Models - Cfd.boiling Water [PDF]

Zitiervorschau

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 As

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 fmix 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