Chemical Reaction Engineering Module Users Guide
Chemical Reaction Engineering Module Users Guide
Engineering Module
User’s Guide
Chemical Reaction Engineering Module User’s Guide
© 1998–2018 COMSOL
Protected by patents listed on www.comsol.com/patents, and U.S. Patents 7,519,518; 7,596,474;
7,623,991; 8,219,373; 8,457,932; 8,954,302; 9,098,106; 9,146,652; 9,323,503; 9,372,673; and
9,454,625. Patents pending.
This Documentation and the Programs described herein are furnished under the COMSOL Software License
Agreement (www.comsol.com/comsol-license-agreement) and may be used or copied only under the terms
of the license agreement.
COMSOL, the COMSOL logo, COMSOL Multiphysics, COMSOL Desktop, COMSOL Server, and
LiveLink are either registered trademarks or trademarks of COMSOL AB. All other trademarks are the
property of their respective owners, and COMSOL AB and its subsidiaries and products are not affiliated
with, endorsed by, sponsored by, or supported by those trademark owners. For a list of such trademark
owners, see www.comsol.com/trademarks.
Version: COMSOL 5.4
Contact Information
Visit the Contact COMSOL page at www.comsol.com/contact to submit general
inquiries, contact Technical Support, or search for an address and phone number. You can
also visit the Worldwide Sales Offices page at www.comsol.com/contact/offices for
address and contact information.
If you need to contact Support, an online request form is located at the COMSOL Access
page at www.comsol.com/support/case. Other useful links include:
CONTENTS |3
The Reaction Engineering Interface 58
Features Nodes Available for the Reaction Engineering Interface . . . . . 66
Initial Values . . . . . . . . . . . . . . . . . . . . . . . . 66
Reaction . . . . . . . . . . . . . . . . . . . . . . . . . . 67
Species . . . . . . . . . . . . . . . . . . . . . . . . . . 71
Reversible Reaction Group . . . . . . . . . . . . . . . . . . . 75
Equilibrium Reaction Group. . . . . . . . . . . . . . . . . . . 77
Species Group . . . . . . . . . . . . . . . . . . . . . . . . 78
Additional Source . . . . . . . . . . . . . . . . . . . . . . 79
Reaction Thermodynamics . . . . . . . . . . . . . . . . . . . 80
Species Activity . . . . . . . . . . . . . . . . . . . . . . . 80
Species Thermodynamics. . . . . . . . . . . . . . . . . . . . 80
Feed Inlet. . . . . . . . . . . . . . . . . . . . . . . . . . 81
Generate Space-Dependent Model . . . . . . . . . . . . . . . . 82
Parameter Estimation . . . . . . . . . . . . . . . . . . . . . 90
Experiment . . . . . . . . . . . . . . . . . . . . . . . . . 91
4 | CONTENTS
Coupling to Other Physics Interfaces . . . . . . . . . . . . . . 117
Adding a Chemical Species Transport Interface and Specifying the
Number of Species . . . . . . . . . . . . . . . . . . . . 117
CONTENTS |5
Theory for the Electrophoretic Transport Interface 160
6 | CONTENTS
Open Boundary . . . . . . . . . . . . . . . . . . . . . . 197
Thin Diffusion Barrier . . . . . . . . . . . . . . . . . . . . 197
Thin Impermeable Barrier . . . . . . . . . . . . . . . . . . 197
Equilibrium Reaction . . . . . . . . . . . . . . . . . . . . 198
Surface Reactions . . . . . . . . . . . . . . . . . . . . . 199
Surface Equilibrium Reaction . . . . . . . . . . . . . . . . . 199
Fast Irreversible Surface Reaction . . . . . . . . . . . . . . . 200
Porous Electrode Coupling . . . . . . . . . . . . . . . . . . 200
Reaction Coefficients . . . . . . . . . . . . . . . . . . . . 201
Electrode Surface Coupling . . . . . . . . . . . . . . . . . . 201
Porous Media Transport Properties. . . . . . . . . . . . . . . 202
Adsorption . . . . . . . . . . . . . . . . . . . . . . . . 204
Partially Saturated Porous Media . . . . . . . . . . . . . . . . 205
Volatilization . . . . . . . . . . . . . . . . . . . . . . . 207
Reactive Pellet Bed . . . . . . . . . . . . . . . . . . . . . 208
Reactions. . . . . . . . . . . . . . . . . . . . . . . . . 211
Species Source. . . . . . . . . . . . . . . . . . . . . . . 212
Hygroscopic Swelling . . . . . . . . . . . . . . . . . . . . 213
Fracture . . . . . . . . . . . . . . . . . . . . . . . . . 213
CONTENTS |7
Electrode Surface Coupling . . . . . . . . . . . . . . . . . . 237
Turbulent Mixing . . . . . . . . . . . . . . . . . . . . . . 238
Reaction . . . . . . . . . . . . . . . . . . . . . . . . . 239
Reaction Sources . . . . . . . . . . . . . . . . . . . . . . 240
Initial Values . . . . . . . . . . . . . . . . . . . . . . . 241
Mass Fraction . . . . . . . . . . . . . . . . . . . . . . . 242
Flux . . . . . . . . . . . . . . . . . . . . . . . . . . . 242
Inflow . . . . . . . . . . . . . . . . . . . . . . . . . . 243
No Flux . . . . . . . . . . . . . . . . . . . . . . . . . 244
Outflow . . . . . . . . . . . . . . . . . . . . . . . . . 244
Symmetry . . . . . . . . . . . . . . . . . . . . . . . . 245
Flux Discontinuity . . . . . . . . . . . . . . . . . . . . . 245
Open Boundary . . . . . . . . . . . . . . . . . . . . . . 246
Equilibrium Reaction . . . . . . . . . . . . . . . . . . . . 246
Surface Equilibrium Reaction . . . . . . . . . . . . . . . . . 247
8 | CONTENTS
The Nernst-Planck-Poisson Equations Interface 261
CONTENTS |9
The Reacting Flow Multiphysics Interface 281
The Reacting Laminar Flow Interface . . . . . . . . . . . . . . 281
The Reacting Flow Coupling Feature . . . . . . . . . . . . . . 282
Physics Interface Features . . . . . . . . . . . . . . . . . . 284
Chapter 6: Thermodynamics
10 | C O N T E N T S
External Thermodynamic System . . . . . . . . . . . . . . . . 309
Exporting and Importing Thermodynamic Systems . . . . . . . . . 313
Species Property . . . . . . . . . . . . . . . . . . . . . . 313
Mixture Property. . . . . . . . . . . . . . . . . . . . . . 319
Equilibrium Calculation . . . . . . . . . . . . . . . . . . . 321
Coupling with the Reaction Engineering and the Chemistry Interfaces . . 326
Evaluating a Property Function in a Physics Interface . . . . . . . . 330
User Defined Species . . . . . . . . . . . . . . . . . . . . 331
References . . . . . . . . . . . . . . . . . . . . . . . . 339
Chapter 7: Glossary
CONTENTS | 11
12 | C O N T E N T S
1
This chapter introduces you to the capabilities of the module. A summary of the
physics interfaces and where you can find documentation and model examples is
also included. The last section is a brief overview with links to each chapter in this
guide.
13
About the Chemical Reaction
Engineering Module
In this section:
Included in these physics interfaces are the kinetic expressions for the reacting systems
and models for the definition of mass transport. A variety of ready-made expressions
are also accessible in order to calculate a system’s thermodynamic and transport
properties.
Like all COMSOL modules, the physics interfaces described in this guide include all
the steps available for the modeling process, which are described in detail in the
COMSOL Multiphysics Reference Manual (see Where Do I Access the
Documentation and Application Libraries?), for example:
Once a model is defined, you can go back and make changes to all the branches listed
above, while maintaining consistency in the other definitions throughout. You can
restart the solver, for example, using the existing solution as an initial guess or even
alter the geometry, while the equations and boundary conditions are kept consistent
through the associative geometry feature. It is also useful to review the Introduction
to the Chemical Reaction Engineering Module included with the module’s
documentation.
While a major focus of this module is on chemical reactors and reacting systems, it is
also extensively used for systems where mass transport is the major component. This
includes unit operations equipment, separation and mixing processes, corrosion,
chromatography, and electrophoresis. The module is also widely used for educational
purposes including courses about chemical engineering, chemical reaction
engineering, electrochemical engineering, biotechnology, and transport phenomena.
The notations and structure in this module were inspired by the book
Transport Phenomena by Bird, Stewart, and Lightfoot. The work by H.
Scott Fogler, Elements of Chemical Reaction Engineering, was also used
as an inspiration.
When one or several physics interfaces are chosen from the Model Wizard (or if you
open the Add Study window), you select an analysis type (stationary, dynamic, or
parametric) and then the modeling interface(s) are available as a node(s) in the Model
Builder along with all the other nodes required for modeling (Definitions, Geometry,
and so forth).
By adding another physics interface, you can account for a phenomenon not previously
described in a model. To do this, right-click a Component node in the Model Builder to
open the Add Physics window. You can do this at any stage during the modeling
process. This action still retains the existing geometry, equations, boundary
conditions, and current solution, which you can build upon for further development
of the model.
The table below lists all the physics interfaces specifically available or enhanced with
this module (in addition to the basic COMSOL Multiphysics license). These physics
interfaces form the backbone of the module and are based on the equations for
chemical reaction kinetics, chemical species transport, as well as for fluid flow and heat
transfer in porous media fundamentals.
Single-Phase Flow
• Materials
• In the Model Builder or Physics Builder click a node or window and then
press F1.
• On the main toolbar, click the Help ( ) button.
• From the main menu, select Help>Help.
• Press Ctrl+F1.
• From the File menu select Help>Documentation ( ).
• Press Ctrl+F1.
• On the main toolbar, click the Documentation ( ) button.
• From the main menu, select Help>Documentation.
Once the Application Libraries window is opened, you can search by name or browse
under a module folder name. Click to view a summary of the model or application and
its properties, including options to open it or its associated PDF document.
To include the latest versions of model examples, from the Help menu
select ( ) Update COMSOL Application Library.
The Reacting Flow in Porous Media Multiphysics Interface section, the Reacting Flow
in Porous Media (rfds) and Reacting Flow in Porous Media (rfcs) interfaces are
described. The theory for these physics interfaces is included.
THERMODYNAMICS INPUT
The chapter Thermodynamics describes different ways of how inputting external
thermodynamics information to a model.
In this chapter:
See also Chemical Species Transport Interfaces for the additional physics interfaces
in the Chemical Species Transport branch ( ).
25
Overview of the Reaction Engineering
and Chemistry interfaces
The Reaction Engineering and Chemistry interfaces share many common features and
functionality for defining species and reactions. The main difference is that the Reaction
Engineering interface always is defined within the context of a (zero dimensional)
reactor model, solving for a number of species dependent variables and an optional
heat balance, whereas the Chemistry interfaces solely defines a number of variable
expressions, based on species and reactions, that can be used when coupling the
Chemistry interface to, for instance, a space-dependent mass transport model.
The Reaction Engineering interface can also be used for defining and simulating
thermodynamic data.
• Rate expressions and heat sources for use in mass and heat balances
• Material property variables (mixture density, diffusivities, viscosity, etc.) for use in
space-dependent transport equations.
The Chemistry interface is also created when the Generate Space-Dependent Model
feature is used in The Reaction Engineering Interface, collecting all mixture variables
and properties for use in a space-dependent model.
For each reaction formula entered, specify the reaction type with a delimiter separating
the two sides of the equation:
Figure 2-1: After adding a Reaction node, enter the chemical reaction formula and
specify the reaction type.
Right-click any Reaction node to disable and enable the corresponding node in the
Model Builder. Creating reaction subsets in this way is a straightforward approach to
investigate the influence of individual reactions on the overall reaction system. Species
that take part only in deactivated reactions are automatically deactivated, as indicated
by the unavailable Species feature nodes. If a reaction is deleted, the interface
automatically deletes those species that take part only in the deleted reactions.
Name/Variable scope
ID number
Figure 2-2: The ID number of each Reaction node is displayed in front of the reaction
formula.
As a general labeling rule, the variable name that refers to the contents of a field
associated with a Reaction node is given by the physics interface Name, followed by the
field name, and ends with the reaction ID number. For example, the contents of the
reaction rate field r for Reaction 1 is assigned the variable name re.r_1.
For an example of how to use variables with scope, see Tank Series with
Feedback Control: Application Library path
Chemical_Reaction_Engineering_Module/Ideal_Tank_Reactors/
tankinseries_control
Figure 2-3: Species nodes are generated automatically as chemical reaction formulas are
entered in the Reaction nodes.
As with Reaction features you can add, remove, or deactivate Species features by
right-clicking a node in the Model Builder. Deactivation of a species automatically
deactivates any reactions in which the species is participating.
In the Reaction Rate section it is possible to alter the definition of the reaction rate of
the species. This will override the settings in the Reaction feature, where the rate is
defined by the stoichiometry of the reaction(s) in which the species participates.
It is also possible to add individual Species nodes with either reactive or non-reactive
species.
The Species type selection has implications in the calculation of thermodynamic and
transport properties.
Bulk species and Surface species, defined per reactor volume and area, respectively, set
the mixture’s physical properties dependent upon its composition. However,
configuring a species as a Solvent sets the physical properties of the reacting fluid
mixture equal to those of the solvent species; specifically, its density, heat capacity,
viscosity, and thermal conductivity. The interface also implements a solute-solvent
approximation for the interaction of species in the fluid and describes the transport
properties accordingly. In material balances this means that the diffusion coefficient is
independent of any of the solute’s concentrations, because every solute only interacts
with solvent molecules, regardless of the concentration. In addition, the convective
term in the flux of species is directly given by the velocity field of the solvent multiplied
by the solute concentration.
The Species type has implications on the Generate Space-Dependent Model procedure
since it determines whether interfaces using a solvent-solute approximation of the
reacting fluid mixture (as in The Transport of Diluted Species Interface) or a full
multicomponent description of the reacting fluid mixture (as in The Transport of
Concentrated Species Interface) should be generated.
Furthermore, the Species type affects the reaction kinetics. Solvent sets the species’
concentration to a constant value (the initial species concentration). The Reacting
Engineering interface does not formulate a mass balance for the solvent species. This
setting corresponds to situations where the solvent does not take part in chemical
reactions at all, or where it reacts but is present in large excess.
Furthermore, it applies that adding (ads) directly after the species name creates a
surface species named speciesname_surf. Additionally, adding (s), (l), or (g) creates
either a solid, liquid, or gaseous species denoted speciesname_solid,
speciesname_liquid speciesname_gas, respectively. Note that the fluid Mixture,
gas or liquid, remains the same despite the presence of any of these species in the
system.
Similar to the labeling rule applying to Reaction nodes, the variable name referring to
the contents of a field associated with a Species node is given by the interface Name,
followed by the field name, and ending with the species name. For example, the
contents of the Rate expression field R for the species roh is assigned the variable name
re.R_roh (for the Reaction Engineering interface). Access the definition of all the
variables used by a specific node by displaying the Equation View node. To display the
node, click the Show button ( ) and select Equation View.
f (2-1)
kj
aA + bB + ... r
xX + yY + ...
kj
For such a reaction set, the reaction rates rj (SI unit: mol/(m3·s)) can be described by
the mass action law:
f –ν r ν
rj = kj ∏ c i ij – k j ∏ c i ij (2-2)
i ∈ react i ∈ prod
f r
Here, k j and k j denote the forward and reverse rate constants, respectively. The
concentration of species i is denoted as ci (SI unit: mol/m3). The stoichiometric
coefficients are denoted νij, and are defined as being negative for reactants and positive
for products. In practice, a reaction seldom involves more than two species colliding
in a reacting step, which means that a kinetic expression is usually of order 2 or less
(with respect to the involved concentrations).
Here, A denotes the frequency factor, n the temperature exponent, E the activation
energy (SI unit: J/mol) and Rg the gas constant, 8.314 J/(mol·K). The
pre-exponential factor, including the frequency factor A and the temperature factor
Tn, is given the units (m3/mol)α − 1/s, where α is the order of the reaction (with
respect to volumetric concentrations).
The default settings for the reaction given by Equation 2-1 and assuming equilibrium,
yields the equilibrium expression in Equation 2-4:
x y
c X c Y ...
K eq = -----------------
a b
- (2-4)
c A c B ...
f
r k
k = --------- (2-5)
K eq
G ( T, P ) = U + PV – TS = H – TS (2-6)
where δQ is (reversible) heat transfer to the fluid and δW is (pressure) work in the
system. The change in Gibbs free energy can be written as
where n is the number of mole of species i in the system. At constant temperature, this
expression can be integrated as a function of pressure
μi
μ dμ = P v̂i dP
Pi
° °
(2-10)
i
o
where v is molar volume and μ i is chemical potential of species at standard state. For
an ideal gas this can be expressed as
o Pi
μ i – μ i = RT ln -----o- (2-11)
P
νi μi = 0 (2-12)
i
o fˆi
μ i – μ i = RT ln ----o- (2-13)
f i
where fˆi is the fugacity of species i in the mixture, and f i is the fugacity of pure species
o
fˆi ν i
o
ν i μ i + RT ln ∏ ----o- = 0
f i
(2-14)
i i
– ΔG orxn νi fˆi ν i
K ( T ) = exp -------------------- =
RT ∏ ai = ∏ ----f o- (2-15)
i i i
o o o
ΔG rxn = ΔH rxn – T ΔS rxn (2-16)
o o
the enthalpy of reaction, ΔH rxn , and entropy of reaction, ΔS rxn , both at a given
temperature T is defined as
T
νi ΔH f, i + T Cp, i dT
o o
ΔH rxn ( T ) = (2-17)
o
i
T C p, i
νi Sabs, i + T ----------
o
ΔS rxn ( T ) = - dT (2-18)
T o
i
o o
Here T is the temperature at standard state. ΔH f, i and S abs, i are the standard
enthalpy of formation and absolute entropy for each species (these data are available in
the COMSOL database).
ACTIVITY
Activity of species, ai is defined by Equation 2-13 as
o
μi – μi
a i = exp ----------------- (2-19)
RT
Activity depends on the choice of an arbitrary standard state. The standard state of
pure species is usually at 105 Pa and for solute in solution is based on hypothetical
molality or amount concentration also referred as infinite dilute behavior.
ci
a i = γ x, i x i = γ c, i --------- (2-20)
c ref
Gas phase
The standard state is the pure species at ideal gas condition, 1 atm and the equilibrium
temperature. Activity of species in mixture is expressed by:
where fˆ i and φ̂ i are fugacity and fugacity coefficient of species i in the mixture.
Liquid phase
The standard state is pure liquid species at 1 atm and equilibrium temperature. The
fugacity of a species in a mixture is given by
fˆ i = y i γ i f i (2-22)
where γi is the activity coefficient of species in the mixture and fi is the fugacity of pure
species at the equilibrium temperature and pressure. The activity is expressed by
yi γi fi
a i = ------------
o
- (2-23)
fi
where fio is the fugacity of pure species at the equilibrium temperature and 1 atm. The
o
ratio f i ⁄ f i is given by
f v̂ i
----i = exp -------- ( P – P i )
sat
(2-24)
o RT
fi
where v̂ i is the partial molar volume of species and Psat is species saturated vapor
pressure. For liquids is weak function of pressure and can be assumed to be 1 unless at
high pressure.
Dilute Solutions
The concentration can in the case of non-ideal mixtures be replaced with the activity.
In these interfaces, the dimensionless activity (ai) depend on species concentration (ci),
activity coefficient (γi) and the standard state concentration (c0s=1 mol/m3).
ci
a i = ------- γ i (2-25)
c 0s
Additionally, an effective species concentration (ce,i) (SI unit: mol/m3) is used in the
reaction rates (Equation 2-2) when activities are utilized.
c e, i = c i γ i = a i c 0s (2-26)
ν
∏ ci i
νi i ∈ prod
K eq = ∏ ci = --------------------------
–ν
(2-27)
i ∏ ci i
i ∈ react
Gas Phase
Inserting Equation 2-21 in Equation 2-15 gives:
φ̂ i P ν i
K ( T ) = K eq ∏ -----------------o-
c sum P
(2-28)
i
Liquid Phase
For low and moderate pressure, the equilibrium constant can be reformulated by
substitution of Equation 2-23 to Equation 2-15 as:
γi νi
K ( T ) = K eq ∏ ----------
c sum
- (2-29)
i
EXAMPLE I
The following short example illustrates how the Reaction Engineering interface and
the Chemistry interface handle equilibrium reactions in the formulation of the material
balance equations.
f (2-30)
k
A r B
k
According to Equation 2-2 the reaction rate (SI unit: mol/(m3·s)) is formulated as:
f r
r = k cA – k cB
f r
RA = – k cA + k cB = –r
f r
RB = k cA – k cB = r
f r
r = k cA – k cB = 0
The relationship between the forward and reverse reaction rates in Equation 2-30 is
given by the following ratio:
f cB
eq k
K = ----r- = ------ (2-31)
k cA
The Reaction Engineering interface also sets up mass balances that are solved. The
general material balances for species A and B, respectively, are:
∂c A
= RA = –r (2-32)
∂t
∂c B
= RB = r (2-33)
∂t
The rate of consumption of species A equals the production rate of species B, as shown
in Equation 2-32 and Equation 2-33.
∂ (c + c ) = 0
(2-34)
∂t A B
eq cB
K = ------ (2-35)
cA
EXAMPLE II
This example shows how equilibrium reactions are considered in the Reaction
Engineering interface using the Equilibrium Species Vector section.
f
k1
A B
f (2-36)
k2
B 2C
k2r
d-
----- c = –r1 (2-37)
dt A
d-
----- c = r1 – r2 (2-38)
dt B
d-
----- c = 2r 2 (2-39)
dt C
Now compare Equation 2-37, Equation 2-38, and Equation 2-39 with the balance
equations that the physics interface sets up for the related chemistry, where the second
reaction is instead an equilibrium reaction:
f (2-40)
k1
A B
K2
eq (2-41)
B 2C
d-
----- c = –r1
dt A
d-
----- ( 2c B + c C ) = 2r 1 (2-42)
dt
2
eq cC
K 2 = ------ (2-43)
cB
MASS BALANCE
The mass balances that are set up in the Reaction Engineering interface are simplified
versions of the general mass transport equation. The main assumption is that the
reactor is perfectly mixed, meaning that any variations in compositions within the
reactor are neglected.
The following reactor types are available in the Reaction Engineering interface (the
Chemistry interfaces does not contain reactor models):
BATCH
In the batch reactor no mass enters or leaves the system. Common for all reactor
models is that reacting fluids in the gas phase are assumed to behave as ideal gases.
Liquid mixtures are assumed to be ideal and incompressible.
d ( ci Vr )
-------------------
- = Vr Ri (2-44)
dt
which takes into account the effect of changing volume. In Equation 2-44, ci (SI unit:
mol/m3) is the species molar concentration, Vr (SI unit: m3) denotes the reactor
volume, and Ri (SI unit: mol/(m3·s)) is the species rate expression.
In Equation 2-45, Cp,i (SI unit: J/(mol·K)) is the species molar heat capacity, T (SI
unit: K) is the temperature, and p (SI unit: Pa) the pressure. On the right-hand side,
Q (SI unit: J/s) is the heat due to chemical reaction, and Qext (SI unit: J/s) denotes
heat added to the system. The heat of reaction is:
Q = –Vr Hj rj
j
where Hj (SI unit: J/mol) is the enthalpy of reaction, and rj (SI unit: mol/(m3·s)) the
reaction rate.
For an incompressible and ideally mixed reacting liquid, the energy balance is:
dT
Vr ci Cp, i -------
dt
- = Q + Q ext (2-46)
i
dc i
-------- = R i
dt
For an ideal reacting gas, the energy balance is given by Equation 2-45. For an
incompressible and ideally mixed reacting liquid, the energy balance is given by
Equation 2-46.
The species mass balances for the CSTR are given by:
d ( ci Vr )
-------------------
dt
- =
vf, m cf, m – vci + Ri Vr (2-47)
m
In Equation 2-47, cf,m (SI unit: mol/m3) is the species molar concentration of the
associated feed inlet stream vf,m (SI unit: m3/s). Vr (SI unit: m3) denotes the reactor
volume and is a function of time.
dV
---------r- =
dt v f, m – v + v p (2-48)
m
In Equation 2-48, vp (SI unit: m3/s) denotes the volumetric production rate. It is
given by Equation 2-49 for ideally mixed liquids and by Equation 2-50 for ideal gases.
Ri Mi
vp = Vr -------------
ρi
-, Ri = νij rj (2-49)
i j
Rg T
v p = ----------- V r
p Ri (2-50)
i
where νij is the stoichiometric coefficient of species i in reaction j, Mi (SI unit: kg/
mol) denotes the species molecular weight, ρi (SI unit: kg/m3) the species density, and
Ri (SI unit: mol/(m3·s)) is the reaction rate of species i.
v f, m ρ f, m
m
v = ------------------------------
- (2-51)
ρ
In contrast, when the model is set to be solved for generic conditions a specific outlet
flow stream can be set.
dT dp
Vr ci Cp, i -------
dt
- = Q + Q ext + V r ------- +
dt vf, mi cf, mi ( hf, mi – hi ) (2-52)
i m i
For an incompressible and ideally mixed reacting liquid, the energy balance is:
dT
Vr ci Cp, i -------
dt
- = Q + Q ext + vf, mi cf, mi ( hf, mi – hi ) (2-53)
i m i
d ( ci Vr )
-------------------- =
dt vf, m cf, m – vci + Ri Vr (2-54)
m
dV
---------r- = 0
dt
dc i
V r -------- = v f, i c f, i – vf, i + vp ci + Ri Vr (2-56)
dt
The volumetric production rate, vp, is defined as in Equation 2-49 and Equation 2-50.
The energy balance is the same as for the CSTR with Constant Mass/Generic reactor
type (Equation 2-53).
SEMIBATCH
In the semibatch reactor, reactants enter the reactor by means of one or several feed
inlet streams.
d ( ci Vr )
-------------------
- = v f, i c f, i + R V r (2-57)
dt i
dV
---------r- =
dt vf, m + vp (2-58)
m
An energy balance over the Semibatch reactor results in the same energy balance
expression as for the CSTR reactor types (Equation 2-53).
PLUG FLOW
In the plug flow reactor the species concentrations and the temperature vary with
position. For a tubular reactor configuration, plug flow assumes concentration and
temperature gradients to only develop in the axial direction but not in the radial
direction of the reactor.
dF
---------i = R i (2-59)
dV
where Fi (SI unit: mol/s) is the species molar flow, V (SI unit: m3) is the reactor
volume, and Ri (SI unit: mol/(m3·s)) denotes the species rate expression.
In order to evaluate the rate expressions Ri, which are functions of the species
concentrations, the physics interface calculates:
Fi
c i = -----
v
where Mi (SI unit: kg/mol) denotes the species molecular weight and ρi (SI unit: kg/
m3) the species density
Rg T
v = -----------
p Fi
i
so that
p Fi
c i = ----------- -------------
Rg T
Fi
i
Neglecting pressure drop, the energy balance for an ideal reacting gas, as well as an
incompressible and ideally mixed reacting liquid is given by:
dT
Fi Cp, i -------
dV
- = Q + Q ext (2-60)
i
Equation 2-60 is similar to the energy balance for the batch reactors (Equation 2-46),
but with a reactor volume dependence instead of a time dependence.
Transport Properties
The Reaction Engineering Interface and The Chemistry Interface can calculate several
transport properties that can be accessed in interfaces in space dependent models.
DIFFUSIVITY
The diffusivity is calculated in terms of binary diffusion coefficients. These are available
for the following fluid mixtures:
3 3
– 22 T ( M A + M B ) ⁄ ( 2 ⋅ 10 M A M B ) –1
D AB = 2.6628 ⋅ 10 ⋅ -------------------------------------------------------------------------------------- ⋅ Ω D (2-61)
pσ A σ B
Here, DAB (SI unit: m2/s) is the binary diffusion coefficient, M (SI unit: kg/mol)
equals the molecular weight, T (SI unit: K) represents the temperature, p (SI unit: Pa)
is the pressure, and σ (SI unit: m) equals the characteristic length of the
Lennard-Jones/ Stockmayer potential. In addition, ΩD is the collision integral, given
by the following equation (Ref. 6 and Ref. 7):
–c2
ΩD = c1 ( T ⁄ εk ) + c 3 [ exp ( – c 4 T ⁄ ε k ) ] + c 5 [ exp ( – c 6 T ⁄ ε k ) ]
– 40 2 2 (2-62)
( 4.748 ⋅ 10 )μ D, A μ D, B
+ c 7 [ exp ( – c 8 T ⁄ ε k ) ] + ----------------------------------------------------------------
2 3 3
-
k b Tε k σ A σ B
where
εA εB
εk = ------------ (2-63)
2
kb
In Equation 2-62, cx are empirical constants, μD is the species dipole moment value
(SI unit: Cm) and ε/kb (SI unit: K) the potential energy minimum value divided by
Boltzmann’s constant. Tabulated data in literature frequently lists values of ε/kb rather
than ε. It should be noted that predefined expressions for binary diffusivities only treat
ideal gas mixtures. Thus, these are applicable as input only for gases at moderate
pressure in multicomponent diffusive transport models.
The binary diffusivity according to Equation 2-61 is also suited for gaseous species in
solvent, simply by setting either the component A or B to the solvent. The binary
diffusion coefficient is in this case equal to the diffusion coefficient of the bulk species.
1⁄2
– 15 ( φ B M B ) T
D AB = 3.7·10 ---------------------------------
- (2-64)
0.6
μ B V b, A
DYNAMIC VISCOSITY
The dynamic viscosities are computed for the following fluid mixtures:
3
T ( M i ⋅ 10 )
– 6 ---------------------------------- –1
μ = 2.669 ⋅ 10 2
ΩD (2-65)
σi
In Equation 2-65, μ (SI unit: Ns/m2) represents the dynamic viscosity, and ΩD is the
dimensionless collision integral given by:
–b2
ΩD = b1 ( T ⁄ ( εi ⁄ kb ) ) + b 3 [ exp ( – b 4 T ⁄ ( ε i ⁄ k b ) ) ]
– 40 4 (2-66)
( 4.998 ⋅ 10 )μ D, i
+ b 5 [ exp ( – b 6 T ⁄ ( ε i ⁄ k b ) ) ] + -------------------------------------------------
2 6
-
k b T ( ε i ⁄ k b )σ i
In Equation 2-65 and Equation 2-66, bx are empirical constants, μD (SI unit: Cm) the
species dipole moment value, and ε/kb (SI unit: K) the potential energy minimum
value divided by Boltzmann’s constant. Tabulated data in literature frequently lists
values of ε/kb rather than ε, and σ (SI unit: Å) is the characteristic length value.
n
μ
--------------------------------------------
-
i
μ = (2-67)
j=n
i = 1 1
1 + ----
x i
x j φ ij
j=1, j ≠ i
In Equation 2-67 and Equation 2-68, xi is the molar composition and μi is computed
with Equation 2-65 for each of the species in the mixture.
– 3 .758
3 – 0.2661 T – T ref
μ = 10 ( μ ref 10 )
–3
+ -------------------- (2-69)
223
where μ (SI unit: Ns/m2) is the dynamic viscosity. As inputs for Equation 2-69, the
physics interface takes the reference viscosity, μref (SI unit: Pa⋅s) at the reference
temperature Tref (SI unit: K).
THERMAL CONDUCTIVITY
The thermal conductivity is calculated as well for some types of fluid mixtures:
μ
k = ----- ( 1.15C p + 0.88R g ) (2-70)
M
where k (SI unit: W/(m·K)) is the thermal conductivity and Cp (SI unit: J/(mol·K))
denotes the molar heat capacity. Equation 2-70 is a function of viscosity, μ, as given by
Equation 2-65. Equation 2-70 is directly used in the case of a solvent; all parameters
being those of the solvent. Without a solvent, however, the following equation is also
used:
1 x i – 1
k = ---
2 xi ki + ----k-i
i=1 i=1
where ki is the thermal conductivity of each species i and xi the molar composition for
each of the species in the mixture.
• Solve mass and energy balances for ideal reactor systems ( reactors, semi- reactors,
CSTRs, and plug flow reactors).
• Evaluate species transport properties as a function of reactor conditions.
• Evaluate mass, energy, and momentum balances for space-dependent models, and
transfer them to COMSOL Multiphysics.
It is possible to read the input files for kinetics, thermodynamic, and transport
properties independently and use these as separate data resources. For example, if a set
of reactions is entered into the Reaction Engineering interface, species thermodynamic
and transport data can be supplied by reading the appropriate input files. The full
functionality of the physics interface is retained even after the import procedure. This
means that all expressions and all data imported into the software are available for
reference and for editing.
• Transport CHEMKIN files that supply data used to compute Transport Properties.
• Thermodynamic CHEMKIN files containing data for Gordon and McBride or
NASA polynomials (Ref. 4). These polynomials are denoted NASA format and
compute the species’ heat capacity, Cp, molar enthalpy, h, and molar entropy, s:
2 3 4
C p, i = R g ( a 1 + a 2 T + a 3 T + a 4 T + a 5 T ) (2-71)
a2 2 a3 3 a4 4 a5 5
h i = R g a 1 T + ------ T + ------ T + ------ T + ------ T + a 6 (2-72)
2 3 4 5
a3 2 a4 3 a5 4
s i = R g a 1 ln T + a 2 T + ------ T + ------ T + ------ T + a 7 (2-73)
2 3 4
- Cp,i (SI unit: J/(mol·K)) denotes the species’ molar heat capacity,
- T (SI unit: K) is the temperature,
- Rg the ideal gas constant, 8.314 J/(mol·K),
- hi (SI unit: J/mol) is the species’ molar enthalpy, and
- si (SI unit: J/(mol·K)) represents its molar entropy at standard state.
From the CHEMKIN files the coefficients a1 to a7 are directly imported into the
corresponding NASA format fields. Coefficients for NASA polynomials are available as
public resources (Ref. 12).
Kinetics CHEMKIN files that can be imported in the Reversible Reaction Group
feature. These consist of reaction kinetics data, such as activation energy (SI unit: J/
mol).
CASE 1
Carbon monoxide is part of a reacting mixture. You want to do several things: make
use of the predefined expressions of species Cp as an input to the heat capacity of the
reacting mixture (Equation 2-71); use the predefined expression h for each species to
To accomplish this, enter the seven coefficients of the NASA format, a1 to a7, into the
appropriate fields, or import a CHEMKIN Import for Species Property thermo input
file.
CASE 2
Carbon monoxide is part of a reacting mixture. You want to make use of the
predefined expressions of species Cp to calculate the heat capacity of the reacting
mixture (Equation 2-98).
Polynomials for Cp are available in the literature (Ref. 2, Ref. 3, and Ref. 5) in the
frequently used form
an T
n–1
Cp = Rg n = 1, …, 5
n
You can directly use the predefined expression for the species’ heat capacity,
Equation 2-71, also given in the Cp field, by supplying coefficients in the a1 to a5 fields.
You also want to use the predefined expressions h for each species to calculate the heat
of reaction (Equation 2-72). An option is to make use of the heat of formation at
standard state (298.15 K) to calculate the coefficient a6. Identifying the coefficient a6
of Equation 2-72 is straightforward. In the NASA polynomial format, the species
molar enthalpy is related to its heat capacity according to
T
h = 0 Cp dT + h ( 0 ) (2-74)
Inserting Equation 2-71 into Equation 2-74, and comparing the result with
Equation 2-72, shows that the term a6 Rg is identified as the species enthalpy of
formation at 0 K, that is, h(0). Evaluate h(0) from the species enthalpy of formation
at standard state temperature, Tstd = 298.15 K, which is given by
h ( T std ) a2 2 a3 3 a4 4 a5 5
a 6 = -------------------- – a 1 T std + ------ T std + ------ T std + ------ T std + ------ T std (2-75)
Rg 2 3 4 5
h = C p ( T – 298 ) + h ( 298 )
2. R.H. Perry and D.W. Green, Perry’s Chemical Engineering Handbook, 7th ed.,
McGraw Hill, 1997.
3. B.E. Poling, J.M. Prausnitz, and J.P. O’Connell, The Properties of Gases and
Liquids, 5th ed., McGraw Hill, 2000.
6. P.D. Neufeld, A.R. Janzen, and R.A. Aziz, “Empirical Equations to Calculate 16 of
the Transport Collision Integrals Ω(l,s)*for the Lennard-Jones (12-6) Potential,” J.
Chem. Phys., vol. 57, pp. 1100–1102, 1972.
7. R.S. Brokaw, “Predicting Transport Properties of Dilute Gases,” Ind. Eng. Process
Design Develop., vol. 8, no. 2, pp. 240–253, 1969.
9. C. R. Wilke, “A Viscosity Equation for Gas Mixtures,” J. Chem. Phys., vol. 18, no.
4, pp. 517–519, 1950.
10. W.R. Gambill, “Fused Salt Thermal Conductivity,” Chem. Eng., vol. 66, p. 129,
1959.
11. L.I. Stiel and G. Thodos, “The Viscosity of Polar Substances in the Dense Gaseous
and Liquid Regions,” AIChE J., vol. 10, p. 26–29, 1964.
Add physics features from the toolbar, or right-click Reaction Engineering to select
features from the context menu. Many of the fields and nodes described in this section
are made available when either a Reaction or a Species (or both) subnode is added to
the Model Builder. Because nodes and subnodes are accessible at any time, and any
change is updated throughout the model, reactions and species are often defined
before the settings described in this section.
The following is a description of the features and fields available on the Settings
window for Reaction Engineering.
SETTINGS
The Label is the default physics interface name.
The Name is used primarily as a scope prefix for variables defined by the physics
interface. Refer to such physics interface variables in expressions using the pattern
<name>.<variable_name>. In order to distinguish between variables belonging to
The default Name (for the first physics interface in the model setup) is re.
EQUATION
This section displays the governing equations according to the selection of reactor
types in the Reactor Type section.
REACTOR TYPE
Select a Reactor type to define the reaction system. The available reactor types are:
Batch, Batch, constant volume, CSTR, constant volume, CSTR, constant mass/generic,
Semibatch, and Plug flow. Each reactor type solves a mass balance based on properties
typical to the type, such as reactor volume, volumetric production, or flow rate:
• BatchBatch: In reactors no mass enters or leaves the system. This type can account
for a variable reactor volume.
• Batch, Constant Volume: Same as Batch but where the reactor volume is assumed
to be constant. As this is the situation for most reacting systems, this is the default
condition.
• CSTR, Constant Mass/Generic: Continuous stirred tank reactors (CSTR) differ
from batch reactors, since these allow species to enter and leave the reactor by means
of feed inlet streams and outlet streams. The system’s volume is allowed to change,
such as in a car engine cylinder or a balloon. This reactor model can be solved either
for a constant mass condition or by selecting a specific outlet flow.
• CSTR, Constant Volume: Same as the CSTR with constant mass/generic reactor
but assumes that the volume is constant during operation.
• Semibatch: Semibatch reactors differ from batch reactors in that they allow reactants
to enter the reactor by means of one or several feed inlet streams.
• Plug Flow: In the plug flow reactor, the species concentrations and the temperature
vary with position. Plug flow in a tubular configuration means that concentration
and temperature gradients develop only in the axial direction, but not in the radial
direction.
MIXTURE
Select a Mixture to specify the phase, Gas or Liquid, of the reaction fluid.
Mixture Density
Mixture density always has two settings available: Automatic or User defined. The
mixture density is transferred to physics interfaces set up by the Generate
Space-Dependent Model feature. The density is compiled for both multicomponent
and solute-solvent solutions.
• Automatic selected for Liquid, considers the liquid as ideal and incompressible. The
liquid mixture density depends on the density of i number of pure species (ρi) and
the species mass fraction (wi).
1
ρ = --------------
wi
------
ρi
i
The mass fraction is given by the species concentration (ci) and the molar mass (Mi).
ci Mi
w i = --------------------
ci Mi
i
• Automatic set for Gas calculates the gas mixture density (ρ) from the concentrations
(ci) and molar masses (Mi) of the mixture species, which are automatically taken
from Species features.
• If a Species Type is set to Solvent and the Mixture is Liquid, the mixture density is the
same as the solvent density as defined in Density in General parameters in the Species
node. When Mixture is Gas, the mixture density is calculated from Equation 2-76
only for the species set as Solvent.
Reactor Pressure
The Mixture - Gas setting displays the Reactor pressure. For all reactor types, except the
Plug flow reactor, select either a pressure computed from the Ideal gas law or from any
other expression using the User defined option. The Batch and Semibatch reactor types
also have the option to keep the reactor pressure Constant during reaction.
For the Plug Flow reactor, the Reactor pressure can be entered in the case of the Mass
Balance section having the Volumetric rate set to Automatic, in which case the User
defined alternative is available and a constant pressure fits the conditions.
The most general description of a mixture is the one that treats the mixture as a
multicomponent solution, where all species interact with each other. A simplified
description, but still a common one, assumes that the solution has a solvent that
dominates the properties of the solution. The solutes in such a solution interact only
with solvent molecules.
Select the Calculate mixture properties check box to enable calculation of mixture
transport properties exported from the Reaction Engineering interface.
From the list for each property, select the in-built Automatic expression or set a User
defined entry. The mixture properties you can transfer for space-dependent models are:
• Heat capacity (cp) (SI unit: J/(kg·K)) (this is available when the Energy Balance is
set to Include)
• Mixture density (ρ) (SI unit: kg/m3) (this is available under Mixture Density
section)
All species properties needed to compute the mixture properties are entered in the
Species Transport Expression or Species Thermodynamic Expression in the Species node.
Transport Properties
In the Predefined dependent species (separated by ‘,’) text field, edit, if necessary, the
species that depends on the composition of the other species according to the
Equilibrium expression in the Reaction node. To minimize the impact of any numerical
errors, it is recommended to set the species with the highest concentration as
dependent species. The default species is set to the leftmost species in the Reaction
formula.
ACTIVITY
Select the Use activity check box to solve for species activities instead of species
concentrations, which is a common approach when non-ideal fluids are modeled.
Activity
DISCRETIZATION
The Uniform scaling of concentration variables check box is not selected by default.
When selected, all concentration variables are scaled using the same scale factor in the
Study node. Enabling uniform scaling can decrease solver time for problems involving
many concentration variables.
MASS BALANCE
This section displays the parameters and expressions used in the mass balance equations
for the available reactor types as selected under Reactor type.
TABLE 2-1: REACTOR TYPE PARAMETERS TO DEFINE UNDER MASS BALANCE
Volumetric Rate
For CSTR, constant volume and Semibatch reactor types, the Volumetric production rate
(vp) is available. For Automatic, the in-built Volumetric production rate expression is
shown. If User defined is selected, the expression can be changed (SI unit: m3/s). For
instance, this enables the setting of zero (0) volumetric production rate, which ignores
volume changes due to reactions.
• For liquid phase reactions, the Automatic expression for the Volumetric production
rate varies with the reaction rate of each species as defined by:
The physics interface automatically inserts the stoichiometric coefficients (νij) and
reaction rate expressions for each species (Ri) that depend on j number of reactions
of a rate (rj), as defined in the Reaction feature node. Furthermore, the values of the
molar mass (Mi) and the species density (ρi) are automatically taken from the
Species features.
• For gas phase reactions, the Automatic expression for the Volumetric production rate
is similarly given by:
Rg T
v p = ----------- V r
p R i, Ri = νij rj (2-78)
i j
For CSTR, constant mass/generic, select Volumetric rate to either Constant mass (the
default) or Generic. The Constant mass selection shows both the expression for
Volumetric production rate (vp) (Equation 2-77 or Equation 2-78) and Volumetric
outlet rate (v), the latter is derived from constant mass flow condition through the
reactor:
vf, m ⋅ ρf, m
m
v = ----------------------------------
- (2-79)
ρ
The mixture density (ρf) of m number of feed inlet streams is determined in the
same way as in the Mixture section. For Generic both the Volumetric rate properties
can be edited (SI unit: m3/s). This means that it is possible to completely control
the volumetric outlet rate from the CSTR.
For the Plug flow reactor, the Volumetric flow rate along the reactor is set. The
Automatic definition computes a variable volumetric flow rate that depends on the
molar flow rate of each species (Fi).
• For liquid phase reactions, the expression is:
Fi Mi
v = ------------
ρi
- (2-80)
i
The default value for p (Reactor pressure) is 1 atm in Equation 2-81 and it is set in
the Mixture section.
Reactor Volume
The section sets the Reactor volume Vr , i.e. the fluid volume in which chemical reaction
takes place. The Batch reactor type can account for a changing volume, thus a
time-dependent volume expression can be entered here.
Two types of CHEMKIN input files can be imported here—Thermo or Transport for
either volumetric or surface species. Click Browse to locate the CHEMKIN file to be
imported, then click Import. The imported data directly enters the NASA format fields
in the Species Thermodynamic Expressions section or the Species Transport Expressions.
section.
Initial Values
The Initial Values node sets the initial values of the dependent variables solved for in
the Reaction Engineering interface.
GENERAL PARAMETERS
For CSTR, constant mass and Semi reactor types, in which the reactor volume changes
with reaction, an Initial system volume Vr0 can be set.
For nonisothermal conditions, i.e. the Energy Balance is set to Include, enter an Initial
temperature T0 for the system. For the Plug flow reactor a corresponding Inlet
temperature T0,in is entered.
For the Plug flow reactor an Inlet molar flow rate table is instead shown. Enter the values
or expressions representing the inlet conditions of each species in the Molar flow rate
(mol/s) column.
The total surface concentration initially is always restricted by the molar amount of
available sites on the reactive surface area. This is defined by the Initial density of surface
sites γ0 (SI unit: mol/m2) entry. Make sure that the sum of the entered initial surface
concentrations do not surpass the value selected here. Note that the actual surface
concentration of each species is the Surface concentration (mol/m^2) multiplied with the
Site occupancy number. The restriction of the total surface concentration is defined as:
REACTOR AREA
This section is displayed when surface species are present in the model. Set the Reactor
area, in other words the total reactive surface area available in the reactor, here.
Reaction
To add a Reaction node ( ) either right-click the Reaction Engineering node or on the
Reaction Engineering toolbar click Reaction.
REACTION FORMULA
Formula
Enter a chemical reaction Formula. Click Apply to make the interface examine the
species taking part in the model’s reactions and automatically add the associated
Species features to the Model Builder.
Reaction Type
Select the Reaction type—Reversible, Irreversible, or Equilibrium—or edit the expression
directly in the Formula field. In the latter case, specify the reaction type with a delimiter
separating the two sides of the equation:
• If the reaction is Reversible or Irreversible, the reaction rate for a species becomes:
dc i
-------- = ν i r (2-83)
dt
REACTION RATE
This section is available when the Reaction type is either Reversible or Irreversible.
The Automatic Reaction rate expression is set up according to the mass action law:
f – ν ij
rj = kj ∏ ci (2-85)
i ∈ react
RATE CONSTANTS
This section applies for Reversible or Irreversible reactions and defines the reaction rate
constants used in the reaction rates.
The SI units of the rate constants are automatically based on the order of the reaction
with respect to the concentrations, as defined in the Reaction formula.
r f
k = k ⁄ K eq0 (2-87)
Thus, in this case, the forward rate constant and equilibrium constant for the reaction
are needed. The Equilibrium constant is edited in the Equilibrium Settings section.
f rf
r r n r
k = A ( T ⁄ T ref ) exp ( – E ⁄ ( R g T ) ) (2-89)
Specify the activation energy and the frequency factor in the Arrhenius expressions to
account for temperature variations. The reference temperature, Tref equals 1 K. The
available fields are based on the Reaction type chosen in the Reaction node. Enter
values or expressions for each of the following (reverse expressions are only available
for reversible reactions):
• Forward frequency factor Af and Reverse frequency factor Ar (SI unit: (mβ/
mol)α − 1/s, where α equals the order of the reaction and β is 3 or 2 for volumetric
or surface reactions, respectively)
• Forward temperature exponent nfand Reverse temperature exponent nr
• Forward activation energy Efand Reverse activation energy Er (SI unit: J/mol)
EQUILIBRIUM SETTINGS
This section is available for equilibrium reactions, and for reversible reactions when the
Specify equilibrium constant check box has been selected.
Equilibrium Expression
For an equilibrium reaction, specify the Equilibrium expression. When the Equilibrium
expression is set to Automatic the following expression is used:
Select User defined from the Equilibrium expression list to instead enter a manually
defined equilibrium expression.
Equilibrium Constant
Specify the Equilibrium constant Keq0 for an equilibrium reaction, or for a reversible
reaction when the Specify equilibrium constant check box has been selected (in the Rate
Constants section).
The Equilibrium constant can either be User defined, or automatically defined when set
to Automatic or Thermodynamics.
Use the Automatic option to compute the equilibrium constant for an ideal system.
This setting requires that the temperature is also solved for by setting Energy Balance
to Included.
The Thermodynamics option is available when all reactions in the interface are
equilibrium reactions, and the interface is fully coupled to a Thermodynamic System
(see Species Matching). Use this setting to automatically compute the equilibrium
constant for an ideal or non-ideal system, dependent on the thermodynamic model
applied for the coupled system.
Using Automatic or Thermodynamics, Keq0 is calculated from the Gibbs free energy of
the reaction. For more details see The Equilibrium Constant and the Automatically
Defined Equilibrium Constants section therein.
Enthalpy of Reaction
The Enthalpy of reaction H (SI unit: J/mol) is calculated by the interface from species
properties and the related stoichiometric coefficients:
Entropy of Reaction
The Entropy of reaction S (SI unit: J/(mol·K)) comes from a similar expression:
Qp Qr
In Equation 2-91 and Equation 2-92, hi and si are the species’ molar enthalpy and
molar entropy, respectively.
Enter these quantities in the Species Thermodynamic Expressions section for the Species
node either by using the predefined polynomial or by providing a custom expression
or constants.
The stoichiometric coefficients, vij, are defined as being negative for reactants and
positive for products. Using Equation 2-91 and Equation 2-92 to equate the Gibbs
free energy of reaction enables the equilibrium constant to be expressed according to
Handling of Equilibrium Reactions.
Qj = –Hj rj
Species
When a Reaction is defined, a Species node ( ) is automatically generated for the
participating reactants and products. This feature enables you to review and enter
species specific information regarding chemical kinetics, thermodynamics and
transport properties.
SPECIES NAME
When a Species node is automatically generated using the Formula text field for the
Reaction node, the Species name is also automatically generated.
For a Species node added individually, enter a Species name in the field and click Apply.
SPECIES TYPE
Select a Species type—Bulk species, Surface species, or Solvent.
Bulk species and Solvent are solved for volumetric concentrations (SI unit: mol/m3),
while Surface species are solved for surface concentration (SI unit: mol/m2). The
compositions for Bulk species and Solvent use the syntax c_speciesname, while Surface
species uses csurf_speciesname_surf.
For the Plug flow reactor only Bulk species and Solvent are allowed.
GENERAL PARAMETERS
The General Parameters section deals with species parameters.
Edit, if necessary, the species Molar mass M and the ionic Charge z of the species in
question.
It is possible to specify the species density ρ when the fluid Mixture is specified as Liquid.
The default value is that of water at 293 K.
Edit either the Rate expression (SI unit: mol/(m3s)), the Surface rate expression (SI
unit: mol/(m2s)), or both. For a bulk species, both expressions appear if surface
reactions are present in the reactor since the reaction of the species can depend both
on bulk reaction R and surface reaction Rads rates. For a surface species, only the
surface reaction rate Rads appears.
The reaction rate is not editable if the species in question participates in an equilibrium
reaction and has been selected as a Predefined dependent species in Equilibrium Species
Vector.
SPECIES CONCENTRATION/ACTIVITY
To account for non-ideality in the fluid mixture adjust the activity coefficient in the
Activity coefficient text field here. The section is only shown if activity instead of
concentration has been chosen in the interface, i.e. the Use activities check box is
selected on the Reaction Engineering interface Settings window
Click to select the Locked concentration/activity check box if the species concentration
or activity should be treated as constant.
The Species enthalpy is by default computed with the NASA format. In this case, enter
the following to compute the species’ heat capacity, Cp (SI unit: J/(mol⋅K)), the molar
enthalpy, h (SI unit: J/mol), and the molar entropy, s (SI unit: J/(mol⋅K)):
Transport Properties
Add the node from the Reaction Engineering toolbar or right-click Reaction Engineering
and add it from the context menu.
For the case when reaction kinetics data are entered manually into the Reaction table
and nonisothermal conditions apply (Energy Balance is set to Include), right-click to
add a Reaction Thermodynamics subnode or select it from the Reaction Engineering
toolbar, Attributes menu. In it, the Enthalpy of Reaction (J/mol) for each reaction can
be specified.
REACTION TABLE
The reversible reactions in the Reaction table are numbered and contain reactants,
products, and kinetic parameters describing the reaction. Use the buttons under the
Reaction table to add and sort the reaction details.
• In general, use the Move Up ( ), Move Down ( ), and Delete ( ) buttons and
the fields under tables to edit the table contents. Or right-click a table cell and select
Move Up, Move Down, or Delete.
• The Add button ( ) adds default reactant, A, and product, B, with a default. Click
the corresponding field to edit the reactant, product, or parameters. After editing
the Species Group node is also updated. It is created together with the reaction
group.
In the case of twenty or more reactions the Disable updating variables during editing
table check box is available. Select this to speed up editing text fields; automatic
updates related to edits do not occur until you click to clear the check box.
Add the node from the Reaction Engineering toolbar or right-click Reaction Engineering
and add it from the context menu.
For the case when reaction kinetics data are entered manually into the Reaction table
and nonisothermal conditions apply (Energy Balance is set to Include), right-click to
add a Reaction Thermodynamics subnode or select it from the Reaction Engineering
toolbar, Attributes menu. In it, the Enthalpy of Reaction (J/mol) for each reaction can
be specified.
• In general, use the Move Up ( ), Move Down ( ), and Delete ( ) buttons and
the fields under tables to edit the table contents. Or right-click a table cell and select
Move Up, Move Down, or Delete.
• The Add button ( ) adds default reactant, A, and product, B, with a default. Click
the corresponding field to edit the reactant, product, or parameters. After editing
the Species Group node is also updated. It is created together with the reaction
group.
• You can save the parameters to a text file to reuse in other models. Click the Save to
File button ( ) and enter a File name in the Save to File dialog box, including the
extension.txt. Click Save to store the parameters in a text file or in a Microsoft Excel
Workbook spreadsheet if the license includes LiveLink™ for Excel®. The information
is saved in space-separated columns in the same order as displayed on screen. When
saving to Excel, an Excel Save dialog box appears where you can specify the sheet and
range and whether to overwrite existing data, include a header, or use a separate
column for units.
• You can import or load data in files from a spreadsheet program, for example, with
the Load from File button ( ) and the Load from File dialog box that appears. Data
must be separated by spaces or tabs. If there is already data in the table, imported
parameters are added after the last row. Move or edit rows as needed. If the license
includes LiveLink™ for Excel® you can also load parameters from a Microsoft Excel
In the case of twenty or more reactions the Disable updating variables during editing
table check box is available. Select this to speed up editing text fields; automatic
updates related to edits do not occur until you click to clear the check box.
Species Group
The Species Group node ( ) contains information on a molecular level about the
volumetric species and the surface species present in the model. The Property for
Volumetric Species or Property for Surface Species tables typically collect parameters
from when importing CHEMKIN transport files.
Add the node from the Reaction Engineering toolbar or right-click Reaction Engineering
and add it from the context menu. This node is also automatically added when either
the Reversible Reaction Group or the Equilibrium Reaction Group are used.
For the case of nonisothermal reactor conditions (Energy Balance is set to Include), a
Species Thermodynamics subnode is automatically created in which the
thermodynamic properties of the species can be specified.
Additional Source
Use the Additional Source node ( ) to add an additional rate expression (SI unit:
mol/m ) and/or an additional volumetric production rate (SI unit: m3/s) to the mass
3
Add the node from the Reaction Engineering toolbar or right-click Reaction Engineering
and add it from the context menu.
Add the Reaction Thermodynamics node from the Reaction Engineering toolbar,
Attributes menu. Alternatively, when the Energy Balance is set to Include, right-click a
Reversible Reaction Group or Equilibrium Reaction Group and select it from the
context menu.
Species Activity
The Species Activity node ( ) creates variables for the activities for all the species/
surface species present in the Species Group parent feature. Edit the Activity coefficient
field in the Species Activity or Surface Species Activity tables by clicking in these.
Species Thermodynamics
The Species Thermodynamics node ( ) creates variables for the enthalpies, entropies,
and heat capacities for all the species/surface species present in the Species Group
parent feature. The purpose is to compute thermodynamic mixture properties and the
heat of reactions.
This node is a subnode to the Species Group node when the Energy Balance is set to
Include.
Feed Inlet
The Feed Inlet ( ) feature is used for adding inlet streams to the reactor.
After adding a Reaction node and setting its Reactor Type to CSTR, constant volume,
CSTR, constant mass/generic, or Semi, add a Feed Inlet node from the Reaction Engineering
toolbar or right-click the Reaction Engineering node to add it from the context menu.
The Feed inlet temperature (Tf) is required as input when nonisothermal conditions are
investigated (Energy Balance is set to Include).
To add this feature, on the Reaction Engineering toolbar click Generate Space-Dependent
Model ( ) or right-click the Reaction Engineering node to add it from the context
menu. Only one Generate Space-Dependent Model node can be added per model file.
1 Select the component to use: 1D axi, 1D, 2D, 2D axi, 3D, or an existing component
within the model builder.
2 Select the physics interfaces.
3 Select a study type.
Figure 2-5: The Generate Space-Dependent Model transfers the properties from the
Reaction Engineering interface to physics interfaces in a1D axi, 1D, 2D, 2D axi, or 3D
component.
If more interfaces are needed in a component after the initial model generation, select
the component and the new interface(s) and click Create/refresh to introduce the
interface into it.
COMPONENT SETTINGS
Select a Component to use. Either specify the space dimension of a new component—
1D, 1Daxi, 2D, 2Daxi, and 3D—or select a component already defined in the Model
Builder.
If desired, physics interfaces for fluid flow, heat transfer, or other features affecting the
reacting system can be added later separately from the physics interfaces created by the
Generate Space-Dependent Model feature.
If necessary, variables and properties in the Chemistry node can be accessed with the
Chemistry node Name (default Chem). An arbitrary reaction rate expression, R, is thus
When surface species are present (that is, when the Species Type is set to Surface species
for at least one species in the reactor), specify how to include surface reactions in the
space-dependent model.
• When Boundaries is selected, the surface reactions will be modeled using a Surface
Reactions interface defined on the boundaries of the geometry.
• When Porous Pellets is selected, the surface reactions will be modeled using a
Reactive Pellet Bed feature added to a Transport of Diluted Species interface or a
Diluted Species in Porous Media interface. In this case the content of the Mass balance
list is restricted to these interfaces.
A corresponding reaction feature is added and set up, in accordance with the reaction
kinetics defined in the Reaction Engineering interface, when clicking the Create/Refresh
button in the Space-Dependent Model Generation section. In Figure 2-7, the surface
reaction kinetics in a Reacting Engineering interface has been implemented in a Reactive
Pellet Bed feature using a Reactions subfeature. Note that the surface reaction rates are
defined by the Chemistry interface.
Heat Transfer
The Reaction Engineering interface can also be set up with time- and space-dependent
energy balance equations. In the model generation process the physics interface
generates expressions used by the Heat Transfer interface, such as the heat generated
by a chemical reaction. It also generates expressions for physical transport properties.
There are two heat transfer interfaces under the Energy balance:
Several expressions for the species density, heat capacity, and thermal conductivity are
available and can be transferred from the Reaction Engineering interface.
The densities are available from the Mixture section in the interface, where the fluid
mixture properties are selected. The density depends on the Species settings and is
computed as follows for:
• Liquid, if assuming the liquid to be ideal and incompressible. The liquid mixture
density depends on the density of i number of pure species, ρi, and the species
weight fraction, wi.
The volume fraction is given by the species concentration, ci, and the molar mass,
Mi.
ci Mi
w i = -------------------- (2-94)
ci Mi
i
• Gas mixture density depends on the concentrations and molar masses of the mixture
species.
ρ = ci Mi (2-95)
i
For mixtures with solvent all values are taken from the species set as solvent.
The heat capacity, cp (SI unit: J/(mol·kg)), of the mixture is calculated by the
species’ molar heat capacity, Cp (SI unit: J/(mol·kg)) according to
C p, i
cp = wi ----------
Mi
- (2-96)
i
where M is the molar mass (SI unit: kg/mol) and wi the weight fraction.
Fluid Flow
The model generation process can be selected to generate a Fluid Flow interface.
Species density (see Equation 2-93, Equation 2-94, and Equation 2-95) and dynamic
viscosity can for some mixture options be transferred from the Reaction Engineering
interface to the Fluid Flow interfaces.
Transport Properties
STUDY TYPE
Select a Study Type, either Stationary or Time Dependent. It is also possible to edit this
choice after the space-dependent model has been generated.
Parameter Estimation
The Parameter Estimation feature ( ) requires the Optimization Module.
Add the node from the Reaction Engineering toolbar or right-click Reaction Engineering
and add it from the context menu. The node and related subnodes have the options
and settings for the definition of parameter estimation calculations based on the
predefined reactor types in the Reaction Engineering interface.
ESTIMATION PARAMETERS
The parameters to be estimated are defined in the Estimation parameters table. The
parameter names are entered in the Parameter column. An Initial value is also required
for each variable. To connect the estimated parameters to the model, make sure that
the same parameter naming is used within the rest of the Reaction Engineering
interface.
Bounds for the parameter values are specified in the Lower bound and Upper bound
columns of the Estimation parameters table. These are applicable when the SNOPT
optimization (Solver Study Steps for Parameter Estimation) solver method is used. The
default solver method for parameter estimation is Levenberg-Marquardt, which
currently does not support parameter bounds. Enter a Scale indicating a typical
magnitude of the control variable. The relative solver tolerances refer to variables
rescaled with respect to this scale, and it may also influence the search pattern of some
optimization solvers.
Right-click the Parameter Estimation node to add one or more Experiment nodes. You
can also select it from the Reaction Engineering toolbar, Attributes menu. The
Experiment node requires the Optimization Module.
EXPERIMENTAL DATA
Experimental data is read into the software by importing comma-separated value files,
CSV files. After import, identifiers display in the Data column of the Experimental Data
table (Figure 2-8). If the CSV files have column headers, these are used to identify the
columns. If the data file does not have column headers, the columns are identified by
a number. It is assumed that the first column corresponds to the independent variable
for time (t). In the Model variables column, the measured data is correlated with model
variables by entering the proper variable names or expressions. The Use column
controls whether to include data in the parameter estimation calculation or not. Data
column entries require associated Model variables to be accounted for the parameter
estimation calculations.
In the Unit column, enter the unit of the imported experimental data.
In the Weight column, enter a strictly positive value. The difference between the value
of the Model variable and the Data column value is squared and multiplied with the
Weight and a factor 0.5 to give the contribution to the total objective for each
measured value.
EXPERIMENTAL PARAMETERS
Each Experiment node can be associated with user-defined parameters during the
parameter estimation. For example, concentration transients can be measured for a
number of different isothermal conditions. Parameter name and Parameter expression
are added from a list. The list is populated with the parameters that are defined in the
Parameters under Global Definitions in the Model Builder.
INITIAL STATES
A given experiment usually determines the initial values for the dependent variables in
the reactor model, such as the initial concentrations for isothermal reactor models and
initial concentrations and temperature for non-isothermal reactor models.
Initial conditions that change from one experiment to another need to be provided in
the Initial States section of the Experiment node. An entry in the Parameter names
column is added by choosing from a list of predefined variable names corresponding
to the initial state of the reactor model.
The Chemistry (chem) interface ( ) is found under the Chemical Species Transport
branch ( ) when adding a physics interface. The Chemistry interface is also created
when the Generate Space-Dependent Model feature is used in the Reaction
Engineering interface, collecting all mixture variables and properties for use in a
space-dependent model.
This physics interface is a tool for generating a set of variables to be used for modeling
chemical species and reactions systems. The variables are generated from species and
reaction properties and can be divided in two categories:
• Rate expressions and heat sources for use in mass and heat balances
• Material property variables (mixture density, diffusivities, viscosity, etc.) for use in
space-dependent transport equations.
Many of the fields and nodes described in this section are only made available when
either a Reaction or a Species (or both) subnode is added to the Model Builder. All
predefined constants and expressions can be overwritten by user-defined expressions.
This makes it possible to go beyond the modeling assumptions set as defaults in this
physics interface.
The following is a description of the features and fields available on the Settings
window for the Chemistry interface.
SETTINGS
The Label is the default physics interface name.
The Name is used primarily as a scope prefix for variables defined by the physics
interface. Refer to such physics interface variables in expressions using the pattern
The default Name (for the first physics interface in the model) is chem.
MODEL INPUTS
This section sets the Temperature and Pressure. In the lists, it is possible to edit the User
defined option to the temperature and pressure from interfaces available in the Model
Builder. The latter choice couples the interfaces.
SPECIES MATCHING
Use the table in this section to specify the concentrations for the variables generated
by the Reaction and the Species features (the reaction rate for example), and to define
mixture properties (transport and thermodynamic properties).
Select Diluted Species in the Mixture type list to use the concentration variables from a
Transport of Diluted Species interface. The same setting should be used for any other
interface solving for species concentrations using a diluted species assumption. Enter
the names of the dependent variables in the Molar concentration column.
Select Concentrated Species in the Mixture type list in order to use the mass fractions
from a Transport of Concentrated Species interface. Enter the names of the dependent
variables in the Mass fraction column.
The name of the reaction rate variables generated by the interface can be seen in the
Reaction Rate column of the table. The syntax of the reaction rate variables depend on
the Species type selected in the Species Settings window. Bulk species and Solvent use
the syntax R_speciesname, while Surface species uses R_speciesname_surf.
Selecting Concentrated species for the Mixture type, the base of syntax for the reaction
rate is changed. to Rw_speciesname.
EXTRA DIMENSION
Select Define variables in extra dimension when the current Chemistry interface is
coupled to a feature on an extra dimension. An example of this is when the
concentrations in the Concentration input column of the Model Inputs, Concentration
corresponds to pellet concentrations from a Reactive Pellet Bed feature in a Transport
of Diluted Species interface.
When selected, generated variables will be defined using concentrations averaged over
the extra dimension. Note that the generated variables are global in order be available
MIXTURE
Select a Mixture to specify the phase, Gas or Liquid, of the reaction fluid.
Mixture Density
Mixture density always has two settings available: Automatic or User defined. The density
is compiled for both multicomponent and solute-solvent solutions.
• Automatic (ideal liquid) selected for Liquid, considers the liquid as ideal and
incompressible. The liquid mixture density depends on the density of i number of
pure species (ρi) and the species weight fraction (wi).
1
ρ = --------------
wi
------
ρi
i
The volume fraction is given by the species concentration (ci) and the molar mass
(Mi).
ci Mi
w i = --------------------
ci Mi
i
• Automatic set for Gas calculates the gas mixture density (ρ) from the concentrations
(ci) and molar masses (Mi) of the mixture species, which are automatically taken
from Species features.
ρ = ci Mi (2-97)
i
• If a Species Type is set to Solvent and the Mixture is Liquid, the mixture density is the
same as the solvent density as defined in Density in General parameters in the Species
node. When Mixture is Gas, the mixture density is calculated from Equation 2-97
only for the species set as Solvent.
ACTIVITY
Select the Use activity check box to solve for species activities instead of species
concentrations, which is a common approach when non-ideal fluids are modeled.
Activity
Transport Properties
REACTION FORMULA
Formula
Enter a chemical reaction Formula. Click Apply to make the interface examine the
species taking part in the model’s reactions and automatically add the associated
Species features to the Model Builder.
Reaction Type
Select the Reaction type — Reversible, Irreversible, or Equilibrium — or edit the
expression directly in the Formula field. In the latter case, specify the reaction type with
a delimiter separating the two sides of the equation:
• If the reaction is Reversible or Irreversible, the reaction rate for a species is:
Ri = νi r (2-98)
REACTION RATE
This section is available when the Reaction type is either Reversible or Irreversible.
The Automatic Reaction rate expression is set up according to the mass action law:
f – ν ij
rj = kj ∏ ci
i ∈ react
RATE CONSTANTS
This section applies for Reversible or Irreversible reactions and defines the reaction rate
constants used in the reaction rates.
The SI units of the rate constants are automatically based on the order of the reaction
with respect to the concentrations, as defined in the Reaction formula.
r f
k = k ⁄ K eq0
Thus, in this case, the forward rate constant and equilibrium constant for the reaction
are needed. The Equilibrium constant is edited in the Equilibrium Settings section.
f rf
r r n r
k = A ( T ⁄ T ref ) exp ( – E ⁄ ( R g T ) )
Specify the activation energy and the frequency factor in the Arrhenius expressions to
account for temperature variations. The reference temperature, Tref equals 1 K. The
available fields are based on the Reaction type chosen in the Reaction node. Enter
• Forward frequency factor Af and Reverse frequency factor Ar (SI unit: (mβ/
mol)α − 1/s, where α equals the order of the reaction and β is 3 or 2 for volumetric
or surface reactions, respectively)
• Forward temperature exponent nf and Reverse temperature exponent nr
• Forward activation energy Ef and Reverse activation energy Er (SI unit: J/mol)
EQUILIBRIUM SETTINGS
This section is available for equilibrium reactions, and for reversible reactions when the
Specify equilibrium constant check box has been selected.
Equilibrium Expression
For an equilibrium reaction, specify the Equilibrium expression. When the Equilibrium
expression is set to Automatic the following expression is used:
νij
∏ ci
i ∈ prod
K eqj = ----------------------------
– ν ij
∏ ci
i ∈ react
Select User defined from the Equilibrium expression list to instead enter a manually
defined equilibrium expression.
Equilibrium Constant
Specify the Equilibrium constant Keq0 for an equilibrium reaction, or for a reversible
reaction when the Specify equilibrium constant check box has been selected (in the Rate
Constants section).
The Equilibrium constant can either be User defined, or automatically defined when set
to Automatic or Thermodynamics.
Use the Automatic option to compute the equilibrium constant for an ideal system.
This settings requires that the Calculate thermodynamic properties in the corresponding
section is selected.
The Thermodynamics option is available when all reactions in the interface are
equilibrium reactions, and the interface is fully coupled to a Thermodynamic System
(see Species Matching). Use this setting to automatically compute the equilibrium
Using Automatic or Thermodynamics, Keq0 is calculated from the Gibbs free energy of
the reaction. For more details see The Equilibrium Constant and the Automatically
Defined Equilibrium Constants section therein.
Enthalpy of Reaction
The Enthalpy of reaction H (SI unit: J/mol) is calculated by the interface from species
properties and the related stoichiometric coefficients:
Hj = ν ij h i – ( – ν ij )h i (2-100)
i ∈ prod i ∈ react
Entropy of Reaction
The Entropy of reaction S (SI unit: J/(mol·K)) comes from a similar expression:
Sj = ν ij s i – ( – ν ij )s i (2-101)
i ∈ prod i ∈ react
In Equation 2-100 and Equation 2-101, hi and si are the species’ molar enthalpy and
molar entropy, respectively.
Enter these quantities in the Species Thermodynamic Expressions section for the Species
node either by using the predefined polynomial or by providing a custom expression
or constants.
The stoichiometric coefficients, vij, are defined as being negative for reactants and
positive for products. Using Equation 2-100 and Equation 2-101 to equate the Gibbs
free energy of reaction enables the equilibrium constant to be expressed according to
Equation 2-100.
Species
When a Reaction is defined, a Species node ( ) is automatically generated for the
participating reactants and products. This feature enables you to review and enter
species specific information regarding chemical kinetics, thermodynamics and
transport properties.
It is also possible to add and define an individual Species node: on the Chemistry toolbar
click Species or right-click the Chemistry node and select it from the context menu.
SPECIES NAME
When a Species node is automatically generated using the Formula text field for the
Reaction node, the Species name is also automatically generated.
For a Species node added individually, enter a Species name in the field and click Apply.
SPECIES TYPE
Select a Species type—Bulk species, Surface species, or Solvent.
Bulk species and Solvent are solved for volumetric concentrations (SI unit: mol/m3),
while Surface species are solved for surface concentration (SI unit: mol/m2). The
compositions for Bulk species and Solvent use the syntax c_speciesname, while Surface
species uses csurf_speciesname_surf.
When Surface species is selected, the corresponding reaction formula introduces (ads)
after the species notation and changes the species’ name to speciesname_surf.
Additionally, the Species node name is updated in a similar fashion.
GENERAL PARAMETERS
The General Parameters section deals with species parameters.
It is possible to specify the species density ρ when the fluid Mixture is specified as Liquid.
The default value is that of water at 293 K.
REACTION RATE
Change the Automatic default setting to User defined to use a species reaction rate other
than the one set up in the associated Reaction node. For individual species, use the User
defined option to set a reaction rate other than zero (that is, non-reactive).
Edit either the Rate expression (SI unit: mol/(m3·s)), the Surface rate expression (SI
unit: mol/(m2·s)), or both. For a bulk species, both expressions appear if surface
reactions are present in the reactor since the reaction of the species can depend both
on bulk reaction R and surface reaction Rads rates. For a surface species, only the
surface reaction rate Rads appears.
ADDITIONAL SOURCE
The Additional Source section is available in order to include additional rate
contribution for the species to the reaction kinetics. When the Additional source check
box is selected, add an Additional rate expression in the text field (SI unit: mol/m3).
SPECIES CONCENTRATION/ACTIVITY
To account for non-ideality in the fluid mixture adjust the activity coefficient in the
Activity coefficient text field here. The section is only shown if activity instead of
concentration has been chosen in the interface, i.e. the Use activities check box is
selected on the Chemistry interface Settings window
Click to select the Locked concentration/activity check box if the species concentration
or activity should be treated as constant.
• For a gas mixture, there are maximum five properties to consider: σ, the
characteristic length (unit: Å) of the Lennard-Jones/Stockmayer potential; ε/kb, the
energy minimum (SI unit: K) of the Lennard-Jones/Stockmayer potential; μD, the
dipole moment (SI unit: Debye); ki, the thermal conductivity of the gas (SI unit:
The Species enthalpy is by default computed with the NASA format. In this case, enter
the following to compute the species’ heat capacity, Cp (SI unit: J/(mol⋅K)), the molar
enthalpy, h (SI unit: J/mol), and the molar entropy, s (SI unit: J/(mol⋅K)):
Any coefficients for the thermodynamic polynomials entered into the alow,k fields
apply to the temperatures in the range Tlo to Tmid; coefficients entered into the ahi,k
fields apply to temperatures in the range Tmid to Thi range. The coefficients can also
be imported in the CHEMKIN Import for Species Property section in the Chemistry node.
Add the node from the Chemistry toolbar or right-click Chemistry and add it from the
context menu.
For the case when reaction kinetics data are entered manually into the Reaction table
and temperature dependent reaction kinetics apply (Calculate Thermodynamic
Properties check box is selected), right-click to add a Reaction Thermodynamics
subnode or select it from the Chemistry toolbar, Attributes menu. In it, the Enthalpy of
Reaction (J/mol) for each reaction can be specified.
REACTION TABLE
The reversible reactions in the Reaction table are numbered and contain reactants,
products, and kinetic parameters describing the reaction. Use the buttons under the
Reaction table to add and sort the reaction details.
• In general, use the Move Up ( ), Move Down ( ), and Delete ( ) buttons and
the fields under tables to edit the table contents. Or right-click a table cell and select
Move Up, Move Down, or Delete.
• The Add button ( ) adds default reactant, A, and product, B, with a default. Click
the corresponding field to edit the reactant, product, or parameters. After editing
the Species Group node is also updated. It is created together with the reaction
group.
• You can save the parameters to a text file to reuse in other models. Click the Save to
File button ( ) and enter a File name in the Save to File dialog box, including the
extension.txt. Click Save to store the parameters in a text file or in a Microsoft Excel
Workbook spreadsheet if the license includes LiveLink™ for Excel®. The information
is saved in space-separated columns in the same order as displayed on screen. When
saving to Excel, an Excel Save dialog box appears where you can specify the sheet and
range and whether to overwrite existing data, include a header, or use a separate
column for units.
• You can import or load data in files from a spreadsheet program, for example, with
the Load from File button ( ) and the Load from File dialog box that appears. Data
must be separated by spaces or tabs. If there is already data in the table, imported
parameters are added after the last row. Move or edit rows as needed. If the license
In the case of twenty or more reactions the Disable updating variables during editing
table check box is available. Select this to speed up editing text fields; automatic
updates related to edits do not occur until you click to clear the check box.
Add the node from the Chemistry toolbar or right-click Chemistry and add it from the
context menu.
For the case when reaction kinetics data are entered manually into the Reaction table
and temperature dependent reaction kinetics apply (Calculate Thermodynamic
Properties check box is selected), right-click to add a Reaction Thermodynamics
subnode or select it from the Chemistry toolbar, Attributes menu. In it, the Enthalpy of
Reaction (J/mol) for each reaction can be specified.
• In general, use the Move Up ( ), Move Down ( ), and Delete ( ) buttons and
the fields under tables to edit the table contents. Or right-click a table cell and select
Move Up, Move Down, or Delete.
• The Add button ( ) adds default reactant, A, and product, B, with a default. Click
the corresponding field to edit the reactant, product, or parameters. After editing
the Species Group node is also updated. It is created together with the reaction
group.
• You can save the parameters to a text file to reuse in other models. Click the Save to
File button ( ) and enter a File name in the Save to File dialog box, including the
extension.txt. Click Save to store the parameters in a text file or in a Microsoft Excel
Workbook spreadsheet if the license includes LiveLink™ for Excel®. The information
is saved in space-separated columns in the same order as displayed on screen. When
saving to Excel, an Excel Save dialog box appears where you can specify the sheet and
range and whether to overwrite existing data, include a header, or use a separate
column for units.
• You can import or load data in files from a spreadsheet program, for example, with
the Load from File button ( ) and the Load from File dialog box that appears. Data
must be separated by spaces or tabs. If there is already data in the table, imported
parameters are added after the last row. Move or edit rows as needed. If the license
includes LiveLink™ for Excel® you can also load parameters from a Microsoft Excel
Workbook spreadsheet. Then an Excel Load dialog box appears where you can specify
the sheet and range and whether to overwrite existing data. It is also possible to
import from a spreadsheet containing a separate column for units.
In the case of twenty or more reactions the Disable updating variables during editing
table check box is available. Select this to speed up editing text fields; automatic
updates related to edits do not occur until you click to clear the check box.
Add the node from the Chemistry toolbar or right-click Chemistry and add it from the
context menu. This node is also automatically added when either the Reversible
Reaction Group or the Equilibrium Reaction Group are used.
CHEMKIN
This section allows import of CHEMKIN transport files. To import the data directly
into the table columns, click Browse to locate the CHEMKIN file to be imported, then
click Import. Note that this section is only shown when transport properties are
computed (Calculate mixture properties check box under the Calculate Transport
Properties is selected).
Reaction Thermodynamics
The Reaction Thermodynamics subnode ( ), the Enthalpy of Reaction (J/mol) of
each reaction can be specified. This node overrides all the automatically calculated
reaction enthalpies as defined in the Species Thermodynamics subnode.
Add the Reaction Thermodynamics node from the Chemistry toolbar, Attributes menu.
Alternatively, when the Calculate Thermodynamic Properties check box is selected on
the Chemistry interface Settings window, right-click a Reversible Reaction Group or
Equilibrium Reaction Group to add the Reaction Thermodynamics subnode.
Species Thermodynamics
The Species Thermodynamics node ( ) creates variables for the enthalpies, entropies,
and heat capacities for all the species/surface species present in the Species Group
parent feature. The purpose is to compute thermodynamic mixture properties and the
heat of reactions.
This node is a subnode to the Species Group node when the Calculate
Thermodynamic Properties check box is selected on the Chemistry interface Settings
window.
The predefined Plug flow reactor type sets up equations for stationary conditions
describing how species molar flow rates (SI unit: mol/s) vary with the reactor volume.
Stationary Plug Flow is the correct study type to use for this reactor type. The reactor
volume is specified in the Stationary Plug Flow study step feature
Note that selecting the Reaction Engineering interface and the Time Dependent study sets
up the Batch, constant volume reactor type by default. The Time Dependent study is also
suitable for all CSTR reactors and the Semibatch reactor. However, if the reactor type is
changed to Plug flow, then the Time Dependent study step needs to be replaced with
a Stationary Plug Flow step. Combining a Reaction Engineering interface with the
Stationary Plug Flow study does not set up a Plug flow reactor type by default. In that
situation, the Plug Flow reactor type needs to be set manually in the Reaction
Engineering interface.
Figure 2-9: The settings for the Optimization Solver are accessed in the Study tree.
www.comsol.com/optimization-module
113
• The Reacting Flow Multiphysics Interface
• The Reacting Flow in Porous Media Multiphysics Interface
Note that The Chemistry and Reaction Engineering Interfaces and The Chemistry
Interface are described in The Chemistry and Reaction Engineering Interfaces chapter.
The Transport of Diluted Species in Porous Media Interface ( ) is tailored for the
modeling of solute transport in porous media. The physics interface supports cases
where either the solid phase substrate is exclusively immobile, or when a gas-filling
medium is also assumed to be immobile.
Most chemical reactions or other type of material processing, such as casting, either
require or produce heat, which in turn affects both the reaction and other physical
processes connected to the system. Heat transfer through conduction and convection,
when flow is laminar, as well as through porous media are supported in this module.
More extensive description of heat transfer, such as in turbulent flow or involving
radiation, require the Heat Transfer Module. The Heat Transfer Module also includes
turbulent flow.
2 Under Chemical Species Transport, navigate to the physics interface to add and
double-click it.
There are other ways to add a physics interface depending on whether you are in the
Model Builder or Add Physics window:
- In the Model Wizard, click Add or right-click and select Add Physics ( ). The
physics interface displays under Added physics interfaces.
- In the Add Physics window, click Add to Component ( ) or right-click and select
Add to Component.
3 Specify the number of species (concentrations or mass fractions) and the names:
- In the Model Wizard, on the Review Physics Interface page under Dependent
Variables.
- In the Add Physics window, click to expand the Dependent Variables section.
- After adding the physics interface, you can also edit this information—click the
node in the Model Builder, then, on the Settings window under Dependent
Variables, specify the information.
4 Continue by adding more physics interfaces and specifying the number of species
(concentrations or mass fractions) that are to be simulated in a mass transport
interface when adding that interface.
5 In the Dependent Variables section, enter the Number of species. To add a single
species, click the Add Concentration button ( ) underneath the table or enter a
value into the Number of species field. Click the Remove Concentration button ( )
underneath the table if required.
6 The Transport of Concentrated Species interface and the Nernst-Planck Equations
interface both need to contain at least two species (the default). Also edit the strings
or names directly in the table. The names must be unique for all species (and all
other dependent variables) in the model.
7 For the Nernst-Planck Equations interface, an electric potential variable is also
required.
When studying mixtures that are not dilute, the mixture and transport properties
depend on the composition, and a different physics interface is recommended. See The
Transport of Concentrated Species Interface for more information.
Fick’s law governs the diffusion of the solutes, dilute mixtures, or solutions, while the
phenomenon of ionic migration is sometimes referred to as electrokinetic flow. The
Transport of Diluted Species interface supports the simulations of chemical species
transport by convection, migration, and diffusion in 1D, 2D, and 3D as well as for
axisymmetric components in 1D and 2D.
In this section:
∂c i
------- + ∇ ⋅ J i + u ⋅ ∇c i = R i (3-1)
∂t
Equation 3-1 in its form above includes the transport mechanisms diffusion and
convection. If Migration in Electric Field is activated (only available in some add-on
products), the migration mechanism will be added to the equation as well. See more
details in the section Adding Transport Through Migration.
The mass flux relative to the mass averaged velocity, Ji (SI unit: mol/(m2·s)), is
associated with the mass balance equation above and used in boundary conditions and
flux computations. The Transport of Diluted Species interface always includes mass
transport due to molecular diffusion. In this case the mass flux Ji defines the diffusive
flux vector
J i = – D ∇c (3-2)
An input field for the diffusion coefficient is available. Anisotropic diffusion tensor
input is supported.
The third term on the left side of Equation 3-1 describes the convective transport due
to a velocity field u. This field can be expressed analytically or obtained from coupling
the physics interface to one that solves for fluid flow, such as Laminar Flow. Note that
all fluid flow interfaces solve for the mass averaged velocity.
On the right-hand side of the mass balance equation (Equation 3-1), Ri represents a
source or sink term, typically due to a chemical reaction or desorption on a porous
matrix. To specify Ri, another node must be added to the Transport of Diluted Species
interface — the Reaction node for example, which includes an input field for specifying
a reaction expression using the variable names of all participating species.
The kinetics of the reaction is so fast that the equilibrium condition is fulfilled at all
times in all space coordinates.
ν
∏ ai i
i ∈ products
K eq = ----------------------------------
–ν
∏ ai i
i ∈ reactants
ci
a i = γ c, i -------
c a0
where ca0 (SI unit: mol/m3) is the standard molarity, and γc,i (dimensionless) an
activity coefficient.
νi
K eq = ∏ ai
i
The Equilibrium Reaction node solves for a reaction rate so that the equilibrium
condition is always fulfilled in the domain.
There are two ways to present a mass balance where chemical species transport occurs
through diffusion and convection. These are the nonconservative and conservative
formulations of the convective term:
∂c
nonconservative: ----- + u ⋅ ∇c = ∇ ⋅ J i + R (3-3)
∂t
∂c
conservative: ----- + ∇ ⋅ ( cu ) = ∇ ⋅ J i + R (3-4)
∂t
and each is treated slightly differently by the solver algorithms. In these equations
Ji (SI unit: mol/(m2·s)) is the diffusive flux vector, R (SI unit: mol/(m3·s)) is a
production or consumption rate expression, and u (SI unit: m/s) is the solvent velocity
field. The diffusion process can be anisotropic, in which case D is a tensor.
If the conservative formulation is expanded using the chain rule, then one of the terms
from the convection part, c∇·u, would equal zero for an incompressible fluid and
would result in the nonconservative formulation above. This is in fact the default
formulation in this physics interface. To switch between the two formulations, click the
Show button ( ) and select Advanced Physics Options.
∂----c-
= ∇ ⋅ Ji + R
∂t
Note: The features below are only available in a limited set of add-on products. For a
detailed overview of which features are available in each product, visit
http://www.comsol.com/products/specifications/
POINT SOURCE
·
A point source is theoretically formed by assuming a mass injection/ejection, Q c (SI
3
unit: mol/(m ·s)), in a small volume δV and then letting the size of the volume tend
to zero while keeping the total mass flux constant. Given a point source strength, q· p,c
(SI unit: mol/s), this can be expressed as
· ·
lim
δV → 0 Qc = qp,c (3-5)
δV
q· p,c test ( c )
·
is added at a point in the geometry. As can be seen from Equation 3-5, Q c must tend
to plus or minus infinity as δV tends to zero. This means that in theory the
concentration also tends to plus or minus infinity.
Observe that “point” refers to the physical representation of the source. A point source
can therefore only be added to points in 3D components and to points on the
symmetry axis in 2D axisymmetry components. Other geometrical points in 2D
components represent physical lines.
LINE SOURCE
·
A line source can theoretically be formed by assuming a source of strength Q l,c (SI
3
unit: mol/(m ·s)), located within a tube with cross section δS and then letting δS tend
to zero while keeping the total mass flux per unit length constant. Given a line source
strength, q· l,c (SI unit: mol/(m·s)), this can be expressed as
·
lim Ql,c = q· l,c (3-6)
δS → 0
δS
As in the point source case, an alternative approach is to assume that mass is injected/
extracted through the surface of a small object. This results in the same mass source,
but requires that effects resulting from the physical object’s volume are neglected.
q· l,c test ( c )
As with a point source, it is important not to mesh too finely around the line source.
For feature node information, see Line Mass Source and Point Mass
Source in the COMSOL Multiphysics Reference Manual.
• The Line Mass Source node is available as two nodes, one for the fluid
flow (Fluid Line Source) and one for the species (Species Line Source).
• The Point Mass Source node is available as two nodes, one for the fluid
flow (Fluid Point Source) and one for the species (Species Point Source).
Note: Migration is only available in a limited set of add-on products. For a detailed
overview of which features are available in each product, visit
http://www.comsol.com/products/specifications/
∂c i
+ ∇ ⋅ ( – D i ∇c i – z i u m, i F c i ∇V + c i u ) = R i (3-7)
∂t
where
The velocity, u, can be a computed fluid velocity field from a Fluid Flow interface or
a specified function of the spatial variables x, y, and z. The potential can be provided
by an expression or by coupling the system of equations to a current balance, such as
the Electrostatics interface. Sometimes it is assumed to be a supporting electrolyte
present, which simplifies the transport equations. In that case, the modeled charged
species concentration is very low compared to other ions dissolved in the solution.
Thus, the species concentration does not influence the solution’s conductivity and the
net charge within the fluid.
The Nernst-Einstein relation can in many cases be used for relating the species mobility
to the species diffusivity according to
Di
u m, i = --------
RT
where R (SI unit: J/(mol·K)) is the molar gas constant and T (SI unit: K) is the
temperature.
Note: In the Nernst-Planck Equations interface, the ionic species contribute to the
charge transfer in the solution. It includes an electroneutrality condition and also
computes the electric potential field in the electrolyte. For more information, see
Theory for the Nernst-Planck Equations Interface. This interface is included in the
Chemical Reaction Engineering Module.
Supporting Electrolytes
In electrolyte solutions, a salt can be added to provide a high electrolyte conductivity
and decrease the ohmic losses in a cell. These solutions are often called supporting
electrolytes, buffer solutions, or carrier electrolytes. The added species, a negative and
a positive ion pair, predominates over all other species. Therefore, the supporting
electrolyte species can be assumed to dominate the current transport in the solution.
In addition, the predominant supporting ions are usually selected so that they do not
react at the electrode surfaces since the high conductivity should be kept through the
process, that is, they should not be electro-active species. This also means that the
concentration gradients of the predominant species in a supporting electrolyte are
usually negligible.
The current density vector is proportional to the sum of all species fluxes as expressed
by Faraday’s law:
i = F zi Ni
i
The electroneutrality condition ensures that there is always a zero net charge at any
position in a dilute solution. Intuitively, this means that it is impossible to create a
current by manually pumping positive ions in one direction and negative ions in the
other. Therefore, the convective term is canceled out to yield the following expression
for the electrolyte current density, where j denotes the supporting species:
–zj um, j F cj ∇φ
2
i = F (3-8)
j
Equation 3-8 is simply Ohm’s law for ionic current transport and can be simplified to
i = – κ ∇φ (3-9)
where κ is the conductivity of the supporting electrolyte. A current balance gives the
current and potential density in the cell
∇⋅i = 0
∇ ⋅ ( – κ ∇φ ) = 0 (3-10)
Equation 3-10 can be easily solved using the Electrostatics or Secondary Current
Distribution interface and, when coupled to the Transport in Diluted Species interface,
the potential distribution shows up in the migration term.
Crosswind Diffusion
Transport of diluted species applications can often result in models with a very high
cell Péclèt number—that is, systems where convection or migration dominates over
diffusion. Streamline diffusion and crosswind diffusion are of paramount importance
to obtain physically reasonable results. The Transport of Diluted Species interface
In some cases, the resulting nonlinear equation system can be difficult to converge.
This can happen when the cell Péclèt number is very high and the model contains
many thin layers, such as contact discontinuities. You then have three options:
CODINA
The Codina formulation is described in Ref. 1. It adds diffusion strictly in the direction
orthogonal to the streamline direction. Compared to the do Carmo and Galeão
formulation, the Codina formulation adds less diffusion but is not as efficient at
reducing over- and undershoots. It also does not work as well for anisotropic meshes.
The advantage is that the resulting nonlinear system is easier to converge and that
underresolved gradients are less smeared out.
Use the Danckwerts condition to specify inlet concentrations to domains where high
reaction rates are anticipated in the vicinity to the inlet (Ref. 2).
Given an inlet concentration ci,0, the Danckwerts inflow boundary condition reads
n ⋅ ( J i + uc i ) = n ⋅ ( uc i, 0 ) (3-11)
Inflow
∂ θc ∂
( i ) + ∂ ( ρc P, i ) + (a v c G, i) + u ⋅ ∇c i =
∂t ∂t ∂t (3-12)
∇ ⋅ [ ( D D, i + D e, i ) ∇c i ] + R i + S i
On the left-hand side of Equation 3-12, the first three terms correspond to the
accumulation of species within the liquid, solid, and gas phases, while the last term
describes the convection due to the velocity field u (SI unit: m/s).
In Equation 3-12 ci denotes the concentration of species i in the liquid (SI unit: mol/
m3), cP, i the amount adsorbed to solid particles (moles per unit dry weight of the
solid), and cG, i the concentration of species i in the gas phase.
The equation balances the mass transport throughout the porous medium using the
porosity εp, the liquid volume fraction θ; the matrix (drained) density, ρ = (1 − εp)ρp,
and the solid phase density ρp.
For saturated porous media, the liquid volume fraction θ is equal to the porosity εp,
but for partially saturated porous media, they are related by the saturation s as θ = sεp.
The resulting gas volume fraction is av = εp − θ = (1-s)εp.
On the right-hand side of Equation 3-12, the first term introduces the spreading of
species due to mechanical mixing resulting from the porous media (dispersion), as well
as from diffusion and volatilization to the gas phase. The dispersion tensor is denoted
DD (SI unit: m2/s) and the effective diffusion by De (SI unit: m2/s).
The last two terms on the right-hand side of Equation 3-12 describe production or
consumption of the species; Ri is a reaction rate expression which can account for
Adsorption
The time evolution of the adsorption, the solute transport to or from the solid phase,
is defined by assuming that the amount of solute adsorbed to the solid, cP,i, is a
function of the concentration in the fluid ci. This implies that the solute concentration
in the liquid and solid phase are in instant equilibrium. The adsorption term can be
expanded to give
∂c P, i ∂c i ∂ε P ∂c i ∂ε P
∂ ρc = ρk P,i – ρc P, i (3-13)
( P, i ) = ρ – ρc P, i ∂ t ∂t
∂t ∂ ci ∂ t ∂t
Volatilization
Volatilization is the process where a solute species in the liquid is transported to the
gas phase due to vaporization. Assuming that the amount of solute in the gas phase,
cG,i, is a linear function of the liquid phase concentration, the volatilization term is
defined as
∂c G, i ∂c i ∂a v ∂c i ∂a v
∂a = av + k G, i c i = a v k G, i + k G, i c i (3-14)
v c G, i ∂ ci ∂ t ∂t ∂ t ∂t
∂t
∂ ε
( c ) + ∂ ( ρc P, i ) + u ⋅ ∇c i = ∇ ⋅ [ ( D D, i + D e, i ) ∇c i ] + R i + S i (3-15)
∂t p i ∂t
The velocity field to be used in the Model Inputs section on the physics
interface can, for example, be prescribed using the velocity field from a
Darcy’s Law or a Brinkman Equations interface.
The average linear fluid velocities ua, provides an estimate of the fluid velocity within
the pores:
u
u a = ----- Saturated
εp
u
u a = ---- Partially saturated
θ
where εp is the porosity and θ = sεp the liquid volume fraction, and s the saturation, a
dimensionless number between 0 and 1.
Figure 3-1: A block of a porous medium consisting of solids and the pore space between the
solid grains. The average linear velocity describes how fast the fluid moves within the pores.
The Darcy velocity attributes this flow over the entire fluid-solid face.
If the conservative formulation is expanded using the chain rule, then one of the terms
from the convection part, ci∇·u, would equal zero for an incompressible fluid and
would result in the nonconservative formulation described in Equation 3-12.
When using the nonconservative formulation, which is the default, the fluid is assumed
incompressible and divergence free: ∇ ⋅ u = 0. The nonconservative formulation
improves the stability of systems coupled to a momentum equation (fluid flow
equation).
To switch between the two formulations, click the Show button ( ) and
select Advanced Physics Options. In the section Advanced Settings select
either Nonconservative form (the default) or Conservative form. The
conservative formulation should be used for compressible flow.
De = DF Free Flow
εp
D e = ----- D L Saturated Porous Media
τL
θ
D e = ----- D L Partially Saturated Porous Media
τL
θ av
D e = ----- D L + ------ k G D G Partially Saturated with Volatilization
τL τG
Here DF, DL, and DG are the single-phase diffusion coefficients for the species diluted
in fluid, pure liquid and gas phases respectively (SI unit: m2/s), and τF, τL, and τG are
the corresponding tortuosity factors (dimensionless).
The tortuosity factor accounts for the reduced diffusivity due to the fact that the solid
grains impede Brownian motion. The interface provides predefined expressions to
–7 ⁄ 3 2 –7 ⁄ 3 2
τL = θ ε , τG = av ε
–5 ⁄ 2 2 –5 ⁄ 2 2
τL = θ ε , τG = av ε
For saturated porous media θ = εp. The fluid tortuosity for the Millington and Quirk
model is
–1 ⁄ 3
τL = εp
–1 ⁄ 2
τL = εp
User defined expressions for the tortuosity factor can also be applied.
Dispersion
The contribution of dispersion to the mixing of species typically overshadows the
contribution from molecular diffusion, except when the fluid velocity is very small.
The spreading of mass, as species travel through a porous medium is caused by several
contributing effects. Local variations in fluid velocity lead to mechanical mixing
referred to as dispersion occurs because the fluid in the pore space flows around solid
particles, so the velocity field varies within pore channels. The spreading in the
direction parallel to the flow, or longitudinal dispersivity, typically exceeds the
transverse dispersivity from up to an order of magnitude. Being driven by the
concentration gradient alone, molecular diffusion is small relative to the mechanical
dispersion, except at very low fluid velocities.
is controlled through the dispersion tensor DD. The tensor components can either be
given by user-defined values or expressions, or derived from the directional
dispersivities.
Using the longitudinal and transverse dispersivities in 2D, the dispersivity tensor
components are (Ref. 9):
2 2
ui uj
D Dii = α L ------- + α T -------
u u
ui uj
D Dij = D Dji = ( α L – α T ) -----------
u
In these equations, DDii (SI unit: m2/s) are the principal components of the
dispersivity tensor, and DDji and DDji are the cross terms. The parameters αL and αT
(SI unit: m) specify the longitudinal and transverse dispersivities; and ui (SI unit: m/
s) stands for the velocity field components.
In order to facilitate modeling of stratified porous media in 3D, the tensor formulation
by Burnett and Frind (Ref. 10) can be used. Consider a transverse isotropic media,
where the strata are piled up in the z direction, the dispersivity tensor components are:
In Equation 3-17 the fluid velocities u, v, and w correspond to the components of the
velocity field u in the x, y, and z directions, respectively, and α1 (SI unit: m) is the
longitudinal dispersivity. If z is the vertical axis, α2 and α3 are the dispersivities in the
transverse horizontal and transverse vertical directions, respectively (SI unit: m).
Setting α2 = α3 gives the expressions for isotropic media shown in Bear (Ref. 9 and
Ref. 11).
Adsorption
As species travel through a porous medium they typically attach to (adsorb), and
detach (desorb) from the solid phase, which slows chemical transport through the
porous medium. Adsorption and desorption respectively reduces or increases species
concentrations in the fluid. The adsorption properties vary between chemicals, so a
plume containing multiple species can separate into components (Ref. 6). The
Adsorption feature includes three predefined relationships to predict the solid
concentrations, cPi from the concentration in the liquid phase, ci:
∂c P
cP = KP c k P = -------- = K P User defined
∂c
c N ∂c P cP
c P = K F -------- k P = -------- = N ------ Freundlich (3-18)
c ref ∂c c
KL c ∂c P K L c Pmax
c P = c Pmax -------------------- k P = -------- = ---------------------------2- Langmuir
1 + KL c ∂c ( 1 + KL c )
• These predefined expressions are adsorption isotherms that describe the amount of
species sorbed to the solid. Defined at equilibrium, the switch between liquid and
solid phases is instantaneous.
Using a Species Source feature, arbitrary expressions can be entered to define, for
example, nonequilibrium and temperature-dependent adsorption laws, including
those set out by Fetter (Ref. 7) and Bear and Verruijt (Ref. 8).
The retardation factor, RF, describes how adsorption slows the solute velocity, uc,
relative to the average linear velocity of the fluid, ua, as in
ρ b ∂c P ua
RF = 1 + ------ -------- = ------
θ ∂c uc
If the contaminant moves at the average linear velocity of the fluid for RF = 1. For
RF > 1, the contaminant velocity is smaller than the fluid velocity owing to residence
time on solids.
Reactions
Chemical reactions of all types influence species transport in porous media. Examples
include biodegradation, radioactive decay, transformation to tracked products,
temperature- and pressure-dependent functions, exothermic reactions, and
endothermic reactions. The reactions represent change in species concentration per
unit volume porous medium per time. Reaction terms are used on the right-hand side
of the governing equation to represent these processes. For reactions in a fluid phase,
multiply the expression by the fluid volume fraction θ. Similarly, solid phase reaction
expressions include the bulk density, ρb, and gas phase reactions include the gas
volume fraction, av.
The fluid flow in a fracture can be modeled using Darcy’s law formulated in a thin
sheet of porous medium (a fracture):
κ
u = --- ∇ t p
μ
Here u is the tangential Darcy velocity, κ is the fracture permeability, μ the fluid’s
dynamic viscosity, and ∇tp is the tangential gradient of the fluid pressure.
The equation to solve for mass transport of species ci in a thin fracture, embedded in
a porous media, is derived from Equation 3-12. The resulting equation is:
∂ρ b c P, i ∂ε p c i
d fr ------------------- + ------------- + ∇ t ⋅ ( D e, i ∇ t c i ) + u ⋅ ∇ t c i = d fr R i + d fr S i + n 0 (3-19)
∂t ∂t
Here dfr is the fracture thickness, cP, i the amount of species adsorbed to (or desorbed
from) the porous matrix (moles per unit dry weight of the solid), εp is the fracture
porosity, and De is the effective diffusivity. The first two terms on the right hand side
represent source terms from reactions, and n0 corresponds to out-of plane flux from
the adjacent porous domain.
In order to arrive at the tangential differential equation, the gradient is split into the
contributions normal and tangential to the fracture:
∇c i = ∇ n c i + ∇ t c i
The normal gradient is defined in the direction normal to the boundary representing
the fracture and the tangential gradient is defined along the boundary. Assuming that
the variations in the normal (thin) direction of the fracture are negligible compared to
those in the tangential direction, the gradient is simplified as:
∇c i = ∇ t c i
See Fracture for more information about the boundary feature solving
Equation 3-19. See The Transport of Diluted Species in Fractures
Interface for more information about the physics interface solving the
equation on boundaries only.
For an example of how to use the Reactive Pellet Bed, see the model
example A Multiscale 3D Packed Bed Reactor, file path
Chemical_Reaction_Engineering_Module/Reactors_with_Porous_Catalysts/
packed_bed_reactor_3d
Macroscale: c
Concentration in
fluid passing through
bed
Microscale:
cpe Concentration
in porous pellet
Inflow
Figure 3-3: Schematic showing the macroscale (bed volume) and the microscale (pellet).
The transport and reaction equations inside the pellets are solved on an extra
dimension attached to the 1D, 2D, or 3D physics interfaces, including axisymmetric
cases.
The equations inside the spherical pellet are solved as a spherical transport equations
on a nondimensional radial coordinate on the domain 0-1. A given pellet size or a
discrete distribution of sizes can be used.
The model equations assume spherical particles (pellets) of a radius rpe. Modeling
assumptions for cylinders, flakes, and user-defined shapes can also be used. Consider
the microscale concentration cpe inside an individual porous pellet or pellets, and the
macro-concentration c in the packed bed gas volume.
• One uniform pellet radius, which can be space dependent f (x, y, z).
• Binary, ternary, and so on, mix up to 5 radii. The user inputs a table with the mix of
sizes (1 mm, 2 mm, for example), and a percentage of each. Different chemical
reactions per pellet size can be specified.
The model equation for the bulk (macroscale) species is, for example:
∂
ε b (c i) + u ⋅ ∇c i + ∇ ⋅ ( – D b, i ∇c i ) = R i (3-20)
∂t
• The dependent variable c for each chemical species i represents the interstitial
concentration (SI unit: mol/m3), that is, the physical concentration based on unit
volume of fluid flowing between the pellets.
• εb is the bed porosity (SI unit: 1). It should be noted that the R term on the right
side is per unit volume of bed, (SI unit: mol/(m3· s)).
0 1
Figure 3-4: Modeling domain in a pellet for dimensional (top) and nondimensional
(bottom) coordinates.
A shell mole balance across a spherical shell at radius rdim (SI unit: m), and a
subsequent variable substitution r = rdim/rpe gives the following governing equation
inside the pellet on the domain 0<r<1:
where
ε pe D
D pe = ------------ .
τ
The available tortuosity models for porous media are the Millington and Quirk
(Ref. 12),
–1 ⁄ 3 4⁄3
τ = ε pe D pe = ε pe D, (3-22)
Bruggeman,
–1 ⁄ 2 3⁄2
τ = ε pe D pe = ε pe D (3-23)
and the Tortuosity model, where the tortuosity expression is entered as user input:
ε pe D
D pe = ------------ (3-24)
τ
These are readily used for both gaseous and liquid fluids along with various types of
pellet shapes. For instance, the first model has been shown to fit mass transport in
soil-vapor and soil-moisture well.
Equation 3-21 can be solved for two types of boundary conditions at the interface
between the pellet surface and the fluid in this feature.
• Continuous concentration: assuming that all resistance to mass transfer to/from the
pellet is within the pellet and no resistance to pellet-fluid mass transfer is on the bulk
fluid side. The concentration in the fluid will thus be equal to that in the pellet pore
where Ni, inward is the molar flux from the free fluid into a pellet and has the unit
moles/(m2· s).
With the film resistance formulation above, the free fluid Equation 3-20 needs to be
amended for flux continuity so that
∂
ε b (c i) + u ⋅ ∇c i + ∇ ⋅ ( – D b, i ∇c i ) = R i – N i,inward S b (3-26)
∂t
where Sb (SI unit: m2/m3) is the specific surface area exposed to the free fluid of the
packed bed (not including the inside of the pores).
For the case of randomly packed spherical pellets, the specific surface area exposed to
the free fluid is (Ref. 3):
3
S b = ------- ( 1 – ε b ) (3-27)
r pe
The mass transfer coefficient in Equation 3-25 can be computed from the fluid
properties and flow characteristics within the porous media. For this, the Sherwood,
Sh, number defined as the ratio between the convective mass transfer coefficient and
the diffusive mass transfer coefficient is often used:
hL
Sh = --------
D
where L is a characteristic length (for spheres, typically the radius), and D is the
diffusion coefficient in the fluid. From the Sherwood number definition, the mass
transfer coefficient can be computed.
Three commonly used empirical expressions for the calculation of the Sherwood
number are the Frössling relation (Ref. 4):
1⁄2 1⁄3
Sh = 2 + 0.552Re Sc , (3-28)
1⁄2 1⁄3
Sh = 0.94Re Sc , (3-30)
All three depend on the Reynolds, Re, and Schmidt, Sc, numbers. The first describing
the fluid flow regime (laminar versus turbulent) and the second, the ratio between the
viscous diffusion rate and the molecular (mass) diffusion rate. In the expressions,
properties such as velocity, u, dynamic viscosity, μ, and density, ρ, of the fluid are
included.
ρuL
Re = -----------
μ
μ
Sc = --------
ρD
SURFACE SPECIES
Surface species corresponds to species bound to the solid interface of the porous pellet,
which is also in contact with the bulk fluid. The surface species hence only exists within
the pellet. The surface species are assumed immobile, and the concentration is only
dependent on the surface reaction rate:
∂c pe,s
= R pe,s (3-31)
∂t
It can be noted that there is no mass flux of surface species within the pellet or across
the pellet outer surface and the bulk fluid.
NONSPHERICAL PARTICLES
For nonspherical pellets (of any shape), the relations above can be applied
approximately by reinterpreting the pellet radius rpe as
V pe
r pe = 3 --------- (3-32)
A pe
A pe
S pe = ---------
V pe
3
S pe = ------- .
r pe
For a packed bed of which the packing has a porosity εb, the specific surface of the bed
will be
A pe
S b = S pe ( 1 – ε b ) = --------- ( 1 – ε b )
V pe
or
3
S b = ------- ( 1 – ε b ) (3-33)
r pe
for any pellet shape. Now rpe and Sb can be calculated for any shape and inserted in
equations Equation 3-21 and Equation 3-26. Some common specific shapes have
automatic support:
Cylinders
For cylindrical shapes, applying Equation 3-32 gives
2
V pe πr cyl L cyl 1
r pe = 3 --------- = 3 -------------------------------------------------
- = -------------------------------
-. (3-34)
A pe 2 2 2 -
2πr cyl L cyl + 2πr cyl ------------ + ------------
3r cyl 3L cyl
It is common practice to assume that the top and bottom surface of cylindrical pellets
have negligible effect on the mass transfer to and from the internals of the pellet, or,
r cyl « L cyl . Equation 3-34 then simplifies to
3
r pe = --- r cyl
2
Flakes
The derivation for a disc-shaped catalyst pellet is exactly the same as for cylindrical
pellets, except that the assumption is reversed about the end surfaces and the envelope
surface: r flake » w flake , where wflake is the thickness of the disc. This gives
3
r pe = --- w flake
2
and
2
S b = -------------- ( 1 – ε b ) .
w flake
SURFACE REACTION
Surface reactions can also be simulated inside the pellet. Surface species are introduced
into pellet by adding them in the Surface Species section of the Reactive Pellet Bed
feature.
A bulk species can take part in both volumetric and surface reactions. The total
reaction rate for a bulk species within a pellet is defined as:
R = R pe + R pe,s ⋅ S b,reac
Here Rpe is the reaction rate for bulk reactions occurring inside the pellet. Rpe,s and
Sb are the reaction rate and the reactive specific surface area for a surface reaction
occurring inside the pellet, on the interface between the solid matrix and the fluid.
HEAT SOURCE
The heat source of endothermic or exothermic reactions inside the pellet needs to be
accounted for in the heat transfer on the bulk level. Thermal equilibrium is assumed in
each pellet, implying that the heat balance is not solved within the pellet. To account
for the heat evolution in the bulk, the source is averaged across the pellet:
Qpe dV
Q bed = V
--------------------- ( 1 – ε b ) W/m3
V pe
If there are multiple pellet sizes i in the bed the heat source computed by summing
over all sizes:
References
1. R. Codina, “A discontinuity-capturing crosswind-dissipation for the finite element
solution of the convection-diffusion equation”, Computer Methods in Applied
Mechanics and Engineering, vol. 110, pp. 325–342, 1993.
3. J.M. Coulson and J.F. Richardson, Chemical Engineering, vol. 2, 4th ed.,
Pergamon Press, Oxford, U.K., 1991.
4. J.M. Coulson and J.F. Richardson, Chemical Engineering, vol. 1, 4th ed.,
Pergamon Press, Oxford, U.K., 1991.
6. D.M. Mackay, D.L. Freyberg, P.V. Roberts, and J.A. Cherry, “A Natural Gradient
Experiment on Solute Transport in a Sand Aquifer: 1. Approach and Overview of
Plume Movement”, Water Resourc. Res., vol. 22, no. 13, pp. 2017–2030, 1986.
10. R.D. Burnett and E.O. Frind, “An Alternating Direction Galerkin Technique for
Simulation of Groundwater Contaminant Transport in Three Dimensions: 2.
Dimensionality Effects”, Water Resour. Res., vol. 23, no. 4, pp. 695–705, 1987.
11. J. Bear, Dynamics of Fluids in Porous Media, Elsevier Scientific Publishing, 1972.
12. R.J. Millington and J.M. Quirk, “Permeability of Porous Solids”, Trans. Faraday
Soc., vol. 57, pp. 1200–1207, 1961.
∂
( ρω i ) + ∇ ⋅ ( ρω i u ) = – ∇ ⋅ j i + R i (3-35)
∂t
where, ρ (SI unit: kg/m3) denotes the mixture density and u (SI unit: m/s) the mass
averaged velocity of the mixture. The remaining variables are specific for each of the
species, i, being described by the mass transfer equation:
The relative mass flux vector ji can include contributions due to molecular diffusion,
mass flux due to migration in an electric field, and thermal diffusion.
Summation of the transport equations over all present species gives Equation 3-36 for
the conservation of mass
assuming that
Q Q Q
ωi = 1 , ji = 0 , Ri = 0
i=1 i=1 i=1
Using the mass conservation equation, the species transport for an individual species,
i, is given by:
ρ ∂ ( ω i ) + ρ ( u ⋅ ∇ )ω i = – ∇ ⋅ j i + R i (3-37)
∂t
Q − 1 of the species equations are independent and possible to solve for using
Equation 3-37. To compute the mass fraction of the remaining species, COMSOL
Multiphysics uses the fact that the sum of the mass fractions is equal to 1:
ω1 = 1 – ωi (3-38)
i=2
D̃ik dk – Di ∇ ln T
T
j i = – ρω i (3-39)
k=1
In Equation 3-39:
˜
• D 2
ik (SI unit: m /s) are the multicomponent Fick diffusivities
where
When an external force applies equally to all species (such as for gravity),
the last two terms disappear.
As can be seen in Equation 3-39 and Equation 3-40, the total diffusive flux for the
species depends on the gradients of all species concentrations, temperature, and
pressure as well as any external force on the individual species.
Using the ideal gas law, p = c·Rg·T, and the definition of the partial pressures, pk = xkp,
the equation can be written as
Q
1
d k = ∇x k + --- ( x k – ω k ) ∇p – ρω k g k + ω k
p ρωl gl (3-41)
l=1
ωk
x k = -------- M (3-42)
Mk
Q
1- ωi
----
M
= ------
Mi
-
i=1
Q
ρ ∂ ( ω i ) + ρ ( u ⋅ ∇ )ω i = ∇ ⋅ ρω i ˜ d + DT ∇ T-
∂t
D ik k i -------
T
+ Ri
k=1
(3-43)
Q
1
d k = ∇x k + --- ( x k – ω k ) ∇p – ρω k g k + ω k
p ρωl gl
l=1
Multicomponent Diffusivities
˜
The multicomponent Fick diffusivities, D ik , are needed to solve Equation 3-43. The
diffusivities are symmetric
˜ ˜
D ik = D ki
and are related to the multicomponent Maxwell-Stefan diffusivities, Dik, through the
following relation (Ref. 2)
xi xk
( adjB i ) jk
j≠i ˜ –D ˜ ,
---------- = – ω i ω k ---------------------------------------- , ( B ) = D kj ij i≠j (3-44)
D ik ˜ i kj
D ij ( adjB i ) jk
j≠i
where (adjBi)jk is the jkth component of the adjoint of the matrix Bi.
Solving for Equation 3-44 leads to a number of algebraic expressions for each of the
components in the multicomponent Fick diffusivity matrix. For two- and
three-component systems, these are implemented and solved directly by COMSOL
˜
Multiphysics. For instance, the component D 12 in a ternary system is given by:
ω1 ( ω2 + ω3 ) ω2 ( ω1 + ω3 ) ω 32
-------------------------------- + -------------------------------- – --------------- -
˜ x 1 D 23 x 2 D 13 x 3 D 12
D 12 = ----------------------------------------------------------------------------------------------
x1 x2 x3
-------------------
- + ------------------- - + ------------------- -
D 12 D 13 D 12 D 23 D 23 D 13
˜ = N –g
D (3-45)
ij ij
˜
where ij are indices in the matrices D and N, and ranges from 1 to the number of
species, Q.
–1
N ij = ( P ) ij (3-46)
ωi ωj ˜
P ij = ------------ – C ij
g
˜
The matrix C in turn is defined as
xi xj
--------- i≠j
˜ D ij
C ij =
– C˜ ik i=j
k≠j
The term g in Equation 3-45 is a scalar value that provides numerical stability and
should be of the same order of magnitude as the multicomponent Maxwell-Stefan
diffusion coefficients. The physics interface therefore defines q as the sum of the
multicomponent Maxwell-Stefan diffusion coefficients:
n–1
n
g =
D ij
i = 1 j = i+1
This definition for g works well in most cases. In rare cases, it might be
necessary to change the value to obtain convergence.
Assuming that the diffusive flux, relative to the mass averaged velocity, is proportional
to the mole fraction gradient, the mass flux is defined as:
m ∇x i
j md, i = – ρ i D i --------- (3-47)
xi
Here ρi is the density, and xi the mole fraction of species i. Using the definition of the
species density and mole fraction
ωi
ρ i = ρω i , x i = ------- M
Mi
Equation 3-47 can be expressed in terms of the mass fractions (ωi) in the manner of
m ∇M
j md, i = – ρD i ∇ω i + ρω i D i ---------
m
M
Using Equation 3-47 together with the Maxwell-Stefan equations, where isobaric and
isothermal conditions have been assumed, the following expression for the
mixture-averaged diffusion coefficients can be derived (Ref. 3):
m 1 – ωi
D i = -------------------------- (3-48)
xk
N
---------
k≠iD
ik
If instead the diffusive flux (relative to the mass averaged velocity) is assumed
proportional to the mass fraction gradient, the mass flux is defined as:
m* ∇ω i m*
j md, i = – ρ i D i ---------- = – ρD i ∇ω i (3-49)
ωi
1 xk xi ωk
k ≠ i --------
D ik 1 – ω i k ≠ i D ik
N N
----------- = - + --------------- --------- (3-50)
m*
Di
i = 1 ρωi ( ud,i + uc ) =
N
j md, i = 0 (3-51)
i = 1 ωi ud,i
N
uc = – (3-52)
Here ud,i is the diffusion velocity resulting from the flux assumption in Equation 3-47
or Equation 3-49. Note that the correction velocity is a constant correction (same for
all species), but varies in space.
Using the correction velocity together with Equation 3-47, the resulting diffusive flux
is
m ∇M Mi
i = 1 ------
m N m
j md, i = – ρD i ∇ω i – ρω i D i --------- + ρω i - D ∇x i (3-53)
M M i
i = 1 Di
m* N m*
j md, i = – ρD i ∇ω i + ρω i ∇ω i (3-54)
T ∇T
j i = ∇ ⋅ D i -------- + ρω i z i u m, i F ∇φ
T
where
T
• D i (SI unit: kg/(m·s)) is the thermal diffusion coefficient
• zi (dimensionless) is the charge number
• um,i the mobility of the ith species, and
• φ (SI unit: V) is the electric potential.
F ∇x i
j md, i = – ρ i D i, kl --------- (3-55)
xi
when assuming that the diffusive flux is proportional to the mole fraction gradient. If
instead assuming that it is proportional to the mass fraction it becomes
F
j md, i = – ρD i, kl ∇ω i (3-56)
F
In the equations above D i, kl represents a general diffusion matrix (SI unit: m2/s)
describing the diffusion of species i into the mixture. This form makes it possible to
use any diffusion coefficient, matrix, or empirical model based on Fick’s law. For
example, in situations when the mass transport is not dominated by diffusion, an
alternative is to use the diffusion coefficients at infinite dilution,
F 0
D i, kk = D i
The mixture diffusion correction described above for the mixture-averaged diffusion
can also be applied in this case. Correspondingly, the resulting diffusive flux is
F ∇M Mi
i = 1 ------
F N F
j md, i = – ρD i, kl ∇ω i – ρω i D i, kl --------- + ρω i - D ∇x (3-57)
M M i, kl i
i = 1 Di, kl ∇ωi
F N F
j md, i = – ρD i, kl ∇ω i + ρω i (3-58)
When using the Fick’s Law approximation, Additional Transport Mechanisms can be
accounted for in the same manner as described above for the mixture-averaged
approximation.
Di
T
= 0
i=1
c c c c
1 Ri – Ri 1 Ri + Ri
R i = --- -------------------------------- max ( ω i, 0 ) + --- -----------------------------------------
- max ( 1 – ω i, 0 ) (3-59)
2 max ( ω , ω ) dl 2 max ( 1 – ω , ω dl )
i i i i
The first term on the right hand side of Equation 3-59 is active if Ric < 0, that is if ωi
is a reactant. The reaction rate contribution, Ri, is equal to the “core” reaction rate,
Ric, as long as ωi > ωidl. As ωi approaches zero, the regularization damps out negative
Ric and for ωi < 0, Ri for reactant ωi is equal to zero.
The second term on the right hand side of Equation 3-59 is active if Ric > 0, that is if
ωi is a product. The reaction rate contribution, Ri, is equal to the “core” reaction rate,
Ric, as long as ωi < 1−ωidl. As ωi approaches one, the regularization damps out
positive Ric and for ωi > 1, Ri for product ωi is equal to zero.
The damping limits, ωidl, should be in an order of magnitude that can be considered
numerical noise for species i. The damping limits are per default set to 1e−6, which is
appropriate for most applications. It can be advantageous to lower some limits when
working with for example catalytic trace species and the limits can sometimes be raised
to gain additional robustness.
2. R.B. Bird, W.E. Stewart, and E.N. Lightfoot, Transport Phenomena, 2nd ed., John
Wiley & Sons, 2005.
3. R.J. Kee, M.E. Coltrin, and P. Glarborg, Chemically Reacting Flow, John Wiley &
Sons, 2003.
The species concentrations are denoted, ci (SI unit: mol/m3), and the potential, φ l
(SI unit: V).
N i = – D i ∇c i – z i u m, i Fc i ∇φ l + uc i = J i + uc i (3-60)
where Di (m2/s) is the diffusion coefficient, zi (1) the corresponding charge, um,i (SI
unit: s·mol/kg) the mobility and u (SI unit: m/s) the velocity vector. Ji denotes the
molar flux relative to the convective transport (SI unit: mol/(m2·s)). For a detailed
description of the theory of these equations and the different boundary conditions, see
Theory for the Transport of Diluted Species Interface.
il = zi Ni (3-61)
i
∇ ⋅ il = Ql (3-62)
where Ql (A/m3) is the electrolyte current source stemming from, for example,
porous electrode reactions. For non-porous electrode domains this source term is
usually zero.
zi ci = 0 (3-63)
i=1
In water-based systems the species H+ and OH- are always present. The auto
ionization reaction for water is
+ -
H 2 O ↔ H + OH (3-64)
2 –6
c H + c OH - = K w × 1 mol dm (3-65)
– 14
where K w ≈ 10 .
Now, the electroneutrality condition, including the two additional species H+ and
OH-, reads
c H + – c OH - + zi ci = 0 (3-66)
i=1
Combining these two equations results in the following algebraic expressions for the
concentrations of H+ and OH-.
2
Σ Σ 2 –6
c H + = – --- + ------ + K w × 1 mol dm (3-67)
2 4
and
2
Σ Σ 2 –6
c OH - = --- + ------ + K w × 1 mol dm (3-68)
2 4
where
Σ = zi ci (3-69)
i=1
to
+ z
z +1 z + [ H ] [ S 00 ]
S k0 ↔ S 00 + H K a, k = -------------------------
z +1
- (3-71)
[ S 10 ]
where z0 is the charge (valence) of species S0 (which has no dissociable protons) and
Ka,j is the acid (equilibrium) constant of the jth dissociation reaction. The brackets [ ]
here represents the species activity. The charge of each species is always deductible from
the index i according to z0+i and will be dropped from now on.
If the proton activity is known, any species Sm may be expressed using any other species
Sl according to
m–l
[H] [ Sl ]
[ S m ] = -------------------------------------
k–l
(3-72)
∏ K a, j
j = k–m+1
if m>l and
k–m
m–l
[ Sm ] = [ H ] [ Sl ] ∏ K a, j
j = k–l+1
(3-73)
if l>m.
Setting m=i and denoting the flux of species i by Ni using equation Equation 3-60,
the mass balance equation for the concentration ci of each subspecies i in the
dissociation chain is
δc
-------i + ∇ ⋅ N i = R eq, i, k – i – R eq, i, k – i + 1 + R i (3-74)
δt
where Req,i,j is the reaction source stemming from the jth dissociation step (with
Req,i,k+1=0), and Ri any additional reaction sources.
The reaction source contributions from the dissociation steps are generally not known,
but may be canceled by taking the sum of all mass balance equations, resulting in
N i = – D∇c i – zu m Fc i ∇φ l + uc i (3-76)
When considering the contribution to the current and the charge balance equation one
2 2
needs to take into account that the squared average charge, z = ( z 0 – ν ) , is not
2 2 2 2 2 2
equal to the “average squared charge”, z = z 0 + 2z 0 ν + ν = z – ν + ν (Ref. 1).
2
i l = … – F ( zD i ∇c i + z u m c i ∇φ l ) (3-77)
The average number of protons removed from the proton typically depends on the
pH. If the average number of removed protons depend only on the pH, the averaged
squared number of protons removed can be written as
2 d ν ν2
ν = – c H+ + (3-78)
d cH+
And from this one can derive the average squared charge according to
2 d ν z2 d z z2
z = – cH+ + = – cH+ + (3-79)
d c H+ d cH+
DIFFUSIVITY-MOBILITY RELATIONS
The Stokes radius r of a molecule is related to the diffusivity according to
kT
r = --------------- (3-80)
6πμD
For small molecules, one frequently uses the Nernst-Einstein relation between the
diffusivity and the mobility
D
u m = -------- (3-81)
RT
For larger molecules, such as proteins, the mobility may instead be calculated based on
the Debye-Hückel-Henry expression (Ref. 2 ) according to
ef ( κr ) Df ( κr )
u m = ------------------------------------ = ------------------------------ (3-82)
6πμF ( 1 + κr ) RT ( 1 + κr )
where κ (1/m) is the Debye parameter, which depends on the ionic strength of the
solution, is defined for ideal solutions as
2 N
2e N A
zi ci
2 2
κ = ----------------- (3-83)
εε 0 kT
i=1
where ε is the dielectric constant of the fluid and ε0 the permittivity of free space.
2
( z should be used if available in the formula above when calculating the ionic
strength).
The function f above is based on a sigmoidal function so that it ranges from 1 for
κr = 0 to 1.5 for κr = ∞ . Note that the Debye-Hückel-Henry expression
approaches the Nernst-Einstein mobility as r → 0 .
For larger molecules (macro ions), where the distance between the charges is large
compared to 1/κ, the Linderstrøm-Lang approximation postulates a smaller
contribution of to the ionic strength so that the z-valent ion behaves as a monovalent
ion with a z-fold concentration. For an assemble of N-M smaller molecules and M
macro ions, the Debye parameter then is defined as
N–M N
2e N A
2
κ = ----------------- abs(z i)c i
2 2
zi ci + (3-84)
εε 0 kT
i=1 i = N–M+1
REFERENCES
1 The Dynamics of Electrophoresis, Mosher, Saville and Thormann, VCH
Verlagsgesellschaft mbH, Weinheim, Germany, 1992.
In this section:
N t, i = – D s , i ∇ t c s , i
where Ds,i (SI unit: m2/s) is the surface diffusion coefficient for species i.
∂c s, i
= – ∇ t ⋅ N t, i + R s , i (3-85)
∂t
where Rs,i (SI unit: mol/(m2·s)) is the sum of all sources due to surface reactions and
adsorption/desorption phenomena.
Of frequent interest for surface reaction kinetics are the fractional surface coverages, θi
(dimensionless), of the species (with index i).
σ i c s, i
θ i = --------------
Γs
(The site occupancy number accounts for the situation when a large species covers
more than one site on the surface.)
For the case of monolayer adsorption, the sum of all fractional coverages of free and
adsorbed sites is unity, and hence the fraction of free sites on the surface, θ*, can be
calculated from:
θ* = 1 – θi
i
The reaction rate in mass basis, rb,k (SI unit: kg/(s·m2)) for species k, is given by:
r k = M k R b, k
with Mk (SI unit: kg/mol) being the molar mass of the species.
Based on this, the species contribution to the bulk growth velocity, vk (SI unit: m/s),
is given by:
r b, k
v b, k = ----------
ρk
r b, tot = r b, k
k
v b, tot = vb, k
k
R b, k
R b, frac, k = ---------------
R b, tot
r b, k
r b, frac, k = -------------
r b, tot
v b, k
v b, frac, k = --------------
v b, tot
The bulk concentration, cb,k (SI unit: mol/m2), for species k is governed by the
equation:
∂c b, k
= R b, k (3-86)
∂t
The bulk concentration in mass basis, mb,k (SI unit: kg/m2) for a species k, can be
derived from:
m b, k = M k c b, k
m b, k
s b, k = -------------
ρk
c b, tot = c b, k
k
m b, tot = mb, k
k
s b, tot = sk
k
FRACTIONAL QUANTITIES
The fractional bulk concentration, bulk mass, and thickness (all dimensionless) are
calculated according to:
c b, k
c b, frac, k = -------------
c b, tot
m b, k
m b, frac, k = ----------------
m b, tot
s b, k
s b, frac, k = -------------
s b, tot
In order to comply with the additional contributions to the mass balance, equations
are added. First, the following terms are added to the right-hand side of Equation 3-85
and Equation 3-86, respectively.
– c s, i ∂ ln ∂ A
∂t
– c b, k ∂ ln ∂ A
∂t
where ∂A is the infinitesimal mesh area segment (area scale factor). The above terms
account for the concentration change due to a fractional area change.
∇ t ⋅ ( c s, i v t, mesh )
∇ t ⋅ ( c b, k v t, mesh )
This convectional term needs often to be stabilized using methods such as streamline
diffusion or isotropic diffusion.
∂c i
+ ∇ ⋅ ( c i u ) = – ∇ ⋅ ( – D i ∇c i – u m, i F c i ∇V ) + R i = – ∇ ⋅ J i + R i (3-87)
∂t
where Ji is the molar flux relative to the convective transport, and Ri (SI unit: mol/
(m3 ·s)) is the reaction term. The velocity, u, is equal to the velocity of the solvent.
This implies that the solute’s contribution to the solvent’s velocity, through shear or
any other forces, is negligible in comparison to the solvent’s contribution to the solute.
Equation 3-87 introduces one variable for the concentration of each of the dissolved
species and the electric potential, V.
Apart from the transport equations, the physics interface assumes that the
electroneutrality condition holds:
zi ci = 0 (3-88)
that is that the net charge in every control volume is zero. This means that Q − 1,
where Q is the number of species present, can be solved for using Equation 3-87. The
remaining species concentration is computed from the electroneutrality condition.
This means that boundary conditions for this species cannot be specified, although it
takes part in the boundary condition descriptions for the current density. Often, the
species chosen to be computed from electroneutrality is the oppositely charged ion, to
the electroactive species, from a supporting electrolyte.
n
∂c i n n
F zi
∂t
+ ∇ ⋅ F
z N
i i = F zi Ri (3-89)
i=1 i=1 i=1
The first term in Equation 3-89 is zero, which can be shown by taking the time
derivative of the electroneutrality condition. The expression under the divergence
operator is the total current density vector, defined by:
i = F z i ( – D i ∇c i – z i u m , i F c i ∇V ) (3-90)
i=1
Notice that no convective term is included in the expression for the current density,
which is also a result of the electroneutrality condition.
∇⋅i = F zi Ri (3-91)
i=1
This equation states the conservation of electric charge and is the one solved for to
compute the electric potential.
Equation 3-87, Equation 3-88, and Equation 3-91 are sufficient for describing the
potential and concentration distribution in an electrochemical cell or in an electrolyte
subjected to an electric field.
A useful observation from Equation 3-90 is that the ionic conductivity, defined in
absence of concentration gradients, is implicitly given by:
κ = F
2
z2i ci um, i
i=1
( –Di ∇ci )
i=1
∇V = ------------------------------------
n
zi um, i Fci
i=1
∇⋅u = 0 (3-92)
∂c i
+ ∇ ⋅ ( – D i ∇c i – z i u m, i Fc i ∇V ) + u ⋅ ∇c i = R i (3-93)
∂t
which can be called the nonconservative formulation and is the default formulation in
the Nernst-Planck Equations interface.
ω i – nojac ( ω i )
ρ ------------------------------------
Δt̃
is added to the left-hand side of the mass fraction equations. Here ρ is the fluid mixture
density, ωi is the mass fraction (dimensionless) of species i, and Δt̃ is the pseudo time
step. Since ωi−nojac(ωi) is always zero, this term does not affect the final solution. It
does, however, affect the discrete equation system and effectively transforms a
nonlinear iteration into a time step of size Δt̃ .
For a description of the pseudo time step term for the Navier-Stokes
equations and the pseudo time step see Pseudo Time Stepping for
Laminar Flow Models and Pseudo Time Stepping in the COMSOL
Multiphysics Reference Manual.
Here, n is the unit normal pointing out of the fluid domain, u is the mass averaged
velocity of the fluid mixture (SI unit: m/s), ji denotes the mass flux of species i relative
to relative to the mixture (typically due to diffusion), and Mi is the species molar mass
(SI unit: kg/mol). Summing the mass balances at the surface, over all species, results
in an effective mixture velocity:
n ⋅ ρu s = rs,i Mi (3-95)
i=1
referred to as the Stefan velocity, here denoted us. To reach Equation 3-95 the fact
that the sum of all mass fractions is one, and that the sum of all relative diffusive fluxes
is zero, was used.
Equation 3-95 implies that surface reactions result in a net flux between the surface
and the domain. A net flux in turn corresponds to an effective convective velocity at
the domain boundary; the Stefan velocity. It should be noted here that when solving
for mass transport inside a fluid domain, an outer boundary of the domain corresponds
to a position just outside of the actual physical wall (on the fluid side). The domain
boundary does not coincide with the physical wall.
In most reacting flow models, the species mass fractions in the fluid domain are solved
for without including the surface concentrations (mol per area) on exterior walls. One
reason for this is that the surface reaction rates are often not known. In this case,
surface reactions can be modeled either by applying a mass flux or prescribing the mass
fraction, or a combination of both, on fluid boundaries adjacent to the reacting
surface. The Stefan velocity on a fluid domain boundary is then defined as the net mass
flux resulting from the boundary conditions applied:
Here, the first term contains contributions from boundary conditions prescribing the
mass flux, while the second contains contributions from boundary conditions
prescribing the mass fractions. Contributions to the Stefan velocity can be added by
selecting Account for Stefan velocity in the Flux or Mass Fraction features in The
Transport of Concentrated Species interface.
n ⋅ j 0, i
i=1
n ⋅ u s = ------------------------------------
- (3-97)
N
ρ 1 – ω i, 0
i
Using a Reacting Flow interface, the Stefan velocity, defined in the manner of
Equation 3-97, is automatically computed and applied on boundaries corresponding
to walls in the coupled fluid flow interface. The Stefan velocity is prescribed in the wall
normal direction on the wall selection.
f (3-98)
kj
aA + bB + ... r
xX + yY + ...
kj
For such a reaction set, the reaction rates rj (SI unit: mol/(m3·s)), can be described by
the mass action law:
f –ν r ν
rj = kj ∏ c i ij – k j ∏ c i ij (3-99)
i ∈ react i ∈ prod
f r
Here, kj and kj denote the forward and reverse rate constants, respectively. The
concentration of species i is denoted ci (SI unit: mol/m3). The stoichiometric
coefficients are denoted νij, and are defined to be negative for reactants and positive
for products. In practice, a reaction seldom involves more than two species colliding
in a reacting step, which means that a kinetic expression is usually of order 2 or less
(with respect to the involved concentrations).
Here, A denotes the frequency factor, n the temperature exponent, E the activation
energy (SI unit: J/mol) and Rg the gas constant, 8.314 J/(mol·K). The
pre-exponential factor, including the frequency factor A and the temperature factor
Tn, is given the units (m3/mol)α − 1/s, where α is the order of the reaction (with
respect to the concentrations).
The interface supports simulation of transport by convection and diffusion in 1D, 2D,
and 3D as well as for axisymmetric components in 1D and 2D. The dependent variable
is the molar concentration, c. Modeling multiple species transport is possible, whereby
the physics interface solves for the molar concentration, ci, of each species i.
SETTINGS
The Label is the default physics interface name.
The Name is used primarily as a scope prefix for variables defined by the physics
interface. Refer to such physics interface variables in expressions using the pattern
<name>.<variable_name>. In order to distinguish between variables belonging to
different physics interfaces, the name string must be unique. Only letters, numbers, and
underscores (_) are permitted in the Name field. The first character must be a letter.
The default Name (for the first physics interface in the model) is tds.
DOMAIN SELECTION
If any part of the model geometry should not partake in the mass transfer model,
remove that part from the selection list.
TRANSPORT MECHANISMS
Mass transport due to diffusion is always included. Use the check boxes available under
Additional transport mechanisms to control other transport mechanisms.
• By default, the Convection check box is selected. Clear the check box to disable
convective transport.
• Select the Migration in electric field check box to activate transport of ionic species in
an electric field. See further the theory section Adding Transport Through
Migration.
• Select the Dispersion in porous media check box to activate the dispersion mechanism
in porous media. See further Dispersion in the theory chapter.
• Select the Volatilization in partially saturated porous media check box to model
volatilization in partially saturated domains. See further Theory for the Transport of
Diluted Species Interface.
The following features are also enabled when selecting the Mass transport in porous
media check box:
• Adsorption
• Fracture
• Partially Saturated Porous Media
• Porous Media Transport Properties
• Reactive Pellet Bed
• When the Crosswind diffusion check box is selected, a weak term that reduces
spurious oscillations is added to the transport equation. The resulting equation
system is always nonlinear. There are two options for the Crosswind diffusion type:
- Do Carmo and Galeão — the default option. This type of crosswind diffusion
reduces undershoots and overshoots to a minimum but can in rare cases give
equation systems that are difficult to fully converge.
- Codina. This option is less diffusive compared to the Do Carmo and Galeão
option but can result in more undershoots and overshoots. It is also less effective
for anisotropic meshes. The Codina option activates a text field for the Lower
gradient limit glim. It defaults to 0.1[mol/m^3)/tds.helem, where tds.helem
is the local element size.
• For both consistent stabilization methods, select an Equation residual. Approximate
residual is the default and means that derivatives of the diffusion tensor components
are neglected. This setting is usually accurate enough and is computationally faster.
If required, select Full residual instead.
INCONSISTENT STABILIZATION
To display this section, click the Show button ( ) and select Stabilization. By default,
the Isotropic diffusion check box is not selected, because this type of stabilization adds
artificial diffusion and affects the accuracy of the original problem. However, this
option can be used to get a good initial guess for underresolved problems.
ADVANCED SETTINGS
To display this section, click the Show button ( ) and select Advanced Physics Options.
Normally these settings do not need to be changed. Select a Convective term—
Nonconservative form (the default) or Conservative form. The conservative formulation
should be used for compressible flow. See Convective Term Formulation for more
information.
DISCRETIZATION
To display all settings available in this section, click the Show button ( ) and select
Advanced Physics Options.
The Compute boundary fluxes check box is activated by default so that COMSOL
Multiphysics computes predefined accurate boundary flux variables. When this option
If the check box is cleared, the COMSOL Multiphysics software instead computes the
flux variables from the dependent variables using extrapolation, which is less accurate
in postprocessing results but does not create extra dependent variables on the
boundaries for the fluxes.
• ndflux_c (where c is the dependent variable for the concentration). This is the
normal diffusive flux and corresponds to the boundary flux when diffusion is the
only contribution to the flux term.
• ntflux_c (where c is the dependent variable for the concentration). This is the
normal total flux and corresponds to the boundary flux plus additional transport
terms, for example, the convective flux when you use the nonconservative form.
Also the Apply smoothing to boundary fluxes check box is available if the previous check
box is checked. The smoothing can provide a more well-behaved flux value close to
singularities.
For details about the boundary fluxes settings, see Computing Accurate Fluxes in the
COMSOL Multiphysics Reference Manual.
The Value type when using splitting of complex variables setting should in most pure
mass transfer problems be set to Real, which is the default. It makes sure that the
dependent variable does not get affected by small imaginary contributions, which can
occur, for example, when combining a Time Dependent or Stationary study with a
frequency-domain study. For more information, see Splitting Complex-Valued
Variables in the COMSOL Multiphysics Reference Manual.
DEPENDENT VARIABLES
The dependent variable name is the Concentration c by default. The names must be
unique with respect to all other dependent variables in the component.
Add or remove species variables in the model and also change the names of the
dependent variables that represent the species concentrations.
Enter the Number of species. Use the Add concentration ( ) and Remove
concentration ( ) buttons as needed.
It applies to one or more diluted species or solutes that move primarily within a fluid
that fills (saturated) or partially fills (unsaturated) the voids in a solid porous medium.
The pore space not filled with fluid contains an immobile gas phase. Models including
a combination of porous media types can be studied.
The main feature nodes are the Porous Media Transport Properties and Partially
Saturated Porous Media nodes, which add the equations for the species concentrations
The physics interface can be used for stationary and time-dependent analysis.
When this physics interface is added, these default nodes are also added to the Model
Builder — Porous Media Transport Properties, No Flux (the default boundary condition),
and Initial Values. Then, from the Physics toolbar, add other nodes that implement, for
example, boundary conditions, reaction rate expressions, and species sources. You can
also right-click Transport of Diluted Species in Porous Media to select physics features
from the context menu.
SETTINGS
The rest of the settings are the same as The Transport of Diluted Species Interface.
FURTHER READING
Transport Properties
The settings in this node are dependent on the check boxes selected under Transport
Mechanisms on the Settings window for the Transport of Diluted Species interface. It
includes only the sections required by the activated transport mechanisms. It has all the
equations defining transport of diluted species as well as inputs for the material
properties.
When the Convection check box is selected, the Turbulent Mixing subnode is available
from the context menu as well as from the Physics toolbar, Attributes menu.
MODEL INPUTS
The temperature model input is always available. Select the source of the Temperature.
For User defined, enter a value or expression for the temperature in the input field. This
input option is always available.
You can also select the temperature solved for by a Heat Transfer interface added to
the model component. These physics interfaces are available for selection in the
Temperature list.
CONVECTION
If transport by convection is active, the velocity field of the solvent needs to be
specified. Select the source of the Velocity field. For User defined, enter values or
expressions for the velocity components in the input fields. This input option is always
available.
You can also select the velocity field solved for by a Fluid Flow interface added to the
model component. These physics interfaces are available for selection in the Velocity
field list.
DIFFUSION
Select an option from the Material list. This selection list can only be used if a material
has been added in the Materials node and if that material has a diffusion coefficient
defined. Else, you need to type in the diffusivity in the User Defined edit field.
Note that multiple species, as well as Migration in Electric fields (described below) is
only available for certain COMSOL Multiphysics add-on products. See details: http:/
/www.comsol.com/products/specifications/.
• Enter a value or expression for the Electric potential V, which is User defined; this
input option is always available.
• Select the electric potential solved by an AC/DC-based interface that has also been
added to the model.
• Select the electric potential defined or solved by Electrochemistry interface that has
been added to the component.
By default the Mobility is set to be calculated based on the species diffusivity and the
temperature using the Nernst-Einstein relation. For User defined, and under Mobility,
select the appropriate scalar or tensor type — Isotropic, Diagonal, Symmetric, or
Anisotropic — and type in the value of expression of the mobility um,c.
Enter the Charge number zc (dimensionless, but requires a plus or minus sign) for each
species.
The temperature (if you are using mobilities based on the Nernst-Einstein relation) is
taken from Model Inputs section.
Note that the migration in electric fields feature is only available in some COMSOL
products. See details: http://www.comsol.com/products/specifications/.
Turbulent Mixing
Use this node to account for the turbulent mixing of the chemical species caused by
the eddy diffusivity. This node should typically be used when the specified velocity field
corresponds to a RANS solution.
The subnode can added from the context menu (right-click the Transport Properties
parent node), as well as from the Physics toolbar, Attributes menu, provided that
Convection is selected as a transport mechanism.
FURTHER READING
See the section About Turbulent Mixing in the CFD Module User’s Guide (this link
is available online or if you have the CFD Module documentation installed).
DOMAIN SELECTION
If there are several types of domains with different initial values defined, it might be
necessary to remove some domains from the selection. These are then defined in an
additional Initial Values node.
INITIAL VALUES
Enter a value or expression for the initial value of the Concentration or concentrations,
ci. This also serves as a starting guess for stationary problems.
Mass-Based Concentrations
Use the Mass-Based Concentrations node to add postprocessing variables for mass-based
concentrations (SI unit: kg/m3) and mass fractions (dimensionless) for all species.
MIXTURE PROPERTIES
The default Solvent density ρsolvent is taken From material. For User defined, enter a
value or expression manually. Define the Molar mass of each species, which is needed
to calculate the mass-based concentration.
Reactions
Use the Reactions node to account for the consumption or production of species
through chemical reactions. Define the rate expressions as required.
DOMAIN SELECTION
From the Selection list, choose the domains on which to define rate expression or
expressions that govern the source term in the transport equations.
Several reaction nodes can be used to account for different reactions in different parts
for the modeling geometry.
REACTION RATES
Add a rate expression R (SI unit: mol/(m3·s)) for species i. Enter a value or expression
in the field. Note that if you have the Chemistry interface available, provided with the
REACTING VOLUME
This section is only available when the Mass Transport in Porous Media property is
available and selected. See http://www.comsol.com/products/specifications/ for
more details on availability.
When specifying reaction rates for a species in porous media, the specified reaction rate
may have the basis of the total volume, the pore volume, or the volume of a particular
phase.
• For Total volume, the reaction expressions in mol/(m3·s) are specified per unit
volume of the model domain (multiplied by unity).
• For Pore volume, the reaction expressions in mol/(m3·s) are specified per unit
volume of total pore space. The reaction expressions will be multiplied by the
domain porosity, εp. (εp equals unity for nonporous domains.)
• For Liquid phase, the reaction expressions in mol/(m3·s) are specified per unit
volume of liquid in the pore space. The expressions will be multiplied by the liquid
volume fraction θ. (θ equals εp for Saturated Porous Media domains).
• For Gas phase, the expressions are multiplied by the gas volume fraction av = εp − θ.
av equals 0 for Saturated Porous Media domains.
FURTHER READING
See the theory chapter on chemical species transport, starting with the section Mass
Balance Equation.
–n ⋅ Ji = 0
Inflow
Use this node to specify all species concentrations at an inlet boundary.
If you want to specify the concentration of a subset of the partaking species, this can
be done by using the Concentration node instead.
For the Electroanalysis interface, this node is available when you select the Convection
check box on the physics interface Settings window.
CONCENTRATION
For the concentration of each species c0,c (SI unit: mol/m3), enter a value or
expression.
CONSTRAINT SETTINGS
To display this section, click the Show button ( ) and select Advanced Physics Options.
You can find details about the different constraint settings in the section Constraint
Reaction Terms in the COMSOL Multiphysics Reference Manual.
FURTHER READING
See the theory chapter in the section Danckwerts Inflow Boundary Condition.
n ⋅ ( – D ∇c ) = 0
Note that the Convection or the Migration in electric field transport mechanisms needs
to be included for this node to be available.
Concentration
This condition node adds a boundary condition for the species concentration. For
example, a c = c0 condition specifies the concentration of species c.
CONCENTRATION
Individually specify the concentration for each species. Select the check box for the
Species to specify the concentration, and then enter a value or expression in the
corresponding field. To use another boundary condition for a specific species, click to
clear the check box for the concentration of that species.
CONSTRAINT SETTINGS
To display this section, click the Show button ( ) and select Advanced Physics Options.
You can find details about the different constraint settings in the section Constraint
Reaction Terms in the COMSOL Multiphysics Reference Manual.
Flux
This node can be used to specify the species flux across a boundary. The prescribed flux
of a species c is defined as:
– n ⋅ ( – D∇c ) = J 0
When the mass transport includes migration of ionic species, the flux is defined as:
– n ⋅ ( – D∇c – zu m Fc∇φ ) = J 0
INWARD FLUX
Select the Species check box for the species to specify and enter a value or expression
for the inward flux in the corresponding field. Use a minus sign when specifying a flux
directed out of the system. To use another boundary condition for a specific species,
click to clear the check box for that species.
External convection
Set Flux type to External convection to prescribe a flux to or from an exterior domain
(not modeled) assumed to include convection. The exterior can for example include a
forced convection to control the temperature or to increase the mass transport. In this
case the prescribed mass flux corresponds to
J0 = kc ( cb – c )
where kc is a mass transfer coefficient and cb is the bulk concentration, the typical
concentration far into the surrounding exterior domain.
Symmetry
The Symmetry node can be used to represent boundaries where the species
concentration is symmetric, that is, where there is no mass flux across the boundary.
Flux Discontinuity
This node represents a discontinuity in the mass flux across an interior boundary:
– n ⋅ [ ( J + uc ) u – ( J + uc ) d ] = N 0 J = – D∇c – zu m Fc∇φ
where the value N0 (SI unit: mol/(m2·s)) specifies the jump in total flux at the
boundary. This can be used to model a boundary source, for example a surface
reaction, adsorption or desorption.
Select the Species check box for the species to specify and enter a value or expression
for the material flux jump in the corresponding field. To use a different boundary
condition for a specific species, click to clear the check box for the flux discontinuity
of that species.
Partition Condition
The Partition Condition node can be used to prescribe the ratio between the
concentration of a solute species in two different immiscible phases. It can for example
be used on interior boundaries separating two liquid phases, a gas-liquid interface, or
on a boundary separating a liquid phase and a solid or porous media. For a species
concentration ci, the ratio between the concentration on the up side and on the down
side of the boundary (ci,u and ci,d respectively) is defined as
c i, u
K i = ---------
c i, d
in terms of a partition coefficient Ki. The up and down side of the selected boundary
is indicated in the Graphics window. The arrows point from the down side into the up
side.
PARTITION COEFFICIENT
Select the Reverse direction check box to reverse the direction of the arrows on the
selected boundaries, and the corresponding definition of the up and down side
concentration.
Use the associated input fields to prescribe the partition coefficient Ki.
Periodic Condition
The Periodic Condition node can be used to define periodicity for the mass transport
between two sets of boundaries. The node prescribes continuity in the concentration
and the mass flux between the “source” and the “destination” side respectively. Note
that these names are arbitrary and does not influence the direction in which mass is
transported. It is dictated by mass transfer equations in the adjacent domains.
The node can be activated on more than two boundaries, in which case the feature tries
to identify two separate surfaces that each consist of one or several connected
boundaries.
For more complex geometries, it might be necessary to add the Destination Selection
subnode, which is available from the context menu (right-click the parent node) as well
as from the Physics toolbar, Attributes menu. With this subnode, the boundaries that
constitute the source and destination surfaces can be manually specified.
FURTHER READING
For an example of using a periodic condition, see this application example:
2D Points
2D Axisymmetry Points not on the symmetry axis and the symmetry axis
3D Edges
SPECIES SOURCE
·
Enter the source strength, q l,c , for each species (SI unit: mol/(m·s)). A positive value
results in species injection from the line into the computational domain, and a negative
value means that the species is removed from the computational domain.
Line sources located on a boundary affect the adjacent computational domains. This
effect makes the physical strength of a line source located in a symmetry plane twice
the given strength.
FURTHER READING
See the section Mass Sources for Species Transport.
SPECIES SOURCE
·
Enter the source strength, q p,c , for each species (SI unit: mol/s). A positive value
results in species injection from the point into the computational domain, and a
negative value means that the species is removed from the computational domain.
Open Boundary
Use this node to set up mass transport across boundaries where both convective inflow
and outflow can occur. On the parts of the boundary where fluid flows into the
domain, an exterior species concentration is prescribed. On the remaining parts, where
fluid flows out of the domain, a condition equivalent to the Outflow node is instead
prescribed.
The direction of the flow across the boundary is typically calculated by a fluid flow
interface and is provided as a model input to the Transport of Diluted Species
interface.
EXTERIOR CONCENTRATION
Enter a value or expression for the Exterior concentration.
Equilibrium Reaction
Use this node to model a reaction where the kinetics is assumed so fast that the
equilibrium condition is fulfilled at all times. The node solves for an additional degree
of freedom (the reaction rate Req) to fulfill the equilibrium condition at all times in all
space coordinates.
If the Apply equilibrium condition on inflow boundaries check box is selected, the
specified inflow concentration values in all active Inflow boundary nodes for the physics
interface are modified to comply with the equilibrium condition.
EQUILIBRIUM CONDITION
The list defaults to Equilibrium constant or select User defined. For either option, the
Apply equilibrium condition on inflow boundaries check box is selected by default.
For Equilibrium constant, enter an Equilibrium constant Keq (dimensionless). Also enter
a value or expression for the Unit activity concentration Ca0 (SI unit: mol/m3).
Selecting Equilibrium constant defines an equilibrium condition based on the
stoichiometric coefficients, the species activities, and the law of mass action.
STOICHIOMETRIC COEFFICIENTS
Enter a value for the stoichiometric coefficientνc (dimensionless). The default is 0. Use
negative values for reactants and positive values for products in the modeled reaction.
Species with a stoichiometric coefficient value of 0 are not affected by the Equilibrium
Reaction node.
Surface Reactions
The Surface Reactions node can be used to account for the species boundary flux due
to chemical reactions occurring on a surface (heterogeneous reactions). For a domain
species participating in a surface reaction, the boundary flux corresponds to the
reaction rate at the surface.
FURTHER READING
For an example of using the Surface Reactions node, see this application example:
The node will set the Rate limiting species concentration to zero at the boundary, and
balance the fluxes of the species participating in the reaction and the current densities
according to the Stoichiometric Coefficients settings.
In the Transport of Concentrated Species interface, the molar sources (or sinks) are
multiplied by the species molar masses to obtain the corresponding mass sources.
Additional Reaction Coefficients subnodes are available from the context menu
(right-click the parent node) as well as from the Physics toolbar, Attributes menu.
Note that if you are also modeling the momentum transport and expect a
non-negligible total mass source or sink, which is often the case in gas diffusion
electrodes, you need to also add a corresponding Porous Electrode Coupling node in
the Fluid Flow interface.
The molar flux or source is proportional to the stoichiometric coefficients and the
current density according to Faraday’s law.
Current densities from Electrode Reaction (iloc, SI unit: A/m2) or Porous Electrode
Reaction nodes (iv, SI unit: A/m3) of any Electrochemistry interface in the model are
available for selection as the Coupled reaction, and user-defined expressions are also
supported.
The flux is proportional to the current densities and the stoichiometric coefficients
according to Faraday’s law as defined by summation over the Reaction Coefficients
subnodes.
Note that if you are also modeling the momentum transport and expect a
non-negligible total mass flux over the boundary, which is often the case for gas
diffusion electrodes, you need to also add a corresponding Electrode Surface Coupling
node in the Fluid Flow interface.
MODEL INPUTS
The temperature model input is always available. Select the source of the Temperature.
For User defined, enter a value or expression for the temperature in the input field. This
input option is always available.
You can also select the temperature solved for by a Heat Transfer interface added to
the model component. These physics interfaces are available for selection in the
Temperature list.
MATRIX PROPERTIES
Use the Porous material list to define a material specifying the matrix properties on the
current selection. By default the Domain material is used.
Specify the Porosity, εp (dimensionless) of the porous matrix. This is by default taken
From material. Select User defined to instead enter a different value.
CONVECTION
If transport by convection is active, the velocity field of the solvent needs to be
specified. Select the source of the Velocity field. For User defined, enter values or
expressions for the velocity components in the input fields. This input option is always
available.
You can also select the velocity field solved for by a Fluid Flow interface added to the
model component. These physics interfaces are available for selection in the Velocity
field list.
DIFFUSION
Select a Fluid material (when available and applicable).
Select the Effective diffusivity model: Millington and Quirk model (the default),
Bruggeman model, Tortuosity model, or User defined. For Tortuosity model, enter a value
for the tortuosity τF,i (dimensionless). The default is 1.
• Enter a value or expression for the Electric potential V, which is User defined; this
input option is always available.
• Select the electric potential solved by an AC/DC-based interface that has also been
added to the model.
• Select the electric potential defined or solved by Electrochemistry interface that has
been added to the component.
By default the Mobility is set to be calculated based on the species effective diffusivity
and the temperature using the Nernst-Einstein relation. For User defined, and under
Mobility, select the appropriate scalar or tensor type — Isotropic, Diagonal, Symmetric,
or Anisotropic — and type in the value of expression of the effective mobility ume,c.
DISPERSION
This section is available when the Dispersion in porous media check box is selected on
the Settings window for the physics interface.
Select the Specify dispersion for each species individually check box to specify the
dispersion tensor DD (SI unit: m2/s) for each species separately. The default is to use
the same dispersion tensor DD for all species.
Select an option from the Dispersion tensor list — User defined (the default) or
Dispersivity. For User defined, use it to specify the dispersion components as
user-defined constants or expressions. Select Isotropic, Diagonal, Symmetric, or
Anisotropic based on the properties of the dispersion tensor.
Select Dispersivity when Convection has been added as the transport mechanism. Specify
the dispersivities (SI unit: m) to define the dispersion tensor DD (SI unit: m2/s)
together with the velocity field u. Select an option from the Dispersivity model list:
Isotropic (the default) or Transverse isotropic based on the properties of the porous
Adsorption
Use this node to model adsorption of the fluid phase species onto the porous media
surface. It is available as a subnode to the Porous Media Transport Properties and the
Partially Saturated Porous Media nodes.
MATRIX PROPERTIES
Use the Porous material list to define a material specifying the matrix properties on the
current selection. By default the Domain material is used.
The density of the porous media is needed when modeling adsorption to the surface
of the porous matrix. By default Density ρ is defined from the domain material.
ADSORPTION
Select a Sorption type — Langmuir (the default), Freundlich, or User defined to specify
how to compute cP, the amount of species sorbed to the solid phase (moles per unit
dry weight of the solid):
• For Langmuir:
KL c ∂c P K L c Pmax
c P = c Pmax -------------------- -------- = ---------------------------
-
1 + KL c ∂c ( 1 + KL c )
2
Enter a Langmuir constant kL,c (SI unit: m3/mol) and an Adsorption maximum
cp,max,c (SI unit: mol/kg):
• For Freundlich:
Enter a Freundlich constant kF,c (SI unit: mol/kg), a Freundlich exponent NF,c
(dimensionless), and a Reference concentration cref,c (SI unit: mol/m3).
∂c P
cP = KP c -------- = K P User defined
∂c
FURTHER READING
See the theory chapter in the section Mass Balance Equation for Transport of Diluted
Species in Porous Media.
MODEL INPUTS
The temperature model input is always available. Select the source of the Temperature.
For User defined, enter a value or expression for the temperature in the input field. This
input option is always available.
You can also select the temperature solved for by a Heat Transfer interface added to
the model component. These physics interfaces are available for selection in the
Temperature list.
Specify the Porosity, εp (dimensionless) of the porous matrix. This is by default taken
From material. Select User defined to instead enter a different value.
SATURATION
Select Saturation or Liquid volume fraction from the list.
For Saturation, enter a value for s (dimensionless) between 0 and 1. The liquid volume
fraction is then computed from the saturation and porosity as θ = sεp.
For Liquid volume fraction, enter a value for θ (dimensionless) between 0 and the value
of porosity. The saturation is then computed from the porosity and the liquid volume
fraction as s = θεp.
Select a Fluid fraction time change: Fluid fraction constant in time (the default), Time
change in fluid fraction, or Time change in pressure head.
• For Time change in fluid fraction, enter dθ/dt (SI unit: 1/s).
• For Time change in pressure head, enter dHp/dt (SI unit: m/s) and a Specific
moisture capacity Cm (SI unit: 1/m).
CONVECTION
If transport by convection is active, the velocity field of the solvent needs to be
specified. Select the source of the Velocity field. For User defined, enter values or
expressions for the velocity components in the input fields. This input option is always
available.
You can also select the velocity field solved for by a Fluid Flow interface added to the
model component. These physics interfaces are available for selection in the Velocity
field list.
DIFFUSION
Select a Liquid material from the list.
Specify the Liquid diffusion coefficient DL,c (SI unit: m2/s). Enter a value or expression
for each of the species in the corresponding input field. The default is 1·10-9 m2/s.
Select the Effective diffusivity model, liquid: Millington and Quirk model (the default),
Bruggeman model, Tortuosity model, or User defined. For Tortuosity model, enter a value
for τL,c (dimensionless). The default is 1.
DISPERSION
This section is available when the Dispersion in porous media check box is selected on
the Settings window for the physics interface. The settings are the same as for Porous
Media Transport Properties.
VOLATILIZATION
This section is available when the Volatilization in partially saturated porous media check
box is selected on the Settings window for the physics interface.
Enter a value for the Volatilization kG,c (dimensionless) for each species.
Volatilization
This feature is available when the Volatilization in partially saturated porous media check
box is selected on the Settings window for the physics interface.
Use this feature to model mass transfer at the boundary due to volatilization. The
feature can be added on boundaries of a Partially Saturated Porous Media domain. In
this case the porous media contains a liquid phase and a gas phase. The species
dissolved in the liquid are assumed to be vaporized at the boundary, and transported
into the surrounding bulk region due to convection and diffusion. The mass transfer
at the boundary is defined as
– n ⋅ J c = – h c ( k G,c c – c Gatm,c )
VOLATILIZATION
Enter a Mass transfer coefficient hc defining the transfer into the surrounding media.
This can be given by boundary layer theory. When assuming that no convective flow
is present in the surrounding, the mass transfer coefficient can be defined from the gas
diffusion coefficient DGc and the thickness of the diffusion layer ds in the manner of
D Gc
h c = ----------
ds
Also give the atmospheric concentration for each species, cGatm,c. The Volatilization
coefficient kG,c for each species are taken from the adjacent Partially Saturated Porous
Media domain.
BED PARAMETERS
Here you can specify the bed porosity, which is the void fraction in the packed bed
structure. Select From densities to calculate the porosity from the bed density and the
individual pellet density. Select User defined to specify the porosity directly.
Depending on the shape selection, equivalent radii or volumes and surface areas will
be required as input. If a size distribution is selected, the volume percentage of each
size is required as input.
Note that different chemical reactions can be specified for each pellet size if a
distribution is specified.
SURFACE SPECIES
In order to add surface species, click the Add button and enter the species name in the
Surface species table. Added surface species are be available inside all pellet types
defined in the Pellet Shape and Size section, but not in the bulk fluid.
For each pellet type, specify the Reactive specific surface area, Sb,reac (SI unit: 1/m),
corresponding to the surface area, per volume, available for surface reactions.
PELLET PARAMETERS
Enter a Pellet porosity εpe (dimensionless) to specify the porosity of the pellet internals.
Select Diffusion model — Millington and Quirk model (the default), Bruggeman model,
Tortuosity model, or User defined to describe the effective correction of the diffusion
coefficient in the pellet. In the case of the Tortuosity model, a value for the tortuosity
τpe within the pellet is required.
Enter also the Diffusion coefficient Dpe,c (SI unit: m2/s). If a User defined diffusion
model is selected, an Effective diffusion coefficient Dpeff,c (SI unit: m2/s) is entered. The
default value is 1·10-9 m2/s in both cases.
PELLET-FLUID SURFACE
For the coupling of concentration between the pellet internals and the surrounding
fluid, two Coupling type options are available:
• Continuous concentration, assuming that all resistance to mass transfer to/from the
pellet is within the pellet and no resistance to pellet-fluid mass transfer is on the bulk
fluid side. The concentration in the fluid will thus be equal to that in the pellet pore
just at the pellet surface: c pe,i = c i . This constraint also automatically ensures flux
continuity between the internal pellet domain and the free fluid domain through
so-called reaction forces in the finite element formulation.
The Film resistance (mass flux) option computes the inward surface flux,
Ni,inward=hDi(ci-cpe,i). hDi is the mass transfer coefficient (SI unit: m/s) and is
calculated with the default Automatic setting from a dimensionless Sherwood number
expression or with User defined mass transfer coefficients.
The Active specific surface area (SI unit: m-1) is required to couple the mass transfer
between the pellets and the bed fluid. Select either the Automatic setting that calculates
the specific surface area from the shape information given above. User defined is also
available for explicit surface area specification.
The Sherwood number expression can be computed from three available expressions:
Frössling, Rosner, and Garner and Keey. The Frössling equation is the default and
probably the most commonly used for packed spheres. All of these are based on the
dimensionless Reynolds, Re, and Schmidt, Sc, numbers, which are computed from
Density and Dynamic viscosity. Select these to be taken either From material or choose
the User defined alternative.
PELLET DISCRETIZATION
The extra dimension in the pellet needs to be discretized into elements. Select a
Distribution — Cubic root sequence (the default), Linear, or Square root sequence. Enter
the Number of elements Nelem.
CONSTRAINT SETTINGS
To display this section, click the Show button ( ) and select Advanced Physics Options.
See the details about the different constraint settings in the section Constraint
Reaction Terms in the COMSOL Multiphysics Reference Manual.
FURTHER READING
Theory for the Reactive Pellet Bed in the Theory section of this manual.
Reactions
The Reactions subfeature to the Reactive Pellet Bed is used to define reaction terms to
the transport within the reactive pellets. The feature also defines the corresponding
averaged heat sources to be applied to heat transport in the bulk fluid.
DOMAIN SELECTION
From the Selection list, choose the domains on which to define rate expression or
expressions that govern source terms in the transport equations.
Several reaction nodes can be used to account for different reactions in different parts
for the modeling geometry.
REACTION RATES
Add a rate expression R (SI unit: mol/(m3·s)) for species i using a value or an
expression. One reaction rate per species and pellet type can be entered.
Note that if you have the Chemistry interface available, provided with the Chemical
Reaction Engineering Module, the reaction rate expressions can be automatically
generated and picked up using the drop-down menu. For an example, see the
application Fine Chemical Production in a Plate Reactor.
Specify the rate expression R (SI unit: mol/(m2·s)) corresponding to the surface
reaction rate of each volumetric species i participating in the surface reaction.
Furthermore, specify the surface reaction rates for the participating surface species in
the corresponding table.
Note that if you have the Chemistry interface available, provided with the Chemical
Reaction Engineering Module, the reaction rate expressions can be automatically
generated and picked up using the drop-down menu.
HEAT SOURCE
Specify the heat source originating from the heat of reaction of the chemical reactions
inside the pellet can be specified. Both heat sources from reactions in the fluid, and
heat sources resulting from surface reactions can be defined. When using several pellet
types, heat sources for each type can be added.
The heat sources are most conveniently picked up from a Chemistry feature that
defines the reaction rate and the heat of reactions. In that case, the Rate expression can
be selected from the drop-down menu. Else it can be set to User Defined.
The defined heat source can be used by a Heat Source feature in any of the heat
transfer interfaces.
Species Source
In order to account for consumption or production of species in porous domains, the
Species Source node adds source terms expressions Si to the right-hand side of the
species transport equations.
DOMAIN SELECTION
From the Selection list, choose the domains on which to define rate expression or
expressions that govern the source term in the transport equations.
If there are several types of domains, with subsequent and different reactions occurring
within them, it might be necessary to remove some domains from the selection. These
are then defined in an additional Species Source node.
SPECIES SOURCE
Add a source term Si (SI unit: mol/(m3·s)) for each of the species solved for. Enter a
value or expression in the field of the corresponding species.
ε hs = β h M m ( c mo – c mo,ref )
where βh is the coefficient of hygroscopic swelling, Mm is the molar mass, cmo is the
moisture concentration, and cmo,ref is the strain-free reference concentration.
It requires a license of either the MEMS Module or the Structural Mechanics Module.
The multiphysics feature will appear automatically if both the Transport of Diluted
Species and the Solid Mechanics interfaces are added to the same component. For the
most current information about licensing, please see See http://www.comsol.com/
products/specifications/.
FURTHER READING
More information about how to use hygroscopic swelling can be found in Hygroscopic
Swelling Coupling section in the Structural Mechanics Module User’s Guide.
More information about multiphysics coupling nodes can be found in the section The
Multiphysics Branch in the COMSOL Multiphysics Reference Manual.
Fracture
Use this node to model mass transport along thin fractures in porous media. The node
assumes that the transport in the tangential direction along the fracture is dominant,
as a result of lower flow resistance.
FRACTURE PROPERTIES
Specify a value for the Fracture thickness dfr.
Specify the Porosity, εp (dimensionless) of the porous matrix. This is by default taken
From material. Select User defined to instead enter a different value.
CONVECTION
Select an option from the Velocity field list to specify the convective velocity along the
fracture. For a consistent model, use a Fracture Flow feature in a Darcy’s Law interface
to compute the fluid flow velocity in the fracture.
For User defined, enter values or expressions for the velocity components in the table
shown.
The settings for the Diffusion, and Dispersion sections are the same as for
Porous Media Transport Properties.
The interface supports simulation of species transport along boundaries in 2D and 3D,
and axisymmetric components in 2D. The dependent variable is the molar
concentration, c. Modeling multiple species transport is possible, whereby the physics
interface solves for the molar concentration, ci, of each species i.
SETTINGS
The Label is the default physics interface name.
The Name is used primarily as a scope prefix for variables defined by the physics
interface. Refer to such physics interface variables in expressions using the pattern
<name>.<variable_name>. In order to distinguish between variables belonging to
different physics interfaces, the name string must be unique. Only letters, numbers, and
underscores (_) are permitted in the Name field. The first character must be a letter.
The default Name (for the first physics interface in the model) is dsf.
BOUNDARY SELECTION
If model geometry includes boundaries that should not be included in the mass
transfer simulation, remove those from the selection list.
TRANSPORT MECHANISMS
Mass transport due to diffusion is always included. Use the Convection check box,
available under Additional transport mechanisms, to control whether to also include
convective transport.
• When the Crosswind diffusion check box is selected, a weak term that reduces
spurious oscillations is added to the transport equation. The resulting equation
system is always nonlinear. There are two options for the Crosswind diffusion type:
- Do Carmo and Galeão — the default option. This type of crosswind diffusion
reduces undershoots and overshoots to a minimum but can in rare cases give
equation systems that are difficult to fully converge.
- Codina. This option is less diffusive compared to the Do Carmo and Galeão
option but can result in more undershoots and overshoots. It is also less effective
for anisotropic meshes. The Codina option activates a text field for the Lower
gradient limit glim. It defaults to 0.1[mol/m^3)/tds.helem, where tds.helem
is the local element size.
• For both consistent stabilization methods select an Equation residual. Approximate
residual is the default and means that derivatives of the diffusion tensor components
are neglected. This setting is usually accurate enough and is computationally faster.
If required, select Full residual instead.
INCONSISTENT STABILIZATION
To display this section, click the Show button ( ) and select Stabilization. By default,
the Isotropic diffusion check box is not selected, because this type of stabilization adds
artificial diffusion and affects the accuracy of the original problem. However, this
option can be used to get a good initial guess for underresolved problems.
DISCRETIZATION
To display all settings available in this section, click the Show button ( ) and select
Advanced Physics Options.
The Value type when using splitting of complex variables setting should in most pure
mass transfer problems be set to Real, which is the default. It makes sure that the
dependent variable does not get affected by small imaginary contributions, which can
occur, for example, when combining a Time Dependent or Stationary study with a
frequency-domain study. For more information, see Splitting Complex-Valued
Variables in the COMSOL Multiphysics Reference Manual.
Add or remove species variables in the model and also change the names of the
dependent variables that represent the species concentrations.
Enter the Number of species. Use the Add concentration ( ) and Remove
concentration ( ) buttons as needed.
FURTHER READING
Boundary, Edge, Point, and Pair Nodes for the Transport of Diluted
Species in Fractures Interface
The Transport of Diluted Species Interface has the following boundary, edge, point,
and pair nodes, listed in alphabetical order, available from the Physics ribbon toolbar
(Windows users), Physics context menu (Mac or Linux users), or by right-clicking to
access the context menu (all users).
• Adsorption • No Flux
• Concentration • Outflow
• Flux • Reactions
• Fracture • Species Source
• Inflow
• Initial Values
MATRIX PROPERTIES
Use the Porous material list to define a material specifying the matrix properties on the
current selection. By default the Domain material is used. The density of the porous
media is needed when modeling adsorption to the surface of the porous matrix. By
default Density ρ is set to from domain material.
adsorption
• For Langmuir:
KL c ∂c P K L c Pmax
c P = c Pmax -------------------- -------- = ---------------------------
-
1 + KL c ∂c ( 1 + KL c )
2
Enter a Langmuir constant kL,c (SI unit: m3/mol) and an Adsorption maximum
cp,max,c (SI unit: mol/kg):
• For Freundlich:
c N ∂c P cP
c P = K F -------- -------- = N -----
- Freundlich
c ref ∂c c
Enter a Freundlich constant kF,c (SI unit: mol/kg), a Freundlich exponent NF,c
(dimensionless), and a Reference concentration cref,c (SI unit: mol/m3).
∂c P ∂
cP = KP c -------- = ( K P c ) User defined
∂c ∂c
Concentration
Use this node to specify the species concentration on a fracture boundary (applied in
points in 2D and along edges in 3D). For example, a c = c0 condition specifies the
concentration of species c.
CONCENTRATION
Individually specify the concentration for each species. Select the check box for the
Species to specify the concentration, and then enter a value or expression in the
corresponding field. To use another boundary condition for a specific species, click to
clear the check box for the concentration of that species.
CONSTRAINT SETTINGS
To display this section, click the Show button ( ) and select Advanced Physics Options.
You can find details about the different constraint settings in the section Constraint
Reaction Terms in the COMSOL Multiphysics Reference Manual.
Flux
This node can be used to specify the species flux across a boundary of a porous fracture
(applied in points in 2D and along edges in 3D). The flux of species c is defined as
n ⋅ ( D e ∇c ) = N 0
INWARD FLUX
Specify the flux of each species individually. To use another boundary condition for a
specific species, click to clear the check box for the mass fraction of that species.
Fracture
Use this node to model mass transport along thin fracture surfaces situated inside
porous or solid material. The node assumes that the transport in the tangential
MATRIX PROPERTIES
Use the Porous material list to define a material specifying the matrix properties on the
current selection. By default the Boundary material is used.
Specify the Porosity, εp (dimensionless) of the porous matrix. This is by default taken
From material. Select User defined to instead enter a different value.
CONVECTION
Select an option from the Velocity field list to specify the convective velocity along the
fracture. For a consistent model, use The Fracture Flow Interface to compute the fluid
flow velocity.
For User defined, enter values or expressions for the velocity components in the table
shown.
The settings for the Diffusion, and Dispersion sections are the same as for Porous Media
Transport Properties.
Inflow
Use this node to specify all species concentrations at a fracture inlet. The condition is
applied in points in 2D and along edges in 3D
If you want to specify the concentration of a subset of the partaking species, this can
be done by using the Concentration node instead.
CONCENTRATION
For the concentration of each species c0,c (SI unit: mol/m3), enter a value or
expression.
CONSTRAINT SETTINGS
To display this section, click the Show button ( ) and select Advanced Physics Options.
You can find details about the different constraint settings in the section Constraint
Reaction Terms in the COMSOL Multiphysics Reference Manual.
FURTHER READING
See the theory chapter in the section Danckwerts Inflow Boundary Condition.
No Flux
This node can be used to specify that the species flux across a boundary of a porous
fracture is zero. The condition is applied in points in 2D and along edges in 3D.
Outflow
Set this condition at fracture outlets where species are transported out of the model
domain by fluid motion. The condition is applied in points in 2D and along edges in
3D. It is assumed that convection is the dominating transport mechanism across
outflow boundaries, and therefore that diffusive transport can be ignored, that is:
n ⋅ ( – D e ∇c ) = 0
Reactions
Use the Reactions node to account for the consumption or production of species
through chemical reactions in the fracture. Define the rate expressions as required.
BOUNDARY SELECTION
From the Selection list, choose the boundaries on which to define rate expression or
expressions that govern the source term in the transport equations.
Several reaction nodes can be used to account for different reactions in different parts
of the fracture.
REACTION RATES
Add a rate expression Ri for species i. Enter a value or expression in the field. Note that
if you have the Chemistry interface available, provided with the Chemical Reaction
REACTING VOLUME
When specifying reaction rates for a species in a fracture, the specified reaction rate may
have the basis of the pore volume of the fracture, or the total volume.
• For Total volume, the reaction expressions in are specified per unit volume of the
fracture. The reaction expressions will be multiplied by the fracture thickness dfr.
• For Pore volume, the reaction expressions in mol/(m3·s) are specified per unit
volume of total pore space in the fracture. The reaction expressions will be
multiplied by the fracture thickness dfr and the fracture porosity, εp.
Species Source
In order to account for consumption or production of species in a fracture, the Species
Source node adds source terms expressions Si to the right-hand side of the species
transport equations.
BOUNDARY SELECTION
From the Selection list, choose the boundaries on which to define expressions that
govern the source term in the transport equations.
If there are several different parts of the fracture, with subsequent and different sources
occurring within them, it might be necessary to remove some boundaries from the
selection. The sources in these can then be defined using an additional Species Source
node.
SPECIES SOURCE
Add a source term Si for each of the species solved for. Enter a value or expression in
the field of the corresponding species.
The physics interface solves for the mass fractions of all participating species. Transport
through convection, diffusion, and migration in an electric field can be included.
The available transport mechanisms and diffusion models differs between various
COMSOL products (see http://www.comsol.com/products/specifications/).
Some examples of what can be studied with this physics interface include:
When this physics interface is added, the following default nodes are also added in the
Model Builder — Transport Properties, No Flux, and Initial Values. Then, from the Physics
toolbar, add other nodes that implement, for example, boundary conditions and
reactions. You can also right-click Transport of Concentrated Species to select physics
features from the context menu.
The Name is used primarily as a scope prefix for variables defined by the physics
interface. Refer to such physics interface variables in expressions using the pattern
<name>.<variable_name>. In order to distinguish between variables belonging to
different physics interfaces, the name string must be unique. Only letters, numbers, and
underscores (_) are permitted in the Name field. The first character must be a letter.
The default Name (for the first physics interface in the model) is tcs.
EQUATION
The basic equation for an individual species i is:
∂ ρω
( i ) + ∇ ⋅ ( ρω i u ) = – ∇ ⋅ j i + R i (3-100)
∂t
The displayed formulation changes depending on the active transport mechanisms and
the selected diffusion model.
TRANSPORT MECHANISMS
The Transport of Concentrated Species interface always accounts for transport due to
diffusion.
The available diffusion models and the additional transport mechanisms differs
between various COMSOL products (see http://www.comsol.com/products/
specifications/).
Diffusion Model
• The Maxwell-Stefan option employs the most detailed diffusion model, but is also the
most computationally expensive. The model is intended for diffusion dominated
models, and requires that the multicomponent Maxwell-Stefan diffusivities of all
component pairs are known. No stabilization is available when selecting this model.
• The Mixture-averaged option is less computationally expensive than the
Maxwell-Stefan model. It is a simpler model that can be used when variations in the
partial pressures and temperature can be assumed to not affect the multicomponent
diffusion. The model includes stabilization but requires the multicomponent
Maxwell-Stefan diffusivities of all component pairs.
• The Fick’s law model is a general model that should be used when the diffusion is
assumed Fickian, or when no multicomponent diffusivities are available. Also, when
molecular diffusion is not the dominating transport mechanism and a robust but
Knudsen Diffusion
For Mixture-averaged and Fick’s law, it is possible to include Knudsen diffusion. This
mechanism accounts for species collisions with the surrounding media, for example,
the pore walls the species pass through. It is also an important component when setting
up a Dusty gas model.
Q T
Di
j i = – ρω i D̃ ik d k – -------- ∇T
T
k=1
˜
where D ik (SI unit: m2/s) are the multicomponent Fick diffusivities, dk (SI unit: 1/
T
m) is the diffusional driving force, T (SI unit: K) is the temperature, and D i (SI
unit: kg/(m·s)) is the thermal diffusion coefficient.
Q
1
d k = ∇x k + --- ( x k – ω k ) ∇p – ρω k g k + ω k
p ρωl gl (3-101)
l=1
zk F
g k = – --------- ∇φ (3-102)
Mk
where zk is the species charge number, F (SI unit: A·s/mol) is Faraday’s constant and
φ (SI unit: V) is the electric potential.
SPECIES
Select the species that this physics interface solves for using the mass constraint in
Equation 3-103 (that is, its value comes from the fact that the sum of all mass fractions
must equal 1). In the From mass constraint list, select the preferred species. To
minimize the impact of any numerical and model introduced errors, use the species
with the highest concentration. By default, the first species is used.
ω1 = 1 – ωi (3-103)
i=2
• There are two consistent stabilization methods available when using the
Mixture-Averaged Diffusion Model or Fick’s Law Diffusion Model—Streamline
diffusion and Crosswind diffusion. Both are active by default.
The Residual setting applies to both the consistent stabilization methods.
Approximate residual is the default setting and it means that derivatives of the
diffusion tensor components are neglected. This setting is usually accurate enough
and computationally faster. If required, select Full residual instead.
• There is one inconsistent stabilization method, Isotropic diffusion, which is available
when using the Mixture-Averaged Diffusion Model or Fick’s Law Diffusion Model.
ADVANCED SETTINGS
To display this section, click the Show button ( ) and select Advanced Physics Options.
Normally these settings do not need to be changed.
Regularization
From the Regularization list, select On (the default) or Off. When turned On, regularized
mass fractions are calculated such that
0 ≤ w i, reg ≤ 1
Diffusion
The Diffusion settings are available for the approximate diffusion models
Mixture-averaged and Fick’s law.
When the Mixture diffusion correction is enabled, a flux correction is added to ensure
that the net diffusive flux is zero. This typically also mean that the solution becomes
less sensitive to the species selected to be computed from the mass constraint in the
Species section. More information on this correction is available in the theory section
Multicomponent Diffusion: Mixture-Averaged Approximation.
The Diffusion flux type list controls the whether the molecular flux is assumed
proportional to the mole fraction or the mass fraction. See Multicomponent Diffusion:
Mixture-Averaged Approximation or Multispecies Diffusion: Fick’s Law
Approximation for information on the diffusive flux formulation.
DISCRETIZATION
To display all settings available in this section, click the Show button ( ) and select
Advanced Physics Options.
For more information about these settings, see the Discretization section under The
Transport of Diluted Species Interface.
In the COMSOL Multiphysics Reference Manual see Table 2-3 for links
to common sections and Table 2-4 to common feature nodes. You can
also search for information: press F1 to open the Help window or Ctrl+F1
to open the Documentation window.
DEPENDENT VARIABLES
Add or remove species in the model and also change the names of the dependent
variables that represent the species concentrations.
Specify the Number of species. There must be at least two species. To add a single
species, click the Add concentration button ( ) under the table. To remove a species,
The species are dependent variables, and their names must be unique with
respect to all other dependent variables in the component.
In the COMSOL Multiphysics Reference Manual see Table 2-3 for links
to common sections and Table 2-4 to common feature nodes. You can
also search for information: press F1 to open the Help window or Ctrl+F1
to open the Documentation window.
• Open Boundary
• Outflow
Transport Properties
The Transport Properties is the main node used to model mass transfer in a fluid
mixture with the Transport of Concentrates species interface. The node adds the
equations governing the mass fractions of all present species, and provides inputs for
the transport mechanisms and for the material properties of the fluid mixture.
The settings in this node are dependent on the check boxes selected under Transport
Mechanisms in the Settings window of the Transport of Concentrated Species
interface.
When the Convection check box is selected, the Turbulent Mixing subnode is available
from the context menu as well as from the Physics toolbar, Attributes menu.
The options available in this feature differs between COMSOL products. (See http:/
/www.comsol.com/products/specifications/).
MODEL INPUTS
Specify the temperature and pressure to be used in the physics interface. The
temperature model input is used when calculating the density from the ideal gas law,
Temperature
Select the source of the Temperature field T:
• Select User defined to enter a value or an expression for the temperature (SI unit: K).
This input is always available.
• If required, select a temperature defined by a Heat Transfer interface present in the
model (if any). For example, select Temperature (ht) to use the temperature defined
by the Heat Transfer in Fluids interface with the ht name.
Absolute Pressure
Select the source of the Absolute pressure p:
• Select User defined to enter a value or an expression for the absolute pressure
(SI unit: Pa). This input is always available.
• In addition, select a pressure defined by a Fluid Flow interface present in the model
(if any). For example, select Absolute pressure (spf) to use the pressure defined in a
Laminar Flow interface with spf as the Name.
DENSITY
Define the density of the mixture and the molar masses of the participating species.
Mixture Density
Select a way to define the density from the Mixture density list — Ideal gas or User
defined:
• For Ideal gas, the density is computed from the ideal gas law in the manner of:
pM
ρ = -----------
Rg T
Here M is the mean molar mass of the mixture and Rg is the universal gas constant.
The absolute pressure, p, and temperature, T, used corresponds to the ones defined
in the Model Inputs section.
• For User defined enter a value or expression for the Mixture density ρ.
CONVECTION
This section is available when the Convection check box is selected.
• Select User defined to enter manually defined values or expressions for the velocity
components. This input is always available.
• Select a velocity field defined by a Fluid Flow interface present in the model (if any).
For example, select Velocity field (spf) to use the velocity field defined by the Fluid
Properties node fp1 in a Single-Phase Flow, Laminar Flow interface with spf as the
Name.
DIFFUSION
Specify the molecular and thermal diffusivities of the present species based on the
selected Diffusion model.
Specify the molecular and thermal diffusivities of the present species based on the
selected Diffusion model.
KNUDSEN DIFFUSION
The Knudsen diffusion transport mechanism accounts for the interaction of the species
with the surroundings (interspecies collisions excluded) — for example, the pore wall
when a species passes through porous media.
Depending on which Diffusion model is selected, either the Fick’s law or the
M
Mixture-averaged diffusion coefficient D i is corrected with the Knudsen diffusion
K
coefficient D i in the following way
MK 1 1 -
-1
Di = --------
- + -------
DM i Di
K
For gases, the Kinetic gas theory is often valid and requires the Mean path length λpath
(SI unit: m). Typically, for transport in porous media, the pore diameter can be entered
here. For other cases, choose User defined.
Settings for the mobilities are used for the Mixture-averaged and Fick’s law transport
models. By default the mobility is set to be calculated based on the species diffusivities
and the temperature using the Nernst-Einstein relation. To manually specify the
mobilities, select User defined for the mobility um,c (SI unit: s·mol/kg) and enter one
value for each species.
Enter the Charge number zc (dimensionless, but requires a plus or minus sign) for each
species.
The temperature (if you are using mobilities based on the Nernst-Einstein relation) is
taken from the Model Inputs section.
The settings in this node are dependent on the check boxes selected under Transport
Mechanisms in the Settings window of the Transport of Concentrated Species
interface.
MODEL INPUTS
Specify the temperature and pressure to be used in the physics interface. The
temperature model input is used when calculating the density from the ideal gas law,
but also when thermal diffusion is accounted for by supplying thermal diffusion
coefficients. The pressure model input is used in the diffusional driving force in
Temperature
Select the source of the Temperature field T:
• Select User defined to enter a value or an expression for the temperature (SI unit: K).
This input is always available.
• If required, select a temperature defined by a Heat Transfer interface present in the
model (if any). For example, select Temperature (ht) to use the temperature defined
by the Heat Transfer in Fluids interface with the ht name.
Absolute Pressure
Select the source of the Absolute pressure p:
• Select User defined to enter a value or an expression for the absolute pressure
(SI unit: Pa). This input is always available.
• In addition, select a pressure defined by a Fluid Flow interface present in the model
(if any). For example, select Absolute pressure (spf) to use the pressure defined in a
Laminar Flow interface with spf as the Name.
MATRIX PROPERTIES
Enter a value or expression for the Porosity, εp (dimensionless), of the porous media.
In order use a porosity defined in a material; specify the material using the Porous
material list, and select From material from the Porosity list.
DENSITY
Use this section to define the mixture density, and to specify the molar masses of the
participating species.
Mixture Density
Select a way to define the density from the Mixture density list—Ideal gas or User
defined:
• For Ideal gas, the density is computed from the ideal gas law in the manner of:
Here M is the mean molar mass of the mixture and Rg is the universal gas constant.
The absolute pressure, p, and temperature, T, used corresponds to the ones defined
in the Model Inputs section.
• For User defined enter a value or expression for the Mixture density ρ.
Molar Mass
Enter a value or expression for the Molar mass Mw for each species. The default value
is 0.032 kg/mol, which is the molar mass of O2 gas.
CONVECTION
This section is available when the Convection check box is selected.
• Select User defined to enter values or expressions for the velocity components. This
input is always available.
• Select a velocity field defined by a Fluid Flow interface present in the model (if any).
For example, select Velocity field (spf) to use the velocity field defined by the Fluid
Properties node fp1 in a Single-Phase Flow, Laminar Flow interface with spf as the
Name.
DIFFUSION
Specify the species molecular and thermal diffusivities in nonporous media in the
manner described for the Transport Properties node.
To account for the effect of porosity in the diffusivities, select an Effective diffusivity
model — Millington and Quirk model, Bruggeman model, Tortuosity model, or No
correction. Using one of the first four models, the effective transport factor, fe, is
defined from the porosity and the fluid tortuosity factor in the manner of:
εp
f e = ----- (3-104)
τF
• Select User defined to enter a value or expression for the electric potential. This input
is always available.
• If required, select an electric potential defined by an AC/DC interface that is
present in the model (if any). For example, select Electric potential (ec) to use the
electric field defined by the Current Conservation node cucn1 in an Electric
Currents interface ec.
Settings for the mobilities are needed for the Mixture-averaged and Fick’s law diffusion
models. By default the mobility is set to be calculated based on the species diffusivities
(adjusted by the Effective diffusivity model in the Diffusion section) using the
Nernst-Einstein relation. To manually specify the mobilities, select User defined for the
mobility um,w (SI unit: s·mol/kg) and enter one value for each species.
Enter the Charge number zc (dimensionless, but requires a plus or minus sign) for each
species.
The temperature (if you are using mobilities based or the Nernst-Einstein relation) is
taken from the Model Inputs section.
The flux is proportional to the current densities and the stoichiometric coefficients
according to Faraday’s law as defined by summation over the Reaction Coefficients
subnodes. The molar fluxes are multiplied by the species molar masses to obtain the
corresponding mass fluxes.
Select Account for Stefan velocity to update the Stefan velocity in accordance with the
mass flux from the electrode reactions. One example that may benefit from this is when
modeling gas diffusion electrodes.
Turbulent Mixing
Use this node to account for the species mixing caused by turbulent flow, typically
when the specified velocity field corresponds to a RANS solution. The node defines
the turbulent mass diffusion from the turbulent kinematic viscosity and a turbulent
Schmidt number.
The Turbulent Mixing feature should not be used together with the
Reacting Flow coupling feature. In this case the turbulent mixing is added
by the coupling feature.
The Turbulent Mixing subnode is available from the context menu (right-click the
Transport Properties parent node) or from the Physics toolbar, Attributes menu.
TURBULENT MIXING
Physics interfaces capable of solving for turbulent fluid flow provide the turbulent
kinematic viscosity, and these appear as options in the Turbulent kinematic viscosity νT
(SI unit: m2/s) list. The list always contains the User defined option that makes it
possible to enter any value or expression.
REACTION RATE
Select a Reaction rate — Automatic (the default), or User defined. Selecting Automatic
the laminar flow reaction rate is computed using the mass action law
For User defined, input a custom expression or constants for the Reaction rate r.
Specify the reaction stoichiometry by entering values for the stoichiometric coefficients
(dimensionless) of each species. Enter negative values for reactants and positive values
for products.
RATE CONSTANTS
When the Use Arrhenius expressions check box is not selected, input custom expressions
or constants for the Forward rate constant kf and Reverse rate constant kr.
When the Use Arrhenius expressions check box is selected, enter values for the following
parameters of the forward and reverse reactions:
TURBULENT FLOW
Note this section is only available when then licensed to the CFD Module (see http:/
/www.comsol.com/products/specifications/).
When the Turbulent-reaction model is set to None, laminar flow is assumed and the
reaction source terms are defined from the reaction stoichiometry and reaction rates
prescribed.
The Eddy-dissipation model also requires an estimation of the turbulent mixing time of
the fluid flow turbulence. When a Fluid Flow interface defining it is present in the
model, it can be selected from the Turbulence time scale list. For example, select
REGULARIZATION
Select Rate expression in order to regularize the individual rate expressions that are
added to each species. If the mass fraction for a reactant species ωi becomes smaller
than its damping limit, ωidl, the rate expression added to species ωi is reduced linearly.
If ω i ≤ 0 for a reactant species, the reaction rate contribution to that species is
completely removed. Similarly, the if the mass fraction for a product species ωj becomes
larger than 1−ωjdl, the rate expression added to that species is damped linearly. If
ω j ≥ 1 for a product species, the reaction rate contribution to that species is
completely removed.
The default value for the damping limit, ωidl, is 1e−6, which is appropriate for most
applications, but can require adjustment when working with for example catalytic trace
species.
Reaction Sources
In order to account for consumption or production of species due to one or more
reactions, the Reaction Sources node adds source terms to the right-hand side of the
species transport equations.
REACTIONS
Add an expression for the reaction mass source, Ri, for each individual species present,
except for the one computed from the mass constraint (see Species). Enter a value or
expression in the field for the corresponding species.
Select the Mass transport to other phases check box if mass is leaving or entering the
fluid as a result of the reactions, for instance due to condensation or vaporization in a
porous matrix. In this case the mass source for the species calculated from the mass
constraint can also be specified. The net mass transfer corresponds to the sum of the
mass sources for all species.
• For Total volume, the reaction mass source expressions are specified per unit volume
of the model domain.
• For Pore volume, the reaction mass source expressions are specified per unit volume
pore space. In this case the reaction mass sources will be multiplied by the domain
porosity εp (εp equals unity for nonporous domains).
Initial Values
The Initial Values node adds initial values for the mass fractions that can serve as an
initial condition for a transient simulation, or as an initial guess for a nonlinear solver.
If required, add additional Initial Values nodes from the Physics toolbar.
INITIAL VALUES
The initial mass fractions can be specified using a number of quantities. Select the type
of input from the Mixture specification list. Select:
Enter a value or expression in the field for each species except for the one computed
from the mass constraint.
• For Ideal gas, also specify the Initial pressure p0 and the Initial Temperature T0. Note
that dependent variables solved for are evaluated to zero for initial values. When
Mass Fraction
The Mass Fraction node adds boundary conditions for the species mass fractions. For
example, the following condition specifies the mass fraction of species i: ωi = ωi,0.
Set the mass fractions of all species except the one computed from the mass constraint.
This ensures that the sum of the mass fractions is equal to one (see Species). This node
is available for exterior and interior boundaries.
MASS FRACTION
Specify the mass fraction for each species individually. Select the check box for the
species to specify the mass fraction, and enter a value or expression in the
corresponding field. To use another boundary condition for a specific species, click to
clear the check box for the mass fraction of that species.
Select Account for Stefan velocity to update the Stefan velocity in accordance with the
prescribed mass fractions. Examples of cases that may benefit from this are for example
when modeling surface reactions or phase change on an exterior boundary.
CONSTRAINT SETTINGS
To display this section, click the Show button ( ) and select Advanced Physics Options.
Flux
The Flux node is available on exterior boundaries and can be used to specify the mass
flux across a boundary. The boundary mass flux for each species is defined in the
manner of:
– n ⋅ j i = j 0, i (3-105)
INWARD FLUX
Specify the Inward flux for each species individually. Select the check box for the species
to prescribe a flux and enter a value or expression for the flux in the corresponding
field. To use another boundary condition for a specific species, click to clear the check
box for the flux of that species. Use a positive value for an inward flux.
External convection
Set Flux type to External convection to prescribe a mass flux to or from an exterior
domain (not modeled) assumed to include convection. The exterior can for example
include a forced convection to control the temperature or to increase the mass
transport. In this case the prescribed mass flux corresponds to
j 0, i = k ω, i ( ω b, i – ω i ) (3-106)
where k ω, i is a mass transfer coefficient and ωb,i is the bulk mass fraction, the typical
mass fraction far into the surrounding exterior domain.
Select Account for Stefan velocity to update the Stefan velocity in accordance with the
prescribed flux. Examples of cases that may benefit from this are for example when
modeling surface reactions or phase change on an exterior boundary.
Inflow
The Inflow node adds a boundary condition for an inflow boundary where one
condition for each species is specified. It is available for exterior boundaries. The
condition can be specified using the following quantities:
A concentration quantity other than the mass fractions can only be used when all
species are defined, as in this boundary condition. The other quantities are
composition dependent and therefore unambiguous only when all species are defined.
For this reason the Mass Fraction node, which allows some species to use a different
boundary condition, only includes inputs for the mass fractions.
INFLOW
Select a Mixture specification:
Enter a value or expression in the field for each species except for the one computed
from the mass constraint.
CONSTRAINT SETTINGS
To display this section, click the Show button ( ) and select Advanced Physics Options.
No Flux
The No Flux node is the default boundary condition available for exterior boundaries.
It should be used on boundaries across which there is no mass flux, typically exterior
solid walls where no surface reactions occur. The condition applied for each species
corresponds to
–n ⋅ ji = 0
Outflow
The Outflow node is the preferred boundary condition at outlets where the species are
to be transported out of the model domain. It is useful, for example, in mass transport
models where it is assumed that convection is the dominating effect driving the mass
– n ⋅ ρω i Dik dk = 0
k
f
– n ⋅ ρD i ∇ω i = 0
Symmetry
The Symmetry node can be used to represent boundaries where the species
concentration is symmetric; that is, there is no mass flux in the normal direction across
the boundary:
– n ⋅ N = – n ⋅ ( ρω i u + j i ) = 0
This boundary condition is identical to the No Flux node, but applies to all species and
cannot be applied to individual species. The Symmetry node is available for exterior
boundaries.
Flux Discontinuity
The Flux Discontinuity node represents a discontinuity in the mass flux across an interior
boundary:
–n ⋅ ( Nd – Nu ) = N0 N = ( ρω i u + j i )
where the value of N0 specifies the size of the flux jump evaluated from the down to
the upside of the boundary.
FLUX DISCONTINUITY
Specify the jump in species mass flux. Use a positive value for increasing flux when
going from the downside to the upside of the boundary. The boundary normal points
Select the Species check boxes to specify a flux discontinuity, and enter a value or
expression for the Flux discontinuity N0 (SI unit: kg/(m2·s)) in the corresponding field,
N0, w1 for example. To use a different boundary condition for a specific species, click
to clear the check box for the flux discontinuity of that species.
Open Boundary
Use the Open Boundary node to set up mass transport across boundaries where both
convective inflow and outflow can occur. Use the node to specify an exterior species
composition on parts of the boundary where fluid flows into the domain. A condition
equivalent to the Outflow node applies to the parts of the boundary where fluid flows
out of the domain. The direction of the flow across the boundary is typically calculated
by a Fluid Flow interface and is entered as Model Inputs.
EXTERIOR COMPOSITION
Enter a value or expression for the species composition. Select:
• Mass fractions (the default) to enter mass fractions (ω0,ω1, for example)
• Mole fractions to enter mole fractions (x0,ω1, for example)
• Molar concentrations (SI unit: mol/m3) to enter molar concentrations (c0,ω1, for
example)
• Number densities (SI unit: 1/m3) to enter number densities (n0,ω1, for example)
and to describe the number of particles per volume n = n0
• Densities (SI unit: kg/m3) to enter densities (ρ0,ω1, for example)
A concentration quantity other than the mass fractions can only be used
when all species are defined.
Equilibrium Reaction
Use this node to model a reaction where the kinetics is assumed so fast that the
equilibrium condition is fulfilled at all times. The node solves for an additional degree
EQUILIBRIUM CONDITION
Selecting Equilibrium constant in the list, the following equilibrium condition based on
the species activities and the law of mass action is used
ν
∏ ai i
i ∈ products
K eq = ----------------------------------
–ν
∏ ai i
i ∈ reactants
where νi are the stoichiometric coefficients and the species activities are defined from
the concentration, ci, and the unit activity concentration ca0.
ci
a i = -------
c a0
Enter a value or expression for the dimensionless Equilibrium constant Keq, and the Unit
activity concentration Ca0.
Select User defined from the list to instead enter a manually defined Equilibrium
expression Eeq.
STOICHIOMETRIC COEFFICIENTS
Enter a value for the stoichiometric coefficientν for all participating species. Use
negative values for reactants and positive values for products in the modeled reaction.
Species with a stoichiometric coefficient value of 0 are not affected by the Equilibrium
Reaction node.
The settings for this node are the same as for Equilibrium Reaction except for the
setting in the section below.
N i = – D i ∇c i – z i u m, i Fc i ∇V + c i u (3-107)
where ci denotes the concentration of species i (SI unit: mol/ m3), Di is the diffusion
coefficient of species i (SI unit: m2/s), u is the solvent velocity (SI unit: m/s), F refers
to Faraday’s constant (SI unit: s·A/mol), V denotes the electric potential (SI unit: V),
zi is the charge number of the ionic species (dimensionless), and um,i its ionic mobility
(SI unit: s·mol/kg).
There are three physics interfaces for this type of transport. In the Nernst-Planck
Equations (npe) interface ( ), found under the Chemical Species Transport
branch ( ) when adding a physics interface, the transport of every charged species is
accounted for, and is solved for in combination with the electroneutrality condition.
The physics interface has the equations, boundary conditions, and rate expression
terms for modeling mass transport of ionic species by convection, diffusion, and
migration, solving for the species concentrations and the electric potential. The main
feature is the Convection, Diffusion, and Migration node, which adds the equation for
the species concentration and the potential and provides an interface for defining
species properties and model inputs.
When this physics interface is added, these default physics nodes are also added to the
Model Builder—Convection, Diffusion and Migration, No Flux (the default boundary
condition for species concentrations), Electric Insulation (the default boundary
SETTINGS
The Label is the default physics interface name.
The Name is used primarily as a scope prefix for variables defined by the physics
interface. Refer to such physics interface variables in expressions using the pattern
<name>.<variable_name>. In order to distinguish between variables belonging to
different physics interfaces, the name string must be unique. Only letters, numbers, and
underscores (_) are permitted in the Name field. The first character must be a letter.
The default Name (for the first physics interface in the model) is npe.
SPECIES
Select the species that the physics interface computes from the electroneutrality
condition in Equation 3-107. That is, its value comes from the fact that the net charge
in every control volume is zero. Select the preferred species in the From
electroneutrality list. To minimize the impact of any numerical errors, use the species
with the highest concentration. By default, the software uses the first species.
DEPENDENT VARIABLES
Add or remove species in the model and also change the names of the dependent
variables that represent the species concentrations.
Specify the Number of species. There must be at least two species. To add a single
species, click the Add species button ( ). To remove a species, select it in the list and
click the Remove species button ( ).Edit the names of the species directly in the
table.
The species are dependent variables, and their names must be unique with
respect to all other dependent variables in the model.
DISCRETIZATION
To display all settings available in this section, click the Show button ( ) and select
Advanced Physics Options.
For more information about these settings, see the Discretization section under The
Transport of Diluted Species Interface.
• When the Crosswind diffusion check box is selected, the Lower gradient limit glim
(SI unit: mol/m4) field defaults to 0.1[mol/m^3)/npe.helem, where npe.helem
is the local element size.
• For both consistent stabilization methods, select an Equation residual. Approximate
residual is the default setting and it means that derivatives of the diffusion tensor
components are neglected. This setting is usually accurate enough and is faster to
compute. If required, select Full residual instead.
DOMAIN NODES
• Convection, Diffusion, and Migration
• Reactions
• Equilibrium Reaction is described for the Transport of Diluted Species interface
• Initial Values
BOUNDARY CONDITIONS
• Concentration • Symmetry
• Flux • Inflow
• No Flux (the default boundary • Outflow
condition for species concentrations)
• Current Density
• Electric Insulation (the default boundary condition for the electric potential)
• Electric Potential (described for the Electrostatics interface in the COMSOL
Multiphysics Reference Manual)
For interior boundaries, continuity in the electric potential is the default boundary
condition. In addition, the following boundary conditions are available on interior
boundaries:
• Current Discontinuity
• Electric Potential (described for the Electrostatics interface in the COMSOL
Multiphysics Reference Manual)
• Select User defined to enter a value or an expression for the temperature (SI unit: K).
This input is always available.
• If required, select a temperature defined by a Heat Transfer interface present in the
component (if any). For example, select Temperature (ht) to use the temperature
defined by the Heat Transfer in Fluids interface with the ht name.
CONVECTION
This section is available when the Convection check box is selected.
The velocity field for the convection of the solvent, which is used in the convective
term, needs to be specified as feature input. Select the source of the Velocity field u:
• Select User defined to enter values or expressions for the velocity components. This
input is always available.
• Select a velocity field defined by a Fluid Flow interface present in the model (if any).
For example, select Velocity field (spf) to use the velocity field defined by the Fluid
Properties node fp1 in a Single-Phase Flow, Laminar Flow interface with spf as the
Name.
DIFFUSION
Select a Material from the list (when applicable and available). The default is None.
Specify the Diffusion coefficient Di for each species. This can be a scalar value for
isotropic diffusion or a tensor describing anisotropic diffusion. Enter the values in the
field (one for each species).
Enter the Charge number zi (dimensionless, but requires a plus or minus sign) for each
species.
Electric Insulation
The Electric Insulation node is the default boundary condition for the electric field. It
provides a boundary condition for the electric field that means that the boundary is
electrically insulated; that is, there is no current flowing across the boundary:
–n ⋅ i = 0 (3-108)
No Flux
The No Flux node is the default boundary condition on exterior boundaries. It should
be used on boundaries across which there is no mass flux, typically exterior solid walls
where no surface reactions occur. The condition applied for each species corresponds
to
–n ⋅ Ji = 0
Initial Values
The Initial Values node adds initial values for the concentrations ci and the electric
potential V that can serve as an initial condition for a transient simulation or as an initial
guess for a nonlinear solver. Add more Initial Values nodes from the Physics toolbar.
INITIAL VALUES
Enter values or expressions for the initial value of the Concentration ci (SI unit: mol/
m3) (except the one that a mass constraint determines) and for the Electric potential V
(SI unit: V). The default concentration values are 0 mol/m3 and the default electric
potential is 0 V.
Reactions
In order to account for consumption or production of species, the Reactions node adds
rate expressions Ri, which appear on the right-hand side of the species transport
equations.
Concentration
The Concentration node adds a boundary condition for the species concentrations: for
example the following condition specifies the concentration of species, ci:
c i = c i, 0
Set the concentration of all species except the one computed from the
electroneutrality condition, which ensures that the net charge is zero. See
Species.
CONCENTRATION
Specify the concentration ci (SI unit: mol/m3) for each species individually. Select the
check box for the species to specify the concentration, and enter a value or expression
in the corresponding field. To use another boundary condition for a specific species,
click to clear the check box for that species’ mass fraction.
CONSTRAINT SETTINGS
To display this section, click the Show button ( ) and select Advanced Physics Options.
Flux
The Flux node can be used to specify the species flux across a boundary. The species
flux, for example at the surface of a electrode, is defined as:
– n ⋅ J i = J 0, i
where n denotes the normal unit vector to the boundary. The flux J0,i can be any
function of concentration and electric potential (and indeed temperature if heat
transfer is included in the model), often provided by a kinetic expression. In
electrochemistry, this function usually depends on concentration and exponentially on
the electric potential through either the Tafel equation or the Butler-Volmer
equation.
Symmetry
The Symmetry node can be used to represents a boundaries where the concentration is
symmetric; that is, there is no species flux in the normal direction across the boundary:
–n ⋅ Ji = 0
Inflow
The Inflow node adds a boundary condition for an inflow boundary, where the
concentration for all species is specified in the manner of c i = c i, 0 .
CONCENTRATION
Specify the concentration ci (SI unit: mol/m3) of each species individually by entering
a value or expression in the corresponding field.
Outflow
The Outflow condition assumes that convection and migration is the dominating
transport mechanism across the outflow boundary, and therefore that diffusive
transport can be ignored, that is:
n ⋅ ( – D ∇c ) = 0
Flux Discontinuity
The Flux Discontinuity node represents a discontinuity in the species flux across an
interior boundary:
– n ⋅ ( N u – N d ) = N 0, i N = ( ci u + Ji )
where the value of N0,i specifies the size of the flux jump, evaluated from the up to the
down side of the boundary. This boundary condition is only available on interior
boundaries.
FLUX DISCONTINUITY
Specify the jump in species flux, using a positive value for increasing flux when going
from the down to the up side of the boundary. The normal direction (nx,ny,nz)
((nr,nz) in axisymmetric models) points in the direction from the downside toward
the upside of an internal boundary and can be plotted for visualization.
Select the check boxes for the species to specify a Flux discontinuity N0 (SI unit: mol/
(m2·s)) and enter a value or expression for the mass flux jump in the corresponding
field. To use another boundary condition for a specific species, click to clear the check
box for that species’ flux discontinuity.
Current Density
The Current Density node provides a boundary condition for the electric field which
makes it possible to set the total current density at a boundary. The condition is
defined as:
–n ⋅ i = i0 (3-109)
Current Discontinuity
The Current Discontinuity node represents a discontinuity in the current density across
an internal boundary:
where the value of i0 specifies the size of the current density jump evaluated from the
up to the down side of the boundary. This boundary condition is only available on
interior boundaries.
CURRENT DISCONTINUITY
2
Specify i0 (SI unit: A/m ), the jump in current density, using a positive value for
increasing current density when going from the up to the down side of the boundary.
The normal direction (nx,ny,nz) ((nr,nz) in axisymmetric models) points in the
direction from the downside toward the upside of an interior boundary and can be
plotted for visualization.
Open Boundary
Use the Open Boundary node to set up mass transport across boundaries where both
convective inflow and outflow can occur. Specify an exterior species concentration on
EXTERIOR CONCENTRATION
Enter a value or expression for the exterior concentration c0,c2 (SI unit: mol/m3).
A ( ) multiphysics node is also added. This node computes the local space charge,
based on the local concentrations and species charges in the Transport of Diluted
Species interface, and adds it to Poisson’s Equation in the Electrostatics interface.
The physics interface can simulate most forms of electrophoresis modes, such as zone
electrophoresis, isothachophoresis, isoelectric focusing and moving boundary
electrophoresis.
The interface supports simulation in 1D, 2D, and 3D as well as for axisymmetric
components in 1D and 2D.
The dependent variables are the electrolyte potential, and the molar concentrations of
the included species, added individual by each species node in the model tree.
SETTINGS
The Label is the default physics interface name.
The Name is used primarily as a scope prefix for variables defined by the physics
interface. Refer to such physics interface variables in expressions using the pattern
<name>.<variable_name>. In order to distinguish between variables belonging to
different physics interfaces, the name string must be unique. Only letters, numbers, and
underscores (_) are permitted in the Name field. The first character must be a letter.
The default Name (for the first physics interface in the model) is el.
OUT-OF-PLANE THICKNESS
For 1D components, enter a Cross sectional area Ac (SI unit: m2) to define
a parameter for the area of the geometry perpendicular to the 1D
component. The value of this parameter is used, among other things, to
automatically calculate the total current from the current density vector.
The analogy is valid for other fluxes. The default is 1 m2.
TRANSPORT MECHANISMS
Mass transport due to diffusion and migration is always included. Use the check boxes
available under Additional transport mechanisms to control other transport
mechanisms.
• By default, the Convection check box is selected. Clear the check box to disable
convective transport.
• The Mass transfer in porous media check box activates functionality specific to species
transport in porous media. When selected, the Porous Matrix Properties node can
be added to a domain to specify the electrolyte volume fraction and tortuosity, and
the Effective Transport Parameter Correction sections are enabled in the species
nodes.
• When the Crosswind diffusion check box is selected, a weak term that reduces
spurious oscillations is added to the transport equation. The resulting equation
system is always nonlinear.
• For both Streamline diffusion and Crosswind diffusion, select an Equation residual.
Approximate residual is the default and means that derivatives of the diffusivity are
neglected. This setting is usually accurate enough and is computationally faster. If
required, select Full residual instead.
INCONSISTENT STABILIZATION
To display this section, click the Show button ( ) and select Stabilization. By default,
the Isotropic diffusion check box is not selected, because this type of stabilization adds
artificial diffusion and affects the accuracy of the original problem. However, this
option can be used to get a good initial guess for under-resolved problems.
ADVANCED SETTINGS
To display this section, click the Show button ( ) and select Advanced Physics Options.
Normally these settings do not need to be changed. Select a Convective term —
Nonconservative form (the default) or Conservative form. The conservative formulation
should be used for compressible flow.
DISCRETIZATION
To display all settings available in this section, click the Show button ( ) and select
Advanced Physics Options.
The Compute boundary fluxes check box is activated by default so that COMSOL
Multiphysics computes predefined accurate boundary flux variables. When this option
is checked, the solver computes variables storing accurate boundary fluxes from each
boundary into the adjacent domain.
If the check box is cleared, the COMSOL Multiphysics software instead computes the
flux variables from the dependent variables using extrapolation, which is less accurate
in postprocessing results but does not create extra dependent variables on the
boundaries for the fluxes.
• <name>.nIl, where <name> is the name of the interface (default is el), set on the
interface top node. This is the normal electrolyte current density.
• <name>.ntflux_<species_name> is the Species name (see Common Settings for
the Species nodes in the Electrophoretic Transport Interface below). This is the
normal total flux for each species.
Also the Apply smoothing to boundary fluxes check box is available if the previous check
box is checked. The smoothing can provide a more well-behaved flux value close to
singularities.
For details about the boundary fluxes settings, see Computing Accurate Fluxes in the
COMSOL Multiphysics Reference Manual.
Regarding the Value type when using splitting of complex variables, see Splitting
Complex-Valued Variables in the COMSOL Multiphysics Reference Manual.
DEPENDENT VARIABLES
The dependent variable name for the electrolyte potential variable is phil by default.
The name of the concentration dependent variables are named as el.xxx, where the el is
the name of the interface as set above, and the xxx string is controlled by the Species
name setting on the individual species nodes.
FURTHER READING
• Ampholyte
• Fully Dissociated Species
• Uncharged Species
• Weak Acid
• Weak Base
Each species node add a dependent variable for the concentration. The initial and
boundary condition, as well as adding additional source reaction terms, for each
species concentration is controlled by adding subnodes to the species nodes:
• Concentration
• Flux
• Inflow
• Initial Concentration
• No Flux
• Outflow
• Species Source
The Ampholyte, Weak Acid and Weak Base are dissociation species and may define an
arbitrary number of dissociation steps. Each dissociation step is defined by it’s pKa (the
acid equilibrium constant) parameter. For the weak bases the pKa refers to the acid
constant of the conjugate acid. Each dissociation step adds one additional subspecies
concentration variable so that the concentration dependent variable represents the sum
of all subspecies, and initial and boundary conditions are defined with respect to this
total concentration.
All species node have a setting for the Species name, which needs to be unique. The
species name is used for naming of all related variables of the species. For species nodes
not defining any subspecies, the concentration variables are named as
<name>.c_<species_name> where <name> is the name of the interface (default is el),
set on the interface top node, and <species_name> is the Species name. For
dissociation species nodes defining multiple subspecies, the concentration nodes are
named as <name>.c<X>_<species_name> where <X> is the integer from 1 up to the
All species except the Uncharged Species carry charge and contribute to the total
electrolyte current which is used in the equation for solving the electrolyte potential.
The Immobile Species check box can be used to lock the concentration of a species, to,
for instance, define the immobile charges in a ion-selective membrane or a gel. When
the check box is enabled the concentration of the species is not added as a dependent
variable to the model; instead the concentration will be set to the value provided in the
Concentration field. The contribution to the electrolyte current for immobile species is
zero.
For dissociation species you may choose to set the transport parameters to be the Same
for all species appearing in the different dissociation steps, or you may use Individual
settings for each subspecies.
Typically the mobilities and diffusivities for small species are related by the
Nernst-Einstein relation, and when this relation is enabled you can choose whether to
specify either the Diffusivity (SI unit: m2/s) or the Mobility (SI-unit: s·mol/kg). The
Debye-Hückel-Henry relation is commonly used for larger molecules, such as proteins.
Note: There are other definitions of the migration transport equations in literature
which use mobilities expressed in m2/(V·s), whereas COMSOL Multiphysics uses
s·mol/kg. To convert mobilities expressed in m2/(V·s) to the corresponding values in
s·mol/kg, you typically divide by the Faraday constant, F_const (96485 C/mol).
The default correction model is Bruggeman, which multiplies the diffusivity and
mobility values by the porosity to the power of 1.5. The porosity of a domain is set by
the Porous Matrix Properties node.
The Debye-Hückel-Henry relation makes use of the ionic strength for calculating the
species mobility from the diffusivity. All charged species contribute to the ionic
strength, either assuming the species contributing to an Ideal solution or by using the
Lindestrøm-Lang assumption. The latter is usually used for macromolecules.
• Ampholyte • No Flux
• Concentration • Outflow
• Current • Porous Matrix Properties
• Current Density • Potential
• Current Source • Protein
• Flux • Species Source
• Fully Dissociated Species • Solvent
• Inflow • Uncharged Species
• Initial Concentration • Weak Acid
• Initial Potential • Weak Base
• Insulation
In the COMSOL Multiphysics Reference Manual, see Table 2-3 for links
to common sections and Table 2-4 for common feature nodes. You can
also search for information: press F1 to open the Help window or Ctrl+F1
to open the Documentation window.
Solvent
The settings of this node are used to define the properties of the aqueous solvent.
If Convection is enabled on the interface top node, you can specify the Velocity field
(m/s) as user defined input using analytical expressions or the velocity field variables
solved for by a separate physics interface.
The Solvent node automatically defines the concentration variables and for protons
(<name>.cH) and hydroxide ions (<name>.cOH), and the corresponding flux
expressions. See Diffusion and Migration Settings for how to set up the transport
parameters for the proton and hydroxide ions.
In the Solvent Properties section you can modify the Dynamic viscosity (Pa·s) and
Relative permittivity (unitless) values. The Built in and default values are applicable to
water. These parameters are used when calculating mobilities according to the
Debye-Hückel-Henry relation in the species nodes.
In the Water Self-Ionization section you can change the default Built in expression for
the Water self-ionization constant, pKw (unitless), to any user defined expression.
See also Common Settings for the Species nodes in the Electrophoretic Transport
Interface and Diffusion and Migration Settings.
Uncharged Species
Use this node to define a species that does not carry any charge, nor is impacted by the
electric field.
See also Common Settings for the Species nodes in the Electrophoretic Transport
Interface and Diffusion and Migration Settings.
Weak Acid
The Weak acid node supports multiple dissociation steps, where the acid of the first
dissociation step is uncharged.
The species may be either Monoprotic, subject to one dissociation step only, or
Polyprotic. For the latter case any Number of dissociation steps larger than one may be
used.
See also Common Settings for the Species nodes in the Electrophoretic Transport
Interface and Diffusion and Migration Settings.
Weak Base
The Weak base node supports multiple dissociation steps, where the base of the last
dissociation step is uncharged.
The species may be either Monoprotic, subject to one dissociation step only, or
Polyprotic. For the latter case any Number of dissociation steps larger than one may be
used.
Note that the pKa refers to the acid constant of the conjugate acid of the weak base.
See also Common Settings for the Species nodes in the Electrophoretic Transport
Interface and Diffusion and Migration Settings.
When using the Equilibrium constants the Base charge in last dissociation step, Z0
(unitless), needs to be set.
See also Common Settings for the Species nodes in the Electrophoretic Transport
Interface and Diffusion and Migration Settings.
Protein
Use the Protein node to define macromolecules. The features of the Protein node are
similar to the Ampholyte node, but with the default settings applicable for larger
molecules.
See also Common Settings for the Species nodes in the Electrophoretic Transport
Interface and Diffusion and Migration Settings.
Current Source
To make this node available, click the Show button ( ) and select Advanced Physics
Options.
Use this node to add a current source in a domain. A current source may appear in a
domain in homogenized porous electrode modeling, but should normally not be used.
Initial Potential
Use this node to specify the Initial Value of the electrolyte potential for the solver.
When using the Total current option in 1D or 2D, the boundary area is based either
on the Cross sectional area (1D) or the Out-of-Plane thickness (2D) properties, set on
the physics interface top node.
Current Density
Use the Current Density node to specify the current density distribution along a
boundary.
This node is typically used to model electrode surfaces where the electrode kinetics
depends on the electrolyte potential.
Note that using this node in 2D or 3D may result in an uneven potential distribution
along the boundary. To mitigate such effects you may use the Current node instead.
Insulation
The Insulation boundary condition describes the walls of a cell or the boundaries of the
cell that do not face an electrode (or a reservoir containing an electrode). The
boundary condition imposes the following equation:
il ⋅ n = 0
Potential
Add the Potential node on a boundary to set a fixed potential. This node is typically
used to model electrode surfaces or boundaries facing an electrolyte reservoir.
The node sets the potential in the electrolyte, φ l, to be equal to the Boundary
electrolyte potential, φ l, bnd (SI unit: V).
CONSTRAINT SETTINGS
To display this section, click the Show button ( ) and select Advanced Physics Options.
This node may be added as a subnode to any species node. See also Common Settings
for the Species nodes in the Electrophoretic Transport Interface.
Initial Concentration
This node may be added as a subnode to any species node. See also Common Settings
for the Species nodes in the Electrophoretic Transport Interface.
This node specifies the initial value for the Concentration, c (mol/m3), of the parent
species. This value serve as the initial condition for a transient simulation. The value
also serves as a start guess for stationary problems.
You can use spatially dependent functions (such as smoothed step functions) available
under Definitions when defining the Concentration expression to specify different
concentrations in different parts of the geometry. You can also use additional Initial
Values node and modify the Selection to set different values for different domains.
Concentration
This node may be added as a subnode to any species node. See also Common Settings
for the Species nodes in the Electrophoretic Transport Interface
This condition node adds a boundary condition for the parent species concentration.
Use the node to, for instance, specify the inlet concentration at the boundary facing a
electrolyte reservoir.
No Flux
This node may be added as a subnode to any species node. See also Common Settings
for the Species nodes in the Electrophoretic Transport Interface.
Flux
This node may be added as a subnode to any species node. See also Common Settings
for the Species nodes in the Electrophoretic Transport Interface.
This node can be used to specify the species inward flux across a boundary. The flux
can represent a flux from or into a much larger surrounding environment, a phase
change, or a flux due to chemical reactions.
INWARD FLUX
Enter a value or expression for the species mass flux J0. Use a minus sign when
specifying a flux directed out of the system.
External convection
Set Flux type to External convection to prescribe a flux to or from an exterior domain
(not modeled) assumed to include convection. The exterior can for example include a
forced convection to control the temperature or to increase the mass transport. In this
case the prescribed mass flux corresponds to
J0 = kc ( cb – c )
where kc is a mass transfer coefficient and cb is the bulk concentration, the typical
concentration far into the surrounding exterior domain.
Inflow
This node is available when you select the Convection check box on the physics interface
Settings window.
The other option, Flux (Danckwerts) can be more stable and fast to solve when high
reaction rates are anticipated in the vicinity of the inlet. Oscillations on the solutions
Outflow
Set this condition at outlets where species are transported out of the model domain by
migration or fluid motion. It is assumed that migration and convection is the
dominating transport mechanism across outflow boundaries, and therefore that
diffusive transport can be ignored, that is:
n ⋅ ( – D ∇c ) = 0
When this physics interface is added, these default nodes are also added to the Model
Builder—Surface Properties, No Flux, and Initial Values. Then, from the Physics toolbar,
add other nodes that implement, for example, boundary conditions. You can also
right-click Surface Reactions to select physics features from the context menu.
SETTINGS
The Label is the default physics interface name.
The Name is used primarily as a scope prefix for variables defined by the physics
interface. Refer to such physics interface variables in expressions using the pattern
<name>.<variable_name>. In order to distinguish between variables belonging to
different physics interfaces, the name string must be unique. Only letters, numbers, and
underscores (_) are permitted in the Name field. The first character must be a letter.
The default Name (for the first physics interface in the model) is sr.
DEPENDENT VARIABLES
Add or remove species and also change the names of the dependent variables that
represent the species concentrations. Note that the names can be changed but the
names of fields and dependent variables must be unique within a model.
Enter the Number of surface species. Use the Add surface concentration ( ) and
Remove surface concentration ( ) buttons as needed. The same number of Surface
concentrations cs, cs2, cs3, … are then listed in the table.
DISCRETIZATION
To display all settings available in this section, click the Show button ( ) and select
Advanced Physics Options.
For more information about these settings, see the Discretization section under The
Transport of Diluted Species Interface.
By default the Compensate for boundary stretching check box is selected for the Surface
Properties node. This section is then used to stabilize the tangential mesh velocity
term.
When the Compensate for boundary stretching check box is cleared (not selected), and
for fixed geometries or moving geometries, the stabilization has no effect.
• Boundary, Edge, Point, and Pair Nodes for the Surface Reactions
Interface
• Theory for the Surface Reactions Interface
Boundary, Edge, Point, and Pair Nodes for the Surface Reactions
Interface
The Surface Reactions Interface has these boundary, edge, point, and pair nodes, listed
in alphabetical order, available from the Physics ribbon toolbar (Windows users),
• Initial Values
• Reactions
• Surface Concentration
• Surface Properties
All other available nodes are described for the Transport of Diluted Species interface.
See Domain, Boundary, and Pair Nodes for the Transport of Diluted Species Interface.
In the COMSOL Multiphysics Reference Manual see Table 2-3 for links
to common sections and Table 2-4 to common feature nodes. You can
also search for information: press F1 to open the Help window or Ctrl+F1
to open the Documentation window.
Surface Properties
Use the Surface Properties node to define the density of sites, the site occupancy
number, and the surface diffusion.
SITES
Enter a value or expression for the Density of sites Γs (SI unit: mol/m2). The default
is 2 x 10-5 mol/m2.
Enter a Site occupancy number σi (dimensionless), indicating how many surface sites a
surface species block upon adsorption.
BULK SPECIES
For each bulk species enter the Molar mass Mi (SI unit: kg/mol) and the Density ρi
(SI unit: kg/m3). The default molar mass is 0.144 kg/mol and the default density is
5320 kg/m3).
Initial Values
The Initial Values node allows the initial value or guess for the surface and bulk
concentrations.
INITIAL VALUES
Based on the number of surface species and number of bulk species entered for the
physics interface under Dependent Variables section, enter values for the same number
of Surface concentration cs, cs2, cs3,... (SI unit: mol/m2) and Bulk concentration cb,
cb2, cb3,... (SI unit: mol/m2) in each field.
Reactions
The Reactions node adds rate expression terms to the species transport equations in
order to account for consumption or production of species due to reactions.
Surface Concentration
Use the Surface Concentration node to set the surface concentrations for one or more
species on an edge (3D components) or a point (2D and 2D axisymmetric
components).
SURFACE CONCENTRATION
Select each species check box as needed and enter a value or expression for each species
concentration, cs,0,cs1,cs2... (SI unit: mol/(m2·s)).
CONSTRAINT SETTINGS
To display this section, click the Show button ( ) and select Advanced Physics Options.
Selecting Laminar Flow under the Chemical Species Transport>Reacting Flow branch of
the Model Wizard or Add Physics windows, a Laminar Flow interface and a Transport of
Concentrated Species interface are added to the Model Builder.
In addition, the Multiphysics node is added, which includes the multiphysics coupling
feature Reacting Flow. The Reacting Flow feature predefines and controls the couplings
between the separate interfaces in order to facilitate easy set up of models.
In this section:
It combines the Laminar Flow, and Transport of Concentrated Species interfaces. The
Reacting Flow multiphysics coupling, which is added automatically, couples fluid flow
and mass transport. The fluid flow can either be free flow or flow in a porous medium.
The species transport supports both a mixture, where the concentrations are of
comparable order of magnitude, and low-concentration solutes in a solvent.
The interface can be used for stationary and time-dependent analysis in 2D, 2D axial
symmetry, and 3D.
The equations solved by the Laminar Flow interface are the Navier-Stokes equations for
conservation of momentum and the continuity equation for conservation of mass. A
Fluid Properties feature is active by default on the entire interface selection. A Fluid and
The Transport of Concentrated Species interface solves for an arbitrary number of mass
fractions. The species equations include transport by convection, diffusion and,
optionally, migration in an electric field.
The Reacting Laminar Flow interface triggers pseudo time stepping for the
flow equations when Use pseudo time stepping for stationary equation form
in the Fluid Flow interface is set to Automatic from physics.
The Name is used primarily as a scope prefix for variables defined by the coupling node.
Refer to such variables in expressions using the pattern <name>.<variable_name>. In
order to distinguish between variables belonging to different coupling nodes or physics
interfaces, the name string must be unique. Only letters, numbers, and underscores (_)
are permitted in the Name field. The first character must be a letter.
The default Name (for the first multiphysics coupling feature in the model) is rf1.
DOMAIN SELECTION
The Reacting Flow coupling is automatically defined on the intersection of the
selections for the coupled interfaces.
The Selection list displays the domains where the coupling feature is active.
COUPLED INTERFACES
This section defines the physics involved in the multiphysics coupling. The Fluid flow
and Species transport lists include all applicable physics interfaces.
• If it is added from the Physics ribbon (Windows users), Physics contextual toolbar
(Mac and Linux users), or context menu (all users), then the first physics interface
of each type in the component is selected as the default.
• If it is added automatically when a multiphysics interface is chosen in the Model
Wizard or Add Physics window, then the two participating physics interfaces are
selected.
You can also select None from either list to uncouple the node from a physics interface.
If the physics interface is removed from the Model Builder, for example Laminar Flow is
deleted, then the Species transport list defaults to None as there is nothing to couple to.
Click the Go to Source buttons ( ) to move to the main physics interface node for
the selected physics interface.
Click the Show or Hide Physics Properties Settings button ( ) to toggle the display of
physics properties settings affecting the coupling feature. When a turbulence model is
used, turbulent mass transfer is automatically accounted for (see the settings in the
Turbulence section below). Using Reacting Flow, the mass transfer treatment at walls
follows that applied for the fluid flow. Therefore the Wall treatment setting is also
displayed when using a turbulence model. For more information on turbulent mass
If a physics interface is deleted and then added to the model again, then
in order to reestablish the coupling, you need to choose the physics
interface again from the Fluid flow or Species transport lists. This is
applicable to all multiphysics coupling nodes that would normally default
to the once present interface. See Multiphysics Modeling Approaches in
the COMSOL Multiphysics Reference Manual.
TURBULENCE
When the fluid flow interface uses a turbulence model, select an option from the Mass
transport turbulence model list — Kays-Crawford, High Schmidt Number, or User-defined
turbulent Schmidt number.
For User-defined turbulent Schmidt number, enter a Turbulent Schmidt number ScT
(dimensionless).
The turbulent mass transfer added to the mass fraction equations is defined as
μT
N i, T = – --------- ∇ω i
Sc T
where μT is the turbulent viscosity defined by the flow interface, and the turbulent
Schmidt number, ScT, depends on the Mass transport turbulence model used.
Note, since the Reacting Flow coupling feature adds the turbulent mass transport, it
should not be combined with a Turbulent Mixing feature (subfeature to Transport
Properties in the Transport of Concentrated Species interface).
In addition, the Multiphysics node is added, which includes the multiphysics coupling
feature. The multiphysics coupling feature controls the coupling between the separate
interfaces in order to facilitate easy set up of models.
In this section:
It combines the Brinkman Equations, and Transport of Diluted Species in Porous Media
interfaces. The Reacting Flow, Diluted Species multiphysics coupling feature, which is
added automatically, couples the fluid flow and mass transport.
The interface can be used for stationary and time-dependent analysis in 2D, 2Daxi and
3D.
The Brinkman Equations interface computes the fluid velocity and pressure fields of
single-phase flow in porous media in the laminar flow regime. A Fluid and Matrix
Properties feature is active by default on the entire interface selection.
The Transport of Diluted Species in Porous Media interface computes the species
concentration in free and porous media, assuming that the species are of solutes,
dissolved in a solvent of significantly higher concentration. The species equations
include transport by convection, diffusion and, optionally, migration in an electric
field.
The interface can be used for stationary and time-dependent analysis in 2D, 2Daxi and
3D.
The Brinkman Equations interface computes the fluid velocity and pressure fields of
single-phase flow in porous media in the laminar flow regime. A Fluid and Matrix
Properties feature is active by default on the entire interface selection.
The Transport of Concentrated Species interface solves for an arbitrary number of mass
fractions in free and porous media. In the current multiphysics interface a Porous Media
Transport Properties feature is active by default on the entire interface selection.
BRINKMAN EQUATIONS
The available physics features for The Brinkman Equations interface are listed in the
Domain, Boundary, Point, and Pair Nodes for the Brinkman Equations Interface
section in the CFD Module User’s Guide.
The functionality for simulating fluid flow in free, porous media is found under
the Fluid Flow branch ( ) when adding a physics interface.
The documentation of the available interfaces and features can be found in CFD
Module User’s Guide in the chapters Single-Phase Flow Interfaces and Porous
Media and Subsurface Flow Interfaces.
289
Modeling Fluid Flow
In this section:
The Creeping Flow Interface ( ), also under the Single-Phase Flow branch ( ), is
used to model flow fluid flows at very low Reynolds numbers, also referred to as Stokes
flow. This typically occurs in fluid systems with high viscosity or small geometrical
length scales (for example in microfluidics and MEMS devices).
The rest of the available physics interfaces for fluid flow are found under the Porous
Media and Subsurface Flow branch ( ).
The Darcy’s Law Interface ( ) is used to model fluid movement through interstices
in a porous medium where a homogenization of the porous and fluid media into a
single medium is done. This physics interface combines the continuity equation and an
equation of state for the pore fluid (or gas), can be used to model low velocity flows,
for which the pressure gradient is the major driving force.
The Brinkman equations extend Darcy’s law to describe the dissipation of the kinetic
energy by viscous shear, similar to the Navier-Stokes equation. Consequently, they are
well suited to transitions between slow flow in porous media, governed by Darcy’s law,
and fast flow in channels described by the Navier-Stokes equations. The equations and
boundary conditions that describe these types of phenomena are also found in The
Free and Porous Media Flow Interface ( ).
More advanced descriptions of fluid flow, such as turbulent and multiphase flow, can
be found in the CFD Module. More extensive descriptions of heat transfer, such as that
in turbulent flow or systems including radiation, can be found in the Heat Transfer
Module. Furthermore, some applications that involve electrochemical reactions and
porous electrodes, particularly in electrochemical power source applications, are better
handled by the Batteries & Fuel Cells Module.
The functionality for simulating heat transfer in free and porous media is found
under the Heat Transfer branch ( ) when adding a physics interface.
The documentation of all features in the Heat Transfer in Porous Media interface
can be found in the Heat Transfer Module User’s Guide, in the chapter Heat
Transfer Interfaces.
293
Modeling Heat Transfer
Available Physics Interfaces
Heat transfer is an important part of chemical reaction engineering. Most chemical
reactions either require or produce heat, which in turn affects both the reactions
themselves and other physical processes connected to the system.
The heat transfer branch includes a number of physics interfaces for the simulation of
heat transfer, particularly in porous media. As with all other physical descriptions
simulated by COMSOL Multiphysics, any description of heat transfer can be directly
coupled to any other physical process. This is particularly important for systems based
on chemical reactions and mass transfer, along with fluid flow.
You can manipulate these physics interfaces with respect to including source and sink
terms or, in other words, account for exothermic or endothermic chemical reactions.
If you also have the Heat Transfer Module there are additional Heat
Transfer interface features, such as radiation, available.
More advanced descriptions of fluid flow, such as turbulent and multiphase flow,
require the CFD Module. In addition, the Heat Transfer Module also includes more
detailed descriptions and tools for simulating energy transport, such as radiation.
Thermodynamics
This chapter describes how you can use the thermodynamics functionality to
define thermo-physical and transport properties. The properties in turn can be used
when simulating chemical reaction systems, or any type of transport model
involving mass transfer, fluid flow, or heat transfer.
In this chapter:
295
Using Thermodynamic Properties
In this section:
2 If you use the Reaction Engineering interface or the Chemistry interface, link
chemical species in these interfaces with the chemical species in the thermodynamic
system. The required property parameters and functions are automatically added
and visualized as nodes under the corresponding Thermodynamic System node.
3 As an alternative to step 2 above, you can also set up functions or constants manually
for a thermodynamic system. You can, for example, create a function for the density
of a fluid for use in a fluid flow interface without having to use the Chemical
Reaction Engineering or the Chemistry interfaces. You can also define mixture
functions that depend on the composition of a mixture, for example for density,
enthalpy, or heat capacity.
4 Both alternatives in steps 3 and 4 above allow you to use the functions defined by a
thermodynamic system in any physics interface in COMSOL Multiphysics to
evaluate properties that depend on model variables such as temperature, pressure,
and composition.
Thermodynamics
To access the functionality for thermodynamic calculations, right-click the Global
Definitions node in the Model Builder tree and select Thermodynamics ( ). Using the
Use User Defined Species to add new species that are not available in COMSOL
database. You can also edit available species in the database.
Thermodynamic System
A thermodynamic system is used to describe properties of pure species and mixtures of
chemical compounds for liquids, gases, liquid-vapor equilibria, and liquid-liquid
equilibria. It specifies the available species and the phases (states of aggregation) that
may present in the modeled system. It also defines and evaluates the functions for
thermodynamic and transport properties of the chemical system, that is, the species
and mixture properties for liquids, gases and phase equilibria.
1 Select System
2 Select Species
3 Selecting a Thermodynamic Model
SELECT SYSTEM
Use the Select System step in the wizard to define the phases in the modeled system.
You can select Gas, Liquid, Vapor-liquid, Vapor-liquid-liquid or Liquid-liquid. The names
of the phases in the Selected system table can be changed by editing the elements in the
Name column. Click the Next button ( ) to proceed to the next step in the wizard.
SELECT SPECIES
Use the Species filter to search among the available species in the COMSOL database or
the User Defined Species. Species can be searched for by typing the name, CAS
number, or the chemical formula. Select one or more species in the list and click the
• Ideal Solution
• Regular Solution
• Extended Regular Solution
• Wilson
• NRTL
• UNIQUAC
• UNIFAC
Use the Advanced options check box in order to manually control the models used for
thermodynamic properties, transport properties and surface tension. When this check
box is selected, all available property models are shown in the property model table.
The available property models are dependent on the phases available in the
thermodynamic system.
When only a gas phase is present, models are available for following properties:
For a single-phase liquid, or for a two-phase liquid-liquid system, there are models for
following properties:
• Liquid volume
• Gas-liquid surface tension
• Liquid thermal conductivity
• Liquid viscosity
Click the Finish button ( ) to exit the Thermodynamic System Wizard and add the
corresponding system under the Thermodynamics node. Note that the default node
label reflects the available phase in the system. For example, when creating a
vapor-liquid system, a node labeled Vapor-Liquid System is added.
Provider:
Contains information about the thermodynamic system such as the version it was
created in and the manager used to create it. For a built-in thermodynamic system the
manager corresponds to COMSOL.
Species:
Lists the species included in the thermodynamic system. You can change the list of
species by adding or removing species in a thermodynamic system. For information on
how to change the list of species, see Add or Remove Species.
Species Properties:
Lists the parameters and function that describe thermodynamics and transport
properties for pure species. Such functions may describe density, heat capacity, thermal
conductivity, or other thermodynamic and transport properties. For more reading, see
Species Property.
Mixture Properties:
Lists the available mixture property functions. Here you can see the mixture models
that are available. Note that you have to have a mixture model defined in order to use
these functions. You can define a mixture property by right-clicking the
Thermodynamic System node and selecting Mixture Property.
Thermodynamic Models:
When a liquid is present, use the Liquid phase model list to select the thermodynamic
model for the this phase.
When a gas is present, use the Gas phase model list to select the thermodynamic model
for this phase. If a liquid phase is also available and the Liquid phase model corresponds
to an equation of state, then the Gas phase model is set to the same model.
Property Models:
Displays the settings for the available individual property models. The property models
can be changed by selecting different values from the combo box in the Model
column. For example, for Gas thermal conductivity you can select Kinetic theory or Ideal
from the corresponding combo box. You can also select the property model in the last
step of the thermodynamic system wizard by selecting the Advanced option check box,
see Selecting a Thermodynamic Model.
Figure 6-7: Entering the Binary Interaction Parameters for the Soave-Redlich-Kwong
model.
When a binary interaction parameter is not available in the database for a pair of
species, the parameter value is set to zero (default value).
Make sure to click the Finish button ( ) button in order for the changes to take
effect.
DEFINE SYSTEM
Selecting Define System takes you to the Select System step in the Thermodynamic
System Wizard. You can select the desired system from the list.
Note that changing the phase(s) in a thermodynamic system that is currently coupled
to a Reaction Engineering interface or a Chemistry interface breaks this coupling. See
the Coupling with the Reaction Engineering and the Chemistry Interfaces section for
how this coupling can be updated.
WARNING INFORMATION
A sanity check is always performed when a Thermodynamic System is created. If any
problems are found, a Warning Information node listing the problems is added under
the package. One example when this occurs is when a parameter required for the
thermodynamic model, typically a binary interaction parameter, is not available in the
database.
1 Download and install COCO, which includes the TEA thermodynamic system
manager. The software is available from
www.cocosimulator.org/index_download.html.
2 Create and configure a thermodynamic system template that handles physical and
thermodynamics calculations needed for your model. If you have already created a
package template earlier, or if an adequate thermodynamic system already exists in
the installation, this step is not needed.
3 Create an External Thermodynamic System node as detailed in the next section.
The installed packages are available in the Select Property Package step of the
Thermodynamic System Wizard.
When adding the package its default label reflects the included phases. For example, a
node labeled Vapor-Liquid System (External) is created when adding an external
package containing a vapor and a liquid phase.
Figure 6-9: Available external thermodynamic systems. Use the Thermodynamic System
Wizard to browse the contents of the installed external thermodynamic systems on your
system. The example shows the packages shipped with the COCO provider.
Species:
Lists the species included in the selected thermodynamic system.
Phases:
Lists the phases included in the selected thermodynamic system for example gas,
liquid, gas-vapor, or liquid-liquid.
Species Properties:
Lists the parameter values or the functions that describe species properties, for example
molar mass, and properties available for pure compounds, such as density as a function
of temperature.
Mixture Properties:
Lists the available functions that describe mixture properties. One example is the
density of a non-ideal mixture as a function of composition.
Figure 6-10: Settings seen when selecting an added External Thermodynamic System.
Note that only the thermodynamic system definition is exported, for example the
underlying species data (from the database), the included phases, and the applied
thermodynamic models. Property values and functions created using the
thermodynamic system are not stored.
Species Property
A Species Property is used to define and compute a pure species property. The available
properties consist of both parameters and functions. Some examples of available
parameters are molar mass, Lennard Jones diameter, and dipole moment. Some
examples of available functions are density, enthalpy, heat capacity, and viscosity. The
property functions created are either dependent on temperature alone, or both on
temperature and pressure. For all property functions, the first order derivative with
respect to temperature and, when applicable, with respect to pressure are automatically
defined.
Right-click the relevant Thermodynamic System node (see Figure 6-8), or the relevant
External Thermodynamic System node, and select Species Property to start the Species
Property Wizard.
SELECT PROPERTIES
First use the Amount base unit list to define the base unit. Select mol or kg.
Use the filter to search among the available properties. Select one or more properties
in the list and click the Add Selected button ( ) to add them to the Selected properties
list.
Click the Next button ( ) to proceed to the next step, selecting the species.
SELECT SPECIES
Select one or more of the species available in the thermodynamic system and use the
Add Selected button ( ) to add them to the Selected species table. One property
function is created for each of the selected species.
Click the Next button ( ) to proceed to the next step in the wizard, which is selecting
the phase for the property.
Definition:
Shows the definition of a property that is defined as a parameter or a function, for
example the name of the parameter or function.
You can use the Parameter name or Function name fields to specify or change the name
of a parameter or a function. The section also provides information about the property
type and the species it is defined for.
Plot Parameters
Available for property functions in order to plot a selected function for a given set of
argument values.
Apply a Lower limit and an Upper limit for each argument, and click the Plot button
( ) to plot the function using the given argument range. You can also click the Create
Plot button ( ) in order to create a plot group, under the Results node.
Properties Window
To see the reference for constant or temperature dependent functions, right-click on
the function and select Properties. This opens the Properties window.
Figure 6-16: Show constant and temperature dependent properties references for species
functions.
Mixture Property
A Mixture Property is used to compute a property function that depends on the
concentration of the species in a thermodynamic system. Some examples of the
available property functions are density, enthalpy, heat capacity, and surface tension.
Apart from the composition, the mixture property functions are also dependent on
temperature and pressure. The first order derivatives with respect to temperature and
pressure are automatically defined.
Right-click the relevant Thermodynamic System node (see Figure 6-8), or the relevant
External Thermodynamic System node, and select Mixture Property to start the Mixture
Property Wizard.
1 Select Properties
2 Select Phase
3 Select Species
4 Mixture Property Overview
SELECT PHASE
Use this list to specify the phase, among the ones available in the system, for the
selected mixture property.
Some properties require that the system consists of two phases, for example surface
tension. A two-phase system may consist of a combination of liquid-vapor or
liquid-liquid phases.
SELECT SPECIES
First select the Species composition base unit to be used for function arguments. Select
Mole fraction or Mass fraction.
Select the species to be included the list. Use the Add All button ( ) to add all species
in the thermodynamic system. It is also possible to select a subset of the available
species. In that case use the Add Selected button ( ) to add species. The Selected
species table is updated as you add species.
Click the Next button ( ) to proceed to the next step in the wizard.
SETTINGS
Selecting a desired Mixture Property node to display its settings window.
Figure 6-18: Mixture properties can be created by right-clicking a Mixture node under a
thermodynamic system.
Equilibrium Calculation
The Equilibrium Calculation functionality is used to compute the resulting equilibrium
conditions for a mixture of a set of species and phases.
Equilibrium calculations are often used for processes with vapor-liquid equilibrium
(VLE), so-called flash calculations. A typical process that requires flash calculations to
be described in a model is a distillation process, where a multiphase feed stream is
separated into a vapor and a liquid product and where the concentrations of the species
in each phase are required. Equilibrium calculations involve combining the
The last three examples are often considered more difficult, since they require energy
balances and relations for computing enthalpy and entropy. The Chemical Reaction
Engineering Module can handle all of the above cases using the equilibrium calculation
functionality. Phase envelopes, bubble point, and dew point can be calculated for any
number of species.
Right-click the relevant Thermodynamic System node (see Figure 6-8), or the relevant
External Thermodynamic System node, and select Equilibrium Calculation to start the
Equilibrium Calculation Wizard.
1 Select species
2 Equilibrium Specifications
3 Equilibrium Function Overview
• The first function type is used to detect whether a phase is present in the system and
includes “exist” in its name.
• The amount function computes the total amount of material in each phase.
• The phase composition functions compute the mass or mole fraction of each species
in each phase, depending on the selected base unit for the equilibrium calculation.
SELECT SPECIES
Select one or more of the species available in the thermodynamic system and use the
Add Selected button ( ) to add them to the Selected species table. Click the Next
button ( ) to proceed to the next step in the wizard.
Select two Equilibrium conditions that define the current equilibrium, for example a
given pressure and a temperature. These equilibrium conditions are used as arguments
in the equilibrium functions, in addition to the composition (overall fractions of
species).
The available equilibrium conditions are: Temperature, Pressure, Phase fraction, Energy
(or Energy of formation), Enthalpy (or Enthalpy of formation), Specific volume, and
Entropy (or Entropy of formation). For chemical reactions, it is recommended to use
Enthalpy of formation, Entropy of formation and Energy of formation, since they account
for heat of reactions.
Selecting Phase fraction as one of the equilibrium conditions activates the Solution type
input field. This can be used to indicate the direction of the desired solution, which is
of great use especially near critical points. The options Undefined, Normal, or Retrograde
define different directions for the search of the solution to the equilibrium equations.
Using Normal means that the derivative of the vapor phase fraction with respect to
temperature (at constant pressure and composition) is kept positive and the derivative
of the vapor phase fraction with respect to pressure (at constant temperature and
For a single species, the critical point is the highest pressure and temperature at which
two phases (liquid and vapor) are distinguishable. However, for some multispecies
systems, the critical point is a point between the dew point and the bubble point. In
this case, the critical point does not represent the maximum pressure or the maximum
temperature of vapor-liquid coexistence. This phenomena is known as retrograde
condensation. This means that under isothermal conditions, when the pressure
decreases, some of the vapor condenses into liquid instead of expanding or vaporizing.
An example of such system is formation of liquid hydrocarbons in a gas reservoir as the
pressure decreases below the dew point pressure. In this case, setting Solution type to
Normal or Retrograde may not be sufficient to distinguish between the two solutions.
Note that the Solution type setting is only available for a built-in thermodynamic
systems. For external thermodynamic systems, the corresponding functionality needs
to be supplied by the thermodynamic software provider. For instance, the
COCO/TEA provider does not support the Normal or Retrograde options. In those
cases, the Solution type should be Undefined.
SETTINGS
Selecting an Equilibrium Calculation node displays the settings including the property
functions, see Figure 6-20.
Examples of species properties that can be automatically created are the molar mass,
the heat capacity, the enthalpy, and the entropy of each species. Parameters and
functions for these properties are created by the package. The Reaction Engineering and
Chemistry interfaces can also be used to define transport properties for the resulting
mixture (all species in the interface). When coupled, the following mixture properties
can be automatically created: heat capacity, density, thermal conductivity, and dynamic
viscosity.
REACTION ENGINEERING
Mixture Properties
You can couple a Reaction Engineering interface with an existing thermodynamic
system in the Reaction Engineering interface’s settings window.
You need to have at least one species defined in the Reaction Engineering interface in
order to couple to a thermodynamic system. You can make this coupling in settings
window for the Reaction Engineering interface by selecting the Thermodynamics check
box in the Mixture Properties section.
Species Matching
The Species Matching section is activated when the Thermodynamics check box is
selected in the Mixture Properties section, see above. The species in the Reaction
Engineering interface can be matched to a species in the thermodynamic system. This
ensures that the arguments in the thermodynamic system functions are correctly
defined.
Use the drop-down lists in the From Thermodynamics column to match each species in
the interface to a species in the coupled thermodynamic system.
For each species matched, the required property parameters and functions are added
under to the corresponding thermodynamic system.
When all species are matched, the interface is considered fully coupled and functions
representing mixture properties, such as density, are also added automatically under
the corresponding thermodynamic system.
Figure 6-22: Matching the species in Reaction Engineering to those in the corresponding
thermodynamic system.
Figure 6-23: Select Calculate mixture properties when coupled to a thermodynamic system.
CHEMISTRY
Mixture Properties
Coupling a Chemistry interface to a Thermodynamic System is initiated by selecting the
Thermodynamics check box, in the same way as for the Reaction Engineering interface
above.
Species Matching
Here you can match the variables for the concentrations and by this calculate mixture
properties (transport and thermodynamic properties). For information on how to
specify the dependent variables to be used, see Species Matching in The Chemistry
Interface documentation.
For each species matched, the species property parameters and functions required by
the Chemistry interface are automatically created and added under the corresponding
thermodynamic system.
When the system is fully matched, the mixture property Zmix is defined as:
∂Z mix ( T, P, n 1, n 2, …, n m )
Z i = ----------------------------------------------------------------------- (6-2)
∂ni T, P, n
j≠i
The definition of partial molar properties can be rewritten using mole fraction
derivatives as:
Z i ( T , P , n 1 , …, n m ) = Z i ( T , P , x 1 , … , x m ) = (6-3)
∂Z mix ( T, P, x 1, …, x m )
Z + ------------------------------------------------------------
mix ( T, P, x 1, …, x m ) ∂xi T, P, x i
m ∂Z mix ( T, P, x 1, …, x m )
– xi -----------------------------------------------------------
∂xi
-
T, P, x
i=1
i
When the system is partially coupled, which means that some but not all species have
been coupled, the mixture property is instead calculated assuming ideal mixing:
m
Z mix = n i Z i ( T, P ) (6-4)
i=1
Figure 6-25: The enthalpy and heat capacity of benzene are evaluated using external
property calculations through a thermodynamics package. In this case, the functions are
temperature dependent properties for a pure component.
The data needed for a user defined species includes both material properties, such as
the molar mass and the vapor pressure, and properties for specific thermodynamic
models or transport models (see Thermodynamic Models and Theory). For instance,
in order to use the UNIFAC Thermodynamic model, you need to define UNIFAC
groups for the new species.
Note that the thermodynamics calculator includes measures to handle missing species
properties, sometimes by applying approximations. For instance, if the Wilson volume
is not available, the liquid volume at normal boiling point is used instead. If data for
this is also missing, it is estimated from the saturated liquid density correlation.
When creating a user defined species it is recommended to add the following common
material properties:
• Molecular mass
EXAMPLE MODEL
For an example of using a User Defined Species see this application example.
Figure 6-26: User Defined Species Wizard; Enter Name and Formula.
If you want to edit a species in the COMSOL database, select the Use species from
database as template check box. Then locate the species to edit in the list. The filter
text field can be used to search among the available species.
Figure 6-27: User Defined Species Wizard; edit an available species by selecting
Use species from database as template.
Constants
Use this table to define material constants such as molar mass, critical temperature, and
standard enthalpy of formation.
Figure 6-29: User Defined Species Wizard; Specify the structure information in the Enter
Parameters step.
Model parameters
Specify the parameters for the thermodynamic models and transport models in use.
DEFINE PROPERTIES
The last step in the wizard is to add temperature dependent properties for the new
species. All temperature dependent properties are defined using cubic polynomials on
the form
Each property can consist of an arbitrary number of temperature intervals, each using
the above form. Click the Add button ( ) under the table for a specific property to
add an interval.
Note that some thermodynamic properties, such as the enthalpy and entropy, of a
species or mixture, are estimated from the ideal gas heat capacity and depends on the
thermodynamic model applied for the system (see Thermodynamic Properties
Definitions).
Figure 6-31: User Defined Species Wizard; Specify temperature dependent properties.
SETTINGS
Selecting a species node under User Defined Species shows the Settings window
including the definitions of all species properties. Properties are categorized into
• Introduction
• Thermodynamic Models
• Selecting the Right Thermodynamic Model
• Thermodynamic Properties Definitions
• Standard Enthalpy of Formation and Absolute Entropy terms
• Reference State
• Transport Properties
• Surface Tension
• References
Introduction
In this chapter, we review the theory behind the thermodynamic properties database
and its functions. The thermodynamic models in the database are available for single
phase, gas or liquid, and phase equilibrium systems for two or more phases such as
vapor-liquid equilibrium (VLE), vapor-liquid-liquid equilibrium (VLLE) and
liquid-liquid equilibrium (LLE).
Thermodynamic Models
In the following sections, the available thermodynamic models are described:
PV
Z ≡ -------- = f v ( V, T ) (6-6)
RT
The equation of state models available in the thermodynamic properties database are:
RT
P = -------- (6-7)
V
As the name suggests, the ideal gas law is only applicable to gases. In fact, its use is
limited to gases at low to moderate pressures.
Peng-Robinson
The classical Peng-Robinson (PR) equation of state Ref. 7 is given by
RT aα
P = ------------- – ------------------------------------
- (6-8)
V – b V 2 + 2bV – b 2
2 2
R Tc
a i = Ω A -------------- (6-9)
Pc
RT c
b i = Ω B ----------- (6-10)
Pc
1 1
Ω A = --- + --- Ω B ( 4 + 10Ω B ) (6-11)
3 3
Ω B = β – ----------------- – 0.03125
63
(6-12)
1024β
2
α i = [ 1 + ( 0.37464 +1.54226ω i – 0.26992ω i2 ) ( 1 – T r, i ) ] (6-13)
For mixtures
xi xj ( 1 – kPR,i,j )
aα = ( aα ) i ( aα ) j (6-14)
i j
xi bi
b = (6-15)
i
The binary interaction parameters (BIPs), kPR, are symmetric with zeros in the
diagonal:
k PR = k PR,i,j (6-16)
k PR,i,j = 0 (6-17)
When binary interaction parameters are missing in the database for a set of species, the
value is set to zero (a warning node is created). The values for critical temperature, Tc,
critical pressure, Pc, and acentric factor, ωi must be specified for all species.
Peng-Robinson (Twu)
For the Twu modificationRef. 8 of the Peng-Robinson model, the alpha function, αi,
is replaced by
( Ni ( Mi – 1 ) ) NM
αi = T + R exp ( L i ( 1 – T r, ii i ) ) (6-18)
The binary interaction parameters kPR, are used for the Twu modification. Acentric
factor, ωi, is not used in this model but critical temperature and critical pressure must
be specified for all species. The species specific fit parameters Li, Mi, Ni can be
determined by fitting the pure species phase equilibrium to the vapor pressure curve.
Soave-Redlich-Kwong
The classical Soave-Redlich-Kwong equation of state Ref. 9 is given by
RT aα
P = ------------- – ----------------------- (6-19)
V – b V(V + b)
with
2 2.5
R Tc
a i = Ω A ------------------ (6-21)
Pc
RT c
b i = Ω B ----------- (6-22)
Pc
1 1
Ω A = --- + --- Ω B ( 3 + 3Ω B ) (6-23)
3 3
1-
-- --2-
1 3 3 1
Ω B = ------ 2 27 – --- (6-24)
27 3
For mixtures
xi xj ( 1 – kSRK,i,j )
a = ai aj (6-25)
i j
xi bi
b = (6-26)
i
The binary interaction parameters, kSRK, are symmetric with zero in the diagonal:
If a value is missing for kSRK,i,j in the database, it is set to zero (a warning node is
created). The Values for critical temperature, Tc, critical pressure, Pc, and acentric
factor, ωi, must be specified for all species. The Soave-Redlich-Kwong equation of state
is a version of Equation 6-19 modified by Soave Ref. 10, where for pure species i, the
alpha function is modified to
2
α i = [ 1 + ( 0.480 +1.574ω i – 0.176ω i2 ) ( 1 – T r, i ) ] (6-29)
2
α i = [ 1 + ( 0.48508 +1.55174ω i – 0.1561ω i2 ) ( 1 – T r, i ) ] (6-30)
which yields
P i, sat
φ̂ i, l = γ i φ i, sat --------------- F i (6-33)
P
where the activity coefficient, γi, describes the non-ideal liquid phase and φ̂ i, sat P i, sat
is the fugacity at the vapor-liquid phase boundary at equilibrium for the pure species
i. The Poynting correction, Fi, describes the pure species fugacity deviation from the
boiling curve and can be expressed as
P V i, l
F i = exp
P i, sat
---------
RT
dp (6-34)
V i, l
F i ≈ exp --------- ( P – P i, sat ) (6-35)
RT
The Poynting correction can often be ignored for moderate pressure. Hence,
Equation 6-33 can be expressed as
If the vapor phase is considered ideal, then φ i, sat = 1 and the above equation reduces
to
P i, sat
φ̂ i, l = γ i --------------- (6-37)
P
This reduction can be selected explicitly in case the vapor phase is not ideal.
Ideal Solution
For an ideal solution the activity coefficient is equal to one, which gives:
ln γ i = 0 (6-39)
Regular Solution
The Scatchard-Hilderbrand equation Ref. 13 for non-polar mixture is
2
V i ( δ i – δ av )
ln γ i = --------------------------------- (6-40)
RT
where Vi is species molar volume and δi is species solubility parameter, and δav is
( xi Vi δi ) (6-41)
i
δ av = ----------------------------
( xi Vi )
i
The volume parameter, Vi, is set equal the liquid volume, Vi,l,b at normal boiling point
which must be specified for all species. The solubility parameter, δi must be specified
for all species and can be estimated from the normal heat of vaporization, ΔHvap,i and
the liquid volume at normal boiling point as below:
ΔH vap, i
δ i ≡ --------------------
- (6-42)
V i, l, b
2
V i ( δ i – δ av )
ln γ i = --------------------------------- + ln ( θ i ) + 1 – θ i (6-43)
RT
where
Vi
θ i = ----------------------- (6-44)
( xi Vi )
i
Wilson
Wilson Ref. 16 derived his activity coefficient model from a consideration of
probabilities of neighboring molecules in a liquid
xΛ (6-45)
-----------------------
j j, i
ln γ i = 1 –
- – ln
xj Λi, j
j
x k Λ j, k
j
k
V w, j λ i, j
Λ i, j = ----------- exp – -------- (6-46)
V w, i T
λ i, i = 0 (6-47)
Λ i, i = 1 (6-48)
NRTL
Renon and Prausnitz (Ref. 17) formulated a three parameter activity coefficient model
that is able to describe liquid-liquid equilibrium; the non-random two-liquid (NRTL)
model:
A i, j
τ i, j = --------- (6-50)
T
G i, j = exp ( – α i, j τ i, j ) (6-51)
The three parameters are Ai,j, Aj,i and αi,j. A more general form is implemented here:
A i, j
τ i, j = --------- + B i, j (6-52)
T
G i, j = exp ( β i, j – α i, j τ i, j ) (6-53)
The binary interaction parameters, Ai,j, are specified in terms of absolute temperature.
The diagonal values are zero and the matrix is non-symmetric. All off-diagonal values
must be specified.
The binary interaction parameters, Bi,j, have values of zero on the diagonal and the
matrix is non-symmetric. For each pair of species at least Ai,j or Bi,j should be specified.
The randomness parameters, αi,j, have values of zero on the diagonal and the matrix
is symmetric. All off-diagonal values must be specified. Alternatively one can set the
more generic form directly specifying parameter βi,j for which diagonal values are zero
and the matrix is non-symmetric. For each pair of species at least αi,j or βi,j should be
specified.
If any value for these parameters is missing in the database, it is set to zero (warning
node is created).
A i, i = 0 (6-54)
B i, i = 0 (6-55)
α i, i = 0 (6-56)
α j, i = α i, j (6-57)
G i, i = 1 (6-59)
β i, i = 0 (6-60)
UNIQUAC
Abrams and Prausnitz followed up with another two-liquid model known as Universal
Quasi Chemical equation (UNIQUAC) Ref. 18 which is formulated in terms of two
activity coefficients:
The first term is the combinatorial part contributes to the Gibbs free energy
originating from size and shape effects as
φi φi
ln γ i, comb = 1 – φ i + ln φ i – --- q i 1 + ---- + ln ----
z
(6-62)
2 θi θ i
and the second term is the residual part from chemical interactions between the
molecules,
j
x j q j τ j, i
x j q j τ j, i
ln γ i, res = q i 1 – ln -------------------------
- –
--------------------------- (6-63)
x j q j x j q j τ k, j
j
j
k
where
ri
φ i = ----------------- (6-64)
xj rj
j
qi
θ i = ----------------- (6-65)
xj qj
j
– Δ E i, j
τ i, j = exp ----------------- (6-66)
T
ΔE i, i = 0 (6-67)
τ i, i = 1 (6-68)
V VDW, i
r i = -----------------------------------
–3
- (6-69)
0.01517 ×10
A VDW, i
q i = -------------------5- (6-70)
2.5 ×10
UNIFAC
The UNIQUAC Functional-group Activity Coefficients (UNIFAC; see Ref. 19) uses
the same equations as UNIQUAC but the parameters are constructed from group
contributions. The model can be used if UNIQUAC parameters are not available for
all species. The activity coefficients are calculated from Equation 6-61. The
combinatorial part follows from equation Equation 6-62, where
ri = ν k, i r k (6-71)
k
qi = ν k, i q k (6-72)
k
where rk and qk are the values for group k in species i, andνk,i is the number of
occurrences of group k in molecule. The residual term in Equation 6-61 is calculated
from a summation over functional groups:
l
x l q l τ l, k
x l q l τ k, l
ln γ k, res = q k 1 – ln --------------------------
- –
----------------------------- (6-74)
x l q l
x l q l τ m, l
l
j
m
where xl and xm are the compositions of functional group l and m in the mixture
xi νm, i
i
x m = ----------------------------------- (6-75)
x i ν k, i
i k
ν m, i
x m = ----------------- (6-76)
ν k, i
k
A k, m
τ k, m = exp – ------------- (6-77)
T
The default group and interaction parameters are those published by the UNIFAC
consortium (Ref. 20 through Ref. 25), with added groups from Balslev and Abildskov
(Ref. 26) but can be modified per package or database. The groups must be specified
for all species. Note that the interaction parameter matrix is sparse, and a package can
only be used if all interaction parameters for all used groups are specified.
ln φ i, l = ln γ i + ln φ i, l, 0 (6-78)
A i, 1
log φ i, l, 1 = A i, 0 + ---------- + A i, 2 T r, i + A i, 3 T r2, i + A i, 4 T r3, i (6-80)
T r, i
A i, 12
log φ i, l, 2 = A i, 10 + ------------- + A i, 11 T r, i + A i, 13 T r3, i + A i, 14 ( P r, i – 0.6 ) (6-81)
T r, i
The enthalpy, entropy and Gibbs free energy can be calculated from Equation 6-104
to Equation 6-113.
For liquids the density is defined as the reciprocal of the liquid volume:
1
ρ l = ------ (6-82)
Vl
Equation of State
When an equation of state is selected as the thermodynamic model, the liquid volume
is calculated using the equation of state model. In the other cases, when the
thermodynamic model is not an equation of state, an equation of state can be explicitly
selected to calculate the liquid volume.
Ideal Mixture
For an ideal mixture the liquid volume is computed from the pure species densities:
xi
V l, m = --------
ρ i, l
(6-83)
i
The pure species densities are evaluated from the built-in database.
COSTALD
Hankinson and Thomson Ref. 27 presented the Corresponding States Liquid Density
(COSTALD) equation as
V -
------------ = V r, ref ( 1 – ωV r, δ ) (6-84)
V mix
V r, ref = 1 + a ( 1 – T r ) 1 / 3 + b ( 1 – T r ) 2 / 3 + c ( 1 – T r ) + d ( 1 – T r ) 4 / 3 (6-85)
e + fT r + gT r2 + hT r3
V r, δ = ---------------------------------------------------- (6-86)
T r – 1.00001
1
V mix = ---
4 xi Vi + 3 xi Vi2 / 3 xi Vi1 / 3 (6-87)
i i i
2
x i T c, i V i
i
T c = ------------------------------------------ (6-88)
V mix
-----
T
T- T < Tc
Tr = c (6-89)
1 T ≥ Tc
ω = xi ωi (6-90)
i
2
V i = 5.385 V VDW, i – 5.1022VVDW 3
+ 79.524VVDW (6-91)
,i ,i
4 5 6
– 99.316VVDW, i + 100.88VVDW, i – 1152.7VVDW, i
V i = V c, i (6-92)
If the COSTALD acentric factor, ωi is not specified then it can be set equal to the
generic acentric factor for species i. Critical temperature, Tc,i must be specified for all
species.
a -1.52816 e -0.296123
b 1.43907 f 0.386914
c -0.81446 g -0.0427258
d 0.190454 h -0.0480645
Rackett
The Rackett equation Ref. 28 computes the liquid density at the saturation point, and
can be used to describe liquid density at any pressure using the assumption that the
liquid is incompressible. The equation and its condition can be expressed as:
x i T c, i ( 1 + ( 1 – Tr )2 ⁄ 7 )
2
Vl = R
-
---------------------------
( MW ) i P c, i Z r, i
xi ( MW )i (6-93)
i i i
xi Tc, i
Tc = (6-94)
i
-----
T
T- T < Tc
Tr = c (6-95)
1 T ≥ Tc
Z r, i = Z c, i (6-96)
For temperature dependent-functions for example liquid Density, the section appears
as
<Density>
<Phase>Liquid</Phase>
<Coefficients>Tlb;a0;a1;a2;a3;Tub</Coefficients>
<Data>Tlb;f(Tlb);Tub;f(Tub)</Data>
<Comment></Comment>
</Density>
where Tlb and Tub are lower and upper bound for temperature. The a1 to a4 are
coefficients for a cubic polynomial as f(T) = a0 + a1*T+ a2*T2+ a3*T3 fitted to the
experimental data. The unit for temperature is K, See Table 6-3 for other units:
TABLE 6-3: UNIT FOR TEMPERATURE DEPENDENT PROPERTIES
PROPERTY UNIT
3
liquid density mol/m
ideal heat capacity J/mol/K
ln vapor pressure 1
heat of vaporization J/mol
vapor viscosity Pa.s
ln liquid viscosity 1
vapor thermal conductivity W/m/K
liquid thermal conductivity W/m/K
surface tension N/m
Database includes all the references for constant and temperature dependent
properties. In the saved thermodynamic system, it is located between
<Comment></Comment>. Alternatively under the Thermodynamic System node, you
can right click on the created species function created and select Properties.
Ideal Gas
The ideal gas law is independent of composition and determines V at given T and P.
Density can be calculated from
1
ρ = ---- (6-97)
V
ln φ̂ i = 0 (6-98)
T
xi Hi, ig, T T
H ig = + C P, i, ig dT (6-99)
ref
ref
i
where Hi,ig,Tref relates the enthalpy of an ideal gas to the enthalpy at the selected
reference state for species i.
T C P, i, ig
- dT – R ln ---------
P
S ig = x i S i, ig, T ref – R ln x i +
T ref
-----------------
T P ref
- (6-100)
i
where Si,ig,Tref is the entropy of an ideal gas to the entropy of species at the selected
reference state.
G ig = H ig – TS ig (6-101)
Equation of State
The equation of state determines V at given x, T, and P. Density can be expressed as
Equation 6-97. The partial fugacity coefficients are derived from
P
ˆ
0 Vi – --------
1 RT
ln φ̂ i = -------- (6-102)
RT P
∂ ln φ̂ i
H = H ig – RT 2 xi -------------
∂T
- (6-103)
i
∂ ln φ̂ i
S = S ig – R xi ln φ̂i + T -------------
∂T
- (6-104)
i
G = G ig – RT xi ln φ̂i (6-105)
i
Heat Capacity
∂H
C P = -------- (6-106)
∂T P,x
∂H
C v = -------- (6-107)
∂T v,x
The relationship between heat capacity at constant pressure and constant volume can
be expressed as:
∂V ∂P
C p – C v = T ------- ------- (6-108)
∂T P,x ∂T v,x
∂V ∂V
Δv = ------- ΔT + ------- ΔP = 0 (6-109)
∂T P,x ∂P T,x
∂V
T -------
∂T P,x
C v = C P + ------------------------ (6-110)
∂V
-------
∂P T,x
Cp – Cv = R (6-111)
Cp
ϒ = ------- (6-112)
Cv
2 ∂ ln γ i
T
H = xi Hi, ig, T ref
+ T ref
C P, i, ig dT – ΔH i, vap – RT --------------
∂T
(6-113)
i
2 ∂ ln φ̂ i, sat
– R T ------------------------
∂T
T C P, i, ig P i, sat
------------------ dT – R ln --------------
S = xi Si, ig, T ref
– R ln x i + T ref
T P ref
- (6-114)
i
∂ ln γ i ∂ ln φ̂ i, sat ΔH i, vap
+ ln γ i + T -------------- + ln φ̂ i, sat + T ------------------------ + R ---------------------
∂T ∂T T
T
G = x i H i, ig, T ref +
T ref
C P, i, ig dT – T S i, ig, T ref +
(6-115)
i
T C P, i, ig Pi, sat
------------------ dT + RT ln x i + ln γ i + ln φ̂ i, sat + ln --------------
-
T ref T P ref
Note that if the vapor phase is ideal, then the saturated fugacity, φ i, sat , contribution
can be ignored.
Other Properties
Partial fugacity is calculated from
where Ui,ig,ref is the enthalpy of an ideal gas to the species enthalpy at the selected
reference state.
φ̂ i, q
K i, p, q = --------- (6-118)
φ̂ i, p
γ i, q
K i, p, q = --------- (6-119)
γ i, p
xi Mi
M mix = (6-120)
i
At reference conditions, the heat of any reaction relates to the heat of formation as
Note that for reacting flow or a heat balance in a reactor when the heat of reaction is
explicitly taken into account, the enthalpy should not include the heat of formation.
xi Sabs, i
SF = S + (6-123)
i
The entropy balance over a process that includes reactions should include either the
entropy of reaction and use S, or use SF without entropy of reaction.
If user needs the absolute value of entropy of formation, it is possible to estimate it by:
S i = ΔS f, i = νj Sabs, j (6-124)
j
where the entropy of formation of species i is calculated from its elemental constituent
j. For example, entropy of formation of ammonia at 298K is
N 2 ( g ) + 3H 2 ( g ) ⇔ 2NH 3 (6-125)
The reference values for enthalpy, Hi,ref, entropy, Si,ref, and internal energy, Ui,ref, are
calculated so that the pure species enthalpy, entropy, and internal energy are equal to
zero at reference conditions.
The reference values for enthalpy, including formation terms, and entropy, including
absolute terms, are calculated such that the corresponding property for pure species i
has a value equal to the specified formation term at the species reference conditions.
Transport Properties
This section includes a variety of models that are available in the thermodynamic
database for thermal conductivity and viscosity as described below:
Ideal
The thermal conductivity correlations is according to:
x i λ i , v + Δλ v , P
λv = (6-127)
i
The pressure correction Δλv,P is calculated from the method of Stiel and Thodos, see
Ref. 29, which is applicable for ρr<3, but is less accurate for H2, strongly polar gases,
and gases with a high degree of hydrogen bonding, such as H2O and NH3
P c2 / 3
Δλ v, P = ----------------------------- A ( exp ( ( Bρ r ) + C ) ) (6-128)
MT c1 / 6 Z c5
1
ρ c = ------ (6-130)
Vc
xi xj Vc, i, j
Vc = (6-131)
i j
xi ωi
ω = (6-133)
i
RT c Z c
P c = ----------------- (6-135)
Vc
xi Mi
M = (6-136)
i
1
V c, i, j = --- ( V c1, /i3 + V c1,/j3 ) 3 (6-137)
8
T c, i, j = T c, i T c, j (6-138)
A B C
ρr < 0.5 A1=2.702×108 B1=0.535 C1=-1
0.5 ≤ ρr < 2.0 A =2.528×108
2 B2=0.670 C2=-1.069
ρr ≥ 2.0 A3=0.574×108 B3=1.155 C3= 2.016
The vapor thermal conductivity correlation must be available for all species. Also
critical volumes, Vc,i, critical temperatures, Tc,i, molecular weights Mi, and acentric
factors ωi must be specified for all species.
Kinetic Theory
Lindsay and Bromley (see Ref. 31) provided an equation for the interaction parameters
of the method of Wassiljewa (see Ref. 32) based on the kinetic theory, to provide
mixture thermal conductivity from pure species values
xλ (6-141)
------------------------
i i, v
λv =
- + Δλ v, P
i
x j λφ i, j
j
2
T + 3
--- T b, i 9
T + --4- T b, i T b, j
η M 3 / 4 2
= --- 1 + ---------- ------- ------------------------------ ----------------------------------------------
1 i, v j
φ i, j (6-142)
4 η j, v M i
T + --- T b, j T + --- T b, i
3 3
2 2
where the pressure correction Δλv,P is calculated from Equation 6-130. Both vapor
thermal conductivity correlation λi,v and the vapor viscosity correlation ηi,v must be
available for all species. In addition, all normal boiling points Ti,b, molecular weights
Mi, critical volumes Vc,i, critical temperatures Tc,i, and acentric factors ωi must be
specified.
Wilke
Wilke, see Ref. 33, based his method for mixture viscosity of the vapor phase on kinetic
theory:
xη (6-143)
----------------------
,
i i v
ηv =
i
x j ψ i, j
j
2
η i, v M j 1 / 4
1 + ---------- --------
η j, v Mi
ψ i, j = ------------------------------------------------------ (6-144)
Mi
8 1 + -------
M j
The vapor viscosity correlation ηi,v must be available for all species. In addition, all
molecular weights Mi must be specified.
Brokaw
Brokaw (see Ref. 34) uses the same basic equation as Wilke (Equation 6-143).
However, Equation 6-144 is replaced by
η i, v
ψ i, j = S i, j A i, j ---------
- (6-145)
η j, v
β i, j RM i, j – RM i0.45 ,j
A i, j = -------------------- 1 + ----------------------------------------------------------------------------------------- (6-146)
RM i, j 1 + RM i0.45
2 ( 1 + RM i, j ) + ------------------------------------------------
,j
( β i, j ) ( 1 + RM i, j )
where
4M i M j 1 / 4 Mi
β i, j = ----------------------------2 , RM i, j = ------- (6-147)
( Mi + M ) Mj
j
The vapor viscosity correlation, ηi,v must be available for all species i. In addition, all
molecular weights Mi must be specified. If Lennard-Jones energy εi (see Ref. 35)
T 2- + 1 --- δ δ
1 + --------
ε i ε j 4 s, i s, j
S i, j = ---------------------------------------------------------------------------- (6-148)
T 1 T 1
1 + ---- + --- δ s, i 1 + ---- + --- δ s, j
εi 4 εj 4
Otherwise,
S i, j = 1 (6-149)
M t P c2 / 3 10 – 7
- --------------------------- ζ
Δη v, P = ----------------------- (6-150)
T c1 / 6 101325 2 / 3
where ξ is calculated from the correlation of Jossi (Ref. 38), which is applicable for ρr<
3.0. It is less accurate for H2, strongly polar gases and gases with a high degree of
hydrogen bonding such as H2O and NH3.
The correction factor is due to using pressure, atm, and viscosity, cP, units in Jossi’s
correlation. It is expressed as:
1
---
4
(ξ + 1) = 1.0230 + 0.23364ρ r + 0.58533 ρ r2 (6-151)
– 0.40758 ρ r3 + 0.093324 ρ r4
ρ
ρ r = ----- (6-152)
ρc
1
ρ c = ------ (6-153)
Vc
x i Z c, i xi Mi
Zc = ,M = (6-155)
t
i i
The values for critical volumes, Vc,i, critical temperatures, Tc,i, critical compressibility
factors, Zc,i and molecular weights, Mi must be specified for all species i.
The high pressure correction is available for the Wilke and Brokaw mixture models.
The vapor viscosity follows from
η v = η v, Wilke + Δη v, P (6-156)
The CH4 viscosity is calculated from Ref. 41, modified by Pedersen and Fredunsland
(Ref. 42) to avoid issues below 91 K where CH4 becomes solid
η CH = η CH 4, 0 + ρ CH 4 η CH , 1 + F 1 η CH , 2 + F 2 η CH , 3 (6-157)
4 4 4 4
where
H+1
F 1 = -------------- (6-158)
2
1–H
F 1 = -------------- (6-159)
2
exp ( ΔT ) – exp ( – ΔT )
H = ------------------------------------------------------------- (6-160)
exp ( ( ΔT ) + exp ( – ΔT ) )
ΔT = T – 91 (T ∈ K) (6-161)
Here, ρCH4 is used in g/cm3; for the mass-mole conversion of ρCH4, a molecular
weight of MCH4 = 16.042568 g/mol is used.
The first density correction for the moderately dense gas is given by
T 2
η CH4, 1 = L 10 – L 11 1.4 – ln ---------- (6-163)
168
L 12 L 15
η CH4, 2 = exp --------- – L 13 exp ( 10 ρ CH4 ) L 14 – ----------- + (6-164)
T T 3 / 2
ρ CH – ρ c, CH4 L 17 L 18
4 - ( ρ CH ) L 16 + --------
------------------------------------ - + --------- – 1
ρ c, CH4 4 T T 2
L 19 L 22
η CH4, 3 = exp --------- – L 20 exp ( 10 ρ CH 4 ) L 21 – -----------
T T 3 / 2 (6-165)
ρ CH 4 – ρ c, CH 4 L L
+ ------------------------------------ ( ρ CH 4 ) L 23 + --------- + --------
25
24
- – 1
ρ c, CH T T2
4
with the values of the parameters L1 through L25 are listed in Table 6-5 below:
TABLE 6-5: METHANE VISCOSITY NUMERICAL COEFFICIENTS.
Here, ρCH4 is used in g/cm3; the critical density is given by ρc,CH4 = 0.16284 g/cm3.
The following equation by McCarty (Ref. 43) is solved for the density of CH4
N 15 N 16 ρ CH7
4 N 18 N 19 ρ CH9
4
+ ρ CH
6 – ---------- + ---------
- + N 17 ------------- + ρ CH
8 ---------- + ---------
- + N 20 -------------
4
T T
2
T 4
T T
2
T
3 N 22 N 23 N 24 N 25 N 26 N 27
+ exp ( – N 21 ρ CH
2 ) ρ CH ---------
2
- – ---------- + ρ CH
3
5 ---------
2
- + ---------- + ρ CH
4
7 ---------
- + ----------
4
4
T T
4
T T
4
T2 T
3
N 28 N 29 N 30 N 31 N 32 N 33 4
+ ρ CH
9 – ---------
2
- – ---------- + ρ CH
4
11 – ---------
2
- + ---------- + ρ CH
3
13 ---------
2
- – ---------- + N 34 T
3
4
T T
4
T T
4
T T
With the viscosity and density of CH4 defined, the viscosity of any mixture, ηm, can be
calculated from the corresponding states principle
–1
Tc ------ Pc 2 / 3 Mm 1 / 2 α
η m = ------------------ 6 ----------------- --------------- ------------- η CH4, P 0, T 0 (6-167)
T c, CH P c, CH M CH α CH
4 4 4 4
P c, CH α CH
P 0 = -----------------4 P -------------4 (6-169)
Pc α
The following mixing rules are used for the critical properties, see Ref. 44:
x i x j β i, j T c, i T c, j
(6-170)
i j
-----------------------------------------------------------
-
Tc =
x i x j β i, j
i j
8 x i x j β i, j T c, i T c, j
(6-171)
i j
--------------------------------------------------------------
-
Pc =
2
xi xj βi, j
i j
T c, i 1 / 3 T c, j 1 / 3 3
β i, j = ---------- + ---------- (6-172)
P c, i P c, j
The parameter α is
–3 0.5173
α = 1 + 7.378 ×10 ρ r1.847 M m (6-173)
ρ
CH 4, P 0, T 0
ρ r = --------------------------
- (6-174)
ρ c, CH4
where
T c, CH
T 0 = ------------------4 T (6-175)
Tc
P c, CH
P 0 = -----------------4 P (6-176)
Pc
–3 0.5173
α CH 4 = 1 + 7.378 ×10 ρ r1.847 M CH 4 (6-177)
–4 2.303 2.303
M m = 1.304 ×10 ( M W – MN ) + MN (6-178)
xi Mi
2
i (6-179)
M W = --------------------
-
MN
xi Mi
MN = (6-180)
i
Note that pure species vapor viscosity correlations ηi,v are not required. However, for
each species i, molecular weight Mi, critical temperature, Tc,i and critical pressure, Pc,i
must be specified.
Ideal
To calculate the mixture liquid thermal conductivity, λl,m, the values of pure liquid
thermal conductivity correlations are mixed ideally
x i λ i, l
λ l, m = (6-181)
i
The pressure dependence is based on the work of Missenard (Ref. 45) where
λ l, m , P
----------------- = 1 + QP r0.7 (6-182)
λ l, m
where Q is correlated as
T
T r = ------ (6-184)
Tc
RT c Z c
P c = ----------------- (6-186)
Vc
x i V c, i
Vc = (6-187)
i
x i T c, i
Tc = (6-188)
i
x i Z c, i
Zc = (6-189)
i
All liquid thermal conductivity correlations must be specified. All values for critical
temperatures, Tc,i, critical volumes, Pc,i and critical compressibility factors, Zc,i must
be specified for all species i.
Power Law
The values of pure liquid vapor thermal conductivity correlations are mixed according
to the following power law
1 xi
-----------
λ l2, m
= --------
λ l, i
(6-190)
i
All liquid thermal conductivity correlations must be specified. The model is valid for
pure compound thermal conductivity values that are no more apart than a factor of 2
(Ref. 46 - Ref. 47). The pressure dependence is introduced using Equation 6-182 –
Equation 6-189. All values for critical temperatures, Tc,i, critical volumes, Pc,i and
critical compressibility factors, Zc,i must be specified for all spices i.
Local Composition
The local composition model by Rowley (Ref. 47) uses an ideal and excess
contribution
xi Mi
ω i = -------------------- (6-193)
xj Mj
j
where Gj,i follows from Equation 6-51. The binary interaction terms follow from
ω i ω i, i λ i, l + ω j ω j, j λ j, l
λ i, j = λ j, i = -------------------------------------------------------- (6-195)
ω i ω i, i + ω j ω j, j
λ i, i = λ i, l (6-196)
and
ω i, i = ω i ( ω i + ω j G j, i ) (6-197)
with ϖi is the composition in the binary mixture of species i and j and the local
composition is equi-molar
M i G j, i
ω i = -------------------------------------------------- (6-198)
M i G j, i + M j G i, j
M i ( ω i ω i, i )λ i, l + M j ( ω j ω j, j )λ j, l
λ i, j = λ j, i = ---------------------------------------------------------------------------------- (6-199)
M i ( ω i ω i, i ) + M j ( ω j ω j, j )
which he found to produce better model predictions in most cases where both the
Local Composition model and Power Law model have trouble. However, the model
is not as generally applicable; for instance, systems containing H2O are not well
described by this model due to the low molecular weight of H2O.
LIQUID VISCOSITY
The following mixture models are available for liquid viscosity.
xi ln ηi, l
ln η l, m = (6-200)
i
ωi ln ηi, l
ln η l, m = (6-201)
i
Ideal
The gas-liquid surface tension is predicted by ideally mixing the pure species
correlations. It is independent of pressure, vapor temperature or composition.
xi, l σi, vl
σ vl = (6-202)
i
where the vapor-liquid surface tension correlations, σi,vl, must be specified for all
species i, and are evaluated at the temperature of the liquid phase.
Winterfeld
Following Winterfeld (Ref. 48), the vapor-liquid surface tension is predicted by mixing
the pure species correlations according to
x i, l x j, l σ i, vl σ j, vl
---------------------------------------------
ρ i, l ρ j, l
-
i j (6-203)
σ vl = -------------------------------------------------------------
-
x i, l 2
--------
ρ i, l
i
References
3. A. Akerberg, CFD analyses of the gas flow inside the vessel of a hot isostatic press,
Master of Science Thesis, KTH School of Industrial Engineering and Management,
Stockholm, Sweden, 2012.
4. M.A. Trebble and P.R. Bishnoi, “Accuracy and consistency comparisons of ten cubic
equations of state for polar and non-polar compounds”, Fluid Phase Equilibria,
29:465–474, 1986.
5. B.E. Poling, J.M. Prausnitz, and J.P. O’Connell. The Properties of Gases and
Liquids. McGraw-Hill, international edition, 2007.
8. C.H. Twu, J.E. Coon, and J.R. Cunningham. “A new generalized alpha function
for a cubic equation of state, Part 1. Peng-Robinson equation”. Fluid Phase
Equilibria, 105:49–59, 1995.
11. M.S. Graboski and T.E. Daubert. “A modified Soave equation of state for phase
equilibrium calculations. 3. systems containing hydrogen”. Industrial &
Engineering Chemistry Process Design and Development, 18(2):300–306, 1979.
12. K.C. Chao and J.D. Seader. “A general correlation of vapor-liquid equilibria in
hydrocarbon mixtures”. AIChE Journal, 7(4):598–605, 1961.
13. J.H. Hildebrand and R.L. Scott. “The solubility of non-electrolytes”. Journal of
Physical Chemistry, 55(4):619–620, 1951.
14. H.G. Grayson and C.W. Streed. “Vapor-liquid equilibria for high temperature,
high pressure hydrogen-hydrocarbon systems”. 6th World Petroleum Congress, 19–
26 June, Frankfurt am Main, Germany, IV:169–180, 1963.
15. J.H. Hildebrand, J.M. Prausnitz, and R.L. Scott. Regular and Related Solutions.
Van Nostrand Reinhold Co., 1970.
16. G.M. Wilson. “Vapor-liquid equilibrium. xi. a new expression for the excess free
energy of mixing”. Journal of the American Chemical Society, 86:127–130, 1964.
18. D.S. Abrams and J.M. Prausnitz. “Statistical thermodynamics of liquid mixtures:
A new expression for the excess Gibbs energy of partly or completely miscible
systems”. AIChE Journal, 21(1):116–128, 1975.
26. K. Balslev and J. Abildskov. “Unifac parameters for four new groups”. Industrial
& Engineering Chemistry Research, 41:2047–2057, 2002.
27. R.W. Hankinson and G.H. Thomson. “A new correlation for saturated densities
of liquid and their mixtures”. AIChE Journal, 25(4):653–663, 1979.
28. H.G. Rackett. “Equation of state for saturated liquids”. Journal of Chemical and
Engineering Data, 15(4):514–517, 1970.
29. L.I. Stiel and G. Thodos. “The thermal conductivity of nonpolar substances in the
dense gaseous and liquid regions”. AIChE Journal, 10(1):26–30, 1964.
31. A.L. Lindsay and L.A. Bromley. “Thermal conductivity of gas mixtures”.
Industrial and Engineering Chemistry, 42(8):1508–1511, 1950.
33. C.R. Wilke. “A viscosity equation for gas mixtures”. The Journal of Chemical
Physics, 18(4):517– 520, 1950.
34. R.S. Brokaw. “Approximate formulas for the viscosity and thermal conductivity of
gas mixtures. ii”.The Journal of Chemical Physics, 42(4):1140–1147, 1965.
35. J.E. Lennard-Jones. “On the determination of molecular fields 1. from the
variation of the viscosity of a gas with temperature”. Proceedings of the Royal Society
of London. Series A, Containing Papers of a Mathematical and Physical Character,
106:441–462, 1924.
36. W.H. Stockmayer. “Second virial coefficients of polar gases”. The Journal of
Chemical Physics, 9:398–402, 1941.
37. F.M. Mourits and F.H.A. Rummens. “A critical evaluation of Lennard-Jones and
Stockmayer potential parameters and of some correlation methods”. Canadian
Journal of Chemical Engineering, 55:3007–3020, 1977.
38. J.A. Jossi, L.I. Stiel, and G. Thodos. “The viscosity of pure substances in the dense
gaseous and liquid phases”. AIChE Journal, 8(1):59–63, 1962.
39. K.S. Pedersen and P.L. Christensen. “Phase Behavior of Petroleum Reservoir
Fluids”. CRC Press/Taylor & Francis Group, 2007.
41. H.J.M. Hanley, W.M. Haynes, and R.D. McCarty. “The viscosity and thermal
conductivity coefficients for dense gaseous and liquid methane”. Journal of Physical
Chemistry, 6(2):597–609, 1977.
42. K.S. Pedersen and A. Fredenslund. “An improved corresponding states model for
the prediction of oil and gas viscosities and thermal conductivities”. Chemical
Engineering Science, 42(1):182–186, 1987.
44. S. Murad and K.E. Gubbins. “Corresponding states correlation for thermal
conductivity of dense fluids”. Chemical Engineering Science, 32(5):499–505, 1977.
46. R.L. Rowley, G.L. White, and M. Chiu. “Ternary liquid mixture thermal
conductivities”. Chemical Engineering Science, 43(2):361–371, 1988.
47. R.L. Rowley. “A local composition model for multicomponent liquid mixture
thermal conductivities”. Chemical Engineering Science, 37(6):897–904, 1982.
48. P.H. Winterfield, L.E. Scriven, and H.T. Davis. “An approximate theory of
interfacial tension of multicomponent systems: Applications binary liquid-vapor
tensions”. AIChE Journal, 24(6):1010–1014, 1978.
50. E.C. Carlson, “Don’t Gamble With Physical Properties For Simulation”,
Chemical Engineering Progress, 92, 10, pp. 35–46, 1996.
Glossary
381
Glossary of Terms
absorption Uptake of a gas into the bulk of a liquid. Gas absorption takes place for
example in the liquid of a scrubber tower where an up-streaming gas is washed by a
down-going flow of a scrubber solution.
Arrhenius rate equation Expression that relates the rate constant of a chemical
reaction to the exponential of the temperature.
batch reactor Reactor characterized by its operation, which means that the reactor
does reaches steady state.
bipolar plate Electrically conducting plate connected to the anode on one side and to
the cathode on the other side in an electrochemical cell.
Butler-Volmer equation Expression that relates the reaction rate of an electron transfer
reaction on an electrode surface to the exponential of the overpotential. The equation
can be derived from the Arrhenius rate equation by accounting for the contribution of
the electric potential to the activation energy.
boundary layer Region in a fluid close to a solid surface. This region is characterized
by large gradients in velocity and is often treated with approximate methods because
it is difficult to geometrically resolve the large gradients found there.
continuous reactor Reactor that operates without interruption. This type of reactor is
characterized by its steady-state operation.
Darcy’s law Equation that gives the velocity vector as proportional to the pressure
gradient. Often used to describe flow in porous media.
electroneutrality condition Condition that states that the sum of charges in a control
volume in an electrolyte should be zero.
Euler flow Flow at high velocities, where incompressibility of the fluid is of importance
whereas the influence of viscous momentum transport is negligible.
Fick’s laws The first law relates the concentration gradients to the diffusive flux of a
solute infinitely diluted in a solvent. The second law introduces the first law into a
differential material balance for the solute.
fully developed laminar flow Laminar flow along a channel or pipe that only has
velocity components in the main direction of the flow. The velocity profile
perpendicular to the flow does not change downstream in the flow.
heterogeneous reaction Reaction that takes place at the interface between two phases.
Maxwell-Stefan equations Set of equations that describe the diffusion of solutes and
solvent in a concentrated solution. In such a solution, the solutes interact with each
other and with the solvent.
monolithic reactor Catalytic reactor made of one single piece of solid material.
Incorporates a catalytic structure in its often porous structure.
Nernst-Planck equation Equation that describes the flux of an ion through diffusion,
convection, and migration in an electric field. The equation is valid for diluted
electrolytes.
Poiseuille’s law Equation that relates the mass rate of flow in a tube as proportional to
the pressure difference per unit length and to the fourth power of the tube radius. The
law is valid for fully developed laminar flow.
specific surface area Internal surface area of a porous structure given in area per unit
volume, which yields the unit one over unit length. Often used to characterize the
structure of porous catalysts.
surface molar flux the tangential flux in the surface dimension as governed by
diffusion according to Fick’s law.
wall function Semi-empirical expression for the anisotropic flow close to a solid surface
used in turbulence models. Often based on negligible variations in pressure gradient
in the direction tangential to the surface.
28 documentation 20
INDEX| 385
transport of diluted species 184 frequency factors 35
Dusty gas model 225 Freundlich exponent 138
386 | I N D E X
M mass action law 34 pair nodes
mass action law, laminar flow 177 electrophoretic transport interface
mass based concentrations (node) 189 268
mass fraction (node) 242 surface reactions 277
mass fractions 117 transport of concentrated species 229
mass transport 150 transport of diluted species 184
Maxwell-Stefan diffusion model 225 parameter estimation, solver 110
Maxwell-Stefan diffusivity matrix 154 partially saturated porous media (node)
mixture-averaged diffusion model 155, 205
226 periodic condition (node)
monolayer adsorption 167 transport of diluted species 195
MPH-files 21 persistence
multicomponent diffusion 153 of thermodynamic packages 313
multiphase flash calculations 309 physics interfaces 15
INDEX| 387
reaction thermodynamics (node) 80, 108 space dimensions 15
Reactions species (node) 71, 102
in porous catalyst pellets 211 species activity (node) 80, 109
reactions (node) species group (node) 78, 108
Nernst-Planck equations 255 species names, syntax for 31
reacting flow, laminar flow 239 species source (node) 273
surface reactions 279 species thermodynamics (node) 80, 109
transport of concentrated species 240 species types 31
transport of diluted species 189 standard settings 19
reactive pellet bed 140 Stiel-Thodos equation 52
reactive pellet bed (node) 208 Stockmayer potential 74, 103
reactor types stoichiometric coefficients 71, 101
batch (constant volume) 44 stratified porous media 136
batch, defined 59 superficial volume averages, porous me-
batch, theory 43 dia 132
CSTR (constant volume), defined 59 Supporting Electrolytes 128
CSTR (constant volume), theory 45 supporting electrolytes
plug flow, definition 59 theory, Nernst-Planck equations 172
plug flow, theory 48 surface concentration (node) 280
semibatch, definition 59 surface equilibrium reaction (node) 199,
semibatch, theory 47 247
removing species 117 surface properties (node) 278
retardation factor 138 surface reactions interface 276
reversible reaction 27, 67, 98 theory 166
reversible reaction group (node) 75, 105 symmetry (node)
RF Module 117 Nernst-Planck equations 257
transport of concentrated species 245
S selecting
transport of diluted species 193
chemical species transport interfaces
113 T Tafel equation 256
space dimensions and physics interfac- technical support, COMSOL 22
es 15 temperature exponents 35
semibatch reactors 47, 59 theory
solute-solvent approximations 31 chemistry interface 34
solvents, association factor of 51 Nernst-Planck equations 172
solver settings, parameter estimation reaction engineering 34
110 surface reactions 166
solver time, discretization and 63 transport of concentrated species in-
Soret effect 158 terface 150
388 | I N D E X
transport of diluted species in porous
media interface 121
transport of diluted species interface
120
thermal conductivity 52
thermal diffusion 158
thermodynamics package 308
adding 309
thickness
fracture 213
thin diffusion barrier (node) 197
Thin Impermeable Barrier 197
tortuosity factors 134
transport mechanisms 230, 234, 237
transport of concentrated species inter-
face 223
theory 150
transport of diluted species in porous
media interface 183
theory 121
transport of diluted species interface
179, 215
theory 120
turbulent mixing (node)
transport of concentrated species 238
transport of diluted species 188
U uniform scaling 63
V variable scope 29
W websites, COMSOL 22
Wilke-Chang equation 50
INDEX| 389
390 | I N D E X