diff --git a/doc/content/examples/fuel_cycle/index.md b/doc/content/examples/fuel_cycle/index.md index 35008286..824d77fe 100644 --- a/doc/content/examples/fuel_cycle/index.md +++ b/doc/content/examples/fuel_cycle/index.md @@ -1,6 +1,6 @@ # Fuel Cycle -This demonstration re-creates the tritium fuel cycle model described by [!cite](Abdou2021). +This demonstration re-creates the tritium fuel cycle model described by [!cite](Abdou2021). ### Generating the Input File @@ -144,8 +144,8 @@ lines are model results. ### Python-based Interactive Script -A python-based graphical interactive script is available at [/test/tests/fuel_cycle/fuel_cycle_gui.py](/scripts/fuel_cycle_gui.py) as a demonstration of the various effects of -individual parameters. To use it, navigate to the scripts directory and run the script. All of the input parameters for the model can be changed by editing the associated entry, +A python-based graphical interactive script is available at [/test/tests/fuel_cycle/fuel_cycle_gui.py](/scripts/fuel_cycle_gui.py) as a demonstration of the various effects of +individual parameters. To use it, navigate to the scripts directory and run the script. All of the input parameters for the model can be changed by editing the associated entry, then clicking the "Run" button. Once the simulation has run, checkboxes will appear for each system-level tritium inventory. Time units can also be adjusted by selecting the appropriate timescale. diff --git a/doc/content/verification_and_validation/val-2a.md b/doc/content/verification_and_validation/val-2a.md index db1382e8..ed77da86 100644 --- a/doc/content/verification_and_validation/val-2a.md +++ b/doc/content/verification_and_validation/val-2a.md @@ -50,9 +50,9 @@ J = 2 A K_r C^2. ## Case and Model Parameters -The beam flux on the upstream side of the sample during the experiment is presented in [val-2a_flux_and_pressure_TMAP4], and only 75 % of the flux remain in the sample. Other case and model parameters used in TMAP8 are listed in [val-2a_set_up_values_TMAP4]. +The beam flux on the upstream side of the sample during the experiment is presented in [val-2a_flux_and_pressure_TMAP4], and only 75% of the flux remains in the sample. Other case and model parameters used in TMAP8 are listed in [val-2a_set_up_values_TMAP4]. -!table id=val-2a_flux_and_pressure_TMAP4 caption=Values of beam flux on the upstream side of the sample during the experiment [!citep](longhurst1992verification). +!table id=val-2a_flux_and_pressure_TMAP4 caption=Values of beam flux on the upstream side of the sample during the experiment [!citep](anderl1985tritium,longhurst1992verification). | time (s) | Beam flux $F$ (atom/m$^2$/s) | | --------- | ---------------------------- | | 0 - 5820 | 4.9$\times 10^{19}$ | @@ -63,7 +63,7 @@ The beam flux on the upstream side of the sample during the experiment is presen | 17678 - 20000 | 0 | !alert warning title=Typo in [!cite](longhurst1992verification) -The times listed in [!cite](longhurst1992verification) for TMAP8 for the start and end times of the beam are not accurate. instead, TMAP8 uses the times directly from [!cite](anderl1985tritium) to better correspond to experimental conditions. +The times listed in [!cite](longhurst1992verification) for TMAP8 for the start and end times of the beam are not accurate. Instead, TMAP8 uses the times directly from [!cite](anderl1985tritium) to better correspond to experimental conditions. !table id=val-2a_set_up_values_TMAP4 caption=Values of material properties. Note that $K_d$ are currently not used in the input file since the upstream and downstream pressure do not noticeably influence the results. | Parameter | Description | Value | Units | Reference | diff --git a/doc/content/verification_and_validation/val-2b.md b/doc/content/verification_and_validation/val-2b.md index f40129a8..70c1e187 100644 --- a/doc/content/verification_and_validation/val-2b.md +++ b/doc/content/verification_and_validation/val-2b.md @@ -6,8 +6,8 @@ This validation problem is taken from [!cite](macaulay1991deuterium) and is part of the validation suite of TMAP4 and TMAP7 [!citep](longhurst1992verification,ambrosek2008verification), which we reproduce here, with some updates. -R.G. Macaulay-Newcombe et al. conducted thermal absorption and desorption experiments, as well as implantation experiments, on wafers of polished beryllium. -Of the several data sets presented, the one modeled here is run 2a1 represented in Figure 2(a) in their publication. The beryllium sample was 0.4 mm thick and had an area of 104 mm$^2$, as illustrated in [val-2b_schematic]. +R.G. Macaulay-Newcombe et al. conducted thermal absorption and desorption experiments, as well as implantation experiments, on wafers of polished beryllium [!citep](macaulay1991deuterium). +Of the several data sets presented, the one modeled here is titled "run 2a1" and is represented in Figure 2(a) in their publication [!citep](macaulay1991deuterium). The beryllium sample was 0.4 mm thick and had an area of 104 mm$^2$, as illustrated in [val-2b_schematic]. It was polished to a mirror finish and then exposed to 13.3 kPa of deuterium at 773 K for 50 hours. It was quickly cooled under a vacuum of about 1 $\mu$Pa. The cooling time constant for the apparatus is taken as 45 minutes, which is consistent with the assumption made in [!citep](longhurst1992verification,ambrosek2008verification). After removing the sample from the charging furnace, it was transferred in the air to a thermal desorption furnace where the temperature was increased from ambient (300 K) to 1073 K at the rate of 3 K/min. This was done under vacuum, and the pressure of the chamber was monitored by the residual gas analysis and calibrated against standard leaks. In that way, the emission rate from the sample could be measured as a function of temperature. The sample pressure and temperature histories are shown in [val-2b_temperature_pressure_history]. @@ -25,7 +25,7 @@ Experimental data from that measurement, given in Figure 2 (a) in [!cite](macaul caption=Pressure and temperature histories. !alert note title=Uncertainty about cooldown duration. -The exact duration of the cooldown period and its temperature history are uncertain. [!cite](macaulay1991deuterium) provides information about 24 hour-cooldown cycles, but it is unclear whether this applies to the charging chamber alone or to the sample as well. In parallel, [!citep](longhurst1992verification,ambrosek2008verification) assumes that the cooldown lasted for 40 minutes. With an assumed cooling time constant for the apparatus of 45 minutes, this did not enable the sample to cool down to the starting temperature of the desorption phase of the experiment (i.e., 300 K). To model this case in TMAP8, we decided to select a cooldown duration that is long enough to bring the temperature of the sample to around 300 K, but did not unnecessarily increase the length of the history since no significant changes happen to the deuterium distribution at 300 K due to slow kinetics. For these reasons, we selected a cooldown duration of 5 hours, as shown in [val-2b_temperature_pressure_history]. +The exact duration of the cooldown period and its temperature history are uncertain. [!cite](macaulay1991deuterium) provides information about 24 hour-cooldown cycles, but it is unclear whether this applies to the charging chamber alone or to the sample as well. In parallel, [!citep](longhurst1992verification,ambrosek2008verification) assume that the cooldown lasted for 40 minutes. With an assumed cooling time constant for the apparatus of 45 minutes, this did not enable the sample to cool down to the starting temperature of the desorption phase of the experiment (i.e., 300 K). To model this case in TMAP8, we decided to select a cooldown duration that is long enough to bring the temperature of the sample to around 300 K, but did not unnecessarily increase the length of the history since no significant changes happen to the deuterium distribution at 300 K due to slow kinetics. For these reasons, we selected a cooldown duration of 5 hours, as shown in [val-2b_temperature_pressure_history]. From Rutherford backscattering measurements made on the samples before charging with deuterium, they deduced that the thickness of the oxide film was 18 nm. This is typical for polished beryllium. The metal is so reactive in air that the film forms almost immediately after any surface oxide removal. On the other hand, it is relatively stable and would only grow slightly when exposed to air between charging and thermal desorption. @@ -33,11 +33,16 @@ This experiment is modeled using a two-segment model in TMAP8 with the segments The diffusivity of deuterium in beryllium was measured by [!cite](abramov1990deuterium). They made measurements on high-grade (99$\%$ pure) and extra-grade (99.8$\%$ pure). The values used here are those for high-grade beryllium, consistent with Dr. Macaulay-Newcombe's measurements of the purity of his samples. -Deuterium transport properties of the BeO are more challenging. First, it is not clear in which state the deuterium exists in the BeO. However, it has been observed [!cite](longhurst1990tritium) that an activation energy of -78 kJ/mol (exothermic reaction) is evident for tritium coming out of neutron-irradiated beryllium in work done by D. L. Baldwin of Pacific Northwest Laboratory. The same value of energy has appeared in other results (can be inferred from Dr. Swansiger's work cited by [!cite](wilson1990beryllium) and by [!cite](causey1990tritium), among others), so one may be justified in using it. Concerning the solubility, measurements reported by [!cite](macaulay1992thermal) and in follow-up conversations indicate about 200 appm of D in BeO after exposure to 13.3 kPa of D$_2$ at 773 K. That suggests a coefficient of only 1.88 x 10$^{18}$ d/m$^3$Pa$^{1/2}$. Since much of the deuterium in the oxide layer will get out during the cool-down process (and because it gives a good fit), the solubility coefficient is taken to be 5 $\times$ 10$^{20}$ d/m$^3$Pa$^{1/2}$. +Deuterium transport properties of the BeO are more challenging. First, it is not clear in which state the deuterium exists in the BeO. +However, it has been observed [!cite](longhurst1990tritium) that an activation energy of -78 kJ/mol (exothermic solution) is evident for tritium coming out of neutron-irradiated beryllium in work done by D. L. Baldwin of Pacific Northwest Laboratory. +The same value of energy has appeared in other results. It can be inferred from Dr. Swansiger's work cited by [!cite](wilson1990beryllium) and by [!cite](causey1990tritium), among others, so one may be justified in using it. +Concerning the solubility, measurements reported by [!cite](macaulay1992thermal) and in follow-up conversations indicate about 200 appm of D in BeO after exposure to 13.3 kPa of D$_2$ at 773 K. +That suggests a coefficient of only 1.88 $\times$ 10$^{18}$ d/m$^3$Pa$^{1/2}$. +Since much of the deuterium in the oxide layer will get out during the cool-down process (and because it gives a good fit), the solubility coefficient is taken to be 5 $\times$ 10$^{20}$ d/m$^3$Pa$^{1/2}$. -Deuterium diffusion measurements in BeO were made by [!cite](fowler1977tritium). They found a wide range of results for diffusivity in BeO depending on the physical form of the material, having measured it for single-crystal, sintered, and powdered BeO. The model in [!citep](longhurst1992verification,ambrosek2008verification) uses one expression for the charging phase and another for the thermal desorption phase, believing that the surface film changed somewhat during the transfer between the two furnaces. For the charging phase diffusivity, the model uses 20 times that for the sintered BeO. Thermal expansion mismatches tend to open up cracks and channels in the oxide layer, so this seems a reasonable value. The same activation energy of 48.5 kJ/mol, is retained, however. For the thermal desorption phase, the diffusivity prefactor of the sintered material (7x10$^{-5}$ m$^2$/sec) and an activation energy of 223.7 kJ/mol (53.45 kcal/mol) are used. These values give good results and lie well within the scatter of Fowler's data. Exposure of the sample to air after heating should have made the oxide more like a single crystal by healing the cracks that may have developed. Diffusivities and solubilities used in the simulation are listed in [val-2b_parameters]. +Deuterium diffusion measurements in BeO were made by [!cite](fowler1977tritium). They found a wide range of results for diffusivity in BeO depending on the physical form of the material, having measured it for single-crystal, sintered, and powdered BeO. The model in [!citep](longhurst1992verification,ambrosek2008verification) uses one expression for the charging phase and another for the thermal desorption phase, believing that the surface film changed somewhat during the transfer between the two furnaces. For the charging phase diffusivity, the model uses 20 times that for the sintered BeO. Thermal expansion mismatches tend to open up cracks and channels in the oxide layer, so this seems a reasonable value. The same activation energy of 48.5 kJ/mol, is retained, however. For the thermal desorption phase, the diffusivity prefactor of the sintered material (7x10$^{-5}$ m$^2$/sec) and an activation energy of 223.7 kJ/mol (53.45 kcal/mol) are used. These values give results lying well within the scatter of Fowler's data. Exposure of the sample to air after heating should have made the oxide more like a single crystal by healing the cracks that may have developed. Diffusivities and solubilities used in the simulation are listed in [val-2b_parameters]. -!table id=val-2b_parameters caption=Model parameter values for the charging and the desorption phases [!citep](longhurst1992verification,ambrosek2008verification). T is the temperature in Kelvin. +!table id=val-2b_parameters caption=Model parameter values for the charging and the desorption phases [!citep](longhurst1992verification,ambrosek2008verification). $T$ is the temperature in Kelvin. | Property of deuterium | Value for charging phase | Value for desorption phase | Units | | --------------------- | ------------------------------------- | ------------------------------------- | ------------------- | | Diffusivity in Be | $8.0 \times 10^{-9} \exp(-4220/T)$ | $8.0 \times 10^{-9} \exp(-4220/T)$ | m$^2$/s | @@ -45,7 +50,7 @@ Deuterium diffusion measurements in BeO were made by [!cite](fowler1977tritium). | Solubility in Be | $7.156 \times 10^{27} \exp(-11606/T)$ | $7.156 \times 10^{27} \exp(-11606/T)$ | at/m$^3$/Pa$^{1/2}$ | | Solubility in BeO | $5.00 \times 10^{20} \exp(9377.7/T)$ | $5.00 \times 10^{20} \exp(9377.7/T)$ | at/m$^3$/Pa$^{1/2}$ | -The model applies 13.3 kPa of D$_2$ for 50 hours and 15 seconds followed by cool down with a 45 minute time constant at 1 $\mu$Pa for five hour ([!cite](longhurst1992verification,ambrosek2008verification) used 40 minutes, but we extend it here to 5 hours to let the temperature go down closer to 300 K - see [val-2b_temperature_pressure_history]). The deuterium concentrations in the sample have a complex distribution that results from first charging the sample and then discharging it during the cool down. This problem is then restarted with different model parameters (see [val-2b_parameters]) to simulate thermal desorption in a $1 \times 10^{-3}$ Pa environment. That begins at 300 K and goes to 1073 K at 3 K/min. +The model applies 13.3 kPa of D$_2$ for 50 hours and 15 seconds followed by cool down with a 45 minute time constant at 1 $\mu$Pa for five hour ([!cite](longhurst1992verification,ambrosek2008verification) used 40 minutes, but we extend it here to 5 hours to let the temperature go down closer to 300 K - see [val-2b_temperature_pressure_history]). The deuterium concentrations in the sample have a complex distribution that results from first charging the sample and then discharging it during the cool down. This problem is then restarted with different model parameters (see [val-2b_parameters]) to simulate thermal desorption in a $1 \times 10^{-3}$ Pa environment that begins at 300 K and goes to 1073 K at 3 K/min. ## Results @@ -55,10 +60,10 @@ The model applies 13.3 kPa of D$_2$ for 50 hours and 15 seconds followed by cool image_name=val-2b_comparison.png style=width:50%;margin-bottom:2%;margin-left:auto;margin-right:auto id=val-2b_comparison - caption=Comparison of TMAP8 calculation with the experimental data, which shows TMAP8's ability to accurately model this validation case. + caption=Comparison of TMAP8 calculation against experimental data, which shows TMAP8's ability to accurately model this validation case. !alert note title=Experimental data from [!cite](macaulay1991deuterium). -The experimental data used in this case comes directly from Figure (2) in [!cite](macaulay1991deuterium). This is in contrast with the data used in [!cite](longhurst1992verification,ambrosek2008verification), which, although the scale of the data is very similar, differs slightly. Note also that the units in Figure (2) from [!cite](macaulay1991deuterium) should be atoms/mm$^2$/s $\times 10^{10}$ instead of atoms/mm$^2$ $\times 10^{10}$, which is corrected in [val-2b_comparison]. +The experimental data used in this case comes directly from Figure (2) in [!cite](macaulay1991deuterium). This is in contrast with the data used in [!cite](longhurst1992verification,ambrosek2008verification), which, although the scale of the data is very similar, differs slightly. Also, note that the units in Figure (2) from [!cite](macaulay1991deuterium) should be atoms/mm$^2$/s $\times 10^{10}$ instead of atoms/mm$^2$ $\times 10^{10}$, which is corrected in [val-2b_comparison]. To verify that the solubility ratio at the interface between the beryllium and its oxide is appropriately modeled, [val-2b_ratio] compares the known solubility ratio with the calculated deuterium concentration ratio at the interface, and they match, as expected. diff --git a/doc/content/verification_and_validation/val-2c.md b/doc/content/verification_and_validation/val-2c.md index a42ce5be..d3ce64e1 100644 --- a/doc/content/verification_and_validation/val-2c.md +++ b/doc/content/verification_and_validation/val-2c.md @@ -9,9 +9,8 @@ Whenever tritium is released into a fusion reactor test cell, it is crucial to c This case models an experiment conducted at Los Alamos National Laboratory at the tritium systems test assembly (TSTA) to study the behavior of tritium once released in a test cell and the efficacy of the emergency tritium cleanup system (ETCS). The experimental set up, described in greater detail in [!cite](Holland1986), can be summarized as such: -the inner walls of an enclosure of volume $V$ was covered with an aluminum foil and then covered in paint with an average thickness $l$, which is then in contact with the enclosure air. -A given amount, T$_2^0$, of tritium, T$_2$, is injected in the enclosure, which initially contained ambient air. -This represents the tritium release. +the inner walls of an enclosure of volume $V$ are covered with an aluminum foil and then covered in paint with an average thickness $l$, which is then in contact with the enclosure air. +A given amount, T$_2^0$, of tritium, T$_2$, is injected in the enclosure, which initially contained ambient air, representing tritium release. A flow rate $f$ through the enclosure represents the air replacement time expected for large test cells. The purge gas is ambient air with 20% relative humidity. A fraction of that amount is diverted through the measurement system to determine the concentrations of chemical species within the enclosure. @@ -36,7 +35,7 @@ Here, $c_i$ represents the concentration of species $i$, and $K^0$ is a constant Second, the different species will permeate in the paint. The elemental tritium species, T$_2$ and HT, have a given solubility $K_S^e$ and diffusivity $D^e$, while the tritiated water, HTO, and water, H$_2$O, have a solubility $K_S^w$ and diffusivity $D^w$. -It is expected that species will initially permeate into the paint and later get released as the purge gas cleans up the enclosure air. +It is expected that the species will initially permeate into the paint and later get released as the purge gas cleans up the enclosure air. The objectives of this case are to determine the time evolution of T$_2$ and HTO concentrations in the enclosure, match the experimental data published in [!cite](Holland1986), and display this comparison with the appropriate error checking (see [val-2c_comparison_T2] and [val-2c_comparison_HTO]). @@ -44,7 +43,7 @@ The objectives of this case are to determine the time evolution of T$_2$ and HTO To model the case described above, TMAP8 simulates a one-dimensional domain with one block to represent the air in the enclosure, and another block to represent the paint. In each block, the simulation tracks the local concentration of T$_2$, HT, HTO, and H$_2$O. -Note that this case can easily be extended to a two- or three-dimensional case, but, consistent with previous analyses, we will maintain the one-dimensional configuration here. +Note that this case can easily be extended to a two- or three-dimensional case, but consistent with previous analyses, we will maintain the one-dimensional configuration here. In the enclosure, to capture the purge gas and the chemical reactions, the concentrations evolve as \begin{equation} \label{eq:enclosure:T2} @@ -64,17 +63,17 @@ where $c_{\text{H}_2\text{O}}^0$ is the concentration of H$_2$O in the incoming In the paint, TMAP8 captures species diffusion through \begin{equation} \label{eq:paint:T2} -\frac{d c_{\text{T}_2}}{dt} = - \nabla D^e \nabla c_{\text{T}_2}, +\frac{d c_{\text{T}_2}}{dt} = \nabla D^e \nabla c_{\text{T}_2}, \end{equation} \begin{equation} \label{eq:paint:HT} -\frac{d c_{\text{HT}}}{dt} = - \nabla D^e \nabla c_{\text{HT}}, +\frac{d c_{\text{HT}}}{dt} = \nabla D^e \nabla c_{\text{HT}}, \end{equation} \begin{equation} \label{eq:paint:HTO} -\frac{d c_{\text{HTO}}}{dt} = - \nabla D^w \nabla c_{\text{HTO}}, +\frac{d c_{\text{HTO}}}{dt} = \nabla D^w \nabla c_{\text{HTO}}, \end{equation} and \begin{equation} \label{eq:paint:H2O} -\frac{d c_{\text{H}_2\text{O}}}{dt} = - \nabla D^w \nabla c_{\text{H}_2\text{O}}. +\frac{d c_{\text{H}_2\text{O}}}{dt} = \nabla D^w \nabla c_{\text{H}_2\text{O}}. \end{equation} At the interface between the enclosure air and the paint, sorption is captured in TMAP8 with Henry's law thanks to the [InterfaceSorption.md] object: @@ -111,8 +110,8 @@ The case and model parameters used in both approaches in TMAP8 are listed in [va | $K^0$ | $\boldsymbol{1.5} \times 2.0 \times 10^{-10}$ | $\boldsymbol{2.8} \times 2.0 \times 10^{-10}$ | m$^3$/Ci/s | Adapted from [!cite](longhurst1992verification) | | $D^e$ | 4.0 $\times 10^{-12}$ | Identical | m$^2$/s | [!cite](Holland1986) | | $D^w$ | 1.0 $\times 10^{-14}$ | Identical | m$^2$/s | [!cite](Holland1986) | -| $K_S^e$ | $\boldsymbol{5.0 \times 10^{-2}} \times 4.0 \times 10^{-19}$ | $\boldsymbol{1.0 \times 10^{-3}} \times 4.0 \times 10^{-19}$ | 1/m$^3$/Pa | Adapted from [!cite](longhurst1992verification) | -| $K_S^w$ | $\boldsymbol{3.5 \times 10^{-4}} \times 6.0 \times 10^{-24}$ | $\boldsymbol{3.0 \times 10^{-4}} \times 6.0 \times 10^{-24}$ | 1/m$^3$/Pa | Adapted from [!cite](longhurst1992verification) | +| $K_S^e$ | $\boldsymbol{5.0 \times 10^{-2}} \times 4.0 \times 10^{19}$ | $\boldsymbol{1.0 \times 10^{-3}} \times 4.0 \times 10^{19}$ | 1/m$^3$/Pa | Adapted from [!cite](longhurst1992verification) | +| $K_S^w$ | $\boldsymbol{3.5 \times 10^{-4}} \times 6.0 \times 10^{24}$ | $\boldsymbol{3.0 \times 10^{-4}} \times 6.0 \times 10^{24}$ | 1/m$^3$/Pa | Adapted from [!cite](longhurst1992verification) | | $t_{injection}$ | N/A | 3 | hr | | ## Results and discussion @@ -132,7 +131,7 @@ It is also possible to perform this optimization with [MOOSE's stochastic tools image_name=val-2c_comparison_TMAP8_Exp_T2_Ci.png style=width:50%;margin-bottom:2%;margin-left:auto;margin-right:auto id=val-2c_comparison_T2 - caption=Comparison of TMAP8 calculations against the experimental data for T$_2$ concentration in the enclosure over time. TMAP8 matches the experimental data well, with an improvement when T$_2$ is injected over a given period rather than immediately. + caption=Comparison of TMAP8 calculations against the experimental data for T$_2$ concentration in the enclosure over time. TMAP8 matches the experimental data well, with an improvement when T$_2$ is injected over a given period rather than immediately. !media comparison_val-2c.py image_name=val-2c_comparison_TMAP8_Exp_HTO_Ci.png diff --git a/doc/content/verification_and_validation/ver-1a.md b/doc/content/verification_and_validation/ver-1a.md index ed564b5a..51d0ec5d 100644 --- a/doc/content/verification_and_validation/ver-1a.md +++ b/doc/content/verification_and_validation/ver-1a.md @@ -7,23 +7,23 @@ This verification case consists of an enclosure containing a finite concentration of atoms which are allowed to thermally diffuse through a SiC layer over time. No Soret effects, solubility, or trapping effects are included. -This is one of the original problems introduced in [!cite](longhurst1992verification) for TMAP4 and adapted in [!cite](ambrosek2008verification) for TMAP7. Note, however, that the verification cases for TMAP4 and TMAP7, although using the exact same set up, use different quantities to verify their implementation (see [ver-1a_comparison_analytical_TMAP4_release_fraction], [ver-1a_comparison_analytical_TMAP7_release_fraction], and [ver-1a_comparison_analytical_TMAP7_flux]). In TMAP8, for completeness, we perform verification on all these quantities and show agreement with analytical solutions from both TMAP4 and TMAP7. +This is one of the original problems introduced in [!cite](longhurst1992verification) for TMAP4 and adapted in [!cite](ambrosek2008verification) for TMAP7. Note, however, that the verification cases for TMAP4 and TMAP7, although using the exact same setup, use different quantities to verify their implementation (see [ver-1a_comparison_analytical_TMAP4_release_fraction], [ver-1a_comparison_analytical_TMAP7_release_fraction], and [ver-1a_comparison_analytical_TMAP7_flux]). In TMAP8, for completeness, we perform verification on all these quantities and show agreement with analytical solutions from both TMAP4 and TMAP7. ## Case Set up -[ver-1a_schematic] shows a schematic of the verification 1a case, along with an illustration of the quantities used for verification. It consists of an enclosure that is pre-charged with a fixed quantity of tritium in gaseous form, and a finite SiC slab. At time t > 0, the tritium is allowed to thermally diffuse through the SiC slab, initially at zero concentration. The surface of the slab in contact with the source is assumed to be in equilibrium with the source enclosure, assuming a temperature-dependent solubility $S$. As a boundary condition, the concentration at the outer surface of the SiC slab is kept at zero for all time. The diffusion of the tritium through the SiC slab is calculated by neglecting Soret and trapping effects. Different aspects of the diffusion are compared against analytical solutions, as described below. The values used to characterize the necessary material properties and the case geometry are provided in [ver-1a_set_up_values]. +[ver-1a_schematic] shows a schematic of the ver-1a case, along with an illustration of the quantities used for verification. It consists of an enclosure that is pre-charged with a fixed quantity of tritium in gaseous form and a finite SiC slab. At time t > 0, the tritium is allowed to thermally diffuse through the SiC slab, initially at zero concentration. The surface of the slab in contact with the source is assumed to be in equilibrium with the source enclosure, assuming a temperature-dependent solubility $S$. As a boundary condition, the concentration at the outer surface of the SiC slab is kept at zero for all time. The diffusion of the tritium through the SiC slab is calculated by neglecting Soret and trapping effects. Different aspects of the diffusion are compared against analytical solutions, as described below. The values used to characterize the necessary material properties and the case geometry are provided in [ver-1a_set_up_values]. !media figures/ver-1a_schematic.png style=width:50%;margin-bottom:2%;margin-left:auto;margin-right:auto id=ver-1a_schematic - caption= Schematic of the verification 1a case and illustration of the quantities used for verification. + caption= Schematic of the ver-1a case and illustration of the quantities used for verification. !table id=ver-1a_set_up_values caption=Values of material properties and case geometry with $R$, the gas constant, as defined in [PhysicalConstants](source/utils/PhysicalConstants.md). | Parameter | Description | Value | Units | | --------- | -------------------------- | -------------------------------------- | ---------- | | $V$ | Enclosure volume | 5.20e-11 | m$^3$ | | $A$ | SiC slab surface area | 2.16e-6 | m$^2$ | -| $P_0$ | Initial enclosure pressure | 1e6 | mPa | +| $P_0$ | Initial enclosure pressure | 1e6 | Pa | | $T$ | Temperature | 2373 | K | | $D$ | Tritium diffusivity in SiC | 1.58e-4 $\times \exp(-308000.0/(R*T))$ | m$^2$/2 | | $S$ | Tritium solubility in SiC | 7.244e22 / $T$ | 1/m$^3$/Pa | @@ -49,10 +49,10 @@ with \phi = \frac{\text{source concentration}}{\text{layer concentration}} = \frac{1}{S k_b T}, \end{equation} -where the "source concentration" is the concentration in the enclosure ($P(t)/k_b/T$), $k_b$ is the Boltzmann constant (as defined in [PhysicalConstants](source/utils/PhysicalConstants.md)), and "layer concentration" is the concentration in the slab at the interface with the enclosure ($S P(t)$). $\phi$ is constant in time. $\alpha_n$ are the roots of +where the "source concentration" is the concentration in the enclosure ($P(t)/k_b/T$), $k_b$ is the Boltzmann constant (as defined in [PhysicalConstants](source/utils/PhysicalConstants.md)), and "layer concentration" is the concentration in the slab at the interface with the enclosure ($S P(t)$). $\phi$ is constant in time, and $\alpha_n$ are the roots of \begin{equation} - \alpha_n = \frac{L}{tan \ \alpha_n}. + \alpha_n = \frac{L}{\tan \ \alpha_n}. \end{equation} @@ -68,7 +68,7 @@ The expression in Eq. (1) of [!cite](longhurst1992verification) writes $-\alpha_ image_name=ver-1a_comparison_analytical_TMAP4_release_fraction.png style=width:50%;margin-bottom:2%;margin-left:auto;margin-right:auto id=ver-1a_comparison_analytical_TMAP4_release_fraction - caption=Comparison of TMAP8 calculation with the analytical solution for the release fraction from [!cite](longhurst1992verification). + caption=Comparison of TMAP8 calculation with the analytical solution for the release fraction on the outer surface of the SiC slab [!cite](longhurst1992verification). ## Verification of the release fraction on the inner surface of the SiC slab (TMAP7) @@ -80,7 +80,7 @@ In [!cite](ambrosek2008verification), i.e., TMAP7, the verification test focuses \end{equation} where $P(t)$ is the pressure at the surface of the enclosure over time. -To derive $FR(t)$, [!cite](ambrosek2008verification) first reference the analytical solution for an analogous heat transfer problem [!cite](Carslaw1959conduction), which provides the solute concentration profile in the membrane as +To derive $FR(t)$, [!cite](ambrosek2008verification) first references the analytical solution for an analogous heat transfer problem [!cite](Carslaw1959conduction), which provides the solute concentration profile in the membrane as \begin{equation} C(x,t) = 2 S P_0 L' \sum_{n=1}^{\infty} \frac{\exp \left(-\alpha_{n}^2 D t\right) \sin (\alpha_{n} x)}{[l(\alpha_{n}^2 + L'^2)+L'] \sin (\alpha_{n} l)}, @@ -91,7 +91,7 @@ where $x=l$ is the position of the surface in contact with the enclosure, and $x \end{equation} and $\alpha_n$ are the roots of \begin{equation} - \alpha_n = \frac{L}{tan (\alpha_n l)}. + \alpha_n = \frac{L}{\tan (\alpha_n l)}. \end{equation} By linking the concentration on the surface of the slab in contact with the enclosure ($x=l$ as used in [!cite](ambrosek2008verification)) with the pressure of the enclosure, Henry's law provides @@ -104,7 +104,7 @@ which leads to \end{equation} !alert warning title=Typos in [!cite](ambrosek2008verification) -1. The units and expression of the solubility provided in [!cite](ambrosek2008verification) have typos. The correct values and units are provided above in [ver-1a_set_up_values]. The value provided in [!cite](ambrosek2008verification) is $S=3.053 \times 10^{\bold{\underline{29}}}$ kg $\cdot$ m $^2$/s$^2$ instead of $S=7.244\times 10^{22} / T = 3.053 \times 10^{\bold{\underline{19}}}$ 1/m$^3$/Pa when $T=2373$ K. +1. The units and expression of the solubility provided in [!cite](ambrosek2008verification) have typos. The correct values and units are provided above in [ver-1a_set_up_values]. The value provided in [!cite](ambrosek2008verification) is $S=3.053 \times 10^{\bold{\underline{29}}}$ kg $\cdot$ m $^2$/s$^2$ instead of $S=7.244\times 10^{22} / T = 3.053 \times 10^{\bold{\underline{19}}}$ 1/m$^3$/Pa when $T=2373$ K. The correct value of $S$ is used in the input file decks in [!cite](longhurst1992verification,ambrosek2008verification). 2. The release fraction in [!cite](ambrosek2008verification) is described as $P(t)/P_0$ in their Eq. (5) but is actually plotted as $1-(P(t)/P_0)$ in Figure 1 of that report. 3. Not a typo, per se, but a potential source of confusion for users is that in [!cite](ambrosek2008verification), $x=0$ represents the surface not exposed to the enclosure (where $c=0$), and $x=l$ represents the surface exposed to the enclosure. In the TMAP8 documentation, we have kept this convention to correspond to the TMAP7 case, but note that the TMAP8 input file fixes $x=0$ at the enclosure surface and $x=l$ for the outer surface. @@ -116,7 +116,7 @@ which leads to image_name=ver-1a_comparison_analytical_TMAP7_release_fraction.png style=width:50%;margin-bottom:2%;margin-left:auto;margin-right:auto id=ver-1a_comparison_analytical_TMAP7_release_fraction - caption=Comparison of TMAP8 calculation with the analytical solution for the release fraction from [!cite](ambrosek2008verification). + caption=Comparison of TMAP8 calculation with the analytical solution for the release fraction on the inner surface of the SiC slab [!cite](ambrosek2008verification). ## Verification of the tritium flux at the outer surface of the SiC slab (TMAP7) @@ -140,7 +140,7 @@ Again, be aware of the typos in [!cite](ambrosek2008verification) and the differ image_name=ver-1a_comparison_analytical_TMAP7_flux.png style=width:50%;margin-bottom:2%;margin-left:auto;margin-right:auto id=ver-1a_comparison_analytical_TMAP7_flux - caption=Comparison of TMAP8 calculation with the analytical solution for the release fraction from [!cite](ambrosek2008verification). + caption=Comparison of TMAP8 calculation with the analytical solution for the tritium flux at the outer surface of the SiC slab [!cite](ambrosek2008verification). ## Input files diff --git a/doc/content/verification_and_validation/ver-1b.md b/doc/content/verification_and_validation/ver-1b.md index 6faa6f8b..09887056 100644 --- a/doc/content/verification_and_validation/ver-1b.md +++ b/doc/content/verification_and_validation/ver-1b.md @@ -3,28 +3,28 @@ # Diffusion Problem with Constant Source Boundary Condition This verification problem is taken from [!cite](longhurst1992verification). Diffusion of tritium through a semi-infinite SiC layer is modeled with a constant -source located on one boundary. No solubility or traping is included. The -concentration as a function of time and position is given by: - +source located on one boundary. No solubility or trapping is included. The +concentration as a function of time and position is given by \begin{equation} -C = C_o \; erfc \left(\frac{x}{2\sqrt{Dt}}\right) +C = C_0 \; erfc \left(\frac{x}{2\sqrt{Dt}}\right), \end{equation} +where $C_0$ the constant source concentration, erfc is the error function, $x$ is the distance from the boundary, $D$ is the diffusion coefficient, and $t$ is the time. Comparison of the TMAP8 results and the analytical solution is shown in [ver-1b_comparison_time] as a function of time at -x = 0.2 mm. For simplicity, both the diffusion coefficient and the initial -concentration were set to unity. The TMAP8 code predictions match very well with -the analytical solution. +$x = 0.2$ mm. For simplicity, both the diffusion coefficient and the initial +concentration were set to unity. The TMAP8 code predictions match +the analytical solution very well. !media comparison_ver-1b.py image_name=ver-1b_comparison_time.png style=width:50%;margin-bottom:2%;margin-left:auto;margin-right:auto id=ver-1b_comparison_time - caption=Comparison of concentration as function of time at x\=0.2m calculated - through TMAP8 and analytically + caption=Comparison of concentration as function of time at $x = 0.2$m calculated + through TMAP8 and analytically. As a second check, the concentration as a function of position at a given time -t = 25s, from TMAP8 was compared with the analytical solution as shown in +$t = 25$ s calculated by TMAP8 was compared with the analytical solution as shown in [ver-1b_comparison_dist]. The predicted concentration profile from TMAP8 is in good agreement with the analytical solution. @@ -33,21 +33,20 @@ good agreement with the analytical solution. style=width:50%;margin-bottom:2%;margin-left:auto;margin-right:auto id=ver-1b_comparison_dist caption=Comparison of concentration as function of distance from the source - at t\=25sec calculated through TMAP8 and analytically + at $t = 25$ s calculated through TMAP8 and analytically. Finally, the diffusive flux ($J$) was compared with the analytic solution where the -flux is proportional to the derivative of the concentration with respect to x and -is given by: - +flux is proportional to the derivative of the concentration with respect to $x$ and +is given by \begin{equation} \label{eq:flux} -J = C_o \; \sqrt{\frac{D}{t\pi}} \; exp \left(\frac{x}{2\sqrt{Dt}}\right) +J = C_0 \; \sqrt{\frac{D}{t\pi}} \; exp \left(\frac{x}{2\sqrt{Dt}}\right). \end{equation} The flux as given by [eq:flux] is compared with values calculated by TMAP8. -The diffusivity, D, and the initial concentration, C$_o$, were both -taken as unity, and the distance, x, was taken as 0.5 in this comparison. -TMAP8 initially under predicts but the results match well subsequently. Comparison +The diffusivity, D, and the initial concentration, $C_0$, were both +taken as unity, and the distance, $x$, was taken as 0.5 in this comparison. +TMAP8 initially underpredicts but the results match well subsequently. Comparison results are shown in [ver-1b_comparison_flux] with a root mean square percentage error of RMSPE = 6.03 %. The error is calculated for $t \geq 10$ s due to infinite value at small $t$. @@ -57,11 +56,11 @@ value at small $t$. style=width:50%;margin-bottom:2%;margin-left:auto;margin-right:auto id=ver-1b_comparison_flux caption=Comparison of flux as function of time at x\=0.5m calculated through - TMAP8 and analytically + TMAP8 and analytically. The oscillations in the permeation graph go away with increasing fineness in the -mesh and in `dt`. +mesh and in the time step `dt`. ## Input files diff --git a/doc/content/verification_and_validation/ver-1c.md b/doc/content/verification_and_validation/ver-1c.md index 8023d701..9edba332 100644 --- a/doc/content/verification_and_validation/ver-1c.md +++ b/doc/content/verification_and_validation/ver-1c.md @@ -3,15 +3,15 @@ # Diffusion Problem with Partially Preloaded Slab This verification problem is taken from [!cite](longhurst1992verification,ambrosek2008verification). Diffusion of tritium through a semi-infinite SiC layer is modeled with an initial -loading of 1 atom/m$^3$ in the first 10 m of a 100 m slab. TMAP4 uses a slab length of 2275 m (the slab length is not specified in the TMAP7 document [!cite](ambrosek2008verification)), however using a smaller slab length was found to not change the results. Additionally, the smaller domain size allows getting a finer simulation mesh for the same computational cost, which improves the agreement between the TMAP8 and analytical calculations. +loading of 1 atom/m$^3$ in the first 10 m of a 100 m slab. TMAP4 uses a slab length of 2275 m (the slab length is not specified in the TMAP7 document [!cite](ambrosek2008verification)); however, using a smaller slab length was found not to change the results. Additionally, the smaller domain size allows getting a finer simulation mesh for the same computational cost, which improves the agreement between the TMAP8 and analytical calculations. Diffusivity is set to 1 m$^2$/s and no trapping is included. The boundary condition on the left-hand side of the slab (at $x=0$ m) is different for the TMAP4 and TMAP7 cases. For TMAP4, an insulating boundary condition is assumed, and the analytical solution is given by [!citep](Carslaw1959conduction): \begin{equation} \label{eq:c_func_4} -C = \frac{C_0}{2} \left[ \text{erf}\bigg(\frac{h-x}{2\sqrt{Dt}}\bigg) + -\text{erf}\bigg(\frac{h+x}{2\sqrt{Dt}}\bigg) \right] +C = \frac{C_0}{2} \left[ \text{erf}\left(\frac{h-x}{2\sqrt{Dt}}\right) + +\text{erf}\left(\frac{h+x}{2\sqrt{Dt}}\right) \right] \end{equation} but for TMAP7, which specifies $C(x=0 \mathrm{\:m})=0$, the analytical solution is [!citep](Carslaw1959conduction) @@ -21,25 +21,24 @@ but for TMAP7, which specifies $C(x=0 \mathrm{\:m})=0$, the analytical solution C = \frac{C_0}{2}\left[2\mathrm{erf}\left(\frac{x}{2\sqrt{Dt}}\right) - \mathrm{erf}\left(\frac{x-h}{2\sqrt{Dt}}\right) - \mathrm{erf}\left(\frac{x+h}{2\sqrt{Dt}}\right)\right] \end{equation} -where $h=10$ m is the thickness of the pre-loaded portion of the layer. +where $h=10$ m is the thickness of the pre-loaded portion of the layer, $C_0$ is the initial concentration in the pre-loaded section of the slab, $\mathrm{erf}$ is the error function, $x$ is the position along the slab, and $D$ is the diffusivity. !alert warning title=Typo in [!cite](longhurst1992verification) The value of $C$ found in [!cite](longhurst1992verification) has a typographical error, $\sqrt{Dt}$ should be at the denominator. [eq:c_func_4] follows the form of [!cite](Carslaw1959conduction). - TMAP4 and TMAP7 verification cases are also evaluated at slightly different locations: TMAP4 verifies the mobile species concentration at three points: -1. a point at the free surface (x = 0 m) -2. a point at the end of the pre-loaded region (x = 10 m) -3. a point in the initially unloaded region (x = 12 m) +1. a point at the free surface ($x = 0$ m) +2. a point at the end of the pre-loaded region ($x = 10$ m) +3. a point in the initially unloaded region ($x = 12$ m) - The comparison of the values calculated with TMAP8 and analytically for the TMAP4 cases is shown in +The comparison of the values calculated with TMAP8 and analytically for the TMAP4 cases is shown in [ver-1c_comparison_time_TMAP4]. TMAP7 verifies the mobile species concentration at similar points -1. a point near the free surface (x = 0.25 m), as the concentration at the surface is specified -2. a point at the end of the pre-loaded region (x = 10 m) -3. a point in the initially unloaded region (x = 12 m) +1. a point near the free surface ($x = 0.25$ m), as the concentration at the surface is specified +2. a point at the end of the pre-loaded region ($x = 10$ m) +3. a point in the initially unloaded region ($x = 12$ m) @@ -50,20 +49,20 @@ The comparison of the values calculated with TMAP8 and analytically for the TMAP image_name=ver-1c_comparison_time_TMAP4.png style=width:50%;margin-bottom:2%;margin-left:auto;margin-right:auto id=ver-1c_comparison_time_TMAP4 - caption=Comparison of concentration as a function of time at x\=0 m, 10 m, and 12 m + caption=Comparison of concentration as a function of time at $x = 0$ m, 10 m, and 12 m calculated with TMAP8 and analytically (TMAP4 cases) !media comparison_ver-1c.py image_name=ver-1c_comparison_time_TMAP7.png style=width:50%;margin-bottom:2%;margin-left:auto;margin-right:auto id=ver-1c_comparison_time_TMAP7 - caption=Comparison of concentration as a function of time at x\=0.25 m, 10 m, and 12 m + caption=Comparison of concentration as a function of time at $x = 0.25$ m, 10 m, and 12 m calculated with TMAP8 and analytically (TMAP7 cases) ## Input files !style halign=left -The input file for this case can be found at [/ver-1c.i], which is also used as test in TMAP8 at [/ver-1c/tests]. The TMAP4 and TMAP7 verification tests use the same input file, +The input file for this case can be found at [/ver-1c.i], which is also used as test in TMAP8 at [/ver-1c/tests]. The TMAP4 and TMAP7 verification tests use the same input file, but different command line arguments for TMAP4. !listing /test/tests/ver-1c/tests line=NeumannBC diff --git a/doc/content/verification_and_validation/ver-1d.md b/doc/content/verification_and_validation/ver-1d.md index 317e5990..9d08a06b 100644 --- a/doc/content/verification_and_validation/ver-1d.md +++ b/doc/content/verification_and_validation/ver-1d.md @@ -8,7 +8,7 @@ This verification problem is taken from [!cite](longhurst1992verification). It m \begin{equation} \label{eqn:diffusion_mobile} - \frac{dC_M}{dt} = - \nabla D \nabla C_M - \text{trap\_per\_free} \cdot \frac{dC_T}{dt}, + \frac{dC_M}{dt} = \nabla D \nabla C_M - \text{trap\_per\_free} \cdot \frac{dC_T}{dt}, \end{equation} \begin{equation} \label{eqn:trapped_rate} @@ -24,7 +24,7 @@ The breakthrough time may have one of two limiting values depending on whether t \begin{equation} \label{eqn:zeta} - \zeta = \frac{\lambda^2 \nu}{\rho D_o} exp \left( \frac{E_d - \epsilon}{kT} \right) + \frac{c}{\rho} + \zeta = \frac{\lambda^2 \nu}{\rho D_0} \exp \left( \frac{E_d - \epsilon}{kT} \right) + \frac{c}{\rho} \end{equation} where @@ -35,7 +35,7 @@ $\nu$ = Debye frequency ($\approx$ $10^{13} \; s^{-1}$) $\rho$ = trapping site fraction -$D_o$ = diffusivity pre-exponential +$D_0$ = diffusivity pre-exponential $E_d$ = diffusion activation energy @@ -47,13 +47,15 @@ $T$ = temperature $c$ = dissolved gas atom fraction -The discriminant for which regime is dominant is the ratio of $\zeta$ to c/$\rho$. If $\zeta$ $\gg$ c/$\rho$ then the effective diffusivity regime applies, and the permeation transient is identical to the standard diffusion transient but with the diffusivity replaced by an effective diffusivity. +The discriminant for which regime is dominant is the ratio of $\zeta$ to c/$\rho$. If $\zeta$ $\gg$ c/$\rho$, then the effective diffusivity regime applies, and the permeation transient is identical to the standard diffusion transient, with the diffusivity replaced by an effective diffusivity \begin{equation} \label{eqn:Deff} D_{eff} = \frac{D}{1 + \frac{1}{\zeta}} \end{equation} +to account for the fact that trapping leads to slower transport. + In this limit, the breakthrough time, defined as the intersection of the steepest tangent to the diffusion transient with the time axis, will be \begin{equation} @@ -63,30 +65,28 @@ In this limit, the breakthrough time, defined as the intersection of the steepes where $l$ is the thickness of the slab and D is the diffusivity of the gas through the material. The permeation transient is then given by - \begin{equation} \label{eqn:Jp} - J_p = \frac{c_o D}{l} \Bigg\{ 1 + 2 \sum_{m=1}^{\infty} \left[ (-1)^m \exp \left( -m^2 \frac{t}{2 \; \tau_{b_e}} \right) \right] \Bigg\} + J_p = \frac{C_0 D}{l} \left\{ 1 + 2 \sum_{m=1}^{\infty} \left[ (-1)^m \exp \left( -m^2 \frac{t}{2 \; \tau_{b_e}} \right) \right] \right\} \end{equation} - [!cite](longhurst2005verification) where $\tau_{b_e}$ is defined in [eqn:tau_be] -In the deep-trapping limit, $\zeta$ $\approx$ c/$\rho$, and no permeation occurs until essentially all the traps have been filled. Then permeation rapidly turns on to its state value. The breakthrough time is given by +In the deep-trapping limit, $\zeta$ $\approx$ c/$\rho$, and no permeation occurs until essentially all the traps have been filled. Then the system quickly reaches steady-state. The breakthrough time is given by \begin{equation} \label{eqn:tau_bd} - \tau_{b_d} = \frac{l^2 \rho}{2 \; c_o \; D} + \tau_{b_d} = \frac{l^2 \rho}{2 \; C_0 \; D} \end{equation} -where $c_o$ is the steady dissolved gas concentration at the upstream (x = 0) side. +where $C_0$ is the steady dissolved gas concentration at the upstream (x = 0) side. Using TMAP8 we examine these two different regimes, one where diffusion is the rate-limiting step, and one where trapping is the rate-limiting step. The upstream-side starting concentration of 0.0001 atom fraction, a diffusivity of 1 $m^2$/s, a trapping site fraction of 0.1, $\lambda^2 = 10^{-15} \; m^2$, and a temperature of 1000 K is considered. ## Diffusion-limited -For the effective diffusivity limit, we selected $\epsilon/k = 100 K$ to give $\zeta = 91.47 c/\rho$. The comparison results are presented in [ver-1d_comparison_diffusion] with a root mean square percentage error of RMSPE = 0.96 % for $t \geq 0.4$ s.. +For the effective diffusivity limit, we selected $\epsilon/k = 100$ K to give $\zeta = 91.47 c/\rho$. The comparison results are presented in [ver-1d_comparison_diffusion] with a root mean square percentage error of RMSPE = 0.96% for $t \geq 0.4$ s. !media comparison_ver-1d.py image_name=ver-1d_comparison_diffusion.png @@ -96,7 +96,7 @@ For the effective diffusivity limit, we selected $\epsilon/k = 100 K$ to give $\ ## Trapping-limited -For the deep trapping limit we took $\epsilon/k = 10000 K$ to give $\zeta = 1.00454 c/\rho$. The comparison results are presented in [ver-1d_comparison_trapping]. +For the deep trapping limit we took $\epsilon/k = 10000$ K to give $\zeta = 1.00454 c/\rho$. The comparison results are presented in [ver-1d_comparison_trapping]. !media comparison_ver-1d.py image_name=ver-1d_comparison_trapping.png @@ -104,14 +104,11 @@ For the deep trapping limit we took $\epsilon/k = 10000 K$ to give $\zeta = 1.00 id=ver-1d_comparison_trapping caption=Permeation transient in a slab subject to strong trapping. - - - ### Notes The trapping test input file can generate oscillations in the solution due to the feedback loop between the diffusion PDE and trap evolution ODE. In order for the oscillations to not take over the simulation, it seems that the ratio of the inverse of the Fourier number must be kept -sufficiently high, e.g. `h^2 / (D * dt)`. Included in this directory are three +sufficiently high, e.g. $h^2 / (D * dt) \gg 1$. Included in this directory are three `png` files that show the permeation for different `h` and `dt` values. They are summarized below: @@ -127,7 +124,9 @@ To keep the oscillations damped, the verification is run with nx=1000 and an ada \begin{equation} C_M(x=0) = \tanh(3t). \end{equation} -This takes the BC to 99.5 % of it's actual value in 3 s, which is a small fraction of the breakthrough time of 500 s. + +This takes the BC to 99.5% of its actual value in 3 s, which is a small fraction of the breakthrough time of 500 s. +It therefore reduces numerical oscillations without affecting the expected solution. ## Input files diff --git a/doc/content/verification_and_validation/ver-1dc.md b/doc/content/verification_and_validation/ver-1dc.md index 7acdbfc1..979362e4 100644 --- a/doc/content/verification_and_validation/ver-1dc.md +++ b/doc/content/verification_and_validation/ver-1dc.md @@ -10,7 +10,7 @@ This problem models permeation through a membrane with a constant source in whic \begin{equation} \label{eqn:diffusion_mobile} - \frac{dC_M}{dt} = - \nabla D \nabla C_M - \text{trap\_per\_free} \cdot \sum_{i=1}^{3} \frac{dC_{T_i}}{dt} , + \frac{dC_M}{dt} = \nabla D \nabla C_M - \text{trap\_per\_free} \cdot \sum_{i=1}^{3} \frac{dC_{T_i}}{dt} , \end{equation} and, for $i=1$, $i=2$, and $i=3$: \begin{equation} @@ -26,7 +26,7 @@ where $C_M$ is the concentrations of the mobile, $C_{T_i}$ is the trapped specie The trapping parameter is defined by \begin{equation} \label{eqn:zeta} - \zeta = \frac{\lambda^2 \nu}{\rho D_o} exp \left( \frac{E_d - \epsilon}{kT} \right) + \frac{c}{\rho}, + \zeta = \frac{\lambda^2 \nu}{\rho D_0} \exp \left( \frac{E_d - \epsilon}{kT} \right) + \frac{c}{\rho}, \end{equation} where @@ -37,7 +37,7 @@ $\nu$ = Debye frequency ($\approx$ $10^{13} \; s^{-1}$) $\rho$ = trapping site fraction -$D_o$ = diffusivity pre-exponential +$D_0$ = diffusivity pre-exponential $E_d$ = diffusion activation energy @@ -60,9 +60,9 @@ Three traps that are relatively weak are assumed to be active in a slab. The tra \begin{equation} \label{eqn:Jp} - J_p = \frac{c_o D}{l} \Bigg\{ 1 + 2 \sum_{m=1}^{\infty} \left[ (-1)^m \exp \left( -m^2 \frac{t}{2 \; \tau_{b_e}} \right) \right] \Bigg\}, + J_p = \frac{C_0 D}{l} \left\{ 1 + 2 \sum_{m=1}^{\infty} \left[ (-1)^m \exp \left( -m^2 \frac{t}{2 \; \tau_{b_e}} \right) \right] \right\}, \end{equation} -where $c_o$ is the steady dissolved gas concentration at the upstream (x = 0) side, $l$ is the thickness of the slab, $D$ is the diffusivity of the gas through the material, and $\tau_{b_e}$, the breakthrough time, is defined as +where $C_0$ is the steady dissolved gas concentration at the upstream (x = 0) side, $l$ is the thickness of the slab, $D$ is the diffusivity of the gas through the material, and $\tau_{b_e}$, the breakthrough time, is defined as \begin{equation} \label{eqn:tau_be} @@ -74,11 +74,11 @@ where $D_{eff}$, the effective diffusivity, is defined as \label{eqn:Deff} D_{eff} = \frac{D}{1 + \sum_{i=1}^3 1 / \zeta_i}, \end{equation} -where $\zeta_i$ is the trapping parameter of trap $i$. The trapping parameters, $\zeta_i$, calculated from [eqn:zeta] for the three traps are 91.47930 $c/\rho$, 61.65009 $c/\rho$, 45.93069 $c/\rho$. +where $\zeta_i$ is the trapping parameter of trap $i$. +The trapping parameters, $\zeta_i$, calculated from [eqn:zeta] for the three traps are 91.47930 $c/\rho$, 61.65009 $c/\rho$, 45.93069 $c/\rho$. !alert warning title=Typo in [!cite](ambrosek2008verification) -The $\zeta_i$ values of the three traps from [!cite](ambrosek2008verification) have a typographical error: They are three orders of magnitude lower than the correct values. however, it dos not impact the final analytical solution. - +The $\zeta_i$ values of the three traps from [!cite](ambrosek2008verification) have a typographical error: They are three orders of magnitude lower than the correct values. However, it does not impact the final analytical solution. ## Results and comparison against analytical solution diff --git a/doc/content/verification_and_validation/ver-1dd.md b/doc/content/verification_and_validation/ver-1dd.md index 2fac53d8..9eddabd5 100644 --- a/doc/content/verification_and_validation/ver-1dd.md +++ b/doc/content/verification_and_validation/ver-1dd.md @@ -6,11 +6,11 @@ This verification problem is taken from [!cite](ambrosek2008verification) and builds on the capabilities verified in [ver-1d](ver-1d.md) and [ver-1dc](ver-1dc.md). The configuration and modeling parameters are the same as in [ver-1d](ver-1d.md), except that, in the current case, there are no traps in the membrane. This case is simulated in [/ver-1dd.i]. -This problem models permeation through a membrane with a constant source in which no trap presented. We solve +This problem models permeation through a membrane with a constant source in which no trap presented. We solve \begin{equation} \label{eqn:diffusion_mobile} - \frac{dC_M}{dt} = - \nabla D \nabla C_M , + \frac{dC_M}{dt} = \nabla D \nabla C_M , \end{equation} where $C_M$ is the concentrations of the mobile, $D$ is the diffusivity of the mobile species, and $t$ is the time. @@ -20,20 +20,20 @@ where $C_M$ is the concentrations of the mobile, $D$ is the diffusivity of the m \begin{equation} \label{eqn:Jp} - J_p = \frac{c_o D}{l} \Bigg\{ 1 + 2 \sum_{m=1}^{\infty} \left[ (-1)^m \exp \left( -m^2 \frac{t}{2 \; \tau_{b_e}} \right) \right] \Bigg\}, + J_p = \frac{C_0 D}{l} \left\{ 1 + 2 \sum_{m=1}^{\infty} \left[ (-1)^m \exp \left( -m^2 \frac{t}{2 \; \tau_{b_e}} \right) \right] \right\}, \end{equation} -where $c_o$ is the steady dissolved gas concentration at the upstream (x = 0) side, $l$ is the thickness of the slab, $D$ is the diffusivity of the gas through the material, and $\tau_{b_e}$, the breakthrough time, is defined as +where $C_0$ is the steady dissolved gas concentration at the upstream (x = 0) side, $l$ is the thickness of the slab, $D$ is the diffusivity of the gas through the material, and $\tau_{b_e}$, the breakthrough time, is defined as \begin{equation} \label{eqn:tau_be} \tau_{b_e} = \frac{l^2}{2 \; \pi^2 \; D_{eff}}, \end{equation} -where $D_{eff}$, the effective diffusivity, is same with $D$ due to no traps in the membrane. +where $D_{eff}$, the effective diffusivity, is equal to $D$ due to the absence of traps in the membrane. ## Results and comparison against analytical solution -The analytical solution for the permeation transient is compared with TMAP8 results in [ver-1dd_comparison_diffusion]. The graphs for the theoretical flux and the calculated flux are in good agreement, with root mean square percentage errors (RMSPE) of RMSPE = 0.14 % for $t \geq 0.01$ s. The breakthrough time calculated from the analytical solution provided by [eqn:tau_be] is 0.05 s, and the breakthrough time from TMAP8 is 0.05 s, so it matches. +The analytical solution for the permeation transient is compared with TMAP8 results in [ver-1dd_comparison_diffusion]. The graphs for the theoretical flux and the calculated flux are in good agreement, with root mean square percentage errors (RMSPE) of RMSPE = 0.14% for $t \geq 0.01$ s. The breakthrough time calculated from the analytical solution provided by [eqn:tau_be] is 0.05 s, and the breakthrough time from TMAP8 is 0.05 s, which matches. !media comparison_ver-1dd.py image_name=ver-1dd_comparison_diffusion.png diff --git a/doc/content/verification_and_validation/ver-1e.md b/doc/content/verification_and_validation/ver-1e.md index 2ec3e74b..d87c830a 100644 --- a/doc/content/verification_and_validation/ver-1e.md +++ b/doc/content/verification_and_validation/ver-1e.md @@ -4,39 +4,39 @@ ## Test Description -This verification problem is taken from [!cite](longhurst1992verification, ambrosek2008verification). In this problem, a composite structure of PyC and SiC is modeled with a constant concentration boundary condition of the free surface of PyC and zero concentration boundary condition on the free surface of the SiC. The steady state solution for the PyC is given as: +This verification problem is taken from [!cite](longhurst1992verification, ambrosek2008verification). In this problem, a composite structure of PyC and SiC is modeled with a constant concentration boundary condition of the free surface of PyC and zero concentration boundary condition on the free surface of the SiC. The steady-state solution for the PyC is given as: \begin{equation} \label{eqn:steady_state_pyc} - C = C_o \left[1 + \frac{x}{l} \left(\frac{a D_{PyC}}{a D_{PyC} + l D_{SiC}} - 1 \right) \right] + C = C_0 \left[1 + \frac{x}{l} \left(\frac{a D_{PyC}}{a D_{PyC} + l D_{SiC}} - 1 \right) \right] \end{equation} while the concentration profile for the SiC layer is given as: \begin{equation} \label{eqn:steady_state_sic} - C = C_o \left(\frac{a+l-x}{l} \right) \left(\frac{a D_{PyC}}{a D_{PyC} + l D_{SiC}} \right) + C = C_0 \left(\frac{a+l-x}{l} \right) \left(\frac{a D_{PyC}}{a D_{PyC} + l D_{SiC}} \right) \end{equation} where $x$ = distance from free surface of PyC - $a$ = thickness of the PyC layer (33 $\mu m$) + $a$ = thickness of the PyC layer (33 $\mu$m) - $l$ = thickness of the SiC layer (63 $\mu m$ in TMAP4 verification case, and 66 $\mu m$ in TMAP7 verification case) + $l$ = thickness of the SiC layer (63 $\mu$m in TMAP4 verification case, and 66 $\mu$m in TMAP7 verification case) - $C_o$ = concentration at the PyC free surface (50.7079 moles/m$^3$) + $C_0$ = concentration at the PyC free surface (50.7079 moles/m$^3$) - $D_{PyC}$ = diffusivity in PyC (1.274 x 10$^{-7}$ m$^2$/sec) + $D_{PyC}$ = diffusivity in PyC (1.274 $\times$ 10$^{-7}$ m$^2$/s) - $D_{SiC}$ = diffusivity in SiC (2.622 x 10$^{-11}$ m$^2$/sec) + $D_{SiC}$ = diffusivity in SiC (2.622 $\times$ 10$^{-11}$ m$^2$/s) The analytical transient solution for the concentration in the SiC side of the composite slab is given as: \begin{equation} \label{eqn:transient} - C = C_o \Bigg\{ \frac{D_{PyC}(l-x)}{l D_{PyC} + a D_{SiC}} - 2 \sum_{n=1}^{\infty} \frac{\sin(a \lambda_n) \sin(k l \lambda_n) \sin \left[k (l-x) \lambda_n \right]}{\lambda_n \left[ a \sin^2(k l \lambda_n) + l \sin^2 (a \lambda_n) \right]} \exp(-D_{PyC} \lambda_n^2 t) \Bigg\} + C = C_0 \left\{ \frac{D_{PyC}(l-x)}{l D_{PyC} + a D_{SiC}} - 2 \sum_{n=1}^{\infty} \frac{\sin(a \lambda_n) \sin(k l \lambda_n) \sin \left[k (l-x) \lambda_n \right]}{\lambda_n \left[ a \sin^2(k l \lambda_n) + l \sin^2 (a \lambda_n) \right]} \exp(-D_{PyC} \lambda_n^2 t) \right\} \end{equation} where @@ -47,15 +47,15 @@ and $\lambda_n$ are the roots of \begin{equation} \label{eqn:roots} - \frac{1}{tan(\lambda a)} + \frac{1}{k \; tan(k l \lambda)} = 0. + \frac{1}{\tan(\lambda a)} + \frac{1}{k \; \tan(k l \lambda)} = 0. \end{equation} -!alert warning title=Typo in [!cite](ambrosek2008verification) -[eqn:roots] for the roots of $\lambda$ is obtained from [!cite](longhurst1992verification). The equation presented in [!cite](ambrosek2008verification) has a typographical error and gives incorrect results. +!alert warning title=Typo in [!cite](longhurst1992verification,ambrosek2008verification) +[eqn:roots] for the roots of $\lambda$ is different from the equations provided in [!cite](longhurst1992verification,ambrosek2008verification), as they both have typographical errors and give different analytical solutions from the one provided in the figures of these reports. ## Results -[ver-1e_comparison_dist] shows the comparison of the TMAP8 calculation and the analytical solution for concentration after steady-state is reached. The plots show the TMAP8 and analytical solution comparisons for both the TMAP4 case ($l$ = 63 $\mu m$) and the TMAP7 case ($l$ = 66 $\mu m$). +[ver-1e_comparison_dist] shows the comparison of the TMAP8 calculation and the analytical solution for concentration after steady-state is reached. The plots show the TMAP8 and analytical solution comparisons for both the TMAP4 case ($l = 63$ $\mu$m) and the TMAP7 case ($l = 66$ $\mu$m). !media comparison_ver-1e.py image_name=ver-1e_comparison_dist.png @@ -63,7 +63,7 @@ and $\lambda_n$ are the roots of id=ver-1e_comparison_dist caption=Comparison of TMAP8 calculation with the analytical solution. Bold text next to the plot curves shows the root mean square percentage error (RMSPE) between the TMAP8 prediction and analytical solution for the TMAP4 and TMAP7 verification cases. -For transient solution comparison, the concentation at a point, which is $x$ $\mu m$ away from the PyC-SiC interface into the SiC layer, is obtained using the TMAP code as well as analytically. [ver-1e_comparison_time] shows comparison of the TMAP calculation with the analytical solution for this transient case. In the TMAP4 case, $x$ = 8 $\mu m$, and in the TMAP7 case $x$ = 15.75 $\mu m$. There is good agreement between TMAP and the analytical solution for both steady state as well as transient cases. In both cases, the root mean square percentage error (RMSPE) is under 0.2 %. +For transient solution comparison, the concentration at a point, which is $x$ $\mu$m away from the PyC-SiC interface into the SiC layer, is obtained using the TMAP code and analytically. [ver-1e_comparison_time] shows comparison of the TMAP calculation with the analytical solution for this transient case. In the TMAP4 case, $x$ = 8 $\mu$m, and in the TMAP7 case $x$ = 15.75 $\mu$m. There is good agreement between TMAP and the analytical solution for both steady state and transient cases. In both cases, the root mean square percentage error (RMSPE) is under 0.2%. !media comparison_ver-1e.py image_name=ver-1e_comparison_time.png @@ -77,7 +77,7 @@ The error is calculated between the TMAP8 and analytical solution values after $ image_name=ver-1e_comparison_time_closeup.png style=width:50%;margin-bottom:2%;margin-left:auto;margin-right:auto id=ver-1e_comparison_time_zoomed - caption=Zoomed-in view of comparison of TMAP8 calculation with the analytical solution for $t$ < 1 s. The analytical solution shows unphysical predictions close to $t$ = 0 s. + caption=Zoomed-in view of comparison of TMAP8 calculation with the analytical solution for $t$ < 1 s. The analytical solution shows unphysical predictions close to $t$ = 0 s, whereas TMAP8's predictions are reasonable. ## Input files diff --git a/doc/content/verification_and_validation/ver-1fa.md b/doc/content/verification_and_validation/ver-1fa.md index 21db27bc..5a991fc5 100644 --- a/doc/content/verification_and_validation/ver-1fa.md +++ b/doc/content/verification_and_validation/ver-1fa.md @@ -14,7 +14,7 @@ This case models a heat conduction problem through a slab with a heat source. Th where $T$ is the temperature, $\rho$ is the density, $C_P$ is the specific heat, $k$ is the thermal conductivity, and $Q$ is the internal volumetric heat generation rate. -One end of the slab is kept at a constant temperature of 300K while the other end acts as an adiabatic surface. +One end of the slab is kept at a constant temperature of 300 K while the other end acts as an adiabatic surface. This case uses internal volumetric heat generation rate of $Q = 1 \times 10^{4}$ W/m$^3$ in the slab. The thickness of slab, $L$, is 1.6 m, and the thermal conductivity is $k = 10$ W/m/K. The surface temperature, $T_s$, on the end with a constant temperature, is 300 K. The material density and specific heat are assumed to be $\rho = 1$ kg/m$^3$ and $C_P = 1$ J/kg/K, respectively. diff --git a/doc/content/verification_and_validation/ver-1fb.md b/doc/content/verification_and_validation/ver-1fb.md index 641a64f7..70146bda 100644 --- a/doc/content/verification_and_validation/ver-1fb.md +++ b/doc/content/verification_and_validation/ver-1fb.md @@ -16,7 +16,7 @@ where $T$ is the temperature, $\rho$ is the density, $C_P$ is the specific heat, The ends of a slab are kept fixed at different temperatures. The temperature distribution in the slab evolves from an initial state to steady-state. -In this case, the thickness of slab, $L$, is 4.0 m, the thermal conductivity is $k =1$ W/m/K, the material density is assumed to be $\rho = 1$ kg/m$^3$, and the specific heat is assumed to be $C_P = 1$ J/kg/K. The fixed surface temperature, $T_0$ and $T_1$, on both ends are defined as 400 K and 300 K, respectively. +In this case, the thickness of slab, $L$, is 4.0 m, the thermal conductivity is $k =10$ W/m/K, the material density is assumed to be $\rho = 1$ kg/m$^3$, and the specific heat is assumed to be $C_P = 10$ J/kg/K. The fixed surface temperature, $T_0$ and $T_1$, on both ends are defined as 400 K and 300 K, respectively. ## Analytical solution @@ -32,8 +32,10 @@ where $x$ is the distance across the slab, $t$ is the time, $\lambda_m$ is a coe \alpha = \frac{k}{\rho C_p}. \end{equation} -!alert warning title=Typo in analytical solution from [!cite](longhurst1992verification) -Both TMAP4 ([!citep](longhurst1992verification)) and TMAP7 ([!citep](ambrosek2008verification)) provide analytical solutions, but they use different equations. At the initial time and steady state, the solution from TMAP7 matches the real thermal distribution, whereas the solution from TMAP4 does not. Thus, TMAP8 select the solution from TMAP7 as the analytical solution. +!alert! warning title=Typo in analytical solution from [!cite](longhurst1992verification) +Both TMAP4 [!citep](longhurst1992verification) and TMAP7 [!citep](ambrosek2008verification) provide analytical solutions, but they use different equations. +At the initial time and steady state, the solution from TMAP7 matches the real thermal distribution, whereas the solution from TMAP4 does not. Thus, TMAP8 selects the solution from TMAP7 as the analytical solution. +!alert-end! ## Results diff --git a/doc/content/verification_and_validation/ver-1g.md b/doc/content/verification_and_validation/ver-1g.md index 3c5f3b09..4e64ac02 100644 --- a/doc/content/verification_and_validation/ver-1g.md +++ b/doc/content/verification_and_validation/ver-1g.md @@ -2,13 +2,13 @@ # A Simple Forward Chemical Reaction -This verification problem is taken from [!cite](longhurst1992verification,ambrosek2008verification). A simple time-dependent chemical reaction given by: +This verification problem is taken from [!cite](longhurst1992verification,ambrosek2008verification). A simple time-dependent chemical reaction given by \begin{equation} A + B \Rightarrow AB \end{equation} -is modeled in a functional enclosure. The reaction rate, R, is positive if the species AB is being produced in the reaction and negative if it is being consumed. The forward rate coefficient, K, for the reaction has no spatial or time dependence. The reaction rate is: +is modeled in a functional enclosure. The reaction rate, $R$, is positive if the species AB is being produced in the reaction and negative if it is being consumed. The forward rate coefficient, $K$, for the reaction has no spatial or time dependence. The reaction rate is: \begin{equation} R = K C_A C_B @@ -23,24 +23,24 @@ R = -\frac{d[C_A]}{dt} = -\frac{d[C_B]}{dt} = \frac{d[C_{AB}]}{dt} = K C_A C_B where $C_{AB}$ is the concentration of species AB. The analytical solution for the concentration of species AB as a function of time ($t$) is given as: \begin{equation} -C_{AB} = C_{B_o} \frac{1 - exp{[K t (C_{B_o} - C_{A_o})]}}{1 - \frac{C_{B_o}}{C_{A_o}}exp{[K t (C_{B_o} - C_{A_o})]}} +C_{AB} = C_{B_0} \frac{1 - exp{[K t (C_{B_0} - C_{A_0})]}}{1 - \frac{C_{B_0}}{C_{A_0}}exp{[K t (C_{B_0} - C_{A_0})]}} \label{eq:conc_AB} \end{equation} -where $C_{A_o}$ and $C_{B_o}$ are the initial concentrations of A and B, respectively. For the special case when $C_{A_o}$ and $C_{B_o}$ are equal, [eq:conc_AB] becomes: +where $C_{A_0}$ and $C_{B_0}$ are the initial concentrations of A and B, respectively. For the special case when $C_{A_0}$ and $C_{B_0}$ are equal, [eq:conc_AB] becomes: \begin{equation} -C_{AB} = C_{A_o} - \frac{1}{\frac{1}{C_{A_o}} + Kt} +C_{AB} = C_{A_0} - \frac{1}{\frac{1}{C_{A_0}} + Kt} \end{equation} For this verification exercise, three cases were considered: (a) the initial concentrations of A and B are equal, (b) the initial concentrations of A and B are different and use the TMAP4 verification case values, and (c) the initial concentrations of A and B are different and use the TMAP7 verification case values. The equal concentration case is the same in TMAP4 and TMAP7, and is used for verification of TMAP8 here. While both TMAP4 and TMAP7 verification cases ((b) and (c) respectively) have the same analytical solution, the different initial conditions produce different $C_{AB}$ values over time, and we replicate both results here. -For case (a), the initial pressures of A and B were 1 $\mu$Pa, and the reaction rate (K) is 4.14 x 10$^3$ $\mu$m$^3$ / atom$\cdot$s. For case (b) the initial pressure of A was same as in case (a) while the initial pressure of B is 0.1 $\mu$Pa as per TMAP4. For case (c) the initial pressure of A was same as in case (a) while the initial pressure of B is 0.5 $\mu$Pa. In all cases, the initial pressures of A and B are first converted to their initial concentrations $C_{A_o}$ and $C_{B_o}$ using the ideal gas law to be used in the TMAP8 simulations and analytical solutions. The initial concentration $C_{i_0}$ of component $i$ is +For case (a), the initial pressures of A and B were 1 $\mu$Pa, and the reaction rate $K$ is 4.14 $\times$ 10$^3$ $\mu$m$^3$ / atom$\cdot$s. For case (b) the initial pressure of A was same as in case (a) while the initial pressure of B is 0.1 $\mu$Pa as per TMAP4. For case (c) the initial pressure of A was same as in case (a) while the initial pressure of B is 0.5 $\mu$Pa. In all cases, the initial pressures of A and B are first converted to their initial concentrations $C_{A_0}$ and $C_{B_0}$ using the ideal gas law to be used in the TMAP8 simulations and analytical solutions. The initial concentration $C_{i_0}$ of component $i$ is \begin{equation} C_{i_0} = 10^{-18} \frac{P_{i_0} N_a}{R T}, \end{equation} -where $P_{i_0}$ is the initial pressure, $N_a$ is Avogardro's constant, $R$ is the gas constant (from [PhysicalConstants](source/utils/PhysicalConstants.md)), and $T$ = 298.15 K (25$\deg$ C) is the temperature. The factor $10^{-18}$ converts the concentration from atoms/m$^3$ to atoms/$\mu$m$^3$. +where $P_{i_0}$ is the initial pressure, $N_a$ is Avogardro's constant, $R$ is the gas constant (from [PhysicalConstants](source/utils/PhysicalConstants.md)), and $T$ = 298.15 K (25$\deg$C) is the temperature. The factor $10^{-18}$ converts the concentration from atoms/m$^3$ to atoms/$\mu$m$^3$. A comparison of the concentration of AB as a function of time is plotted in [ver-1g_comparison_equal_conc] for case (a), and [ver-1g_comparison_diff_conc] for the cases (b) and (c), respectively. The TMAP8 calculations are found to be in good agreement with the analytical solution with the root mean square percentage errors of (a) RMSPE = 0.27 %, (b) RMSPE = 0.22 %, and (c) RMSPE = 0.24 %. diff --git a/test/tests/val-2a/comparison_val-2a.py b/test/tests/val-2a/comparison_val-2a.py index db1fa517..665bb28c 100644 --- a/test/tests/val-2a/comparison_val-2a.py +++ b/test/tests/val-2a/comparison_val-2a.py @@ -54,9 +54,9 @@ def numerical_solution_on_experiment_input(experiment_input, tmap_input, tmap_ou ax = fig.add_subplot(gs[0]) ax.plot(simulation_time_TMAP4/3600, simulation_recom_flux_right_TMAP4, linestyle='-', label=r"TMAP8", c='tab:gray') -ax.plot(experiment_time_TMAP4/3600, experiment_flux_TMAP4, linestyle='--', label=r"experiment", c='k') +ax.plot(experiment_time_TMAP4/3600, experiment_flux_TMAP4, linestyle='--', label=r"Experiment", c='k') -ax.set_xlabel(u'time (hr)') +ax.set_xlabel(u'Time (hr)') ax.set_ylabel(u"Deuterium flux (atom/m$^2$/s)") ax.legend(loc="best") ax.set_ylim(bottom=0) @@ -91,7 +91,7 @@ def numerical_solution_on_experiment_input(experiment_input, tmap_input, tmap_ou ax.plot(coordinate_x, normal_distribution, linestyle='-', label=r"normal distribution", c='k') ax.plot(coordinate_x, piecewise, linestyle='-', label=r"piecewise", c='gray') -ax.set_xlabel(u'depth (m)') +ax.set_xlabel(u'Depth (m)') ax.set_ylabel(u"Source rate (atom/m$^3$/s)") ax.legend(loc="best") ax.set_ylim(bottom=0) diff --git a/test/tests/val-2c/comparison_val-2c.py b/test/tests/val-2c/comparison_val-2c.py index 6eb63127..240dea74 100644 --- a/test/tests/val-2c/comparison_val-2c.py +++ b/test/tests/val-2c/comparison_val-2c.py @@ -163,11 +163,11 @@ def molecules_microns3_to_Ci_m3(x): time_sample, experimental_hto_enclosure_concentration_sampled, tmap8_hto_enclosure_concentration_Ci_sampled = common_times_generator(experimental_time_hto, experimental_hto_enclosure_concentration, tmap8_time,tmap8_hto_enclosure_concentration_Ci) RMSE = np.sqrt(np.mean((tmap8_hto_enclosure_concentration_Ci_sampled-experimental_hto_enclosure_concentration_sampled)**2)) RMSPE = RMSE*100/np.mean(experimental_hto_enclosure_concentration_sampled) -ax.text(0.8,3.9e-5, 'RMSPE = %.2f '%RMSPE+'%',fontweight='bold') +ax.text(0.8,3.9e-5, '(immediate) RMSPE = %.2f '%RMSPE+'%',fontweight='bold') time_sample_delay, experimental_hto_enclosure_concentration_delay_sampled, tmap8_hto_enclosure_concentration_Ci_delay_sampled = common_times_generator(experimental_time_hto, experimental_hto_enclosure_concentration, tmap8_time_delay,tmap8_hto_enclosure_concentration_Ci_delay) RMSE = np.sqrt(np.mean((tmap8_hto_enclosure_concentration_Ci_delay_sampled-experimental_hto_enclosure_concentration_delay_sampled)**2)) RMSPE = RMSE*100/np.mean(experimental_hto_enclosure_concentration_delay_sampled) -ax.text(7.5,1.2e-5, 'RMSPE (delay) = %.2f '%RMSPE+'%',fontweight='bold') +ax.text(7.5,1.2e-5, '(delay) RMSPE = %.2f '%RMSPE+'%',fontweight='bold') # save figure plt.savefig('val-2c_comparison_TMAP8_Exp_HTO_Ci.png', bbox_inches='tight', dpi=300); @@ -200,11 +200,11 @@ def molecules_microns3_to_Ci_m3(x): tmap8_t2_enclosure_concentration_Ci_sampled = sample_solution(experimental_time_t2[1:], tmap8_time, tmap8_t2_enclosure_concentration_Ci) RMSE = np.sqrt(np.mean((tmap8_t2_enclosure_concentration_Ci_sampled-experimental_t2_enclosure_concentration[1:])**2)) RMSPE = RMSE*100/np.mean(experimental_t2_enclosure_concentration[1:]) -ax.text(0.5,1.2e-4, 'RMSPE (delay) = %.2f '%RMSPE+'%',fontweight='bold') +ax.text(0.5,1.2e-4, '(immediate) \n RMSPE = %.2f '%RMSPE+'%',fontweight='bold') tmap8_t2_enclosure_concentration_Ci_delay_sampled = sample_solution(experimental_time_t2[1:], tmap8_time_delay, tmap8_t2_enclosure_concentration_Ci_delay) RMSE = np.sqrt(np.mean((tmap8_t2_enclosure_concentration_Ci_delay_sampled-experimental_t2_enclosure_concentration[1:])**2)) RMSPE = RMSE*100/np.mean(experimental_t2_enclosure_concentration[1:]) -ax.text(19,8e-4, 'RMSPE (delay) = %.2f '%RMSPE+'%',fontweight='bold') +ax.text(19,8e-4, '(delay) RMSPE = %.2f '%RMSPE+'%',fontweight='bold') # save figure plt.savefig('val-2c_comparison_TMAP8_Exp_T2_Ci.png', bbox_inches='tight', dpi=300); diff --git a/test/tests/ver-1a/comparison_ver-1a.py b/test/tests/ver-1a/comparison_ver-1a.py index 8100bc39..bc8f2fcf 100644 --- a/test/tests/ver-1a/comparison_ver-1a.py +++ b/test/tests/ver-1a/comparison_ver-1a.py @@ -210,7 +210,7 @@ def analytical_expression_flux(t, P_0, D, S, V, T, A, l): ax = fig.add_subplot(gs[0]) ax.plot(tmap8_time,tmap8_release_fraction_right,label=r"TMAP8",c='tab:gray') ax.plot(time_analytical,analytical_release_fraction_TMAP4,label=r"Analytical (TMAP4)",c='k', linestyle='--') -ax.set_xlabel(u'Time(s)') +ax.set_xlabel(u'Time (s)') ax.set_ylabel(r"Fractional release") ax.legend(loc="best") ax.set_xlim(left=0) @@ -230,7 +230,7 @@ def analytical_expression_flux(t, P_0, D, S, V, T, A, l): ax = fig.add_subplot(gs[0]) ax.plot(tmap8_time,tmap8_release_fraction_left,label=r"TMAP8",c='tab:gray') ax.plot(time_analytical,analytical_release_fraction_TMAP7,label=r"Analytical (TMAP7)",c='k', linestyle='--') -ax.set_xlabel(u'Time(s)') +ax.set_xlabel(u'Time (s)') ax.set_ylabel(r"Fractional release") ax.legend(loc="best") ax.set_xlim(left=0) @@ -250,7 +250,7 @@ def analytical_expression_flux(t, P_0, D, S, V, T, A, l): ax = fig.add_subplot(gs[0]) ax.plot(tmap8_time,tmap8_flux_right,label=r"TMAP8",c='tab:gray') ax.plot(time_analytical,analytical_flux_TMAP7,label=r"Analytical (TMAP7)",c='k', linestyle='--') -ax.set_xlabel(u'Time(s)') +ax.set_xlabel(u'Time (s)') ax.set_ylabel(r"Flux at outer surface (atoms/m$^2$/s)") ax.legend(loc="best") ax.set_xlim(left=0) diff --git a/test/tests/ver-1b/comparison_ver-1b.py b/test/tests/ver-1b/comparison_ver-1b.py index 0a297460..3d8e8da3 100644 --- a/test/tests/ver-1b/comparison_ver-1b.py +++ b/test/tests/ver-1b/comparison_ver-1b.py @@ -33,7 +33,7 @@ ax.plot(analytical_time,analytical_conc,label=r"Analytical",c='k', linestyle='--') -ax.set_xlabel(u'Time(s)') +ax.set_xlabel(u'Time (s)') ax.set_ylabel(r"Normalized specie concentration") ax.legend(loc="best") ax.set_xlim(left=0) @@ -108,7 +108,7 @@ ax.plot(analytical_time[1:],analytical_flux[1:],label=r"Analytical",c='k', linestyle='--') -ax.set_xlabel(u'Time(s)') +ax.set_xlabel(u'Time (s)') ax.set_ylabel(r"Diffusive flux") ax.legend(loc="best") ax.set_xlim(left=0) diff --git a/test/tests/ver-1c/comparison_ver-1c.py b/test/tests/ver-1c/comparison_ver-1c.py index 437073bf..79cf2bb9 100644 --- a/test/tests/ver-1c/comparison_ver-1c.py +++ b/test/tests/ver-1c/comparison_ver-1c.py @@ -84,7 +84,7 @@ def get_c_analytical_7(x,t): ax.text(20,get_c_analytical_4(12,20), 'RMSPE = %.2f '%RMSPE+'%',fontweight='bold',ha='left',va='top') -ax.set_xlabel(u'Time(s)') +ax.set_xlabel(u'Time (s)') ax.set_ylabel(r"Species concentration (atoms/m$^3$)") ax.legend(loc="best") ax.set_xlim(left=0,right=100) @@ -119,7 +119,7 @@ def get_c_analytical_7(x,t): RMSPE = RMSE*100/np.mean(analytical_conc12_7[idx:]) ax.text(20,get_c_analytical_7(12,20), 'RMSPE = %.2f '%RMSPE+'%',fontweight='bold') -ax.set_xlabel(u'Time(s)') +ax.set_xlabel(u'Time (s)') ax.set_ylabel(r"Species concentration (atoms/m$^3$)") ax.legend(loc="best") ax.set_xlim(left=0,right=100) diff --git a/test/tests/ver-1d/comparison_ver-1d.py b/test/tests/ver-1d/comparison_ver-1d.py index 337496b2..2a3d3f0c 100644 --- a/test/tests/ver-1d/comparison_ver-1d.py +++ b/test/tests/ver-1d/comparison_ver-1d.py @@ -79,7 +79,7 @@ def summation_term(num_terms, time): ax.plot(tmap_time, tmap_perm, label=r"TMAP8", c='tab:gray') -ax.set_xlabel(u'Time(s)') +ax.set_xlabel(u'Time (s)') ax.set_ylabel(u"Permeation (atom/m$^2$s)") ax.legend(loc="best") ax.set_xlim(left=0) @@ -119,7 +119,7 @@ def summation_term(num_terms, time): print(tmap_time[np.argmin(np.abs(tmap_perm-3.1622e18))]) -ax.set_xlabel(u'Time(s)') +ax.set_xlabel(u'Time (s)') ax.set_ylabel(u"Permeation (atom/m$^2$s)") ax.legend(loc="lower right") ax.set_xlim(left=0) diff --git a/test/tests/ver-1dc/comparison_ver-1dc.py b/test/tests/ver-1dc/comparison_ver-1dc.py index 43ba956b..baee3071 100644 --- a/test/tests/ver-1dc/comparison_ver-1dc.py +++ b/test/tests/ver-1dc/comparison_ver-1dc.py @@ -93,7 +93,7 @@ def summation_term(num_terms, time): ax.plot([tau_be, tau_be], [0, 3.15e18], label=r"Analytical breakthrough time", c='tab:brown', linestyle='--') ax.plot([tmap_intercept, 17.5], [0, (17.5 - tmap_intercept) * np.max(tmap_slope)], label=r"Numerical breakthrough time", c='tab:brown') -ax.set_xlabel(u'Time(s)') +ax.set_xlabel(u'Time (s)') ax.set_ylabel(u"Flux (atom/m$^2$s)") ax.legend(loc="best") ax.set_xlim(left=0) diff --git a/test/tests/ver-1dd/comparison_ver-1dd.py b/test/tests/ver-1dd/comparison_ver-1dd.py index d79534eb..97bad40f 100644 --- a/test/tests/ver-1dd/comparison_ver-1dd.py +++ b/test/tests/ver-1dd/comparison_ver-1dd.py @@ -55,7 +55,7 @@ def summation_term(num_terms, time): ax.plot([tmap_intercept, 0.221], [0, (0.221 - tmap_intercept) * np.max(tmap_slope)], label=r"Numerical breakthrough time", c='tab:brown') ax.plot([tau_be, tau_be], [0, 3.15e18], label=r"Analytical breakthrough time", c='tab:brown', linestyle='--') -ax.set_xlabel(u'Time(s)') +ax.set_xlabel(u'Time (s)') ax.set_ylabel(u"Flux (atom/m$^2$s)") ax.legend(loc="best") ax.set_xlim(left=0, right=1.5) diff --git a/test/tests/ver-1e/comparison_ver-1e.py b/test/tests/ver-1e/comparison_ver-1e.py index a98abae2..0e3f1bfd 100644 --- a/test/tests/ver-1e/comparison_ver-1e.py +++ b/test/tests/ver-1e/comparison_ver-1e.py @@ -132,8 +132,9 @@ def get_lambdas(k,l,a): ax.minorticks_on() plt.savefig('ver-1e_comparison_time_closeup.png', bbox_inches='tight', dpi=300) plt.close(fig) -# ============ Comparison of concentration as a function of distance ============ + +# ============ Comparison of concentration as a function of distance ============ fig = plt.figure(figsize=[6.5, 5.5]) gs = gridspec.GridSpec(1, 1) ax = fig.add_subplot(gs[0]) @@ -197,6 +198,8 @@ def get_lambdas(k,l,a): ax.set_xlabel(u'Distance ($\mu$m)') ax.set_ylabel(r"Concentration (moles/m$^3$)") +ax.set_xlim(left=0) +ax.set_ylim(bottom=0) ax.legend(loc="best") plt.grid(visible=True, which='major', color='0.65', linestyle='--', alpha=0.3) diff --git a/test/tests/ver-1fa/comparison_ver-1fa.py b/test/tests/ver-1fa/comparison_ver-1fa.py index 37644e16..295f9f9d 100644 --- a/test/tests/ver-1fa/comparison_ver-1fa.py +++ b/test/tests/ver-1fa/comparison_ver-1fa.py @@ -58,6 +58,8 @@ def get_analytical_solution(x): ax.set_xlabel(u'Distance along slab (m)') ax.set_ylabel(u"Temperature (K)") ax.legend(loc="best") +ax.set_xlim(left=0) +ax.set_xlim(right=max(tmap_x)) ax.set_ylim(bottom=0) plt.grid(visible=True, which='major', color='0.65', linestyle='--', alpha=0.3) RMSE = np.sqrt(np.mean((tmap_temp-analytical_temp)**2) ) diff --git a/test/tests/ver-1fb/comparison_ver-1fb.py b/test/tests/ver-1fb/comparison_ver-1fb.py index 64413f42..dcbdd084 100644 --- a/test/tests/ver-1fb/comparison_ver-1fb.py +++ b/test/tests/ver-1fb/comparison_ver-1fb.py @@ -115,6 +115,8 @@ def get_analytical_solution(x, time_i): ax.set_xlabel(u'Distance along slab (m)') ax.set_ylabel(u"Temperature (K)") ax.legend(loc="best") +ax.set_xlim(left=0) +ax.set_xlim(right=max(tmap_x)) ax.set_ylim(bottom=300) plt.grid(visible=True, which='major', color='0.65', linestyle='--', alpha=0.3) RMSE = np.sqrt(np.mean((tmap_temp[0]-analytical_temp[0])**2) ) # 0.1 seconds