Merge pull request #68 from RemDelaporteMathurin/kschmid_remarks

Reader 1's remarks

bibfile.bib

@@ -10585,3 +10585,32 @@ @article{brown_endfb-viii0_2018
 	pages = {1--142},
 	file = {ScienceDirect Full Text PDF:D\:\\Logiciels\\data_zotero\\storage\\2I9PI5MF\\Brown et al. - 2018 - ENDFB-VIII.0 The 8th Major Release of the Nuclea.pdf:application/pdf;ScienceDirect Snapshot:D\:\\Logiciels\\data_zotero\\storage\\JFGW9MM5\\S0090375218300206.html:text/html},
 }
+
+@techreport{waelbroeck_influence_1984,
+	address = {Germany},
+	title = {Influence of bulk and surface phenomena on the hydrogen permeation through metals},
+	abstract = {We discuss the permeation of hydrogen through metals and alloys such as iron, nickel,
+steels and Inconel wherein H dissolves endothermically from an H2 gas We assume first
+that trapping centers, surface contamination layers, the saturation of the H surface
+coverage and the implantation profile - when energetic ions drive the permeation -
+can be neglected, that a quasi-equilibrium exists between the H atom concentration
+ν in the adsorbed layer and c in the near surface layers and that the H solubility
+and diffusivity are homogeneous in the membrane We evaluate thereafter separately
+the influence of these various effects and identify the parameter domains where appreciable
+corrections result The permeation phenomenon is complex even when these simplifications
+are made: the penetration rate is proportional to the flux of thermal molecules, atoms
+or energetic ions - depending upon the case - which strike the surface; the diffusion
+in the metal is proportional to the gradient of c; the release rate depends on c2;
+the time-dependent diffusion equation includes a double spatial derivative of c Permeation
+can only be fully described when computer codes such as PERI is used Simple analytical
+relations are however obtained in several limiting cases They are the object of this
+report Some of them had already been derived by other authors but they were not shown
+to be part of a single, self consistent permeation model A comparison of predicted
+and experimental results shows that the simplified model describes surprisingly accurately
+the hydrogen exchange between gas and metal solutions (orig/GSCH)},
+	author = {Waelbroeck, F. and Wienhold, P. and Winter, J. and Rota, E. and Bauno, T.},
+	year = {1984},
+	note = {Juel--1966
+INIS Reference Number: 16053332},
+	pages = {201},
+}

chapters/chapter2/model_description.tex

@@ -82,6 +82,8 @@ \subsubsection{Dissociation and recombination fluxes}
     \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}}, $n \in \{1, 2\}$ is the order of the recombination.
+This is the Waelbroeck model \sidecite{waelbroeck_influence_1984}.
+This model may not be valid for all cases \sidecite{schmid_use_2021} and a more extended and more complex model from Pick–Sonnenberg \sidecite{pick_model_1985} could be more generic.
 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).
@@ -230,12 +232,17 @@ \subsection{Interface condition: conservation of chemical potential}\labsec{cons
     \frac{\partial \phi S}{\partial t} &= \nabla \cdot\left(D \nabla \phi S\right) + f \nonumber \\
     &=\nabla \cdot\left( D S \nabla \phi + D \phi \nabla S\right) + f \labeq{diffusion equation changed}
 \end{align}
-
+where $f$ is some source term.
 % Because $\phi$ is computed, the ratio $c_\mathrm{m}/S$ is continuous by default at the material interfaces.
 
 % This second approach is used for instance in the \textit{mass-diffusion} procedure of the Abaqus code \sidecite{smith_abaqusstandard_2009}.
 % This interface model has also been implemented into the current hydrogen transport code FESTIM \sidecite{delaporte-mathurin_finite_2019} using FEniCS \sidecite{alnaes_fenics_2015}.
 
+\refeq{trapped} reads:
+\begin{equation}
+    \frac{\partial c_{\mathrm{t}, i}}{\partial t}=k_i \ \phi \ S \ \left(n_{i}-c_{\mathrm{t}, i}\right)-p_i \ c_{\mathrm{t}, i}
+\end{equation}
+
 After solving \refeq{diffusion equation changed} for $\phi$, $c_m$ can be retrieved by multiplying the solution by $S$.
 
 \section{Heat transfer}
@@ -402,6 +409,18 @@ \subsection{The finite element method: FEniCS}
 After solving \refeq{variational problem matrix form}, the $U_i$ coefficients are known and the approximated solution $u_h$ can be computed.
 Non-linear problems are solved similarly where the solution is approached using Newton's method.
 
+
+The weak formulation of the steady state McNabb \& Foster equation is:
+
+\begin{equation}
+    \int_{\Omega} \nabla c_\mathrm{m} \cdot \nabla v_\mathrm{m}  \ dx = \int_{\Omega} \left( k c_\mathrm{m} (n - c_\mathrm{t}) - p c_\mathrm{t} \right) \ v_\mathrm{t} \ dx \quad \forall \ (v_\mathrm{m}, v_\mathrm{t}) \in \hat{V}
+    \labeq{weak form mcnabb foster}
+\end{equation}
+where $\hat{V}$ is a suitable vector-functionspace.
+
+The transient formulation can be obtained by adding transient terms to \refeq{weak form mcnabb foster} $\int_{\Omega} \frac{c_\mathrm{m} - c_\mathrm{m, old}}{\Delta t} v_\mathrm{m} \ dx$ and $\int_{\Omega} \frac{c_\mathrm{t} - c_\mathrm{t, old}}{\Delta t} v_\mathrm{t} \ dx$ (forward Euler time discretisation).
+$c_\mathrm{m, old}$ and $c_\mathrm{t, old}$ are the previous solutions for the mobile concentration and trapped concentration respectively and $\Delta t$ is the timestep.
+
 \subsection{Main features of \gls{festim}}
 \gls{festim} provides an even higher level of abstraction than \gls{fenics} by providing a user-friendly interface dedicated to H transport and heat transfer problems.
 Users only have to provide inputs such as material properties, traps properties, geometry, solving parameters, without having to dive into the finite element implementation.

chapters/chapter3/monoblocks.tex

@@ -280,7 +280,7 @@ \subsection{Influence of cycling}\labsec{influence of cycling}
 \begin{figure}
     \centering
     \includegraphics[width=\linewidth]{Figures/Chapter3/monoblocks/comparison_inventory_cycling.pdf}
-    \caption{Evolution of the monoblock inventory as a function of the implanted fluence for cycled (solid) and continuous (dashed) exposure.}
+    \caption{Evolution of the monoblock inventory as a function of the implanted fluence for cycled (solid) and continuous (dashed) exposure on a 1D case.}
     \labfig{monoblock inventory cycling}
 \end{figure}
 

chapters/chapter3/monoblocks/parametric_study.tex

@@ -15,6 +15,16 @@
 \subsection{Assumptions and simplifications}
 
 For the sake of simplicity and computational time, Trap 2 (in W) is neglected.
+This trap was neglected for it has the lowest density.
+However this will induce errors (< \SI{80}{\%}) in the inventory computation (see \reffig{neglecting trap 2 error})
+
+\begin{figure}
+    \centering
+    \includegraphics[width=\linewidth]{Figures/Chapter3/monoblocks/neglecting_trap_2_error.pdf}
+    \caption{Evolution of the concentrations of traps 1 and 2 as a function of temperature for a local mobile hydrogen concentration of \SI{1e22}{m^{-3}} and error associated with neglecting trap 2.}
+    \labfig{neglecting trap 2 error}
+\end{figure}
+
 Note that this method could be applied to any set of trapping parameters.
 Continuity of mobile concentration at interfaces between materials is also assumed in order to save computational time (see \refsec{influence of interface conditions}).
 To remain conservative, no recombination on the gaps (toroidal and poloidal) is assumed.

chapters/chapter4/divertor.tex

@@ -69,7 +69,7 @@ \subsubsection{\gls{solps} runs}
 Several \gls{iter} cases were taken from \sidecite{pitts_physics_2019} with \gls{divertor} neutral pressures varying from \SI{1.8}{Pa} to \SI{11.2}{Pa}.
 These \gls{solps} \sidecite{kaveeva_solps-iter_2020} scenarios can be found in the \gls{iter} \gls{imas} database \sidecite{imbeaux_design_2015, park_assessment_2020}.
 The nine simulations used in this work are labelled 122396, 122397, 122398, 122399, 122400, 122401, 122402, 122403 and 122404.
-These have been run in baseline burning plasma conditions ($Q=10$ with \SI{50}{MW} of input power).
+These have been run in baseline detached burning plasma conditions ($Q=10$ with \SI{50}{MW} of input power).
 
 
 \begin{figure}[h!]
@@ -152,6 +152,23 @@ \subsection{Estimation of exposure conditions}
 
 \section{ITER results}
 
+
+\begin{figure*}[h!]
+    \captionsetup[subfigure]{format=plain,singlelinecheck=true}  % needed to center the subcaptions
+    \centering
+    \begin{subfigure}{0.40\linewidth}
+        \includegraphics[width=\linewidth]{Figures/Chapter4/ITER/exposure_conditions_inner.pdf}
+        \caption{\gls{ivt}.}
+    \end{subfigure}%
+    \begin{subfigure}{0.58\linewidth}
+        \includegraphics[width=\linewidth]{Figures/Chapter4/ITER/exposure_conditions_outer.pdf}
+        \caption{\gls{ovt}.}
+    \end{subfigure}
+    \caption{Distribution of exposure conditions (heat flux, particle flux and particle incident energy) along the ITER divertor for several divertor neutral pressures \cite{pitts_physics_2019}.}
+    \labfig{distrib exposure conditions iter divertor}
+\end{figure*}
+
+
 \begin{figure*}[h!]
     \captionsetup[subfigure]{format=plain,singlelinecheck=true}  % needed to center the subcaptions
     \centering
@@ -198,6 +215,7 @@ \section{ITER results}
     \labfig{iter vs time}
 \end{figure}
 
+The exposure conditions vary with respect to the divertor neutral pressure (see \reffig{distrib exposure conditions iter divertor}).
 
 % peak temperature
 Peak temperatures at \glspl{strike point} increased when decreasing the \gls{divertor} neutral pressure (see \reffig{distrib outer target}).

chapters/chapter5/He_transport_in_PFCs/indirect_implantation.tex

@@ -97,7 +97,7 @@ \subsection{Comparison to direct implantation}
 \begin{equation}
     \Gamma = \varphi_\mathrm{imp} \, f(x)
 \end{equation}
-where $\varphi_\mathrm{imp} = \SI{5e25}{m^{-2}.s^{-1}}$ is the implanted helium flux and $f(x)$ is a gaussian distribution centered on $R_p=\SI{1.5}{nm}$ and a width $\sigma=\SI{1.0}{nm}$, which correspond to typical implantation parameters for helium exposure in tokamaks.
+where $\varphi_\mathrm{imp} = \SI{1e23}{m^{-2}.s^{-1}}$ is the implanted helium flux and $f(x)$ is a gaussian distribution centered on $R_p=\SI{1.5}{nm}$ and a width $\sigma=\SI{1.0}{nm}$, which correspond to typical implantation parameters for helium exposure in tokamaks.
 
 When comparing the production of helium from indirect sources with the quantity helium implanted from the plasma $\Gamma$, it appears that the indirect sources are negligible in the exposed region (see \reffig{comparison helium generation}).
 Indirect helium production may become dominant in bulk regions - though may not necessarily be enough to produce helium bubbles.

scripts/he_generation.py

@@ -8,7 +8,7 @@
 gaussian_distribution = (
     1 / (width * (2 * np.pi) ** 0.5) * np.exp(-0.5 * ((x - imp_depth) / width) ** 2)
 )  # m-1
-flux = 5e25  # He/m2/s
+flux = 1e23  # He/m2/s
 
 gaussian_area = np.trapz(gaussian_distribution, x)
 
@@ -39,7 +39,7 @@
     plt.yscale("log")
     matplotx.ylabel_top("He generation \n (m$^{-3}$ s$^{-1}$)")
     plt.xlabel("Depth (nm)")
-    plt.ylim(bottom=1e10, top=1e35)
+    plt.ylim(bottom=1e10, top=1e34)
     matplotx.line_labels(min_label_distance=5)
     plt.tight_layout()
     plt.show()

scripts/trap_occupancy_check.py

@@ -0,0 +1,118 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+k_B = 8.617e-5  # ev/K
+
+
+def trap_occupancy(p_0, E_p, k_0, E_k, c_m, T):
+    p = p_0 * np.exp(-E_p / k_B / T)
+    k = k_0 * np.exp(-E_k / k_B / T)
+    return 1 / (p / k / c_m + 1)
+
+
+T = np.linspace(300, 1200, num=300)
+c_m = 1e22  # m-3
+
+# trap 1
+trap_1_occupancy = trap_occupancy(
+    p_0=1e13, E_p=0.87, k_0=8.96e-17, E_k=0.2, c_m=c_m, T=T
+)
+trap_1_density = 1.1e-3  # at.fr
+trap_1_concentration = trap_1_density * trap_1_occupancy
+
+trap_2_occupancy = trap_occupancy(
+    p_0=1e13, E_p=1.00, k_0=8.96e-17, E_k=0.2, c_m=c_m, T=T
+)
+trap_2_density = 4.0e-4  # at.fr
+trap_2_concentration = trap_2_density * trap_2_occupancy
+
+# equivalent trap
+trap_eq_density = trap_1_density + trap_2_density  # at.fr
+
+trap_eq_occupancy = trap_occupancy(
+    p_0=1e13,
+    E_p=(0.87 * trap_1_density + 1 * trap_2_density) / trap_eq_density,
+    k_0=8.96e-17,
+    E_k=0.2,
+    c_m=c_m,
+    T=T,
+)
+trap_eq_concentration = trap_eq_density * trap_eq_occupancy
+
+# compute difference
+difference = 100 * trap_2_concentration / (trap_1_concentration + trap_2_concentration)
+
+print("Maximum error induced by neglecting trap 2 is {:.0f} %".format(difference.max()))
+
+
+# PLOTTING
+
+fig, axs = plt.subplots(2, 1, sharex=True)
+
+plt.sca(axs[0])
+
+plt.fill_between(T, 0, trap_1_concentration, label="Trap 1", alpha=0.8)
+plt.fill_between(
+    T,
+    trap_1_concentration,
+    trap_1_concentration + trap_2_concentration,
+    label="Trap 2",
+    alpha=0.8,
+)
+# plt.plot(T, trap_eq_concentration, label="Trap eq")
+# plt.annotate("Equivalent trap", (470, 0.0013), color="tab:blue")
+
+plt.annotate("Trap 1", (320, 0.0005), color="white")
+plt.annotate("Trap 2", (320, 0.0012), color="white")
+
+plt.ylabel("Trap concentration (at.fr.)")
+plt.ylim(bottom=0)
+
+
+plt.sca(axs[-1])
+plt.plot(T, difference)
+plt.ylim(bottom=0, top=80)
+
+plt.xlabel("T (K)")
+plt.ylabel("Error (%) \n (neglecting trap 2)")
+plt.tight_layout()
+
+for ax in axs:
+    ax.spines.right.set_visible(False)
+    ax.spines.top.set_visible(False)
+
+plt.show()
+
+
+# 2D map
+
+c_m = np.logspace(20, 23)
+T = np.linspace(300, 1200)
+
+TT, cc = np.meshgrid(T, c_m)
+
+# trap 1
+trap_1_occupancy = trap_occupancy(
+    p_0=1e13, E_p=0.87, k_0=8.96e-17, E_k=0.2, c_m=cc, T=TT
+)
+trap_1_concentration = trap_1_density * trap_1_occupancy
+
+trap_2_occupancy = trap_occupancy(
+    p_0=1e13, E_p=1.00, k_0=8.96e-17, E_k=0.2, c_m=cc, T=TT
+)
+trap_2_concentration = trap_2_density * trap_2_occupancy
+
+
+# compute difference
+difference = 100 * trap_2_concentration / (trap_1_concentration + trap_2_concentration)
+
+CF = plt.contourf(TT, cc, difference, levels=100)
+plt.yscale("log")
+plt.colorbar(label="Error (%)")
+plt.xlabel("T (K)")
+plt.ylabel("Mobile concentration (m$^{-3}$)")
+
+for c in CF.collections:
+    c.set_edgecolor("face")
+
+plt.show()

scripts/zoom_out.py

@@ -0,0 +1,54 @@
+import numpy as np
+import matplotlib as mpl
+import matplotlib.pyplot as plt
+import matplotlib.animation
+import matplotx
+from matplotlib.ticker import FormatStrFormatter
+
+nb_frames_still_beg = 200
+nb_frames_still_end = 100
+
+y_lim_still_beginning = np.ones(nb_frames_still_beg) * 30
+y_lim = np.linspace(30, 750, num=50)
+y_lim_still = np.ones(nb_frames_still_end) * 750
+
+y_lim = np.concatenate([y_lim_still_beginning, y_lim, y_lim_still])
+
+x_lim_still_beginning = np.ones(nb_frames_still_beg) * 2
+x_lim = np.linspace(2, 5, num=50)
+x_lim_still = np.ones(nb_frames_still_end) * 5
+x_lim = np.concatenate([x_lim_still_beginning, x_lim, x_lim_still])
+
+with plt.style.context(matplotx.styles.dufte_bar):
+    fig = plt.figure()
+    ax = fig.add_subplot(1, 1, 1)
+    pos = [-1.5, 0, 1.5]
+    labels = ["calculated", r"80 % error", r"50 % Tritium"]
+    (bar) = plt.bar(pos, [14, 14 * 1.8, 14 * 1.8 * 0.5])
+    plt.xticks(pos, labels)
+    ax.set_ylim(0, y_lim[0])
+    ax.set_xlim(-x_lim[0], x_lim[0])
+    # matplotx.show_bar_values("{:.0f} g")
+    plt.hlines(700, -4, 4, colors="tab:green", linestyles="dashed")
+    plt.annotate("ITER safety limit", xy=(0, 725), color="tab:green", ha="center")
+    plt.gca().yaxis.set_major_formatter(FormatStrFormatter("%d g"))
+
+
+def update(i):
+    ax.set_ylim(0, y_lim[i])
+    ax.set_xlim(-x_lim[i], x_lim[i])
+    if i > nb_frames_still_beg:
+        for child in ax.get_children():
+            if isinstance(child, matplotlib.text.Text):
+                if child.get_text() != "ITER safety limit":
+                    child.set(visible=False)
+    if i > nb_frames_still_beg + 15:
+        plt.xticks([])
+
+
+ani = mpl.animation.FuncAnimation(
+    fig, update, frames=y_lim.size, blit=False, repeat=False
+)
+plt.tight_layout()
+ani.save("animation_divertor_inventories.gif", writer="imagemagick", fps=60)
+# plt.show()