Overview

• Introduction
• Code Development
• Theory and Code Implementation
• View Factor Benchmarks (6)
• Heat Transfer Benchmarks (3)
• Future Work

Introduction

Motivation

• Radiation heat transfer occurs all over the place
• It is a prominent mode of heat transfer when high temperatures are present
• Can be critical for determining realistic temperature distributions (e.g. peak temperature analysis)

Nuclear Applications

• Gap between the fuel and the cladding in LWR fuel pins
• Reactor vessel of VHTRs to the cooling panel of the RCCS
• Corium and other high temperature surfaces in severe accidents
• Surface of heat pipes in terrestrial and space, nuclear reactors

Non-Nuclear Applications

• Continuous annealing furnaces used for recrystallization annealing and surface treatment in cold-rolled strip manufacturing
• Thermal control and thermal protection systems used by spacecraft
• Industrial furnaces
• Etc.

Challenges

• The computation of view factors can be done analytically for only the simplest geometries
• Numerical methods are required for arbitrary geometries
• Simulation of heat diffusion with radiation heat transfer

Goals

To develop an extensible code that can:

1. Calculate view factors for arbitrary geometries
2. Simulate heat diffusion within those geometries while fully resolving all nonlinearities
3. Have the capability of being fully-coupled to other codes that simulate other physical phenomena (e.g. material deformation, neutronics, fluid-dynamics, etc.)

Code Development

History (1/2)

• Began as a weekend side-project of a side-project in the Summer of 2017
• Standalone python code to compute view factors via gaussian and adaptive quadrature rules: thermal_radiation
• Became curious with computational geometry and stumbled across CGAL
• Leveraged CGAL into a MOOSE app and implemented view factor functionality

History (2/2)

• Focused on school and other work
• Revived code in Phoenix in January 2020
• Added blackbody boundary condition capabilities
• Conducted more extensive benchmarking and testing

Radiation Heat Transfer Basics

• The exchange of thermal energy between surfaces via photons
• All surfaces emit photons as a result of oscillations or transitions of the electrons that constitute matter
• The emitted photons encompass a range of wavelengths, and they are distributed in a variety of directions
• Blackbody (idealized opaque, non-reflective) surfaces will be assumed for this thesis

Heat Conduction with Radiation Boundary Conditions

\begin{equation} \rho c_p \frac{\partial T}{\partial t} - \nabla \cdot \left( k \nabla T \right) = q''', \:\:\:\: \forall \vec{r} \in \Omega \end{equation}

\begin{equation} q''_\text{rad} = -k\frac{d T}{d \hat{n}}, \:\:\:\: \forall \vec{r} \in \partial \Omega_\text{rad} \end{equation}

Radiation Boundary

• $\mathcal{G}$ : irradiation incident on the surface
• $\mathcal{E}$ : emissive power of the surface
• $r$ : reflectivity of the surface

Radiation Boundary (cont.)

\begin{equation} r + \alpha + \tau = 1 \end{equation}

• $\alpha$ : absorptivity of the surface
• $\tau$ : transmissivity of the surface

For black, opaque surfaces $r = 0$ and $\tau = 0$.

Blackbody Simplification

\begin{equation} q''_\text{rad}(\vec{r}) = \mathcal{E}^\text{b}(\vec{r}) - \mathcal{G}^\text{b}(\vec{r}), \:\:\:\: \forall \vec{r} \in \partial \Omega_\text{rad} \end{equation}

where, \begin{equation} \mathcal{E}^\text{b}(\vec{r}) = \sigma T(\vec{r})^4 \end{equation}

What about $\mathcal{G}^\text{b}$?

View Factor Theory and Blackbody Irradiation

Infinitesimal Geometric Configuration

Infinitesimal View Factor

\begin{equation} dF_{e \to f} \equiv \frac{\text{diffuse energy leaving } dA_e \text{ directly toward and intercepted by } dA_f}{\text{total diffuse energy leaving } dA_e} \end{equation}

skipping a few details

\begin{equation} dF_{e \to f} = \sphericalangle(\vec{r}_e, \vec{r}_f) \frac{\cos \theta_e \cos \theta_f}{\pi \| \vec{s}_{e \to f} \|^2} dA_f \end{equation}

where $\sphericalangle(\vec{r}_e, \vec{r}_f)$ is the view function.

View Function

\begin{equation} \sphericalangle(\vec{r}_e, \vec{r}_f) = \begin{cases} 0 & \text{if the ray, $\vec{s}_{e \to f}$, is occluded by a surface} \\ 1 & \text{otherwise} \end{cases} \end{equation}

Note that $\sphericalangle(\vec{r}_e, \vec{r}_f) = \sphericalangle(\vec{r}_f, \vec{r}_e)$.

View Factor

\begin{equation} F_{e \to f} = \frac{1}{A_e} \int_{\partial \Omega_e} \int_{\partial \Omega_f} \sphericalangle(\vec{r}_e, \vec{r}_f) \frac{\cos\theta_e \cos\theta_f}{\pi \|\vec{s}_{e \to f}\|^2} dA_f dA_e \end{equation}

where $\partial \Omega_e$ is a finite surface, $e$, and the surface area of that surface is $A_e = | \partial \Omega_e |$.

View Factor Algebra

The reciprocity relationship: \begin{equation} A_e F_{e \to f} = A_f F_{f \to e} \end{equation}

The partial additive relationship: \begin{equation} F_{e \to F} = \sum_{\partial \Omega_f \in \partial \Omega_F} F_{e \to f} \end{equation}

The convex, partial-enclosure constraint: \begin{equation} \sum_{\partial \Omega_f \in \partial \Omega_F} F_{e \to f} \le 1 \end{equation}

The complete additively relationship: \begin{equation} F_{E \to F} = \frac{1}{A_E} \sum_{\partial \Omega_e \in \partial \Omega_E} \sum_{\partial \Omega_f \in \partial \Omega_F} A_e F_{e \to f} \end{equation}

Differential Amount of Diffuse Energy Intercepted

without going into too many details.

\begin{equation} \mathcal{G}^\text{b}(\vec{r}_e) = \sigma \int_{\partial \Omega_f} T(\vec{r}_f)^4 \sphericalangle(\vec{r}_e, \vec{r}_f) \frac{\cos \theta_e \cos \theta_f}{\pi \| \vec{s}_{e \to f} \|^2} dA_f. \end{equation}

• It is impractiacal to approximate $\mathcal{G}^\text{b}(\vec{r}_e)$
• Notice the similarities between the integral kernel for $F_{E \to F}$ and $\mathcal{G}^\text{b}(\vec{r}_e)$
• Perhaps we can exploit this similarity...

Surface-Area Averaging

Surface-Averaging Operator

Consider the element $\partial \Omega^h_e \in \partial \Omega^h_\text{rad}$. Similar to the previous section, $A_e = | \partial \Omega^h_e |$ (i.e. $A_e$ is the area of $\partial \Omega^h_e$)

\begin{equation} \langle \: \cdot \: \rangle_e \equiv \frac{\int_{\partial \Omega^h_e} \cdot \: dA}{A_e} \end{equation}

Total Energy Emitted from an Element via Blackbody Radiation

\begin{align} \mathcal{E}^\text{b}_e & \equiv \int_{\partial \Omega^h_e} \mathcal{E}^\text{b} (\vec{r}) dA \\ & = \int_{\partial \Omega^h_e} \sigma T(\vec{r})^4 dA \nonumber \\ & = \sigma \int_{\partial \Omega^h_e} T(\vec{r})^4 dA \nonumber \\ \mathcal{E}^\text{b}_e & = \sigma A_e \langle T^4 \rangle_e \end{align}

Irradiation of an Element via Blackbody Radiation

The differential, amount of diffuse energy that irradiates an element:

\begin{equation} \mathcal{G}^\text{b}_e(\vec{r}_e) \equiv \sigma \sum_{\partial \Omega_f \in \partial \Omega^h_\text{rad}} \int_{\partial \Omega_f} T(\vec{r}_f)^4 \sphericalangle(\vec{r}_e, \vec{r}_f) \frac{\cos \theta_e \cos \theta_f}{\pi \| \vec{s}_{e \to f} \|^2} dA_f \end{equation}

The total amount of diffuse energy that irradiates an element:

\begin{align} \mathcal{G}^\text{b}_e & \equiv \int_{\partial \Omega_e} \mathcal{G}^\text{b}_e(\vec{r}_e) dA_e \\ & = \sigma \sum_{\partial \Omega_f \in \partial \Omega^h_\text{rad}} \int_{\partial \Omega_f} T(\vec{r}_f)^4 \left( \int_{\partial \Omega_e} \sphericalangle(\vec{r}_e, \vec{r}_f) \frac{\cos \theta_e \cos \theta_f}{\pi \| \vec{s}_{e \to f} \|^2} dA_e \right) dA_f \end{align}

Constant Temperature Assumption

Assumption: $T(\vec{r}_e) = \langle T \rangle_e \:\: \forall \vec{r}_e \in \partial \Omega_e, \: \forall \partial \Omega_e \in \partial \Omega_\text{rad}^h \implies \langle T^4 \rangle_e = \langle T \rangle^4_e$.

Total amount of diffuse energy that irradiates the element: \begin{align} \mathcal{G}^\text{b}_e & = \sigma \sum_{\partial \Omega_f \in \partial \Omega^h_\text{rad}} \langle T \rangle^4_f \int_{\partial \Omega_f}\int_{\partial \Omega_e} \sphericalangle(\vec{r}_e, \vec{r}_f) \frac{\cos \theta_e \cos \theta_f}{\pi \| \vec{s}_{e \to f} \|^2} dA_e dA_f \nonumber \\ & = \sigma \sum_{\partial \Omega_f \in \partial \Omega^h_\text{rad}} \langle T^4 \rangle_f \int_{\partial \Omega_f}\int_{\partial \Omega_e} \sphericalangle(\vec{r}_e, \vec{r}_f) \frac{\cos \theta_e \cos \theta_f}{\pi \| \vec{s}_{e \to f} \|^2} dA_e dA_f \nonumber \\ & = \sigma \sum_{\partial \Omega_f \in \partial \Omega^h_\text{rad}} \langle T^4 \rangle_f A_f F_{f \to e} \nonumber \\ \mathcal{G}^\text{b}_e & = \sum_{\partial \Omega_f \in \partial \Omega^h_\text{rad}} \mathcal{E}^\text{b}_f F_{f \to e} \end{align}

Surface-Area Averaged Blackbody Radiation Boundary Condition

\begin{equation} \left( -k \nabla T (\vec{r}) \right) \cdot \hat{n} = \mathcal{E}^\text{b}_e - \mathcal{G}^\text{b}_e, \:\:\:\: \forall \vec{r} \in \partial \Omega_e \end{equation}

Finite Element Modeling

MOOSE

• The Multiphysics Object-Oriented Simulation Environment
• Language: C++
• License: LGPL
• ~20k commits with ~110 contributors
• Fully-coupled, fully-implicit multiphysics solver.
• Many complementary applications are based on MOOSE (RELAP-7, ​BISON, Rattlesnake, etc.)

Finite Element Prelude

Let $\phi_j(\vec{r})$ be that basis functions. The approximate solution variable is represented as: \begin{equation} T_h(\vec{r}) = \sum_{j = 1}^{N_\text{dofs}} T_j \phi_j(\vec{r}) \end{equation} It follows that: \begin{equation} \nabla T_h(\vec{r}) = \sum_{j = 1}^{N_\text{dofs}} T_j \nabla \phi_j(\vec{r}) \end{equation}

In the remainder of this presentation $h$ will be dropped from $T_h$.

Weak Form of the Heat Diffusion Equation

Multiplying the heat diffusion equation by the test fuction, $\psi_i(\vec{r})$, integrating over that spatial domain, $\Omega$, and appropriately applying Gauss’ divergence theorem:

\begin{equation} \int_\Omega \rho c_p \frac{\partial T}{\partial t} \psi_i \: dV + \int_\Omega \left( k \nabla T \right) \cdot \nabla \psi_i \: dV - \int_{\partial \Omega} \psi_i \left( k \nabla T \right) \cdot \hat{n} \: dA - \int_\Omega q''' \psi_i \: dV = 0 \end{equation}

Weak Form of the Heat Diffusion Equation (cont.)

If $\psi_i$ has local support restricted to an element in $\Omega_e \subset \Omega$:

\begin{equation} \int_{\Omega_e} \rho c_p \frac{\partial T}{\partial t} \psi_i \: dV + \int_{\Omega_e} \left( k \nabla T \right) \cdot \nabla \psi_i \: dV + \int_{\partial \Omega \cap \partial \Omega_e} - \psi_i \left( k \nabla T \right) \cdot \hat{n} \: dA + \int_{\Omega_e} - q''' \psi_i \: dV = 0 \end{equation}

Each of the terms in the above equation is a so called kernel in MOOSE.

\begin{equation} R_i(\vec{T}) = K_i^\text{time}(\vec{T}) + K_i^\text{cond}(\vec{T}) + K_i^\text{bndry}(\vec{T}) + K_i^\text{src}(\vec{T}) \end{equation}

• For radiation boundary conditions, $K_i^\text{bndry}(\vec{T})$, will be nonlinear.
• Implying $R_i(\vec{T})$ is nonlinear.
• How is $R_i(\vec{T})$ solved if it is nonlinear?

Preconditioned Jacobian-free Newton-Krylov (PJFNK)

Consider the residual statement for a general nonlinear system of equation: \begin{equation} \vec{F}(\vec{x}) = 0 \end{equation} The elements of the jacobian matrix for the system are: \begin{equation} J(\vec{x})_{i, j} = \frac{\partial (F)_i}{\partial (x)_j}(\vec{x}) \end{equation} Newton methods (in essence) consists of iteratively updating the iterate solution vector $\vec{x}^n$ in the following fashion: \begin{equation} \vec{x}^{n+1} = \vec{x}^n + \delta \vec{x}^n \end{equation} The newton direction is determined by solving the following linear system: \begin{equation} J(\vec{x}^n) \delta \vec{x}^n = - \vec{F}(\vec{x}^n) \end{equation}

View Factor Computation

Possible Methods

• Direct integration
  • Surface integration. Analytical or numerical.
  • Contour integration. Analytical or numerical.
• Statistical determination (Monte Carlo)
• Special methods
  • View factor algebra.
  • Crossed-strings method.
  • Unit sphere method.
  • Inside sphere method.
• Approximate quadrature methods commonly used in computer graphics
  • The hemicube algorithm.
  • The cubic tetrahedral algorithm.

Surface Integration Utilizing View Factor Algebra

The view factor kernel: \begin{equation} K_{e \leftrightarrow f} \equiv - \sphericalangle(\vec{r}_e, \vec{r}_f) \frac{\left( \hat{n}_e \cdot \vec{s}_{e \to f} \right) \left( \hat{n}_f \cdot \vec{s}_{e \to f} \right)}{\pi \| \vec{s}_{e \to f} \|^4} \end{equation} The view factor residual between finite surfaces $\partial \Omega_e$ and $\partial \Omega_f$: \begin{equation} R_{e \leftrightarrow f} \equiv \int_{\partial \Omega_e} \int_{\partial \Omega_f} K_{e \leftrightarrow f} dA_f dA_e \end{equation} The view factor from $\partial \Omega_e$ to $\partial \Omega_f$ can be computed by: \begin{equation} F_{e \to f} = \frac{R_{e \leftrightarrow f}}{A_e} \end{equation}

Generic Quadrature Rule

\begin{equation} \int_R f(x) dx \approx \sum_{i=1}^n f(x_i) \cdot w_i \end{equation}

• $x_i$ are quadrature points that lie within a reference domain $R$
• $w_i$ are the quadrature weights
• For a quadrature rule to be valid in view factor calculations, every $x_i$ must strictly be within the reference spatial domain and every $w_i > 0$

Change of Integration Domain

\begin{equation} \int_A u(x) dx = \int_R u(m(y)) \, | m'(y) | \, dy \end{equation}

• $R$ is the "reference" domain
• $A$ is the "actual" domain
• $m$ is a mapping between $R$ and $A$
• $R \subset \mathbb{R}^2$ for this application
• $A$ is embedded in $\mathbb{R}^3$
• $m : \mathbb{R}^2 \to \mathbb{R}^3$

View Function Computation

Challenges

• At a low level, resolving the view function involves computational geometry algorithms
• Computational geometry algorithms can occasionally be incorrect when implemented with floating poing numbers

• The Computational Geometry Algorithms Library
• Language: C++
• License: GPL/LGPL
• ~80k commits with ~80 contributors
• Implements all building blocks for the "ray tracing"
• Trivial to switch between exact and standard floating-point computation paradigms

Three Implementations of $\sphericalangle$

• NONE
• BRUTE_FORCE_WITHOUT_BBOX
• BRUTE_FORCE_WITH_BBOX

View Factor Benchmarks

Mesh Invariance

Geometric Parameters

• $a = b = 10$ m
• $c = 9$ m
• $F_{A_1 \to A_2} = F_{A_2 \to A_1} = 0.2285656684427082$

Discretizations

• Hexahedral meshes (HEX)
  • 100 surface elements
  • 1 m² per element
• Triangular prism meshes (PRI)
  • 244 surface elements
  • average element surface area of 0.40984 m²
• Triangular pyramid meshes (PYR)
  • 246 surface elements
  • average element surface area of 0.40650 m²

HEX-HEX Configuration

• $F_{A_1 \to A_2} = F_{A_2 \to A_1} = 0.2285656685503451$
• absolute error: 1.07637 × 10^-10
• relative error: 4.70923 × 10^-8%

PRI-HEX Configuration

PYR-HEX Configuration

PRI-PRI Configuration

• $F_{A_1 \to A_2} = F_{A_2 \to A_1} = 0.2285656684545930$
• absolute error: 1.18846 × 10^-11
• relative error: 5.19968 × 10^-9%

PYR-PRI Configuration

PYR-PYR Configuration

Mesh Refinement

HEX - HEX

HEX - HEX (cont.)

PYR - PYR

PYR - PYR (cont.)

Quadrature Refinement

HEX - HEX

PYR - PYR

Number Type Invariance

Geometric Parameters

• $r_1 = r_3 =2$ cm
• $r_2 =6$ cm
• $r_4 = 4$ cm
• $a = 15$ cm
• $F_{A_1 \to A_2} = 0.1131618692082199$

Discretization

• Triangular pyramid meshes
• Mesh size of 0.25 cm
• Smaller ring
  • 1,478 surface elements (on the surface of interest)
  • Average element surface area of 0.02550 cm²
• Larger ring
  • 3,860 surface elements (on the surface of interest)
  • Average element surface area of 0.02604 cm²

Unoccluded Rings Results

View Function Benchmark

Geometric Parameters

• $r_c =2$ cm
• $r_1 =6$ cm
• $r_2 = 4$ cm
• $h = 15$ cm
• $F_{A_1 \to A_2} = 0.0725555875115452$

Concave Discretization

• Triangular pyramid meshes
• Mesh size of 1.0 cm
• Smaller ring
  • 91 triangular surface elements (on the surface of interest)
  • Average element surface area of 0.414276 cm²
• Larger ring
  • 261 triangular surface elements (on the surface of interest)
  • Average element surface area of 0.385176 cm²
• Cylinder
  • 472 triangular surface elements (on the surface of interest)
  • Average element surface area of 0.399355 cm²

Concave Results

Convex Discretization

• Hexahedra meshes
• Smaller ring
  • 128 quadrangular surface elements (on the surface of interest)
  • Average element surface area of 0.294524 cm²
• Larger ring
  • 192 quadrangular surface elements (on the surface of interest)
  • Average element surface area of 0.523599 cm²
• Cylinder
  • 480 quadrangular surface elements (on the surface of interest)
  • Average element surface area of 0.392699 cm²

Convex Results

Real World Benchmark

Introduction

• Original work by He et al. (2018) as a means for exploring the design space for a component used in continuous annealing furnaces
• Results presented in the report were computed using the monte carlo method

Discretization

• Hexahedra mesh for all entities
• Plate
  • 40 × 40 × 4 elems
  • Physical dimensions change with each example
• W-Tube
  • Consistent physical dimension and mesh for each example

W-Tube Results (1/3)

W-Tube Results (2/3)

W-Tube Results (3/3)

Heat Transfer Benchmarks

Notes

The material properties in every simulation:

• $k = 1 \:\: W/m \cdot K$
• $\rho = 1 \:\: kg/m^3$
• $c_p = 1 \:\: J/kg \cdot K$

All plates are $10 \: m \times 10 \: m \times 1 \: m$

All plates are discretized by $40 \times 40 \times 4$ elementes

Radiation Emission from a Single Plate

Problem Description

• $q''' = 10,000 \:\: W/m^3$
• All faces except for one has adiabatic boundary conditions applied
• One $10 \: m \times 10 \: m$ face is emitting radiation in accordance with the Stefan–Boltzmann law

Temperature at Emitting Boundary

$Q_\text{gen} = A t q'''$

$Q_\text{rad} = A \sigma T_\text{rad, surf}^4$

$\implies T_\text{rad, surf} = \left( \frac{t q'''}{\sigma} \right)^{1/4}$

$T_\text{rad, surf} = 648.0329 \: K$

Analytic Steady State Solution

$-\frac{\partial^2 T}{\partial z^2} = q'''$

left BC: $\left. \frac{\partial T}{\partial z} \right|_{z=0} = 0$

right BC: $T(z=1) = T_\text{rad, surf}$

$\implies T(x, y, z) = T_\text{rad, surf} + \frac{q'''}{2} \left( 1 - z^2 \right)$

$T(x, y, 0) = 5,648.0329 \: K$

Simulations

• 2 simulations:
  • Local Stefan-Boltzmann law
  • Surface-averaged Stefan-Boltzmann law

Emitting Plate Temperature Profile

Emitting Plate Line Sample

Parallel and Directly Opposed Plates with Symmetric Heat Generation

Problem Description

• $q''' = 10,000 \:\: W/m^3$ in one plate
• $q''' = 0 \:\: W/m^3$ in the other plate
• All faces have adiabatic boundary conditions except for the faces facing eachother

Complete View

Plate with Heat Generation Adiabatic Face

Plate with Heat Generation Emitting Face

Linesamples of Plate with Heat Generation

Plate with Heat Generation Adiabatic Face

Plate with Heat Generation Emitting Face

Linsamples of Plate without Heat Generation

Parallel and Directly Opposed Plates with Asymmetric Dirichlet Boundary Condition

Problem Description

• $q''' = 0 \:\: W/m^3$ in both plates
• Dirichlet boundary conditions are applied at the top ($500 \: K$) and bottom ($273.15 \: K$) of the left plate
• All other faces have adiabatic boundary conditions except for the faces facing eachother

Complete View

Emitting Plate's Radiation Face

Receiving Plate's Radiation Face

Future Work

• Include additional methods to compute view factors
• Improve the efficiency of the view factor computation by parallelizing the quadrature computation
• Clustering of view factors could also be incorporated
• Extend boundary conditions to include gray-diffuse surfaces
• Additional studies could be conducted to verify, validate and compare Phoenix