Implementation of Monte Carlo method
This page detail the implementation of Monte Carlo method, see [Modest] and [Ozturk] for more details.
The code are in project : github.com/feelpp/hygrothermie
In folder src/radiative-heat-transfert (github.com/feelpp/hygrothermie/tree/model/src/radiative-heat-transfert)
1. Algorithm in pseudo-code
For this stage, I choose three algortihms (principle difference the order of execution of different element)
1.1. Element
This section present different function use to global algorithm.
1.1.1. Generate Source
\begin{algorithm}
\caption{Source Point}
\begin{algorithmic}
\FUNCTION{IsInTriangle}{$\mathcal{T}$, $Q$}
\STATE $\mathcal{T}$ triangle of origin point $P$ and directional vectors $\mathbf{e_1}$ and $\mathbf{e_2}$
\STATE $u \leftarrow \frac{(\mathbf{e_2} \cdot \mathbf{e_2}) (\vec{PQ} \cdot \mathbf{e_1}) - (\mathbf{e_1} \cdot \mathbf{e_2}) (\vec{PQ} \cdot \mathbf{e_2}) }{(\mathbf{e_1} \cdot \mathbf{e_1})(\mathbf{e_2} \cdot \mathbf{e_2})-(\mathbf{e_1} \cdot \mathbf{e_2})^2}$
\STATE $v \leftarrow \frac{(\mathbf{e_1} \cdot \mathbf{e_1}) (\vec{PQ} \cdot \mathbf{e_2}) - (\mathbf{e_1} \cdot \mathbf{e_2}) (\vec{PQ} \cdot \mathbf{e_1}) }{(\mathbf{e_1} \cdot \mathbf{e_1})(\mathbf{e_2} \cdot \mathbf{e_2})-(\mathbf{e_1} \cdot \mathbf{e_2})^2}$
\RETURN $(0 \leq u \leq 1)(0 \leq v \leq 1)$
\ENDFUNCTION
\FUNCTION{SourcePoint}{$\mathcal{T}$}
\STATE $\mathcal{T}$ triangle of origin point $P$ and directional vectors $\mathbf{e_1}$ and $\mathbf{e_2}$
\STATE $\xi_1$ and $\xi_2$ random number in $[0,1]$
\STATE $\phi \leftarrow 2\pi \xi_1$
\STATE $r \leftarrow R \sqrt{\xi_2}$
\STATE $Q \leftarrow \begin{pmatrix} x_C + r cos \phi \\ y_C + r sin \phi \end{pmatrix}$
\WHILE{ IsInTriangle($\mathcal{T}$, $Q$)}
\STATE $\xi_1$ and $\xi_2$ random number in $[0,1]$
\STATE $\phi \leftarrow 2\pi \xi_1$
\STATE $r \leftarrow R \sqrt{\xi_2}$
\STATE $Q \leftarrow \begin{pmatrix} x_C + r cos \phi \\ y_C + r sin \phi \end{pmatrix}$
\ENDWHILE
\RETURN $Q$
\ENDFUNCTION
\end{algorithmic}
\end{algorithm}
1.1.2. Generate Ray
\begin{algorithm}
\caption{Ray}
\begin{algorithmic}
\FUNCTION{Ray}{$\mathcal{T}$, $s$}
\STATE $\mathcal{T}$ triangle of origin point $P$ and directional vectors $\mathbf{e_1}$ and $\mathbf{e_2}$
\STATE $s$ source point
\STATE $\phi_n \leftarrow $ incident angle of $normal(\mathcal{T})$
\STATE $\theta_n \leftarrow $ azimuthal angle of $normal(\mathcal{T})$
\STATE $\phi_{random} \leftarrow$ random number in $[0,\pi]$
\STATE $\theta_{random} \leftarrow$ random number in $[0,2\pi]$
\STATE $\phi \leftarrow \phi_n + \phi_{random}$
\STATE $\theta \leftarrow \theta_n + \theta_{random}$
\RETURN $\phi$, $\theta$
\ENDFUNCTION
\end{algorithmic}
\end{algorithm}
1.1.3. Visibility
\begin{algorithm}
\caption{Visibility}
\begin{algorithmic}
\FUNCTION{Visibility}{$\mathcal{T}_1$, $\mathcal{T}_2$}
\STATE $\mathbf{n_1} \leftarrow$ $normal(\mathcal{T}_1)$
\STATE $\mathbf{n_2} \leftarrow$ $normal(\mathcal{T}_2)$
\RETURN $\left( \mathbf{n_1} \cdot \mathbf{n_2} < 0 \right)$
\ENDFUNCTION
\end{algoritmic}
\end{algorithm}
1.1.4. Interception
\begin{algorithm}
\caption{Interception}
\begin{algorithmic}
\FUNCTION{Interception}{$r$,$\mathcal{T}$}
\STATE $\mathcal{T}$ triangle of origin point $P$ and directional vectors $\mathbf{e_1}$ and $\mathbf{e_2}$
\STATE $r$ is ray, compose by origin point $O$ and directional vector $\mathbf{v}$
\STATE $div = \left( \mathbf{v} \wedge \mathbf{e_2} \right) \cdot \mathbf{e_1} \ne 0$
\IF{$div \ne 0$}
\STATE $u = \frac{ \vec{OP} \cdot \left( \mathbf{e_2} \wedge \mathbf{v} \right)}{div}$
\IF{$0 \leq u \leq 1$}
\STATE $t = \frac{ \vec{OP} \cdot \left( \mathbf{e_1} \wedge \mathbf{v} \right)}{div}$
\IF{$0 \leq t \leq 1$}
\RETURN 1
\ENDIF
\ENDIF
\ENDIF
\RETURN 0
\ENDFUNCTION
\end{algorithmic}
\end{algorithm}
1.2. First Strategy
Here, for all element \(e^i_k\) of surface \(A_i\), we generate once a set of ray and for all element \(e^i_l\), we browse the set of ray to know wich ray is intercept.
\begin{algorithm}
\caption{ViewFactor}
\begin{algorithmic}
\FUNCTION{ViewFactor1}{$A_i, A_j, N_s, N_r$}
\STATE $A_i$ is surface decomposed by finite element $e^i_k$
\STATE $A_j$ is surface decomposed by finite element $e^j_l$
\STATE $N_s$ number of source point of one element
\STATE $N_r$ number of ray from one source point
\STATE $Area_i \leftarrow 0$ : area of surface $A_i$
\STATE $\Sigma_2 \leftarrow 0$
\FOR{$e^i_k$ in set of element of surface $A_i$}
\FOR{$1 \leq i_s \leq N_s$}
\STATE generate a source point $s$
\FOR{$1 \leq i_r \leq N_r$}
\STATE generate a ray $r$ from source point $s$
\ENDFOR
\ENDFOR
\STATE $\Sigma_1 \leftarrow 0$
\FOR{$e^j_l$ in set of element of surface $A_j$}
\IF{$e^j_k$ visible from $e^i_l$}
\STATE $n \leftarrow 0$ : number of ray intercept
\FOR{ray $r$ in set of ray}
\IF{$r$ is intercept by $e^j_k$}
\STATE $n \leftarrow n + 1$
\ENDIF
\ENDFOR
\STATE $F_{local} \leftarrow \frac{n}{N_s N_r}$
\STATE $\Sigma_1 \leftarrow \Sigma_1 + F_{local}$
\ENDIF
\ENDFOR
\STATE $Area_i \leftarrow Area_i + Area(e^i_k)$
\STATE $\Sigma_2 \leftarrow \Sigma_2 + Area(e^i_k) \Sigma_1$
\ENDFOR
\STATE $F_{global} \leftarrow \frac{1}{Area(A_i)} \Sigma_2$
\RETURN $F_{global}$
\ENDFUNCTION
\end{algorithmic}
\end{algorithm}
Complexity : \(O\left( K (L+1) N_s N_r \right) = O\left( K L N_s N_r \right)\), with \(K\) number of element of \(A_i\), \(L\) number of element of \(A_j\), \(N_s\) number of source point and \(N_r\) number of ray from one source point.
Complexity more precise : \(K \times \left(N_s \times complexity(\mathbf{SourcePoint}) \times N_r complexity(\mathbf{Ray}) + L \times complexity(\mathbf{Visibility}) \times N_r \times N_s \times complexity(\mathbf{Interception}) \right)\)
Disadvantage : demand more memories. (not really important)
1.3. Second Strategy
Here, for all element \(e^i_k\) of surface \(A_i\), we generate once a set of ray and for all ray, we browse the set of element \(e^j_l\) to know wich element intercept the ray.
\begin{algorithm}
\caption{ViewFactor}
\begin{algorithmic}
\FUNCTION{ViewFactor2}{$A_i, A_j, N_s, N_r$}
\STATE $A_i$ is surface decomposed by $K$ finite element $e^i_k$
\STATE $A_j$ is surface decomposed by $L$ finite element $e^j_l$
\STATE $N_s$ number of source point of one element
\STATE $N_r$ number of ray from one source point
\STATE $Area_i \leftarrow 0$ : area of surface $A_i$
\STATE $\Sigma_2 \leftarrow 0$
\FOR{$e^i_k$ in set of element of surface $A_i$}
\STATE $n \leftarrow \{n_l\}_{1\leq l \leq L}$ is a list of zeros
\FOR{$1 \leq i_s \leq N_s$}
\STATE generate a source point $s$
\FOR{$1 \leq i_r \leq N_r$}
\STATE generate a ray $r$ from source point $s$
\STATE $l \leftarrow 0$
\WHILE{$r$ isn't intercept by $e^j_l$ and $l \leq L$}
\STATE $l \leftarrow l+1$
\ENDWHILE
\IF{$l \leq L$}
\STATE $n_j \leftarrow n_j+1$
\ENDIF
\ENDFOR
\ENDFOR
\STATE $Area_i \leftarrow Area_i + Area(e^i_k)$
\STATE ($F_{e^i_k \rightarrow e^j_l} \leftarrow \frac{n_j}{N_s N_r}$)
\STATE $F_{e^i_k \rightarrow A_j} \leftarrow \frac{\Sigma(n_j)}{N_s N_r} \mbox{ or } = \Sigma(F_{e^i_k \rightarrow e^j_l} )$
\STATE $\Sigma_2 \leftarrow \Sigma_2 + Area(e^i_k) F_{local} $
\ENDFOR
\STATE $F_{global} \leftarrow F_{global} + \frac{1}{Area(A_i)} \Sigma_2$
\RETURN $F_{global}$
\ENDFUNCTION
\end{algorithmic}
\end{algorithm}
| The interception must contain the visibility between \(e^i_k\) and \(e^j_l\). |
Complexity : worst case \(O\left( KN_sN_r \frac{L(L+1)}{1} \right)\), with \(K\) number of element of \(A_i\), \(L\) number of element of \(A_j\), \(N_s\) number of source point and \(N_r\) number of ray from one source point.
1.4. Third Strategy
Here, we take all couple of element \(e^i_l\) and \(e^j_k\) and we calculate the correspondant view factor.
\begin{algorithm}
\caption{Calculation of View Factor}
\begin{algorithmic}
\FUNCTION{ViewFactor3}{$A_i$, $A_j$, $N_s$, $N_r$}
\STATE $A_i$ is surface decomposed by finite element $e^i_k$
\STATE $A_j$ is surface decomposed by finite element $e^j_l$
\STATE $N_s$ number of source point of one element
\STATE $N_r$ number of ray from one source point
\STATE $Area_i \leftarrow 0$ : area of surface $A_i$
\STATE $\Sigma_2 \leftarrow 0$
\FOR{$e^i_k$ in set of element of surface $A_i$}
\STATE $\Sigma_1 \leftarrow 0$
\FOR{$e^j_l$ in set of element of surface $A_j$}
\IF{$e^i_k$ and $e^j_l$ are visible}
\STATE $n \leftarrow 0$
\FOR{$1 \leq i_s \leq N_s$}
\STATE generate source point $s$
\FOR{$1 \leq i_r \leq N_r$}
\STATE generate ray $r$
\IF{$r$ is intercept by $e^j_l$}
\STATE $n \leftarrow n + 1$
\ENDIF
\ENDFOR
\ENDFOR
\STATE $F_{e^i_k \rightarrow e^j_l} \leftarrow \frac{n}{N_s N_r}$
\ENDIF
\STATE $\Sigma_1 \leftarrow \Sigma_1 + F_{e^i_k \rightarrow e^j_l}$
\ENDFOR
\STATE $Area_i \leftarrow Area_i + Area(e^i_k)$
\STATE $\Sigma_2 \leftarrow \Sigma_2 + Area(e^i_k) \Sigma_1$
\ENDFOR
\STATE $F_{global} \leftarrow \frac{1}{Area(A_i)} \Sigma_2$
\RETURN $F_{global}$
\ENDFUNCTION
\end{algorithmic}
\end{algorithm}
Complexity : \(O\left( K L N_s N_r \right)\), with \(K\) number of element of \(A_i\), \(L\) number of element of \(A_j\), \(N_s\) number of source point and \(N_r\) number of ray from one source point.
Complexity more precise : \(K \times L \times N_s \times complexity(\mathbf{SourcePoint}) \times N_r \times complexity(\mathbf{Ray}) \times complexity(\mathbf{Interception})\)
2. Structutre of code
(not definitive)
-
rht_option.cpp : the option of radiative-heat-transfert application
-
rht_main.cpp : the main script of radiative-heat-transfert application
-
geometrical.hpp : contains some function and class to resolve geometrical problem
-
radiative_heat_transfert.hpp : contains the class wich calculate view factor
3. Detail of geometrical.hpp
3.1. Class
This script contains some geometrical class :
-
Point : represent Point in three dimension
-
Vector : represent Vector in three dimension
-
VectorUnit (sub-class of Vector) : represent unit vector, this class inherited the cartesian coordonated of Vector and have polar coordonated
-
S_phi : incident angle
-
S_theta : azimuthal angle
-
-
Square : represent a square, build with 4 vertex and 1 direction (+1 or -1) to have the direction of emission or absorption, it have the arguments:
-
S_side : side of square
-
S_R : ray of square (the ray of circumscribed circle)
-
S_midlle : middle of square
-
S_normal : directional normal of surface wich contain the square, the sense of direction give the side of emission (or absorption) of surface
-
Some overload of Point and Vector is implement (rule of caculation : addition, substraction, division …) adn function of vectorial and scalar produce in three dimension.
3.2. Function
There are function :
-
generate direction(VectorUnit<T> normal) : generate a VectorUnit (direction of emissive ray) around a Vector normal (normal of square) with distribution (for more detail see the section Distribution of rays) :
-
\(\phi = \pi \xi_1\) with \(\xi_1\) follow a uniform distribution between \([0,1\)]
-
\(\theta = 2\pi \xi_2\) with \(\xi_2\) follow a uniform distribution between \([0,1\)]
-
-
IsInSquare(Square<T> square, Point<T> point) : say if point is in square (for more detail see In Square ?)
-
IsIntercept(Point<T> source, VectorUnit<T> direction, Square<T> square) : say if the ray of source source and direction direction is intercept by square (for more detail see Interception of ray)
-
Point<T> generate_source(Square<T> square) : generate a source point of square (fore more detail see Distribution of point source of ray)
4. Detail of radiative_heat_transfert.hpp
radiative_heat_transfert.hpp contains the class :
-
local_view_factor(Square<T> S1, Square<T> S2, int N_s, int N_r) : calculate the view factor between two square with :
-
N_s number of source point in emissive surface
-
N_r number of emissive ray of one source point
-
References
-
[Incropera] Fundamentals of Heat and Mass Transfer, Theodore L. Bergman, Adrienne S. Lavine, Frank P. Incropera, David P. DeWitt, Wiley.
-
[Kumaran1996] Kumaran, M. Kumar, Heat, Air and Moisture Transfer in Insulated Envelope Parts, Final Report, Volume 3.
-
[Ozturk] Abdurrahman Ozturk, Implementation of View Factor Model and Radiative Heat Transfer Model in MOOSE. University of South Carolina, Download PDF
-
[Modest] M.F. Modest, Radiative Heat Transfet. Elsevier Science, 2013, ISBN 9780123869906, Chapter 8
-
[Visionary] Visionaray: A Cross-Platform Ray Tracing Template Library, Zellmann, Stefan and Wickeroth, Daniel and Lang, Ulrich, Proceedings of the 10th Workshop on Software Engineering and Architectures for Realtime Interactive Systems (IEEE SEARIS 2017), 2017
-
[Visionary::github] Visionaray, Zellman, github.com/szellmann/visionaray