Solved with COMSOL Multiphysics 4.3b
Boiling Water
Note: This model requires the CFD Module and the Heat Transfer Module.
Introduction
Boiling flow is an example of phase transition of the first kind and can be initiated by
raising the temperature of a liquid above its saturation temperature. It is possible to
accomplish this in many ways; two of the most common ways are by applying an
external heat flux to a solid surface in contact with a liquid or by reducing the pressure
in the surrounding environment. There are three distinct regimes that characterize
boiling induced by a heated surface: nucleate, transition, and film. The boiling regime
depends on the excess temperature and the magnitude of the applied heat flux to the
surface. The excess temperature is defined as the difference between the temperature
of the solid surface and the saturation temperature of the liquid.
This model shows how to solve a boiling flow problem with the Laminar Two-Phase
Flow, Phase Field user interface.
NUCLEATE BOILING
For water at atmospheric pressure, nucleate boiling occurs for an excess temperature
between 5 °C and 30 °C. The heat flux required to sustain this surface temperature is
between 104 W/m2 and 106 W/m2. In this regime nucleation sites, which form due
to imperfections in the solid surface, separate and expand due to mass transfer from
the liquid phase to the vapor phase. The bubbles begin to rise due to buoyancy forces
and may merge with other isolated bubbles to form jets or columns. Because the
detaching bubbles are quickly replaced by liquid, the effective thermal conductivity of
the fluid layer close to surface is high. This means that even though the applied heat
flux is very high, the excess temperature remains low. Nucleate boiling is most efficient
at the critical heat flux (see Figure 1). At this point liquid can still rapidly replace the
vapor produced, in effect continuously wetting the surface. This results in a heat
©2013 COMSOL
1 |
B O I L I N G WA T E R
Solved with COMSOL Multiphysics 4.3b
transfer coefficient of more than 104 W/(m2·K), much higher than any heat transfer
coefficient that occurs due to convection alone.
Convection
Nucleate
Film
Transition
Critical heat flux
Leidenfrost point
Onset of boiling
Figure 1: Boiling curve for water at atmospheric pressure (101,325 Pa). The same surface
heat flux can result in different excess temperatures depending on the boiling regime (Ref.
2).
TR A N S I T I O N B O I L I N G
To generate excess temperatures above 30 °C the heat flux absorbed must actually
decrease below the critical heat flux. This is because a thin layer of low thermal
conductivity vapor begins to insulate the solid surface from the liquid, reducing the
heat transfer coefficient between the fluid and solid surface. The region is referred to
as transition boiling because there are continuous local oscillations between nucleate
and film boiling regimes.
FILM BOILING
To sustain film boiling, the excess temperature and applied heat flux must be greater
than the Leidenfrost point. For water at atmospheric pressure, this corresponds to a
heat flux of around 4·104 W/m2 and an excess temperature of around 120 °C. A layer
of vapor continuously insulates the hot surface from the liquid and periodic shedding
of vapor bubbles from the vapor/liquid interface occurs. Unlike nucleate and
2 |
B O I L I N G WA T E R
©2013 COMSOL
Solved with COMSOL Multiphysics 4.3b
transition boiling, film boiling is relatively stable in nature and easier to simulate.
Beyond the Leidenfrost point, an increase in the applied heat flux leads directly to an
increase in excess temperature. At very high excess temperatures, conduction is not the
only mechanism of heat transfer from the hot surface to the liquid, and radiation effects
must be considered.
Figure 1 shows how the applied heat flux varies with the excess temperature for the
three main boiling regimes. At heat fluxes and excess temperatures below the onset of
boiling, natural convection occurs.
Model Definition
PROBLEM STATEMENT
The geometry and initial conditions for a boiling flow problem are shown in Figure 2.
Initially, there are pockets of vapor in each of the two heated cavities at the base of a
tank of water. This leads to stable film boiling when a heat flux above the Leidenfrost
point is applied.
A surface heat flux of 105 W/m2 is applied, which is above the Leidenfrost point. An
excess temperature on the surface of around 600 K is expected (1000 K total
temperature) in order to sustain film boiling.
©2013 COMSOL
3 |
B O I L I N G WA T E R
Solved with COMSOL Multiphysics 4.3b
Water vapor, 1 atm
Initial vapor pocket
Water
q
q
Figure 2: Initially two heated cavities are filled with a pocket of vapor. An inward heat
flux, q, is applied to the cavity surfaces.
PROBLEM FORMULATION
The physical behavior of a boiling flow is driven by the interface dynamics. The
governing partial differential equations are the standard Navier-Stokes and
convection/conduction equations. The boundary conditions are, however, rather
complicated because the interface between the liquid and vapor is moving. First, the
boiling flow problem is formulated with the exact equations and boundary conditions.
Then, a series of approximations are made so that the problem can be solved on a fixed
mesh where the interface is tracked by the phase field equation.
Domain Equations
The velocity field and pressure for the liquid phase are described by the incompressible
Navier-Stokes equations:
4 |
B O I L I N G WA T E R
©2013 COMSOL
Solved with COMSOL Multiphysics 4.3b
∂u L
T
ρ L ---------- + ρ L ( u L ⋅ ∇ )u L = ∇ ⋅ [ – p L I + μ L ( ∇u L + ( ∇u L ) ) ] + ρ L g
∂t
(1)
∇ ⋅ uL = 0
(2)
where ρL is the fluid density (kg/m3), uL is the fluid velocity (m/s) and μL is the
viscosity (Pa·s). Subscript L denotes liquid phase. For the vapor phase, the
compressible Navier-Stokes equations are solved:
∂u v
ρ v ---------- + ρ v ( u v ⋅ ∇ )u v =
∂t
(3)
2
∇ ⋅ – p v I + η v ( ∇u v + ( ∇u v ) ) – --- μ v ( ∇ ⋅ u )I + ρ v g
3
T
∂ρ v
-------- + ∇ ⋅ ( ρ v u v ) = 0
∂t
(4)
where subscript v denotes liquid phase. The heat equation is solved only in the vapor
phase:
∂T v
ρ v C p ---------- + ρ v C p ( u v ⋅ ∇ )T v = ∇ ⋅ κ v ∇T v
∂t
(5)
where Cp is the specific heat capacity (J/(kg·K)) and κv is the thermal conductivity of
the vapor (W/(m·K)). The heat conduction equation is only solved in the vapor phase
because the temperature at the liquid/vapor interface is set to the saturation
temperature. This results in a constant temperature throughout the liquid phase and
so it is not necessary to solve the heat equation there.
Boundary Conditions
The boundary conditions for a boiling flow model are rather complicated. It is
important to realize that the interface velocity, liquid velocity and vapor velocity are
not necessarily equal:
·
m
u int = u L – ------ n
ρL
Where n is the unit normal vector to the interface directed from the liquid phase to
the vapor phase. The natural boundary condition on the interface for the vapor phase
is:
©2013 COMSOL
5 |
B O I L I N G WA T E R
Solved with COMSOL Multiphysics 4.3b
ρv
·
n ⋅ ρ v u v = m  1 – ------ + ( n ⋅ ρ v u L )

ρ L
(6)
·
where m is the rate of vaporization (kg/m2·s). There is a lot of physics in Equation 6
and a short discussion is necessary. It is obvious that if there is no phase change,
·
m = 0 , so uv = uL = uint and the velocity is continuous across the interface.
Furthermore, if the density of the liquid and vapor phase is the same, the first term on
the right hand side is zero and the velocity is continuous across the interface. When
phase change is occurring, and the liquid density is not equal to the vapor density, the
first term on the right hand side results in an inward flow velocity normal to the
interface. The second term on the right hand side results in a velocity in the outward
direction normal to the interface. This leads to a discontinuity in the velocity field
across the interface, the velocity on the liquid in the outward direction normal from
the vapor phase side being higher than the vapor side.
There are three forces that act on the liquid at the interface and so the natural
boundary condition for the liquid is
T
·
n ⋅ [ – p L I + μ L ( ∇u L + ( ∇u L ) ) ] = m ( u L – u v ) + σκn
T 2
+ n ⋅ – p v I + μ v ∇u v + ( ∇u v ) – --- μ v ( ∇ ⋅ u )I
3
(7)
which results from a force balance on the interface. The first term on the right hand
side represents a reaction force due to the acceleration of the vapor away from the
liquid surface. This force is nearly always negligible. The second term is surface tension
and the last term on the right hand side is the sum of the pressure and viscous forces
acting on the liquid from the vapor. The mass flux leaving the liquid surface leads to
an increase in pressure of the vapor. The pressure exerts a force on the liquid surface
and the vapor region begins to expand. The presence of the surface tension force leads
to a discontinuity in pressure across the interface.
In the energy equation, the temperature at the interface is fixed at the saturation
temperature, which may be a function of pressure:
T = T sat ( p )
(8)
The mass flux leaving the interface can then be evaluated from the conductive heat
flux:
6 |
B O I L I N G WA T E R
©2013 COMSOL
Solved with COMSOL Multiphysics 4.3b
Mw
·
m = –  -------------  n ⋅ κ v ∇T v
ΔH vl
(9)
where Mw is the molecular weight of the vapor (kg/mol) and ΔHvl is the enthalpy of
vaporization (J/mol). This approximation is obtained by neglecting the kinetic
energies and work due to viscous forces, see Ref. 4.
Equation 1 through Equation 5 along with boundary conditions Equation 6 through
Equation 8 represent a complete description of a boiling flow. It is possible, in
principle, to solve these equations using the arbitrary Lagrangian-Eulerian (ALE)
method, but no topological changes can occur. By making some suitable assumptions,
it is possible to solve the problem on a fixed mesh with the Laminar Two-Phase Flow,
Phase Field user interface.
PROBLEM FORMULATION—PHASE FIELD METHOD
When moving to a fixed mesh, a number of approximations must be made. They
essentially involve rewriting the boundary conditions above as volumetric sources or
sinks. The equations governing the interface dynamics of a two-phase flow can be
described by the Cahn-Hilliard equation. As stated in Ref. 1, the equation for the
phase field variable is modified to allow for the change of phase:
γλ
∂φ
· V f, v V f, L
+ u ⋅ ∇φ – m δ  ---------- + ----------- = ∇ ⋅ -----2- ∇ψ
∂t
ρv
ρL
ε
where φ is the dimensionless phase field variable such that – 1 ≤ φ ≤ 1 . The volume
fraction of the vapor phase is Vf,v, and the volume fraction of the liquid is Vf,L. The
quantity λ is the mixing energy density (N) and ε is a capillary width that scales with
the thickness of the interface (m). These two parameters are related to the surface
tension coefficient via
2 2λ
σ = ----------- --3 ε
and γ is the mobility (m3·s/kg). 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. The
equation governing ψ is:
2
2
ψ = – ∇ ⋅ ε ∇φ + ( φ – 1 )φ
©2013 COMSOL
7 |
B O I L I N G WA T E R
Solved with COMSOL Multiphysics 4.3b
The quantity, δ (1/m) is a smoothed representation of the interface between the two
phases. In the Laminar Two-Phase Flow, Phase Field user interface, it is defined as
∇φ
δ = 6V f ( 1 – V f ) ---------2
The momentum equation includes surface tension effects as a volumetric body force:
ρ
∂u
T
+ ρ ( u ⋅ ∇ )u = ∇ ⋅ [ – p I + μ ( ∇u + ( ∇u ) ) ] + ρg + G∇φ
∂t
where G is the chemical potential (Pa). The continuity equation is modified to account
for the phase change from liquid to vapor:
1
1
·
∇ ⋅ u = m δ  ------ – ------
 ρ v ρ L
The problem now is how to define an expression for the rate of phase change.
Equation 9 cannot be used because the peak temperature gradient does not coincide
with the interface. This leads to a substantial underestimate of the mass flux emanating
from the surface. Instead, the mass flux can be approximated by the following
expression:
( T – T sat )
Mw
·
m = –  ------------- n ⋅ κ v ∇T v ≈ Cρ L ------------------------- ΔH vl
T sat
(10)
where C is a constant (m/s). This expression is analogous to specifying a heat transfer
coefficient on an external boundary of a heat transfer problem. The mass flux appears
in the energy equation as:
·
m δ ΔH vl
∂T
ρC p ------- + ρC p ( u ⋅ ∇ )T = ∇ ⋅ κ∇T – ---------------------∂t
Mw
(11)
The combination of Equation 10 and Equation 11 naturally forces the interface
temperature to the saturation temperature. Practically, choose C to be large enough
that the temperature at the interface remains at the saturation temperature but not so
large that numerical instabilities result.
Physical Properties
The two fluids considered are water and water vapor. The saturation temperature of
the liquid is 373 K (100 °C) and the pressure in the surrounding environment is
101325 Pa. The liquid density is 1000 kg/m3, and the viscosity is 1·10−3 Pa·s. The
8 |
B O I L I N G WA T E R
©2013 COMSOL
Solved with COMSOL Multiphysics 4.3b
vapor density is given by the ideal gas law, and a viscosity of 4·10−5 Pa·s is specified.
This is slightly higher than the actual viscosity of water vapor but results in a more
stable implementation. The surface tension coefficient is 0.0588 N/m. The thermal
conductivity and specific heat capacity are computed as functions of the volume
fraction of the two phases:
κ = ( κ L – κ v )V f, L + κ v
C p = ( C p, L – C p, v )V f, L + C p, v
The walls of the heater are allowed to radiate to free space and to each other.
Results and Discussion
As the water begins to boil, the surface area (length in 2D) of the interface begins to
increase. The surface area can be computed by integrating the expression for δ over the
modeling domain:
As =
 δ dV
V
©2013 COMSOL
9 |
B O I L I N G WA T E R
Solved with COMSOL Multiphysics 4.3b
Because the vapor is continuously being released and detaching from the surface, the
surface area of the interface oscillates as pockets of vapor are produced expand and
then penetrate the free surface:
Figure 3: Plot of interface length versus time. The shedding frequency of the vapor bubbles
is periodic.
The structure of the free surface is examined at various time steps in Figure 4. After
about 0.5 seconds, bubbles of vapor begin to shed from the cavities. The bubbles of
vapor rise through the liquid and penetrate the free surface. Boiling, by nature is
extremely sensitive to the input parameters. For example, lowering the applied heat
flux will not lead to stable film boiling and the size and period with which the vapor
bubbles are shed is strongly dependent on the latent heat and surface tension
coefficient.
10 |
B O I L I N G WAT E R
©2013 COMSOL
Solved with COMSOL Multiphysics 4.3b
The average temperature on the heated surface is around 1080 °C, 700 °C above the
boiling point. This qualitatively agrees with the boiling curve (see Figure 1).
t=0s
t=0.52s
t=0.64s
t=0.92s
t=1.27s
t=1.53s
Figure 4: Plot of fluid volume fractions at different time steps.
You can use the following two definitions to compute the total free energy of the
system:
F = σA s
F =
  --2- λ ∇φ
1
2
(12)
2
λ
2
+ --------2- ( φ – 1 )  dV

4ε
(13)
Any difference in the two values indicates that the Laminar Two-Phase Flow, Phase
Field user interface is not exactly conserving the total free energy in the system.
Practically, this means that “mass loss” is occurring. Figure 5 plots the total free energy
for the two different methods. Initially, the two methods agree very well, but as the
periodic shedding of vapor bubbles begins, deviations begin to occur. The “mass
losses” mainly occur when the vapor bubbles penetrate the free surface. In these
©2013 COMSOL
11 |
B O I L I N G WA T E R
Solved with COMSOL Multiphysics 4.3b
instances very sharp interfaces are formed, and there is inadequate mesh resolution to
accurately capture this effect.
Figure 5: Plot of total free energy using Equation 12 (blue) and Equation 13 (green).
limitations of the Model
There were a number of assumptions used in this model and their effect is important
to understand:
• The approximate expression for the rate of vapor generation generally results in a
rate of vaporization that is lower than expected.
• The isotropic artificial diffusion in the heat transfer equation results in a higher than
expected transport of energy from the wall to the vapor surface. This results in an
unphysical increase in the rate of vaporization. Additionally, this makes the solution
highly mesh dependent.
• The model is solved in the 2D xy-plane which results in an unphysical
implementation of the surface tension force. As a result, only qualitative information
about the topological evolution of the vapor bubbles can be extracted from the
model.
• The model neglects radiation from the hot surface to the liquid interface. In reality,
this effect enhances the rate of vaporization on the surface.
12 |
B O I L I N G WAT E R
©2013 COMSOL
Solved with COMSOL Multiphysics 4.3b
References
1. Y. Sun and C. Beckermann, “Diffuse interface modeling of two-phase flows based
on averaging: mass and momentum equations,” Physica D, vol. 198, pp. 281–308,
2004.
2. F.P. Incropera and D.P. DeWitt, Fundamentals of Heat and Mass transfer, John
Wiley & Sons, Inc., 2002.
3. G. Son and V.K. Dhir, “A Level Set Method for Analysis of Film Boiling on an
Immersed Solid Surface,” Numerical Heat Transfer, Part B, vol. 52, pp. 153–177,
2007.
4. D. Jamet, “Diffuse interface models in fluid mechanics,”
http://pmc.polytechnique.fr/mp/GDR/docu/Jamet.pdf.
Model Library path: CFD_Module/Multiphase_Tutorials/boiling_water
Modeling Instructions
MODEL WIZARD
1 Go to the Model Wizard window.
2 Click the 2D button.
3 Click Next.
4 In the Add physics tree, select Fluid Flow>Multiphase Flow>Two-Phase Flow, Phase
Field>Laminar Two-Phase Flow, Phase Field (tpf).
5 Click Add Selected.
6 In the Add physics tree, select Heat Transfer>Heat Transfer in Fluids (ht).
7 Click Add Selected.
8 Click Next.
9 Find the Studies subsection. In the tree, select Preset Studies for Selected
Physics>Transient with Initialization.
10 Click Finish.
©2013 COMSOL
13 |
B O I L I N G WA T E R
Solved with COMSOL Multiphysics 4.3b
GEOMETRY 1
Rectangle 1
1 In the Model Builder window, under Model 1 right-click Geometry 1 and choose
Rectangle.
2 In the Rectangle settings window, locate the Size section.
3 In the Width edit field, type 0.04.
4 In the Height edit field, type 0.05.
5 Locate the Position section. In the x edit field, type -0.02.
6 Click to expand the Layers section. In the table, enter the following settings:
Layer name
Thickness (m)
Layer 1
0.025
7 Click the Build Selected button.
Polygon 1
1 In the Model Builder window, right-click Geometry 1 and choose Polygon.
2 In the Polygon settings window, locate the Coordinates section.
3 In the x edit field, type -0.0095
-0.0125 -0.0025 -0.0055 -0.0095.
4 In the y edit field, type 0 -0.0075 -0.0075 0 0.
5 Click the Build Selected button.
Circle 1
1 Right-click Geometry 1 and choose Circle.
2 In the Circle settings window, locate the Size and Shape section.
3 In the Radius edit field, type 0.0075.
4 Locate the Position section. In the x edit field, type -0.0075.
5 In the y edit field, type -0.01.
6 Click the Build Selected button.
7 Click the Zoom Extents button on the Graphics toolbar.
Compose 1
1 Right-click Geometry 1 and choose Boolean Operations>Compose.
2 Select the objects pol1 and c1 only.
3 In the Compose settings window, locate the Compose section.
4 In the Set formula edit field, type pol1+pol1*c1.
14 |
B O I L I N G WAT E R
©2013 COMSOL
Solved with COMSOL Multiphysics 4.3b
5 Click the Build Selected button.
Copy 1
1 Right-click Geometry 1 and choose Transforms>Copy.
2 Select the object co1 only.
3 In the Copy settings window, locate the Displacement section.
4 In the x edit field, type 0.015.
5 Click the Build Selected button.
Form Union
6 In the Model Builder window, under Model 1>Geometry 1 right-click Form Union and
choose Build Selected.
The model geometry is now complete and should look like the figure below.
GLOBAL DEFINITIONS
Parameters
1 In the Model Builder window, right-click Global Definitions and choose Parameters.
2 In the Parameters settings window, locate the Parameters section.
©2013 COMSOL
15 |
B O I L I N G WA T E R
Solved with COMSOL Multiphysics 4.3b
3 In the table, enter the following settings:
Name
Expression
Description
Mw
0.018[kg/mol]
Molecular weight of water
L
42000[J/mol]/Mw
Latent heat of vaporization
Cp_l
4200[J/(kg*K)]
Heat capacity, liquid
Cp_g
1840[J/(kg*K)]
Heat capacity, gas
T_sat
373[K]
Boiling point, 1atm
k_l
0.63[W/m/K]
Thermal conductivity, liquid
q
1E5[W/m^2]
Inward heat flux
p0
101325[Pa]
Atmospheric pressure
C
0.03[m/s]
Constant in expression for
mass flux
DEFINITIONS
Variables 1
1 In the Model Builder window, under Model 1 right-click Definitions and choose
Variables.
2 In the Variables settings window, locate the Variables section.
3 In the table, enter the following settings:
16 |
Name
Expression
kappa
(k_l-k_g)*tpf.Vf2+k_g
Thermal conductivity
Cp
(Cp_l-Cp_g)*tpf.Vf2+Cp_g
Specific heat
k_g
8.3154E-05[W/m/K/
K]*T-7.4556E-03[W/m/K]
Thermal conductivity of
vapor
rho1
(p+p0)*Mw/8.314[J/mol/K]/
T
Density of vapor
q1
q+0.1*q*sin(2*pi*0.2*t[1/
s])
Heat flux, cavity 1
q2
q-0.1*q*sin(2*pi*0.2*t[1/
s])
Heat flux, cavity 2
mdot
C*tpf.rho2*(T-T_sat)*flc2
hs((T-T_sat)[1/K],0.001)/
T_sat
Expression for rate of
vaporization
delta
6*tpf.Vf2*(1-tpf.Vf2)*0.5
*sqrt(phipfx^2+phipfy^2+e
ps)
Interface delta function
B O I L I N G WAT E R
Description
©2013 COMSOL
Solved with COMSOL Multiphysics 4.3b
Name
Expression
Description
phi_source
-mdot*delta*(tpf.Vf1/
tpf.rho1+tpf.Vf2/
tpf.rho2)
Source term in the level
set equation
usource
mdot*delta*(1/tpf.rho1-1/
tpf.rho2)
Source in continuity
equation
Qs
-mdot*L*delta
Heat sink in energy
equation
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
1 In the Model Builder window’s toolbar, click the Show button and select Stabilization
in the menu.
2 In the Model Builder window, under Model 1 click Laminar Two-Phase Flow, Phase
Field.
3 In the Laminar Two-Phase Flow, Phase Field settings window, click to expand the
Consistent Stabilization section.
4 Find the Navier-Stokes equations subsection. Clear the Streamline diffusion check
box.
5 Clear the Crosswind diffusion check box.
6 Click to expand the Inconsistent Stabilization section. Select the Isotropic diffusion
check box.
Fluid Properties 1
1 In the Model Builder window, expand the Laminar Two-Phase Flow, Phase Field node,
then click Fluid Properties 1.
2 In the Fluid Properties settings window, locate the Fluid 1 Properties section.
3 From the ρ1 list, choose User defined. In the associated edit field, type rho1.
4 From the μ1 list, choose User defined. In the associated edit field, type 4e-5.
5 Locate the Fluid 2 Properties section. From the ρ2 list, choose User defined. In the
associated edit field, type 1000.
6 From the μ2 list, choose User defined. In the associated edit field, type 1e-3.
7 Locate the Surface Tension section. From the Surface tension coefficient list, choose
User defined. In the σ edit field, type 0.0588.
8 Locate the Phase Field Parameters section. In the εpf edit field, type 4e-4.
9 In the χ edit field, type 10.
©2013 COMSOL
17 |
B O I L I N G WA T E R
Solved with COMSOL Multiphysics 4.3b
Initial Interface 1
1 In the Model Builder window, under Model 1>Laminar Two-Phase Flow, Phase Field
click Initial Interface 1.
2 Select Boundaries 4 and 22–25 only.
Initial Values 2
1 In the Model Builder window, right-click Laminar Two-Phase Flow, Phase Field and
choose Initial Values.
2 Select Domains 1, 4, and 6 only.
3 In the Initial Values settings window, locate the Initial Values section.
4 Click the Fluid 2 button.
Gravity 1
1 Right-click Laminar Two-Phase Flow, Phase Field and choose Gravity.
2 In the Gravity settings window, locate the Domain Selection section.
3 From the Selection list, choose All domains.
Outlet 1
1 Right-click Laminar Two-Phase Flow, Phase Field and choose Outlet.
2 Select Boundary 5 only.
3 In the Model Builder window’s toolbar, click the Show button and select Advanced
Physics Options in the menu.
Weak Contribution 1
1 Right-click Laminar Two-Phase Flow, Phase Field and choose the domain setting
More>Weak Contribution.
2 In the Weak Contribution settings window, locate the Domain Selection section.
3 From the Selection list, choose All domains.
4 Locate the Weak Contribution section. In the Weak expression edit field, type
test(psi)*phi_source+test(p)*usource.
H E A T TR A N S F E R I N F L U I D S
1 In the Model Builder window, under Model 1 click Heat Transfer in Fluids.
2 In the Heat Transfer in Fluids settings window, locate the Physical Model section.
3 Select the Surface-to-surface radiation check box.
4 Click to expand the Consistent Stabilization section. Clear the Streamline diffusion
check box.
18 |
B O I L I N G WAT E R
©2013 COMSOL
Solved with COMSOL Multiphysics 4.3b
5 Click to expand the Inconsistent Stabilization section. Select the Isotropic diffusion
check box.
Heat Transfer in Fluids 1
1 In the Model Builder window, expand the Heat Transfer in Fluids node, then click Heat
Transfer in Fluids 1.
2 In the Heat Transfer in Fluids settings window, locate the Model Inputs section.
3 In the pA edit field, type p0.
4 From the u list, choose Velocity field (tpf/fp1).
5 Locate the Heat Conduction, Fluid section. From the k list, choose User defined. In
the associated edit field, type kappa.
6 Locate the Thermodynamics, Fluid section. From the ρ list, choose User defined. In
the associated edit field, type tpf.rho.
7 From the Cp list, choose User defined. In the associated edit field, type Cp.
8 From the γ list, choose User defined.
Initial Values 1
1 In the Model Builder window, under Model 1>Heat Transfer in Fluids click Initial Values
1.
2 In the Initial Values settings window, locate the Initial Values section.
3 In the T edit field, type T_sat.
Heat Source 1
1 In the Model Builder window, right-click Heat Transfer in Fluids and choose Heat
Source.
2 In the Heat Source settings window, locate the Domain Selection section.
3 From the Selection list, choose All domains.
4 Locate the Heat Source section. In the Q edit field, type Qs.
Heat Flux 1
1 Right-click Heat Transfer in Fluids and choose Heat Flux.
2 Select Boundaries 1–3, 11, 18, 20, and 21 only.
3 In the Heat Flux settings window, locate the Heat Flux section.
4 Click the Inward heat flux button.
5 In the h edit field, type 10.
6 In the Text edit field, type T_sat.
©2013 COMSOL
19 |
B O I L I N G WA T E R
Solved with COMSOL Multiphysics 4.3b
Heat Flux 2
1 Right-click Heat Transfer in Fluids and choose Heat Flux.
2 Select Boundaries 6–8, 10, and 12 only.
3 In the Heat Flux settings window, locate the Heat Flux section.
4 In the q0 edit field, type q1.
Heat Flux 3
1 Right-click Heat Transfer in Fluids and choose Heat Flux.
2 Select Boundaries 13–15, 17, and 19 only.
3 In the Heat Flux settings window, locate the Heat Flux section.
4 In the q0 edit field, type q2.
Surface-to-Surface Radiation 1
1 Right-click Heat Transfer in Fluids and choose Surface-to-Surface
Radiation>Surface-to-Surface Radiation.
2 Select Boundaries 6–8, 10, 12–15, 17, and 19 only.
3 In the Surface-to-Surface Radiation settings window, locate the Radiation Settings
section.
4 From the Radiation direction list, choose Positive normal direction.
5 Locate the Ambient section. In the Tamb edit field, type 0.
6 Locate the Surface Emissivity section. From the ε list, choose User defined. In the
associated edit field, type 1.
Outflow 1
1 Right-click Heat Transfer in Fluids and choose Outflow.
2 Select Boundary 5 only.
MESH 1
In the Model Builder window, under Model 1 right-click Mesh 1 and choose Free
Triangular.
Size
1 In the Model Builder window, under Model 1>Mesh 1 click Size.
2 In the Size settings window, locate the Element Size section.
3 Click the Custom button.
4 Locate the Element Size Parameters section. In the Maximum element size edit field,
type 1.25e-3.
20 |
B O I L I N G WAT E R
©2013 COMSOL
Solved with COMSOL Multiphysics 4.3b
Size 1
1 In the Model Builder window, under Model 1>Mesh 1 right-click Free Triangular 1 and
choose Size.
2 In the Size settings window, locate the Geometric Entity Selection section.
3 From the Geometric entity level list, choose Domain.
4 Select Domains 3–6 only.
5 In the Size settings window, locate the Element Size section.
6 Click the Custom button.
7 Locate the Element Size Parameters section. Select the Maximum element size check
box.
8 In the associated edit field, type 5e-4.
Size 2
1 Right-click Free Triangular 1 and choose Size.
2 In the Size settings window, locate the Geometric Entity Selection section.
3 From the Geometric entity level list, choose Boundary.
4 Select Boundary 4 only.
5 Locate the Element Size section. Click the Custom button.
6 Locate the Element Size Parameters section. Select the Maximum element size check
box.
7 In the associated edit field, type 1e-3.
©2013 COMSOL
21 |
B O I L I N G WA T E R
Solved with COMSOL Multiphysics 4.3b
8 Click the Build All button.
STUDY 1
Step 2: Time Dependent
1 In the Model Builder window, expand the Study 1 node, then click Step 2: Time
Dependent.
2 In the Time Dependent settings window, locate the Study Settings section.
3 In the Times edit field, type range(0,0.01,2).
4 In the Model Builder window, right-click Study 1 and choose Show Default Solver.
5 Expand the Study 1>Solver Configurations node.
Solver 1
1 In the Model Builder window, expand the Study 1>Solver Configurations>Solver 1
node, then click Time-Dependent Solver 1.
2 In the Time-Dependent Solver settings window, click to expand the Absolute Tolerance
section.
3 In the Tolerance edit field, type 5e-4.
4 In the Model Builder window, right-click Study 1 and choose Compute.
22 |
B O I L I N G WAT E R
©2013 COMSOL
Solved with COMSOL Multiphysics 4.3b
RESULTS
Volume Fraction (tpf)
The default plot shows the phase field function at the end of the simulated time
interval. To visualize the initial phase field function (Figure 2), plot the function at
t = 0 instead.
1 In the 2D Plot Group settings window, locate the Data section.
2 From the Time list, choose 0.
3 Click the Plot button.
4 Click the Zoom Extents button on the Graphics toolbar.
Similarly, you can reproduce the plots in Figure 4 by choosing the corresponding
time values from the Time list.
To reproduce the plot in Figure 3, do the following:
Derived Values
1 In the Model Builder window, under Results right-click Derived Values and choose
Integration>Surface Integration.
2 In the Surface Integration settings window, locate the Selection section.
3 From the Selection list, choose All domains.
4 Locate the Expression section. In the Expression edit field, type delta.
5 Click the Evaluate button.
1D Plot Group 6
1 In the Model Builder window, right-click Results and choose 1D Plot Group.
2 Right-click 1D Plot Group 6 and choose Table Graph.
3 Right-click Results>1D Plot Group 6>Table Graph 1 and choose Plot.
4 In the Model Builder window, click 1D Plot Group 6.
5 In the 1D Plot Group settings window, locate the Plot Settings section.
6 Select the x-axis label check box.
7 In the associated edit field, type Time (s).
8 Select the y-axis label check box.
9 In the associated edit field, type Interface length (m).
To reproduce the plot in Figure 5, do the following:
©2013 COMSOL
23 |
B O I L I N G WA T E R
Solved with COMSOL Multiphysics 4.3b
Derived Values
1 In the Model Builder window, under Results right-click Derived Values and choose
Integration>Surface Integration.
2 In the Surface Integration settings window, locate the Selection section.
3 From the Selection list, choose All domains.
4 Locate the Expression section. In the Expression edit field, type delta*tpf.sigma.
5 Click the Evaluate button.
6 In the Model Builder window, right-click Derived Values and choose
Integration>Surface Integration.
7 In the Surface Integration settings window, locate the Selection section.
8 From the Selection list, choose All domains.
9 Locate the Expression section. In the Expression edit field, type
0.5*tpf.lam*(phipfx^2+phipfy^2)+tpf.lam*(phipf^2-1)^2/
(4*tpf.epsilon_pf^2).
10 Click the Evaluate button.
1D Plot Group 7
1 In the Model Builder window, right-click Results and choose 1D Plot Group.
2 Right-click 1D Plot Group 7 and choose Table Graph.
3 In the Table Graph settings window, locate the Data section.
4 From the Table list, choose Table 2.
5 Click to expand the Legends section. Select the Show legends check box.
6 From the Legends list, choose Manual.
7 Click the Plot button.
8 In the table, enter the following settings:
Legends
sigma A<sub>s</sub>
9 In the Model Builder window, right-click 1D Plot Group 7 and choose Table Graph.
10 In the Table Graph settings window, locate the Data section.
11 From the Table list, choose Table 3.
12 Locate the Legends section. Select the Show legends check box.
13 From the Legends list, choose Manual.
14 Click the Plot button.
24 |
B O I L I N G WAT E R
©2013 COMSOL
Solved with COMSOL Multiphysics 4.3b
15 In the table, enter the following settings:
Legends
integral f<sub>mix</sub> dV
16 In the Model Builder window, click 1D Plot Group 7.
17 In the 1D Plot Group settings window, click to expand the Legend section.
18 From the Position list, choose Lower right.
©2013 COMSOL
25 |
B O I L I N G WA T E R
Solved with COMSOL Multiphysics 4.3b
26 |
B O I L I N G WAT E R
©2013 COMSOL