Merge pull request #59 from RemDelaporteMathurin/yanns_remarks_ch2

Yanns remarks Chapter 2

chapters/abstract.tex

@@ -5,11 +5,11 @@ \chapter*{Abstract}
 As a radioactive isotope of hydrogen, tritium can represent a nuclear safety hazard and its inventory in the reactors materials must be controlled.
 In ITER, the tritium in-vessel safety limit is \SI{700}{g}.
 
-The tritium inventory of the ITER divertor was estimated using numerical modelling.
-The FESTIM code was developed to simulate hydrogen transport in tungsten monoblocks.
+The tritium inventory of the ITER divertor was numerically estimated.
+To this end, the FESTIM code was developed to simulate hydrogen transport in tungsten monoblocks.
 A parametric study was performed varying the exposure conditions (surface temperature and surface hydrogen concentration) and a behaviour law was extracted.
 This behaviour law provided a rapid way of estimating a monoblock inventory for a given exposure time and for given surface concentration and temperature.
-This behaviour law was then used and interfaced with output data from the edge-plasma code SOLPS-ITER.
+This behaviour law was then used and interfaced with output data from the edge-plasma code SOLPS-ITER in order to estimate the hydrogen inventory of the whole ITER divertor.
 Under conservative assumptions, the total hydrogen inventory (deuterium and tritium) was found to be well below the ITER tritium safety limit, reaching $\approx \SI{14}{g}$ after 25 000 pulses of \SI{400}{s}.
 
 To investigate the influence of helium exposure on these results, a helium bubble growth model was developed.

chapters/chapter2/model_description.tex

@@ -3,7 +3,7 @@
 \chapter{Model description}\labch{Chapter2}
 \label{Chapter2} % For referencing the chapter elsewhere, use \ref{Chapter2} 
 \section{Introduction}
-This Chapter describes the mathematical models needed to simulate H transport in tokamaks plasma facing components.
+This Chapter describes the mathematical models used to simulate H transport in tokamaks plasma facing components.
 The numerical scheme to solve these equations and an introduction to the finite element method alongside with its implementation in the \gls{festim} code are also presented.
 Finally, the verification \& validation of \gls{festim} is performed to guaranty its reliability.
 
@@ -35,7 +35,7 @@ \subsection{Macroscopic Rate Equations model}
 % This effect will enhance the diffusion when particles diffuse from a hot region to a cold region.
 % It should be noted however that the Soret coefficient is fairly difficult to measure.
 
-The spatio-temporal evolution of these concentrations is commonly described by the following reaction-diffusion system:
+The spatio-temporal evolution of these concentrations is commonly described by the following reaction-diffusion system \sidecite{mcnabb_new_1963}:
 
 \begin{equation}
     \frac{\partial c_\mathrm{m}}{\partial t}=\nabla \cdot J+\Gamma-\sum \frac{\partial c_{\mathrm{t}, i}}{\partial t}
@@ -55,10 +55,10 @@ \subsection{Macroscopic Rate Equations model}
 At steady state (i.e.\ $\frac{\partial c_\mathrm{m}}{\partial t} = 0$ and $\frac{\partial c_{\mathrm{t}, i}}{\partial t} = 0$), the mobile H concentration is independent of $c_{\mathrm{t}, i}$.
 \refeq{trapped} can be rewritten as:
 \begin{equation}
-    c_{\mathrm{t}, i} = n_i \frac{1}{\frac{p}{k c_\mathrm{m}} + 1}
+    c_{\mathrm{t}, i} = n_i \frac{1}{\frac{p_i}{k_i c_\mathrm{m}} + 1}
     \labeq{steady state ct}
 \end{equation}
-The quantity $(p / (k c_\mathrm{m}) + 1)^{-1}$ determines the trap occupancy. 
+The quantity $(p_i / (k_i c_\mathrm{m}) + 1)^{-1}$ determines the trap occupancy. 
 As it approaches one (high mobile concentration, low detrapping rate or high trapping rate), the trapped concentration approaches the trap density.
 As it approaches zero (high detrapping rate, low mobile concentration or low trapping rate), the trapped concentration approaches zero.
 
@@ -68,20 +68,20 @@ \subsection{Boundary conditions}
 
 \subsubsection{Dissociation and recombination fluxes}
 
-The concentration gradient can also be constrained on the boundaries (see \refeq{neuman robin bc}).
+The concentration gradient can be constrained on the boundaries (see \refeq{neuman robin bc}).
 
 \begin{equation}
     - D(T)\nabla c_\mathrm{m} \cdot \mathbf{n} = f(x, t) \quad \text { on } \partial \Omega
     \labeq{neuman robin bc}
 \end{equation}
-where $D(T) = D_0 \exp(\frac{-E_D}{k_B T}) $ is the diffusion coefficient in \si{m^2.s^{-1}}, $T$ is the local temperature in \si{K}, $\mathbf{n}$ is the boundary normal vector and $\partial \Omega$ is the domain boundary.
+where $D(T) = D_0 \exp(\frac{-E_D}{k_B T}) $ is the diffusion coefficient in \si{m^2.s^{-1}}, $\mathbf{n}$ is the boundary normal vector and $\partial \Omega$ is the domain boundary.
 
-$f$ can also be expressed as a molecular recombination flux:
+$f$ can be expressed as a molecular recombination flux:
 \begin{equation}
     - D(T)\nabla c_\mathrm{m} \cdot \mathbf{n} = K_r(T) c_\mathrm{m}^n \quad \text { on } \partial \Omega
     \labeq{recombination flux}
 \end{equation}
-where $K_r(T) = K_{r_0} \exp(\frac{-E_{K_r}}{k_B T}) $ is the recombination coefficient expressed in \si{m^{3n-2}.s^{-1}}, $\mathbf{n}$ is the boundary normal vector and $n \in \{1, 2\}$ is the order of the recombination.
+where $K_r(T) = K_{r_0} \exp(\frac{-E_{K_r}}{k_B T}) $ is the recombination coefficient expressed in \si{m^{3n-2}.s^{-1}}, $n \in \{1, 2\}$ is the order of the recombination.
 In a metal, $n=2$ and in a non-metallic liquid, $n=1$.
 Recombination occurs when hydrogen particles located at the surface of the material combine with other particles (which can be other hydrogen particles) and are no longer bonded with the metal surface.
 It can happen both in presence of a vacuum or when the metal is in contact with a fluid (gas or liquid).
@@ -91,7 +91,7 @@ \subsubsection{Dissociation and recombination fluxes}
     - D(T)\nabla c_\mathrm{m} \cdot \mathbf{n} = K_d(T) P \quad \text { on } \partial \Omega
     \labeq{dissociation flux}
 \end{equation}
-where $K_d(T) = K_{d_0} \exp(\frac{-E_{K_d}}{k_B T}) $ is the dissociation coefficient expressed in \si{m^{-3}.Pa^{-1}} and $\mathbf{n}$ is the boundary normal vector.
+where $K_d(T) = K_{d_0} \exp(\frac{-E_{K_d}}{k_B T}) $ is the dissociation coefficient expressed in \si{m^{-3}.Pa^{-1}}.
 Dissociation is the opposite process of recombination and occurs when particles in the surrounding atmosphere or fluid reach the metal surface and are adsorbed.
 These particles can then reach the bulk and diffuse in the metal.
 
@@ -108,8 +108,11 @@ \subsubsection{Dissociation and recombination fluxes}
 \subsubsection{Analytical simplification for an implanted source of H} \labsec{triangle model}
 
 Plasma implantation of hydrogen particles can be modelled with a volumetric source.
-Typically, the depth of the implantation profile is a few nanometres depending on the incident particles energy and incident angle.
+Typically, the depth of the implantation profile $R_p$ is a few nanometres depending on the incident particles energy and incident angle.
 These profiles can be simulated by Monte Carlo codes like \gls{srim} \sidecite{ziegler_srim_2010} and have a gaussian-like shape  (see \reffig{srim_implantation_profile_example}).
+% In order to accurately model this source term, the size of the cells constituing the mesh (the spatial discretisation of the domain) must be less than a nanometre.
+% This can be done easily in 1D but is very complicated in higher dimensions, especially when simulating centimetre-sized components.
+This volumetric source term can be simplified into a Dirichlet boundary condition (i.e.\ enforcing the mobile particle concentration at the exposed surface).
 
 \begin{figure}
     \centering
@@ -118,10 +121,6 @@ \subsubsection{Analytical simplification for an implanted source of H} \labsec{t
     \labfig{srim_implantation_profile_example}
 \end{figure}
 
-% In order to accurately model this source term, the size of the cells constituing the mesh (the spatial discretisation of the domain) must be less than a nanometre.
-% This can be done easily in 1D but is very complicated in higher dimensions, especially when simulating centimetre-sized components.
-This volumetric source term can be simplified into a Dirichlet boundary condition (i.e.\ enforcing the mobile particle concentration at the exposed surface).
-
 Let us consider a volumetric source term of hydrogen $\Gamma = \varphi_\mathrm{imp} \; f(x)$ where $f$ is a narrow Gaussian distribution.
 The concentration profile of mobile particles can be approximated by a triangular shape \sidecite{schmid_diffusion-trapping_2016} (see \reffig{recomb sketch}).
 
@@ -133,7 +132,7 @@ \subsubsection{Analytical simplification for an implanted source of H} \labsec{t
 \end{figure*}
 
 The concentration profile is therefore maximum at $x=R_p$.
-The expression of $c_\mathrm{max}$ can be obtained by expressing the flux balance at equilibrium:
+The expression of maximum concentration value $c_\mathrm{max}$ can be obtained by expressing the flux balance at equilibrium:
 
 \begin{equation}
     \varphi_\mathrm{imp} = \varphi_\mathrm{recomb} + \varphi_\mathrm{bulk}
@@ -191,7 +190,7 @@ \subsubsection{Analytical simplification for an implanted source of H} \labsec{t
 \begin{equation}
     \tau = \frac{R_p \sum R_i \, n_i}{8 \varphi_\mathrm{imp}}
 \end{equation}
-In this expression, $R_i = (p / (k c_\mathrm{max}) + 1)^{-1}$ represents the maximum filling ratio of the trap $i$ and $n_i$ is the trap density.
+In this expression, $R_i = (p_i / (k_i \, c_\mathrm{max}) + 1)^{-1}$ represents the maximum filling ratio of the trap $i$ and $n_i$ is the trap density.
 When $t \gg \tau$, $c_\mathrm{max}(t) \approx \frac{R_p \varphi_\mathrm{imp}}{D} + \sqrt{\frac{\varphi_\mathrm{imp}}{K_r}}$
 
 \subsection{Interface condition: conservation of chemical potential}\labsec{conservation of chemical potential}
@@ -269,41 +268,41 @@ \section{Heat transfer}
     -\lambda \nabla T \cdot \mathbf{n} = f(x, y, z, t) \quad \text { on } \partial \Omega
     \labeq{neumann bc T}
 \end{equation}
-where $\lambda$ is the thermal conductivity in \si{W.m^{-1}.K^{-1}}, $\mathbf{n}$ is the boundary normal vector and $\partial \Omega$ is the domain boundary.
+where $\lambda$ is the thermal conductivity in \si{W.m^{-1}.K^{-1}}.
 
 Finally, to model a convective heat flux when the surface is in contact with a fluid (e.g.\ cooling pipes, natural convection...), a Robin boundary condition needs to be employed (see \refeq{convective bc T}).
 \begin{equation}
     -\lambda \nabla T \cdot \mathbf{n} = h (T - T_\mathrm{ext}) \quad \text { on } \partial \Omega
     \labeq{convective bc T}
 \end{equation}
-where $h$ is the heat transfer coefficient in \si{W.m^{-2}.K^{-1}}, $T_\mathrm{ext}$ is the fluid temperature in \si{K} and $\partial \Omega$ is the domain boundary.
+where $h$ is the heat transfer coefficient in \si{W.m^{-2}.K^{-1}}, $T_\mathrm{ext}$ is the fluid temperature in \si{K}.
 The heat transfer coefficient can be dependent on the temperature and the flow characteristics.
 It is obtained by computing the Nusselt number from correlations linking it to the Reynolds number of the flow and the Prandtl number of the fluid \sidecite{poirier_correlations_2016} (e.g.\ Dittus-Boetler, Sieder-State, Gnielinski...).
 Once the Nusselt number is known, the heat transfer coefficient $h$ reads:
 \begin{equation}
     h = \frac{\lambda \, \textit{Nu}}{L}
 \end{equation}
-with $\lambda$ the thermal conductivity of the fluid in \si{W.m^{-1}.K^{-1}}, $\textit{Nu}$ the Nusselt number and $L$ the characteristic length in \si{m}.
+with $\textit{Nu}$ the Nusselt number and $L$ the characteristic length in \si{m}.
 
 
 \section{Implementation}
 
 
-The models described in this Section can be hard to solve analytically for complex problems (complex geometries, transients, combined boundary conditions, etc.).
+The models described in the previous sections can be hard to solve analytically for complex problems (complex geometries, transients, combined boundary conditions, coupling, etc.).
 The code \gls{festim} \sidecite{delaporte-mathurin_finite_2019} was therefore developed in order to numerically solve these coupled equations.
 
 \subsection{The finite element method: FEniCS}
 \gls{festim} is based on the Finite Element Method to solve this set of differential equations and boundary conditions.
 Several finite element libraries are available open-source (deal.II \sidecite{arndt_dealii_2021}, MFEM \sidecite{kolev_tzanio_modular_2010}, MOOSE \sidecite{permann_moose_2020}, FreeFEM++ \sidecite{hecht_new_2012}...).
-The open-source python/C++ package \gls{fenics} \sidecite{alnaes_fenics_2015} was employed.
+The open-source python/C++ package \gls{fenics} \sidecite{alnaes_fenics_2015} was employed for it provides a user-friendly python interface.
 The finite element method is a versatile tool that can solve any partial differential equation on an arbitrary geometry in 1D, 2D or 3D.
 The main advantage of this method compared to the finite difference method is the simplicity of its application to complex geometries and unstructured meshes.
 Indeed, implementing a finite difference scheme for such a problem would be tedious and extra care must be taken for mistakes in the implementation could result in losses in efficiency and accuracy of the numerical solution.
 
 This section illustrates the finite element method applied to Poisson's equation \sidecite{logg_automated_2012}.
 
 The mathematical problem can be described by \refeq{example poisson} where $u$ is the unknown to be solved governed by Poisson's equation.
-The problem is constrained by boundary conditions defined on $\partial \Omega$, the boundary of the domain.
+The problem is constrained by boundary conditions defined on $\partial \Omega$, the boundary of the domain $\Omega$.
 The value of $u$ is prescribed on the subset $\Gamma_D$ (Dirichlet boundary condition) whereas the value of the normal derivative of $u$ is prescribed on the remaining boundary $\Gamma_N$ (Neumann boundary condition) (see \reffig{fe problem sketch}).
 
 \begin{subequations}

chapters/chapter2/verification_and_validation/Parametric_optimisation/methodology.tex

@@ -15,17 +15,17 @@
 In \refeq{cost function}, $f(\textbf{x})$ can be weighted by coefficients $\alpha_i$ in order to have a better fit on specific regions of the spectrum.
 %Note that this cost function could as well be a root mean square error or any type of residual.
 
-The parametric optimisation problem can now be solved by finding the minimum of the cost function $f$.
-The global optimisation routine is illustrated on \reffig{diagramm}.
-\begin{figure}
+% The global optimisation routine is illustrated on \reffig{diagram}.
+\begin{figure*}
     \centering
     \includegraphics[width=\linewidth]{Figures/Chapter3/Parametric_optimisation/algorithm diagram.pdf}
-    \caption{Diagram of the embedding of FESTIM within a parametric optimisation routine based on SciPy \cite{virtanen_scipy_2020}.}
-    \labfig{diagramm}
-\end{figure}
+    \caption{Diagram of the parametric optimisation routine.}
+    \labfig{diagram}
+\end{figure*}
 
-A comparative study of the several optimisation algorithms which can be employed has been made.
-These algorithms require the user to give an initial set of parameters called \textit{initial guess} and evaluate the cost function with several parameters sets until the convergence criterion is reached.
+The parametric optimisation problem can now be solved by finding the minimum of the cost function $f$.
+A comparative study of the several optimisation algorithms which can be employed has been performed.
+These algorithms require the user to give an initial set of parameters called \textit{initial guess} and evaluate the cost function with several parameters sets until the convergence criterion is reached (see \reffig{diagram}).
 As in \sidecite{drexler_model-based_2019}, the Python package SciPy \sidecite{virtanen_scipy_2020} will be employed.
 
 Four minimisation algorithms have been benchmarked against a test case.

chapters/chapter2/verification_and_validation/analytical_verification.tex

@@ -39,21 +39,22 @@ \subsubsection{Case 1: H transport (\gls{mes})} \labsec{analytical}
 \end{equation}
 
 In our case, we choose the trapping and detrapping rates $k$ and $p$, the concentration $c_0$ and the temperature $T$ so that $\zeta \gg \frac{c_\mathrm{m}}{n}$.
+This condition is equivalent to having the trap filling ratio $(p / (k \, c_\mathrm{m}) + 1)^{-1} \ll 1$.
 This is known as the \textit{effective diffusivity regime} where the diffusion is almost identical to the case where there are no traps.
 In this regime, the governing equations are identical as a pure diffusion regime and are therefore easy to solve analytically.
 
 The coefficient $D$ is then replaced by an effective diffusion coefficient:
 \begin{equation}
     D_\mathrm{eff} = \frac{D}{1+\frac{1}{\zeta}}
 \end{equation}
-The particle flux at the background surface ($x=l$) is expressed in $\si{H.m^{-2}.s^{-1}}$ and finally defined in \sidecite{longhurst_verification_2005} by:
+The particle flux at the background surface ($x=l$) is expressed in $\si{H.m^{-2}.s^{-1}}$ and finally defined in \cite{longhurst_verification_2005} by:
 \begin{equation}
     \varphi(t) = \frac{c_0 D}{l}\bigg[1+2\sum_{m=1}^{\infty}(-1)^m \exp\bigg(-m^2\frac{\pi^2 \:D_\mathrm{eff} \: t}{l^2}\bigg)\bigg]
 \labeq{flux analytical}
 \end{equation}
-Note: The infinite sum has been truncated at $m=10000$.
+In \refeq{flux analytical}, the infinite sum has been truncated at $m=10000$.
 
-All the parameters are defined in \reftab{parameters analytical verification}.
+All the parameters used in this verification case are defined in \reftab{parameters analytical verification}.
 These parameters have been chosen for the sake of verification and do not necessarily represent realistic conditions as verification is a mathematical exercise.
 \begin{table}
     \centering
@@ -159,7 +160,7 @@ \subsubsection{Case 2: H transport (\gls{mms})} \labsec{mms}
 \begin{equation}
     \begin{cases}
     f = \cos(t) - \sin(t) - 2D \\
-    g_1 = \nu_1 c_{{t,1}_D} - \nu_m c_{m_D} ( n_1 - c_{{t,1}_D}) - \sin(t)
+    g_1 = p_1 c_{{t,1}_D} - k_1 c_{m_D} ( n_1 - c_{{t,1}_D}) - \sin(t)
     \end{cases}
     \labeq{sources}
 \end{equation}