Micro Fluidics Module Users Guide
Micro Fluidics Module Users Guide
Users Guide
Contact Information
Visit the Contact COMSOL page at [Link]/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 [Link]/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 [Link]/support/case.
Other useful links include:
Support Center: [Link]/support
Product Download: [Link]/product-download
Product Updates: [Link]/support/updates
Discussion Forum: [Link]/community
Events: [Link]/events
COMSOL Video Gallery: [Link]/video
Support Knowledge Base: [Link]/support/knowledgebase
Part number: CM021901
C o n t e n t s
Chapter 1: Introduction
About the Microfluidics Module
14
About Microfluidics . . . . . . . . . . . . . . . . . . . . . . 14
About the Microfluidics Module . . . . . . . . . . . . . . . . . 15
The Microfluidics Module Physics Interface Guide . . . . . . . . . . 15
Coupling to Other Physics Interfaces . . . . . . . . . . . . . . . 18
The Microfluidics Module Study Capabilities by Physics Interface . . . . . 19
Common Physics Interface and Feature Settings and Nodes. . . . . . . 20
The Liquids and Gases Materials Database . . . . . . . . . . . . . 21
Where Do I Access the Documentation and Application Libraries? . . . . 21
Overview of the Users Guide
25
28
31
36
44
CONTENTS
|3
Heat Transfer . . . . . . . . . . . . . . . . . . . . . . . . 53
Coupling to Other Physics Interfaces . . . . . . . . . . . . . . . 54
Modeling Rarefied Gas Flows
56
60
82
4 | CONTENTS
100
102
104
106
108
108
109
C h a p t e r 4 : M u l t i p h a s e F l ow , Two - P h a s e F l ow
Interfaces
The Laminar Two-Phase Flow, Level Set and Laminar
Two-Phase Flow, Phase Field Interfaces
The Laminar Two-Phase Flow, Level Set Interface
112
. . . . . . . . .
112
114
Domain, Boundary, Point, and Pair Nodes for the Laminar and
Turbulent Flow, Two-Phase, Level Set and Phase Field Interfaces. . .
115
Wall. . . . . . . . . . . . . . . . . . . . . . . . . . .
116
Fluid Properties . . . . . . . . . . . . . . . . . . . . . .
118
Gravity
. . . . . . . . . . . . . . . . . . . . . . . . .
120
Initial Values. . . . . . . . . . . . . . . . . . . . . . . .
121
Initial Interface . . . . . . . . . . . . . . . . . . . . . . .
121
122
Domain, Boundary, Edge, Point, and Pair Nodes for the Laminar
CONTENTS
|5
124
Free Deformation . . . . . . . . . . . . . . . . . . . . .
126
126
Initial Values. . . . . . . . . . . . . . . . . . . . . . . .
127
Fixed Mesh . . . . . . . . . . . . . . . . . . . . . . . .
127
127
Fluid-Fluid Interface . . . . . . . . . . . . . . . . . . . . .
128
129
Wall-Fluid Interface . . . . . . . . . . . . . . . . . . . . .
130
Navier Slip . . . . . . . . . . . . . . . . . . . . . . . .
130
132
132
135
Phase Initialization . . . . . . . . . . . . . . . . . . . . .
135
Numerical Stabilization
. . . . . . . . . . . . . . . . . . .
136
137
138
138
139
140
144
145
6 | CONTENTS
148
. . . .
149
149
Initial Values. . . . . . . . . . . . . . . . . . . . . . . .
150
Inlet . . . . . . . . . . . . . . . . . . . . . . . . . . .
150
Initial Interface . . . . . . . . . . . . . . . . . . . . . . .
151
No Flow . . . . . . . . . . . . . . . . . . . . . . . . .
151
152
Domain, Boundary, and Pair Nodes for the Phase Field Interface . . . .
153
153
Initial Values. . . . . . . . . . . . . . . . . . . . . . . .
155
Inlet . . . . . . . . . . . . . . . . . . . . . . . . . . .
155
Initial Interface . . . . . . . . . . . . . . . . . . . . . . .
155
Wetted Wall . . . . . . . . . . . . . . . . . . . . . . .
156
157
157
159
159
160
161
162
162
162
. . . . . . . . . . .
164
164
165
. . . . . . . . . . . . . . . . . .
165
166
C h a p t e r 6 : Po r o u s M e d i a F l ow I n t e r f a c e s
The Darcys Law Interface
168
Domain, Boundary, Edge, Point, and Pair Nodes for the Darcys
Law Interface . . . . . . . . . . . . . . . . . . . . . .
169
171
Mass Source
. . . . . . . . . . . . . . . . . . . . . . .
172
Initial Values. . . . . . . . . . . . . . . . . . . . . . . .
172
Pressure . . . . . . . . . . . . . . . . . . . . . . . .
172
Mass Flux. . . . . . . . . . . . . . . . . . . . . . . . .
173
Inlet . . . . . . . . . . . . . . . . . . . . . . . . . . .
173
Symmetry . . . . . . . . . . . . . . . . . . . . . . . .
173
CONTENTS
|7
No Flow . . . . . . . . . . . . . . . . . . . . . . . . .
174
Flux Discontinuity . . . . . . . . . . . . . . . . . . . . .
174
Outlet . . . . . . . . . . . . . . . . . . . . . . . . . .
174
176
178
179
Forchheimer Drag . . . . . . . . . . . . . . . . . . . . .
179
Mass Source
. . . . . . . . . . . . . . . . . . . . . . .
180
Volume Force . . . . . . . . . . . . . . . . . . . . . . .
180
Initial Values. . . . . . . . . . . . . . . . . . . . . . . .
181
Fluid Properties . . . . . . . . . . . . . . . . . . . . . .
181
183
Domain, Boundary, Point, and Pair Nodes for the Free and Porous
Media Flow Interface . . . . . . . . . . . . . . . . . . .
185
186
Volume Force . . . . . . . . . . . . . . . . . . . . . . .
186
Forchheimer Drag . . . . . . . . . . . . . . . . . . . . .
187
Initial Values. . . . . . . . . . . . . . . . . . . . . . . .
187
187
189
189
. . . . . . . . . . . . . .
8 | CONTENTS
184
Fluid Properties . . . . . . . . . . . . . . . . . . . . . .
191
. . . . . . . . . . . . . . . .
191
192
193
194
194
196
Domain, Boundary, Edge, Point, and Pair Nodes for the Slip Flow
Interface. . . . . . . . . . . . . . . . . . . . . . . .
Fluid
198
. . . . . . . . . . . . . . . . . . . . . . . . . .
199
200
Initial Values. . . . . . . . . . . . . . . . . . . . . . . .
201
Slip Wall . . . . . . . . . . . . . . . . . . . . . . . . .
201
Periodic Condition . . . . . . . . . . . . . . . . . . . . .
202
Symmetry . . . . . . . . . . . . . . . . . . . . . . . .
202
Continuity . . . . . . . . . . . . . . . . . . . . . . . .
202
203
203
205
C h a p t e r 8 : C h e m i c a l S p e c i e s Tr a n s p o r t I n t e r f a c e s
The Transport of Diluted Species Interface
208
211
213
Transport Properties . . . . . . . . . . . . . . . . . . . .
214
Turbulent Mixing . . . . . . . . . . . . . . . . . . . . . .
216
Initial Values. . . . . . . . . . . . . . . . . . . . . . . .
217
Mass-Based Concentrations . . . . . . . . . . . . . . . . . .
217
Reactions. . . . . . . . . . . . . . . . . . . . . . . . .
218
No Flux . . . . . . . . . . . . . . . . . . . . . . . . .
219
Inflow . . . . . . . . . . . . . . . . . . . . . . . . . .
219
Outflow . . . . . . . . . . . . . . . . . . . . . . . . .
220
Concentration . . . . . . . . . . . . . . . . . . . . . . .
220
Flux . . . . . . . . . . . . . . . . . . . . . . . . . . .
221
Symmetry . . . . . . . . . . . . . . . . . . . . . . . .
221
Flux Discontinuity . . . . . . . . . . . . . . . . . . . . .
221
CONTENTS
|9
Periodic Condition . . . . . . . . . . . . . . . . . . . . .
222
223
Open Boundary . . . . . . . . . . . . . . . . . . . . . .
224
224
225
Equilibrium Reaction . . . . . . . . . . . . . . . . . . . .
225
226
226
Reaction Coefficients . . . . . . . . . . . . . . . . . . . .
227
227
227
229
Volatilization . . . . . . . . . . . . . . . . . . . . . . .
230
231
Species Source. . . . . . . . . . . . . . . . . . . . . . .
233
Hygroscopic Swelling . . . . . . . . . . . . . . . . . . . .
233
235
236
237
239
239
. . . . . . . . . . . . . . .
240
242
Supporting Electrolytes . . . . . . . . . . . . . . . . . . .
243
Crosswind Diffusion
10 | C O N T E N T S
222
. . . . . . . . . . . . . . . . . . . .
244
245
245
Convection . . . . . . . . . . . . . . . . . . . . . . . .
247
248
Diffusion . . . . . . . . . . . . . . . . . . . . . . . . .
249
Dispersion . . . . . . . . . . . . . . . . . . . . . . . .
250
Adsorption . . . . . . . . . . . . . . . . . . . . . . . .
251
Reactions. . . . . . . . . . . . . . . . . . . . . . . . .
252
253
References . . . . . . . . . . . . . . . . . . . . . . . .
258
Chapter 9: Glossary
Glossary of Terms
260
CONTENTS
| 11
12 | C O N T E N T S
Introduction
This guide describes the Microfluidics Module, an optional add-on package that
extends the COMSOL Multiphysics modeling environment with customized
physics interfaces for microfluidics.
This chapter introduces you to the capabilities of this 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.
In this chapter:
About the Microfluidics Module
Overview of the Users Guide
13
About Microfluidics
The field of microfluidics evolved as engineers and scientists explored new avenues to
exploit the fabrication technologies developed by the microelectronics industry. These
technologies enabled complex micron and submicron structures to be integrated with
electronic systems and batch fabricated at low cost. Mechanical devices fabricated using
these technologies have become known as microelectromechanical systems (MEMS),
whilst fluidic devices are commonly referred to as microfluidic systems or
lab-on-a-chip devices. A proper description of these microsystems usually requires
multiple physical effects to be incorporated.
At the microscale different physical effects become important to those dominant at
macroscopic scales. Properties that scale with the volume of the system (such as inertia)
become comparatively less important than those that scale with the surface area of the
system (such as viscosity and surface tension). Fluid flow is therefore usually laminar
and chemical migration is often limited by diffusion. Electrokinetic effects become
important as the electric double layers present at interfaces in the system interact with
external applied fields. As systems are further miniaturized, the mean free path of the
fluid can become comparable to the size of the system and rarefaction effects become
important. At moderate Knudsen numbers (the Knudsen number is the ratio of the
mean free path to the system size), it is still possible to use the Navier Stokes equations
to solve the flow, however special slip boundary conditions are required.
14 |
CHAPTER 1: INTRODUCTION
15
Microfluidics Module are listed in the table. The functionality of the COMSOL
Multiphysics base package is given in the COMSOL Multiphysics Reference Manual.
In the COMSOL Multiphysics Reference Manual:
Studies and Solvers
The Physics Interfaces
For a list of all the core physics interfaces included with a COMSOL
Multiphysics license, see Physics Interface Guide.
PHYSICS INTERFACE
ICON
TAG
SPACE
DIMENSION
tds
all dimensions
Transport of Diluted
Species in Porous Media
tds
all dimensions
Laminar Flow1
spf
3D, 2D, 2D
axisymmetric
Creeping Flow
spf
3D, 2D, 2D
axisymmetric
tpf
3D, 2D, 2D
axisymmetric
Fluid Flow
Single-Phase Flow
Multiphase Flow
Two-Phase Flow, Level Set
Laminar Two-Phase
Flow, Level Set
16 |
CHAPTER 1: INTRODUCTION
PHYSICS INTERFACE
ICON
TAG
SPACE
DIMENSION
tpf
3D, 2D, 2D
axisymmetric
3D, 2D, 2D
axisymmetric
time dependent
tpfmm
br
3D, 2D, 2D
axisymmetric
Darcys Law
dl
all dimensions
fp
3D, 2D, 2D
axisymmetric
slpf
3D, 2D, 2D
axisymmetric
Level Set
ls
all dimensions
Phase Field
pf
all dimensions
Rarefied Flow
Slip Flow
Mathematics
Moving Interface
1 This physics interface is included with the core COMSOL package but has added
17
Pressure-driven flow
Chemical Reactions in two-phase flow
Dielectrophoresis
Electroosmotic flow
Electrophoresis
Electrothermal flow
Electrowetting
Magnetophoresis
Mass transport using diffusion, migration, and convection
Slip Flow
DEVICES
Lab-on-a-chip devices
Microfluidic channels
Microreactors
Micromixers
MEMS heat exchangers
Nonmechanical pumps and valves
18 |
CHAPTER 1: INTRODUCTION
When using the axisymmetric interfaces, the horizontal axis represents the
r direction and the vertical axis the z direction. The geometry must be
created in the right half-plane (that is, only for positive r).
DEPENDENT
VARIABLES
PRESET
STUDIES*
tds
tds
slpf
u, v, w, p, T
Laminar Flow
spf
u, v, w, p
Creeping Flow
spf
u, v, w, p
tpf
u, v, w, p,
TIME DEPENDENT
NAME
STATIONARY
PHYSICS INTERFACE
TABLE 1-2: MICROFLUIDICS MODULE DEPENDENT VARIABLES AND PRESET STUDY AVAILABILITY
Slip Flow
FLUID FLOW>SINGLE-PHASE FLOW
19
TABLE 1-2: MICROFLUIDICS MODULE DEPENDENT VARIABLES AND PRESET STUDY AVAILABILITY
tpf
u, v, w, p, ,
tpfmm
u, v, w, p
PRESET
STUDIES*
TRANSIENT WITH INITIALIZATION
DEPENDENT
VARIABLES
TIME DEPENDENT
NAME
STATIONARY
PHYSICS INTERFACE
Brinkman Equations
br
u, v, w, p
Darcys Law
dl
u, v, w, p
fp
u, v, w, p
Level Set
ls
Phase Field
pf
MATHEMATICS
20 |
CHAPTER 1: INTRODUCTION
The COMSOL Multiphysics Reference Manual describes the core physics interfaces
and functionality included with the COMSOL Multiphysics license. This book also has
instructions about how to use COMSOL and how to access the electronic
Documentation and Help content.
21
displays information about that feature (or click a node in the Model Builder followed
by the Help button (
). This is called topic-based (or context) help.
To open the Help window:
In the Model Builder, Application Builder, or Physics Builder click a node or
window and then press F1.
On any toolbar (for example, Home, Definitions, or Geometry), hover the
mouse over a button (for example, Add Physics or Build All) and then
press F1.
From the File menu, click Help (
).
) button.
).
22 |
CHAPTER 1: INTRODUCTION
) button.
):
) Applications
23
24 |
COMSOL website
[Link]
Contact COMSOL
[Link]/contact
Support Center
[Link]/support
Product Download
[Link]/product-download
Product Updates
[Link]/support/updates
Discussion Forum
[Link]/community
Events
[Link]/events
[Link]/video
[Link]/support/knowledgebase
CHAPTER 1: INTRODUCTION
To help you navigate through this guide, see the Contents, Glossary, and Index.
MODELING IN MICROFLUIDICS
The Microfluidic Modeling chapter familiarizes you with modeling procedures useful
when working with this module. Topics include Dimensionless Numbers in
Microfluidics, Modeling Microfluidic Fluid Flows, Modeling Coupled Phenomena in
Microfluidics, and Modeling Rarefied Gas Flows.
THE FLUID FLOW BRANCH INTERFACES
There are several fluid flow interfaces available. The various types of momentum
transport that you can simulate includes laminar and creeping flow, multiphase
two-phase flow, and flow in porous media. Every section describes the applicable
physics interfaces in detail and concludes with the underlying interface theory.
Single-Phase Flow
Single-Phase Flow Interfaces chapter describes the Laminar Flow and Creeping Flow
interfaces.
O V E R V I E W O F T H E U S E R S G U I D E
25
Rarefied Flow
Rarefied Flow Interface chapter describes the Slip Flow interface, including the
underlying theory.
T H E C H E M I C A L S P E C I E S TR A N S P O R T B R A N C H I N T E R F A C E S
26 |
CHAPTER 1: INTRODUCTION
Microfluidic Modeling
This chapter gives an overview of the physics interfaces available for modeling
microfluidic flows and provides guidance on choosing the appropriate physics
interface for a specific problem.
In this chapter:
Physics and Scaling in Microfluidics
Dimensionless Numbers in Microfluidics
Modeling Microfluidic Fluid Flows
Modeling Coupled Phenomena in Microfluidics
Modeling Rarefied Gas Flows
27
28 |
Volume
LENGTH
SCALING
L3
2
Area
Inertial Forces
L3
Viscous Forces
Laplace Pressure
L1
CONSEQUENCES
Capillary Force
Permeability of Porous
media
L2
L2
L1
L1
Laminar flows make mixing particularly difficult, so mass transport is often diffusion
limited. The diffusion time scales as L2, but even in microfluidic systems diffusion is
often a slow process. This has implications for chemical transport and hence reactions
within microfluidic systems. Chemical Transport and Reactions describes how to
model diffusion-based transport and chemical reactions in microfluidics.
As the length scale of the flow becomes comparable to intermolecular length scale
more complex kinetic effects become important. For gases the ratio of the molecular
mean free path to the flow geometry size is given by the Knudsen number (Kn).
Clearly, Kn scales as 1/L. For Kn < 0.01 fluid flow is usually well described by the
Navier-Stokes equations with no-slip boundary conditions. In the slip flow regime
(0.01 < Kn < 0.1) appropriate slip boundary conditions can be used with the
Navier-Stokes equations to describe the flow away from the boundary (see Slip Flow).
At Knudsen numbers above 0.1 a fully kinetic approach is required. Flows in this
regime can be modeled using tools from the Molecular Flow Module.
29
The Free Molecular Flow and Transitional Flow interfaces are available in
the Molecular Flow Module.
30 |
Dimensionless Numbers in
M i c r o f lui di c s
Information about the dominant physics in a microfluidics problem is contained in
many of the dimensionless numbers that are used to characterize the flow. When using
the finite element method, dimensionless numbers defined on the element or cell
level can contain important information about the numerical stability of the problem
(the term cell is carried over from the finite volume method in this context). This
section includes information about the dimensionless numbers that are relevant to
microfluidic flows.
In this section:
Dimensionless Numbers Important for Solver Stability
Other Dimensionless Numbers
31
number as an example oscillations can occur when the cell Peclet number is greater
than one in the following circumstances:
A Dirichlet boundary condition can lead to a solution containing a steep gradient
near the boundary, forming a boundary layer. If the mesh cannot resolve the
boundary layer, this creates a local disturbance.
A space-dependent initial condition that the mesh does not resolve can cause a local
initial disturbance that propagates through the computational domain.
A small initial diffusion term close to a nonconstant source term or a nonconstant
Dirichlet boundary condition can result in a local disturbance.
In theory the grid can be refined to bring the cell Reynolds or Peclet number below
one, although this is often impractical for many problems. Several stabilization
techniques are included, which enable problems with larger cell Reynolds or Peclet
numbers to be solved. At the crudest level additional numerical diffusion can be added
to the problem to improve its stability. This is achieved by selecting Isotropic Diffusion
under Inconsistent Stabilization for any physics interface.
) is clicked
32 |
TABLE 2-2: IMPORTANT DIMENSIONLESS NUMBERS RELATED TO THE STABILITY OF VARIOUS MICROFLUIDIC
PROBLEMS.
DIMENSIONLESS
NUMBER
SYMBOL
GLOBAL
DEFINITION
CELL
DEFINITION
NOTES
Reynolds
number
Re
vL
Re = ----------
c
vh
Re = ---------2
[Link]
Re < 10
Re > 1
Peclet number
(Heat Flow)
Pe
2
= ( c p ) .
[Link]
Stabilization required:
c
Pe > 1
Peclet number
(Mass
Transport)
Pe
vL
Pe = ------D
c
vh
Pe = -------2D
Ratio of concentration
convection to diffusion.
Stabilization required:
c
Pe > 1
Symbol definitions: v is the characteristic velocity or cell velocity, r is the fluid density, L is the characteristic
length scale for the problem, is the fluid viscosity, h is the element size, cp is the heat capacity at constant
pressure, is the thermal conductivity, and D is the diffusion constant. For anisotropic diffusion or conductivity an appropriate average is computed. The variables names assume the default name spf.
33
SYMBOL
Bond number/
Etvs number
Bo
Capillary
number
Ca
Dukhin number
Du
Knudsen
number
Kn
Kn = ---L
Mach number
Ma
v
Ma = --a
Eo
DEFINITION
NOTES
2
Bo
gL
= ------------
Eo
v
Ca = -----
K
Du = ------B
K
Symbol definitions: is the fluid density, g is the body force acceleration (usually the acceleration
due to gravity), L is the characteristic length scale for the problem,
is the surface tension coef
B
ficient, is the fluid viscosity, v is the characteristic velocity, K is the surface conductivity, K is
the bulk conductivity, is the mean free path, is the thermal diffusivity (=/(cp)), cp is the
heat capacity at constant pressure, is the thermal conductivity, and T is the absolute temperature
(T is the characteristic temperature difference).
34 |
SYMBOL
DEFINITION
Marangoni
number
Mg
Ohnesorge
/Laplace
/Surataman
numbers
Oh
Oh = --------------L
La
NOTES
Su
La
1
= ----------2Su
Oh
Weber number
We
v L
We = ------------
Symbol definitions: is the fluid density, g is the body force acceleration (usually the acceleration
due to gravity), L is the characteristic length scale for the problem, is the surface tension coef
B
ficient, is the fluid viscosity, v is the characteristic velocity, K is the surface conductivity, K is
the bulk conductivity, is the mean free path, is the thermal diffusivity (=/(cp)), cp is the
heat capacity at constant pressure, is the thermal conductivity, and T is the absolute temperature
(T is the characteristic temperature difference).
35
36 |
Single-Phase Flow
The Fluid Flow>Single-Phase Flow branch (
) when adding a physics interface includes
the Laminar and Creeping Flow interfaces.
The Laminar Flow Interface (
) is used primarily to model slow-moving flow in
environments without sudden changes in geometry, material distribution, or
temperature. The Navier-Stokes equations are solved without a turbulence model.
laminar flow typically occurs at Reynolds numbers less than 1000. By default the flow
is compressible (see Figure 2-1).
The Creeping Flow Interface (
) uses the same equations as the Laminar Flow
interface with the additional assumption that the contribution of the inertia term is
negligible. This is often referred to as Stokes flow and is appropriate for use when
viscous flow is dominant, which is often the case in microfluidics applications.
Creeping flow applies when the Reynolds number is much less than one. A creeping
flow problem is significantly simpler to solve than a laminar flow problem so it is best
to make this assumption explicitly if it applies. By default the flow is compressible (see
below).
By selecting the Neglect Inertial Form (Stokes Flow) check-box, on the
Settings window for Laminar Flow quickly convert a Laminar Flow
interface into a Creeping Flow interface.
For both physics interfaces several additional options are available.
SHALLOW CHANNEL APPROXIMATION
Often you might want to simplify long, narrow channels by modeling them in 2D. The
Use Shallow Channel Approximation option is useful as it includes a drag term to
approximate the added affects given by thinness of the gap between one set of
boundaries in comparison to the others.
COMPRESSIBLE FLOW
Compressible flow (for speeds of less than Mach 0.3) is the default option for both
physics interfaces. For compressible flow it is important that the density and any mass
balances are well defined throughout the domain. Choosing to model incompressible
flow simplifies the equations to be solved and decreases solution times. Most gas flows
should be modeled as compressible flows; however, liquid flows can usually be treated
as incompressible.
37
NON-NEWTONIAN FLOW
The physics interfaces also allow for easy definition of non-Newtonian fluid flow
through access to the dynamic viscosity in the Navier-Stokes [Link] fluid can
be modeled using the power law and Carreau models or by means of any expression
that describes the dynamic viscosity appropriately.
Both physics interfaces also include a feature to compute a laminar velocity profile at
arbitrarily shaped inlets and outlets, which makes models much easier to set up.
Figure 2-1: The Settings window for Laminar Flow. Model compressible or
non-compressible flow, or laminar and Stokes flow. Combinations are also possible.
Multiphase Flow
The Multiphase Flow branch (
) enables the modeling of multiphase flows. These
physics interfaces are included:
The Laminar Two-Phase Flow, Level Set Interface (
Two-Phase Flow, Level Set branch (
)).
38 |
The Two-Phase Flow interfaces add surface tension forces (including the Marangoni
effect) automatically at the two fluid interface(s). A library of surface tension
coefficients between some common substances is available.
For problems involving topological changes (for example, jet breakup), use either the
Level Set or Phase Field interfaces. These techniques use an auxiliary function (the
level set and phase field functions, respectively) to track the location of the interface,
which is necessarily diffuse. The Level Set interface usually produces a more accurate
representation of surface tension forces, and is recommended for use in smaller scale
problems with lower velocities, where the surface tension is dominant. The phase field
method is physically motivated and is usually more numerically stable than the level set
method. It is can also be extended to more phases and is compatible with
fluid-structure interactions (requires the MEMS or the Structural Mechanics Module).
Switch between the Two-Phase Flow, Level Set interface and the Two-Phase
Flow, Phase Field interface by selecting one or other from a list available in
both interfaces. This is useful if you are not sure which provides the best
solution.
The moving mesh method represents the interface as a boundary condition along a
line or surface in the geometry. Because the physical thickness of phase boundaries is
usually very small, for most practical meshes The Laminar Two-Phase Flow, Moving
Mesh Interface describes the two-phase boundary the most accurately. However, it
cannot accommodate topological changes in the boundary.
For all the Two-Phase Flow interfaces, compressible flow is possible to model at speeds
of less than 0.3 Mach. You can also choose to model incompressible flow by simplifying
the equations to be solved. Stokes law is also an option.
In each physics interface, the density and viscosity are specified for both fluids. You can
easily use non-Newtonian models for any of the two fluids, based on the power law,
the Carreau model, or using an arbitrary user-defined expression.
It is often advantageous to use more than one of these techniques to solve a problem
for example, a level set model for jet breakup could be checked prior to breakup by a
39
moving mesh model to ensure that the surface tension is captured accurately by the
diffuse interface.
) when adding an
Under the Mathematics>Moving Interface branch (
interface, The Level Set Interface (
) and The Phase Field
Interface (
) are available in a form uncoupled with the equations of
flow. These interfaces can be used to model phenomena in which other
factors dominate over the fluid flow, such as some forms of phase
separation.
Similarly, under the Mathematics>Moving Mesh branch ( ), Moving
Mesh Interface ( ) is available (and described in the COMSOL
Multiphysics Reference Manual) to treat a range of other problems
involving moving meshes, such as the transport of solid particles in a fluid
domain.
40 |
equations are similar in form to the Navier Stokes equations. Figure 2-2 shows the
Settings window for the Brinkman Equations interface.
For low Reynolds number flows in which other terms in the Brinkman equations are
necessary the Neglect Inertial Term (Stokes-Brinkman Flow) feature can be selected
to neglect the inertial term in the equations (this is selected by default). The flow is
treated as incompressible by default, but compressible flow can be enabled by selecting
Compressible Flow (Ma<0.3). When using the compressible flow feature the fluid density
must be defined as a function of the local pressure.
Figure 2-2: The Settings window for Brinkman Equations. Model compressible or
non-compressible flow and Stokes flow. Combinations are also possible.
FREE A N D PO ROUS MEDIA
41
It should be noted that if the porous medium is large in comparison to the free
channel, and the results in the region of the interface are not of interest, then it is
possible to manually couple a Fluid Flow interface to the Darcys Law interface. This
makes the model computationally cheaper.
NAME
COMPRESSIBILITY
Laminar Flow
spf
Compressible
flow (Ma<0.3)
None
Creeping Flow
spf
Compressible
flow (Ma<0.3)
Stokes Flow
42 |
PHYSICS INTERFACE
LABEL
NAME
MULTIPHASE
FLOW MODEL
COMPRESSIBILITY
NEGLECT
INERTIAL TERM
(STOKES FLOW)
Laminar, Two-Phase
Flow, Level Set
tpf
Incompressible
flow
None
Laminar, Two-Phase
Flow, Phase Field
tpf
Incompressible
flow
None
Laminar, Two-Phase
Flow, Moving Mesh
tpfmm
not applicable
Compressible
flow (Ma<0.3)
None
NAME
COMPRESSIBILITY
NEGLECT INERTIAL
TERM
PORE SIZE
Darcy's Law
dl
n/a
n/a
Brinkman
Equations
br
Incompressible
flow
Yes Stokes-Brinkman
High permeability
and porosity, faster
flow
fp
Incompressible
flow
Not selected
High permeability
and porosity, fast
flow
43
44 |
of their concentration, the concentration of other species and other model parameters,
such as the temperature.
Figure 2-3: The Settings window for Transport of Diluted Species, with Convection selected
as the Transport Mechanism by default.
S T A B I L I Z A T I O N S E T T I N G S F O R D I L U T E D S P E C I E S TR A N S P O R T
For some laminar flow problems it can be useful to change the settings for the dilute
) and select
species transport stabilization. To do this, click the Show button (
Stabilization. A Consistent stabilization section is now visible in the Transport of Diluted
species settings, and in some cases it can be desirable to change the settings for the
crosswind diffusion. By default the Crosswind diffusion type is set to Do Carmo and
Galeo. This type of crosswind diffusion reduces undershoot and overshoot to a
minimum but can in rare cases give equations systems that are difficult to fully
converge. The alternative option, Codina is less diffusive and so should be used if the
species transport is highly convective, or if convergence problems occur. This option
can result in more undershoot and overshoot and is also less effective for anisotropic
meshes. The Codina option activates a text field for the Lower gradient limit glim
45
tensor components are neglected. This setting is usually accurate enough and is faster
to compute. If required, select Full residual instead.
Electrohydrodynamics
Electrohydrodynamics is a general term describing phenomena that involve the
interaction between solid surfaces, ionic solutions, and applied electric and magnetic
fields. Electrohydrodynamics is frequently employed in microfluidic devices to
manipulate fluids and move particles for sample handling and chemical separation.
Electrokinetics refers to a range of fluid flow phenomena involving electric fields.
These include electroosmosis, electrothermal effects, electrophoresis and
dielectrophoresis. Electrosmosis describes the motion of fluids induced by the forces
on the charged EDLs at the surfaces of the fluid. Electrothermal effects occur in a
conductive fluid where the temperature is modified by Joule heating from an AC
electric field. This creates variations in conductivity and permittivity and thus Coulomb
and dielectric body forces. Electrophoresis and dielectrophoresis describe the motion of
charged and polarized particles in a non-uniform AC or DC applied field.
Magnetohydrodynamics refers to fluid flow phenomena involving magnetic fields.
Magnetophoresis is the motion of diamagnetic particles in a nonuniform magnetic field
and is commonly used for magnetic bead separation.
Table 2-7 summarizes these categories. Although these examples describe specific
multiphysics couplings, COMSOL Multiphysics is not limited to these casesfor
example, it is possible to include both magnetic and electric fields to simulate
electromagnetophoresis.
TABLE 2-7: ELECTROHYDRODYNAMIC PHENOMENA
TYPE OF FIELD/FORCE
DC
AC
46 |
Electroosmosis
AC electroosmosis
Electrophoresis /
Dielectrophoresis
AC Electrophoresis
/AC
Dielectrophoresis
Electrothermal
DC
AC
Magnetophoresis
Magnetic Field
Force on suspended particles
When a polar liquid (such as water) and a solid surface (such as glass or a
polymer-based substrate) come into contact, charge transfer occurs between the
surface and the electrolytic solution. At finite temperature the charges on the surface
are not screened perfectly by the ions in the liquid and a finite thickness electric double
layer (EDL) or Debye layer develops. Electroosmosis is the process by which motion
is induced in a liquid due to the body force acting on the EDL in an electric field.
A complete model of the system includes the space charge layer explicitly. The electric
potential is the solution of a nonlinear partial differential equation, the
Poisson-Boltzmann equation, which can be solved by coupling an Electrostatics
interface to a Transport of Diluted Species interface (with migration enabled) for the
ion species. The software then computes the forces on the fluid and a further coupled
Laminar Flow or Creeping Flow interface is used to compute the overall fluid flow. In
practice this approach is only possible for nanoscale channelsas typically the EDL
thickness is 1-10 nm.
The Poisson Boltzmann equation is sometimes linearizedthis is referred to as the
Debye-Hckel approximation which applies when ze k B T , where is the potential
at the surface of the moving volume of fluid (the zeta potential). At room temperature
this corresponds to the limit 26 mV. Note that there is a layer of immobile ions
trapped adjacent to the surface (the Stern Layer which is of order one hydrated ion
radius thick) with an associated volume of immobile fluid; this means that is not
simply the wall potential. is usually determined experimentally from electroosmotic
flow measurements.
The Poisson-Boltzmann equation has a characteristic length scalethe Debye length,
D:
47
D =
k B T
-------------------2 2
2z e c
u eo = ----- E
where E is the applied electric field and is the liquids dynamic viscosity. From this
equation, the electroosmotic mobility is naturally defined as:
eo = ----
(2-1)
The Laminar Flow and Creeping Flow interfaces include a wall boundary condition
option for an Electroosmotic Velocity boundary condition. This enables you to specify
the external electric field (which can be manually specified or coupled from an
Electrostatics interface) and the electroosmotic mobility to define an electroosmotic
flow.
Further details of the theory of electroosmosis can be found in Ref. 1 and Ref. 2.
AC ELECTROOSMOSIS
Because an alternating electric field does not generate a net force on the EDL, AC
electroosmosis is not used for fluid transport in microfluidics. However, the
back-and-forth movements an AC field generates are useful for mixing purposes. To
model AC electroosmotic flow when the frequency of the electric field is sufficiently
low, the same approaches can be taken as for DC electroosmosis. An Electrostatics
interface should still be employed to calculate the electric field, but when this is
coupled into other physics interfaces the AC dependence should be explicitly added
using an expression. With increasing frequency, AC electroosmosis becomes less
important so this approach is valid for most practical examples.
48 |
ELECTROPHORESIS
49
DEP
( t ) = 2 m a K ( m , p ) ( E ( t ) E ( t ) )
(2-2)
where m and p are the complex permittivities of the medium and the particle
respectively, a is the radius of the particles equivalent homogeneous sphere and
K(m,p) is the Clausius-Mossotti function. The complex permittivity, *, for an
isotropic homogeneous dielectric is
= i ---
where is the electric permittivity, is the electrical conductivity, and is the angular
field frequency. The Clausius-Mossotti function is given by:
p m
K ( m, p ) = --------------------- p + 2 m
which depends on the particles complex permittivity, p, and that of the medium, m.
Taking the time average of Equation 2-2 gives the time averaged force:
F
DEP
( t ) = 2 m r 0 Re ( K ( m, p ) ) ( E rms E rms )
By balancing this force with the Stokes drag force the dielectrophoretic velocity is
obtained:
u
DEP
m r 0 Re ( K ( m, p ) ) ( E rms E rms )
= --------------------------------------------------------------------------------------------3
50 |
for the dot product of the electric field with itself. Finally the dielectrophoretic velocity
term can be added to the velocity field to compute the particle velocity field.
MAGNETOPHORESIS
MAP
m r 0 K ( m, p ) ( H ext H ext )
= ------------------------------------------------------------------------------3
51
The gradient of the field can then be computed term by term from the spatial
derivatives of an expression for the dot product of the magnetic field with itself. Finally
the magnetophoretic velocity term can be added to the velocity field to compute the
particle velocity field.
ELECTROTHERMALLY-DRIVEN FLOW
2
1 + ( )
where is the conductivity, is the fluids permittivity, is the angular frequency of
the electric field, and = / is the fluids charge-relaxation time. The electric-field
vector E contains the amplitude and direction of the AC electric field but not its
instantaneous value.
Because of the heating, and are temperature dependent, and their gradients are
functions of the temperature gradient: = (/T)T and = (/T)T. With
water, for example, the relative change rates for the permittivity and the conductivity
are (1/)(/T) = 0.004 1/K and (1/)(/T) = 0.02 1/K, respectively.
ELECTROWETTING
The contact angle of a two-fluid interface with a solid surface is determined by the
balance of the forces at the contact point. The equilibrium contact angle, 0, is given
by Youngs equation:
s1 + 12 cos 0 = s2
Here s1 is the surface energy per unit area between fluid 1 and the solid surface, s2
is the surface energy per unit area between fluid 2 and the solid surface, and 12 is the
surface tension at the interface between the two fluids.
In electrowetting the balance of forces at the contact point is modified by the
application of a voltage between a conducting fluid and the solid surface. For many
52 |
applications the solid surface consists of a thin dielectric deposited onto a conducting
layer; this is often referred to as Electrowetting on Dielectric (EWOD). In this case,
the capacitance of the dielectric layer dominates over the double layer capacitance at
the solid-liquid interface (Ref. 3). The energy stored in the capacitor formed between
the conducting liquid and the conducting layer in the solid reduces the effective
surface energy of the liquid to which the voltage is applied. For the case when a voltage
difference occurs between fluid 1 and the conductor beyond the dielectric Youngs
equation is modified as follows:
2
V
s1 ---------- + cos ew = s2
2d f
12
Here is the permittivity of the dielectric, V is the potential difference applied, and df
is the dielectric thickness. This equation can be rewritten as
2
V
cos ew = cos 0 + -----------------2 12 d f
(2-3)
Heat Transfer
It is often necessary to consider the effect of heat flow and temperature dependent
material properties (such as the viscosity or the surface tension) in microfluidic systems.
The COMSOL Multiphysics base package includes physics interfaces for heat transfer
in incompressible fluids and solids. The Heat Transfer Module is required in addition
to this module to model heat flow in compressible fluids and to include viscous heating
terms. This section describes how to couple heat transfer to microfluidics models using
the functionality available in the COMSOL Multiphysics base package and the
Microfluidics Module.
53
The Heat Transfer interface includes the equations for heat transfer in both fluid and
solid domains. To model a system consisting of both solid and liquid domains add a
single Heat Transfer interface and then add nodes for Heat Transfer in Solids and Heat
Transfer in Fluids with selections corresponding to the solid and fluid domains
respectively. For a fluid domain the convective flow is coupled into domain by selecting
the velocity field from a corresponding Laminar Flow or Creeping Flow interface.
The fluid domain models the heat transfer including both convection and conduction
according to the equation:
T
C p ------- + ( u )T = ( q ) + Q
t
where, is the fluid density (SI unit: kg/m3), Cp is the specific heat capacity at
constant pressure (SI unit: J/(kgK)), T is absolute temperature (SI unit: K), u is the
velocity vector (SI unit: m/s), q is the heat flux by conduction (SI unit: W/m2), and
Q contains heat sources other than viscous heating (SI unit: W/m3). For a solid
domain the following equation applies:
T
C p ------- = ( q ) + Q
t
When using the heat transfer in solids and in fluids domain properties the software
automatically includes the heat transfer between the solid and the fluid. It is important
to include the correct thermal boundary conditions on regions where fluid flows into
or out of a domain, typically a Temperature condition for an inlet and an Outflow
condition for an outlet.
For a detailed discussion of the fundamentals of heat transfer, see Ref. 4.
More details of the heat transfer capabilities can be found in The Heat
Transfer Interfaces in the COMSOL Multiphysics Reference Manual.
54 |
55
56 |
Continuum
flow
Slip flow
Transitional
flow
Free molecular
flow
Figure 2-4: A plot showing the main fluid flow regimes for rarefied gas flows. Different
regimes are separated by lines of constant Knudsen numbers. The number density of the
gas is normalized to the number density of an ideal gas at a pressure of 1 atmosphere and
a temperature of 0 C (n0).
Slip Flow
In the slip flow regime, the Navier-Stokes equations can still be used with modified
boundary conditions to account for rarefaction effects close to the wall. A layer of
rarefied gas with a size similar to the mean free path develops close to the wallthis is
termed the Knudsen layer. The Navier-Stokes equations are not applicable in this
layer, but the flow outside the layer can be described by extrapolating the bulk gas flow
towards the wall and applying Maxwells slip boundary condition at the wall (Ref. 5).
For thermal flows the von Smoluchowski temperature jump boundary condition must
also be applied (Ref. 5).
To model flows in the slip flow regime use the Slip Flow interface. That physics
interface can be used to model isothermal and non-isothermal flows with or without
explicit modeling of thermal processes in adjacent solid domains. To model isothermal
flows use the Temperature boundary on all external model surfaces, with the
temperature set equal to the fluid temperature. The Slip Wall and External Slip Wall
boundary conditions are used to model slip on interior and exterior model boundaries
respectively. These boundary conditions include thermal creep or transpiration, viscous
57
slip, and the von Smoluchowski temperature jump. The slip coefficients can be
specified directly or by means of Maxwells model. For further details see The Slip Flow
Interface and Theory for the Slip Flow Interface.
The Temperature boundary condition is described for the Heat
Transfer interface in the COMSOL Multiphysics Reference Manual.
Slip Velocity
58 |
59
60 |
When the Laminar Flow interface is added, the following default nodes are also added
in the Model BuilderFluid Properties, Wall (the default boundary condition is No slip)
and Initial Values. Other nodes, that implement, for example, boundary conditions and
volume forces, can be added from the Physics toolbar or from the context menu
displayed when right-clicking Laminar Flow.
SETTINGS
Compressibility
Fully compressible flow can be simulated by selecting the Compressible flow (Ma<0.3)
option.
61
(3-1)
where is the fluids dynamic viscosity, u is the velocity field, and dz is the channel
thickness. This term represents the resistance that the parallel boundaries impose on
the flow; however, it does not account for any changes in velocity due to variations in
the cross-sectional area of the channel.
The following dependent variables (fields) are defined for this physics interfacethe
Velocity field u and its components, and the Pressure p.
If required, the names of the field, component, and dependent variable may be edited.
Editing the name of a scalar dependent variable changes both its field name and the
dependent variable name. If a new field name coincides with the name of another field
of the same type, the fields share degrees of freedom and dependent variable names. A
new field name must not coincide with the name of a field of another type, or with a
component name belonging to some other field. Component names must be unique
within a model except when two fields share a common field name.
62 |
ADVANCED SETTINGS
By default, the Neglect inertial term (Stokes flow) check box is selected. If unchecked,
the inertial terms are included in the computations.
63
DISCRETIZATION
By default, the Creeping Flow interface uses P2+P1 elements. Contrary to general
laminar and turbulent single-phase flow simulations employing purely linear P1+P1
elements, P2+P1 elements are well suited for Creeping flow simulations.
CONSISTENT STABILIZATION
This check box is selected by default and should remain selected for optimal
performance. The consistent stabilization method does not perturb the original
transport equation.
The Laminar Flow Interface
Theory for the Single-Phase Flow Interfaces
64 |
Flow Continuity
Pipe Connection1
Fluid Properties
Initial Values
Inlet
Symmetry
Volume Force
Open Boundary
Wall
Outlet
1 A feature that may require an additional license
Fluid Properties
The Fluid Properties node adds the momentum and continuity equations solved by the
physics interface, except for volume forces which are added by the Volume Force
feature. The node also provides an interface for defining the material properties of the
fluid.
MODEL INPUTS
Fluid properties, such as density and viscosity, can be defined through user inputs,
variables or by selecting a material. For the latter option, additional inputs, for example
temperature and/or pressure, may be required to define these properties.
65
Temperature
By default, the single-phase flow interfaces are set to model isothermal flow. If a Heat
Transfer interface is included in the component, the temperature field may
alternatively be selected from this physics interface. All physics interfaces have their
own tags (Name). For example, if a Heat Transfer in Fluids interface is included in the
component, the Temperature (ht) option is available for T.
Absolute Pressure
This input appears when a material requires the absolute pressure as a model input.
The absolute pressure is used to evaluate material properties, but it also relates to the
value of the calculated pressure field. There are generally two ways to calculate the
pressure when describing fluid flow: either to solve for the absolute pressure or for a
pressure (often denoted gauge pressure) that relates to the absolute pressure through
a reference pressure.
The choice of pressure variable depends on the system of equations being solved. For
example, in a unidirectional incompressible flow problem, the pressure drop over the
modeled domain is probably many orders of magnitude smaller than the atmospheric
pressure, which, when included, may reduce the stability and convergence properties
of the solver. In other cases, such as when the pressure is part of an expression for the
gas volume or the diffusion coefficients, it may be more convenient to solve for the
absolute pressure.
The default Absolute pressure pA is p+pref where p is the dependent pressure variable
from the Navier-Stokes or RANS equations, and pref is from the user input defined at
the physics interface level. When pref is non zero, the physics interface solves for a
gauge pressure. If the pressure field instead is an absolute pressure field, pref should be
set to 0.
The Absolute pressure field can be edited by clicking Make All Model Inputs Editable
(
) and entering the desired value in the input field.
FLUID PROPERTIES
Density
If density variations with respect to pressure are to be included in the computations,
the flow must be set to compressible (at the physics interface level).
66 |
Dynamic Viscosity
The Dynamic viscosity describes the relationship between the shear rate and the shear
stresses in a fluid. Intuitively, water and air have low viscosities, and substances often
described as thick (such as oil) have higher viscosities.
Using the built-in variable for the shear rate magnitude, [Link], makes it possible to
define arbitrary expressions of the dynamic viscosity as a function of the shear rate.
For laminar flow, the Non-Newtonian power law may be used to model the viscosity of
a non-Newtonian fluid. The following model parameters are required for the
Non-Newtonian power law:
Fluid consistency coefficient m
Flow behavior index n
Volume Force
The Volume Force node specifies the volume force F on the right-hand side of the
momentum equation, and may, for example, be used to incorporate the effects of
gravity in a component.
u
2
T
+ ( u )u = pI + ( u + ( u ) ) --- ( u )I + F
3
t
67
If several volume-force nodes are added to the same domain, then the sum of all
contributions are added to the momentum equation.
Initial Values
The initial values serve as initial conditions for a transient simulation or as an initial
guess for a nonlinear solver in a stationary simulation. Note that for a transient
compressible-flow simulation employing a material for which the density depends on
the pressure (such as air), discontinuities in the initial values trigger pressure waves
even when the Mach number is small. The pressure waves must be resolved and this
puts a restriction on the time step.
INITIAL VALUES
Initial values or expressions should be specified for the Velocity field u and the Pressure
p.
Wall
The Wall node includes a set of boundary conditions describing fluid-flow conditions
at stationary, moving and leaking walls.
BOUNDARY CONDITION
Leaking Wall
Slip
Electroosmotic Velocity
Sliding Wall
Slip Velocity
Moving Wall
1 The default
No Slip
No slip is the default boundary condition for a stationary solid wall for laminar flow.
The condition prescribes u = 0, that is, the fluid at the wall is not moving.
68 |
Slip
The Slip option prescribes a no-penetration condition, un=0. It is implicitly assumed
that there are no viscous effects at the slip wall and hence, no boundary layer develops.
From a modeling point of view, this can be a reasonable approximation if the main
effect of the wall is to prevent fluid from leaving the domain.
Sliding Wall
The Sliding wall boundary condition is appropriate if the wall behaves like a conveyor
belt; that is, the surface is sliding in its tangential direction. A velocity is prescribed at
the wall and the boundary itself does not have to actually move relative to the reference
frame.
For 3D components, values or expressions for the Velocity of sliding wall uw should
be specified. If the velocity vector entered is not in the plane of the wall, COMSOL
Multiphysics projects it onto the tangential direction. Its magnitude is adjusted to
be the same as the magnitude of the vector entered.
For 2D components, the tangential direction is unambiguously defined by the
direction of the boundary. For this reason, the sliding wall boundary condition has
different definitions in different space dimensions. A single entry for the Velocity of
the tangentially moving wall Uw should be specified in 2D.
Moving Wall
For an arbitrary wall movement, the condition u = uw may be prescribed. In this case,
the components of the Velocity of moving wall uw should be specified.
Specifying this boundary condition does not automatically cause the associated wall to
move. An additional Moving Mesh interface needs to be added to physically track the
wall movement in the spatial reference frame.
Leaking Wall
This boundary condition may be used to simulate a wall where fluid is leaking into or
leaving the domain with the velocity u = ul through a perforated wall. The
components of the Fluid velocity ul on the leaking wall should be specified.
Electroosmotic Velocity
When an electric field drives a flow along the boundary, the components for the Electric
field E along with the Electroosmotic mobility eo should be defined. The Built-in
expression for the Electroosmotic mobility requires values or expressions for the Zeta
potential and the Relative permittivity r.
69
Slip Velocity
In the microscale range, the flow condition at a boundary is seldom strictly no slip or
slip. Instead, the boundary condition is something in between, and there is a Slip
velocity at the boundary. Two phenomena account for this velocity: noncontinuum
effects and the flow induced by a thermal gradient along the boundary. The
components of Velocity of moving wall: uw should be specified. Zero values are used for
a stationary wall.
When the Use viscous slip check box is selected, the default Slip length Ls is User defined.
Another value or expression may be entered if the default value is not applicable. For
Maxwells model values or expressions for the Tangential momentum accommodation
coefficient av and the Mean free path should be specified. Tangential accommodation
coefficients are typically in the range of 0.85 to 1.0 and can be found in G. Kariadakis,
A. Beskok, and N. Aluru, Microflows and Nanoflows, Springer Science and Business
Media, 2005.
When the Use thermal creep check box is selected, a thermal creep contribution with
Thermal slip coefficient T is activated. Thermal slip coefficients are typically between
0.3 and 1.0 and can be found in G. Kariadakis, A. Beskok, and N. Aluru, Microflows
and Nanoflows, Springer Science and Business Media, 2005.
CONSTRAINT SETTINGS
Physics Options.
Inlet
This condition should be used on boundaries for which there is a net flow into the
domain. To obtain a numerically well posed problem, it is advisable to also consider
the Outlet condition(s) when specifying an Inlet condition. For example, if the
pressure is specified at the outlet, the velocity may be specified at the inlet, and vice
versa. Specifying the velocity vector at both the inlet and the outlet may cause
convergence difficulties.
70 |
BOUNDARY CONDITION
The available Boundary condition options for an inlet are Velocity, Laminar inflow, Mass
flow, and Pressure. After selecting a Boundary Condition from the list, a section with the
The Normal inflow velocity is specified as u = nU0, where n is the boundary normal
pointing out of the domain and U0 is the normal inflow speed.
The Velocity field option sets the velocity vector to u = u0. The components of the inlet
velocity vector u0 should be defined for this choice.
PRESSURE CONDITIONS
This option specifies the normal stress which in most cases is approximately equal to
the pressure. If the reference pressure pref, defined at the physics interface level, is
equal to 0, the value of the Pressure p0, at the boundary, is the absolute pressure.
Otherwise, p0 is the relative pressure at the boundary.
The Suppress backflow option adjusts the inlet pressure locally in order to prevent
fluid from exiting the domain through the boundary. If suppress backflow is
deselected, the inlet boundary can become an outlet depending on the pressure field
in the rest of the domain.
Flow direction controls in which direction the fluid enters the domain.
- For Normal flow it prescribes zero tangential velocity component.
- For User defined, an Inflow velocity direction du (dimensionless) should be
specified. The magnitude of du does not matter, only the direction. du must
point into the domain.
LAMINAR INFLOW
This boundary condition is applicable when the fluid enters the domain from a long
pipe or channel, in which the laminar flow profile is fully developed. The normal stress
at the inlet is determined from the flow conditions at the entrance to a fictitious
channel of length Lentr appended to the boundary. The inflow can be specified by the
Average velocity Uav, the Flow rate V0, or the Entrance pressure pentr.
The Entrance length Lentr should be significantly greater than 0.06ReD, where Re is
the Reynolds number and D is the inlet length scale (hydraulic diameter), in order that
the flow can adjust to a fully developed laminar profile.
71
The Constrain outer edges to zero option forces the laminar profile to go to zero at the
bounding points or edges of the inlet channel. Otherwise the velocity is defined by the
boundary condition of the adjacent boundary in the computational domain. For
example, if one end of a boundary with a Laminar inflow condition connects to a slip
boundary, the laminar profile will have a maximum at that end.
MASS FLOW
The mass flow at an inlet can be specified by the Mass flow rate, the Pointwise mass flux,
the Standard flow rate, or the Standard flow rate (SCCM).
72 |
Standard density, for which the Standard molar volume Vm should be specified.
Standard pressure and temperature, for which the Standard pressure Pst and the
Standard temperature Tst should be defined.
For 2D components, the Channel thickness dbc is used to define the area across which
the mass flow occurs. This setting is not applied to the whole model. Line or surface
integrals of the mass flow over the boundary evaluated during post-processing or used
in integration coupling operators do not include this scaling automatically. Such results
should be appropriately scaled when comparing them with the specified mass flow.
(standard cubic centimeters per minute) without the requirement to specify units.
Here, the dimensionless Number of SCCM units Qsccm should be specified.
CONSTRAINT SETTINGS
Outlet
This condition should be used on boundaries for which there is a net outflow from the
domain. To obtain a numerically well posed problem, it is advisable to also consider
the Inlet condition(s) when specifying an Outlet condition. For example, if the velocity
is specified at the inlet, the pressure may be specified at the outlet, and vice versa.
Specifying the velocity vector at both the inlet and the outlet may cause convergence
difficulties. Selecting appropriate outlet conditions for the Navier-Stokes equations is
a non-trivial task. Generally, if there is something interesting happening at an outflow
boundary, the computational domain should be extended to include this
phenomenon.
BOUNDARY CONDITION
The available Boundary condition options for an outlet are Pressure, Laminar outflow,
and Velocity.
73
PRESSURE CONDITIONS
This option specifies the normal stress which in most cases is approximately equal to
the pressure. The tangential stress component is set to zero. If the reference pressure
pref, defined at the physics interface level, is equal to 0, the value of the Pressure p0, at
the boundary, is the absolute pressure. Otherwise, p0 is the relative pressure at the
boundary.
The Normal flow option changes the no tangential stress condition to a no tangential
velocity condition. This forces the flow to exit (or enter) the domain perpendicularly
to the outlet boundary.
The Suppress backflow check box is selected by default. This option adjusts the outlet
pressure in order to prevent fluid from entering the domain through the boundary.
VE L O C I T Y
This boundary condition is applicable when the flow exits the domain into a long pipe
or channel, at the end of which a laminar flow profile is fully developed. The normal
stress at the outlet is determined from the flow conditions at the end of a fictitious
channel appended to the outlet boundary. The outflow can be specified by the Average
velocity Uav, the Flow rate V0, or the Exit pressure pexit.
The Exit length Lexit should be significantly greater than 0.06ReD, where Re is the
Reynolds number and D is the outlet length scale (hydraulic diameter), in order that
the flow can adjust to a fully developed laminar profile.
The Constrain outer edges to zero option forces the laminar profile to go to zero at the
bounding points or edges of the outlet channel. Otherwise the velocity is defined by
the boundary condition of the adjacent boundary in the computational domain. For
example, if one end of a boundary with a Laminar outflow condition connects to a slip
boundary, the laminar profile will have a maximum at that end.
CONSTRAINT SETTINGS
74 |
Symmetry
The Symmetry boundary condition prescribes no penetration and vanishing shear
stresses. The boundary condition is a combination of a Dirichlet condition and a
Neumann condition:
pI + ( u + ( u ) T ) 2
--- ( u )I n = 0
u n = 0,
u n = 0,
( pI + ( u + ( u ) T ) )n = 0
for the compressible and incompressible formulations. The Dirichlet condition takes
precedence over the Neumann condition, and the above equations are equivalent to
the following equation for both the compressible and incompressible formulations:
u n = 0,
K ( K n )n = 0
K = ( u + ( u ) T )n
BOUNDARY SELECTION
For 2D axial symmetry, a boundary condition does not need to be defined for the
symmetry axis at r = 0. The software automatically provides a condition that prescribes
ur = 0 and vanishing stresses in the z direction and adds an Axial Symmetry node that
implements these conditions on the axial symmetry boundaries only.
CONSTRAINT SETTINGS
Physics Options.
Open Boundary
The Open Boundary conditions describes boundaries in contact with a large volume of
fluid. Fluid can both enter and leave the domain on boundaries with this type of
condition.
BOUNDARY CONDITIONS
The Boundary condition options for open boundaries are Normal stress and No viscous
stress.
Normal Stress
The Normal stress f0 condition implicitly imposes p f0 .
75
No Viscous Stress
The No Viscous Stress condition specifies vanishing viscous stress on the boundary. This
condition does not provide sufficient information to fully specify the flow at the open
boundary and must at least be combined with pressure constraints at adjacent points.
The No viscous stress condition prescribes:
( u + ( u ) T ) 2
--- ( u )I n = 0
3
( u + ( u ) T )n = 0
for the compressible and the incompressible formulations. This condition can be useful
in some situations because it does not impose any constraint on the pressure. A typical
example is a model with volume forces that give rise to pressure gradients that are hard
to prescribe in advance. To make the model numerically stable, this boundary
condition should be combined with a point constraint on the pressure.
Boundary Stress
The Boundary Stress node adds a boundary condition that represents a general class of
conditions also known as traction boundary conditions.
BOUNDARY CONDITION
The Boundary condition options for the boundary stress are General stress, Normal
stress, and Normal stress, normal flow.
General Stress
When General stress is selected, the components for the Stress F should be specified.
The total stress on the boundary is set equal to the given stress F:
pI + ( u + ( u ) T ) 2
--- ( u )I n = F
3
( pI + ( u + ( u ) T ) )n = F
for the compressible and the incompressible formulations.
This boundary condition implicitly sets a constraint on the pressure that for 2D flows is
u n
p = 2 ---------- n F
n
76 |
(3-2)
Normal Stress
Normal Stress is described for the Open Boundary node.
3
T
n ( pI + ( u + ( u ) T ) )n = f 0 ,
tu = 0
tu = 0
(3-3)
Physics Options.
If Normal Stress, Normal Flow is selected as the Boundary condition, then to Apply
reaction terms on all dependent variables, the All physics (symmetric) option should be
77
This section is available when Incompressible flow is selected for Compressibility under
the Physical Model section for the physics interface.
A value or expression should be specified for the Pressure difference, psrc pdst. This
pressure difference can, for example, drive the fully developed flow in a channel.
To set up a periodic boundary condition both boundaries must be selected in the
Periodic Flow Condition node. COMSOL Multiphysics automatically assigns one
boundary as the source and the other as the destination. To manually set the
destination selection, a Destination Selection subnode is available from the context
menu (by right-clicking the parent node) or from the Physics toolbar, Attributes menu.
All destination sides must be connected.
CONSTRAINT SETTINGS
Physics Options.
Pipe Connection
This feature is available with a license for the Pipe Flow Module. For details see Pipe
Connection the in the Pipe Flow Module Users Guide.
Flow Continuity
The Flow Continuity condition is suitable for pairs where the boundaries match; it
prescribes that the flow field is continuous across the pair.
A Wall subnode is added by default and it applies to the parts of the pair boundaries
where a source boundary lacks a corresponding destination boundary and vice versa.
The Wall feature can be overridden by any other boundary condition that applies to
exterior boundaries. By right-clicking the Flow Continuity node, additional Fallback
feature subnodes can be added.
78 |
The relative pressure value is set by specifying the Pressure p0. Or, if the reference
pressure pref defined at the physics interface level is equal to zero, p0 represents the
absolute pressure.
CONSTRAINT SETTINGS
Physics Options.
The source Mass flux, q p should be specified. A positive value results in mass being
ejected from the point into the computational domain. A negative value results in mass
being removed from the computational domain.
79
The Line Mass Source feature is available for all dimensions, but the applicable selection
differs between the dimensions.
MODEL DIMENSION
2D
Points
2D Axisymmetry
3D
Edges
SOURCE STRENGTH
The source Mass flux, q l should be specified. A positive value results in mass being
ejected from the line into the computational domain and a negative value means that
mass is removed from the computational domain.
Line sources located on a boundary affect the adjacent computational domains. This,
for example, has the effect that a line source located on a symmetry plane has twice the
given strength.
80 |
81
82 |
------ + ( u ) = 0
t
(3-4)
u
------- + ( u )u = [ pI + ] + F
t
(3-5)
T
T p
C p ------- + ( u )T = ( q ) + :S ---- ------- ------ + ( u )p + Q
t
T p t
(3-6)
where
is the density (SI unit: kg/m3)
u is the velocity vector (SI unit: m/s)
p is pressure (SI unit: Pa)
is the viscous stress tensor (SI unit: Pa)
F is the volume force vector (SI unit: N/m3)
Cp is the specific heat capacity at constant pressure (SI unit: J/(kgK))
T is the absolute temperature (SI unit: K)
q is the heat flux vector (SI unit: W/m2)
Q contains the heat sources (SI unit: W/m3)
S is the strain-rate tensor:
1
S = --- ( u + ( u ) T )
2
The operation : denotes a contraction between tensors defined by
a:b =
anm bnm
(3-7)
n m
83
(3-8)
The dynamic viscosity, (SI unit: Pas), for a Newtonian fluid is allowed to depend on
the thermodynamic state but not on the velocity field. All gases and many liquids can
be considered Newtonian. Examples of non-Newtonian fluids are honey, mud, blood,
liquid metals, and most polymer solutions. With the Microfluidics Module, you can
model flows of non-Newtonian fluids using the predefined power law and Carreau
models, which describe the dynamic viscosity for non-Newtonian [Link]
commonly used constitutive relations are Fouriers law of heat conduction and the
ideal gas law.
In theory, the same equations describe both laminar and turbulent flows. In practice,
however, the mesh resolution required to simulate turbulence with the Laminar Flow
interface makes such an approach impractical.
There are several books where derivations of the Navier-Stokes equations
and detailed explanations of concepts such as Newtonian fluids can be
found. See, for example, the classical text by Batchelor (Ref. 3) and the
more recent work by Panton (Ref. 4).
Many applications describe isothermal flows for which Equation 3-6 is decoupled from
Equation 3-4 and Equation 3-5.
2D AXISYMMETRIC FORMULATIONS
84 |
= 0
u = 0
Compressible Flow
The equations of motion for a single-phase fluid are the continuity equation:
------ + ( u ) = 0
t
(3-9)
u
2
------- + u u = p + ( u + ( u ) T ) --- ( u )I + F
t
3
(3-10)
These equations are applicable for incompressible as well as compressible flow with
density variations.
85
boundary into the domain changes as the Mach number passes through unity. Hence,
the number of boundary conditions required to obtain a numerically well posed system
must also change. The compressible formulation of the laminar and turbulent
interfaces uses the same boundary conditions as the incompressible formulation, which
implies that the compressible interfaces are not suitable for flows with Mach number
larger than or equal to one.
The practical Mach number limit is lower than one, however. The main reason is that
the numerical scheme (stabilization and boundary conditions) of the Laminar Flow
interface does not recognize the direction and speed of pressure waves. The fully
compressible Navier-Stokes equations do, for example, start to display very sharp
gradients already at moderate Mach numbers. But the stabilization for the single-phase
flow interface does not necessarily capture these gradients. It is impossible to give an
exact limit where the low Mach number regime ends and the moderate Mach number
regime begins, but a rule of thumb is that the Mach number effects start to appear at
Ma = 0.3. For this reason the compressible formulation is referred to as Compressible
flow (Ma<0.3) in COMSOL Multiphysics.
Incompressible Flow
When the temperature variations in a flow are small, a single-phase fluid can often be
assumed incompressible; that is, is constant or nearly constant. This is the case for all
liquids under normal conditions and also for gases at low velocities. For constant ,
Equation 3-9 reduces to
u = 0
(3-11)
u
T
+ ( u )u = [ pI + ( u + ( u ) ) ] + F
t
(3-12)
86 |
numbers, viscous forces dominate and tend to damp out all disturbances, which leads
to laminar flow. At high Reynolds numbers, the damping in the system is very low
giving small disturbances the possibility to grow by nonlinear interactions. If the
Reynolds number is high enough, the flow field eventually ends up in a chaotic state
called turbulence.
Observe that the Reynolds number can have different meanings depending on the
length scale and velocity scale. To be able to compare two Reynolds numbers, they
must be based on equivalent length and velocity scales.
The Fluid Flow interfaces automatically calculate the local cell Reynolds number
Rec = |u|h/(2) using the element length h for L and the magnitude of the velocity
vector u for the velocity scale U. This Reynolds number is not related to the character
of the flow field, but to the stability of the numerical discretization. The risk for
numerical oscillations in the solution increases as Rec grows. The cell Reynolds
number is a predefined quantity available for visualization and evaluation (typically it
is available as: [Link]).
= =
1
--- :
2
anm bnm
n m
87
= ()
The Laminar Flow interfaces have the following predefined models to prescribe a
non-Newtonian viscositythe power law and the Carreau model.
POWER LAW
= m n 1
(3-13)
where m and n are scalars that can be set to arbitrary values. For n > 1, the power law
describes a shear thickening (dilatant) fluid. For n < 1, it describes a shear thinning
(pseudoplastic) fluid. A value of n equal to one gives the expression for a Newtonian
fluid.
Equation 3-13 predicts an infinite viscosity at zero shear rate for n < 1. This is however
never the case physically. Instead, most fluids have a constant viscosity for shear rates
smaller than 10-2 s-1 (Ref. 16). Since infinite viscosity also makes models using
Equation 3-13 difficult to solve, COMSOL Multiphysics implements the power law
model as
= mmax ( , min ) n 1
(3-14)
where min is a lower limit for the evaluation of the shear rate magnitude. The default
value for min is 10-2 s-1, but can be given an arbitrary value or expression using the
corresponding text field.
CARREAU MODEL
The Carreau model defines the viscosity in terms of the following four-parameter
expression
(n 1)
---------------- = + ( 0 inf ) [ 1 + ( ) 2 ] 2
(3-15)
where is a parameter with the unit of time, 0 is the zero shear rate viscosity, inf is
the infinite shear-rate viscosity, and n is a dimensionless parameter. This expression is
able to describe the viscosity for most stationary polymer flows.
88 |
(3-16)
where g is the gravity vector. A further simplifications is often possible. Since g can be
written in terms of a potential, , it is possible to write Equation 3-16 as:
F = ( 0 ) + g
The first part can be canceled out by splitting the true pressure, p, into a hydrodynamic
component, P, and a hydrostatic component, 0. Equation 3-11 and
Equation 3-12 are expressed in terms of the hydrodynamic pressure P = p + 0:
u = 0
(3-17)
u
0 ------- + ( 0 u )u = P + ( ( u + ( u ) T ) ) + g
t
(3-18)
To obtain the Boussinesq approximation in this form, enter the expression for g for
the Wall feature.
In practice, the shift from p to P can be ignored except where the pressure appears in
boundary conditions. The pressure that is specified at boundaries is the hydrodynamic
pressure in this case. For example, on a vertical outflow or inflow boundary, the
hydrodynamic pressure is typically a constant, while the true pressure is a function of
the vertical coordinate.
The system that Equation 3-17 and Equation 3-18 form has its limitations. The main
assumption is that the density fluctuations must be small; that is, /0 << 1. There are
also some more subtle constraints that, for example, make the Boussinesq
approximation unsuitable for systems of very large dimensions. An excellent discussion
of the Boussinesq approximation and its limitations appears in Chapter 14 of Ref. 10.
89
SLIP
The Slip condition assumes that there are no viscous effects at the slip wall and hence,
no boundary layer develops. From a modeling point of view, this is a reasonable
approximation if the important effect of the wall is to prevent fluid from leaving the
domain. Mathematically, the constraint can be formulated as:
u n = 0,
( pI + ( u + ( u ) T ) )n = 0
The no penetration term takes precedence over the Neumann part of the condition
and the above expression is therefore equivalent to
u n = 0,
K ( K n )n = 0
K = ( u + ( u ) T )n
expressing that there is no flow across the boundary and no viscous stress in the
tangential direction.
The boundary condition for is n = 0 .
SL ID IN G WAL L
The Sliding Wall boundary condition is appropriate if the wall behaves like a conveyor
belt; that is, the surface is sliding in its tangential direction. The wall does not have to
actually move in the coordinate system.
In 2D, the tangential direction is unambiguously defined by the direction of the
boundary, but the situation becomes more complicated in 3D. For this reason, this
boundary condition has slightly different definitions in the different space
dimensions.
For 2D and 2D axisymmetric components, the velocity is given as a scalar Uw and
the condition prescribes
u n = 0,
u t = Uw
where t = (ny, nx) for 2D and t = (nz, nr) for axial symmetry.
For 3D components, the velocity is set equal to a given vector uw projected onto
the boundary plane:
u w ( n u w )n
u = -------------------------------------------- u w
u w ( n u w )n
The normalization makes u have the same magnitude as uw even if uw is not exactly
parallel to the wall.
90 |
S L I P VE L O C I T Y
In the microscale range, the flow at a boundary is seldom strictly no slip or slip.
Instead, the boundary condition is something in between, and there is a slip velocity
at the boundary. Two phenomena account for this velocity: violation of the continuum
hypothesis for the viscosity and flow induced by a thermal gradient along the
boundary.
The following equation relates the viscosity-induced jump in tangential velocity to the
tangential shear stress along the boundary:
1
u = --- n, t
= ------------------------2
v
--------------
v
where is the fluids dynamic viscosity (SI unit: Pas), v represents the tangential
momentum accommodation coefficient (TMAC) (dimensionless), and is the
molecules mean free path (SI unit: m). The tangential accommodation coefficients are
typically in the range of 0.85 to 1.0 and can be found in Ref. 15.
A simpler expression for is
= -----Ls
where Ls, the slip length (SI unit: m), is a straight channel measure of the distance
from the boundary to the virtual point outside the flow domain where the flow profile
extrapolates to zero. This equation holds for both liquids and gases.
Thermal creep results from a temperature gradient along the boundary. The following
equation relates the thermally-induced jump in tangential velocity to the tangential
gradient of the natural logarithm of the temperature along the boundary:
u = T --- t log T
where T is the thermal slip coefficient (dimensionless) and is the density of the fluid.
The thermal slip coefficients range between 0.3 and 1.0 and can be found in Ref. 15.
Combining the previous relationships results in the following equation:
91
Ls
u u w, t = ------ n, t + T ------- t T
T
Relate the tangential shear stress to the viscous boundary force by
n, t = K ( n K )n
where the components of K are the Lagrange multipliers that are used to implement
the boundary condition. Similarly, the tangential temperature gradient results from the
difference of the gradient and its normal projection:
t T = T ( n T )n
Most solid surfaces acquire a surface charge when brought into contact with an
electrolyte. In response to the spontaneously formed surface charge, a charged
solution forms close to the liquid-solid interface. This is known as an electric double
layer. If an electric field is applied to the fluid, this very narrow layer starts to move
along the boundary.
It is possible to model the fluids velocity near the boundary using the
Helmholtz-Smoluchowski relationship between the electroosmotic velocity u and the
applied electric field:
u = eo E t
where eo is the electroosmotic mobility and Et is the fluid electric field tangential to
the wall.
Built-in Expression
Use the Electroosmotic mobility eo (SI unit: m2/(sV)) Built-in expression to
compute the electroosmotic mobility from:
eo = r 0 --
92 |
(3-19)
An inlet requires specification of the velocity components. The most robust way to do
this is to prescribe a velocity field using a Velocity condition.
A common alternative to prescribing the complete velocity field is to prescribe a
pressure and all but one velocity component. The pressure cannot be specified
pointwise since this is a mathematically over-constraining. Instead the pressure can be
specified via a stress condition:
u n
p + 2 --------- = F n
n
(3-20)
93
u t
-------- = 0
n
which is what the Normal stress condition does. Vanishing tangential stress becomes a
less well posed inlet condition as the Reynolds number increases. The Pressure
condition in the Inlet feature therefore requires a flow direction to be prescribed which
provides a well-posed condition independent of Reynolds number.
OUTLET CONDITIONS
The most common approach is to prescribe a pressure via a normal stress condition on
the outlet. This is often accompanied by a vanishing tangential stress condition:
u t
-------- = 0
n
where ut/n is the normal derivative of the tangential velocity field. It is also possible
to prescribe ut to be zero. The latter option should be used with care since it can have
a significant effect on the upstream solution.
The elliptic character of the Navier-Stokes equations mathematically permit specifying
a complete velocity field at an outlet. This can however be difficult to apply in practice.
The reason being that it is hard to prescribe the outlet velocity so that it at each point
is consistent with the interior solution. The adjustment to the specified velocity then
occurs across an outlet boundary layer. The thickness of this boundary layer depends
on the Reynolds number; the higher the Reynolds number, the thinner the boundary
layer.
ALTERNATIVE FORMULATIONS
94 |
Laminar Inflow
In order to prescribe a fully developed inlet velocity profile, this boundary condition
adds a weak form contribution and constraints corresponding to unidirectional flow
perpendicular to the boundary. The applied condition corresponds to the situation
shown in Figure 3-1: a fictitious domain of length Lentr is assumed to be attached to
the inlet of the computational domain. The domain is an extrusion of the inlet
boundary, which means that laminar inflow requires the inlet to be flat. The boundary
condition uses the assumption that the flow in this fictitious domain is fully developed
laminar flow. The wall boundary conditions for the fictitious domain is inherited
from the real domain, , unless the option to constrain outer edges or endpoints to
zero is selected in which case the fictitious walls are no-slip walls.
pentr
Lentr
Figure 3-1: An example of the physical situation simulated when using the Laminar
inflow boundary condition. is the actual computational domain while the dashed
domain is a fictitious domain.
If an average inlet velocity or inlet volume flow is specified instead of the pressure,
COMSOL Multiphysics adds an ODE that calculates a pressure, pentr, such that the
desired inlet velocity or volume flow is obtained.
Laminar Outflow
In order to prescribe an outlet velocity profile, this boundary condition adds a weak
form contribution and constraints corresponding to unidirectional flow perpendicular
to the boundary. The applied condition corresponds to the situation shown in
Figure 3-2: assume that a fictitious domain of length Lexit is attached to the outlet of
the computational domain. The domain is an extrusion of the outlet boundary, which
means that laminar outflow requires the outlet to be flat. The boundary condition uses
the assumption that the flow in this fictitious domain is fully developed laminar flow.
The wall boundary conditions for the fictitious domain is inherited from the real
95
domain, , unless the option to constrain outer edges or endpoints to zero is selected
in which case the fictitious walls are no-slip walls.
pexit
Lexit
Figure 3-2: An example of the physical situation simulated when using the Laminar
outflow boundary condition. is the actual computational domain while the dashed
domain is a fictitious domain.
If the average outlet velocity or outlet volume flow is specified instead of the pressure,
the software adds an ODE that calculates pexit such that the desired outlet velocity or
volume flow is obtained.
Mass Flow
The Mass Flow boundary condition constrains the mass flowing into the domain across
an inlet boundary. The mass flow can be specified in a number of ways.
POINTWISE MASS FLUX
The pointwise mass flux sets the velocity at the boundary to:
mf
u = ------- n
The mass flow rate boundary condition sets the total mass flow through the boundary
according to:
dbc ( u n ) dS = m
where dbc (only present in the 2D Cartesian axis system) is the boundary thickness
normal to the fluid-flow domain and m is the total mass flow rate.
In addition to the constraint on the total flow across the boundary, the tangential
velocity components are set to zero on the boundary
96 |
un = 0
(3-21)
The standard flow rate boundary condition specifies the mass flow as a standard
volumetric flow rate. The mass flow through the boundary is set by the equation:
( u n ) dS = Q sv
dbc ------ st
where dbc (only present in the 2D component Cartesian axis system) is the boundary
thickness normal to the fluid-flow domain, st is the standard density, and Qsv is the
standard flow rate. The standard density is defined by one of the following equations:
Mn
st = -------Vn
p st M n
st = ----------------RT st
where Mn is the mean molar mass of the fluid, Vn is the standard molar volume, pst is
the standard pressure, R is the universal molar gas constant, and Tst is the standard
temperature.
Equation 3-21 or Equation 3-22 is also enforced for compressible and incompressible
flow, respectively, ensuring that the normal component of the viscous stress and the
tangential component of the velocity are zero at the boundary.
No Viscous Stress
For this module, and in addition to the Pressure, No Viscous Stress boundary
condition, the viscous stress condition sets the viscous stress to zero:
( u + ( u ) T ) 2
--- ( u )I n = 0
3
( ( u + ( u ) T ) )n = 0
using the compressible and the incompressible formulation, respectively.
97
The condition is not a sufficient outlet condition since it lacks information about the
outlet pressure. It must hence be combined with pressure point constraints on one or
several points or lines surrounding the outlet.
This boundary condition is numerically the least stable outlet condition, but can still
be beneficial if the outlet pressure is nonconstant due to, for example, a nonlinear
volume force.
3
( pI + ( u + ( u ) T ) )n = f 0 n
using the compressible and the incompressible formulation, respectively.
This implies that the total stress in the tangential direction is zero. This boundary
condition implicitly sets a constraint on the pressure which for 2D flows is
u n
p = 2 ---------- + f 0
n
(3-22)
(3-23)
( pI + ( u + ( u ) T ) )n = p 0 n
(3-24)
98 |
(3-25)
pI + ( u + ( u ) T ) 2
--- ( u )I n = p 0 n
3
( pI + ( u + ( u ) T ) )n = p 0 n
p p
0
(3-26)
(3-27)
3
T
n ( pI + ( u + ( u ) T ) )n = p 0
p p
0
(3-28)
together with the tangential condition in Equation 3-25, or, a general flow direction
is prescribed.
99
T pI + ( u + ( u ) T ) 2
--- ( u )I n = p 0 ( r n )
ru
u
3
(r n)
T ( pI + ( u + ( u ) T ) )n = p
ru
0 u
p p
0
(3-29)
du
u ( u r u )r u = 0, r u = ------------du
The > option is used with suppress backflow to have u n 0 or u r u 0 .
See Inlet, Outlet, Open Boundary, and No Viscous Stress for the individual node
settings. Note that some modules have additional theory sections describing options
available with that module.
V 0
Q = qp
(3-30)
100 |
For this alternative approach, effects resulting from the physical object volume, such
as drag and fluid displacement, need to be neglected.
The weak contribution
q p test ( p )
is added to a point in the geometry. As can be seen from Equation 3-30, Q must tend
to plus or minus infinity as V tends to zero. This means that in theory the pressure
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.
The finite element representation of Equation 3-30 corresponds to a finite pressure in
a point with the effect of the point source spread out over a region around the point.
The size of the region depends on the mesh and on the strength of the source. A finer
mesh gives a smaller affected region, but also a more extreme pressure value. It is
important not to mesh too finely around a point source since the resulting pressure can
result in unphysical values for the density, for example. It can also have a negative effect
on the condition number for the equation system.
LINE SOURCE
A line source can theoretically be formed by assuming a source of strength Q (SI unit:
kg/(m3s)), located within a tube with cross-section area 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 (SI unit: kg/(ms)), this can be expressed as
lim
S 0
= q l
(3-31)
101
For feature node information, see Line Mass Source and Point Mass
Source in the COMSOL Multiphysics Reference Manual.
For the Reacting Flow in Porous Media, Diluted Species interface, which
is available with the CFD Module, Chemical Reaction Engineering
Module, or Batteries & Fuel Cells Module, these shared physics nodes are
renamed as follows:
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).
102 |
STREAMLINE DIFFUSION
For strongly coupled systems of equations, the streamline diffusion method must be
applied to the system as a whole rather than to each equation separately. These ideas
were first explored by Hughes and Mallet (Ref. 7) and were later extended to Galerkin
least-squares (GLS) applied to the Navier-Stokes equations (Ref. 8). This is the
streamline diffusion formulation that COMSOL Multiphysics supports. The time-scale
tensor is the diagonal tensor presented in Ref. 9.
Streamline diffusion is active by default because it is necessary when convection is
dominating the flow.
The governing equations for incompressible flow are subject to the Babuska-Brezzi
condition, which states that the shape functions (basis functions) for pressure must be
of lower order than the shape functions for velocity. If the incompressible
Navier-Stokes equations are stabilized by streamline diffusion, it is possible to use
equal-order interpolation. Hence, streamline diffusion is necessary when using
first-order elements for both velocity and pressure. This applies also if the model is
solved using geometric multigrid (either as a solver or as a preconditioner) and at least
one multigrid hierarchy level uses linear Lagrange elements.
CROSSWIND DIFFUSION
Crosswind diffusion can also be formulated for systems of equations, and when applied
to the Navier-Stokes equations it becomes a shock-capturing operator. COMSOL
Multiphysics supports the formulation in Ref. 8 with a shock capturing viscosity of the
Hughes-Mallet type Ref. 7.
Incompressible flows do not contain shock waves, but crosswind diffusion is still useful
for introducing extra diffusion in sharp boundary layers and shear layers that otherwise
would require a very fine mesh to resolve.
Crosswind diffusion is active by default as it makes it easier to obtain a solution even if
the problem is fully resolved by the mesh. Crosswind diffusion also enables the iterative
solvers to use inexpensive presmoothers. If crosswind diffusion is deactivated, more
expensive preconditioners must be used instead.
103
ISOTROPIC DIFFUSION
Stationary Solver
In the stationary case, a fully coupled, damped Newton method is applied. The initial
damping factor is low since a full Newton step can be harmful unless the initial values
are close to the final solution. The nonlinear solver algorithm automatically regulates
the damping factor in order to reach a converged solution.
For advanced models, the automatically damped Newton method might not be robust
enough. A pseudo time-stepping algorithm can then be invoked. See Pseudo Time
Stepping for Laminar Flow Models.
Time-Dependent Solver
In the time-dependent case, the initial guess for each time step is (loosely speaking) the
previous time step, which is a very good initial value for the nonlinear solver. The
automatic damping algorithm is then not necessary. The damping factor in the
Newton method is instead set to a constant value slightly smaller than one. Also, for
the same reason, it suffices to update the Jacobian once per time-step.
104 |
It is seldom worth the extra computational cost to update the Jacobian more than once
per time step. For most models it is more efficient to restrict the maximum time step
or possibly lower the damping factor in the Newton method.
LINEAR SOLVER
The linearized Navier-Stokes equation system has saddle point character, unless the
density depends on the pressure. This means that the Jacobian matrix has zeros on the
diagonal. Even when the density depends on the pressure the equation system
effectively shares many numerical properties with a saddle point system.
For small 2D and 3D models, the default solver suggestion is a direct solver. Direct
solvers can handle most nonsingular systems and are very robust and also very fast for
small models. Unfortunately, they become slow for large models and their memory
requirement scales as somewhere between N1.5and N2, where N is the number of
degrees of freedom in the model. The default suggestion for large 2D and 3D models
is therefore the iterative GMRES solver. The memory requirement for an iterative
solver optimally scales as N.
Geometric Multigrid (GMG) is used to accelerate GMRES. GMG needs smoothers
but the saddle point character of the linear system restricts the number of applicable
smoothers. The choices are further restricted by the anisotropic meshes frequently
encountered in fluid-flow problems. Pointwise smoothers, such as SOR, are not very
efficient on anisotropic meshes.
The efficiency of the smoothers is highly dependent on the numerical stabilization.
Iterative solvers perform at their best when both Streamline Diffusion and Crosswind
Diffusion are active.
The default smoother for P1+P1 elements is SCGS. This is an efficient and robust
smoother specially designed to solve saddle point systems on meshes that contain
anisotropic elements. The SCGS smoother works well even without crosswind
diffusion. SCGS can sometimes work for higher-order elements, especially if Method in
the SCGS settings is set to Mesh element lines. But there is no guarantee for this, so the
default smoother for P2+P1 elements and P3+P2 elements is an SOR Line smoother.
SOR Line handles mesh anisotropy but does not formally address the saddle point
character. It does, however, function in practice provided that streamline diffusion and
crosswind diffusion are both active.
A different kind of saddle point character can arise if the equation system contains
ODE variables. Some advanced boundary conditions can add equations with such
variables. These variables must be treated with the Vanka algorithm. SCGS includes an
105
option to invoke Vanka. Models with higher-order elements must apply SCGS or use
the Vanka smoother. The latter is the default suggestion for higher-order elements, but
it does not work optimally for anisotropic meshes.
TIME-DEPENDENT SOLVERS
The default time-dependent solver for Navier-Stokes is the BDF method with
maximum order set to two. Higher BDF orders are not stable for transport problems
in general nor for Navier-Stokes in particular.
BDF methods have been used for a long time and are known for their stability.
However, they can have severe damping effects, especially the lower-order methods.
Hence, if robustness is not an issue, a model can benefit from using the generalized-
method instead. Generalized- is a solver which has properties similar to those of the
second-order BDF solver but it is much less diffusive.
Both BDF and generalized- are per default set to automatically adjust the time step.
While this works well for many models, extra efficiency and accuracy can often be
gained by specifying a maximum time step. It is also often beneficial to specify an initial
time step to make the solver progress smoothly in the beginning of the time series.
In the COMSOL Multiphysics Reference Manual:
Time-Dependent Solver
Multigrid, Direct, Iterative, SCGS, SOR Line, and Vanka
Stationary Solver
( u )u = [ pI + ( u + ( u ) ) ] + F
(3-32)
Solving Equation 3-32 requires a starting guess that is close enough to the final
solution. If no such guess is at hand, the fully transient problem can be solved instead.
This is, however, a rather costly approach in terms of computational time. An
intermediate approach is to add a fictitious time derivative to Equation 3-32:
T
u nojac ( u )
--------------------------------- + ( u )u = [ pI + ( u + ( u ) ) ] + F
t
106 |
where t is a pseudo time step. Since unojac(u) 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 step of size t of a time-dependent
solver.
Pseudo time stepping is not active per default. The pseudo time step t can be chosen
individually for each element based on the local CFL number:
h
t = CFL loc ------u
where h is the mesh cell size. A small CFL number means a small time step. It is
practical to start with a small CFL number and gradually increase it as the solution
approaches steady state.
If the automatic expression for CFLloc is set to the built-in variable CFLCMP. The
automatic setting then suggests a PID regulator for the pseudo time step in the default
solver. The PID regulator starts with a small CFL number and increases CFLloc as the
solution comes closer to convergence.
The default manual expression is
1.3 min ( niterCMP, 9 ) +
if ( niterCMP > 20, 9 1.3 min ( niterCMP 20, 9 ), 0 ) +
if ( niterCMP > 40, 90
(3-33)
0)
The variable niterCMP is the nonlinear iteration number. It is equal to one for the first
nonlinear iteration. CFLloc starts at 1.3 and increases by 30% each iteration until it
reaches 1.3 9 10.6 . It remains there until iteration number 20 at which it starts to
increase until it reaches approximately 106. A final increase after iteration number 40
then takes it to 1060. Equation 3-33 can for some advanced flows increase CFLloc too
slowly or too quickly. CFLloc can then be tuned for the specific application.
For details about the CFL regulator, see Pseudo Time Stepping in the
COMSOL Multiphysics Reference Manual.
107
d x
dx
= F t, x,
dt
dt2
where x is the position of the particle, m the particle mass, and F is the sum of all forces
acting on the particle. Examples of forces acting on a particle in a fluid are the drag
force, the buoyancy force, and the gravity force. The drag force represents the force
that a fluid exerts on a particle due to a difference in velocity between the fluid and the
particle. It includes the viscous drag, the added mass, and the Basset history term.
Several empirical expressions have been suggested for the drag force. One of those is
the one proposed by Khan and Richardson (Ref. 11). That expression is valid for
108 |
spherical particles for a wide range of particle Reynolds numbers. The particle
Reynolds number is defined as
u u p 2r
Re p = -----------------------------
where u is the velocity of the fluid, up the particle velocity, r the particle radius, the
fluid density, and the dynamic viscosity of the fluid. The empirical expression for the
drag force according to Khan and Richardson is
2
-0.31
F = r u u p ( u u p ) [ 1.84Re p
+ 0.293Re p0.06 ]
3.45
109
110 |
111
112 |
C H A P T E R 4 : M U L T I P H A S E F L O W , TW O - P H A S E F L O W I N T E R F A C E S
Initial Interface, and Initial Values. Then, from the Physics toolbar, add other nodes that
implement, for example, boundary conditions and volume forces. You can also
right-click Laminar Two-Phase Flow, Level Set to select physics features from the context
menu.
Except where included below, the Laminar Two-Phase Flow, Level Set has the same
sections and settings as the Laminar Flow interface.
SETTINGS
The Laminar Two-Phase Flow, Level-Set interface uses a level set method to track the
fluid-fluid interface.
This physics interface changes to a Laminar Two-Phase Flow, Phase Field interface if the
Multiphase flow model selected is Two-phase flow, phase field. See The Laminar
Two-Phase Flow, Phase Field Interface for details.
The Compressibility defaults to Incompressible flow (constant density flow). Select
Compressible flow (Ma<0.3) to use the compressible formulation of the Navier-Stokes
equations.
Select the Neglect inertial term (Stokes flow) check box to model flow at very low
Reynolds numbers for which the inertial term in the Navier-Stokes equations can be
neglected. The physics interface then solves the linear Stokes equations instead. The
Stokes equations are valid for creeping flow, which can occur in microfluidics and
MEMS devices, where the flow speed or length scales are very small.
Enter a Reference pressure level pref (SI unit: Pa). The default value is 1[atm].
DEPENDENT VA RIA BLES
T H E L A M I N A R T W O - P H A S E F L O W , L E V E L S E T A N D L A M I N A R TW O - P H A S E F L O W , P H A S E F I E L D I N T E R F A C E S
113
Domain, Boundary, Point, and Pair Nodes for the Laminar and
Turbulent Flow, Two-Phase, Level Set and Phase Field Interfaces
Theory for the Level Set and Phase Field Interfaces
114 |
C H A P T E R 4 : M U L T I P H A S E F L O W , TW O - P H A S E F L O W I N T E R F A C E S
The main node is the Fluid Properties feature, which adds the Navier-Stokes equations
and the phase field equations and provides an interface for defining the properties of
the fluids and the surface tension.
When this physics interface is added, the following default nodes are also added in the
Model BuilderFluid Properties, Wall (with a default No slip boundary condition),
Initial Interface, and Initial Values. Then, from the Physics toolbar, add other nodes that
implement, for example, boundary conditions and volume forces. You can also
right-click Laminar Two-Phase Flow, Phase Field to select physics features from the
context menu.
Except for the Physical Model section, the settings and model examples are the same as
for The Laminar Two-Phase Flow, Level Set Interface.
PHYSICAL MODEL
The Laminar Two-Phase Flow, Phase Field interface uses a phase field method to track
the fluid-fluid interface.
Domain, Boundary, Point, and Pair Nodes for the Laminar and
Turbulent Flow, Two-Phase, Level Set and Phase Field Interfaces
Theory for the Level Set and Phase Field Interfaces
Domain, Boundary, Point, and Pair Nodes for the Laminar and
Turbulent Flow, Two-Phase, Level Set and Phase Field Interfaces
The Laminar Two-Phase Flow, Level Set and Laminar Two-Phase Flow, Phase Field
Interfaces has these domain, boundary, point, and pair nodes, which are available from
the Physics ribbon toolbar (Windows users), Physics context menu (Mac or Linux
users), or right-click to access the context menu (all users).
In general, to add a node, go to the Physics toolbar, no matter what
operating system you are using. Subnodes are available by clicking the
parent node and selecting it from the Attributes menu.
In addition to the surface level boundary conditions listed below, the pressure point
constraint boundary condition is also available and should be used in addition to the
boundary conditions listed when the pressure level is not specified by any other
boundary condition.
T H E L A M I N A R T W O - P H A S E F L O W , L E V E L S E T A N D L A M I N A R TW O - P H A S E F L O W , P H A S E F I E L D I N T E R F A C E S
115
Initial Values
Gravity
Wall
Initial Interface
The following nodes (listed in alphabetical order) are described for the Laminar Flow
interface:
Flow Continuity
Inlet
Symmetry
Outlet
Volume Force
Wall
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.
Wall
The Wall node represents the walls in a two-phase flow simulation. The Wall feature
for this physics interface includes slip, no slip, and slip velocity boundaries as well as
sliding, leaking, moving, and wetted walls. The No slip boundary condition is the
default boundary condition. These boundary conditions are described for the Laminar
Flow interface.
The Wetted wall and Moving, wetted wall boundary conditions are described in this
section. See the Laminar Flow interface for the other settings (Wall).
Wetted Wall
The Wetted wall boundary condition is suitable for walls in contact with the fluid-fluid
interface. If this boundary condition is used, the fluid-fluid interface can move along
the wall. For applications where the interface is fixed on the wall, the no slip condition
is suitable.
116 |
C H A P T E R 4 : M U L T I P H A S E F L O W , TW O - P H A S E F L O W I N T E R F A C E S
The implementation of the wetted wall boundary condition differs depending on the
method used to track the fluid-fluid interface (the level set method or the phase field
method).
Level Set For the level set method, this boundary condition enforces the
F fr = --- u
where is the slip length. For numerical calculations it is suitable to set = h, where
h is the mesh element size. The boundary condition does not set the tangential velocity
component to zero; however, the extrapolated tangential velocity component is 0 at a
distance outside the wall (see Figure 4-1).
Finally, the boundary condition adds the following weak boundary term:
The boundary term results from the partial integration of the surface tension force in
the momentum equation. Define the contact angle w (that is, the angle between the
fluid interface and the wall). Figure 4-1 illustrates the definition of the contact angle.
Fluid 1
Wall
Wall
Fluid 2
Figure 4-1: Definition of the contact angle at interface/wall contact points (left) and
an illustration of the slip length (right).
With the level set method, define the following two properties for the wetted wall:
Enter a value or expression for the Contact angle w. The default is pi/2 (/2) rad.
Enter a value or expression for the Slip length (SI unit: m). The default is h, which
is the variable for the local mesh element size h.
T H E L A M I N A R T W O - P H A S E F L O W , L E V E L S E T A N D L A M I N A R TW O - P H A S E F L O W , P H A S E F I E L D I N T E R F A C E S
117
Phase Field The motion of the interface on the boundary due to advection is zero and
so the no slip boundary condition, u = 0, is used in the momentum equation. The
following boundary condition defines the contact angle between Fluid 2 and the wall:
2
n = cos ( w )
(4-1)
where w is the user-defined contact angle and n is the unit vector normal to the wall.
The phase field help variable is assigned the boundary condition:
n -----2- = 0
(4-2)
With the phase field method, enter a value or expression for the Contact angle w. The
default value is pi/2 (/2) rad.
Fluid Properties
The Fluid Properties node adds the flow equations and the equations for the level set
variable or the phase field variable (where applicable). For the level set and phase field
interfaces, specify the fluid that each domain is initially filled with using the Initial
Values node.
FLUID 1 AND FLUID 2 PROPERTIES
Specify the Density and the Dynamic viscosity of Fluid 1 and Fluid 2 in the same way as
in a single-phase flow interface (see Fluid Properties).
Care should be taken when using the Domain Material setting for the
material properties for Fluid 1 and Fluid 2.
The material properties are obtained from the domain irrespective of the location of
the interface. If two different materials are selected in domains 1 and 2, with the phase
boundary initially coincident with the domain boundary, the model has convergence
118 |
C H A P T E R 4 : M U L T I P H A S E F L O W , TW O - P H A S E F L O W I N T E R F A C E S
issues once the phase boundary moves away from the domain boundary. This is
because a density discontinuity and a viscosity discontinuity occurs at the boundary
separating the two fluids. For this reason selecting the material directly is
recommended when setting the material properties for Fluid 1 and Fluid 2.
The fluid defined as Fluid 1 affects the wetting characteristics on wetted walls. See the
Wall node for details.
S U R F A C E TE N S I O N
This section is not available for the turbulent versions of these physics interfaces.
Select the Neglect surface tension in momentum equation check box to neglect surface
tension.
Select a Surface tension coefficient (SI unit: N/m):
To use a predefined expression, select Library coefficient, liquid/gas interface or
Library coefficient, liquid/liquid interface. Then select an option from the list that
displays below (for example, Water/Air, Glycerol/Air and so forth).
For User defined enter a value or expression for the surface tension coefficient (SI
unit: N/m).
LEVEL SET PARAMETERS
This section is available with The Laminar Two-Phase Flow, Level Set Interface and
when turbulence is active.
Specify the parameters for the level set method:
Enter a value or expression for the Reintialization parameter (SI unit: m/s). The
default value is 1 m/s. Use an approximate value for the maximum speed occurring
in the flow. Because the flow speed is not always known in advance, an initial
computation might be needed to find a proper value for .
Enter a value or expression for the Parameter controlling interface thickness ls
(SI unit: m). The default expression is [Link]/2, which corresponds to half of
the maximum mesh element size in the model. In general, the results are optimal if
the default value is used for this parameter.
PHASE FIELD PARAMETERS
This section is available with The Laminar Two-Phase Flow, Phase Field Interface and
when turbulence is active.
Define the parameters for the phase field method.
T H E L A M I N A R T W O - P H A S E F L O W , L E V E L S E T A N D L A M I N A R TW O - P H A S E F L O W , P H A S E F I E L D I N T E R F A C E S
119
Enter a value or expression for the Parameter controlling interface thickness pf (SI unit:
m). The default expression is [Link]/2, which corresponds to half of the maximum
mesh element size in the model. In general, simply use the default value for optimal
accuracy.
Enter a value or expression for the Mobility tuning parameter (SI unit: ms/kg). The
default is 1 ms/kg, which is a good starting point for most models. This parameter
determines the time scale of the Cahn-Hilliard diffusion and thereby also governs the
diffusion-related time scale of the interface. Keep the parameter value high enough
to maintain a constant interface thickness yet low enough not to damp the convective
motion. A too high mobility value can also lead to excessive diffusion of droplets.
EXTERNAL FREE ENERGY
This section is available with The Laminar Two-Phase Flow, Phase Field Interface and
when turbulence is active.
The expression for the external free energy must be manually
differentiated with respect to and then entered into the derivative
of external free energy field for f (SI unit: J/m3).
Add a source of external free energy as an external force term in the flow equation (the
Fext term in Equation 3-2). The external free energy fext is a user-defined free energy.
In most cases, the external free energy is zero.
Gravity
The Gravity node adds the force Fg in Equation 3-2, which is equal to g, where g is
the gravity vector.
GRAVITY
Enter the coordinates for the Gravity vector g (SI unit: m/s2).
For 3D and 2D axisymmetric components, the default value is -g_const in the z
component. g_const is the acceleration of gravity (a predefined physical constant).
For 2D components, the default value is -g_const in the y component. g_const is
the acceleration of gravity (a predefined physical constant).
120 |
C H A P T E R 4 : M U L T I P H A S E F L O W , TW O - P H A S E F L O W I N T E R F A C E S
Initial Values
The Initial Values node adds initial values for the flow variables that serve as an initial
condition for a transient simulation. Also specify which of the fluids initially occupies
the domain selection.
IN IT IA L VA LUES
Enter values or expressions for the initial value of the Velocity field u and for the
Pressure p. The default values are all 0.
For the Level Set and Phase Field interfaces, click the Fluid 1 button or the Fluid 2
button under Fluid initially in domain to prescribe the fluid initially in the current
domain selection. Add another Initial Values node to specify the fluid initially in
another domain.
) study is being used, for the
If the Transient with Phase Initialization (
initialization to work it is crucial that there are two Initial Value nodes and
one Initial Interface node. One of the Initial Values nodes should use Fluid
initially in domain: Fluid 1and the other Fluid initially in domain: Fluid 2. The
Initial Interface node should have all interior boundaries where the
interface is initially present as selection. If the selection of the Initial
Interface node is empty, the initialization fails. See Phase Initialization.
Initial Interface
Use the Initial Interface node to define the initial position as a boundary condition on
interior boundaries. During the initialization step, this boundary condition sets the
level set function to 0.5 or the phase field function to 0. Transient simulations of the
fluid flow treats the boundary as an interior boundary.
If the Transient with Phase Initialization (
) study is being used, it is
crucial that there are two Initial Value nodes and one Initial Interface
node for the initialization to work. One of the Initial Values nodes should
use Fluid initially in domain: Fluid 1and the other Fluid initially in domain:
Fluid 2. The Initial Interface node should have all interior boundaries
where the interface is initially present as selection. If the selection of the
Initial Interface node is empty, the initialization fails. See Phase
Initialization.
T H E L A M I N A R T W O - P H A S E F L O W , L E V E L S E T A N D L A M I N A R TW O - P H A S E F L O W , P H A S E F I E L D I N T E R F A C E S
121
This physics interface deforms the mesh within the domains on either side of the two
fluid interface to track its movement.
The Compressibility defaults to Compressible flow (Ma<0.3), which uses a compressible
formulation of the Navier-Stokes equations. Select Incompressible flow to use the
incompressible (constant density) formulation.
122 |
C H A P T E R 4 : M U L T I P H A S E F L O W , TW O - P H A S E F L O W I N T E R F A C E S
If flow is occurring at very low Reynolds numbers the inertial term in the
Navier-Stokes equations can be neglected and the linear Stokes equations can be solved
on the domain. This flow type is referred to as creeping flow or Stokes flow and can
occur in microfluidics and MEMS devices, where the flow length scales are very small.
To make this approximation select the Neglect inertial term (Stokes flow) check box,
which significantly improves the solver speed.
Enter a Reference pressure level pref (SI unit: Pa). The default value is 1[atm].
FREE DEFORMATION SETTINGS
Select a Mesh smoothing typeWinslow (the default), Laplace, Hyperelastic, or Yeoh. For
the Yeoh mesh smoothing type, also specify a Stiffening factor (default: 100). See
Smoothing Methods in the COMSOL Multiphysics Reference Manual for more
information.
FRAME SETTINGS
The Material frame coordinates setting enables the names of the space coordinates for
tracking the mesh movement in the x, y, and z directions to be changed. The default
names are the coordinates of the spatial frame in uppercase letters.
The Geometry shape order setting controls the order of polynomials used for
representing the geometry shape in the spatial frame. The same order is used for
Lagrange shape functions defining the mesh position in domains where Free
displacement has been activated.
DEPENDENT VA RIA BLES
The dependent variables (field variables) are for the Velocity field and Pressure. The
names can be changed in the corresponding fields, but the names of fields and
dependent variables must be unique within a model.
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.
T H E L A M I N A R TW O - P H A S E F L O W , M O V I N G M E S H I N T E R F A C E
123
Domain, Boundary, Edge, Point, and Pair Nodes for the Laminar
Two-Phase Flow, Moving Mesh Interface
Smoothing Methods in the COMSOL Multiphysics Reference
Manual
Domain, Boundary, Edge, Point, and Pair Nodes for the Laminar
Two-Phase Flow, Moving Mesh Interface
The Laminar Two-Phase Flow, Moving Mesh Interface has these domain, boundary,
edge, point, and pair nodes available from the Physics ribbon toolbar (Windows users),
Physics context menu (Mac or Linux users), or right-click to access the context menu
(all users)..
In general, to add a node, go to the Physics toolbar, no matter what
operating system you are using. Subnodes are available by clicking the
parent node and selecting it from the Attributes menu.
Boundary conditions (and edge and point nodes) must be specified for
the moving mesh on all external boundaries of the model. It can also be
beneficial to set boundary conditions on internal boundaries - depending
on the problem and the geometry. The boundary conditions are of two
forms: a Prescribed Mesh Displacement and a Prescribed Mesh Velocity.
124 |
C H A P T E R 4 : M U L T I P H A S E F L O W , TW O - P H A S E F L O W I N T E R F A C E S
Navier Slip
Fluid-Fluid Interface
Fixed Mesh
Free Deformation
Wall-Fluid Interface
Initial Values
These nodes are described for the Moving Mesh interface in the COMSOL
Multiphysics Reference Manual (listed in alphabetical order):
Prescribed Deformation
These nodes are described for the Laminar Flow interface (listed in alphabetical order):
No Viscous Stress
Flow Continuity
Fluid Properties
Symmetry
Inlet
Volume Force
Open Boundary
Wall
Outlet
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.
T H E L A M I N A R TW O - P H A S E F L O W , M O V I N G M E S H I N T E R F A C E
125
Free Deformation
The Free Deformation default node applies the equations for the moving mesh to the
domains selected, to allow deformation of the mesh. Domains either side of a
fluid-fluid interface must have this option selected.
INITIAL DEFORMATION
Enter values or expressions for each coordinate based on the space dimension for the
Initial mesh displacement (dx0, dy0, dz0) (SI unit: m). The default is for no mesh
displacement.
126 |
C H A P T E R 4 : M U L T I P H A S E F L O W , TW O - P H A S E F L O W I N T E R F A C E S
CONSTRAINT SETTINGS
Initial Values
The Initial Values nod adds initial values for the velocity field and the pressure that can
serve as an initial condition for a transient simulation or as an initial guess for a
nonlinear solver.
IN IT IA L VA LUES
Enter values or expressions for the initial value of the Velocity field u (SI unit: m/s) and
for the Pressure p (SI unit: Pa). The default values are 0.
Fixed Mesh
The Fixed Mesh node specifies that the mesh does not move in the selected domains. A
fixed mesh should only be used on domains that are not adjacent to a fluid-fluid
interface.
All the elements on the domain boundary are also fixed. Care should be
taken to prevent mesh movement in an undesirable manner in an adjacent
free deformation domain.
T H E L A M I N A R TW O - P H A S E F L O W , M O V I N G M E S H I N T E R F A C E
127
(partially) warped inside-out. In this case introducing extra boundaries with explicit
deformation inside the domains can help. Using a quadrilateral mesh (which is
typically stiffer as it deforms) is another possibility. Also generate a new mesh for the
region covered by the deformed mesh and let the solver continue by deforming the
new mesh; see Remeshing a Deformed Mesh in the COMSOL Multiphysics
Reference Manual.
When using Geometry shape order larger than 1 in the Moving Mesh and Deformed
Geometry interfaces, the mesh moving techniques often produce elements with
distorted shapes. The default value of 1 is appropriate for most situations, provided
a sufficiently fine mesh is employed.
The measure of mesh quality does not capture these distorted shapes
because it is computed from the positions of the corners of the mesh
element (ignoring midsize nodes, for example).
P RE S C R I B E D M E S H VE L O C I T Y
By default the Prescribed x velocity, Prescribed y velocity, and Prescribed z velocity check
boxes are selected. The available fields are based on the space dimension of the model.
Enter values or expressions for each of the Prescribed mesh velocity fields vx, vy, and vz
(SI unit: m/s). Click to clear the check boxes as needed.
CONSTRAINT SETTINGS
Fluid-Fluid Interface
The Fluid-Fluid Interface node defines the initial position of a fluid-fluid interface and
includes equations to track the evolution of the interface. The Wall-Fluid Interface
subnode is available from the context menu (right-click the parent node) or from the
Physics toolbar, Attributes menu.
128 |
C H A P T E R 4 : M U L T I P H A S E F L O W , TW O - P H A S E F L O W I N T E R F A C E S
S U R F A C E TE N S I O N
The default Surface tension coefficient (SI unit: N/m) is User defined. Or select Library
coefficient, liquid/gas interface or Library coefficient, liquid/liquid interface.
For Library coefficient, liquid/gas interface select an option from the listWater/Air,
Acetone/Air, Acetic acid/Air, Ethanol/Air, Ethylene glycol/Ethylene glycol vapor, Diethyl
ether/Air, Glycerol/Air, Heptane/Nitrogen, Mercury/Mercury vapor, or Toluene/Air.
For Library coefficient, liquid/liquid interface select an option from the list
Benzene/Water, 20C, Corn oil/Water, 20C, Ether/Water, 20C, Hexane/Water, 20C,
Mercury/Water, 20C, or Olive oil/Water, 20C.
MASS FLUX
The mass flux setting specifies the mass transfer across the boundary, due to processes
such as boiling. The default Mass Flux Mf (SI unit: kg/(m2 s)) is User defined, with a
value of 0.
The default Surface tension coefficient (SI unit: N/m) is User defined. Or select Library
coefficient, liquid/gas interface or Library coefficient, liquid/liquid interface.
For Library coefficient, liquid/gas interface select an option from the listWater/Air,
Acetone/Air, Acetic acid/Air, Ethanol/Air, Ethylene glycol/Ethylene glycol vapor, Diethyl
ether/Air, Glycerol/Air, Heptane/Nitrogen, Mercury/Mercury vapor, or Toluene/Air.
For Library coefficient, liquid/liquid interface select an option from the list
Benzene/Water, 20C, Corn oil/Water, 20C, Ether/Water, 20C, Hexane/Water, 20C,
Mercury/Water, 20C, or Olive oil/Water, 20C.
T H E L A M I N A R TW O - P H A S E F L O W , M O V I N G M E S H I N T E R F A C E
129
Wall-Fluid Interface
The Wall Fluid Interface subnode is available from the context menu (right-click the
Fluid-Fluid Interface and External Fluid Interface parent node) or from the Physics
toolbar, Attributes menu. This feature applies the forces necessary to maintain the
contact angle at the wall.
The Wall-Fluid Interface node should only be used with the Navier Slip
boundary condition, which should be applied to the adjacent wall(s). If
another boundary condition is used with this feature the normal for the
wall is not correct on the wall-fluid interface, leading to an incorrect force
on the contact line.
WA LL -F LU ID IN TE RFA CE
Select an option from the Specify contact angle listDirectly (the default) or Through
Youngs equation.
For Directly enter a Contact angle w (SI unit: rad). The default is /2 radians.
For Through Youngs equation enter values or expressions for Phase 1-Solid surface
energy density s1 (SI unit: J/m2) and Phase 2-Solid surface energy density s2 (SI
unit: J/m2).
The contact angle w is defined between the fluid-fluid interface and the
surface of the wall adjacent to phase 1.
Navier Slip
The Navier Slip boundary condition is suitable for walls in contact with the fluid-fluid
interface. If this boundary condition is used in combination with a Prescribed
Displacement perpendicular to the wall only, the fluid-fluid interface can move along
the wall. For applications where the interface is fixed on the wall, the Wall-No slip
condition together with a zero Prescribed Displacement both parallel and perpendicular
to the wall can be used.
The boundary condition enforces the slip condition u nwall = 0 and adds a frictional
force of the form
130 |
C H A P T E R 4 : M U L T I P H A S E F L O W , TW O - P H A S E F L O W I N T E R F A C E S
F fr = --- u
where is the slip length. By default = h/5, where h is the mesh element size. Strictly
speaking the slip length should be set to a fixed physically motivated value (which is
typically a few 10 s of nm) and the mesh should be fine enough to resolve distances of
the order of the slip length. For many practical purposes the above approximation is
sufficient. This boundary condition does not set the tangential velocity component to
zero; however, the extrapolated tangential velocity component is 0 at the distance
outside the wall (this is illustrated for the wetted wall boundary condition in
Figure 4-2).
For axisymmetric components when the Swirl flow property is active, there is an option
to set a value for the Out-of-plane wall velocity, vw. In this case, an additional constraint
is imposed on the component of the velocity:
u = vw
NAVIER SLIP
T H E L A M I N A R TW O - P H A S E F L O W , M O V I N G M E S H I N T E R F A C E
131
u
T
+ ( u )u = [ p I + ( u + u ) ] + F g + F st + F ext + F
t
u = 0
(4-3)
(4-4)
If the level set method is used to track the interface, it adds the following equation:
----+ u = ( 1 ) ----------
(4-5)
132 |
C H A P T E R 4 : M U L T I P H A S E F L O W , TW O - P H A S E F L O W I N T E R F A C E S
= 1 + ( 2 1 )
and the dynamic viscosity is given by
= 1 + ( 2 1 )
where 1 and 2 are the constant densities of Fluid 1 and Fluid 2, respectively, and 1
and 2 are the dynamic viscosities of Fluid 1 and Fluid 2, respectively. Here, Fluid 1
corresponds to the domain where < 0.5 , and Fluid 2 corresponds to the domain
where > 0.5 .
Further details of the theory for the level set method are in Ref. 1.
USING THE PHASE FIELD METHOD
If the phase field method is used to track the interface, it adds the following equations:
+ u = -----2-
t
(4-6)
2
f ext
2
2
= + ( 1 ) + -----
(4-7)
where the quantity (SI unit: N) is the mixing energy density and (SI unit: m) is a
capillary width that scales with the thickness of the interface. These two parameters are
related to the surface tension coefficient, (SI unit: N/m), through the equation
2 2
= ----------- --3
and is related to through =2 where is the mobility tuning parameter (set to 1
by default). The volume fraction of Fluid 2 is computed as
V f = min ( max ( [ ( 1 + ) 2 ], 0 ), 1 )
where the min and max operators are used so that the volume fraction has a lower limit
of 0 and an upper limit of 1. The density is then defined as
= 1 + ( 2 1 )V f
and the dynamic viscosity according to
= 1 + ( 2 1 )V f
133
where 1 and 2 are the densities and 1 and 2 are the dynamic viscosities of Fluid 1
and Fluid 2, respectively.
The mean curvature (SI unit: 1/m) can be computed by entering the following
expression:
G
= 2 ( 1 + ) ( 1 ) ---
where G is the chemical potential defined as:
( 1 ) f
G = 2 + -------------------+ -----
2
Details of the theory for the phase field method are in Ref. 2.
F O R C E TE R M S
The four forces on the right-hand side of Equation 4-3 are due to gravity, surface
tension, a force due to an external contribution to the free energy (using the phase field
method only), and a user-defined volume force.
F st = ( ( I ( nn ) ) )
For a derivation of this formulation, see Appendix A in Ref. 3. In the weak formulation
of the momentum equation, it is possible to move the divergence operator, using
integration by parts, to the test functions for the velocity components.
The -function is approximated by a smooth function according to
= 6 ( 1 )
134 |
C H A P T E R 4 : M U L T I P H A S E F L O W , TW O - P H A S E F L O W I N T E R F A C E S
f
F st = G ------
where G is the chemical potential (SI unit: J/m3) defined in The Equations for the
Phase Field Method and f is a user-defined source of free energy.
Phase Initialization
If the study type Transient with Phase Initialization is used in the model, the level set
or phase field variable is automatically initialized. For this study, two study steps are
created, Phase Initialization and Time Dependent. The Phase Initialization step solves
for the distance to the initial interface >, Dwi. The Time Dependent step then uses the
initial condition for the level set function according to the following expression:
1 0 = ----------------------D
1 + e wi
in domains initially filled with Fluid 1 and
135
1
0 = ------------------------D
1 + e wi
in domains initially filled with Fluid 2.
Correspondingly, for the phase field method the following expressions are used:
D wi
0 = tanh ----------
2
in Fluid 1 and
D wi
0 = tanh ----------
2
in Fluid 2. The initial condition for the help variable is 0 = 0. These expressions are
based on the analytical solution of the steady state solution of Equation 4-5,
Equation 4-6, and Equation 4-7 for a straight, non-moving interface.
For the initialization to work it is crucial that there are two Initial Values
nodes and one Initial Interface node. One of the Initial Values nodes
should use Fluid initially in domain: Fluid 1and the other Fluid initially in
domain: Fluid 2. The Initial Interface node should have all interior
boundaries where the interface is initially present as selection. If the
selection of the Initial Interface node is empty, the initialization fails.
Numerical Stabilization
Four types of stabilization methods are available for the flow (Navier-Stokes) and
interface (level set or phase field) equations. Two are consistent stabilization
methodsStreamline diffusion and Crosswind diffusionand two are inconsistent
Isotropic diffusion and Anisotropic diffusion.
136 |
C H A P T E R 4 : M U L T I P H A S E F L O W , TW O - P H A S E F L O W I N T E R F A C E S
137
138 |
C H A P T E R 4 : M U L T I P H A S E F L O W , TW O - P H A S E F L O W I N T E R F A C E S
(4-8)
The dynamic viscosity (SI unit: Pas) is allowed to depend on the thermodynamic
state but not on the velocity field.
y = y ( X , Y, t )
The original, undeformed, mesh is referred to as the material frame (or reference
frame) whilst the deformed mesh is called the spatial frame. COMSOL Multiphysics
also defines geometry and mesh frames, which are coincident with the material frame
for this physics interface.
For the Two-Phase Flow, Moving Mesh interface the fluid flow equations (along with
other coupled equations such as electric fields or chemical species transport) are solved
in the spatial frame in which the mesh is perturbed. The movement of the phase
boundary is therefore accounted for in these interfaces.
T H E O R Y F O R T H E TW O P H A S E F L O W M O V I N G M E S H I N T E R F A C E
139
The boundary conditions applied at an interface between two immiscible fluids, fluid 1
and fluid 2 (see Figure 4-2), are given by (Ref. 1):
1
1
u 1 = u 2 + ------ ------ M f n i
1 2
(4-9)
n i 2 = n i 1 + f st
(4-10)
Mf
u mesh = u 1 n i --------- n i
(4-11)
where u1 and u2 are the velocities of the fluids 1 and 2 respectively, umesh is the
velocity of the mesh at the interface between the two fluids, ni is the normal of the
interface (outward from the domain of fluid 1), 1 and 2 are the total stress tensors in
domains 1 and 2 respectively, fst is the force per unit area due to the surface tension
and Mf is the mass flux across the interface (SI unit: kg/(m2s)).
Figure 4-2: Definition of fluid 1 and fluid 2 and the interface normal.
140 |
C H A P T E R 4 : M U L T I P H A S E F L O W , TW O - P H A S E F L O W I N T E R F A C E S
The tangential components of Equation 4-9 enforce a no-slip condition between the
fluids at the boundary. In the absence of mass transfer across the boundary,
Equation 4-9 and Equation 4-11 ensure that the fluid velocity normal to the boundary
is equal to the velocity of the interface. When mass transfer occurs these equations
result from conservation of mass and are easily derived in the frame where the
boundary is stationary.
The components of the total stress tensor, uv, represent the uth component of the
force per unit area perpendicular to the v-direction. n=nv uv (using the summation
convention) is therefore interpreted as the force per unit area acting on the boundary
- in general this is not normal to the boundary. Equation 4-10 therefore expresses the
force balance on the interface between the two fluids.
Two boundary conditions (Equation 4-9 and Equation 4-10) are
necessary to couple the two domains as there are two separate sets of
Navier-Stokes equations, one set for each of the domains.
The force due to the surface tension is given by the following expression:
f st = ( s n i )n i s
(4-12)
F n = ---R
The quantity sni in the first term on the right-hand side of Equation 4-12 is related
to the mean curvature, , of the surface by the equation =sni. In two dimensions
the mean curvature = 1/R so sni = 1/R. The first term in Equation 4-12 is
therefore the normal force per unit area acting on the boundary due to the surface
tension.
T H E O R Y F O R T H E TW O P H A S E F L O W M O V I N G M E S H I N T E R F A C E
141
The tangential force per unit area, Ft, acting on the interface in Figure 4-3 can be
obtained from the force balance along the direction of s in the limit s 0:
F t s = ( + ) cos ( ) ( + )
( + )
F t ----------------------------s
This is equivalent to the second term on the right of Equation 4-12.
To obtain additional insight into the boundary condition, it is helpful to re-write
Equation 4-10 as
T
n i ( ( p 1 p 2 )I 1 ( u 1 ( u 1 ) ) + 2 ( u 2 ( u 2 ) ) ) = ( s n i )n i s
assuming Newtonian fluids with viscosities 1 and 2 for fluids 1 and 2 respectively and
that p1 and p2 are the pressures in the respective fluids adjacent to the boundary. This
equation expresses the equality of two vector quantities. It is instructive to consider the
components perpendicular and tangential to the boundary. In the direction of the
boundary normal
T
p 1 p 2 + n i ( 2 ( u 2 ( u 2 ) ) 1 ( u 1 ( u 1 ) ) ) n i = ( s n i )n(4-13)
i
whereas in the tangential direction, ti,
T
n i ( 2 ( u 2 ( u 2 ) ) 1 ( u 1 ( u 1 ) ) ) t i = s
142 |
C H A P T E R 4 : M U L T I P H A S E F L O W , TW O - P H A S E F L O W I N T E R F A C E S
(4-14)
It is often the case that the viscosity of fluid 1 is significantly greater than that of fluid 2
(for example for a liquid-vapor interface). In this case, the viscosity terms in the total
stress from fluid 2 can be neglected and Equation 4-10 becomes:
n i 1 = p ext n i + f st
(4-15)
The outer fluid (fluid 2) now enters the equation system only through the pressure
term and the system can be represented by a domain consisting solely of the fluid 1
domain with an expression (or constant value) for the external pressure, pext, in the
fluid 2 domain. Since the velocity in fluid 2 does not affect fluid 1, fluid 2 does not
need to be explicitly modeled and Equation 4-9 can be dropped. Equation 4-11 is
retained with Mf=0.
T H E O R Y F O R T H E TW O P H A S E F L O W M O V I N G M E S H I N T E R F A C E
143
Figure 4-4: Forces per unit length acting on a fluid-fluid interface at a three phase
boundary with a solid wall. The surface tension force per unit length, , is balanced by a
reaction force per unit length at the surface, Fn, and by the forces generated by the surface
energies of the two phases at the interface: s1 and s2.
In equilibrium, the surface tension forces and the normal restoring force from the
surface are in balance at a constant contact angle (c), as shown in Figure 4-4. This
equilibrium is expressed by Young's equation, which considers the components of the
forces in the plane of the surface:
cos ( c ) + s1 = s2
(4-16)
where is the surface tension force between the two fluids, s1 is the surface energy
density on the fluid 1solid interface and s2 is the surface energy density on the
fluid 2solid interface.
There is still debate in the literature as to precisely what occurs in non-equilibrium
situations (for example, drop impact) when the physical contact angle deviates from
the contact angle specified by a simple application of Young's equation. A simple
approach, is to assume that the unbalanced part of the in plane Young Force acts on
the fluid to move the contact angle towards its equilibrium value (Ref. 2). COMSOL
144 |
C H A P T E R 4 : M U L T I P H A S E F L O W , TW O - P H A S E F L O W I N T E R F A C E S
Figure 4-5: Diagram showing normal and tangential vectors defined on the interface
between two fluids and a solid surface in 3D (left) and 2D (right). The following vectors
are defined: ni, the fluid-fluid interface unit normal; ns, the solid surface unit normal;
tc, the tangent to the three-phase contact line; and mi and ms, the two unit binormals,
which are defined as mi=tcni and ms=tcns, respectively.
T H E O R Y F O R T H E TW O P H A S E F L O W M O V I N G M E S H I N T E R F A C E
145
3. W. Ren and E. Weinan, Boundary Conditions for the Moving Contact Line
Problem, Physics of Fluids, vol. 19, p. 022101, 2007.
4. W. Ren and D. Hu, Continuum Models for the Contact Line Problem, Physics of
Fluids, vol. 22, p. 102103, 2010.
146 |
C H A P T E R 4 : M U L T I P H A S E F L O W , TW O - P H A S E F L O W I N T E R F A C E S
147
T he L e v e l S e t In t erfac e
The Level Set (ls) interface (
), found under the Mathematics>Moving Interface
branch (
) when adding an interface, is used to track moving interfaces in fluid-flow
models by solving a transport equation for the level set function. Simulations using the
Level Set interface are always time dependent since the position of an interface almost
always depends on its history.
The main node is the Level Set Model feature, which adds the level set equation and
provides an interface for defining the level set properties and the velocity field.
When this physics interface is added, the following default nodes are also added in the
Model BuilderLevel Set Model, No Flow (the default boundary condition) and Initial
Values. Then, from the Physics toolbar, add other nodes that implement, for example,
boundary conditions. You can also right-click Level Set to select physics features from
the context menu.
SETTINGS
The dependent variable (field variable) is the Volume fraction of fluid 2 phi. The name
can be changed but the names of fields and dependent variables must be unique within
a model.
Conservative and Non-Conservative Form
Domain, Boundary, and Pair Nodes for the Level Set Interface
Theory for the Level Set Interface
Theory for the Level Set and Phase Field Interfaces
148 |
Domain, Boundary, and Pair Nodes for the Level Set Interface
The Level Set Interface has the following domain, boundary and pair nodes described.
Initial Interface
No Flow
Initial Values
Outlet1
Inlet
Symmetry1
Boundary conditions for axial symmetry boundaries are not required. For
the symmetry axis at r = 0, the software automatically provides a suitable
boundary condition and adds an Axial Symmetry node that is valid on the
axial symmetry boundaries only.
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.
+ u = ( 1 ) ----------
t
and provides the options to define the associated level set parameters and the velocity
field.
LEVEL SET PARAMETERS
Enter a value or expression for the Reinitialization parameter (SI unit: m/s). The
default is 1 m/s.
149
Enter a value or expression for the Parameter controlling interface thickness els
(SI unit: m). The default expression is [Link]/2, which means that the value is half
of the maximum mesh element size in the region through which the interface passes.
CONVECTION
Enter values or expressions for the components (u, v, and w in 3D, for example) of the
Velocity field u (SI unit: m/s). The applied velocity field transports the level set
function through convection.
Initial Values
Use the Initial Values node to set up the initial conditions.
INITIAL VALUES
In order to be able to use the automatic initialization, select an option from the Domain
initially: listInside interface (the default), Outside interface, or Specify level set function
explicitly if automatic initialization is not required.
) study is being used, for the
If the Transient with Phase Initialization (
initialization to work it is crucial that there are two Initial Values nodes
and one Initial Interface node. One of the Initial Values nodes should use
Domain initially: Inside interface and the other Domain initially: Outside
interface. The Initial Interface node should have all interior boundaries
where the interface is initially present as selection. If the selection of the
Initial interface node is empty, the initialization fails.
See Initializing the Level Set Function.
Inlet
The Inlet node adds a boundary condition for inlets (inflow boundaries). At inlets a
value of the level set function must be specified. Typically set to either 0 or 1.
SETTINGS
Enter a value for the Level set function value . The value must be in the range from 0
to 1, and the default is 0.
150 |
Initial Interface
The Initial Interface node defines the boundary as the initial position of the interface
= 0.
If the Transient with Initialization (
) study is being used, for the
initialization to work it is crucial that there are two Initial Values nodes
and one Initial Interface node. One of the Initial Values nodes should use
Domain initially: Inside interface and the other Domain initially: Outside
interface. The Initial Interface node should have all interior boundaries
where the interface is initially present as selection. If the selection of the
Initial interface node is empty, the initialization fails.
See Initializing the Level Set Function.
No Flow
The No Flow node adds a boundary condition that represents boundaries where there
is no flow across the boundary. This is the default boundary condition.
151
This interface defines the dependent variables (fields) Phase field variable and Phase
field help variable . If required, edit the name, but dependent variables must be
unique within a model.
Conservative and Non-Conservative Forms
Domain, Boundary, and Pair Nodes for the Phase Field Interface
Theory for the Phase Field Interface
152 |
Domain, Boundary, and Pair Nodes for the Phase Field Interface
Compared to the single-phase flow interfaces, the two-phase flow
interfaces include two additional boundary conditionsthe Wall
boundary conditions Wetted wall and Moving wetted wall.
The Phase Field Interface includes the following domain, boundary, and pair nodes,
listed in alphabetical order, available from the Physics ribbon toolbar (Windows users),
Physics context menu (Mac or Linux users), or right-click to access the context menu
(all users).
In general, to add a node, go to the Physics toolbar, no matter what
operating system you are using. Subnodes are available by clicking the
parent node and selecting it from the Attributes menu.
Initial Interface
Initial Values
Symmetry1
Inlet
Wetted Wall
Outlet1
1
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.
153
Define the following phase field parameters. Enter a value or expression for the:
Surface tension coefficient (SI unit: N/m).
Parameter controlling interface thickness epf (SI unit: m). The default expression is
[Link]/2, which means that the value is half of the maximum mesh element size
in the region through which the interface passes.
Mobility tuning parameter (SI unit: ms/kg). The default is 1 ms/kg, which is a
good starting point for most models. This parameter determines the time scale of
the Cahn-Hilliard diffusion and it thereby also governs the diffusion-related time
scale for the interface.
Keep the parameter value large enough to maintain a constant interface
thickness but still low enough to not damp the convective motion. A too
high mobility can also lead to excessive diffusion of droplets.
EXTERNAL FREE ENERGY
Add a source of external free energy to the phase field equations. This modifies the last
term on the right-hand side of the equation:
2
2
2
f
= + ( 1 ) + -----
The external free energy f (SI unit: J/m3) is a user-defined free energy. In most cases,
the external free energy can be set to zero. Manually differentiate the expression for
the external free energy with respect to and then enter it into the -derivative of
external free energy field f .
CONVECTION
Enter values or expressions for the components (u, v, and w in 3D, for example) of the
Velocity field u (SI unit: m/s). The applied velocity field transports the phase field
variables through convection.
154 |
Initial Values
The Initial Values node adds initial values for the phase field variable and the phase field
help variable that can serve as initial conditions for a transient simulation.
IN IT IA L VA LUES
Enter initial values or expressions for the Phase field variable and the Phase field help
variable . The default values are 0.
Inlet
The Inlet feature node adds a boundary condition for inlets (inflow boundaries). At
inlets a volume fraction Vf must be specified, typically either 0 or 1. Mathematically
this boundary condition imposes
= 2V f 1 and n -----2- = 0
INLET
Specify a value of the Inlet volume fraction Vf. The value must be in the range from 0
to 1 and the default is 0.
Initial Interface
The Initial Interface node defines the boundary as the initial position of the interface
= 0.
If the Transient with Phase Initialization (
) study is being used, for the
initialization to work it is crucial that there are two Initial Values nodes and
one Initial Interface node. One of the Initial Values nodes is set to phipf
= 1 and the other to phipf = -1. The Initial Interface node should have
all interior boundaries where the interface is initially present as selection.
If the selection of the Initial Interface node is empty, the initialization
fails.
155
Wetted Wall
The Wetted Wall node is the default boundary condition representing wetted walls.
Along a wetted wall the contact angle for the fluid, w, is specified, and across it the
mass flow is zero. This is prescribed by
2
n = cos ( w )
in combination with
n -----2- = 0
WE TTE D WA LL
Enter a value or expression for the Contact angle w. The default value is /2 rad.
156 |
157
Figure 5-1: An example of two domains divided by an interface. In this case, one of the
domains consists of two parts. Figure 5-2 shows the corresponding level set representation.
Figure 5-2: A surface plot of the level set function corresponding to Figure 5-1.
The physics interface solves Equation 5-1 in order to move the interface with the
velocity field u:
+ u = ( 1 ) ----------
t
(5-1)
The terms on the left-hand side give the correct motion of the interface, while those
on the right-hand side are necessary for numerical stability. The parameter, ,
determines the thickness of the region where varies smoothly from zero to one and
is typically of the same order as the size of the elements of the mesh. By default, is
constant within each domain and equals the largest value of the mesh size, h, within
158 |
(5-2)
the volume (area for 2D problems) bounded by the interface should be conserved if
there is no inflow or outflow through the boundaries. To obtain exact numerical
conservation, switch to the conservative form
+ ( u ) = ( 1 ) ----------
(5-3)
159
1 0 = ----------------------D
1 + e wi
in domains initially inside the interface. Here, inside refers to domains where <0.5
and outside refers to domains where >0.5.
For the initialization to work it is crucial that there are two Initial Values
nodes and one Initial Interface node. One of the Initial Values nodes
should use Domain initially: Inside interface and the other Domain initially:
Outside interface. The Initial Interface node should have all interior
boundaries where the interface is initially present as selection. If the
selection of the Initial interface node is empty, the initialization fails.
n = ---------
(5-4)
= 0.5
= 0.5
(5-5)
These variables are available in the physics interface as the interface normal and mean
curvature.
160 |
161
162 |
F ( , , T ) =
--2-
1
2
+ f ( , T ) dV =
ftot dV
where is a measure of the interface thickness. Equation 5-6 describes the evolution
of the phase field parameter:
f tot
f tot
+ ( u ) =
(5-6)
where ftot (SI unit: J/m3) is the total free energy density of the system, and u
(SI unit: m/s) is the velocity field for the advection. The right-hand side of
Equation 5-6 aims to minimize the total free energy with a relaxation time controlled
by the mobility (SI unit: m3s/kg).
The free energy density of an isothermal mixture of two immiscible fluids is the sum
of the mixing energy and elastic energy. The mixing energy assumes the
Ginzburg-Landau form:
2
2
2
1
f mix ( , ) = --- + --------2- ( 1 )
2
4
where is the dimensionless phase field variable, defined such that the volume fraction
of the components of the fluid are (1+ )/2 and (1 )/2. The quantity
(SI unit: N) is the mixing energy density and (SI unit: m) is a capillary width that
scales with the thickness of the interface. These two parameters are related to the
surface tension coefficient, (SI unit: N/m), through the equation
2 2
= ----------- --3
(5-7)
The PDE governing the phase field variable is the Cahn-Hilliard equation:
+ u = G
t
(5-8)
where G (SI unit: Pa) is the chemical potential and (SI unit: m3s/kg) is the mobility.
The mobility determines the time scale of the Cahn-Hilliard diffusion and must be
large enough to retain a constant interfacial thickness but small enough so that the
convective terms are not overly damped. In COMSOL Multiphysics the mobility is
determined by a mobility tuning parameter that is a function of the interface thickness
= 2. The chemical potential is:
163
2
( 1)
G = + ---------------------2
(5-9)
+ u = -----2-
t
(5-10)
= + ( 1 )
(5-11)
+ u = -----2-
t
Using the conservative phase field form, exact numerical conservation of the integral
of is obtained. However, the non-conservative form is better suited for numerical
calculations and usually converges more easily. The non-conservative form, which is
the default form, only conserves the integral of the phase field function approximately,
but this is sufficient for most applications.
f
2
2
= + ( 1 ) + -----
(5-12)
164 |
165
G = ------2
and the surface tension force F = G .
The mean curvature (SI unit: 1/m) of the interface can be computed by entering the
following expression:
G
= 2 ( 1 + ) ( 1 ) ---
166 |
167
Enter a Reference pressure level pref (SI unit: Pa). The default value is 1[atm].
DEPENDENT VARIABLES
The dependent variable (field variable) is the Pressure. The name can be changed but
the names of fields and dependent variables must be unique within a model.
168 |
DISCRETIZATION
The Compute boundary fluxes check box is not activated by default. 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, COMSOL 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.
Also the Apply smoothing to boundary fluxes check box is available if the previous check
box is checked. The smoothing can provide a better 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 transport 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.
Domain, Boundary, Edge, Point, and Pair Nodes for the Darcys Law
Interface
Theory for the Darcys Law Interface
Physical Constants in the COMSOL Multiphysics Reference Manual
Domain, Boundary, Edge, Point, and Pair Nodes for the Darcys Law
Interface
The Darcys Law Interface has the following domain, boundary, edge, point, and pair
nodes, These nodes available from the Physics ribbon toolbar (Windows users), Physics
T H E D A R C Y S L A W I N T E R F A C E
169
context menu (Mac or Linux users), or right-click to access the context menu (all
users).
In general, to add a node, go to the Physics toolbar, no matter what
operating system you are using. Subnodes are available by clicking the
parent node and selecting it from the Attributes menu.
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.
DOMAIN
Initial Values
Mass Source
B O U N D A R Y, E D G E , A N D PO I N T
The following nodes (listed in alphabetical order) are available on exterior boundaries:
The relevant physics interface condition at interior boundaries is continuity:
n ( 1 u1 2 u2 ) = 0
The continuity boundary condition ensures that the pressure and mass flux are
continuous. In addition, the Pressure boundary condition is available on interior
boundaries.
170 |
(6-1)
u = --- p
(6-2)
FLUID PROPERTIES
Select the Fluid material to use for the fluid properties. Select Domain material (the
default) to use the material defined for the domain. Select another material to use that
materials properties for the fluid.
Density
The default Density (SI unit: kg/m3) uses values From material based on the Fluid
material selection.
For User defined enter another value or expression. The default is 0 kg/m3.
For Ideal gas it uses the ideal gas law to describe the fluid. In this case, specify the
thermodynamics properties. Select a Gas constant typeSpecific gas constant Rs (the
default) or Mean molar mass Mn (SI unit: J/(molK)). For Mean molar mass the
universal gas constant R = 8.314 J/(molK) is used as the built-in physical constant.
For both properties, the defaults use values From material. For User defined enter
another value or expression.
Dynamic Viscosity
Select a Dynamic viscosity (SI unit: Pas). The default uses values From material as
defined by the Fluid material selected. For User defined the default is 0 Pas.
MATRIX PROPERTIES
Select the material to use as porous matrix. Select Domain material from the Porous
material list (the default) to use the material defined for the porous domain. Select
another material to use that materials properties.
The default Porosity p (a dimensionless number between 0 and 1) uses the value From
material, defined by the Porous material selected. For User defined the default is 0.
T H E D A R C Y S L A W I N T E R F A C E
171
The default Permeability (SI unit: m2) uses the value From material, as defined by the
Porous material selected. For User defined select Isotropic to define a scalar value or
Diagonal, Symmetric, or Anisotropic to define a tensor value and enter another value or
expression in the field or matrix.
Mass Source
The Mass Source node adds a mass source Qm, which appears on the right-hand side of
the Darcys Law equation (Equation 6-8, the equation for porosity).
( ) + ( u )
= Qm
t
(6-3)
MASS SOURCE
Enter a value or expression for the Mass source Qm (SI unit: kg/(m3s)). The default is
0 kg/(m3s).
Initial Values
The Initial Values node adds an initial value for the pressure that can serve as an initial
condition for a transient simulation or as an initial guess for a nonlinear solver.
INITIAL VALUES
Enter a value or expression for the initial value of the Pressure p (SI unit: Pa). The
default value is 0 Pa.
Pressure
Use the Pressure node to specify the pressure on a boundary. In many cases the
distribution of pressure is known, giving a Dirichlet condition p = p0 where p0 is a
known pressure given as a number, a distribution, or an expression involving time, t,
for example.
PRESSURE
Enter a value or expression for the Pressure p0(SI unit: Pa). Enter a relative pressure
value in p0 (SI unit: Pa).
CONSTRAINT SETTINGS
172 |
Mass Flux
Use the Mass Flux node to specify the mass flux into or out of the model domain
through some of its boundaries. It is often possible to determine the mass flux from
the pumping rate or from measurements. With this boundary condition, positive
values correspond to flow into the model domain:
n --- p = N 0
where N0 is a value or expression for the specified inward (or outward) Darcy flux.
MASS FLUX
Enter a value or expression for the Inward mass flux N0. A positive value of N0
represents an inward mass flux whereas a negative value represents an outward mass
flux. The units are based on the geometric entity: Boundaries: (SI unit: kg/(m2s)),
Edges (SI unit: kg/(ms), and Points (SI unit: kg/s)).
Inlet
The Inlet node adds a boundary condition for the inflow (or outflow) perpendicular
(normal) to the boundary:
n --- p = U 0
where U0 is a value or expression for the specified inward (or outward) Darcy velocity.
A positive value of the velocity U0 corresponds to flow into the model domain whereas
a negative value represents an outflow.
INLET
Enter a value or expression for the Normal inflow velocity U0 (SI unit: m/s). A positive
value of U0 represents an inflow velocity. A negative value represents an outflow
velocity.
Symmetry
The Symmetry node describes a symmetry boundary. The following condition
implements the symmetry condition on an axis or a flow divide:
n --- p = 0
T H E D A R C Y S L A W I N T E R F A C E
173
No Flow
The No Flow node is the default boundary condition stating that there is no flow across
impervious boundaries. The mathematical formulation is:
n --- p = 0
Flux Discontinuity
Use the Flux Discontinuity node to specify a mass flux discontinuity through an interior
boundary. The condition is represented by the following equation:
n ( u 1 u 2 ) = N 0
In this equation, n is the vector normal (perpendicular) to the interior boundary, is
the fluid density, u1 and u2 are the Darcys velocities in the adjacent domains (as
defined in Equation 6-7) and N0 is a specified value or expression for the flux
discontinuity.
u = --- p
(6-4)
Enter a value or expression for the Inward mass flux N0 (SI unit: kg/(m2s)). A positive
value of N0 represents a mass flux discontinuity in the opposite direction to the normal
vector of the interior boundary.
Outlet
The Outlet node adds a boundary condition for the outflow (or inflow) perpendicular
(normal) to the boundary:
174 |
n --- p = U 0
where U0 is a specified value or expression for the outward (or inward) Darcy velocity.
A positive value of the velocity U0 corresponds to flow out of the model domain
whereas a negative value represents an inflow.
OUTLET
Enter a value or expression for the Normal outflow velocity U0 (SI unit: m/s). A positive
value of U0 represents an outflow velocity whereas a negative value represents an
inflow velocity.
T H E D A R C Y S L A W I N T E R F A C E
175
2
1-
T
--- ( u + ( u ) ) --- ( u )I
3
p
SETTINGS
176 |
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 br.
PHYSICAL MODEL
This node specifies the properties of the Brinkman Equations interface, which describe
the overall type of fluid flow model.
Compressibility
By default the physics interface uses the Incompressible flow formulation of the
Brinkman equations to model constant density flow. Alternatively, select Compressible
flow (Ma<0.3) from the Compressibility list if there are small variations in the density,
typically dependent on the temperature (non-isothermal flow). For compressible flow
modeled with the Brinkman Equations interface, the Mach number must be below
0.3.
The following dependent variables (fields) are defined for this physics interfacethe
Velocity field u (SI unit: m/s) and its components, and the Pressure p (SI unit: Pa).
Domain, Boundary, Point, and Pair Nodes for the Brinkman Equations
Interface
Pseudo Time Stepping for Laminar Flow Models and Pseudo Time
Stepping in the COMSOL Multiphysics Reference Manual
Theory for the Brinkman Equations Interface
177
Mass Source
Forchheimer Drag
Volume Force
Initial Values
Fluid Properties
The following nodes(listed in alphabetical order) are described for the Laminar Flow
interface :
Flow Continuity
Inlet
No Viscous Stress
Symmetry
Outlet
Wall
Open Boundary
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.
178 |
The default Fluid material uses the Domain material (the material defined for the
domain). Select another material as needed.
Both the default Density (SI unit: kg/m3) and Dynamic viscosity (SI unit: Pas) use
values From material based on the Fluid material selection. For User defined enter
another value or expression. In this case, the default is 0 kg/m3 for the density and
0 Pas for the dynamic viscosity. The dynamic viscosity describes the relationship
between the shear stresses and the shear rate in a fluid. Intuitively, water and air have
a low viscosity, and substances often described as thick, such as oil, have a higher
viscosity. Non-Newtonian fluids have a viscosity that is shear-rate dependent. Examples
of non-Newtonian fluids include yoghurt, paper pulp, and polymer suspensions.
PO RO US MATR IX PRO PER TIES
The default Porous material uses the Domain material (the material defined for the
domain) for the porous matrix. Select another material as needed.
Both the default Porosity p (a dimensionless number between 0 and 1) and
Permeability (SI unit: m2) use values From material as defined by the Porous material
selection. For User defined select Isotropic, Diagonal, Symmetric, or Anisotropic based on
the characteristics of the thermal conductivity, and enter another value or expression.
The components of a permeability in the case that it is a tensor (xx, yy, and so on,
representing an anisotropic permeability) are available as [Link], [Link],
and so on (using the default name br).
Forchheimer Drag
The Forchheimer Drag subnode is available from the context menu (right-click the Fluid
and Matrix Propertiesparent node) or from the Physics toolbar, Attributes menu. While
the drag of the fluid on the porous matrix in the basic Brinkman equations is
proportional to the flow velocity, (Darcys law drag), the Forchheimer drag is
proportional to the square of the fluid velocity. The latter term accounts for an inertial
179
turbulent drag effect that comes into play for fast flows through large pores. Adding
the Forchheimer term takes into account all drag contributions that the Ergun
equation covers.
FORCHHEIMER DRAG
Enter a value for the Forchheimer coefficient F (SI unit: kg/m4). The default is
0 kg/m4.
Mass Source
The Mass Source node adds a mass source (or mass sink) Qbr to the right-hand side of
the continuity equation: Equation 6-10. This term accounts for mass deposit and/or
mass creation in porous domains. The physics interface assumes that the mass
exchange occurs at zero velocity.
( ) + ( u ) = Q br
t p
(6-5)
DOMAIN SELECTION
Enter a value or expression for the Source term Qbr (SI unit: kg/(m3s)). The default
is 0 kg/(m3s).
Volume Force
Use the Volume Force node to specify the force F on the right-hand side of
Equation 6-11. It then acts on each fluid element in the specified domains. A common
application is to include gravity effects.
- u
u
---+ ( u ) ----- =
p t
p
(6-6)
Q br
1
1
T
2
- u + F
p + ----- ( u + ( u ) ) --- ( u )I + -------p
3
p2
VO L U M E F O R C E
180 |
Initial Values
The Initial Values node adds initial values for the velocity field and the pressure that can
serve as an initial condition for a transient simulation or as an initial guess for a
nonlinear solver.
IN IT IA L VA LUES
Enter initial values or expressions for the Velocity field u (SI unit: m/s) and the Pressure
p (SI unit: Pa). The default values are 0 m/s and 0 Pa, respectively.
Fluid Properties
The Fluid Properties node adds the momentum and continuity equations to solve for
free flow in non-porous domains. The node also provides an interface for defining the
material properties of the fluid.
MODEL INPUTS
Fluid properties, such as density and viscosity, can be defined through user inputs,
variables or by selecting a material. For the latter option, additional inputs, for example
temperature and/or pressure, may be required to define these properties.
Temperature
By default, the single-phase flow interfaces are set to model isothermal flow. Hence,
the Temperature is User defined and defaults to 293.15 K. If a Heat Transfer interface
is included in the component, the temperature may alternatively be selected from this
physics interface. All physics interfaces have their own tags (Name). For example, if a
Heat Transfer in Fluids interface is included in the component, the Temperature (ht)
option is available.
Absolute Pressure
This input appears when a material requires the absolute pressure as a model input.
The absolute pressure is used to evaluate material properties, but it also relates to the
value of the calculated pressure field. There are generally two ways to calculate the
pressure when describing fluid flow: either to solve for the absolute pressure or for a
pressure (often denoted gauge pressure) that relates to the absolute pressure through
a reference pressure.
The choice of pressure variable depends on the system of equations being solved. For
example, in a unidirectional incompressible flow problem, the pressure drop over the
modeled domain is probably many orders of magnitude smaller than the atmospheric
181
pressure, which, when included, may reduce the stability and convergence properties
of the solver. In other cases, such as when the pressure is part of an expression for the
gas volume or the diffusion coefficients, it may be more convenient to solve for the
absolute pressure.
The default Absolute pressure pA is p+pref where p is the dependent pressure variable
from the Navier-Stokes equations, and pref is from the user input defined at the physics
interface level. When pref is non zero, the physics interface solves for a gauge pressure.
If the pressure field instead is an absolute pressure field, pref should be set to 0.
The Absolute pressure field can be edited by clicking Make All Model Inputs Editable
(
) and entering the desired value in the input field.
FLUID PROPERTIES
182 |
Compressibility
By default the physics interface uses the Incompressible flow formulation of the
Navier-Stokes and Brinkman equations to model constant density flow. If required,
select Compressible flow (Ma<0.3) from the Compressibility list, to account for small
183
The following dependent variables (fields) are defined for this physics interfacethe
Velocity field u (SI unit: m/s) and its components, and the Pressure p (SI unit: Pa).
Domain, Boundary, Point, and Pair Nodes for the Free and Porous
Media Flow Interface
Theory for the Free and Porous Media Flow Interface
Domain, Boundary, Point, and Pair Nodes for the Free and Porous
Media Flow Interface
The Free and Porous Media Flow Interface has the following domain, boundary,
point, and pair nodes, listed in alphabetical order, available from the Physics ribbon
toolbar (Windows users), Physics context menu (Mac or Linux users), or right-click to
access the context menu (all users).
In general, to add a node, go to the Physics toolbar, no matter what
operating system you are using. Subnodes are available by clicking the
parent node and selecting it from the Attributes menu.
184 |
Fluid Properties
Forchheimer Drag
Initial Values
Volume Force
Mass Source
The following nodes (listed in alphabetical order) are described for the Laminar Flow
interface:
No Viscous Stress
Flow Continuity
Inlet
Symmetry
Outlet
Wall
Open Boundary
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.
Fluid Properties
Use the Fluid Properties node to define the fluid material, density, and dynamic
viscosity.
FLUID PROPERTIES
The default Fluid material uses the Domain material (the material defined for the
domain). Select another material as needed.
The default Density (SI unit: kg/m3) uses values From material based on the Fluid
material selection. For User defined enter another value or expression. The default is
0 kg/m3.
The Dynamic viscosity (SI unit: Pas) uses values From material based on the Fluid
material selection. For User defined enter another value or expression. The default is
0 Pas.
185
Choose domains from the Selection list, to solve for porous media flow governed by
the Brinkman equations. In the domains not selected, the Free and Porous Media Flow
interface solves for laminar flow governed by the Navier-Stokes (or Stokes) equations.
POROUS MATRIX PROPER TIES
The default Porous material uses the Domain material (the material defined for the
domain) for the porous matrix. Select another material as needed.
Porosity
The default Porosity p (a dimensionless number between 0 and 1) uses values From
material as defined by the Porous material selection. For User defined enter another
value or expression. The default is 0.
Permeability
The default Permeability br (SI unit: m2) uses values From material as defined by the
Porous material selection. For User defined select Isotropic, Diagonal, Symmetric, or
Anisotropic from the list and then enter other values or expressions. The components
of a permeability in the case that it is a tensor (xx, yy, and so on, representing an
anisotropic permeability) are available as [Link], [Link], and so on (using
the default name fp). The defaults is 0 m2.
Source Term
Enter a value or expression for an optional mass source (or sink) Source term Qbr (SI
unit: kg/(m3s)). This term accounts for mass deposit and mass creation within
domains. The physics interface assumes that the mass exchange occurs at zero velocity.
Volume Force
The Volume Force node specifies the force F on the right-hand side of the
Navier-Stokes or Brinkman equations, depending on whether the Porous Matrix
186 |
Properties node is active for the domain. Use it, for example, to incorporate the effects
of gravity in a model.
VO L U M E F O R C E
Forchheimer Drag
The Forchheimer Drag subnode is available from the context menu (right-click the
Porous Matrix Properties parent node) or from the Physics toolbar, Attributes menu.
It can be used on the domain selection that corresponds to the porous medium. For
the Brinkman equations the drag of the fluid on the porous matrix is proportional to
the flow velocity, in the same way as for Darcys law. Add a Forchheimer drag,
proportional to the square of the fluid velocity, as needed.
FORCHHEIMER DRAG
Initial Values
The Initial Values node adds initial values for the velocity field and the pressure that can
serve as an initial condition for a transient simulation or as an initial guess for a
nonlinear solver.
IN IT IA L VA LUES
Enter initial values or expressions for the Velocity field u (SI unit: m/s) and for the
Pressure p (SI unit: Pa). The default values are 0 m/s and 0 Pa, respectively.
The default Boundary condition for the wall is Slip velocity. Enter values or expressions
for the components of the Velocity of moving wall uw (SI unit: m/s).
187
188 |
u = --- p
(6-7)
In this equation, (SI unit: m2) denotes the permeability of the porous medium,
(SI unit: kg/(ms)) the dynamic viscosity of the fluid, p (SI unit: Pa) the pressure, and
u (SI unit: m/s) the Darcy velocity. The Darcys Law interface combines Darcys law
with the continuity equation:
( ) + ( u )
= Qm
t
(6-8)
In the above equation, (SI unit: kg/m3) is the density of the fluid, (dimensionless)
is the porosity, and Qm (SI unit: kg/(m3s)) is a mass source term. Porosity is defined
as the fraction of the control volume that is occupied by pores. Thus, the porosity can
vary from zero for pure solid regions to unity for domains of free flow.
If the Darcys Law interface is coupled to an energy balance, then the fluid density can
be a function of the temperature, pressure, and composition (for mixture flows). For
gas flows in porous media, the relation is given by the ideal gas law:
pM
= --------RT
(6-9)
T H E O R Y F O R T H E D A R C Y S L A W I N T E R F A C E
189
where R= 8.314 J/(molK) is the universal gas constant, M (SI unit: kg/mol) is the
molecular weight of the gas, and T (SI unit: K) is the temperature.
190 |
In porous domains, the flow variables and fluid properties are defined at any point
inside the medium by means of averaging of the actual variables and properties over a
certain volume surrounding the point. This control volume must be small compared
to the typical macroscopic dimensions of the problem, but it must be large enough to
contain many pores and solid matrix elements.
Porosity is defined as the fraction of the control volume that is occupied by pores.
Thus, the porosity can vary from zero for pure solid regions to unity for domains of
free flow.
The physical properties of the fluid, such as density and viscosity, are defined as
intrinsic volume averages that correspond to a unit volume of the pores. Defined this
way, they present the relevant physical parameters that can be measured experimentally,
and they are assumed to be continuous with the corresponding parameters in the
adjacent free flow.
191
(6-10)
- u
u
---+ ( u ) ----- =
p t
p
(6-11)
Q br
1
2
T
- u + F
p + ----- ( u + ( u ) ) --- ( u )I 1 + -------p
3
p2
In these equations:
(SI unit: kg/(ms)) is the dynamic viscosity of the fluid
u (SI unit: m/s) is the velocity vector
(SI unit: kg/m3) is the density of the fluid
p (SI unit: Pa) is the pressure
p is the porosity
(SI unit: m2) is the permeability tensor of the porous medium, and
Qbr (SI unit: kg/(m3s)) is a mass source or mass sink
Influence of gravity and other volume forces can be accounted for via the force term
F (SI unit: kg/(m2s2)).
When the Neglect inertial term (Stokes-Brinkman) check box is selected, the term
(u )(u/p) on the left-hand side of Equation 6-11 is disabled.
The mass source, Qbr, accounts for mass deposit and mass creation within the domains.
The mass exchange is assumed to occur at zero velocity.
192 |
The Forchheimer drag option, F (SI unit: kg/m4), adds a viscous force proportional
to the square of the fluid velocity, FF = F|u|u, to the right-hand side of
Equation 6-11.
In case of a flow with variable density, Equation 6-10 and Equation 6-11 must be
solved together with the equation of state that relates the density to the temperature
and pressure (for instance the ideal gas law).
For incompressible flow, the density stays constant in any fluid particle, which can be
expressed as
( ) + u = 0
t p
and the continuity equation (Equation 6-10) reduces to
u = Q br
193
194 |
195
T he S li p Flo w In t erface
The Slip Flow (slpf) interface (
), found under the Rarefied Flow branch (
) when
adding a physics interface, is used to model thermal and isothermal flows within the
slip flow regime. In the slip flow regime, the Navier-Stokes equations can be used to
model the flow of the gas, except within a thin layer of rarefied gas adjacent to the walls
(known as the Knudsen layer). The effect of the Knudsen layer on the continuum part
of the flow can be modeled by means of modified boundary conditions for the Navier
Stokes equations. Thermal effects are also important in this regime, with effects such
as thermal creep or transpiration often playing a significant role. For this reason the
Slip Flow interface includes the heat flow equations. Typically slip flow applies at
Knudsen numbers between 0.01 and 0.1.
When this physics interface is added, these default nodes are also added to the Model
BuilderFluid, External Slip Wall, and Initial Values. Then, from the Physics toolbar, add
other nodes that implement, for example, boundary conditions and volume forces.
You can also right-click Slip Flow to select physics features from the context menu.
SETTINGS
Select the Neglect inertial term (Stokes flow) check box to model flow at very low
Reynolds numbers where the inertial term in the Navier-Stokes equations can be
neglected. When this option is checked COMSOL Multiphysics solves the linear
Stokes equations for the fluid flow. The Stokes flow or creeping flow regime frequently
196 |
applies in microfluidic devices, where the flow length scales are very small. Enter a
Reference pressure level pref (SI unit: Pa). The default value is 1[atm].
If you have the Heat Transfer Module, additional check boxes are
available. See the Heat Transfer Module Users Guide for information.
DEPENDENT VA RIA BLES
197
Domain, Boundary, Edge, Point, and Pair Nodes for the Slip Flow
Interface
Theory for the Slip Flow Interface
In the COMSOL Multiphysics Reference Manual:
Pseudo Time Stepping for Laminar Flow Models
Handling Frames in Heat Transfer
Domain, Boundary, Edge, Point, and Pair Nodes for the Slip Flow
Interface
The Slip Flow Interface has these domain, boundary, edge, point, and pair nodes
available from the Physics ribbon toolbar (Windows users), Physics context menu (Mac
or Linux users), or right-click to access the context menu (all users)..
In general, to add a node, go to the Physics toolbar, no matter what
operating system you are using. Subnodes are available by clicking the
parent node and selecting it from the Attributes menu.
Periodic Condition
Slip Wall
Fluid
Symmetry
Initial Values
198 |
H E A T TR A N S F E R I N F L U I D S S U B M E N U
These nodes and one subnode are described for the Heat Transfer interface in the
COMSOL Multiphysics Reference Manual (listed in alphabetical order):
Boundary Heat Source
Heat Flux
Diffuse Surface
Heat Source
Temperature
Thermal Insulation
Thin Layer
Outflow
Translational Motion
These nodes are described for the Laminar Flow interface in this guide (listed in
alphabetical order):
Boundary Stress
Inlet
Volume Force
Open Boundary
Wall
Outlet
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.
Fluid
The Fluid node prescribes its domains to be a fluid by adding momentum and energy
transport equations.
199
The default uses the Thermal conductivity k (SI unit: W/(mK)) From material. For User
defined, select Isotropic, Diagonal, Symmetric, or Anisotropic based on the characteristics
of the thermal conductivity and enter another value or expression in the field or matrix.
The thermal conductivity describes the relationship between the heat flux
vector q and the temperature gradient T as in q = kT which is
Fouriers law of heat conduction. Enter this quantity as power per length
and temperature.
THERMODYNAMICS, FLUID
The default Dynamic viscosity (SI unit: Pas) is taken From material. Or select
Non-Newtonian power law, Non-Newtonian Carreau model, or User defined. For User
defined it uses a built-in variable for the shear rate magnitude, [Link], which makes it
possible to define arbitrary expressions of the dynamic viscosity.
Enter the components for the Velocity of moving wall uw (SI unit: m/s). the defaults
are 0 m/s
Enter a value or expression for the Wall temperature Tw (SI unit: K). the default is
293.15 K.
Select the Slip coefficientsMaxwells model (the default) or User defined.
200 |
Maxwells Model
For Maxwells Model enter values or expressions for the Tangential momentum
accommodation coefficient av (dimensionless). The tangential accommodation
coefficients are typically in the range of 0.85 to 1.0 and can be found in Ref. 4.
User defined
For User defined select a Mean free path definitionStandard (the default) or Alternate
(Shapirov). Then enter values for the following. Note that the defaults are different
based on the Mean free path definition chosen.
Thermal slip coefficient T. The defaults are 0.97 for Standard and 1.1 for Alternate
(Shapirov).
Viscous slip coefficient S. The defaults are 0.89 for Standard and 1 for Alternate
(Shapirov).
Temperature jump coefficient T. The defaults are 1.73 for Standard and 1.95 for
Alternate (Shapirov).
Theoretical and experimental values of these coefficients for various gas
surface combinations are available in Ref. 5. Note that it is more
convenient to use the Alternative (Shapirov) mean free path definition
when using this data.
Initial Values
The Initial Values node adds initial values for the velocity field, pressure, and
temperature. For the transient solver, these values define the state of the problem at
the initial time step; for the stationary solver, they serve as a starting point for the
nonlinear solver.
IN IT IA L VA LUES
Enter values or expressions for the initial estimate the solver uses for the Velocity field u
(SI unit: m/s), Pressure p (SI unit: Pa), and Temperature T (SI unit: K).
Slip Wall
Use the Slip Wall node to specify a wall with slip adjacent to a solid region in which the
heat transfer equations are solved. A wall velocity can be specified along with the slip
coefficients.
201
SLIP WALL
Enter the components for the Velocity of moving wall uw (SI unit: m/s). Select the Slip
coefficientsMaxwells model (the default) or User defined. These settings are the same
Periodic Condition
The Periodic Condition node applies periodic boundary conditions for both the fluid
flow and/or the heat transfer equations, as appropriate.
CONSTRAINT SETTINGS
Symmetry
The Symmetry node applies symmetry conditions for both the fluid flow flow and/or
the heat transfer equations, as appropriate.
CONSTRAINT SETTINGS
Continuity
The Continuity node can be added to pairs. It prescribes that the temperature field is
continuous across the pair. Continuity is only suitable for pairs where the boundaries
match.
Since it is not usual to use an assembly to represent the fluid flow domains, Continuity
is only appropriate between two domains in which Heat Transfer in Solids applies.
CONSTRAINT SETTINGS
202 |
4 T g
3
(7-1)
where uslip is the slip velocity, n is the boundary normal, is the viscous stress tensor,
is the viscosity of the gas, is its density, and Tg is its temperature. The factor G
(which has dimensions of length) is given by:
2 av
G = ---------------
av
where is the mean free path and av is the tangential momentum accommodation
coefficient (for a model in which the surface reflects some molecules diffusely and
203
some secularly, this is equivalent to the fraction of molecules which are reflected
diffusely).
A slightly different definition of the mean free path is used here: in terms
of Maxwells original parameter, l, the mean free path is = 2l/3. In
Maxwells original paper, he uses a coordinate based notation, the vector
notation is available in Ref. 2.
In Equation 7-1 the left-hand term represents the phenomena of viscous slip, whilst
the right-hand term produces thermal slip or transpiration. Maxwell was aware that the
temperature of the gas was not necessarily equal to that of the wall but formulated the
boundary condition using only the gas temperatures. When the wall temperature is
included, the following equations are obtained (Ref. 2 and Ref. 3):
T g
T w = T g T n T g
where Tw is the wall temperature, s is the viscous slip coefficient, T is the thermal
slip coefficient, and T is the temperature jump coefficient. Within a generalized
Maxwells model the three coefficients s, T, and T are given by:
2 av
s = --------------av
3
T = --4
2 a v 2
(7-2)
where is the thermal conductivity of the gas. The mean free path can be computed
from the gas properties using the following equation (Ref. 3):
2
= ---------- c
8p
8RT 1 2
= -------
c = ------------
M n
consequently:
204 |
12
(7-3)
----------
2p
In some cases an alternative definition of the mean free path, , is used (as, for
example, in Ref. 5):
' = -------
2
It is possible to use this definition of the mean free path when entering user defined
values for the slip coefficients.
Values of the parameters given in Equation 7-2, are available in Ref. 5,
although care should be taken with the slip and temperature jump
coefficients to adjust for differences between the authors definition of
the mean free path and that used in COMSOL Multiphysics.
Also note that the formulation assumes that the ideal gas laws apply. This
assumption is implicit in the derivation of the above equations (Ref. 3)
and at the level of approximation of the equations is reasonable.
205
206 |
Species interface that is available with the basic COMSOL Multiphysics license.
In addition, the Transport of Diluted Species in Porous Media interface is available.
Both physics interfaces are found under the Chemical Species Transport
branch (
).
In this chapter:
The Transport of Diluted Species Interface
The Transport of Diluted Species in Porous Media Interface
Theory for the Transport of Diluted Species Interface
207
T he T r a ns po r t of D i l u t ed Sp eci es
Interface
The Transport of Diluted Species (tds) interface (
), found under the Chemical Species
), is used to calculate the concentration field of a dilute solute in
a solvent. Transport and reactions of the species dissolved in a gas, liquid or solid can
be handled with this interface. The driving forces for transport can be diffusion by
Ficks law, convection when coupled to a flow field, and migration, when coupled to
an electric field.
Transport branch (
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.
Some features are only available in a limited set of add-on products. For
a detailed overview of which features are available in each product, visit
[Link]
SETTINGS
If any parts of the model geometry should not partake in the mass transfer model,
remove that part from the selection list.
TR A N S P O R T M E C H A N I S M S
Diffusion is always included. By default, the Convection check box is selected under
Additional transport mechanisms.
208 |
C H A P T E R 8 : C H E M I C A L S P E C I E S TR A N S P O R T I N T E R F A C E S
Note: Not all additional transport mechanisms listed below are available in all
products. For detail see [Link]
Select the Migration in electric field check box to activate the migration transport of
ionic species. See further the theory section Adding Transport Through Migration.
Select the Adsorption in porous media check box to activate the adsorption of solutes
in porous media. See further Adsorption.
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.
CONSISTENT STABILIZATION
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 Galeothe 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 options is less diffusive compared to the Do Carmo and Galeo
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)/[Link], where [Link]
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
T H E TR A N S P O R T O F D I L U T E D S P E C I E S I N T E R F A C E
209
ADVANCED SETTINGS
The Compute boundary fluxes check box is activated by default so that COMSOL
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, COMSOL 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.
The flux variables affected in the interface are
ndflux_c (where c is the dependent variable for the concentration) 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) is the normal
total flux and corresponds to the boundary flux plus additional transport terms, for
example, the convective flux when you use the non-conservative 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.
210 |
C H A P T E R 8 : C H E M I C A L S P E C I E S TR A N S P O R T I N T E R F A C E S
The dependent variable name is 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 (
concentration (
) buttons as needed.
) and Remove
FURTHER READING
T H E TR A N S P O R T O F D I L U T E D S P E C I E S I N T E R F A C E
211
This interface includes free and porous media flow with immobile and mobile phases,
including diffusion, convection, dispersion, adsorption, and volatilization in porous
media. It supports cases where either the solid phase substrate is exclusively immobile,
or when a gas-filling medium is also assumed to be immobile.
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,
provide an interface for defining the properties of the porous media, as well as
additional properties governing adsorption, volatilization, dispersion and diffusion,
and the velocity field to model convection.
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
BuilderPorous 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
212 |
C H A P T E R 8 : C H E M I C A L S P E C I E S TR A N S P O R T I N T E R F A C E S
Web link:
[Link]
-transport-sorbing-solute-490
Concentration
Periodic Condition
Flux
Flux Discontinuity
Inflow
Reactions
Initial Values
Species Source
Symmetry
Mass-Based Concentrations
No Flux
Outflow
Transport Properties
Volatilization
T H E TR A N S P O R T O F D I L U T E D S P E C I E S I N T E R F A C E
213
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.
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 sub-node is available
from the context menu as well as from the Physics toolbar, Attributes menu. Note that
this feature is only available in some COMSOL products. See details:
[Link]
MODEL INPUTS
214 |
C H A P T E R 8 : C H E M I C A L S P E C I E S TR A N S P O R T I N T E R F A C E S
When the Migration in electric field check box is selected on the Settings window for
Transport of Diluted Species, select the source of the electric potential field and,
optionally, temperature.
Enter values or expressions for the Electric potential V, which is User defined; this
input option is always available.
Select the electromagnetic field solved by an AC/DC-based interface that has also
been added to the model. This works like the velocity field model input described
above.
The default mobility model is the Nernst-Einstein relation. This also requires a
temperature input, which will be exposed in the same manner as the electric
potential described above.
Note that the migration in electric fields feature is only available in some COMSOL
products. See details: [Link]
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.
Enter the Diffusion coefficient Dc for each species. This can be a scalar value for isotropic
diffusion or a tensor describing anisotropic diffusion. Select the appropriate tensor
type Isotropic, Diagonal, Symmetric, or Anisotropic that describes the diffusion
transport, and then enter the values in the corresponding element (one value for each
species).
Note that multiple species, as well as Migration in Electric fields (described below) is
only available for certain COMSOL add-on products. See details:
[Link]
MIGRATION IN ELECTRIC FIELD
This section is available when the Migration in electric field check box is selected. 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 typeIsotropic, Diagonal, Symmetric, or
Anisotropicand 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.
T H E TR A N S P O R T O F D I L U T E D S P E C I E S I N T E R F A C E
215
Specify the temperature (if you are using mobilities based on the Nernst-Einstein
relation) and electric field in the Model Inputs section.
EXAMPLE MODELS
Web link:
[Link]
Web link:
[Link]
Turbulent Mixing
Note that the Turbulent Mixing node is only available for some product add-ons. See
details: [Link]
This subnode is available from the context menu (right-click the Transport Properties
parent node) as well as from the Physics toolbar, Attributes menu, if Convection is
selected as a transport mechanism. Use this node to account for the turbulent mixing
caused by the eddy diffusivity. An example is when the specified velocity field
corresponds to a RANS solution.
TU R B U L E N T M I X I N G P A R A M E T E R S
Some physics interfaces provide the turbulent kinematic viscosity, and these appear as
options in the Turbulent kinematic viscosity T list. The list always contains the User
defined option where any value or expression can be entered.
The default Turbulent Schmidt number ScT is 0.71 (dimensionless).
216 |
C H A P T E R 8 : C H E M I C A L S P E C I E S TR A N S P O R T I N T E R F A C E S
FURTHER READING
See the section About Turbulent Mixing (this link is available if you have the CFD
Module documentation installed)
Turbulent Mixing of a Trace Species: Application Library path
CFD_Module/Single-Phase_Tutorials/turbulent_mixing
Web link:
[Link]
2727
Initial Values
This node specifies the initial values for the concentration of each species. These serve
as an initial guess for a stationary solver or as initial conditions for a transient
simulation.
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.
IN IT IA L VA LUES
Enter a value or expression for the initial value of the Concentration or concentrations,
ci. This also serves as a start 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 the 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.
T H E TR A N S P O R T O F D I L U T E D S P E C I E S I N T E R F A C E
217
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/(m3s)), 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 Engineering Module, the reaction rate expressions can be
automatically generated and picked up using the drop-down menu. See example in the
application Fine Chemical Production in a Plate Reactor as linked below.
R E A C T I N G VO L U M E
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 non-porous domains, the settings of the Reacting Volume section has no
impact.
Note that the Reacting Volume section is only available in the products that provide
the Transport of Diluted Species in Porous Media interface. See details:
[Link]
For Total volume the reaction expressions in mol/(m3s) are specified per unit volume
of the model domain (multiplied by unity).
For Pore volume the reaction expressions in mol/(m3s) are specified per unit volume
of total pore space. The reaction expressions will be multiplied by the domain porosity,
p. (p equals unity for non-porous domains).
For Liquid phase the reaction expressions in mol/(m3s) 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.
218 |
C H A P T E R 8 : C H E M I C A L S P E C I E S TR A N S P O R T I N T E R F A C E S
FURTHER READING
See the theory chapter on chemical species transport, starting with the section Mass
Balance Equation.
Fine Chemical Production in a Plate Reactor: Application Library
path Chemical_Reaction_Engineering_Module/
Reactors_with_Mass_and_Heat_Transfer/plate_reactor
Web link:
[Link]
eactor-8589
No Flux
This node is the default boundary condition on exterior boundaries. It represents
boundaries where no mass flows in or out of the boundaries. Hence, the total flux is
zero.
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.
B O U N D A R Y C O N D I T I O N TY P E
This section in the settings is only available for some products. Search for Inflow on
the page: [Link] for more details on
availability.
The option Concentration constraint constrains the concentration values on the
boundary by the use of point-wise constraints. The other option, Flux (Danckwerts) can
be more stable and fast to solve when high reaction rates are anticipated in the vicinity
T H E TR A N S P O R T O F D I L U T E D S P E C I E S I N T E R F A C E
219
of the inlet. Oscillations on the solutions can also be avoided in such cases. The latter
condition uses a flux boundary condition based on the velocity across the boundary
and the concentration values. See further details in the theory section
CONSTRAINT SETTINGS
Outflow
This node is not available if Diffusion only is included in the model.
Set this condition at outlets where species are transported out of the model domain by
fluid motion. 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 c ) = 0
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 (SI unit: mol/m3) 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
220 |
C H A P T E R 8 : C H E M I C A L S P E C I E S TR A N S P O R T I N T E R F A C E S
Flux
This node can be used to specify the total species flux across a boundary. The total flux
of species c is defined accordingly:
n ( cu Dc zu m Fc ) = N 0
where N0 is an arbitrary user-specified flux expression (SI unit: mol/(m2s)). For
example, N0 can represent a flux from or into a much larger surrounding environment,
a phase change, or a flux due to chemical reactions. N0 can also be a function of the
concentration and the electric potential (if the mass transport includes migration of
ionic species).
When diffusion is the only transport mechanism present, the flux condition is extended
to include a mass transfer term to describe flux into a surrounding environment:
n ( Dc ) = N 0 + k c ( c b c )
where kc is a mass transfer coefficient (SI unit: m/s), and cb is the concentration
(SI unit: mol/m3) in the surroundings of the modeled system (the bulk
concentration). The mass transfer coefficient (to be specified) is often given by
boundary-layer theory.
INWARD FLUX
This is used to individually specify the flux of each species. To use another boundary
condition for a specific species, click to clear the check box for the mass fraction of that
species.
Note: Use a minus sign when specifying a flux leaving the system.
Symmetry
The Symmetry node can be used to represent boundaries where the species
concentration is symmetric, that is, where there is no mass flux in the normal direction
across the boundary.
This boundary condition is identical to that of the No Flux node, but applies to all
species and cannot be applied to individual species.
Flux Discontinuity
This node represents a discontinuity in the mass flux across an interior boundary:
T H E TR A N S P O R T O F D I L U T E D S P E C I E S I N T E R F A C E
221
n ( Nd Nu ) = N0
N = ( cu Dc zu m Fc )
where the value N0 (SI unit: mol/(m2s)) specifies the jump in flux at the boundary.
This can be used to model a boundary source, for example a surface reaction,
adsorption or desorption.
FLUX DISCONTINUITY
In this section the jump in species flux (or surface source) is specified.
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.
Periodic Condition
The Periodic Condition node can be used to define periodicity or antiperiodicity
between two boundaries. The node can be activated on more than two boundaries, in
which case the feature tries to identify two separate surfaces that can each consist of
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
222 |
C H A P T E R 8 : C H E M I C A L S P E C I E S TR A N S P O R T I N T E R F A C E S
For the Reacting Flow in Porous Media, Diluted Species interface, which is available
in some add-on products, the Line Mass Source node is available in two versions, one
for the fluid flow (Fluid Line Source) and one for the species (Species Line Source).
SELECTION
The Line Mass Source feature is available for all dimensions, but the applicable selection
differs between the dimensions.
MODEL DIMENSION
2D
Points
2D Axisymmetry
3D
Edges
SPECIES SOURCE
Enter the source strength, q l,c , for each species (SI unit: mol/(ms)). 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
T H E TR A N S P O R T O F D I L U T E D S P E C I E S I N T E R F A C E
223
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.
Point sources located on a boundary or on an edge affect the adjacent computational
domains. This has the effect, for example, that the physical strength of a point source
located in a symmetry plane is twice the given strength.
FURTHER READING
Open Boundary
This feature is only available in a limited set of add-on products. See
[Link] for more details on availability.
Use this node to set up mass transport across boundaries where both convective inflow
and outflow can occur. Use this boundary condition to specify an exterior species
concentration 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 provided as a model input to the Transport of Diluted Species
interface.
EXTERIOR CONCENTRATION
Enter a Layer thickness ds (SI unit: m). The default is 0.005 m (5 mm). Enter a
Diffusion coefficient Ds,c (SI unit: m2/s). The default is 0.
224 |
C H A P T E R 8 : C H E M I C A L S P E C I E S TR A N S P O R T I N T E R F A C E S
Equilibrium Reaction
This feature is only available in a limited set of add-on products. See
[Link] for more details on availability.
Use this node to model an equilibrium reaction in a domain. This feature is available
with two species or more.
The equilibrium reaction is defined by the relation between the chemical activities of
the chemical species participating in the reaction (the equilibrium condition), and the
stoichiometry of the reaction.
The node solves for an additional degree of freedom (the reaction rate) 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). The default
is 1. Enter a value or expression for the Unit activity concentration Ca0 (SI
unit: mol/m3). The default is 1 x 10-3 mol/m3. Equilibrium constant creates an
equilibrium condition based on the stoichiometric coefficients, the species activities,
and the law of mass action.
For User defined enter an Equilibrium expression Eeq (dimensionless).
T H E TR A N S P O R T O F D I L U T E D S P E C I E S I N T E R F A C E
225
STOICHIOMETRIC COEFFICIENTS
Enter a value for the stoichiometric coefficientc (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.
226 |
C H A P T E R 8 : C H E M I C A L S P E C I E S TR A N S P O R T I N T E R F A C E S
Reaction Coefficients
Add this node to the Electrode-Electrolyte Interface Coupling and Porous Electrode
Coupling features to define molar fluxes and sources based on electrode current
densities in an Electrochemistry interface.
The molar flux or source is proportional to the stoichiometric coefficients and the
current density according to Faradays law.
All current densities from the Electrode Reaction (iloc, SI unit: A/m2) or the Porous
Electrode Reaction (iv, SI unit: A/m3) nodes are available for selection as the Coupled
reaction, and user-defined expressions are also supported.
Enter the Number of participating electrons nm (dimensionless) and the Stoichiometric
coefficient vc (dimensionless) as explained in the Electrode Reaction documentation or
the theory section.
Use multiple subnodes to couple to multiple reactions.
T H E TR A N S P O R T O F D I L U T E D S P E C I E S I N T E R F A C E
227
MATRIX PROPERTIES
Select an option from the Porous material list. The default is Domain material.
By default the Porosity, p (dimensionless) is taken From material. For User defined
enter a different value. The default is 0.3.
When the Adsorption in porous media check box is selected on the Settings window for
the physics interface, the default Density (SI unit: kg/m3) is taken From material. For
User defined enter a different value. The default is 1400 kg/m3.
DIFFUSION
This section is available when the Adsorption in porous media check box is selected on
the Settings window for the physics interface.
Select a Sorption typeLangmuir (the default), Freundlich, or User defined.
For Langmuir enter a Langmuir constant kL,c (SI unit: m3/mol) and an Adsorption
maximum cp,max,c (SI unit: mol/kg). The defaults are 0 m3/mol and 0 mol/kg,
respectively.
For Freundlich enter a Freundlich constant kF,c (dimensionless), a Freundlich exponent
NF,c (dimensionless), and Reference concentration cref,c (SI unit: m3/mol).
For User defined enter an Adsorption isotherm kP,c (SI unit: m3/kg).
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.
228 |
C H A P T E R 8 : C H E M I C A L S P E C I E S TR A N S P O R T I N T E R F A C E S
Select an option from the Dispersion tensor listUser 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 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 listIsotropic
(the default) or Transverse isotropic based on the properties of the porous media. For
isotropic porous media specify the longitudinal and transverse dispersivities. For
transverse isotropic porous media specify the longitudinal, horizontal transverse, and
vertical transverse dispersivities.
FURTHER READING
See the theory chapter in the section Transport of Diluted Species in Porous Media.
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).
T H E TR A N S P O R T O F D I L U T E D S P E C I E S I N T E R F A C E
229
DIFFUSION
This section is available when the Adsorption 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.
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.
VO L A T I L I Z A T I O N
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 only available in a limited set of add-on products. See
[Link] for more details on availability.
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 the boundary condition to model a thin layer through which mass is transported
by volatilization only. To set up the node, specify the layer thickness and the
atmospheric concentration of each species in the thin layer for each transported
species.
230 |
C H A P T E R 8 : C H E M I C A L S P E C I E S TR A N S P O R T I N T E R F A C E S
VO L A T I L I Z A T I O N
Enter a Layer thickness ds and the atmospheric concentration for each species. The Gas
diffusion coefficient DG,c (SI unit: m2/s) and the Volatilization coefficient kG,c
(dimensionless) for each species are taken from the adjacent Partially Saturated Porous
Media domain.
Here you can specify the bed porosity which is the void fraction between the pellets in
the packed bed structure. Either you can specify the densities of the bed and one pellet,
Density option, or specify the porosity directly Porosity option.
PELL ET SIZE D IST RIBUTIO N
Different pellet sizes can be modeled in the same bed. Select a Pellet size distribution
Uniform size (the default), Two sizes, Three sizes, Four sizes, or Five sizes to select up to
five different particle sizes. Then based on this choice enter the following as applicable:
Radius rpe (SI unit: m). The default is 1 x 10-3 m.
Radii and volume percentages: If several sizes are selected, the radius of each size, and
its volume percentage of the total pellet volume must be added.
Note that different chemical reactions can be specified for each size of the pellet.
PELL ET PARA MET ERS
T H E TR A N S P O R T O F D I L U T E D S P E C I E S I N T E R F A C E
231
Enter also the Diffusion coefficient Dpe,c (SI unit: m2/s). In the User defined case an
Effective diffusion coefficient Dpeff,c (SI unit: m2/s) is entered. The default is 1 x
10-9 m2/s in both cases.
P E L L E T- F L U I D S U R F A C E
The pellet surface to fluid coupling has two Coupling type options available:
Continuous concentration, cpei=ci, assuming no resistance to the mass transport
between pellet and fluid. This options assumes that all the resistance against mass
transfer is inside the pellet pores.
Film resistance (mass flux) when the above assumption does not hold and the mass
transfer of also limited by the convective and diffusional properties of the fluid
surrounding the pellet. This option employs a film resistance theory as described in
the section Theory for the Reactive Pellet Bed.
The Film resistance (mass flux) option computes the inward surface flux,
Ni,inward=hDi(ci-cpei). hDi is the mass transfer coefficient (SI unit: m/s) and can be
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 in most
cases apply for spherical dumped pellets or otherwise use the User defined option.
The Sherwood number expression can be computed from three available expressions:
Frssling, Rosner, and Garner and Keey. The Frssling equation is the default and
probably the most commonly used for packed spheres. All 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
DistributionCubic root sequence (the default), Linear, or Square root sequence. Enter
the Number of elements Nelem.
CONSTRAINT SETTINGS
232 |
C H A P T E R 8 : C H E M I C A L S P E C I E S TR A N S P O R T I N T E R F A C E S
FURTHER READING
Theory for the Reactive Pellet Bed in the Theory section of this manual.
For an application using the Reactive Pellet Bed feature, see
A Multiscale 3D Packed Bed Reactor: Application Library path
Chemical_Reaction_Engineering_Module/Reactors_with_Porous_Catalysts/
packed_bed_reactor_3d
Web link:
[Link]
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 S (SI unit: mol/(m3s)) for each of the species solved for. Enter a
value or expression in the field of the corresponding species.
Hygroscopic Swelling
The Hygroscopic Swelling multiphysics coupling node (
) is used for moisture
concentration coupling between the Solid Mechanics interface and either the
Transport of Diluted Species or Transport of Diluted Species in Porous Media
interfaces.
Hygroscopic swelling is an effect of internal strain caused by changes in moisture
content. This strain can be written as
hs = h ( c mo c mo,ref )
T H E TR A N S P O R T O F D I L U T E D S P E C I E S I N T E R F A C E
233
More information about how to use hygroscopic swelling can be found in Hygroscopic
Swelling Coupling section in the Structural Mechanics module Users Guide.
More information about multiphysics coupling nodes can be found in the section The
Multiphysics Node in the COMSOL Multiphysics Reference Manual.
234 |
C H A P T E R 8 : C H E M I C A L S P E C I E S TR A N S P O R T I N T E R F A C E S
T he o r y f o r t he Tran sp ort of D i l u t ed
S pe c i e s Inte r f a c e
The Transport of Diluted Species Interface provides a predefined modeling
environment for studying the evolution of chemical species transported by diffusion
and convection as well as migration due to an electric field. The physics interface
assumes that all species present are dilute; that is, that their concentration is small
compared to a solvent fluid or solid. As a rule of thumb, a mixture containing several
species can be considered dilute when the concentration of the solvent is more than
90 mol%. Due to the dilution, mixture properties such as density and viscosity can be
assumed to correspond to those of the solvent.
When studying mixtures that are not dilute, the mixture and transport properties
depend of the composition, and a different physics interface is recommended. See The
Transport of Concentrated Species Interface for more information.
Ficks 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:
Mass Balance Equation
Convective Term Formulation
References
Crosswind Diffusion
The Transport of Diluted Species in Porous Media Interface theory is also included.
T H E O R Y F O R T H E TR A N S P O R T O F D I L U T E D S P E C I E S I N T E R F A C E
235
(8-1)
Equation 8-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.
ci is the concentration of the species (SI unit: mol/m3)
Di denotes the diffusion coefficient (SI unit: m2/s)
Ri is a reaction rate expression for the species (SI unit: mol/(m3s))
u is the velocity vector (SI unit: m/s)
The flux vector N (SI unit: mol/(m2s)) is associated with the mass balance equation
above and used in boundary conditions and flux computations. For the case where the
diffusion and convection are the only transport mechanisms, the flux vector is defined
as
N i = D c + u c
(8-2)
If Migration in Electric Fields is activated, the flux vector is amended with the
migration term as show in the section Adding Transport Through Migration.
The first term on the left-hand side of Equation 8-1 corresponds to the accumulation
(or indeed consumption) of the species.
The second term accounts for the diffusive transport, accounting for the interaction
between the dilute species and the solvent. A user input field for the diffusion
coefficient is available. Anisotropic diffusion coefficient tensor input is supported.
The third term on the left hand side of Equation 8-1 describes the convective transport
due to a velocity field u. This field can be expressed analytically or be obtained from
coupling this physics interface to one that computes fluid flow, such as Laminar Flow.
On the right-hand side of the mass balance equation (Equation 8-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
236 |
C H A P T E R 8 : C H E M I C A L S P E C I E S TR A N S P O R T I N T E R F A C E S
interfacethe Reaction node, which has a field for specifying a reaction equation using
the variable names of all participating species.
ai i
i products
K eq = ---------------------------------
ai i
i reactants
ai
i
T H E O R Y F O R T H E TR A N S P O R T O F D I L U T E D S P E C I E S I N T E R F A C E
237
The Equilibrium Reaction node solves for a reaction rate so that the equilibrium
condition is always fulfilled in the domain.
c,i is set to unity when the Equilibrium constant is selected on the
Settings window. For non-unity activity coefficients, a user defined
equilibrium condition can be used.
EQUILIBRIUM REACTIONS AND INFLOW BOUNDAR Y CONDITIONS
238 |
C H A P T E R 8 : C H E M I C A L S P E C I E S TR A N S P O R T I N T E R F A C E S
solving for a stationary problem with all non-equilibrium reaction rates set to zero.
Manual scaling of the reaction rate dependent variables is needed in this study step.
(8-3)
c
conservative: ----- + ( cu ) = ( D c ) + R
t
(8-4)
and each is treated slightly differently by the solver algorithms. In these equations D
(SI unit: m2/s) is the diffusion coefficient, R (SI unit: mol/(m3s)) 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, cu, would equal zero for an incompressible fluid and
would result in the non-conservative formulation above. This is in fact the default
formulation in this physics interface and ensures that nonphysical source terms do not
emerge from a solution for the flow field. To switch between the two formulations,
) and select Advanced Physics Options.
click the Show button (
T H E O R Y F O R T H E TR A N S P O R T O F D I L U T E D S P E C I E S I N T E R F A C E
239
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
[Link]
POINT SOURCE
V 0
Qc =
q p,c
(8-5)
is added at a point in the geometry. As can be seen from Equation 8-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.
The finite element representation of Equation 8-5 corresponds to a finite
concentration at a point with the effect of the point source spread out over a region
around the point. The size of the region depends on the mesh and on the strength of
240 |
C H A P T E R 8 : C H E M I C A L S P E C I E S TR A N S P O R T I N T E R F A C E S
the source. A finer mesh gives a smaller affected region, but also a more extreme
concentration value. It is important not to mesh too finely around a point source since
this can result in unphysical concentration values. It can also have a negative effect on
the condition number for the equation system.
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/(ms)), this can be expressed as
lim
S 0
Ql,c =
q l,c
(8-6)
For feature node information, see Line Mass Source and Point Mass
Source.
T H E O R Y F O R T H E TR A N S P O R T O F D I L U T E D S P E C I E S I N T E R F A C E
241
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
[Link]
(8-7)
where
ci (SI unit: mol/ m3) denotes the concentration of species i
Di (SI unit: m2/s) is the diffusion coefficient of species i
u (SI unit: m/s) is the fluid velocity
F (SI unit: As/mol) refers to Faradays constant
V (SI unit: V) denotes the electric potential
zi (dimensionless) is the charge number of the ionic species, and
um,i (SI unit: mols/kg) its ionic mobility
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 ElectrostaticsSecondary Current Distribution interface. Sometimes it is assumed to
be a supporting electrolyte, which simplifies the transport equations.
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/(molK)) is the molar gas constant and T (SI unit: K) the
temperature.
242 |
C H A P T E R 8 : C H E M I C A L S P E C I E S TR A N S P O R T I N T E R F A C E S
Note: With regards to migration, the assumption in the Transport of Diluted Species
Interface is that the concentrations are so low that each ionic species does not
contribute to a net charge in the solution. If this assumption does not hold, the
physics interface called Nernst-Planck Equations needs to be used. This latter
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.
Modeling and solving for a supporting electrolyte in the Electrostatics or Secondary
Current Distribution interfaces will give a potential distribution that drives the
migration in the Transport of Diluted Species Interface.
The current density vector is proportional to the sum of all species fluxes as expressed
by Faradays 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:
T H E O R Y F O R T H E TR A N S P O R T O F D I L U T E D S P E C I E S I N T E R F A C E
243
i = F
z j u m , j F c j
2
(8-8)
Equation 8-8 is simply Ohms law for ionic current transport and can be simplified to
i =
(8-9)
where is the conductivity of the supporting electrolyte. A current balance gives the
current and potential density in the cell
i = 0
which in combination with Equation 8-9 yields:
( ) = 0
(8-10)
Equation 8-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 very high cell
Pclt numberthat 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
provides two crosswind diffusion options using different formulations. Observe that
crosswind diffusion makes the equation system nonlinear even if the transport
equation is linear.
DO CARMO AND GALEO
244 |
C H A P T E R 8 : C H E M I C A L S P E C I E S TR A N S P O R T I N T E R F A C E S
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 Galeo
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
under-resolved gradients are less smeared out.
The following equations for the concentrations, ci, describe the transport of solutes in
a variably saturated porous medium for the most general case, when the pore space is
primarily filled with liquid but also contain pockets or immobile gas:
( c ) + ( c ) + (a c ) + u c =
i
i
t
t b P, i t v G, i
[ ( D D, i + D e, i ) c i ] + R i + S i
(8-11)
On the left-hand side of Equation 8-11, 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).
T H E O R Y F O R T H E TR A N S P O R T O F D I L U T E D S P E C I E S I N T E R F A C E
245
In Equation 8-11 ci denotes the concentration of species i in the liquid (SI unit:
mol/m3), cP, i the amount adsorbed to (or desorbed from) 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 bulk (or drained) density, b = (1 p),
and the solid phase density (SI unit: kg/m3).
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 = sp.
The resulting gas volume fraction is av = p = (1-s)p.
On the right-hand side of Equation 8-11, the first term introduces the spreading of
species due to mechanical mixing () as well as from diffusion and volatilization to the
gas phase. The 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 8-11 describe production or
consumption of the species; Ri is a reaction rate expression which can account for
reactions in the liquid, solid, or gas phase, and Si is an arbitrary source term, for
example due to a fluid flow source or sink.
In order to solve for the solute concentration of species i, ci, the solute mass sorbed to
solids cP,i and dissolved in the gas-phase cG,i are assumed to be functions of ci.
Expanding the time-dependent terms gives
c
( i ) + ( b c P, i ) + (a v c G, i) =
t
t
t
(8-12)
c i
p
( + b k P, i + a v k G, i )
+ ( 1 k G, i )c i ( P c P, i k G, i c i )
t
t
t
where kP,i = cP,i/ci is the adsorption isotherm and kG,i = cG,i/ci is the linear
volatilization. Equation 8-11 can then be written as
c i
p
+ ( 1 k G, i )c i ( P c P, i k G, i c i )
+ u c i
t
t
t
(8-13)
= [ ( D D + D e ) c i ] + R i + S i
( + b k P, i + a v k G, i )
246 |
C H A P T E R 8 : C H E M I C A L S P E C I E S TR A N S P O R T I N T E R F A C E S
c i
p
+ ( c i P c P, i )
+ u c i =
t
t
[ ( D D, i + F, i D F, i ) c i ] + R i + S i
( p + b k P, i )
(8-14)
Convection
Convection describes the movement of a species, such as a pollutant, with the bulk
fluid velocity. The velocity field u corresponds to a superficial volume average over a
unit volume of the porous medium, including both pores and matrix. This velocity is
sometimes called Darcy velocity, and defined as volume flow rates per unit cross
section of the medium. This definition makes the velocity field continuous across the
boundaries between porous regions and regions with free flow.
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
Darcys 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 = ----p
Saturated
u
u a = ---
Partially saturated
where p is the porosity and = sp the liquid volume fraction, and s the saturation, a
dimensionless number between 0 and 1.
T H E O R Y F O R T H E TR A N S P O R T O F D I L U T E D S P E C I E S I N T E R F A C E
247
Figure 8-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.
c
( i ) + ( b c P, i ) + (a v c G, i) + uc i =
t
t
t
[ ( D D, i + D e, i ) c i ] + R i + S i
(8-15)
If the conservative formulation is expanded using the chain rule, then one of the terms
from the convection part, ciu, would equal zero for an incompressible fluid and
would result in the non-conservative formulation described in Equation 8-11.
When using the non-conservative formulation, which is the default, the fluid is
assumed incompressible and divergence free: u = 0. The non-conservative
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 Non-conservative form (the default) or Conservative form. The
conservative formulation should be used for compressible flow.
248 |
C H A P T E R 8 : C H E M I C A L S P E C I E S TR A N S P O R T I N T E R F A C E S
Diffusion
The effective diffusion in porous media, De, depends on the structure of the porous
material and the phases involved. Depending on the transport of diluted species occurs
in free flow, saturated or partially saturated porous media, the effective diffusivity is
defined as:
De = DL
Free Flow
p
D e = ----- D L
L
D e = ----- D L
L
av
7 3 2
7 3 2
F =
5 2 2
5 2 2
, G = av
For saturated porous media = p. The fluid tortuosity for the Millington and Quirk
model is
F = p
1 3
1 2
User defined expressions for the tortuosity factor can also be applied.
T H E O R Y F O R T H E TR A N S P O R T O F D I L U T E D S P E C I E S I N T E R F A C E
249
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.
D Dii
ui
uj
= 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.
250 |
C H A P T E R 8 : C H E M I C A L S P E C I E S TR A N S P O R T I N T E R F A C E S
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:
2
u
v
w
D Lxx = 1 ------- + 2 ------- + 3 ------u
u
u
v
u
w
D Lyy = 1 ------- + 2 ------- + 3 ------u
u
u
w
u
v
D Lzz = 1 ------- + 3 ------- + 3 ------u
u
u
uv
D Lxy = D Lyx = ( 1 2 ) ------u
uw
D Lxz = D Lzx = ( 1 3 ) -------u
vw
D Lyz = D Lzy = ( 1 3 ) -------u
(8-16)
In Equation 8-16 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 physics
interface predefines three relationships to predict the solid concentrations, cPi from the
concentration in the liquid phase, ci:
T H E O R Y F O R T H E TR A N S P O R T O F D I L U T E D S P E C I E S I N T E R F A C E
251
cP = KP c
cP = KF c
K L c Pmax c
c P = ------------------------1 + KL c
c P
-------- = ( K P c )
c
c
c P
-------- = NK F c N 1
c
K L c Pmax
c P
-------- = --------------------------2
c
( 1 + KL c )
User defined
Freundlich (Ref. 3)
(8-17)
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
252 |
C H A P T E R 8 : C H E M I C A L S P E C I E S TR A N S P O R T I N T E R F A C E S
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.
Outflow
Macroscale: c
Concentration in
fluid passing through
bed
cpe
Microscale:
Concentration
in porous pellet
Inflow
Figure 8-3: Schematic showing the macroscale (bed volume) and the microscale (pellet)
T H E O R Y F O R T H E TR A N S P O R T O F D I L U T E D S P E C I E S I N T E R F A C E
253
The transport and reaction equations inside the pellets is done with an extra dimension
feature 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 non-dimensional radial coordinate on the domain 0-1. Different pellet radii and
even uneven radius distributions can be used.
The model equations assume spherical particles of a radius rpe. Consider the microscale
concentration cpe inside an individual porous pellet or particle, and the
macro-concentration c in the packed bed gas volume.
The pellet radius input can be:
one uniform pellet radius, that 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
t
(8-18)
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
hand side is per unit volume of bed, (SI unit: mol/(m3 s)).
Looking inside a pellet: Assuming no concentration variations in the space-angle (, )
direction, but only in the radial (r) direction of the spherical particle allows a
spherically symmetric reaction-diffusion transport equation inside the pellet. If rdim
(SI unit: m) is the spatial radial coordinate in the pellet, and rpe is the pellet radius, the
non-dimensional coordinate r=rdim/rpe can defined. The modeling domain on r goes
from 0 to 1.
254 |
C H A P T E R 8 : C H E M I C A L S P E C I E S TR A N S P O R T I N T E R F A C E S
0 (center)
Figure 8-4: Modeling domain in a pellet for a dimensional coordinate (top) and
non-dimensional coordinate (bottom.)
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:
2 2
2
2
4N r r pe pe (c pe, i) + ( r D pe, i c pe, i ) = r r pe R pe, i
(8-19)
where
N is the number of pellets per unit volume of bed.
Equally as in Equation 8-18, cpe is the interstitial (physical) species concentration in
moles/m3 fluid volume element inside the pore channel,
Rpe is the reaction rate in moles/(m3 s) of particle volume. It should be stressed
that the user input of R is per unit volume of pellet.
The effective diffusion coefficient in Equation 8-18 and Equation 8-19 depend on the
porosity pe, tortuosity and the physical gas diffusivity D of in the porous particle
generally as
pe D
D pe = ------------ .
The available tortuosity models for porous media are the Millington and Quirk (Ref.
12),
= pe
1 3
D pe = pe
43
D,
(8-20)
Bruggemann,
T H E O R Y F O R T H E TR A N S P O R T O F D I L U T E D S P E C I E S I N T E R F A C E
255
= pe
1 2
D pe = pe
32
(8-21)
and the Tortuosity model, where the tortuosity expression is entered as user input:
pe D
D pe = -----------
(8-22)
These are readily used for both gaseous and liquid fluids along with various types of
particle shapes. For instance, the first model has been shown to fit mass transport in
soil-vapor and soil-moisture well.
Equation 8-18 can be solved for two types of boundary conditions at the interface
between the pellet surface and the fluid in this feature.
Continuous concentrations: assuming that all resistance to mass transfer is within the
pellet and no resistance to pellet-fluid mass transfer on 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 pellet system and the free fluid system though so-called reaction forces
in the finite element formulation.
Film resistance (mass flux): The flux of mass across the pellet-fluid interface into the
pellet is possibly rate determined on the bulk fluid side. The resistance is expressed
in terms of a film mass transfer coefficient, hDi, such that:
N i,inward = h D, i ( c ii c pe, i ) ,
(8-23)
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, the free fluid Equation 8-18 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
t
(8-24)
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 particles, the specific surface area exposed
to the free fluid is (Ref. 3):
3
S b = ------- ( 1 b )
r pe
256 |
C H A P T E R 8 : C H E M I C A L S P E C I E S TR A N S P O R T I N T E R F A C E S
(8-25)
The mass transfer coefficient in Equation 8-23 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 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 Frssling relation (Ref. 4):
Sh = 2 + 0.552Re
12
Sc
13
(8-26)
which was measured on particles in the size region 1 mm, the Rosner relation (Ref. 5)
Sh = Sc
0.4
( 0.4Re
12
+ 0.2Re
23
),
(8-27)
12
Sc
13
(8-28)
Sc = -------D
T H E O R Y F O R T H E TR A N S P O R T O F D I L U T E D S P E C I E S I N T E R F A C E
257
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. 325342, 1993.
2. P.V. Danckwerts, Continuous flow systems: Distribution of residence times,
Chem. Eng. Sci., vol. 2, no. 1, 1953.
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.
5. D.E Rosner, Transport Processes in Chemically Reacting Flow Systems, ISBN-13:
978-1483130262, Butterworth-Heinemann, 1986
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. 20172030, 1986.
7. C.W. Fetter, Contaminant Hydrogeology, Prentice Hall, 1999.
8. J. Bear and A. Verruijt, Modeling Groundwater Flow and Pollution, D. Reidel
Publishing, 1994.
9. J. Bear, Hydraulics of Groundwater, McGraw-Hill, 1979.
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. 695705, 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. 12001207, 1961.
13. I. Langmuir, Chemical Reactions at Low Temperatures, J. Amer. Chem. Soc.,
vol. 37, 1915.
258 |
C H A P T E R 8 : C H E M I C A L S P E C I E S TR A N S P O R T I N T E R F A C E S
Glossary
This Glossary of Terms contains finite element modeling terms specific to the
Microfluidics Module and its applications. For mathematical terms as well as
geometry and CAD terms specific to the COMSOL Multiphysics software and its
documentation, see the glossary in the COMSOL Multiphysics Reference Manual.
To find references in the documentation set where you can find more information
about a given term, see the index.
259
Glossary of Terms
absorption (gas) 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.
adsorption (molecule) Attachment of a molecule or atom to a solid surface.
Adsorption involves a chemical bond between the adsorbed species and the surface.
ALE See arbitrary Lagrangian-Eulerian method.
arbitrary Lagrangian-Eulerian (ALE) method A technique to formulate equations in a
mixed kinematical description. An ALE referential coordinate system is typically a mix
between the material (Lagrangian) and spatial (Eulerian) coordinate systems.
biosensor A general term for sensor devices that either detect biological substances or
gL
Bo = -------------
where is the density, g is the body force acceleration, L is the characteristic length,
and is the surface tension coefficient. The density can refer to a density difference
when considering buoyancy.
Brinkman equations A set of equations extending Darcys law in order to include the
surface tension for a two phase flow. At low capillary numbers flow in porous media is
dominated by surface tension. The capillary number, Ca, is given by:
260 |
CHAPTER 9: GLOSSARY
v
Ca = -----
where is the viscosity, v is the characteristic velocity of the problem, and is the
surface tension coefficient.
continuum flow Fluid flow which is well described by approximating the liquid as a
medium.
creeping flow Models the Navier-Stokes equations without the contribution of the
inertia term. This is often referred to as Stokes flow and is appropriate for use when
viscous flow is dominant, such as in very small channels or microfluidic applications.
crosswind diffusion A numerical technique for stabilization of the numeric solution to
a convection-dominated PDE by artificially adding diffusion perpendicular to the
direction of the streamlines.
Darcys Law Equation that gives the velocity vector as proportional to the pressure
K
Du = ------B
K
where K is the surface conductivity and KB is the bulk conductivity. For low Dukhin
numbers, surface conductivities can be neglected when modeling electrokinetic flows.
G L O S S A R Y O F TE R M S
261
electric double layer (EDL) Sometimes referred to as the Debye layer. At the contact
of a solid and a polar fluid (such as water), the solid acquires an electric charge. This
charge attracts ions within the fluid, and a narrow fluid layer of opposite charge, the
Stern layer, forms on the boundary. In addition, adjacent to the Stern layer, a wider
layer with the same charge as in the Stern layer forms in the fluid. Together, the Stern
layer and the wider layer (called the diffuse or Gouy-Chapman layer) form the electric
double layer. Due to the close distance between the charges, the Stern layer is fixed on
the surface, but the more distant diffuse layer can move.
electrohydrodynamics A general term describing phenomena that involve the
interaction between solid surfaces, ionic solutions, and applied electric and magnetic
fields. It is frequently used in microfluidic devices to manipulate fluids and move
particles for sample handling and chemical separation.
electrokinetics Study of the motion of charged particles under an applied electric field
in moving substances such as water.
electrokinetic flow Transport of fluid or charged particles within a fluid by means of
field.
electrothermal effects These effects occur in a conductive fluid where the
262 |
CHAPTER 9: GLOSSARY
electrothermal flow Fluid flow resulting from an applied nonuniform AC electric field
on a fluid. The Joule heating changes the fluids electrical properties locally, and that
effect, together with the power gradient of the AC electric field, results in fluid motion.
electrowetting The electrowetting effect describes the change in solid-electrolyte
contact angle that occurs when a potential difference is applied between the solid and
the electrolyte.
electrowetting-on-dielectric (EWOD) A form of electrowetting in which a thin
insulating layer separates the conducting solid surface from the electrolyte.
Etvs number See Bond number.
Eulerian frame A frame of reference with its coordinate axes fixed in space.
Ficks laws The first law relates the concentration gradients to the diffusive flux of low
concentration solute diluted in a solvent. The second law introduces the first law into
a differential material balance for the solute.
fluid-structure interaction (FSI) When a flow affects the deformation of a solid object
smaller than the mean free path (Knudsen number, Kn>10). In the free molecular flow
regime the gas molecules collide with the walls of the geometry much more frequently
than they collide with themselves.
fully developed laminar flow Laminar flow along a channel or pipe that has velocity
components only in the main direction of the flow. The velocity profile perpendicular
to the flow does not change downstream in the flow.
Gouy-Chapman layer See electric double layer (EDL).
Hagen-Poiseuille equation See Poiseuilles law.
Helmholtz-Smoluchowski equation Gives the velocity of a parallel electroosmotic flow
G L O S S A R Y O F TE R M S
263
gas flow is, in other words, the mean free path of the gas molecules compared to the
length scale of the flow. The following equation defines the Knudsen number Kn
where is the mean free path of the molecules and L is a length scale characteristic to
the flow.
Kn = ---L
Knudsen layer A layer of rarefied fluid flow that occurs within a few mean free paths
of the walls in a gas flow. The continuum Navier-Stokes equations break down in this
layer.
Lagrangian frame A frame of reference with its coordinate axes fixed in a reference
L
La = ---------2
where is the fluid density, is the surface tension coefficient, L is the characteristic
length scale of the problem, and is the viscosity. The Laplace number is directly
related to the Ohnesorge number, Oh, through the equation La=1/Oh2.
Mach number Ratio of the convective speed, v, to the speed of sound in the medium,
a. The Mach number, Ma, is defined by the equation:
v
Ma = --a
magnetohydrodynamics Fluid flow phenomena involving magnetic fields.
magnetophoresis Migration of magnetic or paramagnetic particles suspended in a
fluid as a result of an applied magnetic field. The field must be non-uniform in the case
of paramagnetic particles.
264 |
CHAPTER 9: GLOSSARY
gradients.
Marangoni number Ratio of thermal surface tension forces to viscous forces. The
Marangoni number, Mg, is given by:
d LT
Mg = -------- -----------dT
where is the surface tension coefficient, T is the absolute temperature (T is the
characteristic temperature difference), L is the characteristic length, and is the
thermal diffusivity (=/(cp)), where cp is the heat capacity at constant pressure, is
the thermal conductivity, and is the density of the fluid.
microfluidics Study of the behavior of fluids at the microscale. Also refers to MEMS
fluidic devices.
migration The transport of charged species in an electrolyte due to the electric force
diffusion, convection, and migration in an electric field. The equation is valid for
diluted electrolytes.
Ohnesorge number A dimensionless number relating the inertial and surface tension
forces to the viscous forces. Used to describe the breakup of liquid jets and sheets: at
low Ohnesorge and Reynolds numbers the Rayleigh instability occurs; at high
Ohnesorge and Reynolds numbers atomization occurs. The Ohnesorge number, Oh,
is given by:
Oh = --------------L
G L O S S A R Y O F TE R M S
265
where is the fluid density, is the surface tension coefficient, L is the characteristic
length scale of the problem, and is the viscosity. The Ohnesorge number is directly
related to the Laplace number, La, through the equation Oh=1/La1/2.
Peclet number A dimensionless number describing the ratio of convection to
diffusion in a fluid flow. The Peclet number, Pe, can be used to describe concentration
diffusion or heat diffusion. It is given by:
vL
Pe = ------D
where v is the convective velocity, L is the flow length scale, and D is the diffusion
constant.
Poiseuilles 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.
Reynolds number A dimensionless number classifying how laminar or turbulent a flow
is. The Reynolds number Re is a measure of the relative magnitude of the flows
viscous and inertial forces. It is defined by the following equation where is the fluid
density, is the dynamic viscosity, is its kinematic viscosity, v is a velocity
characteristic to the flow, and L is a length scale characteristic to the flow.
vL
vL
Re = ----------- = ------
slip flow Fluid flow that occurs when the Knudsen number, Kn, is in the range
266 |
CHAPTER 9: GLOSSARY
the medium including both pores and matrix. They are sometimes called Darcy
velocities, defined as volume flow rates per unit cross section of the medium.
surface tension Surface tension is a property of the surface of a liquid that allows it to
resist an external force. Equivalently it can be through of as the energy per unit area of
the liquid surface. It is caused by asymmetries in the cohesive forces between molecules
at the surface of the liquid.
Surataman number See Laplace number.
transitional flow Fluid flow that occurs when the Knudsen number, Kn, is in the range
0.1<Kn<10. In this regime the flow is so rarefied that continuum equations break
down completely. However collisions between the molecules are still important, so free
molecular flow is not applicable.
Weber number A dimensionless number describing the ratio of inertial forces to
v L
We = ------------
where is the density of the fluid, v is the characteristic velocity of the flow, L is the
characteristic length scale, and is the surface tension coefficient.
zeta potential The potential of at the interface between the electric double layer and
G L O S S A R Y O F TE R M S
267
268 |
CHAPTER 9: GLOSSARY
I n d e x
2D axisymmetric models
AC electric fields 52
Boussinesq approximation 89
AC electroosmosis 48
AC/DC Module 51
creeping flow 64
Carreau model 88
electroosmosis 48
CFL number
electroosmotic velocity 93
settings 63
electrowetting lens 53
theory 107
laminar flow 63
114
15, 44
dielectrophoresis 50
magnetophoresis 51
common settings 20
complex permittivity 50
compressible flow 86
boundary conditions
concentration (node)
continuum flow 56
convection 247
124
coupling interfaces 46
INDEX|
269
electrokinetic flows 46
electrokinetics 46
electroosmosis 47
condition 69
theory 189
theory 92
electrophoresis 46, 49
Debye length 47
Debye-Hckel approximation 47
dielectrophoresis (DEP) 49
emailing COMSOL 23
entrance length 71
Equilibrium Reaction
discontinuous Galerkin
documentation 21
exit length 74
domain nodes
270 | I N D E X
flux (node)
127
theory 194
tion 103
Internet resources 21
tion) 76
isotropic diffusion 32
Joule heating 52
heat transfer 54
Knudsen layer 57
Knudsen number 14, 34, 56
92
Henrys function 49
Hygroscopic Swelling 233
I
lab-on-a-chip devices 14
laminar flow
Reynolds number, and 86
laminar flow interface 60
theory 82
laminar inflow (inlet boundary condition)
71
INDEX|
271
terface 112
interface 114
MEMS, definition 14
Laplace number 35
microfluidics
definition 14
modeling techniques 28
stabilization techniques 32
fluid flow 80
modeling
electrohydrodynamics 46
line source
microfluidics 28
44
local
MPH-files 23
multiphase flow
Mach number 34
magnetohydrodynamics 46
magnetophoresis 46, 51
Newtonian model 88
135
mass sources
fluid flow 100
272 | I N D E X
Marangoni number 35
non-Newtonian fluids 84
stress condition) 77
point nodes
bers and 31
O
Ohnesorge number 35
124
single-phase flow 75
outflow (node)
point source
single-phase flow 73
pair nodes
183
124
194
(node) 227
229
126
Peclet number 31
settings 63
theory 107
pseudoplastic fluids 88
INDEX|
273
numbers and 32
reactions (node)
73
standard settings 20
Stern layer 47
Stokes equations 63
bers and 31
streamline diffusion
fluid flow 102
selecting
dia 192
dia 247
Surataman number 35
laminar flow 60
theory 82
theory 90
slip
Maxwells boundary condition 57
slip flow 57
slip flow interface 196
theory 203
slip velocity, wall boundary condition 70
theory 91
theory 90
274 | I N D E X
93
laminar flow 82
[Link] variable 87
[Link] variable 67
single-phase flow 82
Weber number 35
websites, COMSOL 24
Youngs equation 52
235
variables
level set interface 160
niterCMP 107
phase field interface 165
[Link] 87
viscous slip, wall boundary condition 70
viscous stress tensors, theory 87
viscous stress, theory 97
volume averages 191
volume force (node) 67
Brinkman equations 180
free and porous media flow 186
W wall (node)
single-phase flow 68
INDEX|
275
276 | I N D E X