Heat Transfer Module: User's Guide
Heat Transfer Module: User's Guide
User’s Guide
Heat Transfer Module User’s Guide
© 1998–2021 COMSOL
Protected by patents listed on www.comsol.com/patents, or see Help>About COMSOL Multiphysics on the
File menu in the COMSOL Desktop for a less detailed lists of U.S. Patents that may apply. Patents pending.
This Documentation and the Programs described herein are furnished under the COMSOL Software License
Agreement (www.comsol.com/comsol-license-agreement) and may be used or copied only under the terms
of the license agreement.
COMSOL, the COMSOL logo, COMSOL Multiphysics, COMSOL Desktop, COMSOL Compiler,
COMSOL Server, and LiveLink are either registered trademarks or trademarks of COMSOL AB. All other
trademarks are the property of their respective owners, and COMSOL AB and its subsidiaries and products
are not affiliated with, endorsed by, sponsored by, or supported by those trademark owners. For a list of such
trademark owners, see www.comsol.com/trademarks.
Version: COMSOL 6.0
Contact Information
Visit the Contact COMSOL page at www.comsol.com/contact to submit general inquiries
or search for an address and phone number. You can also visit the Worldwide Sales Offices
page at www.comsol.com/contact/offices for address and contact information.
If you need to contact Support, an online request form is located at the COMSOL Access
page at www.comsol.com/support/case. Other useful links include:
Chapter 1: Introduction
Chapter 2: Notations
Symbols 46
CONTENTS |3
Moisture Transport Variables 80
Predefined Variables . . . . . . . . . . . . . . . . . . . . . 80
Moist Air Properties . . . . . . . . . . . . . . . . . . . . . 82
Domain Moisture Fluxes . . . . . . . . . . . . . . . . . . . . 83
Boundary Moisture Fluxes . . . . . . . . . . . . . . . . . . . 85
Domain Moisture Source. . . . . . . . . . . . . . . . . . . . 86
Mass Balance . . . . . . . . . . . . . . . . . . . . . . . . 86
4 | CONTENTS
Plotting and Evaluating Results in Layered Materials 118
Plotting Along and Through the Layered Material . . . . . . . . . 118
Built-in Operators for Evaluation in the Layered Material . . . . . . 119
CONTENTS |5
References 161
6 | CONTENTS
Theory for Heat Transfer in Building Materials 208
CONTENTS |7
Discrete Ordinates Method Implementation in 2D . . . . . . . . . 278
P1 Approximation Theory . . . . . . . . . . . . . . . . . . 280
Radiation in Absorbing-Scattering Media Theory . . . . . . . . . . 283
Polychromatic Radiation . . . . . . . . . . . . . . . . . . . 285
Radiative Beam in Absorbing Media Theory . . . . . . . . . . . . 285
Rosseland Approximation Theory . . . . . . . . . . . . . . . 286
8 | CONTENTS
Heat Transfer Coefficients — External Forced Convection . . . . . . 335
Heat Transfer Coefficients — Internal Forced Convection . . . . . . 338
Using the Heat and Mass Transfer Analogy for the Evaluation of
Moisture Transfer Coefficients . . . . . . . . . . . . . . . 339
References 355
CONTENTS |9
Versions of the Moisture Transport Physics Interface . . . . . . . . 369
Benefits of the Different Moisture Transport Interfaces . . . . . . . 370
Additional Physics Options . . . . . . . . . . . . . . . . . . 370
Settings for the Heat Transfer Interface . . . . . . . . . . . . . 371
Feature Nodes for the Heat Transfer Interface . . . . . . . . . . 374
Settings for the Heat Transfer in Shells Interface . . . . . . . . . . 377
Settings for the Moisture Transport Interface . . . . . . . . . . . 382
Feature Nodes for the Moisture Transport Interface . . . . . . . . 384
10 | C O N T E N T S
The Lumped Thermal System Interface 412
Feature Nodes for the Lumped Thermal System Interface . . . . . . 413
CONTENTS | 11
Cross Section . . . . . . . . . . . . . . . . . . . . . . . 461
Fluid . . . . . . . . . . . . . . . . . . . . . . . . . . 462
Fluid (Porous Medium) . . . . . . . . . . . . . . . . . . . 467
Geothermal Heating . . . . . . . . . . . . . . . . . . . . 472
Heat Source . . . . . . . . . . . . . . . . . . . . . . . 473
Immobile Fluids (Porous Medium) . . . . . . . . . . . . . . . 477
Initial Values . . . . . . . . . . . . . . . . . . . . . . . 479
Initial Values (Radiative Beam in Absorbing Medium Interface) . . . . . 480
Initial Values (Radiation in Participating Media and Radiation in
Absorbing-Scattering Media Interfaces) . . . . . . . . . . . . 481
Irreversible Transformation . . . . . . . . . . . . . . . . . . 483
Isothermal Domain . . . . . . . . . . . . . . . . . . . . . 486
Liquid Water (Moist Porous Medium) . . . . . . . . . . . . . . 488
Moist Air (Heat Transfer Interface) . . . . . . . . . . . . . . . 489
Moist Air (Moist Porous Medium) . . . . . . . . . . . . . . . 492
Moist Porous Medium . . . . . . . . . . . . . . . . . . . . 492
Opacity (Surface-to-Surface Radiation Interface) . . . . . . . . . . 496
Optically Thick Participating Medium . . . . . . . . . . . . . . 497
Out-of-Plane Heat Flux . . . . . . . . . . . . . . . . . . . 499
Out-of-Plane Radiation . . . . . . . . . . . . . . . . . . . 501
Participating Medium (Radiation in Participating Media Interface) . . . . 503
Pellets (Porous Medium) . . . . . . . . . . . . . . . . . . . 509
Pellet-Fluid Interface (Porous Medium) . . . . . . . . . . . . . 511
Phase Change Material. . . . . . . . . . . . . . . . . . . . 512
Porous Matrix (Porous Medium, Moist Porous Medium) . . . . . . . 516
Porous Medium . . . . . . . . . . . . . . . . . . . . . . 520
Pressure Work . . . . . . . . . . . . . . . . . . . . . . 526
Radiative Source . . . . . . . . . . . . . . . . . . . . . . 527
Shape Memory Alloy . . . . . . . . . . . . . . . . . . . . 529
Solid . . . . . . . . . . . . . . . . . . . . . . . . . . 533
Thermal Damage . . . . . . . . . . . . . . . . . . . . . . 537
Thermal Dispersion . . . . . . . . . . . . . . . . . . . . . 541
Thermoelastic Damping . . . . . . . . . . . . . . . . . . . 543
Thickness. . . . . . . . . . . . . . . . . . . . . . . . . 544
Translational Motion . . . . . . . . . . . . . . . . . . . . 545
Viscous Dissipation . . . . . . . . . . . . . . . . . . . . . 546
12 | C O N T E N T S
Boundary Features 548
Boundary Heat Source. . . . . . . . . . . . . . . . . . . . 550
Continuity . . . . . . . . . . . . . . . . . . . . . . . . 553
Continuity (Radiation in Participating Media and Radiation in
Absorbing-Scattering Media Interfaces) . . . . . . . . . . . . 554
Continuity on Interior Boundary (Radiation in Participating Media and
Radiation in Absorbing-Scattering Media Interfaces) . . . . . . . 554
Deposited Beam Power . . . . . . . . . . . . . . . . . . . 555
Diffuse Mirror (Surface-to-Surface Radiation Interface) . . . . . . . 557
Diffuse Surface (Surface-to-Surface Radiation Interface) . . . . . . . 558
External Temperature (Thin Layer, Thin Film, Fracture) . . . . . . . 565
Fracture (Heat Transfer Interface) and Porous Medium (Heat Transfer
in Shells Interface) . . . . . . . . . . . . . . . . . . . . 566
Harmonic Perturbation . . . . . . . . . . . . . . . . . . . 571
Heat Flux. . . . . . . . . . . . . . . . . . . . . . . . . 572
Heat Source (Heat Transfer in Shells Interface) . . . . . . . . . . 577
Heat Source (Thin Layer, Thin Film, Fracture) . . . . . . . . . . . 579
Incident Intensity (Radiation in Participating Media and Radiation in
Absorbing-Scattering Media Interfaces) . . . . . . . . . . . . 582
Incident Intensity (Radiative Beam in Absorbing Medium Interface) . . . 584
Inflow . . . . . . . . . . . . . . . . . . . . . . . . . . 586
Initial Values (Thin Layer, Thin Film, Fracture, and Heat Transfer in
Shells Interface) . . . . . . . . . . . . . . . . . . . . . 587
Initial Values (Surface-to-Surface Radiation Interface) . . . . . . . . 590
Isothermal Domain Interface . . . . . . . . . . . . . . . . . 591
Layer Opacity (Surface-to-Surface Radiation Interface) . . . . . . . 595
Line Heat Source on Axis . . . . . . . . . . . . . . . . . . 596
Local Thermal Nonequilibrium Boundary. . . . . . . . . . . . . 597
Lumped System Connector . . . . . . . . . . . . . . . . . . 598
Opaque Surface (Surface-to-Surface Radiation Interface) . . . . . . . 599
Opaque Surface (Radiation in Participating Media and Radiation in
Absorbing-Scattering Media Interfaces) . . . . . . . . . . . . 602
Opaque Surface (Radiative Beam in Absorbing Medium Interface) . . . 606
Open Boundary . . . . . . . . . . . . . . . . . . . . . . 607
Outflow . . . . . . . . . . . . . . . . . . . . . . . . . 608
Periodic Condition (Heat Transfer Interface) . . . . . . . . . . . 609
Periodic Condition (Radiation in Participating Media and Radiation in
Absorbing-Scattering Media Interfaces) . . . . . . . . . . . . 610
CONTENTS | 13
Phase Change Interface . . . . . . . . . . . . . . . . . . . 611
Phase Change Interface, Exterior . . . . . . . . . . . . . . . . 614
Prescribed Radiosity (Surface-to-Surface Radiation Interface) . . . . . 617
Radiation Group (Surface-to-Surface Radiation Interface) . . . . . . 621
Semitransparent Surface (Radiation in Participating Media and
Radiation in Absorbing-Scattering Media Interfaces) . . . . . . . 623
Semitransparent Surface (Surface-to-Surface Radiation Interface). . . . 628
Surface-to-Ambient Radiation (Heat Transfer Interface) . . . . . . . 632
Symmetry (Heat Transfer Interface). . . . . . . . . . . . . . . 634
Symmetry (Radiation in Participating Media and Radiation in
Absorbing-Scattering Media Interfaces) . . . . . . . . . . . . 635
Temperature . . . . . . . . . . . . . . . . . . . . . . . 635
Thermal Contact . . . . . . . . . . . . . . . . . . . . . . 637
Thermal Insulation . . . . . . . . . . . . . . . . . . . . . 642
Thin Film (Heat Transfer Interface) and Fluid (Heat Transfer in
Shells Interface) . . . . . . . . . . . . . . . . . . . . . 643
Thin Layer (Heat Transfer Interface) and Solid (Heat Transfer in
Shells Interface) . . . . . . . . . . . . . . . . . . . . . 647
Transparent Surface (Radiative Beam in Absorbing Medium Interface) . . 653
14 | C O N T E N T S
Shell Continuity (Heat Transfer Interface) and Continuity (Heat
Transfer in Shells Interface) . . . . . . . . . . . . . . . . 688
Surface-to-Ambient Radiation (Thin Layer, Thin Film, Fracture, and
Heat Transfer in Shells Interface) . . . . . . . . . . . . . . 689
Symmetry (Thin Layer, Thin Film, Fracture, and Heat Transfer in
Shells Interface) . . . . . . . . . . . . . . . . . . . . . 692
Temperature (Thin Layer, Thin Film, Fracture, and Heat Transfer in
Shells Interface) . . . . . . . . . . . . . . . . . . . . . 694
Thermal Insulation (Thin Layer, Thin Film, Fracture, and Heat
Transfer in Shells Interface) . . . . . . . . . . . . . . . . 696
Thin Rod . . . . . . . . . . . . . . . . . . . . . . . . . 698
CONTENTS | 15
Chapter 7: The Moisture Transport Features
16 | C O N T E N T S
The Nonisothermal Flow, Laminar Flow and Turbulent Flow
Interfaces . . . . . . . . . . . . . . . . . . . . . . . 783
The Conjugate Heat Transfer, Laminar Flow and Turbulent Flow
Interfaces . . . . . . . . . . . . . . . . . . . . . . . 784
Settings for Physics Interfaces and Coupling Features . . . . . . . . 785
Coupling Feature . . . . . . . . . . . . . . . . . . . . . . 786
Physics Interface Features . . . . . . . . . . . . . . . . . . 786
Preset Studies . . . . . . . . . . . . . . . . . . . . . . . 787
CONTENTS | 17
The Thermoelectric Effect Interface 800
About The Thermoelectric Effect Interface . . . . . . . . . . . . 800
Settings for Physics Interfaces and Coupling Features . . . . . . . . 801
Coupling Features . . . . . . . . . . . . . . . . . . . . . 802
Physics Interface Features . . . . . . . . . . . . . . . . . . 803
18 | C O N T E N T S
Coupling Feature . . . . . . . . . . . . . . . . . . . . . . 826
Index 865
CONTENTS | 19
20 | C O N T E N T S
1
Introduction
This guide describes the Heat Transfer Module, an optional package that extends
the COMSOL Multiphysics® modeling environment with customized physics
interfaces for the analysis of heat transfer.
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.
21
About the Heat Transfer Module
In this section:
Heat transfer is involved in almost every kind of physical process, and can in fact be the
limiting factor for many processes. Therefore, its study is of vital importance, and the
need for powerful heat transfer analysis tools is virtually universal. Furthermore, heat
transfer often appears together with, or as a result of, other physical phenomena.
The modeling of heat transfer effects has become increasingly important in product
design including areas such as electronics, automotive, and medical industries.
Computer simulation has allowed engineers and researchers to optimize process
efficiency and explore new designs, while at the same time reducing the need for costly
experimental trials.
22 | CHAPTER 1: INTRODUCTION
How the Heat Transfer Module Improves Your Modeling
The Heat Transfer Module has been developed to greatly expand upon the base
capabilities available in COMSOL Multiphysics. The module supports all fundamental
mechanisms including conductive, convective, and radiative heat transfer. Using the
physics interfaces in this module along with the inherent multiphysics capabilities of
COMSOL Multiphysics, you can model a temperature field in parallel with other
physics — a versatile combination increasing the accuracy and predicting power of your
models.
This book introduces the basic modeling process. The different physics interfaces are
described and the modeling strategy for various cases is discussed. These sections cover
different combinations of conductive, convective, and radiative heat transfer. This
guide also reviews special modeling techniques for thin layers, thin shells, participating
media, and out-of-plane heat transfer. Throughout the guide the topics and examples
increase in complexity by combining several heat transfer mechanisms and also by
coupling these to physics interfaces describing fluid flow — conjugate heat transfer.
Most of the examples involve multiple heat transfer mechanisms and are often coupled
to other physical phenomena, for example, fluid dynamics, moisture transport, or
electromagnetics. The authors developed several state-of-the art examples by
reproducing examples that have appeared in international scientific journals. See
Where Do I Access the Documentation and Application Libraries?.
Moisture Transport
Turbulent Flow
24 | CHAPTER 1: INTRODUCTION
PHYSICS INTERFACE ICON TAG SPACE AVAILABLE PRESET STUDY
DIMENSION TYPE
Fluid Flow
Single-Phase Flow
Turbulent Flow
26 | CHAPTER 1: INTRODUCTION
PHYSICS INTERFACE ICON TAG SPACE AVAILABLE PRESET STUDY
DIMENSION TYPE
Turbulent Flow
Thin Structures
28 | CHAPTER 1: INTRODUCTION
PHYSICS INTERFACE ICON TAG SPACE AVAILABLE PRESET STUDY
DIMENSION TYPE
30 | CHAPTER 1: INTRODUCTION
PHYSICS INTERFACE ICON TAG SPACE AVAILABLE PRESET STUDY
DIMENSION TYPE
DOMAIN
CHEMICAL SPECIES TRANSPORT>MOISTURE TRANSPORT
Laminar Flow(2) — u, p, φ √ √
Turbulent flow, Algebraic — u, p, G, √ √
yPlus(2) yPlus, φ
Turbulent flow, L-VEL(2) — u, p, G, √ √
uPlus, φ
Turbulent Flow, k-ε(2) — u, p, k, ep, φ √ √
(2)
Turbulent Flow, Low Re k-ε — u, p, k, ep, G, √ √
φ
Porous Media Flow(2) — u, p, φ √ √
FLUID FLOW>SINGLE-PHASE FLOW
32 | CHAPTER 1: INTRODUCTION
TABLE 1-1: HEAT TRANSFER MODULE DEPENDENT VARIABLES AND PRESET STUDY AVAILABILITY
1
PHYSICS INTERFACE NAME DEPENDENT PRESET STUDIES
VARIABLES
DOMAIN
Turbulent Flow, Low Re k-ε spf u, p, k, ep, G √ √
FLUID FLOW>NONISOTHERMAL FLOW
Laminar Flow(2) — u, p, T √ √
Turbulent flow, Algebraic — u, p, G, √ √
yPlus(2) yPlus, T
Turbulent flow, L-VEL(2) — u, p, G, √ √
uPlus, T
Turbulent Flow, k-ε(2) — u, p, k, ep, T √ √
(2)
Turbulent Flow, Low Re k-ε — u, p, k, ep, G, √ √
T
HEAT TRANSFER
Laminar Flow(2) — u, p, T √ √
Turbulent flow, Algebraic — u, p, G, √ √
yPlus(2) yPlus, T
Turbulent flow, L-VEL(2) — u, p, G, √ √
uPlus, T
Turbulent Flow, k-ε(2) — u, p, k, ep, T √ √
(2)
Turbulent Flow, Low Re k-ε — u, p, k, ep, G, √ √
T
DOMAIN
HEAT TRANSFER>RADIATION
Joule Heating(2) — T, V √ √ √
Lumped Thermal System lts not applicable √ √
HEAT TRANSFER>THIN STRUCURES
34 | CHAPTER 1: INTRODUCTION
TABLE 1-1: HEAT TRANSFER MODULE DEPENDENT VARIABLES AND PRESET STUDY AVAILABILITY
1
PHYSICS INTERFACE NAME DEPENDENT PRESET STUDIES
VARIABLES
DOMAIN
Moist Air(2) — T, φ √ √
(2)
Moist Porous Media — T, φ √ √
(2)
Building Materials — T, φ √ √
HEAT TRANSFER>HEAT AND MOISTURE TRANSPORT>HEAT AND MOISTURE FLOW
Laminar Flow(2) — u, p, T, φ √ √
Turbulent flow, Algebraic — u, p, G, √ √
yPlus(2) yPlus, T, φ
Turbulent flow, L-VEL(2) — u, p, G, √ √
uPlus, T, φ
Turbulent Flow, k-ε(2) — u, p, k, ep, T, √ √
φ
Turbulent Flow, Low Re k-ε(2) — u, p, k, ep, G, √ √
T, φ
Porous Media Flow(2) — u, p, T, φ √ √
Heat Transfer in Porous Media ht T √ √ √ √
Local Thermal Nonequilibrium ht T √ √ √ √
Heat Transfer in Packed Beds ht T √ √ √ √
Bioheat Transfer ht T √ √ √ √
(2)
Thermoelectric Effect — T √ √ √ √
1
Custom studies are also available based on the physics interface.
2 Multiphysics interfaces.
36 | CHAPTER 1: INTRODUCTION
Opening Topic-Based Help
The Help window is useful as it is connected to the features in the COMSOL Desktop.
To learn more about a node in the Model Builder, or a window on the Desktop, click
to highlight a node or window, then press F1 to open the Help window, which then
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.
• In the Model Builder or Physics Builder click a node or window and then
press F1.
• On the main toolbar, click the Help ( ) button.
• From the main menu, select Help>Help.
• Press Ctrl+F1.
• From the File menu select Help>Documentation ( ).
• Press Ctrl+F1.
• On the main toolbar, click the Documentation ( ) button.
• From the main menu, select Help>Documentation.
Once the Application Libraries window is opened, you can search by name or browse
under a module folder name. Click to view a summary of the model or application and
its properties, including options to open it or its associated PDF document.
To include the latest versions of model examples, from the Help menu
select ( ) Update COMSOL Application Library.
38 | CHAPTER 1: INTRODUCTION
CONTACTING COMSOL BY EMAIL
For general product information, contact COMSOL at [email protected].
40 | CHAPTER 1: INTRODUCTION
After the establishment of the heat balance equation from the energy conservation laws
in Foundations of the General Heat Transfer Equation, the various versions of the heat
equation solved in COMSOL Multiphysics are presented in the following sections:
Then the theory related to multiphysics interfaces is described in the section Theory
for the Heat Transfer Multiphysics Couplings.
Finally, topics related to specific features or variables are treated in the sections Theory
for Thermal Contact, Theory for Heat Transfer in Moist Air, Out-of-Plane Heat
Transfer, Convective Heat Transfer Correlations, Equivalent Thermal Conductivity
Correlations, Temperature Dependence of Surface Tension, Heat Flux and Heat
Balance, and Frames for the Heat Transfer Equations.
The sections The Heat Transfer in Solids Interface, The Heat Transfer in Fluids
Interface, and The Heat Transfer in Solids and Fluids Interface discuss modeling heat
transfer in solids and fluids.
The particular case of heat transfer in moist air and building materials is considered in
the sections The Heat Transfer in Moist Air Interface and The Heat Transfer in
Building Materials Interface.
The section The Bioheat Transfer Interface discusses modeling heat transfer within
biological tissue using the Bioheat Transfer interface.
The sections The Heat Transfer in Shells Interface, The Heat Transfer in Films
Interface, and The Heat Transfer in Fractures Interface describe the physics interfaces
which are suitable for solving thermal conduction, convection, and radiation problems
in layered materials defined on boundaries.
The section The Lumped Thermal System Interface describes the modeling of heat
transfer in a system using a thermal network representation.
Finally, the sections The Moisture Transport in Air Interface, The Moisture Transport
in Porous Media Interface, and The Moisture Transport in Building Materials
Interface describe the modeling of moisture transfer in a porous medium through
moisture storage, vapor diffusion and capillary moisture flows; or in air, through
convection and diffusion.
42 | CHAPTER 1: INTRODUCTION
The chapter The Nonisothermal Flow and Conjugate Heat Transfer Interfaces
describes the multiphysics versions of both the Nonisothermal Flow Laminar Flow and
Turbulent Flow interfaces found under the Fluid Flow branch, which are identical to
the Conjugate Heat Transfer interfaces. Each section describes the applicable physics
interfaces in detail and concludes with the underlying theory.
The section The Heat Transfer with Surface-to-Surface Radiation Interface describes
the predefined multiphysics interface used to model heat transfer by conduction,
convection, and radiation in a transparent media.
The section The Heat Transfer with Radiation in Participating Media Interface
describes the predefined multiphysics interface used to model heat transfer by
conduction, convection, and radiation in semitransparent media. The radiative
intensity equations are approximated by the Discrete Ordinates Method or the P1
Approximation. When no emission should be considered, see the section The Heat
Transfer with Radiation in Absorbing-Scattering Media Interface.
The section The Heat Transfer with Radiative Beam in Absorbing Media Interface
describes the predefined multiphysics interface used to model heat transfer by
conduction, convection, and radiation in semitransparent media. The Beer-Lambert
law is used for the approximation of the radiative intensity.
The section The Thermoelectric Effect Interface describes the predefined multiphysics
interface used to model the Peltier-Seebeck-Thomson effect.
The section The Heat and Moisture Transport Interfaces describes the predefined
multiphysics interfaces used to model coupled heat and moisture transport either in
building materials, by taking into account heat and moisture storage, latent heat
effects, and liquid and convective transport of moisture; or in moist air by convection
and diffusion of moisture and heat.
The section The Moisture Flow Interfaces describes the predefined multiphysics
interfaces used to model moisture transport in air by laminar flow (in free and porous
media) and turbulent flows (in free media).
The section The Heat and Moisture Flow Interfaces describes the predefined
multiphysics interfaces used to model heat transfer and moisture transport in air by
laminar and turbulent flows, and in porous media by laminar flows.
44 | CHAPTER 1: INTRODUCTION
2
Notations
This chapter introduces the notations used in the remaining of the guide. The
notations are listed in alphabetical order and grouped in two tables for Latin and
Greek symbols, respectively.
For each entry the SI unit and a short description are given.
45
Symbols
LATIN SYMBOLS
46 | CHAPTER 2: NOTATIONS
NOTATION SI UNIT DESCRIPTION
SYMBOLS | 47
NOTATION SI UNIT DESCRIPTION
48 | CHAPTER 2: NOTATIONS
NOTATION SI UNIT DESCRIPTION
SYMBOLS | 49
NOTATION SI UNIT DESCRIPTION
50 | CHAPTER 2: NOTATIONS
NOTATION SI UNIT DESCRIPTION
SYMBOLS | 51
NOTATION SI UNIT DESCRIPTION
52 | CHAPTER 2: NOTATIONS
NOTATION SI UNIT DESCRIPTION
SYMBOLS | 53
NOTATION SI UNIT DESCRIPTION
rh m Hydraulic radius
rl m Rod radius
rpe m Average pellets radius
Rs J/(kg·K) Specific gas constant
Rt,s K·m2/W Surface thermal resistance
Rt K/W Thermal resistance
Sb 1/m Specific surface area
sd m Vapor diffusion equivalent air layer thickness
S N/m2 Second Piola-Kirchhoff stress tensor
S V/K Seebeck coefficient
ScT dimensionless Turbulent Schmidt number
Si dimensionless Unit vector of discrete direction in space, i-th component
(angular space discretization)
Sp dimensionless Sparrow number
T K Temperature
T0 K Equilibrium temperature
T’ K Complex amplitude of harmonic perturbation
+
T dimensionless Dimensionless temperature
Tamb K Ambient temperature
Tamb, d K Ambient temperature, downside
Tamb, u K Ambient temperature, upside
Tb K Arterial blood temperature
Td K Temperature, downside
td, c s Damage time, cryogenic analysis
Td, c K Damage temperature, cryogenic analysis
td, h s Damage time, hyperthermia analysis
Td, h K Damage temperature, hyperthermia analysis
Text K External temperature
Text, d K Out-of-plane external temperature, downside
Text, u K Out-of-plane external temperature, upside
Text, z K Out-of-plane external temperature, 1D
Tf K Temperature, fluid phase
54 | CHAPTER 2: NOTATIONS
NOTATION SI UNIT DESCRIPTION
SYMBOLS | 55
GREEK SYMBOLS
56 | CHAPTER 2: NOTATIONS
NOTATION SI UNIT DESCRIPTION
SYMBOLS | 57
NOTATION SI UNIT DESCRIPTION
58 | CHAPTER 2: NOTATIONS
NOTATION SI UNIT DESCRIPTION
SYMBOLS | 59
60 | CHAPTER 2: NOTATIONS
3
61
Heat Transfer Variables
In this section:
• Predefined Variables
• Global Variables for Heat and Energy Balance
• Domain Fluxes
• Out-of-Plane Domain Fluxes
• Boundary Fluxes (Heat Transfer Interface)
• Internal Boundary Heat Fluxes
• Domain Heat Sources
• Boundary Heat Sources
• Line and Point Heat Sources
• Moist Air Variables
Predefined Variables
This section lists some predefined variables that are available to evaluate heat fluxes,
sources, and integral quantities used in heat and energy balance. All the variable names
begin with the physics interface name (the prefix). By default the Heat Transfer
interface prefix is ht, and the Heat Transfer in Shells interface prefix is htlsh. As an
example, you can access the variable named tflux using ht.tflux (as long as the
physics interface is named ht).
TABLE 3-1: PREDEFINED HEAT BALANCE VARIABLES
Global variables
dEiInt Total Accumulated Heat Rate Global
ntfluxInt Total Net Heat Rate Global
WstrInt Total Stress Power Global
heatBalance Heat Balance Global
Domain variables
tflux Total Heat Flux (Heat Transfer Domains
Interface)
Global variables
dEi0Int Total Accumulated Energy Rate Global
ntefluxInt Total Net Energy Rate Global
WInt Total Work Source Global
energyBalan Energy Balance Global
ce
Domain variables
teflux Total Energy Flux (Heat Transfer Domains
Interface)
dflux Conductive Heat Flux (Heat Domains
Transfer Interface)
thflux Total Enthalpy Flux (Heat Transfer Domains
Interface)
Boundary variables
nteflux Normal Total Energy Flux (Heat Boundaries
Transfer Interface)
ndflux Normal Conductive Heat Flux Boundaries
(Heat Transfer Interface)
nthflux Normal Total Enthalpy Flux (Heat Boundaries
Transfer Interface)
q0 Inward Heat Flux (Heat Transfer Boundaries
Interface)
nteflux_u Internal normal total energy flux, Interior boundaries
nteflux_d upside and downside (Internal
Boundary Energy Fluxes)
ndflux_u Internal normal conductive heat flux, Interior boundaries
ndflux_d upside and downside (Internal
Boundary Heat Fluxes)
nthflux_u Internal normal enthalpy flux, upside Interior boundaries
nthflux_d and downside (Internal Boundary
Energy Fluxes)
TABLE 3-3: PREDEFINED HEAT SOURCE VARIABLES
Global variables
QInt Total Heat Source Global
Domain variables
Global variables
Tref Reference temperature Global
Domain variables
alpha Damage indicator (see Theory for Domains
Bioheat Transfer)
alphanecr Instant tissue necrosis indicator (see Domains
Theory for Bioheat Transfer)
theta_d Fraction of damage (see Theory for Domains
Bioheat Transfer)
theta_d_sm Fraction of total damage (see Theory Domains
for Bioheat Transfer)
T_dp Dew Point Temperature Domains
T_eq Equivalent Temperature Domains
psat Saturation Pressure Domains
phi Relative Humidity Domains
Lv Latent Heat of Evaporation Domains
Boundary variables
h Heat transfer coefficient (see Exterior boundaries
Convective Heat Transfer
Correlations)
d
dEiInt =
dt Ω ρEi dω
TOTAL NET HEAT RATE
The total net heat rate, ntfluxInt, is the integral of Total Heat Flux (Heat Transfer
Interface) over all external boundaries. In the case of a fluid domain, it reads:
This indicates the sum of incoming and outgoing total heat flux through the system.
d
dEi0Int =
dt Ω ρEi0 dω
where the total internal energy, Ei0, is defined as
u⋅u
E i0 = E i + ------------
2
This indicates the sum of incoming and outgoing total energy flux through the system.
T 2
K = μ ( ∇u + ( ∇u ) ) – --- μ ( ∇ ⋅ u )I ⋅ n
3
and K ⋅ u ∂Ω int
is the jump of viscous stress at interior boundaries.
HEAT BALANCE
According to Equation 4-200, the following equality between COMSOL Multiphysics
variables holds:
This is the most general form that can be used for time-dependent models. At
steady-state the formula is simplified. The accumulated heat rate equals zero, so the
total net heat rate (the sum of incoming and outgoing heat rates) should correspond
to the heat and work sources:
The sign convention used in COMSOL Multiphysics for QInt is positive when energy
is produced (as for a heater) and negative when energy is consumed (as for a cooler).
For WstrInt, the losses that heat up the system are negative and the gains that cool
down the system are positive.
For stationary models with convection by an incompressible flow, the heat balance
becomes:
ntfluxInt = QInt
∂Ω ext
ρuE i ⋅ n dσ – ∂Ω ext
k∇T ⋅ n dσ = Q Int
At steady state, and without any additional heat source (QInt equal to zero), the
integral of the net energy flux on all boundaries of the flow domain, ntefluxInt,
vanishes. On the other hand, the corresponding integral of the net heat flux does not,
in general, vanish. It corresponds instead to the losses from mass and momentum
equations, such as WstrInt for pressure work and viscous dissipation in fluids. Hence,
energy is the conserved quantity, not heat.
Domain Fluxes
On domains the fluxes are vector quantities. The definition can vary depending on the
active physics nodes and selected properties.
See Radiative Heat Flux (Heat Transfer Interface) to evaluate the radiative
heat flux.
For solid domains — for example, the solid and biological tissue domains — the total
heat flux is defined as:
tflux = dflux
For fluid domains (for example, Fluid), the total heat flux is defined as:
dflux = – k eff ∇T
For heat transfer in fluids with turbulent flow, keff = k + kT, where kT is the turbulent
thermal conductivity.
For heat transfer in porous media, keff is the effective conductivity computed from the
solid and fluid conductivities.
For heat transfer in building materials, a latent heat source due to evaporation is
included in the conductive heat flux variable:
When considering thermoelectric effect, the thermoelectric heat flux is also included
in the diffusive flux, see Equation 4-150.
cflux = ρuE i
The convective heat flux may be oriented in the opposite direction of the velocity
field’s direction, when the internal energy, Ei, has a negative value. This happens when
the sensible enthalpy (variation from reference enthalpy Href) is negative. Href is set to
0 J/kg at pref (1 atm), cref (0 mol/m3), and Tref (293.15 K) in COMSOL
Multiphysics. See Thermodynamic Description of Heat Transfer for details.
thflux = ρuH 0
u⋅u
H 0 = H + ------------
2
• In 2D:
4 4
upside: rflux_u = ε u σ ( T amb, u – T )
4 4
downside: rflux_d = ε d σ ( T amb, d – T )
• In 1D:
4 4
rflux_z = ε z σ ( T amb, z – T )
• In 2D:
upside: q0_u = h u ( T ext, u – T )
• In 1D:
q0_z = h z ( T ext, z – T )
• In 2D:
upside: tflux_u = q0_u + rflux_u
• In 1D:
Frames for the Heat Transfer Equations for a description of spatial and
material frames.
q0 = h ( T ext – T )
where Text is the external temperature defined in the Heat Flux feature.
rflux = ε ( G – e b ( T ) ) + q r, net
ndflux_u = uflux_spatial ( T )
ndflux_d = dflux_spatial ( T )
Frames for the Heat Transfer Equations for a description of spatial and
material frames.
nthflux_u = up ( thflux ) ⋅ un
• Q’s, which are the heat sources added by the Heat Source (described for the Heat
Transfer interface) and Electromagnetic Heating (described for the Joule Heating
interface in the COMSOL Multiphysics Reference Manual) features.
• Qmet, which is the metabolic heat source added by the Bioheat feature.
• Qdmg, which is the cooling source added by the Irreversible Transformation feature.
• Qr, which is the radiative heat source added by the Heat Transfer with Radiation in
Participating Media, Heat Transfer with Radiation in Absorbing-Scattering Media,
and Heat Transfer with Radiative Beam in Absorbing Media multiphysics features.
• Qgeo, which is the geothermal heat source added by the Geothermal Heating
feature.
• Qevap, which is the latent heat source added by the Heat and Moisture coupling
node.
• Qb, which is the boundary heat source added by the Boundary Heat Source
boundary condition.
• Qsh, which is the boundary heat source added by the Electromagnetic Heating
condition (described for the Joule Heating interface in the COMSOL Multiphysics
Reference Manual).
• Qs, which is the boundary heat source added by a Layer Heat Source subfeature of
a thin layer, see Heat Source (Thin Layer, Thin Film, Fracture).
The sum of the point heat sources is available in a variable called Qptot (SI unit: W).
When the presence of water vapor is accounted for in the model, other temperatures
may be considered, depending on vapor pressure.
where phi is the Relative Humidity variable. See Saturation State for the definition of
saturation pressure psat as a function of temperature. See also Saturation Pressure for
the definition of the variable psat.
EQUIVALENT TEMPERATURE
The equivalent temperature is obtained by adiabatically condensing all the water vapor
of a sample of air with initial vapor pressure pv. In this process, the latent heat decrease
due to total removal of the vapor is balanced by a increase of the sensible heat and
temperature.
phi ⋅ psat
T_eq = T + -----------------------------
γ Teq
where phi is the Relative Humidity, and γTeq (SI unit: Pa/K) is the psychrometer
constant, defined in Ref. 1 by:
where p is the total pressure, Cp, a is the heat capacity at constant pressure of dry air at
temperature T, Lv is the latent heat of evaporation at temperature T (see Latent Heat
of Evaporation), and Ma and Mv are the molar mass of dry air and water vapor,
respectively.
See also Saturation Pressure for the definition of the variable psat.
Figure 3-1: Relation between dry bulb, dew point, equivalent, and wet bulb temperatures.
SATURATION PRESSURE
The variable psat is defined by:
psat = fpsat ( T )
where T is the temperature. See Functions for the definition of the function fpsat.
RELATIVE HUMIDITY
The variable phi is defined by:
Lv = lv ( T )
where T is the temperature. See Functions for the definition of the function Lv.
pv = phi ⋅ psat
psat
csat = ------------
RT
HUMIDITY RATIO
The humidity ratio xvap is obtained from the vapor and air concentrations cv and ca,
the absolute pressure pA and temperature T by:
cv Mv
x vap = ----- ⋅ --------
ca Ma
where R is the universal gas constant, Mv is the water vapor molar mass, and Ma is the
dry air molar mass.
pA
c a = -------- – c v
RT
cv Mv
ω v = ----------------------------------
ca Ma + cv Mv
wc _ l = M _ v ⋅ cl
wc _ v = M _ v ⋅ cv
MOISTURE CONTENT
The moisture content wcVar is expressed in different ways depending on the domain
type:
• In moist air:
w ( φw ) = ρg ωv
Lv = lv ( T )
where T is the temperature. See Functions for the definition of the function Lv.
cflux = – δ p φ w ∇ ( p sat ( T ) )
In hygroscopic porous media and building materials domains, the first term of the
definition of dflux corresponds to liquid capillary flux. This term is numerically
handled like a diffusion flux, and is therefore included in the variable dflux.
tflux = cflux+dflux
ncflux_u = up ( cflux ) ⋅ un
ndflux_u = up ( dflux ) ⋅ un
ntflux_u = up ( tflux ) ⋅ un
Mass Balance
The mass balance variable, massBalance, is defined as
where:
• the variation of total moisture content per unit time in the domain, dwcInt, is:
d
dwcInt =
dt Ω w ( φ w ) dω
• the total net moisture rate, ntfluxInt, is the integral of the total moisture flux over
all external boundaries:
ntfluxInt =
∂Ω ext
tflux ⋅ n dσ
• the total moisture source, GInt, is the integral of the moisture sources over the
domains:
GInt =
Ω Gtot dω
T = T0 on ∂Ω
–n ⋅ q = q0 on ∂Ω
where
The inward heat flux, q0, is often a sum of contributions from different heat transfer
processes (for example, radiation and convection). The special case q0 = 0 is called
thermal insulation.
A common type of heat flux boundary conditions is one for which q0 = h·(Text − T),
where Text is the temperature far away from the modeled domain and the heat transfer
coefficient, h, represents all the physics occurring between the boundary and “far
away.” It can include almost anything, but the most common situation is that h
The CFD Module and the Heat Transfer Module contain a set of
correlations for convective heat flux and heating. See Heat Transfer and
Fluid Flow Coupling.
Several categories of boundary condition exist in heat transfer. Table 3-6 gives the
overriding rules for these groups.
A\B 1 2 3 4 5 6
1-Temperature X X X
2-Thermal Insulation X X
3-Heat Flux X X
4-Boundary heat source
5-Thin Layer X X
6-Isothermal Domain Interface X
• Locate the line that corresponds to the A group (see above the definition of the
groups). In the table above only the first member of the group is displayed.
• Locate the column that corresponds to the group of B.
• If the corresponding cell is empty, A and B contribute. If it contains an X, B
overrides A.
Example
Consider a boundary where Heat Flux is applied. Then a Symmetry boundary condition
is applied on the same boundary afterward.
In the heat transfer interfaces, the entire physics (equations and variables) are defined
in the spatial frame. When a moving mesh is detected, the user inputs for certain
features are defined in the material frame and are converted so that all the
corresponding variables contain the value in the spatial frame.
This subsection contains the list of all heat transfer nodes and the corresponding
definition frame:
Spatial: The inputs are entered by the user and defined in the spatial frame. No
conversion is done.
Material/(Spatial): For these physics nodes, select from the Material type list to decide
if the inputs are defined in the material or spatial frame. The default definition frame
is the material frame, which corresponds to Solid in the Material type list.
(Material)/Spatial: For these physics nodes, select from the Material type list to decide
if the inputs are defined in the material or spatial frame. The default definition frame
is the spatial frame, which corresponds to Nonsolid in the Material type list.
The volume reference temperature defines the density in the reference geometry that
should match with the geometry in the material frame. It is a model input of all the
features defined in the material frame with an input field for the density. Following
Table 3-8, Table 3-9, and Table 3-10, this concerns the following features: Solid,
Porous Medium, Biological Tissue, Building Material, Shape Memory Alloy, Thin
Layer (Heat Transfer Interface) and Solid (Heat Transfer in Shells Interface), Fracture
(Heat Transfer Interface) and Porous Medium (Heat Transfer in Shells Interface), and
Thin Rod.
In this section:
• Consistent Stabilization
• Inconsistent Stabilization
Consistent Stabilization
This section contains two consistent stabilization methods: streamline diffusion and
crosswind diffusion. These are consistent stabilization methods, which means that they
do not perturb the original transport equation.
STREAMLINE DIFFUSION
Streamline diffusion is active by default and should remain active for optimal
performance for heat transfer in fluids or other applications that include a convective
or translational term.
CROSSWIND DIFFUSION
Streamline diffusion introduces artificial diffusion in the streamline direction. This is
often enough to obtain a smooth numerical solution provided that the exact solution
of the heat equation does not contain any discontinuities. At sharp gradients, however,
undershoots and overshoots can occur in the numerical solution. Crosswind diffusion
addresses these spurious oscillations by adding diffusion orthogonal to the streamline
direction — that is, in the crosswind direction.
By default there is no isotropic diffusion. To add isotropic diffusion, select the Isotropic
diffusion check box. The field for the tuning parameter δid then becomes available. The
default value is 0.25; increase or decrease the value of δid to increase or decrease the
amount of isotropic diffusion.
• Stabilization Techniques
• Selection Information
All these options make it possible to build a coupling in different ways. Even if the use
of the predefined multiphysics coupling interfaces — Nonisothermal Flow and Conjugate
Heat Transfer — is the preferred choice, other alternatives can be of interest in
particular cases. This section describes the possibility for coupling heat transfer and
fluid flow interface and lists the advantages and limitations of each approach.
In this section:
See Domain, Boundary, Pair, and Point Nodes for Single-Phase Flow in the CFD
Module User’s Guide for a description of the nodes associated to these interfaces.
It is common to start a model with a single physics (for example, fluid flow), then
implement the second one (for example, heat transfer). Then adding the Nonisothermal
Flow multiphysics feature realizes the bidirectional coupling. It also redefines the
consistent stabilization so that the multiphysics coupling effects are accounted for in
the numerical stabilization.
With the Heat Transfer Module or the CFD Module, it accounts for the turbulence in
the coupling. In particular, it modifies the effective thermal conductivity and
implements thermal wall functions if the fluid flow model requires them. Those
modifications affect the implementation of several heat transfer features. It allows to
include work done by pressure changes and viscous dissipation, and Boussinesq
approximation is supported. Finally some physics features are updated when the
Nonisothermal Flow multiphysics feature is active. In particular, the Interior Fan and
Screen fluid-flow features are updated to account for the multiphysics coupling.
Note that the physics interface settings may not be optimal for the numerical treatment
of the coupling when the multiphysics feature is added afterward.
In addition, the heat transfer and fluid flow interfaces are set up with optimal interface
settings: the discretization order of the heat transfer interface is the same as the one
used for the fluid flow interface, and the pseudo time stepping is activated in both
interfaces.
Note that you can do a gradual implementation of the model: It is possible to start
from these multiphysics interfaces and to disable the multiphysics feature or one of the
physics in a first step and then reactivate them when the first step is validated.
The boundary temperature variable called ht.Tvar describes the wall temperature.
When the wall has a nonconstant temperature across its thickness, this variable contains
the average value between the temperatures of the two sides of the wall. The actual
definition of ht.Tvar depends on the model configuration.
The following list includes existing boundary temperature variables that are available
depending on the model configuration:
• T: general temperature variable that coincides with the wall temperature in most
cases
• TWall_u: upside wall temperature defined by a Wall or an Interior Wall feature with
turbulence only if a Fluid feature is defined on the upside of the wall.
• TWall_d: downside wall temperature defined by a Wall or an Interior Wall feature
with turbulence only if a Fluid feature is defined on the downside of the wall.
• Tu: temperature on the upside of the boundary.
• Td: temperature on the downside of the boundary.
• TExtFace: external temperature of an external boundary defined by a thermally
thick boundary condition.
• TuWF: temperature of the fluid in the turbulent boundary layer near the wall only if
the Fluid feature is defined on the upside of the Wall feature.
• TdWF: temperature of the fluid in the turbulent boundary layer near the wall only if
the Fluid feature is defined on the downside of the Wall feature.
Depending on the turbulence model selected for the flow, wall functions are used or
not:
The following sections summarize the definitions of the temperature variables for the
abovementioned configurations.
Fluid
Solid Fluid
Solid Fluid
Tu=Td=down(T) Tu=Td=TWall_d
Fluid Solid
Tu=Td=Twall_d
Td=TWall_d Tu=up(T)
TdWF=down(T)
Thermally thick layer boundary
Fluid Solid
Tu=Td=Twall_d Tu=Td=up(T)
Tu=Td=TWall_d Tu=Td=TWall_d
Td=TWall_d Tu=TWall_u
TdWF=down(T) TuWF=up(T)
In this section:
∂T
ρC p +∇⋅q = Q
∂t
q = – k ∇T
T = T0
–n ⋅ q = q0
In its basic form, the density, ρ, heat capacity, Cp, thermal conductivity, k, heat sources,
Q, constraint temperatures, T0, and heat fluxes, q0, are all constant, which leads to a
Different nonlinear solvers are also provided for these kinds of problems.
Linear Solver
For small number of degrees of freedom, the direct PARDISO solver is used. It is
known to be robust and fast for small-sized problems.
For larger models, the linear iterative GMRES solver with multigrid preconditioner is
used. In most cases, SOR is the presmoother and postsmoother. This solver is memory
effective and fast for large models.
When the model contains settings that lead to a system matrix with 0 on the diagonal
(for example, Lagrange multipliers for weak constraints), SOR cannot be used and is
replaced by Vanka, which is usually slower and uses more memory.
Several options are available to tune the linear solver settings. This paragraph focuses
only on the most commonly used ones.
• When the convergence graph of GMRES shows a slow down every 50 iterations, the
Number of iteration before restart parameter (default value of 50) should be
increased — doubled for example. This may also increase the memory consumption.
• Increasing the Number of iteration in the Multigrid settings, and in the presmoother
and postsmoother nodes improves the quality of the preconditioner and
convergence of GMRES.
• Since an excessive difference between two multigrid levels can affect the
convergence, lowering the Mesh coarsening factor in the Multigrid settings can help
convergence.
• Consider creating the multigrid level meshes manually if the automatic coarsening
method fails or leads to poor quality meshes.
Nonlinear Solver
• Another physics interface is solved together with heat transfer. The dependent
variables of the Heat Transfer interface are placed in a separate segregated group.
• Radiation in participating media using the Discrete ordinates method defines a large
number of dependent variables (up to 168 for a single wavelength), which are placed
in segregated groups. The number of dependent variables per segregated group and
the nonlinear method settings depend on the Performance index parameter available
in the heat transfer interface settings in the Participating Media Settings section.
When multiple wavelength are considered, the variables relative to distinct
wavelengths are not mixed together in the segregated groups.
• The Thermal Damage subfeature (added under Biological Tissue feature) defines an
additional variable alpha that is placed in a dedicated segregated group.
A Newton nonlinear method with a constant damping factor is set by default, with the
following settings:
• Using the Automatic highly nonlinear (Newton) option forces to start the
computation with a very low damping factor and increases it carefully. Alternatively
a low constant damping factor can be used. The damping factor ranges between 0
and 1. A constant damping factor equal to 0.1 is a very low value and should be
robust but slow to converge. For low values of the damping factor, it is thus usually
needed to increase the number of nonlinear iterations. If the nonlinear solver is
unstable with such a damping factor then the automatic option should be used
• Using a constant damping factor equal to 1 for linear problems. The linearity is
determined at the beginning of the resolution and indicated in the Log section of the
solver window.
• Providing a good initial value is an asset for computational speed.
• In the convergence area, the fully coupled solver has a better convergence rate than
the segregated solver.
• Using minimal Jacobian update option avoid to spend time in Jacobian
computation. This is suited for linear models and models with mild nonlinearities.
When a Thermal Damage subfeature is present under Biological Tissue feature, particular
settings for the time-dependent solver are used to efficiently compute the damage
indicators:
• The Absolute Tolerance of the scaled damage indicator variable is set to 1, meaning
that these variable are neglected in the error estimate.
• The damaged tissue indicator, α, is solved with an iterative Jacobi method.
• If the Adaptive mesh refinement option is selected in the study settings, the error
indicator is set to ∇θ d, sm ⋅ ∇θ d, sm , where θd, sm is the smoothed indicator of
necrotic tissue (the fraction of necrotic tissue, θd, is discontinuous in general).
• If the Temperature threshold option is used in the Biological Tissue feature, the instant
necrosis indicator, alphanecr, is placed in the Previous Solution step. This setting
avoids wrong detection of irreversible damage due to nonlinear iterations that may
go through a state where the damage criteria is met and then converge to a solution
where the damage criteria is no longer met. It uses a direct linear solver. The default
nonlinear method is the Newton method with constant damping factor.
When the Irreversible Transformation feature is active, similar settings are used:
The default solver settings for transient heat transfer and moisture transport define the
maximal number of nonlinear iterations to 5. If this is not sufficient, it is recommended
• An implicit way is to define a lower relative tolerance in the study settings. When the
relative tolerance is lowered, the absolute tolerance should be reduced in the same
proportion. The default Relative tolerance is set to 0.01 for Time-Dependent studies,
and to 0.1 for Stationary studies.
• The most explicit way is to define a maximum time step. This is an appropriate
option when the same maximum time step is relevant for the entire simulation.
Otherwise, it is possible to include times of interest in the Times field of the
time-dependent study and to use the Intermediate option in the Time Stepping
settings.
• Lastly you can control the time step by triggering an event when a particular
condition is meet (see the documentation about The Events Interface in the
COMSOL Multiphysics Reference Manual). This advanced method can be
efficient when the other simpler methods are not applicable.
It is also recommended to inspect the solver log and check the default scaling of
dependent variables in case of convergence failure. In case of incorrect automatic
scaling, consider using Manual settings in the Dependent Variable attribute node.
• Surface-to-surface radiation makes the Jacobian matrix of the discrete model partly
filled as opposed to the usual sparse matrix. The additional nonzero elements in the
matrix appear in the rows and columns corresponding to the radiosity degrees of
freedom. It is therefore common practice to keep the element order of the radiosity
variable, J, low. By default, linear Lagrange elements are used irrespective of the
shape-function order specified for the temperature. When you need to increase the
MULTIPHYSICS MODELS
Unless the model contains a multiphysics node that defines a coupling between a Heat
Transfer interface and another interface (see Multiphysics Couplings below), each
physics interface defines default solver settings that are merged.
The Heat Transfer interfaces always define a dedicated segregated group that uses a
linear solver optimized for the heat transfer equations. For strongly coupled models, it
may be efficient to merge two (or more) segregated steps. In this case, a unique linear
solver must be chosen for the fully coupled solver or the new segregated group.
Time-dependent settings from different physics interfaces may compete. When the
different settings are merged the strictest one is kept.
MULTIPHYSICS COUPLINGS
When a Heat Transfer interface is coupled with another physics interface through a
multiphysics coupling feature, additional predefined default settings are loaded. The
next two paragraphs describes the subtleties of the Nonisothermal Flow,
Electromagnetic Heating, Heat and Moisture Transport, and Moisture Flow
interfaces.
Nonisothermal Flow
The Nonisothermal Flow multiphysics coupling controls the solver settings for the
flow and the temperature dependent variables.
When it assumes a weak coupling between the flow and the heat interfaces (typically
no gravity force accounted for in the flow interface), the default solver contains
When a strong coupling is assumed (gravity force accounted for in the flow interface),
the default solver merges the temperature, pressure, and velocity groups. In this case,
the solver corresponds to the default solver of the heat transfer interface. This is meant
to suits well for nonisothermal flows in which natural convection dominates. See
Default Linear Settings for Heat Transfer and Moisture Transport interfaces, Default
Nonlinear Settings for Heat Transfer Interfaces, and Default Time Settings for Heat
Transfer Interfaces for details.
Electromagnetic Heating
The Electromagnetic Heating multiphysics interfaces (Joule heating, Laser Heating,
Induction Heating, and Microwave Heating) define default settings that solve the
temperature and the electromagnetic fields using a coupled step. It can be the fully
coupled nonlinear solver if there is no additional variable to solve for, otherwise it is a
segregated step containing the temperature and the electromagnetic variables.
However when radiation in participating media or damage variable are solved they are
placed in a separate group as described above.
The default solver merges the temperature and relative humidity groups. The solver
corresponds to the default solver of the moisture transport interface. See Default
Linear Settings for Heat Transfer and Moisture Transport interfaces, Default
Nonlinear Settings for Moisture Transport Interfaces, and Default Time Settings for
Moisture Transport Interfaces for details.
Moisture Flow
The Moisture Flow multiphysics coupling controls the solver settings for the flow
(velocity and pressure) and the relative humidity variables.
When a strong coupling is assumed (gravity force accounted for in the flow interface),
the default solver merges the relative humidity, pressure, and velocity groups. In this
case, the solver corresponds to the default solver of the moisture transport interface.
See Default Linear Settings for Heat Transfer and Moisture Transport interfaces,
Default Nonlinear Settings for Moisture Transport Interfaces, and Default Time
Settings for Moisture Transport Interfaces for details.
Therefore, dedicated tools are available to plot and evaluate the results in the 1D extra
dimension.
• The Layered Material Slice plot displays the variables along the layered material, for
specified through-thickness locations in the 1D extra dimension.
• The Through Thickness plot displays the variables in the thickness direction, at
specified points on the boundary
• xdimTag is the extra dimension tag. For example, it may be slmat1_xdim, when a
single layer material (slmat1) is selected in the Layered material list of the Thin Layer
node. This tag can be deduced from the Selection column of the Equation View
subnode of the node applied on the thin layer, by clicking the Show More Options
button ( ) on the Model Builder tool bar and selecting Equation View.
• xd is the coordinate in the extra dimension. It varies from 0 to slmat1.th, which
is the total thickness of the layered material. By convention, xd=0 corresponds to
the downside of the boundary where the thin layer is defined, whereas
xd=slmat1.th corresponds to its upside. Upside and downside settings can be
visualized by plotting the global normal vector (nx, ny, nz), that always points from
downside to upside. See Tangent and Normal Variables in the COMSOL
Multiphysics Reference Manual. Note that the normal vector (ht.nx, ht.ny,
ht.nz) may be oriented differently.
• expr is the quantity to be evaluated at the point xd. For example, it can be set to T
to evaluate the temperature. There are others postprocessing variables defined on
the extra dimension that can be found in the Equation View subnode of the node
applied on the thin layer.
Extra
dimension
Computational
domain Along
the layer
xd
Figure 3-2: Schematic representation of a 2D geometry with a thin layer composed of three
layers, with an evaluation of the results along the layer at the coordinate xd.
• compTag is the component tag. In most cases, this tag is comp1. It is possible to
check it in the Properties window of the component node (display it by
right-clicking on the node and selecting Properties).
• x0 and y0 are the coordinates of the point in the base geometry that belongs to the
boundary linked with the extra dimension. Note that these are 2D coordinates from
the global coordinate system and not curvilinear coordinates.
• expr is the quantity to be evaluated at the point (x0, y0). For example, it can be set
to T to evaluate the temperature.
( x0, y0)
Figure 3-3: Schematic representation of a 2D geometry with a thin layer composed of three
layers, with an evaluation of the results through the layer at the point (x0,y0).
This operator can be used to define manually a through thickness plot in the layered
material. For all dimensions, the section is represented in a line graph under a 1D plot
group. In order to use this, the Dataset selected in the Data section of the 1D plot
group has to select the extra dimension as component. One method is to duplicate the
default Solution node under the Datasets node, set the Component of the new node to
the extra dimension, and use this new dataset into the 1D plot group, with the domains
of the extra dimension selected.
• featTag is the tag of the physics node defining the operator xdintopall. For
example, for a Thin Layer node in the Heat Transfer interface, it may be ht.sls1.
• expr is the quantity to be integrated. For example, it can be set to T to evaluate the
temperature.
This chapter describes how to evaluate the bulk temperature in 3D and 2D models,
using the built-in geometry variables and operators. The nonisothermal flow is
supposed to take place along a tube, with any orientation and curvature, even though
the bulk temperature may be defined in other geometries. It contains the following
sections:
( u ⋅ τ )T ds
S
--------------------------------
⊥
Tb =
u ⋅ τ ds
S⊥
where S ⊥ is a cross-sectional surface of the tube where the flow takes place, u is the
flow velocity field, τ is the tangential direction of the flow, and T is the temperature
field.
The bulk temperature can be used to evaluate the heat transfer coefficient h at tube
walls:
qw
h = --------------------
Tw – Tb
See ball, circle, disk, and sphere in the COMSOL Multiphysics Reference
Manual for details about the operators.
In the evaluation of the bulk temperature, the radius R should be chosen large enough
to encompass the tube cross section, but not too large regarding numerical precision.
If the evaluation disk overlaps other domains, the evaluation may be restricted to the
domain of interest by adding a test on the dom variable.
The evaluation of the tangential velocity over the cross section requires to access the
vector tangent to the edge, in the computational domain. However, the tangent
vectors t1 and t2 are defined on boundaries only. By applying the general extrusion
operator genext1 on the edge, with a mapping to the domain, the tangential vector
genext1(t1) is defined all over the cross section, and can be used in the diskint
operator expression.
Tb =
diskint(R,if(dom==i,ut*T,0),tol)/diskint(R,if(dom==i,ut,0),tol)
In order to integrate over the cross section of the flow, the cross section line is
parametrized as follows:
S ⊥ = { ( x + sn x, y + sn y ), s ∈ [ 0, 1 ] }
where (nx,ny) is the normal vector to the boundary, and hence the vector tangent to
the cross section.
The at2 operator, which is the spatial evaluation operator in 2D, should be used to
evaluate the expression expr on this parametrization of the cross section.
In practice, the integral is performed over a larger interval [-R,R], where R should be
chosen large enough to encompass the cross section of the tube, but not too large
regarding numerical precision. The isnan operator can be used to ensure that the
solution is defined at the prescribed coordinates, and as in 3D models, the dom variable
may be used to identify the domain of interest.
See isinf and isnan in the COMSOL Multiphysics Reference Manual for
details about the operators.
As in 3D models, the evaluation of the tangential velocity over the cross section
requires to access the vector tangent to the boundary, in the computational domain.
However, the tangent vector t is defined on boundaries only. By applying the general
extrusion operator genext1 on the boundary, with a mapping to the domain, the
Tb = A/B
where
A = integrate(if(isnan(at2(x+s*nx,y+s*ny,ut*T)),0,at2(x+s*nx,y+s*
ny,ut*T)),s,-R,R,tol)
B = integrate(if(isnan(at2(x+s*nx,y+s*ny,ut*T)),0,at2(x+s*nx,y+s*
ny,ut)),s,-R,R,tol)
In this section:
When no special mention is added, the term temperature stands for the
dry bulb temperature. See Moist Air Variables for the definition of the dry
bulb temperature, dew point temperature, and relative humidity.
In the Ambient Properties node ( ), you can define ambient variables to be available
as inputs from several features: the temperature Tamb, the absolute pressure pamb, the
Three options are available for the definition of the Ambient data:
• When User defined (the default) is selected, the Temperature Tamb, the Absolute
pressure pamb, the Relative humidity φamb, the Wind velocity vamb, the Precipitation
rate P0,amb, the Clear sky noon beam normal irradiance Isn,amb, and the Clear sky
noon diffuse horizontal irradiance Ish,amb should be specified directly.
• When either Meteorological data (ASHRAE 2013) or Meteorological data (ASHRAE 2017)
is selected, the ambient variables are computed from monthly and hourly averaged
measurements, made over several years at weather stations worldwide. See
Processing of ASHRAE Data for more information. Further settings for the choice
of the location, time, and ambient conditions are needed; and additional input fields
are displayed underneath.
Location
In this section you can set the location by choosing among more than 8000 weather
stations worldwide. Two options are available for the selection of the Weather station:
• When From list is selected, click on the Set Weather Station button to open the
Weather Station browser that allows you to select a Region, a Country, and a Station.
A single country may be available for more than one region selection if it
has stations spread over different regions. For example, United States of
America is available in the Country list when either North America,
Eurasia, or Oceania is selected in the Region list.
Time
The Date and Local time should be set by entering values or expressions in the Day,
Month, Hour, Minute, and Second fields of the two tables.
If On is selected in the Specify year list, a value or expression for the Year should also be
set. As the data are given as averages over several past years, this input is only used for
the detection of leap years, in order to interpolate the data over the months.
For temporal studies, these inputs define the start time of the simulation. By default,
the Update time from solver check box is selected, and the time is then automatically
updated with the time from the solver to evaluate the variables by interpolation of the
measured data. Clear this check box to manually set the time update.
See Processing of ASHRAE Data for more information about the data.
Ambient Conditions
Based on the measured data, several conditions are available for the Temperature, the
Dew point temperature, the Wind speed, and the Precipitation rate. The formula for each
condition is recalled in Table 3-11, Table 3-12, and Table 3-13. The Average
conditions correspond to weighted means of the measured data, whereas the other
conditions are obtained by applying standard or modified deviations (Low, High, and
User defined coefficient for deviation conditions), user defined corrections, or wind
correlations to the average conditions; or by taking the minimum or maximum of the
measured data (Lowest and Highest conditions). More information about these
definitions can be found in Ambient Variables and Conditions.
TABLE 3-11: TEMPERATURE CONDITIONS
CONDITION DEFINITION
CONDITION DEFINITION
CONDITION DEFINITION
CONDITION DEFINITION
The conditions set for Temperature and Dew point temperature should be
consistent in order to keep the temperature larger than the dew point
temperature. However, all settings combinations are available, and the
relative humidity is majored by 1 when necessary.
The sum of the Clear sky noon beam normal irradiance and the Clear sky
noon diffuse horizontal irradiance is available through the postprocessing
variable ht.Is_amb, defined as the Ambient solar irradiance.
Figure 3-4: Computation of weighted mean from frequencies of observations for the
diurnal fluctuations of temperature.
These values are used for the definition of different conditions, as detailed in Ambient
Variables and Conditions.
All the monthly averaged observations except the solar irradiance are supposed to be
made at the middle of each month. This time depends on the number of days in the
month:
• Months with 31 days (January, March, May, July, August, October, December):
data at the 16th at noon
• Months with 30 days (April, June, September, November): data at the 16th at
midnight
• Months with 29 days (February, leap years): data at the 15th at noon
• Months with 28 days (February, other years): data at the 15th at midnight
Finally, the solar irradiance observations are made at the 21st of each month at noon.
Depending on the number of days in the month, this date corresponds to 68% (for
months with 31 days), 70% (for months with 30 days), or 75% (for February) of the
month. The leap years are not considered and the 21st of February always corresponds
to 75% of this month.
• The annual fluctuation of the dew point temperature, the relative humidity, the
wind speed, and the direct and diffuse solar irradiances.
• The annual and diurnal fluctuation of the temperature.
In all cases, the interpolation is of second order, with continuous first-order derivative.
CONDITIONS OF TEMPERATURE
• Average:
T amb = T station
• Low:
T amb = T station – σ T, station
• High:
T amb = T station + σ T, station
• Lowest:
T amb = min ( T station )
• Highest:
• <Tstation> (SI unit: K) is the weighted mean of the observed values of temperature
at the station.
• σT,station (SI unit: K) is the standard deviation of the observed values of
temperature at the station.
• Tstation (SI unit: K) is the set of the observed values of temperature at the station.
• cσ (dimensionless) is a user-defined multiplicative coefficient applied to σT,station.
• ΔT (SI unit: K) is a user-defined additive correction applied to <Tstation>.
All these conditions are illustrated on Figure 3-5 for the variation of temperature over
1 day at New York/John F. Ke, on the 1st of June.
Figure 3-5: Comparison of ambient conditions for the temperature at New York/John F.
Ke, on the 1st of June, with ASHRAE Weather Data Viewer 5.0 (©2013 ASHRAE,
www.ashrae.org. Used with permission.)
Additional conditions are defined from observed couples of temperature and wind
speed and direction values:
1
ΔT wind = --- max ( ΔT ws, station, ΔT wd, station )
2
where ΔTws,station (SI unit: K) and ΔTwd,station (SI unit: K) are respectively the
maximal variations of observed values of temperature correlated with a set of wind
speed and direction observed values.
The heating and cooling wind correlations are illustrated on Figure 3-6 for the
variation of temperature over 1 day, at New York/John F. Ke, on the 1st of June.
Figure 3-6: Comparison of heating and cooling wind correlations for the temperature at
New York/John F. Ke, on the 1st of June, with ASHRAE Weather Data Viewer 5.0
(©2013 ASHRAE, www.ashrae.org. Used with permission.)
• Low:
DPT amb = DPT station – σ DPT, station
• High:
• Lowest:
DPT amb = min ( DPT station )
• Highest:
where:
• <DPTstation> (SI unit: K) is the weighted mean of the observed values of dew point
temperature at the station.
• σDPT,station (SI unit: K) is the standard deviation of the observed values of dew
point temperature at the station.
• DPTstation (SI unit: K) is the set of the observed values of dew point temperature
at the station.
All these conditions are illustrated on Figure 3-7 for the variation of the dew point
temperature over 1 year at New York/John F. Ke.
Figure 3-7: Comparison of the ambient conditions for the dew point temperature at New
York/John F. Ke, with ASHRAE Weather Data Viewer 5.0 (©2013 ASHRAE,
www.ashrae.org. Used with permission.)
• Low:
v amb = v station – σ v, station
• High:
v amb = v station + σ v, station
• Lowest:
v amb = min ( v station )
• Highest:
where:
• <vstation> (SI unit: m/s) is the weighted mean of the observed values of wind
velocity at the station.
• σv,station (SI unit: m/s) is the standard deviation of the observed values of wind
velocity at the station.
• vstation (SI unit: m/s) is the set of the observed values of wind velocity at the
station.
Figure 3-8: Comparison of the ambient conditions for the wind speed at New York/John
F. Ke, with ASHRAE Weather Data Viewer 5.0 (©2013 ASHRAE, www.ashrae.org.
Used with permission.)
• Low:
P 0, amb = P 0, station – σ P
0, station
• High:
P 0, amb = P 0, station + σ P
0, station
• Lowest:
P 0, amb = min ( P 0, station )
• Highest:
• <P0,station> (SI unit: m/s) is the weighted mean of the observed values of
precipitation rate at the station.
• σP0,station (SI unit: m/s) is the standard deviation of the observed values of
precipitation rate at the station.
• P0,station (SI unit: m/s) is the set of the observed values of precipitation rate at the
station.
PRESSURE
p amb = p station
where pstation (SI unit: Pa) is the observed value of absolute pressure at the station.
Only a single value is available, so this data does not vary with time.
RELATIVE HUMIDITY
The relative humidity φ amb (dimensionless) is computed from the temperature Tamb
and the dew point temperature DPTamb with the following relation:
Figure 3-9: Diurnal fluctuations of relative humidity for different ambient conditions at
New York/John F. Ke, on the 1st of June, with ASHRAE Weather Data Viewer 5.0
(©2013 ASHRAE, www.ashrae.org. Used with permission.)
MOISTURE CONTENT
The moisture content xvap,amb (dimensionless) is computed from the temperature
Tamb, the absolute pressure pamb, and the relative humidity φ amb with the following
relation:
where psat(Tamb) is the saturation pressure of vapor at Tamb, and Mv and Ma are the
molar masses of water vapor and dry air.
where Isn,station (SI unit: W/m3) and Ish,station (SI unit: W/m3) are respectively the
observed values of the clear sky noon beam normal and horizontal diffuse solar
irradiances.
Figure 3-10 illustrates the evolution of ambient solar irradiance for New York/John F.
Ke, over the year.
Figure 3-10: Decomposition of solar irradiance into normal and horizontal irradiance at
New York/John F. Ke, with ASHRAE Weather Data Viewer 5.0 (©2013 ASHRAE,
www.ashrae.org. Used with permission.)
As an example, consider a plate of concrete with a cold bottom wall at temperature T0,
placed in a hot environment with an air flux at temperature Tamb:
Figure 3-12: Integral of convective and radiative heat fluxes (W/m) along the top
boundary, for two values of emissivity, as a function of the temperature difference.
Whereas the boundary radiative and convective heat fluxes are of the same order for
temperature difference up to 500K when the emissivity is low (ε=0.1), radiation
becomes the dominant mode of heat transfer even for a small temperature difference
when the emissivity is high (ε=0.9). Note that the convective heat flux decrease
observed for high temperature gradients is related to the fact that the velocity, not the
mass flow rate, is prescribed at the air inlet.
See Heat Transfer Variables for the definition of the variables q0 and
rflux giving access to the convective and radiative heat fluxes on
boundaries.
Two configurations are considered regarding the functionalities available for the
modeling of radiative heat transfer:
SURFACE-TO-AMBIENT SURFACE-TO-SURFACE
Description
SURFACE-TO-AMBIENT SURFACE-TO-SURFACE
Description
In this section:
This is of particular interest when the geometry contains inlets that are fed by channels
that are not represented in the geometry.
– n ⋅ q = ρΔHu ⋅ n
where:
T
ΔH T = Tustr Cp dT (3-2)
and
p
1
ΔH p = pustr --ρ- ( 1 – αp T ) dp (3-3)
Equation 3-1 expresses the fact that the normal conductive heat flux at the inflow
boundary is proportional to the flow rate and enthalpy variation between the upstream
conditions and inlet conditions.
There is another classical case where this term cancels out: when the fluid is modeled
as an ideal gas. Indeed, in this case,
1
α p = ----
T
When the pressure contribution to the enthalpy is neglected, the boundary condition
reads:
T
k∇T ⋅ n = ρ
Tustr Cp dT u ⋅ n (3-4)
When advective heat transfer dominates at the inlet (large flow rates), the temperature
gradient, and hence the heat transfer by conduction, in the normal direction to the
inlet boundary is very small. So in this case, Equation 3-4 imposes that the enthalpy
variation is close to zero. As Cp is positive, the Inflow boundary condition requires
T=Tustr to be fulfilled. So, when advective heat transfer dominates at the inlet, the
Inflow boundary condition is almost equivalent to a Dirichlet boundary condition that
prescribes the upstream temperature at the inlet.
Conversely, when the flow rate is low or in the presence of large heat sources or sinks
next to the inlet, the conductive heat flux cannot be neglected. In addition, the inlet
temperature has to be adjusted to balance the energy brought by the flow at the inlet
and the energy transferred by conduction from the interior, as described by
Equation 3-4. This makes it possible to observe a realistic upstream feedback due to
thermal conduction from the inlet surroundings.
In this section:
The heat sinks are composed of a rectangular base and an array of pin or straight fins,
as shown on Figure 3-14.
border offset
base
All entities are fully parameterized, making them easy to use as parts in industrial
models where heat sinks are used for cooling. For example, you can control the
number, the shape, the dimensions, and the placement of the fins on the base. In
addition, fillet, chamfer, and notch transformations can be applied to the fins, and
parameter checks are applied to ensure that the values set in the Input Parameters
section are valid. Finally, the fins can be defined as solids or as boundaries for
computational efficiency.
In the Heat Sink — Pin Fins part, all the fins are pins with the same dimension, whereas
the outer and inner fins (in the y direction) can have distinct dimensions in the Heat
Sink — Dissimilar Border Pins part. You may use the Heat Sink — Straight Fins part to
define a heat sink made of only straight fins.
By default, the base of the heat sink is positioned at the origin of the xy-plane. You can
apply a Displacement and a Rotation to this configuration in the Position and Orientation
of Output section.
n_fins_x=3
n_fins_y=4 Y_fins_top_2
Y_fins_top
X_fins_top
Z_fins
Z_base o_y
o_x
Y_fins_bottom_2
Y_fins_bottom X_fins_bottom
Y_base X_base
Figure 3-15: Fins and base parameters in Heat Sink — Parameterized part
TABLE 3-18: DEFAULT TETRAHEDRAL MESH WITH 3D FINS AND SHELL FINS
SHELL=0 SHELL=1
STEP
A step can be defined in the x direction at the center on the base. The parameter
step_width specifies the number of filled gaps from middle to border by the step, as
shown on Figure 3-16. This option is not available for shell fins (shell=1).
step_height step_height
step_width=1 step_width=1
FILLET
Finally, a fillet transformation can be applied at the top and bottom of the fins, as
shown on Figure 3-17. The fillet transformation at the bottom of the fins is not
available for shell fins (shell=1).
fillet_top=1
fillet_bottom=1
In addition, notch and chamfer transformations can be applied to the fins, as shown
on Figure 3-18.
notch=1 chamfer=1
notch_height chamfer_height
notch_width chamfer_width
Figure 3-18: Notch and chamfer parameters in Heat Sink — Straight Fins part
REFERENCES | 161
162 | CHAPTER 3: MODELING WITH THE HEAT TRANSFER MODULE
4
This chapter details the theory of the physics interfaces, multiphysics couplings,
and features found under the Heat Transfer branch ( ).
In this chapter:
163
• Theory for Radiation in Participating Media
• Theory for Moisture Transport
• Theory for the Heat Transfer Multiphysics Couplings
• Theory for Thermal Contact
• Out-of-Plane Heat Transfer
• Convective Heat Transfer Correlations
• Nucleate Pool Boiling Correlation
• Equivalent Thermal Conductivity Correlations
• Temperature Dependence of Surface Tension
• Heat Flux and Heat Balance
• Frames for the Heat Transfer Equations
• References
In this section:
The internal energy, EΩ (SI unit: J), is an extensive state function of these three
variables. It measures the amount of energy in the system excluding kinetic energy and
potential energy from external applied forces and is the subject of conservation laws
more detailed in The Heat Balance Equation section. To fit with the finite element
method solved by COMSOL Multiphysics, specific quantities per unit mass are
preferred:
SΩ VΩ
S = --------- ν = ---------
MΩ MΩ
The specific internal energy, E (SI unit: J/kg), is then a function of specific entropy,
S, and specific volume, ν, related to EΩ by:
1
E ( S, ν ) = --------- E Ω ( S Ω, V Ω, M Ω )
MΩ
For a solid, the specific internal energy, E(S, F), is a function of entropy and
deformation gradient, F.
Internal energy is related to the enthalpy, H, via the following for a fluid:
p
H = E + ---
ρ
1
H = E – --------------------- P:F
det(F)ρ
Compared to the internal energy, the enthalpy also includes the pressure-volume
potential energy, p ⁄ ρ, necessary for instance in volume expansion after an isobaric
transformation.
FIRST-ORDER PARAMETERS
The variations of internal energy correspond to variations of entropy and volume
according to:
∂E ∂E
dE = dS + dν
∂S ν ∂ν S
∂E ∂E
T = p = – (4-1)
∂ S ν ∂ν S
dE = T dS – p dν
∂E ∂E
T = P = det(F)ρ (4-2)
∂S F ∂F S
1
dE = T dS + --------------------- P: dF
det(F)ρ
Here, the counterpart of the fluid pressure is the first Piola-Kirchhoff stress tensor, P.
SECOND-ORDER PARAMETERS
Second-order parameters correspond to second partial derivatives of the specific
internal energy and provide a various number of thermodynamic coefficients. These
are usually given as material properties of the domain material. Among them, the heat
capacity at constant pressure and the coefficient of thermal expansion are most often
provided. For a fluid, these are
T 1
C p = --------------- α p = ------------------- (4-3)
∂T ∂T
ν
∂ S ν ∂ν S
–1
– 1 ∂T
α = F
T
C p = ---------------- (4-4)
∂T ∂F S
∂ S F
The heat capacity at constant pressure and coefficient of thermal expansion are related
to the enthalpy, seen as a function of T and p (or P), according to:
∂H = C ∂H = ν ( 1 – α T )
∂T p p ∂p T p
SENSIBLE ENTHALPY
The enthalpy can then be retrieved from Cp and αp (or α) by:
r1
H = H ref + r ∇r H ( r ) ⋅ dr
0
(4-5)
P
11
P
22
P
33
r =
p
or r = P
T 12
P
23
P
13
T
The starting point, r0, is the value of r at reference conditions, that is, pref (1 atm) and
Tref (293.15 K) for a fluid. The ending point, r1, is the solution returned after
simulation. In theory any value can be assigned to the enthalpy at reference conditions,
Href (Ref. 2), and COMSOL Multiphysics sets it to 0 J/kg by default. The integral in
1
α p = ----
T
therefore
∂H = 0
∂p T
∂H = 0
∂p T
∂H = ( C ∂ω v
p, v – C p, a ) ---------
-
∂T p ∂c v
and Equation 4-6 should be modified to include an integral over the vapor
concentration:
cv
T ∂ω v ∂H
H = H ref + Tref ( C p, v – C p, a ) ---------- ⋅ dT +
∂c v cref ∂ cv ⋅ dcv (4-7)
k xx k xy k xz
k = k xy k yy k yz
k xz k yz k zz
INTEGRAL FORM
The first law of thermodynamics states that the variations of macroscopic kinetic
energy, KΩ, and internal energy, EΩ, of a domain Ω are caused either by the mechanical
power of forces applied to the system, Pext, or by exchanged heat rate, Qexch (2.3.53
in Ref. 4):
dE Ω dK Ω
----------- + ----------- = P ext + Q exch (4-9)
dt dt
Mass and momentum balance are needed to complete the description of the system.
The mechanical laws, either for solids or fluids, generate the following balance
equation between variation of kinetic energy, KΩ, stress power, Pstr, and power of
applied forces, Pext (2.3.64 in Ref. 4):
dK Ω
----------- + P str = P ext (4-10)
dt
This equation involves quantities of the macroscopic level where the variation of the
kinetic energy due to some forces applied to it reflects a sensible displacement. In
COMSOL Multiphysics, the Solid Mechanics or Single-Phase Flow interfaces are
examples of physics interfaces that simulate the macroscopic level described by
Equation 4-10.
Combining Equation 4-9 and Equation 4-10 yields the so-called heat balance
equation (2.3.65 in Ref. 4):
dE Ω
----------- = P str + Q exch (4-11)
dt
This time, the equation involves quantities of the microscopic level (exchanged heat
rate, Qexch, and internal energy, EΩ) more concerned with the atomic vibrations and
similar microscopic phenomena that are felt as heat. The presence of the stress power,
Pstr, in both Equation 4-10 and Equation 4-11 stands for the fact that such power is
LOCALIZED FORM
In this paragraph, the different terms of Equation 4-11 are more detailed to obtain the
localized form of the heat balance equation.
EΩ = Ω E dm
Note that by conservation of mass, the variation of internal energy in time is:
dE Ω dE dE
----------- =
dt Ω d t dm = Ω ρ d t dv
In these last relations, ρ is the density, and dv denotes an elementary volume of Ω.
Contrary to the constant elementary mass, dm, the elementary volume changes by
expansion or contraction of the domain. Recall that the derivation operator d ⁄ dt
under the integrals is in the material frame (see Time Derivative in the Frames for the
Heat Transfer Equations section).
Stress Power
The stress power, derived from the Continuum Mechanics theory, is defined by
(2.3.59 in Ref. 4):
P str = Ω ( σ:D ) dv
where σ is the Cauchy stress tensor and D is the strain rate tensor. The operation “:”
is a contraction and can in this case be written on the following form:
Note that in fluid mechanics, the Cauchy stress tensor is divided into a static part for
the pressure, p, and a symmetric deviatoric part, τ, as in:
so that Pstr becomes the following sum of pressure-volume work and viscous
dissipation:
P str = Ω p ( ∇ ⋅ u ) dv – Ω ( τ: ∇u ) dv
Equivalently, the stress power can also be expressed as:
----------------- P: dv
1 dF
P str = Ω det (F) d t
Exchanged Heat
Finally, the exchanged heat rates, Qexch, account for thermal conduction (see Fourier’s
Law at Equation 4-8), radiation and potentially additional heat sources. Joule heating
and exothermic chemical reactions are such examples of domain heat source. The
different kinds of exchanged heat are summarized by the equality below:
dE
Ω ρ d t dv + ∂Ω ( q ⋅ n ) ds + ∂Ω ( qr ⋅ n ) ds = Ω ( σ:D ) dv + Ω Q dv (4-13)
dE
ρ + ∇ ⋅ ( q + q r ) = σ:D + Q (4-14)
dt
∂E
ρ + ρu ⋅ ∇E + ∇ ⋅ ( q + q r ) = σ:D + Q (4-15)
∂t
This verbally means that variations of internal energy in time are balanced by
convection of internal energy, thermal conduction, radiation, dissipation of mechanical
See Frames for the Heat Transfer Equations for more details about the use
of material and spatial frames in the Heat Transfer interfaces.
• Conservation of mass
• Conservation of linear momentum
• Conservation of angular momentum
The equations corresponding to each of them are recalled below in Table 4-1. For
more details about the theory of Solid and Fluid Mechanics, see the Structural
Mechanics Module User’s Guide and CFD Module User’s Guide.
TABLE 4-1: CONSERVATION OF MASS AND MOMENTUM
Conservation of ρ 0 = ρdet(F) ∂ρ
Mass + ∇ ⋅ ( ρu ) = 0
∂t
Conservation of du ∂u
Linear Momentum ρ = ∇ ⋅ σ + Fv ρ + ρ ( u ⋅ ∇ )u = ∇ ⋅ σ + F v
dt ∂t
Conservation of T T
σ = σ σ = σ
Angular Momentum
When modeling a heat transfer problem with one of the Heat Transfer interfaces, the
aforementioned laws needs to be respected. For example, the velocity field, u,
provided in the energy equation and responsible for convection in a fluid, should
satisfy the continuity equation here below in order to avoid unphysical results.
∂ρ
+ ∇ ⋅ ( ρu ) = 0
∂t
∂T
ρC p ------- + u trans ⋅ ∇ T + ∇ ⋅ ( q + q r ) = – αT: + Q
dS
(4-16)
∂t dt
For a steady-state problem the temperature does not change with time and the terms
with time derivatives disappear.
The first term on the right-hand side of Equation 4-16 is the thermoelastic damping
and accounts for thermoelastic effects in solids:
dS
Q ted = – αT: (4-17)
dt
It should be noted that the d ⁄ dt operator is the material derivative, as described in the
Time Derivative subsection of Material and Spatial Frames.
∂T ∂p
ρC p ------- + u ⋅ ∇ T + ∇ ⋅ ( q + q r ) = α p T + u ⋅ ∇p + τ: ∇u + Q (4-18)
∂t ∂t
• the Cauchy stress tensor, σ, is split into static and deviatoric parts as in:
σ = – pI + τ
1 ∂ρ
α p = – ---
ρ∂T
for ideal gases, the thermal expansion coefficient takes the simpler form αp = 1 ⁄ T
• p is the pressure (SI unit: Pa)
• τ is the viscous stress tensor (SI unit: Pa)
• Q contains heat sources other than viscous dissipation (SI unit: W/m3)
For a steady-state problem the temperature does not change with time and the terms
with time derivatives disappear.
The first term of the right-hand side of Equation 4-18 is the work done by pressure
changes and is the result of heating under adiabatic compression as well as some
thermoacoustic effects. It is generally small for low Mach number flows.
Q vd = τ: ∇u (4-20)
DANCKWERTS CONDITION
The application of the Danckwerts condition on the enthalpy allows to express the
normal conductive heat flux at the inlet boundary as proportional to the flow rate ρu
and the enthalpy variation ΔH between the upstream conditions and inlet conditions:
k ∇T ⋅ n = ρΔHu ⋅ n (4-21)
The enthalpy variation between the upstream conditions and inlet conditions, ΔH,
depends in general both on the difference in temperature and in pressure, and is
defined as:
T pA
1
ΔH = Tupstream Cp dT + pupstream --ρ- ( 1 – αp T ) dp (4-22)
In the unexpected case of a velocity field corresponding to an outgoing flow across the
inlet boundary, a zero conductive flux condition is applied to avoid a unphysical
conductive flux condition:
k ∇T ⋅ n = 0, n⋅u≥0
T – T upstream = 0
For low flow rates or in the presence of large heat sources or sinks next to the inlet, the
conductive heat flux cannot be neglected. The first integral in Equation 4-22 has for
effect to adjust the inlet temperature to balance the energy brought by the flow at the
inlet and the energy transferred by conduction from the interior.
In addition to the cases where the upstream and inlet absolute pressures are equal, this
term may be neglected when the work due to pressure changes is not included in the
energy equation, or when the fluid is modeled as an ideal gas (in this case the
coefficient of thermal expansion is the inverse of the temperature).
∂T ∂p
ρC p ------- + u ⋅ ∇ T + ∇ ⋅ ( q + q r ) = + u ⋅ ∇p + τ: ∇u + Q H + Q (4-23)
∂t ∂t
involving the same variables as in Theory for Heat Transfer in Fluids, and the
additional QH term, accounting for the diffusive flux of thermal enthalpy due to the
rate of change of air and vapor in moist air. It is defined as (Ref. 5):
Q H = – ( C p, v – C p, a )g w ⋅ ∇T
where:
• Cp,v is the specific heat capacity at constant pressure of vapor (SI unit: J/(kg·K))
• Cp,a is the specific heat capacity at constant pressure of air (SI unit: J/(kg·K))
• gw is the vapor flux by diffusion (SI unit: kg/(m2·s))
Humidity
This part defines the different definitions of humidity in the moist air theory.
MOISTURE CONTENT
The moisture content (also called mixing ratio or humidity ratio) is defined as the
ratio of water vapor mass, mv, to dry air mass, ma:
mv pv Mv
x vap = -------- = --------------- (4-24)
ma pa Ma
where pv is the water vapor partial pressure, pa is the dry air partial pressure, and Ma
and Mv are the molar mass of dry air and water vapor, respectively. The moisture
content represents a ratio of mass, and it is thus a dimensionless number.
RELATIVE HUMIDITY
The relative humidity of an air mixture is expressed as follows:
where pv is the water vapor partial pressure and psat is the saturation pressure of water
vapor.
According to Dalton’s law, the total pressure of a mixture of gases is the sum of all the
partial pressures of each individual gas; that is, p = pv + pa where pa is the dry air
partial pressure.
The relative humidity formulation is often used to quantify humidity. However, for the
same quantity of moisture content, the relative humidity changes with temperature
and pressure, so in order to compare different values of φ w it has to be at the same
temperature and pressure conditions. Then the thermodynamical properties of moist
air can be deduced through the mixture formula described below.
The relative humidity is useful to study the condensation as it defines the boundary
between the liquid phase and the vapor phase. In fact, when the relative humidity φ w
reaches unity, it means that the vapor is saturated and that water vapor may condense.
SPECIFIC HUMIDITY
The specific humidity is defined as the ratio of water vapor, mv, to the total mass,
mtot = mv + ma:
mv
ω = ------------ (4-26)
m tot
When the water vapor only accounts for a few percent in the total mass,
the moisture content and the specific humidity are very close: xvap ≈ ω
(only for low values). For larger values of ω, the two quantities are more
precisely related by:
ω
x vap = -------------
1–ω
CONCENTRATION
The concentration is defined by:
nv
c v = ------ (4-27)
V
p sat ( T )
c sat = -------------------- (4-28)
RT
Saturation State
The saturation state is reached when the relative humidity reaches one. It means that
the partial pressure of the water vapor is equal to the saturation pressure (which also
depends on the temperature).
From Ref. 36, the saturation pressure can be defined using the following expression:
T – 273.15 [ K ]
7.5 --------------------------------------
T – 35.85 [ K ]
p sat ( T ) = 610.7 [ Pa ] ⋅ 10 (4-29)
PRELIMINARY DEFINITIONS
Molar Fraction
The molar fraction of dry air, Xa, and the molar fraction of water vapor, Xv, are defined
such as:
na pa p – φ w p sat
X a = ---------- = ------ = --------------------------- (4-30)
n tot p p
nv pv φ w p sat
X v = ---------- = ------ = ------------------ (4-31)
n tot p p
where:
From Equation 4-30 and Equation 4-31, the following relation holds:
Xa + Xv = 1
c v RT x vap p
φ w = -------------------- = ---------------------------------------------------- (4-32)
p sat ( T ) Mv
p sat ( T ) -------- + x vap
Ma
MIXTURE PROPERTIES
The thermodynamical properties are built through a mixture formula. The expressions
depend on dry air properties and pure steam properties and are balanced by the mass
fraction.
Density
According to the ideal gas law, the mixture density ρm expression is defined as follows:
p
ρ m = -------- ( M a X a + M v X v ) (4-33)
RT
where Ma and Mv are the molar mass of dry air and water vapor, respectively, and Xa
and Xv are the molar fraction of dry air and water vapor, respectively.
The ideal gas assumption sets the compressibility factor and the
enhancement factor to unity. In fact, the accuracy lost by this assumption
is small as the pure steam represents a small fraction.
Mm = Xa Ma + Xv Mv
and where Cp, a and Cp, v are the heat capacity at constant pressure of dry air and
steam, respectively.
Dynamic Viscosity
According to Ref. 37 and Ref. 38, the dynamic viscosity is defined as:
Xi μi
μm = ---------------------------- (4-35)
i = a, v X j ϕ ij
j a v
where ϕ ij is given by
1 1 2
--- ---
μi 2 Mj 4
1 + ----- -------
μj Mi
ϕ ij = -----------------------------------------------
1
-
---
Mi 2
8 1 + -------
M j
Here, μa and μv are the dynamic viscosity of dry air and steam, respectively.
Thermal Conductivity
According to Ref. 38 and Ref. 37, the thermal conductivity of the mixture is defined
similarly:
Xi ki
km = ---------------------------- (4-36)
i = a, v X j ϕ ij
j a v
where ka and kv are the thermal conductivity of dry air and steam, respectively.
The valid temperature range is 200 K < T < 1200 K for dry air properties and
273.15 K < T < 873.15 K for steam properties.
FUNCTIONS
The following functions are defined and can be used as feature parameters as well as in
postprocessing. Here, feature stands for fluid or porous, depending on whether the
function is defined in the Fluid or in the Porous Medium feature:
x vap p
c v = -----------------------------------------
x Mv
+ -------- RT
vap M a
∂T
ρCp + ∇ ⋅ q = ρ b C p, b ω b ( T b – T ) + Q met (4-37)
∂t
For a steady-state problem the temperature does not change with time and the terms
with time derivatives disappear.
Damaged Tissue
Add a Thermal Damage subnode under the Biological Tissue node to calculate tissue
damage.
Temperature Threshold
In the first form of damage integral, tissue necrosis occurs in four cases:
• When the temperature exceeds the hyperthermia damage temperature Td, h for
more than a certain time period td, h,
• When the temperature falls below the cryogenic damage temperature Td, c for more
than a certain time period td, c,
• Instantly after the temperature exceeds the hyperthermia necrosis temperature
Tn, h,
• Instantly after the temperature falls below the cryogenic necrosis temperature Tn, c.
For the first two cases, the damaged tissue indicator, α, defined either by
t
1
α = ----------
t d, h 0 ϕ d, h dt
for hyperthermia analysis, or by
t
1
α = ---------
t d, c 0 ϕd, c dt
for cryogenic analysis, with
1 if T > T d, h 1 if T < T d, c
ϕ d, h ( t ) = ϕ d, c ( t ) =
0 otherwise 0 otherwise
is the ratio of the period of time when T > Td, h to the time limit td, h, or the ratio of
the period of time when T < Td, c to the time limit td, c. It gives an indication of
For the last two cases, the necrosis time indicator, αnecr, defined either by
t
α necr = 0 ϕn, h dt
for hyperthermia analysis, or by
t
α necr = 0 ϕ n, c d t
for cryogenic analysis, with
1 if T > T n, h 1 if T < T n, c
ϕ n, h ( t ) = ϕ n, c ( t ) =
0 otherwise 0 otherwise
evaluates the period of time when T > Tn, h or the period of time when T < Tn, c. If
αnecr > 0, the tissue is necrotic because it already reached the necrosis temperatures
Tn, h or Tn, c at some time step of the simulation. Hence, the fraction of necrotic tissue
due to immediate necrosis is equal to 1 if αnecr > 0 and 0 otherwise.
Combining all cases, the overall fraction of necrotic tissue, θd, is equal to:
1 if α necr > 0
θd = (4-38)
min (α, 1) otherwise
Arrhenius Kinetics
The second form of damage integral is applicable only for hyperthermia processes and
provides the degree of tissue injury, α, based on the polynomial Arrhenius equation:
– ΔE
∂α -----------
------- = ( 1 – α ) n Ae RT
∂t
Here, A is the frequency factor (SI unit: 1/s), and ΔE is the activation energy for the
irreversible damage reaction (SI unit: J/mol). The parameters A and ΔE are dependent
on the type of tissue and have been characterized for liver tissues by Jacques et others
(Ref. 7) to be A = 7.39 ⋅ 1039 s–1 and ΔE = 2.577 ⋅ 105 J/mol. See Ref. 8, Ref. 9, and
Ref. 10 for the characterization of these parameters for prostate, skin, and fat. See also
Ref. 11 and Ref. 12 for more references on biological tissues material properties.
θ d = min(max ( α, 0 ), 1) (4-39)
Thermal Properties
The material properties of the damaged tissue are redefined to take into account the
influence of tissue injury. If ρd, Cp, d, and kd denote the density, heat capacity at
constant pressure, and thermal conductivity of the necrotic tissue, respectively, then
two effective quantities are defined:
In these equalities, θd takes one of the two definitions given above in Equation 4-38
or Equation 4-39 according to the integral form chosen.
Heat Source
A cooling or heating source is associated with the reaction leading to damage of tissue.
Depending on the damage integral model, this source is expressed as follows:
• Energy absorption:
∂θ d
Q = – ρL ----------
∂t
∂T s
ρ s C p, s --------- + ∇ ⋅ q s = Q s
∂t
and for a fluid domain where pressure work and viscous dissipation are neglected,
Equation 4-18 becomes:
∂T f
ρ f C p, f --------- + ρ f C p, f u f ⋅ ∇T f + ∇ ⋅ q f = Q f
∂t
The local thermal equilibrium hypothesis assumes equality of temperature in both fluid
and solid phases:
Tf = Ts = T (4-40)
When this hypothesis can be assumed, the heat transfer equation for porous media is
derived from the mixture rule on energies appearing in solid and fluid heat transfer
equations (see Ref. 13). This rule applies by multiplying the equation of the solid
domain by the solid volume fraction, θs, multiplying the fluid equation by the porosity,
εp, and summing resulting equations.
The theory for this hypothesis is detailed in the Local Thermal Equilibrium section
below. Otherwise, the Local Thermal Nonequilibrium section describes the theory for
modeling heat transfer in porous media using two temperatures.
In the particular case of a packed bed of pellets, and under some assumptions about
thermal conductivities ratio, the local thermal nonequilibrium model can be tuned to
capture more precisely some thermal effects within the pellets. This is done by
replacing the macroscale heat equation for the solid phase by an ordinary differential
In the case of conduction in porous plates, Ref. 33 provides criteria based on the
dimensionless Sparrow number, Sp, to indicate if temperature equilibrium is still valid
or if a nonequilibrium point of view should be preferred. In Ref. 34, the influence of
the Darcy number, Da, and the ratio of phase conductivities is examined for transient
heat transfer in packed beds. The Sparrow and Darcy numbers are defined by:
2
h sf L κ
Sp = --------------- Da = -----2-
k eff r h d
where:
• hsf is the interstitial heat transfer coefficient between solid and fluid phases (SI unit:
W/(m2·K))
• L is the plate layer thickness (SI unit: m)
• keff is the equivalent thermal conductivity of the porous medium (SI unit:
W/(m·K))
• rh is the hydraulic radius (SI unit: m)
• κ is the permeability (SI unit: m2)
• d is the average particle diameter (SI unit: m)
In the situations described in Ref. 33 and Ref. 34, small values of Sp (less than 100 or
500) and large values of Da (from order of magnitude 10-7) indicate discrepancies of
temperature in each phase. However, in general, assessing the validity of local thermal
equilibrium assumption remains not straightforward in specific situations. The Local
∂T
( ρC p ) eff ------- + ρ f C p, f u ⋅ ∇T + ∇ ⋅ q = Q (4-41)
∂t
q = – k eff ∇T (4-42)
• εp is the porosity.
• θs is the solid matrix volume fraction.
• ρs is the solid matrix density.
• Cp,s is the solid matrix heat capacity at constant pressure.
• keff is the effective thermal conductivity (a scalar or a tensor if the thermal
conductivity is anisotropic).
• q is the conductive heat flux.
• u is the velocity field, either an analytic expression or computed from a Fluid Flow
interface. It should be interpreted as the Darcy velocity, that is, the volume flow rate
per unit cross sectional area. The average linear velocity (the velocity within the
For a steady-state problem the temperature does not change with time, and the terms
with time derivatives of Equation 4-41 disappear.
The effective thermal conductivity of the solid-fluid system, keff, is related to the
conductivity of the solid, ks, and to the conductivity of the fluid, kf, and depends in a
complex way on the geometry of the medium. In Ref. 13, Ref. 14, and Ref. 18, the
following models are proposed:
• If the heat conduction occurs in parallel in the solid and the fluid, then the effective
thermal conductivity is the weighted arithmetic mean of the conductivities kf and
ks:
k eff = θ s k s + ε p k f
This volume average model provides an upper bound for the effective thermal
conductivity.
• If the heat conduction takes place in series, with all of the heat flux passing through
both solid and fluid, then the effective thermal conductivity is the weighted
harmonic mean of the conductivities kf and ks:
1- θs εp
-------- = ----- + -----
k eff ks kf
This reciprocal average model provides a lower bound for the effective thermal
conductivity.
• A good estimate is given by the weighted geometric mean of kf and ks, as long as kf
and ks are not too different from each other
θs εp
k eff = k s ⋅ k f
• When the porous medium is composed of solid spherical inclusions in a fluid phase,
the effective thermal conductivity can be expressed as:
2k f + k s – 2 ( k f – k s )θ s
k eff = k f ----------------------------------------------------------
2k f + k s + ( k f – k s )θ s
• In the case of a wrapped screen, the effective thermal conductivity of the solid-fluid
system can be defined as:
k f + k s – ( k f – k s )θ s
k eff = k f ---------------------------------------------------
k f + k s + ( k f – k s )θ s
• Finally, when the porous medium is made of sintered metal fibers, the effective
thermal conductivity of the solid-fluid system is defined as:
2 4ε p θ s k f k s
k eff = ε p k f + θ 2s k s + --------------------------
kf + ks
Note that when kf and ks are equal, all these models give the same effective thermal
conductivity. In addition, the effective conductivity is equal to kf when the porosity is
1, and to ks when the porosity is 0.
∂T s
θ s ρ s C p, s --------- + ∇ ⋅ q s = q sf ( T f – T s ) + θ s Q s (4-43)
∂t
q s = – θ s k s ∇T s
∂T f
ε p ρ f C p, f -------- + ρ f C p, f u ⋅ ∇T f + ∇ ⋅ q f = q sf ( T s – T f ) + ε p Q f (4-44)
∂t
In these expressions:
θ s ρ s C p, s T s + ε p ρ f C p, f T f
T = -----------------------------------------------------------------
θ s ρ s C p, s + ε p ρ f C p, f
q sf = S b h sf
The specific surface area, Sb (SI unit: 1/m), for a bed packed with spherical particles
of average diameter dpe is:
6θ s
S b = ---------
d pe
1 d pe d pe
------- = -------------- + ---------
h sf k f Nu βk s
1⁄3 0.6
Nu = 2.0 + 1.1Pr Re p
The Prandtl number, Pr, and particle Reynolds number, Rep, are defined by:
μC p, f d pe ρ f u
Pr = -------------- Re p = ------------------------
kf μ
As in the standard local thermal nonequilibrium model, two temperatures, one for the
fluid phase and one for the pellets, are solved for. By assuming that the thermal
conductivity of the pellets is much smaller than the one of the fluid phase, heat
conduction among different pellets is neglected, and a 1D microscale equation for
temperature conduction along the radial coordinate in the pellets is defined to replace
the macroscale heat transfer equation. It is coupled to the macroscale heat transfer
equation in the fluid phase, either by assuming temperature continuity, or a convective
heat flux, at the outer surface of the pellets.
This model provides higher accuracy than the standard local thermal nonequilibrium
model for highly nonlinear problems, in particular when combustion occurs within the
pellets.
∂T f
ε p ρ f C p, f -------- + ρ f C p, f u ⋅ ∇T f + ∇ ⋅ q f = Q pe, f + ε p Q f (4-45)
∂t
q f = – ε p k f ∇T f
In these expressions:
∂T pe ∂ ∂T pe
( ρC p ) pe, eff ------------ + ------------------2 ----- – r k pe, eff ------------- = Q pe
1 2
(4-46)
∂t ∂ r ∂r
( rr pe )
where
• (ρCp)pe,eff is the effective volumetric heat capacity at constant pressure inside the
pellet filled with fluid (SI unit: J/(kg·K)), defined by
• ρpe is the density of the porous matrix within the pellet (SI unit: kg/m3)
• Cp,pe is the volumetric heat capacity at constant pressure of the porous matrix within
the pellet (SI unit: J/(kg·K))
• εpe is the porosity within the pellet (dimensionless)
• Tpe is the temperature along the radial coordinate of the pellet (SI unit: K)
• rpe is outer radius of the pellet (SI unit: m)
• kpe,eff is the effective thermal conductivity of the pellet filled with fluid (SI unit:
W/(m·K)), defined by
k pe, eff = ( 1 – ε pe )k pe + ε pe k f
• kpe is the thermal conductivity of the porous matrix within the pellet (SI unit:
W/(m·K))
• Qpe is any heat source applied to the pellet (SI unit: W/m3)
The heat flux at the outer surface of the pellet, applied in the microscale pellet
equation, is
with hpe,f the interstitial heat transfer coefficient (SI unit: W/(m2·K)), defined as in
the standard local thermal nonequilibrium model.
The corresponding heat exchange term applied in the macroscale heat transfer
equation of the fluid is:
Q pe, f = – S b q pe r=1
with Sb the specific surface area of the pellets. For a packed bed of spherical pellets, it
is defined as
6 ( 1 – εp )
S b = -----------------------
d pe
When continuity of temperatures at the fluid-pellet interface can be assumed, i.e. when
T pe r=1 = Tf
∂T pe
= – -------- k pe, eff -------------
1
q pe r=1 r pe ∂r r = 1
∂T
( ρC p ) eff ------- + ( ρ g C p, g u g + ρ l C p, l u l ) ⋅ ∇T + ∇ ⋅ q = Q + Q evap (4-47)
∂t
q = – k eff ∇T (4-48)
• ρg and ρl (SI unit: kg/m3) the moist air and liquid water densities.
• Cp,g and Cp,l (SI unit: J/(kg·K)) the moist air and liquid water heat capacities at
constant pressure.
See Theory for Moisture Transport in Porous Media for the definition of gw and glc.
• Qevap (SI unit: W/m3) the heat source (or sink) due to water phase change, defined
as:
Q evap = L v G evap
See Theory for Moisture Transport in Porous Media for the definition of Gevap.
• Matter exists in a certain phase, liquid, solid or gas, at equilibrium under a given
condition of pressure and temperature. The phase transition occurs when the
temperature and pressure of the system come across a critical condition at which the
material changes phase. The critical value of temperature at which the transition
occurs for a given reference pressure can be computed through statistical physics or
tabulated from experiments. In the following, we assume that the phase change is
driven by the temperature, neglecting all other factors that may affect phase change,
as equilibrium time, or purity of the material.
• At the macroscopic and mesoscopic level, phase change costs or releases a certain
amount of energy which is called latent heat of phase transition, which is an input
of the model. Because there is a balance between sensible and latent heat, phase
change induces temperature variations.
The dedicated features for phase change simulation that are available in the Heat
Transfer module handle models were it is assumed that the different phases are non
miscible so that each phase is located in a distinct domain. The simulation predicts the
position of the phase change interface and the associated thermal effects. The
following sections present the two numerical methods available to do the modeling:
It applies to the case of phase change between a solid (or immobile fluid) and a fluid
phase, with a sharp transition between the phases, and no topology changes.
The Phase Change Interface, Exterior boundary condition is applicable on the exterior
boundaries adjacent to a solid, fluid, or porous medium domain.
The latent heat of phase change is taken into account through the definition of the
phase change velocity.
The Stefan condition defines the phase change interface velocity vn as follows:
q ⋅ n mesh
v n = -------------------------
ρs Ls → f
with L s → f the latent heat of phase change from solid to fluid (SI unit: J/kg), ρs the
solid density (SI unit: kg/m3), and q the conductive heat flux jump across the interface
(SI unit: W/m2), defined as
q = – k s ∇T s + k f ∇T f
with ks the solid thermal conductivity (SI unit: W/(m·K)), kf the fluid thermal
conductivity (SI unit: W/(m·K)), and Ts, Tf the temperatures on each side of the
interface.
When the density changes during the phase change, the volume change is
compensated by a velocity at the interface in the fluid phase. The normal fluid velocity
is expressed as:
ρf – ρs
v f = ---------------- v n
ρf
The phase change velocity vn appearing in Stefan condition is relative to the solid
position. In case of translation of the solid (using Translational Motion subfeature
under Solid feature for instance), as in a continuous casting process, the solid
translation velocity us contributes to vn to describe the interface velocity relative to the
spatial frame dx/dt:
dx-
------ ⋅ n = vn + us ⋅ n
dt
Similarly, the same term is added to the normal fluid velocity definition relative to the
spatial frame:
It applies well to materials showing a mushy zone around the phase change interface.
And it allows topology changes of the phase change interface.
The Phase Change Material domain condition should be used on both phases
domains, to solve the heat equation after specifying the properties of the phase change
material according to the apparent heat capacity formulation.
Instead of adding a latent heat L in the energy balance equation exactly when the
material reaches its phase change temperature Tpc, it is assumed that the
transformation occurs in a temperature interval between Tpc − ΔT ⁄ 2 and Tpc + ΔT ⁄ 2.
In this interval, the material phase is modeled by a smoothed function, θ, representing
the fraction of phase before transition, which is equal to 1 before Tpc − ΔT ⁄ 2 and to
0 after Tpc + ΔT ⁄ 2. The density, ρ, and the specific enthalpy, H, are expressed by:
ρ = θρ 1 + ( 1 – θ )ρ 2
1
H = --- ( θρ 1 H 1 + ( 1 – θ )ρ 2 H 2 )
ρ
∂H
Cp =
∂T
1 dα m
C p = --- ( θ 1 ρ 1 C p, 1 + θ 2 ρ 2 C p, 2 ) + ( H 2 – H 1 ) -----------
ρ dT
Figure 4-2: Phase indicators, phase change temperature, and transition interval.
The mass fraction, αm, is defined from ρ1, ρ2 and θ according to:
1 θ2 ρ2 – θ1 ρ1
α m = --- ------------------------------
2 ρ
1
C eq = --- ( θ 1 ρ 1 C p, 1 + θ 2 ρ 2 C p, 2 )
ρ
dα m
C L ( T ) = ( H 2 – H 1 ) -----------
dT
dα m
C L ( T ) = L -----------
dT
ΔT ΔT
T pc + -------- T pc + -------- dα
2 2 m
C ( T ) dT = L ----------- dT = L
ΔT L ΔT
T pc – -------- T pc – -------- dT
2 2
The latent heat, L, can depend on the absolute pressure but should not
depend on the temperature.
Finally, the apparent heat capacity, Cp, used in the heat equation, is given by:
1
C p = --- ( θ 1 ρ 1 C p, 1 + θ 2 ρ 2 C p, 2 ) + C L
ρ
Cp θ 1 ρ 1 C p, 1 + θ 2 ρ 2 C p, 2
γ = ------- = --------------------------------------------------------
Cv C p, 1 C p, 2
θ 1 ρ 1 ------------ + θ 2 ρ 2 ------------
γ1 γ2
k = θ1 k1 + θ2 k2
ρ = θ1 ρ1 + θ2 ρ2
To satisfy energy and mass conservation in phase change models, particular attention
should be paid to the density in time simulations. When the fluid density is not
constant over time, for example, dependent on the temperature, the transport velocity
field and the density must be defined so that mass is conserved locally.
ρ = ρ1 = ρ2
H = θH 1 + ( 1 – θ )H 2
The apparent heat capacity, Cp, used in the heat equation, is given by:
dα m
C p = ( θ 1 C p, 1 + θ 2 C p, 2 ) + L -----------
dT
θ2 – θ1
α m = -----------------
2
∂T
( ρC p ) eff ------- + ∇ ⋅ q = Q (4-49)
∂t
which is derived from Equation 4-15, considering the building material as a porous
medium in local thermal equilibrium in which the following mixing rules apply:
• (ρCp)eff (SI unit: J/(m3·K)) is the effective volumetric heat capacity at constant
pressure, defined to account for both solid matrix and moisture properties:
( ρC p ) eff = ρ s C p, s + w ( φ w )C p, w
where ρs (SI unit: kg/m3) is the dry solid density, Cp,s (SI unit: J/(kg·K)) is the dry
solid specific heat capacity, w ( φ w ) (SI unit: kg/m3) is the water content given by a
moisture storage function, and Cp,w (SI unit: J/(kg·K)) is the water heat capacity at
constant pressure.
• keff (SI unit: W/(m·K)) is the effective thermal conductivity, defined as a function
of the solid matrix and moisture properties:
bw ( φ w )
k eff = k s 1 + --------------------
ρs
where ks (SI unit: W/(m·K)) is the dry solid thermal conductivity and b
(dimensionless) is the thermal conductivity supplement.
This definition neglects the contribution due to the volume fraction change of the
moist air.
The heat source due to moisture content variation is expressed as the vapor diffusion
flow multiplied by latent heat of evaporation:
L v δ p ∇( φ w p sat )
jωt
T ( t ) = T 0 + T'e (4-51)
where T0 is the equilibrium temperature that verifies the steady-state heat transfer
equation and may come from the solution of a previous study, T′ is the complex
amplitude of the harmonic perturbation around T0, and ω is the angular frequency,
related to the ordinary frequency, f, according to
ω = 2πf
Note: The amplitude, T′, is complex-valued since it includes the phase term ejϕ.
From the temperature decomposition in Equation 4-51, and according to the heat
transfer equation in Equation 4-16, heat transfer by conduction in solids is then
governed by:
where Q′ejωt is the harmonic perturbation in domain around an average heat source,
Q. Removing the terms of the steady-state heat transfer equation satisfied by T0, and
simplifying by ejωt, this reduces to:
∂k
jωρ 0 C p, 0 T' + ∇ ⋅ – k 0 ∇T' – T' ∇T = Q' (4-53)
∂ T T 0 0
Here, ρ0, Cp, 0, and k0 denote the density, heat capacity at constant pressure, and
thermal conductivity, evaluated at T0, that is: ρ(T0), Cp(T0), and k(T0), respectively.
When the linearized heat transfer equation, such as Equation 4-52 or Equation 4-53,
can still describe the model accurately, the problem becomes steady-state in the
frequency domain, therefore computationally less expensive than a time-dependent
simulation. An automatic linearization process is performed by COMSOL
Multiphysics so that no additional action is needed from the user to get these
equations, even in the presence of temperature-dependent coefficients, in domains and
boundaries. Only the expressions of the material properties and other parameters, as
functions of the temperature, are required as for usual nonlinear modeling.
Recalling Equation 4-18 given previously in the Theory for Heat Transfer in Fluids
section, without pressure-volume work and viscous dissipation, the equation to be
solved reduces to:
dT
ρC p +∇⋅q = Q
dt
dT
mC p
dt
+ S ( n ⋅ q ) ds = V Q dv (4-54)
where the domain mass and the heat capacity at constant pressure are
1
m = V ρ dv C p = -----
m V ρCp dv
Isothermal domain 3
Isothermal domain 2
Isothermal domain 1
THERMAL INSULATION
The Thermal insulation condition prevents any heat transfer between both adjacent
domains.
CONTINUITY
The Continuity condition ensures equal temperature at both sides of the interface.
VENTILATION
The Ventilation condition is used for cases when an isothermal domain is considered
fluid and has an adjacent domain containing the same fluid. An opening lets the fluid
going from one domain to another with a determined mass flux, denoted by ϕ d → u
or ϕ u → d , respectively, along or opposite to the geometrical normal vector. The
Ventilation condition is written
–nd ⋅ qd = –h ( Tu – Td ) (4-56)
THERMAL CONTACT
When an isothermal domain is considered solid and is adjacent to another solid,
thermal contact occurs and is characterized by a given thermal resistance, Rt. At the
interface, the condition Thermal contact reads
Tu – Td
– n d ⋅ q d = – -------------------- (4-57)
R t, s
In this section:
In addition, standalone physics interfaces are available for the modeling of heat transfer
by conduction, convection and radiation in thin structures:
Either the Solid, Fluid, or Porous Medium feature is available by default in each of these
interfaces.
All these functionalities have in common the fact that the thin domains they model are
lumped into boundaries (for Thin Layer, Thin Film and Fracture) or 3D edges (for
Thin Rod).
Figure 4-4: Modeling a copper wire as a domain (top) requires a denser mesh compared to
modeling it as a boundary with a conductive layer (bottom).
An additional 1D segmented line represents the thickness of the thin structure. The
number of mesh points for each interval of the extra dimension is set to 2 by default.
The normal gradient is the projection of the gradient operator onto the normal
vector, n, of the boundary representing the thin structure. This is mathematically
expressed for any scalar field T as:
∇ n T = ( ∇T ⋅ n )n
The tangential gradient removes the normal component from the gradient
operation, so that only tangential components remain. This is mathematically
expressed for any scalar field T as:
∇ t T = ∇T – ( ∇T ⋅ n )n
The gradient operator is then split into a tangential part and a normal part:
∇T = ∇ t T + ∇ n T (4-58)
∇T = ∇ t T or ∇T = ∇ n T
Equation 4-58 is valid for flat layered shells. For the curved ones, the gradient
expression should account for the surface area and the scaling of each layer. The
gradient in the product geometry of a curved layered shell with variable thickness can
be written as:
∂X –1 T
∇T = --------- ⋅ ( ∇t T + ∇n T )
∂X r
with
∂X-
= I + z + --- ( – 1 + z off ) ∇tn
d
--------
∂X r 2
where
∂X – 1 T ∇n T
∇T = --------- ⋅ ∇ t T + ------------
∂X r lsc
∂X-
= I + z + --- ( – 1 + z off ) ∇t( lscn )
d
--------
∂X r 2
∂X
ASF = --------
-
∂X r
q = – k ∇T
which is Fourier’s law of heat conduction (see also The Physical Mechanisms
Underlying Heat Transfer).
The tensor components can be specified in the local coordinate system of the
boundary, which is defined from the geometric tangent and normal vectors. The local
x direction, ex, loc, is the surface tangent vector t1, and the local z direction, ez, loc, is
the normal vector n. Their cross product defines the third orthogonal direction such
that:
x, loc
e = t1
y, loc
e = n × t1
z, loc
e = n
From this, a transformation matrix between the local coordinate system and the global
coordinate system can be constructed in the following way:
The thermal conductivity tensor in the local coordinate system, kbnd, is then expressed
as
k bnd = AkA T
• The general formulation, using the Extra Dimension tool to solve the equations
into the boundaries and through the thin structure’s thickness
• The thermally thin approximation, a lumped formulation assuming that heat
transfer mainly follows the tangential direction of the thin structure
• The thermally thick approximation, a lumped formulation assuming that heat
transfer is dominant in the direction normal to the thin structure
They all derive from the energy equation established in Equation 4-15, and recalled
here below:
∂E
ρ + ρu ⋅ ∇E + ∇ ⋅ ( q + q r ) = – ( σ:D ) + Q
∂t
GENERAL FORMULATION
The general formulation uses the Extra Dimension tool to solve the equations through
the thin structure’s thickness. The thin structure has its domain represented by the
product space between the lumped boundary and the additional dimension for the
thickness. Applying the split of the gradient operator given earlier at Equation 4-58,
the energy equation becomes
∂E s
ρ + ρu ⋅ ( ∇ t E s + ∇ n E s ) + ∇ ⋅ ( q + q r ) = – ( σ:D ) + Q (4-59)
∂t
The ∇t operator is the tangential derivative in the thin structure boundary, and the ∇n
operator is the derivation operator along the extra dimension which is normal to the
thin structure (see Tangential and Normal Gradients). The subscript s appended on E
(and T in the following) is here to recall that this variable lives in the product space of
the thin structure.
q = –k ( ∇t Ts + ∇n Ts ) (4-60)
Here, ds is the length of the extra dimension, or equivalently the thickness of the thin
structure, and Tu and Td are the temperature at the upside and the downside of the
thin structure.
∇T = ∇ t T
This assumption is often valid for thin structures that are good thermal conductors
compared to the adjacent domains, and/or with fast convection along the tangential
direction.
∂E
ds ρ + d s ρu ⋅ ∇ t E + ∇ t ⋅ ( q s + q r ) = – d s ( σ:D ) + d s Q + q 0 (4-62)
∂t
q s = – d s k∇ t T (4-63)
where ds is the layer thickness (SI unit: m). The heat source Q is a density distributed
in the layer while q0 is the received out-of-plane heat flux.
In 2D, Equation 4-62 and Equation 4-63 have an additional factor, dz,
to account for the out-of-plane thickness.
q0 = n ⋅ q
In this coupling relation, the outgoing heat flux n ⋅ q leaves the domain and is received
in the source term q0 by the adjacent thin layer modeled as a boundary. From the point
of view of the domain, and neglecting thermoelastic effects, the following heat source
is received from the thin structure:
∂E
–n ⋅ q = ds Qs – ds ρ – ( d s ρu ⋅ ∇ t E ) – ∇ t ⋅ ( q s + q r ) (4-64)
∂t
∇T = ∇ n T
This assumption is often valid for thin structures that are thermally resistive compared
to the adjacent domains.
∂E
ds ρ + d s ρu ⋅ ∇ n E + ∇ n ⋅ ( q s + q r ) = – d s ( σ:D ) + d s Q + q 0 (4-65)
∂t
q s = – d s k∇ n T (4-66)
where ds is the layer thickness (SI unit: m). The heat source Q is a density distributed
in the layer while q0 is the received out-of-plane heat flux.
In 2D, Equation 4-65 and Equation 4-66 have an additional factor, dz,
to account for the out-of-plane thickness.
q0 = n ⋅ q
In this coupling relation, the outgoing heat flux n ⋅ q leaves the domain and is received
in the source term q0 by the adjacent thin layer modeled as a boundary. From the point
of view of the domain, and neglecting thermoelastic effects, the following heat source
is received from the thin structure:
∂E
–n ⋅ q = ds Qs – ds ρ – ( d s ρu ⋅ ∇ n E ) – ∇ n ⋅ ( q s + q r ) (4-67)
∂t
T d – 2T 1 ⁄ 2 + T u
∇ n ⋅ ( – d s k∇ n T ) ≈ – k s ------------------------------------------
ds
which can be seen as the sum of two contributive sources on the upside and on the
downside of the boundary that compensate:
Tu – Td Td – Tu
– k s -------------------- – k s --------------------
ds ds
∂T ∂T 1 ⁄ 2 1 ∂T u 1 ∂T d
ρC p ------- ≈ ρC p --------------- = --- ρC p ---------- + --- ρC p ----------
∂t ∂t 2 ∂t 2 ∂t
leading to:
ds Q d s ∂T d Tu – Td
– n d ⋅ q d = ----------- – ρC p ------ ---------- + u d ⋅ n d T d – – k s -------------------- – ( q r, d ⋅ n d ) (4-68)
2 2 ∂t ds
ds Q d s ∂T u Td – Tu
– n u ⋅ q u = ----------- – ρC p ------ ---------- + u u ⋅ n u T u – – k s -------------------- – ( q r, u ⋅ n u ) (4-69)
2 2 ∂t ds
Equations for all supported types of medium are presented in the next sections, Thin
Layer, Thin Film, Fracture, and Thin Rod.
Td Tu
Thin Layer boundary
Figure 4-5: Upside and downside temperatures at a thin layer applied on an interior
boundary. The thin layer is represented by the gray domain.
Downside domain
of the boundary
T = Td
Tu = TextFace
Figure 4-6: Upside and downside temperatures at a thin layer applied on an exterior
boundary.
Thin Layer
Thin layers of solid materials can be considered as boundaries when their thickness is
significantly smaller than the typical lengths of the adjacent domains.
GENERAL FORMULATION
With this formulation, multiple sandwiched layers with different material properties
and thicknesses can be modeled. An additional 1D segmented line represents the
multiple layers in the thin structure. In this extra dimension, the governing equation
is derived from Equation 4-59 to give:
∂T s
ρ si C p, si --------- + ∇ t ⋅ q si = Q si (4-70)
∂t
Td = ( Ts )L = 0
Tu = ( Ts )L = d
s
See Thin Layer (Heat Transfer Interface) and Solid (Heat Transfer in
Shells Interface) with Layer type set as General or more information about
the boundary feature solving Equation 4-70 and Equation 4-71.
The thermally thin approximation is derived from Equation 4-62 to Equation 4-64.
Inside the thin layer, the heat equation becomes:
∂T
d s ρC p, s + ∇t ⋅ qs = ds Q s + q0 (4-72)
∂t
q s = – d s k ∇t T (4-73)
The heat source Qs is a density distributed in the layer while q0 is the received
out-of-plane heat flux.
In 2D, Equation 4-72 and Equation 4-73 have an additional factor, dz,
to account for the out-of-plane thickness.
From the point of view of the domain, the following heat source, derived from
Equation 4-64, is received from the layer:
∂T
– n ⋅ q = d s Q s – d s ρC p, s ------- – ∇ t ⋅ q s (4-74)
∂t
ρ i d s, i
i
ρ = ----------------------
-
d s, i
i
• The homogenized heat capacity is defined as
C p , i d s, i
i
C p = ----------------------------
-
d s, i
i
• The homogenized thermal conductivity is defined as
( roti ⋅ A )
T
⋅ k i ⋅ ( rot i ⋅ A )d s, i
i
--------------------------------------------------------------------------------------
-
k =
d s, i
i
with roti the rotation tensor of layer i:
cos ( θ ) sin ( θ ) 0
rot i = – sin ( θ ) – cos ( θ ) 0
0 0 1
and A the transformation matrix between the local and global coordinates systems:
See Thin Layer (Heat Transfer Interface) and Solid (Heat Transfer in
Shells Interface) with Layer type set as Thermally thin approximation for
more information about the boundary feature solving Equation 4-74. See
The Heat Transfer in Shells Interface for more information about the
physics interface solving Equation 4-72.
ds
R s = ------
ks
The heat flux across the thermally thick structure is derived from Equation 4-67 and
gives
1 ∂T d Tu – Td 1
– n d ⋅ q d = – --- d s ρ s C p, s ---------- – k s -------------------- + --- d s Q s (4-75)
2 ∂t ds 2
1 ∂T u Td – Tu 1
– n u ⋅ q u = – --- d s ρ s C p, s ---------- – k s -------------------- + --- d s Q s (4-76)
2 ∂t ds 2
where the u and d subscripts refer to the upside and downside of the layer, respectively.
dsj
d tot = (4-77)
j=1
d tot
k tot = -----------------
n
- (4-78)
d sj
-------
k sj
j 1
where n is the number of layers.
See Thin Layer (Heat Transfer Interface) and Solid (Heat Transfer in
Shells Interface) with Layer type set as Thermally thick approximation for
more information about the boundary feature solving Equation 4-75 and
Equation 4-76.
Thin Film
Thin films of fluid can be considered as boundaries of thickness significantly smaller
than the typical lengths of the overall model.
GENERAL FORMULATION
With this formulation, heat transfer is modeled in the whole film, including its
thickness. An additional 1D segmented line represents the thickness in the thin film.
In this extra dimension, the governing equation is derived from Equation 4-59 to give:
∂T s
ρC p --------- + ρC p u ⋅ ( ∇tT s + ∇nT s ) + ∇ t ⋅ q f = Q f (4-79)
∂t
q f = – k ( ∇t T s + ∇n T s ) (4-80)
Td = ( Ts )L = 0
Tu = ( Ts )L = d
f
See Thin Film (Heat Transfer Interface) and Fluid (Heat Transfer in
Shells Interface) with Thin film model set as General for more information
about the boundary feature solving Equation 4-79 and Equation 4-80.
∂T
d f ρC p ------- + u ⋅ ∇ t T + ∇ t ⋅ q f = d f Q f + q 0 (4-81)
∂t
q f = – d f k∇ t T (4-82)
where df is the film thickness (SI unit: m). The heat source Qf is a density distributed
in the layer while q0 is the received out-of-plane heat flux.
In 2D, Equation 4-72 and Equation 4-73 have an additional factor, dz,
to account for the out-of-plane thickness.
From the point of view of the domain, the following heat source, derived from
Equation 4-64, is received from the layer:
See Thin Film (Heat Transfer Interface) and Fluid (Heat Transfer in
Shells Interface) with Thin film model set as Thermally thin approximation
for more information about the boundary feature solving Equation 4-81.
See The Heat Transfer in Films Interface for more information about the
physics interface solving Equation 4-83.
Fracture
When fractures occur in porous media, fluid flow tends to move faster than in the bulk
medium. The transport of heat occurs faster in the fractures that in the surrounding
medium, so in this sense, heat transfer in fractures filled with fluids is more similar to
a highly conductive layer than to a thin thermally resistive layer.
The mass transport in fractures can be modeled as Darcy’s law in a thin sheet of porous
medium:
κ
u = --- ∇ t p
μ
where u is the tangential Darcy’s velocity (SI unit: m/s), κ is the fracture permeability
(SI unit: m2), μ the fluid’s dynamic viscosity (SI unit: Pa⋅s), and ∇tp the tangential
gradient of the fluid’s pressure.
Typically, Darcy’s Law with tangential derivatives is solved to compute mass transport,
so in addition to the fluid properties, the fracture should define its own permeability
(or hydraulic conductivity in case the fluid is water), porosity, and fracture thickness.
For heat transfer in fractures, the fracture also needs to define the density of the porous
sheet, heat capacity, and thermal conductivity. The effective thermal conductivity of
the fracture must be adjusted to the fracture porosity and thermal conductivity of the
fluid. In rocks and geological formations, the fracture might also contain highly
conductive material, different than the bulk porous matrix.
The equation to solve for computing heat transfer in fractures is derived from
Equation 4-62 to Equation 4-64 and using the procedure detailed in Theory for Heat
Transfer in Porous Media to apply the mixture rule on solid and fluid internal energies.
The resulting equations are:
q fr = – d fr k eff ∇ t T (4-85)
Here (ρCp)eff is the effective heat capacity at constant pressure of the fracture-fluid
volume, ρ is the fluid’s density, Cp is the fluid’s heat capacity at constant pressure, qfr
is the conductive heat flux in the fracture-fluid volume, keff is the effective thermal
conductivity of the fluid-fracture mixture, and Q is a possible heat source.
From the point of view of the domain, the following heat source, derived from
Equation 4-64, is received from the fracture:
∂T
– n ⋅ q = d fr Q 0 – d fr ( ρC p ) eff – d fr ρC p u ⋅ ∇ t T – ∇ t ⋅ q fr (4-86)
∂t
Thin Rod
The Thin Rod feature is similar to Thin Layer (Heat Transfer Interface) and Solid
(Heat Transfer in Shells Interface) with Layer type set as Thermally thin approximation.
It provides a lumped heat transfer model to model thermally thin rods as edges.
∂T
S ( R ) Q ds = Al Ql – Al ρl Cp, l + ∂ t – ∇t ⋅ ql (4-87)
ql = – Al kl ∇ t T (4-88)
with
2
A l = πr l
The analogy with electrical circuits is made by considering the temperature difference
across the network components as the effort variable, and the heat rate through the
network component as the flow variable.
By solving for heat balance within the circuit, the temperatures at the nodes in the
circuit can be determined.
This section presents the underlying theory of the Lumped Thermal System interface,
and the resulting heat transfer equations that hold.
In this section:
For each two-port component, the temperature difference across the component, ΔT,
is expressed as:
ΔT = T p2 – T p1
Radiative Thermal ΔT -P P
Resistor – ------------
R rad
Heat Pipe ΔT -P P
– --------
R
Expression in the other two-port components
Thermal Capacitor ∂ΔT -P P
– C -----------
∂t
Heat Rate Source Psrc 0 P
Thermoelectric P1,src+P2,src P1,src P2,src
Module
Note that for thermal resistors components, the heat rate at port p1 has the same sign
as the temperature difference ΔT across the component. In practice, when the
temperature at port p2 is higher than the temperature at port p1 (ΔT>0), applying
conductive heat transfer between these two ports is equivalent to apply a heat source
(positive heat rate Pp1) at port p1.
For the heat pipe component, the port p1 should correspond to the evaporator side
(hot side), and the port p2 should correspond to the condenser side (cold side).
Therefore, the temperature difference ΔT across the component should be negative in
For the heat rate source component, the heat rate is fully applied on port p2.
For the thermoelectric module component, the heat rates P1,src and P2,src can be
expressed in different ways. See Theory for the Thermoelectric Module Component
for more details.
ΔT
P = – -------- (4-89)
R
The expression of the thermal resistance used in Equation 4-89 depends on the
geometric configuration.
Plane Shell
When considering steady conduction through a plane shell of surface area A (SI unit:
m²), thickness L (SI unit: m), and constant thermal conductivity k (SI unit: W/(m·K)),
the thermal resistance R is:
L
R = -------
kA
Cylindrical Shell
When considering steady conduction through a cylindrical shell of inner radius ri and
outer radius ro (SI unit: m), height H (SI unit: m), and constant thermal conductivity
k (SI unit: W/(m·K)), the thermal resistance R is:
ro
R = ---------------- ln -----
1
2πkH r i
ΔT overall
P = – -------------------------
R tot
R tot = Ri
i=1
with
where ΔToverall = Tn+1 − T1. Note that by analogy with electrical circuits, the elements
placed in series in the circuit share the heat rate (flow variable).
with
where ΔToverall = T2 − T1. Note that by analogy with electrical circuits, the elements
placed in parallel in the circuit share the temperature difference (effort variable).
E bi – J i
P i = --------------------
Ri
1 – εi
R i = -------------
Ai εi
with Ai the surface area (SI unit: m²) and εi the surface emissivity (dimensionless). See
Theory for Surface-to-Surface Radiation for more details.
This allows a network representation of the radiative heat rate with Ebi and Ji as nodes.
The surface resistance tends toward 0 for a large surface or a surface with large
emissivity, for which the blackbody approximation holds.
When considering the surface of index i as part of an enclosure, the net radiative heat
rate is expressed as:
n n
E bi ( J i – J j )
Pi = A i F ij ( J i – J j ) = ------------------------------
R ij
-
j=1 j=1
1
R ij = -------------
A i F ij
TWO-SURFACE ENCLOSURE
In the network representation of the radiative heat transfer in a two-surface enclosure,
the surface resistance and the space resistance are connected in a serial way. Therefore
the total radiative resistance is expressed as follows:
1 – ε1 1 1 – ε2
R 12 = -------------- + ---------------- + --------------
A 1 ε 1 A 1 F 12 A 2 ε 2
For specific configurations, the surface areas can be evaluated from the geometric
dimensions and the view factors are equal to 1. Therefore the expression above
simplifies (see Ref. 21 for details).
Concentric Spheres
The surface areas A1 and A2 verify
A1 r 2
------- = ----1-
A2 r 2
4 4
A1 σ ( T1 – T2 )
P 12 = ----------------------------------------2-
1- 1 – ε2 r1
---- + -------------- -----
ε1 ε2 r2
A1 = A2
4 4
A1 σ ( T1 – T2 )
P 12 = ------------------------------------
1- ----
---- 1
+ -–1
ε1 ε2
A1 r
------- = ----1-
A2 r2
4 4
A1 σ ( T1 – T2 )
P 12 = --------------------------------------
1- 1 – ε2 r1
---- + -------------- -----
ε1 ε2 r2
THERMAL CAPACITANCE
The thermal capacitance C (SI unit: J/K) is a measure of how much heat a body can
store. It is defined as:
C = VρC p = mC p
with V the volume (SI unit: m3), ρ the density (SI unit: kg/m3), Cp the heat capacity
at constant volume (SI unit: J/(kg·K)), and m the mass (SI unit: kg).
hA
T – T∞ – -------- t
- = e C
-------------------
T0 – T∞
where:
This approximation, referred as the lumped thermal capacitance model (see Ref. 41),
holds when h, C, and A are constant and the gradients of temperature within the body
are expected to be smaller than the gradients of temperature between the body and the
surrounding. It happens for example when the thermal contact between the solid and
the fluid is poor, or when the solid is a good thermal conductor.
where L is a relevant length scale of the body, and k is its thermal conductivity.
The lumped thermal capacitance model is usually assumed to be valid when Bi<0.1.
Regarding the network representation of thermal systems, the Biot number may be
used to determine how many nodes should be included, assuming an homogeneous
temperature distribution in the corresponding domains.
DEVICE DESCRIPTION
Although various geometric configurations are available, a heat pipe includes a vapor
channel delimited by a solid wall, with a porous wick in between, see Figure 4-9 below.
The working fluid flows from the hot side to the cold side of the heat pipe under vapor
state in the core channel, and under liquid state on the way back through the porous
wick by capillary action. It evaporates when leaving the wick to the core channel on
the hot side, named evaporator side, and condensates when entering the wick on the
cold side, named condenser side.
The heat is transferred by conduction through the wall, and by conduction and
convection in the wick and in the vapor channel. The latent heat absorbed and released
by the evaporation and condensation of the working fluid makes heat pipes very
effective heat transfer devices, with large effective thermal conductivities.
By making the analogy with electrical circuits, the total thermal resistance R (SI unit:
W/K) is expressed as:
where
1
R a = ------------------------------------------------------------------
1 1 1
------------ + -------------------- + -------------------
R v, a R wick, a R wall, a
with
P
ΔT = – ----
R
See Ref. 18 for details about the network representation of heat pipes and the
expressions of effective thermal conductivity of the wick used to express its thermal
resistance.
DEVICE DESCRIPTION
They are composed of thermocouples, each consisting of p-type and n-type
semiconductors, connected electrically in series and thermally in parallel and
sandwiched between two high thermally conductive but low electrically conductive
ceramic plates, as described on Figure 4-11 below.
See Theory for the Thermoelectric Effect Interface for details about the Peltier effect.
• Hot junction:
2
Re I ΔT ∂T hot
P hot = – SIT hot + ------------ + -------- – C hot --------------- (4-90)
2 R ∂t
• Cold junction:
2
Re I ΔT ∂T cold
P cold = SIT cold + ------------ – -------- – C cold ----------------- (4-91)
2 R ∂t
where
For a steady-state problem the temperature does not change with time and the heat
storage terms disappear.
The Peltier effect is a cooling effect at the hot junction, and a heating effect at the cold
junction.
R e = R e, p + R e, n
where Re,p and Re,n (SI unit: Ω) are the electrical resistances of the p-type and n-type
semiconductors.
1- 1 1
--- = ---- + ----
R Rp Rn
where Rp and Rn (SI unit: K/W) are the thermal resistances of the p-type and n-type
semiconductors.
S = Sp – Sn
NETWORK REPRESENTATION
The following network representation corresponds to Equation 4-90 and
Equation 4-91:
See Ref. 18 and Ref. 19 for details about the network representation of thermoelectric
modules.
Removed
Heat, Q
Qmax
Temperature Difference, ΔT
ΔTmax
Qmax is the maximum amount of heat that the thermoelectric cooler can remove,
when there is no temperature difference between the two sides of the module.
The linear curves from (ΔT=0, Q=Qmax) to (ΔT=ΔTmax, Q=0) define the removed
heat Q as:
ΔT
Q = Q max – Q max ⋅ ------------------
ΔT max
Note that this relation is valid when the temperature difference satisfies:
0 ≤ ΔT ≤ ΔT max
Then the powers applied at the cold and hot sides are respectively:
For cases where the current through the device varies, the performance parameters
should be made a function of the current.
General Model
The removed heat may be defined as a more general function of the temperature
difference and other parameters such as the temperature at the hot side.
εn2σT4
ρdG
P P
The total incoming radiative flux at P is called irradiation and denoted G. The total
diffuse outgoing radiative flux at P is called radiosity and denoted J. This radiosity is
the sum of diffusively reflected and emitted radiation:
According to the Stefan-Boltzmann law, eb(T) is the power radiated across all
wavelengths and depends on the forth power of the temperature:
e b ( T ) = n 2 σT 4
The net inward radiative heat flux, q, is then given by the difference between the
irradiation and the outgoing radiation (radiosity and specular reflected radiation):
q = G – ( J + ρs G )
q = ( 1 – ρ s )G – J (4-93)
Using Equation 4-92 and Equation 4-93, J can be eliminated and a general expression
is obtained for the net inward heat flux into the opaque body based on G and T.
q = ( 1 – ( ρ d + ρ s ) )G – εe b ( T ) (4-94)
Most opaque bodies also behave as ideal gray bodies, meaning that the absorptivity and
emissivity are equal, and the reflectivity ρd+ρs is therefore obtained from the following
relation:
α = ε = 1 – ( ρd + ρs ) (4-95)
q = ε ( G – eb ( T ) ) (4-96)
Ju = ρd,uGu + εunu2σT4
εunu2σT4
Gu ρs,uGu
ρd,uGu
τdGd
P P
Gd
Figure 4-14: Upside and downside incoming irradiation (left), upside outgoing radiosity
(right). The downside outgoing radiosity is defined in a similar way.
The total incoming radiative flux at P is called irradiation, and is denoted Gu on the
upside and Gd on the downside. The total diffuse outgoing radiative flux at P is called
radiosity and denoted Ju on the upside and Jd on the downside. This radiosity is the
sum of diffusively reflected radiation, emitted radiation and transmitted radiation
coming from the other side of the semitransparent layer:
J u = ρ d, u G u + ε u e b, u ( T u ) (4-97)
J d = ρ d, d G d + ε d e b, d ( T d ) (4-98)
The net inward radiative heat fluxes on the upside and downside, qu and qd, are then
given by the difference between the irradiation and the radiosity:
q u = ( 1 – ρ s, u – τ u )G u – J u (4-99)
q d = ( 1 – ρ s, d – τ d )G d – J d (4-100)
Bodies are considered to behave as ideal gray bodies, meaning that the absorptivity and
emissivity are equal, and the reflectivity ρs is therefore obtained from the following
relation:
ε u + ρ d , u = 1 – ρ s, u – τ u (4-101)
ε d + ρ d, d = 1 – ρ s, d – τ d (4-102)
q u = ε u ( G u – e b, u ( T u ) ) (4-103)
q d = ε d ( G d – e b, d ( T d ) ) (4-104)
q = ε u ( G u – e b, u ( T u ) ) + ε d ( G d – e b, d ( T d ) ) (4-105)
Incident rays which angle of incidence (measured between the ray and the normal to
the surface) is higher than the critical angle are not transmitted, regardless the
transmittance of the surface. They contribute to total reflection instead. Hence the
directional transmissivity coefficient can be defined as
τ if θ ≤ θ c
τ(θ) =
0 if θ > θ c
ρs ( θ ) + τ ( θ ) = 1 – ( ε + ρd )
we can establish
ρ s if θ ≤ θ c
ρs ( θ ) =
ρ s + τ if θ > θ c
2πn 2 C 1
e b, λ ( λ, T ) = ------------------------------ (4-106)
C2
-------
λ e – 1
5 λT
where:
• the two constants C1 (SI unit: W·m2/sr) and C2 (SI unit: m·K) are given by
2 hc 0
C 1 = hc 0 C 2 = ---------
kB
The integral of eb, λ(λ, T) over a spectral band represents the power radiated on the
spectral band and is defined by
λ2 ∞
λ eb, λ ( λ, T ) dλ
1
= F λ1 T → λ2 T 0 eb, λ ( λ, T ) dλ
where F λ1 T → λ 2 T is the fractional blackbody emissive power,
0 eb, λ ( λ, T ) dλ
Recall the Stefan-Boltzmann law that computes the power radiated across all
wavelengths:
∞
0 eb, λ ( λ, T ) dλ = e b ( T ) = n 2 σT 4
λ2
λ e b , λ ( λ , T ) dλ = F λ T → λ T e b ( T )
1
1 2
Notice that:
F λ 1 T → λ2 T = F 0 → λ2 T – F 0 → λ 1 T and F 0 → ∞ = 1
The figure below shows the value of F 0 → λT for different values of λT.
The assumption that the surface emissivity is independent of the radiation wavelength
is often valid when most of the radiative power is concentrated on a relatively narrow
spectral band. This is likely the case when the radiation is emitted by a surface at
temperatures in limited range.
It is interesting to notice that about 97% of the radiated power from a blackbody at
5800 K is at wavelengths of 2.5 µm or shorter, and 97% of the radiated power from a
blackbody at 700 K is at wavelengths of 2.5 µm or longer (see Figure 4-17).
Many problems have a solar load, but the peak temperatures are below 700 K.
For each surface, properties are then described in terms of a solar absorptivity and an
emissivity.
Solar irradiation,
λ < 2.5 µm
Reradiation to
surroundings,
λ > 2.5 µm
By splitting the bands at the default of 2.5 μm, the fraction of absorbed solar radiation
on each surface is defined primarily by the solar absorptivity.
Emissivity
Wavelength
Figure 4-19: Solar and ambient spectral band approximation of the surface emissivity by
a constant per band emissivity.
The Heat Transfer Module supports constant surface properties per spectral bands and
to adjust spectral intervals endpoints.
Emissivity
λ1 λ2 λ3 Wavelength
The multiple spectral bands approach is used in cases when the surface properties vary
significantly over the bands of interest.
e3
θ
e2
ϕ
e1
• Emissivity
ε tot ( θ, ϕ, T, x ) = f ε ( θ, ϕ ) + ε ( T, x )
• Transmissivity
τ tot ( θ, ϕ, T, x ) = f τ ( θ, ϕ ) + τ ( T, x )
• Specular reflectivity
ρ s, tot ( θ, ϕ, T, x ) = 1 – ρ d – ε tot ( θ, ϕ, T, x ) – τ tot ( θ, ϕ, T, x )
q = ( 1 – ρ s )G – J
J = ρ d G + εe b ( T )
Here
The irradiation, G, at a given point is split into three contributions according to:
where:
• εamb is the ambient emissivity; (SI unit: 1), a dimensionless number in the range
0 ≤ εamb ≤ 1. εamb < 1 means that part of the energy coming from radiative bodies
is not absorbed by the ambient air.
ε + ρd + ρs = 1 ,
ε + ρd + ρs + τ = 1 ,
J = ( ρd ( λ, T )G ( λ ) + ε ( λ, T )eb, λ ( λ, T ) ) dλ
0
where
• ε(λ, T) and ρd(λ, T) are is the hemispherical spectral surface emissivity and diffuse
reflectivity, dimensionless quantities in the range [0,1]. Diffuse-spectral surface
corresponds to a surface properties are dependent on the radiation wavelength and
surface temperature.
• T is the surface temperature (SI unit: K).
• eb, λ(λ, T) is the blackbody hemispherical emissive power (SI unit: W/(m3·sr))
defined in Equation 4-106.
The Surface-to-Surface Radiation Interface assumes that the surface emissivity and
opacity properties are constant per spectral band. It defines N spectral bands (N = 2
when solar and ambient radiation model is used),
J = Ji
i=1
J i = ρ d, i G i + ε i e b ( T )
• External radiation sources on Bi with q0, s, i and Ps, i the external radiation source
heat flux and heat rate, respectively, over Bi:
λi
G ext, i = λ = λ i–1
G ext ( λ ) dλ = F ext, i ( i s )q 0, s, i
or
λi
G ext, i = λ = λ i–1
G ext ( λ ) dλ = F ext, i ( i s )P s, i
When the external source fractional emissive power corresponds to the one of a
blackbody at Text, external radiation sources on Bi can be defined from the external
radiation source heat flux, q0, s, and heat rate, Ps, over all wavelengths:
G ext, i = F ext, i F λi – 1 T → λi T ( i s )q 0, s
or
G ext, i = F ext, i F λi – 1 T → λi T ( i s )P s
ε i + ρ d , i + ρ s, i = 1
ε i + ρ d , i + ρ s, i + τ i = 1
The quantities Gm and Famb in Equation 4-108 are not strictly view factors in the
traditional sense. Instead, Famb is the view factor of the ambient portion of the field
of view, which is considered to be a single boundary with constant radiosity
J amb = e b ( T amb )
On the other hand, Gm is the integral over all visible points of a differential view factor,
multiplied by the radiosity of the corresponding source point. In the discrete model,
think of it as the product of a view factor matrix and a radiosity vector. This is, however,
not necessarily the way the calculation is performed.
Samb
S′
J′
P′
r
n′ n
Jamb
Samb
( – n′ ⋅ r ) ( n ⋅ r )
Gm = S′ -------------------------------------
πr
4
- J′ ds
The heat flux that arrives from P′ depends on the local radiosity J′ projected onto P.
The projection is computed using the normal vectors n and n′ along with the vector
r, which points from P to P′.
The ambient view factor, Famb, is determined from the integral of the surrounding
surfaces S′, here denoted as F′:
( – n′ ⋅ r ) ( n ⋅ r )
F amb = 1 – F ′ = 1 – S′ -------------------------------------
πr
4
- ds
The two last equations plug into Equation 4-107 to yield the final equation for
irradiative flux.
The equations used so far apply to the general 3D case. 2D geometries result in simpler
integrals. For the 2D case, the resulting equations for the mutual irradiation and
ambient view factor are
( – n′ ⋅ r ⊥ ) ( n ⋅ r ⊥ )
F amb = 1 – S ' --------------------------------------------
2r
3
- ds
⊥
⊥
where the integral over S⊥′ denotes the line integral along the boundaries of the 2D
geometry.
• Hemicube
• Direct area integration
• Ray shooting
View factors are always calculated directly from the mesh, which is a
polygonal representation of the geometry. To improve the accuracy of the
radiative heat transfer simulation, the mesh must be refined rather than
raising the element order.
cos θ
------------
2
4πr
cos θ-
-----------
2πr
SOLAR POSITION
The Sun is the most common example of an external radiation source. The position of
the Sun is necessary to determine the direction of the corresponding external radiation
source. The direction of sunlight (zenith angle and the solar elevation) is automatically
computed from the latitude, longitude, time zone, date, and time using similar a
method as described in Ref. 20. The estimated solar position is accurate for a date
between year 2000 and 2199, due to an approximation used in the Julian Day calendar
calculation.
The zenith angle, θs, and azimuth angle, ϕ s , of the Sun are converted into a direction
vector is = (isx, isy, isz) in Cartesian coordinates assuming that the north, the west, and
the up directions correspond to the x, y, and z directions, respectively, in the model.
The relation between θs, ϕ s , and is is given by:
i sx = – cos ( ϕ s ) sin ( θ s )
i sy = sin ( ϕ s ) sin ( θ s )
i sz = – cos ( θ s )
Select between the hemicube and the direct area integration methods also in axial
symmetry. Their settings work the same way as in 3D.
Let I(Ω) denote the radiative intensity traveling in a given direction, Ω. Different kinds
of interactions are observed:
• Absorption: The medium absorbs a fraction of the incident radiation. The amount
of absorbed radiation is κI(Ω), where κ is the absorption coefficient.
• Emission: The medium emits radiation in all directions. The amount of emitted
radiative intensity is equal to κIb, where Ib is the blackbody radiation intensity.
• Scattering: Part of the radiation coming from a given direction is scattered in other
directions. The scattering properties of the medium are described by the scattering
phase function φ ( Ω′, Ω ) , which gives the probability that a ray coming from one
direction Ω′ is scattered into the direction Ω. The phase function φ ( Ω′, Ω ) satisfies:
1
------
4π 4π φ ( Ω′, Ω ) dΩ′ = 1
σs
Ω ⋅ ∇I ( Ω ) = κI b ( T ) – βI ( Ω ) + ------
4π 4π I ( Ω′ )φ ( Ω′, Ω ) dΩ′ (4-115)
• I(Ω) is the radiative intensity at a given position following the Ω direction (SI unit:
W/(m2·sr))
• Ib(T) is the blackbody radiative intensity (SI unit: W/(m2·sr)), defined as
2 4
n r σT
I b ( T ) = ----------------- (4-116)
π
φ ( μ0 ) = 1 + an Pn ( μ0 )
n=1
1 dk 2 k
P k ( x ) = -----------
k
((x – 1) )
k
2 k! d x
2
1 1–η
φ ( μ 0 ) = ---- ⋅ ------------------------------------------------
3⁄2
-
K 2
( 1 + η – 2ημ 0 )
where – 1 < η < 1 is the anisotropy parameter and K is defined as follows to produce
a normalized phase function:
2
1 1–η
K = ------ ⋅
4π ------------------------------------------------
( 1 + η
2
– 2ημ )
3⁄2
- dΩ
4π 0
INCIDENT RADIATION
A quantity of interest is the incident radiation, denoted G (SI unit: W/m2), and
defined by
G = 4π I ( Ω ) dΩ
1–ε
I ( Ω ) = εI b ( T ) + ----------- q r, out for all Ω such that n ⋅ Ω < 0
π
where
2 4
n r σT
I b ( T ) = -------------------- (4-117)
π
q r, out = n ⋅ Ω > 0 I ( Ω ) ( n ⋅ Ω ) dΩ
For black walls ε = 1. Thus I(Ω) = Ib(T).
The net heat flux corresponding to the balance between the energy received by the
surface and emitted by the surface is defined by
q r, net = q r, in – q r, out
This net heat flux accounts for the radiation coming from the domain where the RTE
is solved for, the radiation coming from the exterior is not modeled and hence not
accounted for.
ρd τd
I ( Ω ) = εI b ( T ) + ------ q r, out + ----- q r, in, ext + τ s I ext
π π
where
The net heat flux corresponding to the balance between the energy received by the
surface and emitted by the surface is defined by
For semitransparent surfaces, the net heat flux accounts for the radiation coming from
the domain and from the exterior.
qr = 4π I ( Ω )Ω dΩ
Heat flux divergence can be defined as a function of G and T (see Ref. 23):
Q r = ∇ ⋅ q r = κ ( G – 4πI b ( T ) )
In order to couple radiation in participating media, radiative heat flux is taken into
account in addition to conductive heat flux. Recalling Equation 4-18, the heat transfer
equation reads:
∂T ∂p
ρC p ------- + u ⋅ ∇ T + ∇ ⋅ ( q + q r ) = α p T + u ⋅ ∇p + τ: ∇u + Q
∂t ∂t
∂T ∂p
ρC p ------- + u ⋅ ∇ T + ∇ ⋅ q = κ ( G – 4nσT 4 ) + α p T + u ⋅ ∇p + τ: ∇u + Q
∂t ∂t
Radiative intensity is defined for any direction Ω, because the angular space is
continuous. In order to handle the radiative intensity equation numerically, the angular
space is discretized.
4π I ( Ω ) dΩ ≈ wj Ij
j=1
n
σs
S i ⋅ ∇I i = κI b ( T ) – βI i + ------
4π w j I j φ ( S j, S i )
j=1
1–ε
I i, bnd = εI b ( T ) + ----------- q r, out for all S i such that n ⋅ S i < 0 (4-118)
π
with
q r, out = wj Ij ( n ⋅ Sj )
n ⋅ Sj > 0
In Equation 4-118, the first term in the right hand side is the emitted radiative
intensity, while the second term represents the reflected radiative intensity.
SEMITRANSPARENT SURFACE
The node Semitransparent Surface (Radiation in Participating Media and Radiation in
Absorbing-Scattering Media Interfaces) defines the following boundary condition for
radiative intensities I1, I2, …, In:
ρd τd
I i, bnd = εI b ( T ) + ------ q r, out + ----- q r, out, ext + τ s I i, ext (4-119)
π π
q r, out = wj Ij ( n ⋅ Sj )
n ⋅ Sj > 0
and
In Equation 4-119, the right hand side is composed of the emitted radiative intensity,
the reflected radiative intensity, and the radiative intensity transmitted in a diffuse and
specular way.
RADIATIVE SOURCE
The node Radiative Source accounts for a directional power density Ii in the radiative
transfer equation:
n
σs
S i ⋅ ∇I i = κI b ( T ) – βI i + ------
4π w j I j φ ( S j, S i ) + q i .
j=1
Assuming that a model is invariant in the z direction, the radiative transfer equation in
two directions, Si+ and Si-, for the discrete ordinates method (DOM) reads:
n
σs
S i + ⋅ ∇I i + = κI b ( T ) – βI i+ + ------
4π w j I j φ ( S j, S i )
+
j=1
n
σs
S i- ⋅ ∇I i- = κI b ( T ) – βI i- + ------
4π wj Ij φ ( Sj, Si ) -
j=1
Halving the sum of the two equations above and using I i+ = I i- = I˜i ( I i- = I i+ in
2D) yields
n
σs
S˜ i ⋅ ∇I˜i = κI b ( T ) – βI˜i + ------ wj Ij ( φ ( Sj, Si ) + φ ( Sj, Si ) )
+ -
8π
j=1
j=1
or
n⁄2
σs
S˜ i ⋅ ∇I˜i = κI b ( T ) – βI˜i + ------ wj Ij φ ( Sj , Si ) + wj Ij φ ( Sj , Si )
+ + + + - - - -
8π
j=1
n⁄2
σs
+ ------
8π wj Ij φ ( Sj , Si ) + wj Ij φ ( Sj , Si )
+ + + - - - - +
j=1
In addition
n⁄2
σs
S˜ i ⋅ ∇I˜i = κI b ( T ) – βI˜i + ------ ˜ ˜
wj Ij φ ( Sj, Si )
˜ ˜ (4-120)
4π
j=1
with
S i, 1
S˜ i = S
i, 2
0
with w˜ i = 2w i .
Using results from Equation 4-120 and Equation 4-121 the DOM is formulated in
2D using only radiative intensities, I˜i , on half of the 3D DOM directions, S˜ i , except
for the scattering term. In other expressions than the scattering term, the z component
of the radiative intensities Ii and of the discrete directions Ωi can be ignored (or set to
zero) and the weight wi, multiplied by 2.
P1 Approximation Theory
The P1 approximation is available as a radiation discretization method in The Radiation
in Participating Media Interface.
• The media is optically thick media: τ >>1, where τ is the optical thickness defined by
the integral of absorption coefficient, κ, along a typical optical path:
s
τ = 0 κ ds
• The scattering is linear isotropic.
From a computational point of view this approximation has a limited impact because
it introduces only one additional degree of freedom for the incident radiation G (SI
unit: W/m2), which is a scalar quantity and adds a heat source or sink to the
temperature equation to account for radiative heat transfer contributions. This
method, however, fails to accurately represent cases where the radiative intensity
propagation dominates over its diffusivity or where the scattering effects cannot be
described by a linear isotropic phase function.
σs
Ω ⋅ ∇I ( Ω ) = κI b ( T ) – βI ( Ω ) + ------
4π 4π I ( Ω′ )φ ( Ω′, Ω ) dΩ′
where
1
D P1 = ----------------------------------------
3κ + σ s ( 3 – a 1 )
Q r = κ ( G – 4πI b ) (4-123)
When scattering is modeled as isotropic, a1=0 and the P1 diffusion coefficient reduces
to
1
D P1 = ----------------------
3κ + 3σ s
n ⋅ ( – D P1 ∇G ) = – q r, net
where qr, net is the net radiative heat flux at the boundary.
In addition Qr, defined by Equation 4-123, is added as an heat source in the heat
transfer equation:
∂T ∂p
ρC p ------- + u ⋅ ∇ T + ∇ ⋅ q = Q r + α p T + u ⋅ ∇p + τ: ∇u + Q
∂t ∂t
OPAQUE SURFACE
The Opaque Surface (Radiation in Participating Media and Radiation in
Absorbing-Scattering Media Interfaces) boundary condition defines a boundary
opaque to radiation and defines the incident intensity on a boundary:
The Opaque Surface feature accounts for the net radiative heat flux, qr, net, in the heat
balance.
Gray Wall
The radiative heat flux at the boundary depends on the surface emissivity, ε:
ε
q r, net = -------------------- ( 4πI b, w – G )
2(2 – ε)
with
n2 σT4
I b, w = I b = -----------------
π
Black Wall
The radiative heat flux at the boundary expression simplifies to
1
n ⋅ ( – D P1 ∇G ) = – --- ( 4πI b, w – G ) ,
2
with
n 2 σT 4
I b, w = I b = -----------------
π
SEMITRANSPARENT SURFACE
The Semitransparent Surface (Radiation in Participating Media and Radiation in
Absorbing-Scattering Media Interfaces) boundary condition defines a boundary that
emits radiation, and reflects one part of the incident intensity while the remaining is
transmitted diffusively.
The net radiative heat flux, qr, net, accounted for in the heat balance, is defined as
2 1 – ρd
q r, net = --------------- ( πεI b, w + τ d G ) – ------------------------ G
1 + ρd 2 ( 1 + ρd )
where ρd is the surface diffuse reflectivity, and τd is the surface diffuse transmissivity.
1
q r, net = --- ( 4πI ext – G )
2
which defines the heat radiative heat flux and also contributes to G boundary
condition:
n ⋅ ( – D P1 ∇G ) = – q r, net
RADIATIVE SOURCE
The node Radiative Source accounts for a power density Q in the P1 approximation:
∇ ⋅ ( – D P1 ∇G ) = – Q r + Q .
The radiative intensity I(Ω) at a given position following the Ω direction is the solution
of the radiative transfer equation with no emission term (see Ref. 23):
σs
Ω ⋅ ∇I ( Ω ) = – βI ( Ω ) + ------
4π 4π I ( Ω′ )φ ( Ω′, Ω ) dΩ′ (4-124)
G = 4π I ( Ω ) dΩ
If the Discrete Ordinates Method (DOM) is used for the approximation of
Equation 4-124, G is computed as
G= ωi Ii
i=1
and
N
σs
S i ⋅ ∇I i = – βI i + ------
4π ω j I j φ ( S j, S i )
j=1
where
∇ ⋅ ( – D P1 ∇G ) = – κ G , (4-125)
ε
n ⋅ ( – D P1 ∇G ) = -------------------- G ,
2(2 – ε)
In the Discrete Ordinates Method (DOM), the Radiative Transfer Equation is solved
for each spectral band:
N
σ s, k
S i ⋅ ∇I i, k = κ k I b ( T )FEP k ( T ) – β k I i, k + ----------
4π ω j I j, k φ k ( S j, S i ) ,
j=1
where FEP is the fractional emissive power, also presented in the section Wavelength
Dependence of Surface Properties.
If the P1 Approximation Theory is used instead, the incident radiation Gk for the
spectral band k is solution of the following equation,
∇ ⋅ ( – D P1, k ∇G k ) = – Q r, k , (4-126)
where
1
D P1, k = ----------------------------------------------------
3κ k + σ s, k ( 3 – a 1, k )
• The radiative beam in the absorbing medium is collimated and each beam
propagates always in the same direction.
• The light experiences no refraction, reflection, or scattering within the material
itself.
• There is no significant emission of the material in the wavelength range of the
incident light. This applies well to laser beams, whose wavelength is in general much
shorter than the one of the radiation emitted by the medium.
In these conditions, the radiative intensity Ii (SI unit: W/m2) of the ith beam through
the material decreases as the beam propagates and is absorbed by the medium. This is
described by the Beer-Lambert Law equation:
ei
--------- ⋅ ∇I i = – κI i
ei
where ei is the orientation of the ith beam, and κ is the absorption coefficient (SI unit:
m–1) of the medium.
The radiative heat source Qr (SI unit: W/m3), corresponding to the energy deposited
by the radiative beam, is defined by:
κIi
Qr =
.
i
s
τ = 0 κ ds
From a computational cost point of view this approximation has a limited impact
because it does not introduce any extra degree of freedom to the heat equation.
In this case, the radiative heat flux can be evaluated by (Ref. 23):
4π
q r, λ = – ------ ∇i b, λ
βλ
4σ
q r = – ---------- ∇ ( n 2 T 4 )
3β R
16n 2 σT 3
k R = ------------------------
3β R
and
16n 2 σT 3
q r = – ------------------------ ∇T
3β R
Rosseland mean extinction coefficient, βR, can be defined from the absorption and
scattering coefficients of the media from βR=κ+σs.
Two formulations are available, depending on the magnitude of moist air density
variations with moisture content.
The Moisture Transport in Air Interface solves for the following equation, in which
the moisture content variation is expressed through the transport of vapor
concentration, cv, by convection and binary diffusion in air:
∂c v
M v -------- + M v u ⋅ ∇c v + ∇ ⋅ g w = G (4-128)
∂t
g w = – M v D ∇c v
c v = φ w c sat
Mv cv
ω v = --------------
ρg
The Heat Transfer in Moist Air Interface solves for the following equation:
∂ω v
ρ g ---------- + ρ g u ⋅ ∇ω v + ∇ ⋅ g w = G (4-129)
∂t
g w = – ρ g D ∇ω v
∂ ( cv ⁄ ρg )
M v ρ g ----------------------- + M v ρ g u ⋅ ∇( c v ⁄ ρ g ) + ∇ ⋅ g w = G
∂t
g w = – M v ρ g D ∇( c v ⁄ ρ g )
By neglecting the spatial and temporal variations of ρg in the equations above, the
equation of the diluted species formulation, Equation 4-128, is obtained.
• Adsorption: the vapor is fixed under molecular films on the pores walls. Equilibrium
vapor pressure is reached in the pores under thermodynamic conditions that are
significantly different than those necessary in a free medium. Thus saturation
conditions are changed, and the relative humidity is significantly changed as well,
becoming a relevant variable for the characterization of moisture content.
• Capillarity: when the porous matrix is wetted by liquid water, there is a pressure
difference between the liquid and gas (moist air) phases, due to the curvature of the
wetting interfaces in each pore. The pressure difference, called capillary pressure,
The relative humidity and capillary pressure are called state variables: they describe the
thermodynamic state of water in the pore space under two phases in equilibrium. They
are related through Kelvin’s law, which is a relation of thermodynamical origin. As long
as equilibrium is satisfied, the state variables are continuous between two media (at the
interface between two porous media with different porosities, or between a porous and
a free medium), and Kelvin’s law is valid on both sides of the interface. In opposition,
the liquid saturation or the moisture content depend on the porous matrix.
Note that capillary transport of liquid water should be considered in both regions.
εp sl ρl + εp ρ ωv ( 1 – sl ) = w ( φw )
g
where the liquid saturation, sl, describes the amount of liquid water within the pores,
and εp (dimensionless) is the porosity.
∂w ( φ w )
-------------------
- + ρ g u g ⋅ ∇ω v + ∇ ⋅ g w + u l ⋅ ∇ρ l + ∇ ⋅ g lc = G (4-130)
∂t
M v φ w c sat ( T )
ω v = -----------------------------------
ρg
• Deff (SI unit: m2/s) is the effective vapor diffusion coefficient in the porous
medium, computed from the diffusion coefficient in a free medium, and accounting
for the porosity and tortuosity of the porous medium.
• ul (SI unit: m/s) is the liquid water velocity field, defined from the absolute
pressure gradient by the Darcy’s law as:
κ rl κ
u l = – ----------- ∇p A
μl
• κrl (dimensionless) is the relative liquid water permeability, that may be defined as
a function of the liquid saturation.
• κ (SI unit: m2) is the porous medium permeability.
• G (SI unit: kg/(m3⋅s)) is a moisture source (or sink). See the Moisture Source node.
The moisture source due to evaporation may be obtained from the transport equation
for vapor:
∂ [ εp ρg ωv ( 1 – sl ) ]
G evap = ---------------------------------------------- + ρ g u g ⋅ ∇ω v + ∇ ⋅ g w
∂t
Q evap = L v G evap
For a steady-state problem, the relative humidity does not change with time and the
first term in the moisture transport equation disappears.
Starting from the equation for the transport of moisture in porous media (using the
diffusion model for capillary flux):
∂w ( φ w ) ∂w ( φ w )
- + ρ g u g ⋅ ∇ω v + ∇ ⋅ g w + u l ⋅ ∇ρ l + ∇ ⋅ – D w -------------------
------------------- - ∇φ w = G
∂t ∂φ w
let’s make the assumption that there is no variation of the ambient total pressure
(pA=1 atm). In this case, by defining the gas and liquid velocities ug and ul as Darcy
velocities, the convective fluxes for vapor and liquid water vanish:
∂w ( φ w ) ∂wφ w
- + ∇ ⋅ g w – D --------------
------------------- - ∇φ w = G
∂t w ∂φ w
∂w ( φ w )
ξ = --------------------
∂φ w
∂φ w
ξ ---------- + ∇ ⋅ ( g w – D w ξ ∇φ w ) = G
∂t
Ma Mv pv
g w = – ρ g D eff ∇ω v = – ρ g D eff ∇ ------------------------------------------2- ------
( xa Ma + xv Mv ) A p
By making the approximation that xaMa+xvMv is constant, and that moist air is an
ideal gas, we have:
M v D eff Ma
g w = – ------------------- ------------------------------------------2- ∇p v
RH2 O ( x M + x M )
a a v v
M v D eff Ma
δ p = ------------------- ------------------------------------------2-
RH2 O ( x M + x M )
a a v v
∂φ w
ξ ---------- + ∇ ⋅ ( – ξ D w ∇φ w – δ p ∇ ( φ w p sat ( T ) ) ) = G (4-131)
∂t
This equation models the moisture transfer as the sum of the capillary moisture flux:
∂w
– D w ∇ ( w ( φ w ) ) = – D w ---------- ∇φ w = – ξD w ∇φ w
∂φ w
δ p ∇p v ( T ) = δ p ∇ ( φ w p sat ( T ) )
• Theory for the Nonisothermal Flow and Conjugate Heat Transfer Interfaces
• Theory for the Moisture Flow Interface
• Theory for the Thermoelectric Effect Interface
• Theory for the Building Materials Version of the Heat and Moisture Transport
Interface
• Theory for the Moist Air Version of the Heat and Moisture Transport Interface
• Theory for the Moist Porous Media Version of the Heat and Moisture Transport
Interface
• Theory for the Heat and Moisture Flow Interfaces
• Theory for the Electromagnetic Heating Interfaces
• Theory for the Thermal Stress Interface
See Theory for the Single-Phase Flow Interfaces and Theory for the Turbulent Flow
Interfaces in the CFD Module User’s Guide for a description of the theory related to
laminar and turbulent single-phase flow interfaces.
Other situations where the density might vary includes chemical reactions, for instance
where reactants associate or dissociate.
The Nonisothermal Flow and Conjugate Heat Transfer interfaces contain the fully
compressible formulation of the continuity and momentum equations:
∂ρ
------ + ∇ ⋅ ( ρu ) = 0
∂t
(4-132)
∂u
ρ ------- + ρu ⋅ ∇u = – ∇p + ∇ ⋅ τ + F
∂t
where
It also solves the heat equation, which for a fluid is given in Equation 4-18 by
∂T ∂p
ρC p ------- + u ⋅ ∇ T + ∇ ⋅ ( q + q r ) = α p T + u ⋅ ∇p + τ: ∇u + Q
∂t ∂t
1 ∂ρ
α p = – ---
ρ∂T
• Q contains heat sources other than viscous heating (SI unit: W/m3)
∂p
Q p = α p T ------ + u ⋅ ∇ p
∂t
Q vd = τ: ∇u
are not included by default because they are usually negligible. These terms can,
however, be added by selecting corresponding check boxes in the Nonisothermal Flow
feature.
The physics interface also supports heat transfer in solids (Equation 4-16):
∂T
ρC p ------- + u trans ⋅ ∇ T + ∇ ⋅ ( q + q r ) = Q ted + Q
∂t
where Qted is the thermoelastic damping heat source (SI unit: W/(m3)). This term is
not included by default but must be added by selecting the corresponding check box.
Equations for compressible turbulence are derived using the Favre average. The Favre
˜
average of a variable T is denoted T and is defined by
˜ ρT
T = -------
ρ
where the bar denotes the usual Reynolds average. The full field is then decomposed as
˜ ˜ ˜ u˜
∂ ρE
˜ u i u i ρu i ″u i ″ ∂ ρu˜ ˜ + u i i ˜ ρu i ″u i ″
+ ----------- + -------------------- + j H ----------- + u j -------------------- = (4-133)
∂t 2 2 ∂ xj 2 2
∂ – q – ρu ″H″ τ u ″ ρu j ″u i ″u i ″ ∂ ( u˜ ( – ρu ″u ″ ) )
+ ij i – ----------------------------- + τ
∂ xj j j 2 ∂ x j i ij i j
∂T
q j = – λ ------- (4-134)
∂x j
∂u i ∂u j 2 ∂u k
τ ij = μ -------- + -------- – --- μ --------- δ ij
∂x j ∂x i 3 ∂x k
is the laminar, viscous stress tensor. Notice that the thermal conductivity is denoted λ.
The modeling assumptions are in large part analogous to those for incompressible
turbulence modeling. The stress tensor
– ρu i''u'' j
˜
T ∂u˜i ∂u˜ j 2 ∂u k 2
– ρu i ″u j ″ = ρτ ij = μ T -------- + -------- – --- μ --------- δ ij – --- ρkδ ij (4-135)
∂x j ∂x i 3 ∂x k 3
1
ρk = --- ρu i ″u i ″ (4-136)
2
The correlation between uj″ and H″ in Equation 4-133 is the turbulent transport of
heat. It is modeled analogously to the laminar conductive heat flux
˜ μ T C p ∂T ˜
T ∂T
ρu j ″H″ = q j = – λ T ------- = – -------------- ------- (4-137)
∂x j Pr T ∂x j
τ ij u i ″
ρu j ″u i ″u i ″
----------------------------
-
2
ρu j ″u i ″u i ″ μ T ∂k
τ ij u i ″ – ----------------------------- = μ + ------ ------- (4-138)
2 σ k ∂x j
Inserting Equation 4-134, Equation 4-135, Equation 4-136, Equation 4-137 and
Equation 4-138 into Equation 4-133 gives
˜ ˜ ˜ u˜
∂ ρE ˜ ui ui ∂ ρu˜ ˜ + u i i
+ ----------- + k + j H ----------- + k = (4-139)
∂ t 2 ∂ xj 2
∂ q μ T ∂k ∂ u˜
– j – q j + μ + ------ ------- +
T T
( ( τ + ρτ ij ) )
∂ xj σ k ∂x j ∂ x j i ij
The Favre average can also be applied to the momentum equation, which, using
Equation 4-135, can be written
∂ ρu˜ ∂ ρu˜ ˜ ∂p ∂ τ T
( i) + ( j u i ) = – ------- + ( + ρτ ij ) (4-140)
∂t ∂ xj ∂x j ∂ x j ij
Taking the inner product between u˜ i and Equation 4-140 results in an equation for
the resolved kinetic energy, which can be subtracted from Equation 4-139 with the
following result:
∂ (ρ(E ˜ ˜
+ k ) ) + ∂ ( ρu˜ j ( E + k ) ) = (4-141)
∂t ∂ xj
∂ũ j ∂ μ T ∂k
– q j – q j + μ + ------ ------- + ∂ ( u˜i ( τ ij + ρτ ij ) )
T T
–p +
∂ xj ∂ xj σ k ∂x j ∂ xj
˜
∂ ( ρE ˜ ∂ũ j ∂
) + ∂ ( ρu˜ j E ) = – p
T ∂ u˜ T
+ ( – qj – qj ) + ( ( τ + ρτ ij ) ) (4-142)
∂t ∂ xj ∂ xj ∂ xj ∂ x j i ij
τ ij = τ˜ij + τ ij ″
Since
τ˜ij » τ ij ″
τ ij ≈ τ˜ij
and consequently
˜
˜
∂ ( ρE ˜ ∂ũ j ∂ ∂T
) + ∂ ( ρu˜ j E ) = – p ∂ ( u˜ tot )
+ ( λ + λ T ) ------- + τ̃ (4-143)
∂t ∂ xj ∂ xj ∂ xj ∂x j ∂ x j i ij
where
˜
tot ∂u˜ i ∂u˜ j 2 ∂u k
τ̃ ij = ( μ + μ T ) -------- + -------- – --- --------- δ ij
∂x j ∂x i 3 ∂x k
˜ ˜ ˜ ˜
∂T ∂T ∂ λ λ ------
∂T 1 ∂ρ ˜ ∂p ˜ ∂p ˜ ∂u i
ρC p ------- + ρC p u˜ j ------- + – ( + T ) - = – --- ------˜- T ------ + u j ------- + τ ij --------
∂t ∂x j ∂ x j ∂x j ρ ∂T ∂t ∂x j ∂x j
p
Turbulent Conductivity
Kays-Crawford This is a relatively exact model for PrT, while still quite simple. In
Ref. 29, it is compared to other models for PrT and found to be a good approximation
for most kinds of turbulent wall bounded flows except for turbulent flow of liquid
metals. The model is given by
1 0.3C p μ T 0.3C p μ T 2 λ –1
Pr T = ----------------- + ----------------------- – ----------------------- 1 – exp – ------------------------------------------ (4-144)
2Pr T∞ λ Pr T∞ λ 0.3C μ Pr
p T T∞
where the Prandtl number at infinity is PrT∞ = 0.85 and λ is the conductivity.
100λ
Pr T∞ = 0.85 + -------------------------------
-
C p μRe ∞0.888
where Re∞, the Reynolds number at infinity must be provided either as a constant or
as a function of the flow field. This is entered in the Model Inputs section of the Fluid
feature.
The heat flux between the fluid with temperature Tf and a wall with temperature Tw, is:
ρC p u τ ( T w – T f )
q wf = -----------------------------------------
+
-
T
where ρ is the fluid density, Cp is the fluid heat capacity, and uτ is the friction velocity.
T+ is the dimensionless temperature and is given by (Ref. 31):
where in turn
+ δ w ρ C μ1 / 2 k + 10
δ w = ------------------------------- δ w1 = --------------
μ Pr 1 / 3
+ κ Cp μ
δ w2 = 10 10 ---------- Pr = -----------
Pr T λ
Pr T κ
β = 15Pr 2 / 3 – ---------- 1 + ln 1000 ----------
2κ Pr T
λ is the thermal conductivity, and κ is the von Karman constant equal to 0.41.
The distance between the computational fluid domain and the wall, δw, is always hw/2
for automatic wall treatment where hw is the height of the mesh cell adjacent to the
wall. hw/2 is almost always very small compared to any geometrical quantity of interest,
at least if a boundary layer mesh is used. For wall function, δw is at least hw/2 and can
be bigger if necessary to keep δw+ higher than 11.06. The computational results should
be checked so that the distance between the computational fluid domain and the wall,
δw, is everywhere small compared to any geometrical quantity of interest. The distance
δw is available for evaluation on boundaries.
+
[ H0 ]- = 0 (4-145)
• If direction is Along normal vector, the outlet temperature Tavg is defined by:
down ( u ⋅ nρC p ) dS T avg =
down ( u ⋅ nρCp T ) dS
Γ Γ
• If the direction is opposite to normal vector, the outlet temperature Tavg is defined
by:
up ( u ⋅ nρC p ) dS T avg =
up ( u ⋅ nρCp T ) dS
Γ Γ
See Theory for the Single-Phase Flow Interfaces and Theory for the Turbulent Flow
Interfaces in the CFD Module User’s Guide for a description of the theory related to
laminar and turbulent single-phase flow interfaces.
The density and viscosity of air might also vary due to temperature variations. See
Theory for the Nonisothermal Flow and Conjugate Heat Transfer Interfaces for details
on this other kind of dependency.
The Moisture Flow interface contains the fully compressible formulation of the
continuity and momentum equations. For laminar flow they read:
where
It also solves the equation for moisture transport in air, given in Equation 4-128 by
∂c v
M v --------- + M v u ⋅ ∇c v + ∇ ⋅ g w = G
∂t
where
When the spatial and temporal variations of the vapor concentration induce significant
variations of the moist air density, the moisture content variation is expressed through
the transport of vapor mass fraction, ωv, defined as:
Mv cv
ω v = --------------
ρg
with
g w = – ρ g D ∇ω v
ρ g u Stefan = n ⋅ g w
The Moisture Flow interface allows to use the Stefan velocity as a leakage velocity at
walls in the fluid flow equations.
ρg cv
c˜v = -----------
ρg
where the bar denotes the usual Reynolds average. The full field can be decomposed
as the sum of the Favre average and the Favre fluctuation:
c v = c˜v + c v ″
∂c v ″
u j ″ -----------
∂x j
This flux is modeled using a gradient based assumption, where the additional transport
is related to the turbulent viscosity νT through the turbulent Schmidt number ScT.
The turbulent diffusivity DT is defined by
νT
D T = ----------
Sc T
1 0.3μ T Sc μ T Sc 2 μ –1
Sc T = ----------------- + ----------------------- – 0.3 ------------- 1 – exp – ----------------------------------------- (4-148)
2Sc T∞ μ Sc T∞ μ 0.3μ Sc Sc
T T∞
where the Schmidt number at infinity is ScT∞ = 0.85 and the Schmidt number is
defined as
ν
Sc = ----
D
Assuming that the turbulent heat and moisture transfer in the near-wall region are
analogous, the same type of wall functions used for the temperature (Ref. 31) is also
applicable for the moisture transport. The moisture transfer wall function is formulated
as a function of the turbulent Schmidt number, instead of the corresponding Prandtl
number.
The moisture flux at the lift-off position between the air with vapor mass fraction ωv,a
and a wall with vapor mass fraction ωv,w, is:
ρ g u τ ( ω v, w – ω v, a )
g wf = -----------------------------------------------
-
ωw +
(Ref. 31):
where in turn
+ δ w ρ C μ1 / 2 k + 10 -
δ w = ------------------------------- δ w1 = -------------
μ Sc 1 / 3
+ κ μ
δ w2 = 10 10 ---------- Sc = --------
Sc T ρD
Sc T κ
β = 15Sc 2 / 3 – ---------- 1 + ln 1000 ----------
2κ Sc T
The distance between the computational fluid domain and the wall, δw, is always hw/2
for automatic wall treatment where hw is the height of the mesh cell adjacent to the
wall. hw/2 is almost always very small compared to any geometrical quantity of interest,
at least if a boundary layer mesh is used. For wall function, δw is at least hw/2 and can
be bigger if necessary to keep δw+ higher than 11.06. The computational results should
be checked so that the distance between the computational fluid domain and the wall,
δw, is everywhere small compared to any geometrical quantity of interest. The distance
δw is available for evaluation on boundaries.
+
[ cv ]- = 0 (4-149)
• If direction is Along normal vector, the outlet vapor concentration cavg is defined
by:
down ( ρ g u ⋅ n ) dS c avg =
down ( ρg u ⋅ ncv ) dS
Γ Γ
• If the direction is opposite to normal vector, the outlet vapor concentration cavg is
defined by:
up ( ρ g u ⋅ n ) dS c avg =
up ( ρg u ⋅ ncv ) dS
Γ Γ
Historically, the thermoelectric effect is known by three different names, reflecting its
discovery in experiments by Seebeck, Peltier, and Thomson. The Seebeck effect is the
conversion of temperature differences into electricity, the Peltier effect is the
conversion of electricity to temperature differences, and the Thomson effect is heat
produced by the product of current density and temperature gradients. These effects
are thermodynamically related by the Thomson relations:
P = ST
where P is the Peltier coefficient (SI unit: V), S is the Seebeck coefficient
(SI unit: V/K), T is the temperature (SI unit: K), and μTh is the Thomson coefficient
(SI unit: V/K). These relations show that all coefficients can be considered different
descriptions of one and the same quantity. The COMSOL formulation primarily uses
the Seebeck coefficient. The Peltier coefficient is also used as an intermediate variable,
but the Thomson coefficient is not used.
When simulating the thermoelectric effect, the following fluxes are the quantities of
interest:
Thermoelectric efficiency is measured by the figure of merit Z (SI unit: 1/K), defined
as:
S2 σ
Z = ----------
k
Some other quantities of relevance are the electric field E and the Joule heat source Q:
E = – ∇V
Q = J⋅E
∂T
ρC p +∇⋅q = Q
∂t
∇ ⋅ J = Qj
where ρ is the density, Cp the heat capacity, and Qj is the current source.
Seebeck Effect
The Seebeck effect is described as the conversion of temperature gradient into electric
current. The contribution of the Seebeck effect is defined as a current contribution
J Se = – σS∇T
Peltier Effect
The Peltier effect is described as the conversion of t electric current in heat source or
sink. It is defined as an heat source contribution
Q Pe = – P∇ ⋅ J
Thomson Effect
The Thomson effect defines the heat source induced by a current in presence of a
temperature gradient in thermoelectric material. The heat source is defined by
Q Th = – μ Th J ⋅ ∇T
This contribution is obtained again by developing the divergence of the q term in the
heat equation when q is defined following Equation 4-150. This time consider the
term −TJ ⋅ ∇S. Assuming that S is function of T, then:
dS
– TJ ⋅ ∇S = – T J ⋅ ∇T = – μ Th J ⋅ ∇T
dT
∂T
( ρC p ) eff ------- + ∇ ⋅ ( – k eff ∇T – L v δ p ∇( φ w p sat ) ) = Q (4-152)
∂t
∂φ w
ξ ---------- + ∇ ⋅ ( – ξ D w ∇φ w – δ p ∇( φ w p sat ) ) = G (4-153)
∂t
See Theory for Heat Transfer in Building Materials and Theory for Moisture Transport
in Building Materials for more details about the notations.
The heat transfer and moisture transport equations are coupled through:
• the definition of the effective thermal properties accounting for moisture content,
• the definition of the latent heat source from the evaporation source in domains,
• the definition of the vapor saturation in function of the temperature.
Theory for the Moist Air Version of the Heat and Moisture Transport
Interface
The Moist Air version of the Heat and Moisture multiphysics coupling implements the
following equations in domains for heat and moisture transport:
∂T
ρ g C p ------- + ρ g C p u ⋅ ∇T + ∇ ⋅ q = Q (4-154)
∂t
q = – k∇T (4-155)
Q = – ( C p, v – C p, a )g w (4-156)
∂c v
M v -------- + M v u ⋅ ∇c v + ∇ ⋅ g w = G (4-157)
∂t
g w = – M v D ∇c v (4-158)
c v = φ w c sat (4-159)
The heat transfer and moisture transport equations are coupled through:
• the definition of the effective thermal properties accounting for vapor concentration
in moist air,
• the definition of the diffusive flux of enthalpy (Equation 4-156),
• the definition of the latent heat source from the evaporation source on surfaces
(Equation 4-160),
• the definition of the vapor saturation in function of the temperature
(Equation 4-159).
Theory for the Moist Porous Media Version of the Heat and Moisture
Transport Interface
The Moist Porous Media version of the Heat and Moisture multiphysics coupling
implements the following equations in domains for heat and moisture transport:
∂T
( ρC p ) eff ------- + ( ρ g C p, g u g + ρ l C p, l u l ) ⋅ ∇T + ∇ ⋅ q = Q + Q evap (4-161)
∂t
q = – k eff ∇T (4-162)
Q = – [ ( C p, v – C p, a )g w + C p, l g lc ] (4-163)
∂w ( φ w )
-------------------- + ρ g u g ⋅ ∇ω v + ∇ ⋅ g w + u l ⋅ ∇ρ l + ∇ ⋅ g lc = G (4-165)
∂t
g w = – ρ g D eff ∇ω v (4-166)
∂ [ εp ρg ωv ( 1 – sl ) ]
G evap = ---------------------------------------------- + ρ g u g ⋅ ∇ω v + ∇ ⋅ g w (4-167)
∂t
The heat transfer and moisture transport equations are coupled through:
• the definition of the effective thermal properties accounting for moisture content in
liquid and gas phases,
• the liquid water saturation and velocity field used to defined the convective heat
flux,
• the diffusive flux of enthalpy in moist air and the liquid capillary flux
(Equation 4-163)
• the definition of the latent heat source from the evaporation source in domains
(Equation 4-164),
• the definition of the vapor saturation in function of the temperature.
They all have in common the multiphysics coupling feature Electromagnetic Heating,
which adds weak contributions due to resistive losses in the domains and boundaries,
and shares the temperature variable with the electromagnetics interfaces. The
underlying theory can be found in the AC/DC Module User’s Guide, RF Module
User’s Guide, and Wave Optics Module User’s Guide.
Both require the Structural Mechanics Module. They have in common the use of the
Thermal Expansion multiphysics coupling that models temperature dependence of the
strain tensor and thermoelastic damping. For more details about the underlying
theory, see the Structural Mechanics Module User’s Guide.
The heat fluxes at the upside and downside boundaries depend on the temperature
difference according to the relations:
– n d ⋅ ( – k d ∇T d ) = – h ( T u – T d ) + rQ b
– n u ⋅ ( – k u ∇T u ) = – h ( T d – T u ) + ( 1 – r )Q b
masp,u asp,u
Y
The joint conductance h has three contributions: the constriction conductance, hc,
from the contact spots, the gap conductance, hg, due to the fluid at the interstitial
space, and the radiative conductance, hr:
h = hc + hg + hr
SURFACE ASPERITIES
The microscopic surface asperities are characterized by the average height σu, asp and
σd, asp and the average slope mu, asp and md, asp. The RMS values σasp and masp are
(4.16 in Ref. 1):
σ asp = σ u2, asp + σ d2, asp m asp = m u2, asp + m d2, asp
m asp p 0.95
h c = 1.25k contact ------------- -------
σ asp H c
Here, Hc is the microhardness of the softer material, p is the contact pressure, and
kcontact is the harmonic mean of the contacting surface conductivities:
2k u k d
k contact = -------------------
ku + kd
-----------------------------------
1
p- p ( 1 + 0.071c 2 )
------ = -----------------------------------------------------
-
Hc σ asp
c2
c
1 1.62 -----------
- m
σ 0 asp
The coefficients c1 and c2 are the Vickers correlation coefficient and size index,
respectively, and σ0 is equal to 1 µm. For materials with Brinell hardness between 1.30
and 7.60 GPa, c1 and c2 are given by the correlation below (4.16.1 in Ref. 1):
c1 HB HB 2 HB 3
------- = 4.0 – 5.77 -------- + 4.0 -------- – 0.61 --------
H0 H0 H0 H0
m asp 0.94
h c = 1.54k contact ------------- ------------------------
2p
σ asp mE contact
Here, Econtact is an effective Young’s modulus for the contact interface, satisfying
(4.16.3 in Ref. 1):
1 - 1 – ν u2 1 – ν d2
------------------ = --------------- + ---------------
E contact Eu Ed
where Eu and Ed are the Young’s moduli of the two contacting surfaces and νu and νd
are the Poisson’s ratios.
GAP CONDUCTANCE
The gap conductance due to interstitial fluid cannot be neglected for high fluid
thermal conductivity or high contact pressure. The parallel-plate gap gas correlation
assumes that the interstitial fluid is a gas and defines hg by:
kg
h g = ------------------
Y + Mg
Here kg is the gas conductivity, Y denotes the mean separation thickness (see
Figure 4-23), and Mg is the gas parameter equal to:
kB Tg
M g = αβΛ Λ = ------------------------
2
-
2πD p g
The mean separation thickness, Y, is a function of the contact pressure, p. For low
values of p near 0 Pa, Y goes to infinity since no contact occur. For high values of p —
greater than Hc ⁄2 in the Cooper-Mikic-Yovanovich model and greater than Hc ⁄4 in
the Mikic elastic model — Y reduces to 0 meaning that the contact is considered as
perfect.
RADIATIVE CONDUCTANCE
At high temperatures, above 600°C, radiative conductance needs to be considered.
The gray-diffuse parallel plate model provides the following formula for hr:
εu εd
h r = ----------------------------------- σ ( T u3 + T u2 T d + T u T d2 + T d3 )
εu + εd – εu εd
εu εd
h r ( T u – T d ) = ----------------------------------- σ ( T u4 – T d4 )
εu + εd – εu εd
εu εd
h r ( T d – T u ) = ----------------------------------- σ ( T d4 – T u4 )
εu + εd – εu εd
THERMAL FRICTION
The friction heat, Qb, is partitioned into rQb and (1 − r)Qb at the contact interface. If
the two bodies are identical, r and (1 − r) would be 0.5 so that half of the friction heat
goes to each surface. However, in the general case where the two bodies are made of
different materials, the partition rate might not be 0.5. The Charron’s relation (Ref. 2)
defines r as:
1 ρ u C p, u k u
r = --------------- ξd = -------------------------
1 + ξd ρ d C p, d k d
1 ρ d C p, d k d
( 1 – r ) = --------------- ξu = -------------------------
1 + ξu ρ u C p, u k u
Thermal Contact
REFERENCES
1. A. Bejan and others, Heat Transfer Handbook, John Wiley & Sons, 2003.
qup
qdown
The reduced geometry does not include all the boundaries of the original 3D
geometry. For example, the reduced geometry does not represent the upside and
downside surfaces of the plate in Figure 4-24 as boundaries.
Equation Formulation
2D GEOMETRIES
In 2D geometries, the temperature is assumed to be constant in the out-of-plane
direction (z direction with default spatial coordinate names). The equation for heat
transfer in solids, Equation 4-16, and in fluids, Equation 4-18, are replaced by:
∂T
d z ρC p ------- + ρC p d z u ⋅ ∇T + ∇ ⋅ q = d z Q + q 0 (4-169)
∂t
Here dz is the thickness of the domain in the out-of-plane direction. Here, the
conductive heat flux, q, becomes
q = – d z k∇T
1D AXISYMMETRIC GEOMETRIES
In 1D axisymmetric geometries, the temperature is assumed to be constant in the
out-of-plane direction (z direction with default spatial coordinate names) in addition
to the axisymmetry ( φ -coordinate with default spatial coordinate names). The
equation for heat transfer in solids, Equation 6-15 is replaced by
∂T
( 2πrd z )ρC p ------- + ∇ ⋅ q = ( 2πrd z )Q + q 0 (4-170)
∂t
where dz is the thickness of the domain in the z direction. The equation for heat
transfer in fluids, Equation 6-5, is replaced by
∂T
( 2πrd z )ρC p ------- + ( 2πrd z )ρC p u ⋅ ∇T + ∇ ⋅ q = ( 2πrd z )Q + q 0 (4-171)
∂t
q = – ( 2πrd z )k∇T
1D GEOMETRIES
In 1D geometries, the temperature is assumed to be constant in the radial direction.
The equation for heat transfer in solids, Equation 6-15 is replaced by
∂T
A c ρC p ------- + ∇ ⋅ q = A c Q + q 0 (4-172)
∂t
where Ac is the cross section of the domain in the plane perpendicular to the 1D
geometry. The equation for heat transfer in fluids, Equation 6-5, is replaced by
∂T
A c ρC p ------- + A c ρC p u ⋅ ∇T + ∇ ⋅ q = A c Q + q 0 (4-173)
∂t
Out-of-plane flux conditions would apply to the exterior boundaries of the domain if
the 1D geometry was seen as a cylinder. With the geometry reduction process, this
heat flux condition is mathematically expressed using the cross section perimeter, Pc,
as in:
q 0 = P c q 0, z
where q0, z is the heat flux density distributed along the cross section perimeter.
The second approach is the most accurate if the geometry or the external flow is
complicated. The Heat Transfer Module includes the Conjugate Heat Transfer
predefined multiphysics coupling and the CFD Module includes the Nonisothermal
Flow predefined multiphysics coupling for this purpose. However, such a simulation
can become costly, both in terms of computational time and memory requirement.
The first method is simple, yet powerful and efficient. The convective heat flux on the
boundaries in contact with the fluid is then modeled as being proportional to the
temperature difference across a fictitious thermal boundary layer. Mathematically, the
heat flux is described by the equation
– n ⋅ q = h ( T ext – T )
where h is a heat transfer coefficient and Text the temperature of the external fluid far
from the boundary.
The main difficulty in using heat transfer coefficients is in calculating or specifying the
appropriate value of the h coefficient. That coefficient depends on the fluid’s material
properties, and the surface temperature — and, for forced convection, also on the
fluid’s flow rate. In addition, the geometrical configuration affects the coefficient. The
Heat Transfer interface has built-in functions for the heat transfer coefficients. For
most engineering purposes, the use of such coefficients is an accurate and numerically
efficient modeling approach.
In this section:
External
Internal
Laminar Flow
Turbulent Flow
The difference between natural and forced convection is that in the forced convection
an external force such as a fan creates the flow. In natural convection, buoyancy forces
induced by temperature differences together with the thermal expansion of the fluid
drive the flow.
Heat transfer books generally contain a large set of empirical and theoretical
correlations for h coefficients. This module includes a subset of them. The expressions
are based on the following set of dimensionless numbers:
where
Further, GrL refers to the Grashof number, which is the squared ratio of the viscous
time scale to the buoyancy time scale multiplied by the Reynolds number.
3
gρ ( ρ ext – ρ s )L
Gr L = -----------------------------------------
2
μ
where g is the acceleration of gravity, ρs denotes the density of the hot surface, ρext
equals the free stream density, L is the length scale, μ represents the fluid’s dynamic
viscosity, and ρ its density.
When the density depends on temperature only, in dry air for example, the Grashof
number is well approximated by using the fluid’s coefficient of thermal expansion αp
(SI unit: 1/K):
3
gα p ( T s – T ext )L
Gr L = ----------------------------------------------
2
(μ ⁄ ρ)
where g is the acceleration of gravity, Ts denotes the temperature of the hot surface,
Text equals the free stream temperature, L is the length scale, μ represents the fluid’s
dynamic viscosity, and ρ its density.
1 ∂ρ
α p = – --- -------
ρ ∂T p
1
α p = ----
T
The transition from laminar to turbulent flow occurs at a Gr value of 109; the flow is
turbulent for larger values.
The Rayleigh number, Ra, is another indicator of the regime. It is similar to the
Grashof number except that it accounts for the thermal diffusivity: Ra = Pr Gr. A small
value of the Ra number indicates that the conduction dominates. It such case using
heat transfer coefficients to model convective heat transfer is not relevant. Instead,
modeling the fluid as immobile is likely to be accurate.
VERTICAL WALL
Figure 4-26: Schematic representation of geometry and parameters for the heat transfer
coefficient correlation applied to natural convection on a vertical wall.
The following correlations correspond to equations 9.26 and 9.27 in Ref. 21:
1⁄4
--- k- 0.67Ra L
L 0.68 + ---------------------------------------------------------
- if Ra L ≤ 10 9
9 / 16 4 / 9
1 + -------------------
0.492k
μC p
h = (4-174)
1⁄6 2
k 0.387Ra L
---- 0.825 + ------------------------------------------------------------
9 / 16 8 / 27
- if Ra L > 10 9
L 1 + -------------------
0.492k
μC p
2
gα p ρ C p T – T ext L 3
Ra L = -------------------------------------------------------- (4-175)
kμ
gρC p ρ ext – ρ s L 3
Ra L = ----------------------------------------------- (4-176)
kμ
All material properties are evaluated at (T + Text) ⁄ 2, except ρs, which is evaluated at
the wall temperature, T, and g is the acceleration of gravity equal to 9.81 m/s2. This
INCLINED WALL
Figure 4-27: Schematic representation of geometry and parameters for the heat transfer
coefficient correlation applied to natural convection on an inclined wall.
The following correlations correspond to equations 9.26 and 9.27 in Ref. 21 (the
same as for a vertical wall):
--- k- 0.67 ( cos ϕRa L ) 1 / 4
L 0.68 + ---------------------------------------------------------- if Ra L ≤ 10 9
9 / 16 4 / 9
1 + 0.492k
-------------------
μC p
h = (4-177)
1⁄6 2
k 0.387Ra L
---- 0.825 + ------------------------------------------------------------
/ /
- if Ra L > 10 9
L 1 + 0.492k -------------------
9 16
8 27
μC p
where the length of the wall, L, is a correlation input and ϕ is the tilt angle (the angle
between the wall and the vertical direction; ϕ = 0 for vertical walls). These
correlations are valid for −60° < ϕ < 60° and 104≤ RaL ≤ 1013.
The definition of the Raleigh number, RaL, is analogous to the one for vertical walls
and is given by the following:
2
gα p ρ C p T – T ext L 3
Ra L = -------------------------------------------------------- (4-178)
kμ
For turbulent flow, 1 is used instead of cos ϕ in the expression for h, because this gives
better accuracy (see Ref. 41).
According to Ref. 21, correlations for inclined walls are only satisfactory
for the top side of a cold plate or the down face of a hot plate. Hence,
these correlations are not recommended for the bottom side of a cold face
and for the top side of a hot plate.
Ra L > 10 9
All material properties are evaluated at (T + Text) ⁄ 2, except ρswhich is evaluated at the
wall temperature, T. The laminar-turbulent transition at RaL = 109 is handled by the
use of a smoothed Heaviside function with a transition of size 108, which corresponds
to 10% of the threshold value.
Figure 4-28: Schematic representation of geometry and parameters for the heat transfer
coefficient correlation applied to natural convection on the top surface of an horizontal
plate.
k-
---
1⁄4
0.54Ra L if 10 4 ≤ Ra L ≤ 10 7
L
h = (4-180)
k-
---
1⁄3
0.15Ra L if 10 7 ≤ Ra L ≤ 10 11
L
k 1⁄4
h = ---- 0.27Ra L if 10 5 ≤ Ra L ≤ 10 10 (4-181)
L
RaL is given by Equation 4-175 or Equation 4-176, and L, the characteristic length
(defined as area/perimeter, see Ref. 41) is a correlation input. The material data are
evaluated at (T + Text) ⁄ 2, except ρswhich is evaluated at the wall temperature, T.
When the density depends only on temperature, the conditions ρ < ρext and ρ ≥ ρext
can be replaced by T > Text and T ≤ Text respectively.
Figure 4-29: Schematic representation of geometry and parameters for the heat transfer
coefficient correlation applied to natural convection on the bottom surface of an
horizontal plate.
Figure 4-30: Schematic representation of geometry and parameters for the heat transfer
coefficient correlation applied to natural convection on a long horizontal cylinder.
The following correlations correspond to equations 9.34 in Ref. 21. It is validated for
RaD ≤ 1012.
1⁄6 2
k 0.387Ra D
h = ---- 0.6 + ---------------------------------------------------------
⁄ ⁄
- (4-182)
D 1 + 0.559
9 16 8 27
---------------
Pr
2
gα p ρ C p T – T ext D 3
Ra D = ---------------------------------------------------------
kμ
gρC p ρ ext – ρ s D 3
Ra D = ------------------------------------------------
kμ
The material data are evaluated at (T + Text) ⁄ 2, except ρswhich is evaluated at the wall
temperature, T.
Figure 4-31: Schematic representation of geometry and parameters for the heat transfer
coefficient correlation applied to natural convection on a sphere.
The following correlations correspond to equation 9.35 in Ref. 21. It is validated for
RaD ≤ 1011 and Pr ≥ 0.7.
1⁄4
k 0.589Ra D
h = ---- 2 + ------------------------------------------------------
⁄ ⁄
- (4-183)
D 1 + 0.469
9 16 4 9
---------------
Pr
2
gα p ρ C p T – T ext D 3
Ra D = ---------------------------------------------------------
kμ
gρC p ρ ext – ρ s D 3
Ra D = ------------------------------------------------
kμ
The material data are evaluated at (T + Text) ⁄ 2, except ρswhich is evaluated at the wall
temperature, T.
Figure 4-32: Schematic representation of geometry and parameters for the heat transfer
coefficient correlation applied to natural convection on a vertical thin cylinder.
The following correlation corresponds to equation 7.83 in Ref. 41. It evaluates the
heat transfer coefficient on the lateral surface of the thin cylinder. The horizontal disks
(top and bottom) should be treated as horizontal plates.
7Ra H Pr 1⁄4
4 ( 272 + 315Pr )H
h = ----- --- ----------------------------------- + -----------------------------------------------
k 4
H 3 5 ( 20 + 21Pr ) 35 ( 64 + 63Pr )D
where D is the cylinder diameter, H is the cylinder height, and RaH is given by
3
gα p T – T ext H
Ra H = -------------------------------------------
kμ
gρC p ρ ext – ρ s H 3
Ra H = ------------------------------------------------
kμ
The material data are evaluated at (T + Text) ⁄ 2, except ρswhich is evaluated at the wall
temperature, T.
If the thermal boundary layer thickness δT is much smaller than the cylinder diameter
D, the curvature does not play any significant role, and vertical wall correlations should
be used. In practice, the (δT < D) criterion requires
1
– ---
D- 4
---- > Ra H
H
Figure 4-33: Schematic representation of geometry and parameters for the heat transfer
coefficient correlation applied to natural convection in a parallel-plate narrow chimney.
k 1
h = ----- ------ Ra L (4-184)
H 24
where the plate distance, L, and the chimney height, H, are correlation inputs
(equation 7.96 in Ref. 41). RaL is given by Equation 4-175 or Equation 4-176. The
material data are evaluated at (T + Text) ⁄ 2.
Figure 4-34: Schematic representation of geometry and parameters for the heat transfer
coefficient correlation applied to natural convection in a circular narrow chimney.
k 1
h = ----- ---------- Ra D
H 128
Figure 4-35: Schematic representation of geometry and parameters for the averaged heat
transfer coefficient correlation applied to forced convection on an horizontal plate.
1⁄2
0.3387Pr 1 / 3 Re L
k- ------------------------------------------------------
2 --- - if Re L ≤ 5 ⋅ 10 5
L 0.0468- 2 / 3 1 / 4
------------------
h = 1 + (4-185)
Pr
2 ---- Pr 1 / 3 ( 0.037Re 4 ⁄ 5 – 871 ) if Re > 5 ⋅ 10 5
k
L L L
where Pr = μCp ⁄ k and ReL = ρUL ⁄ μ. The plate length, L, and the exterior velocity,
U, are correlation inputs. The material data are evaluated at (T + Text) ⁄ 2.
Figure 4-36: Schematic representation of geometry and parameters for the local heat
transfer coefficient correlation applied to forced convection on an horizontal plate.
k 1⁄2
--- 0.332Pr 1 / 3 Re x if Re x ≤ 5 ⋅ 10 5
x
h = (4-186)
k 4⁄5
--- 0.0296Pr 1 / 3 Re x if Re x > 5 ⋅ 10 5
x
where Pr = μCp ⁄ k and Rex = ρUx ⁄ μ. The correlation inputs are x, the position along
the plate, and U, the exterior velocity. The material data are evaluated at (T + Text) ⁄ 2.
To avoid division by zero when the position along the plate is located at the origin
point (x = 0), the implementation replaces k ⁄ x by k ⁄ max(x, √ε) where ε is the
floating-point relative accuracy.
Figure 4-37: Schematic representation of geometry and parameters for the averaged heat
transfer coefficient correlation applied to forced convection on a cylinder in cross-flow.
1⁄2
k 0.62Re D Pr 1 / 3 Re D 5 ⁄ 8 4 ⁄ 5
2 ⁄ 3 1 ⁄ 4
- 1 + --------------------
h = ---- 0.3 + ---------------------------------------------- 282000 if Re D Pr ≥ 0.2 (4-187)
D 1 + 0.4
--------
-
Pr
where Pr = μCp ⁄ k and ReD = ρUD ⁄ μ. The cylinder diameter, D, and the exterior
velocity, U, are correlation inputs. The material data are evaluated at (T + Text) ⁄ 2.
SPHERE
Figure 4-38: Schematic representation of geometry and parameters for the averaged heat
transfer coefficient correlation applied to forced convection on a sphere.
0.4 μ 1 ⁄ 4
h = ---- 2 + ( 0.4Re D + 0.06Re D )Pr -----
k 1⁄2 2⁄3
(4-188)
D μs
where Pr = μCp ⁄ k and ReD = ρUD ⁄ μ. The sphere diameter, D, and the exterior
velocity, U, are correlation inputs. All material data are evaluated at Text except μs,
which is evaluated at the wall temperature, T.