From 7b6d2af1d14fa244e0d56e9d3978ec5554494c24 Mon Sep 17 00:00:00 2001 From: Humrickhouse Date: Wed, 21 Sep 2022 11:23:03 -0400 Subject: [PATCH] Added SolubilityRatioMaterial to handle dissimilar solid material interfaces ('linked segments' in legacy TMAP parlance). Also added problem val-2d (Refs #12) though it isn't working yet, and the digitized experiment data and python script needed to make the comparison figure. --- include/materials/SolubilityRatioMaterial.h | 35 +++ src/materials/SolubilityRatioMaterial.C | 48 ++++ test/tests/val-2b/comparison_val-2b.py | 43 ++++ test/tests/val-2b/val-2b.i | 235 ++++++++++++++++++++ 4 files changed, 361 insertions(+) create mode 100644 include/materials/SolubilityRatioMaterial.h create mode 100644 src/materials/SolubilityRatioMaterial.C create mode 100644 test/tests/val-2b/comparison_val-2b.py create mode 100644 test/tests/val-2b/val-2b.i diff --git a/include/materials/SolubilityRatioMaterial.h b/include/materials/SolubilityRatioMaterial.h new file mode 100644 index 00000000..88d54ae0 --- /dev/null +++ b/include/materials/SolubilityRatioMaterial.h @@ -0,0 +1,35 @@ +//* This file is part of the MOOSE framework +//* https://www.mooseframework.org +//* +//* All rights reserved, see COPYRIGHT for full restrictions +//* https://github.com/idaholab/moose/blob/master/COPYRIGHT +//* +//* Licensed under LGPL 2.1, please see LICENSE for details +//* https://www.gnu.org/licenses/lgpl-2.1.html + +#pragma once + +#include "InterfaceMaterial.h" + +/** + * Calculates the jump in concentration across an interface + * based on the ratio of solubilities + */ +class SolubilityRatioMaterial : public InterfaceMaterial +{ +public: + static InputParameters validParams(); + + SolubilityRatioMaterial(const InputParameters & parameters); + +protected: + virtual void computeQpProperties() override; + + const std::string _solubility_primary_name; + const std::string _solubility_secondary_name; + const MaterialProperty & _solubility_primary; + const MaterialProperty & _solubility_secondary; + const VariableValue & _concentration_primary; + const VariableValue & _concentration_secondary; + MaterialProperty & _jump; +}; diff --git a/src/materials/SolubilityRatioMaterial.C b/src/materials/SolubilityRatioMaterial.C new file mode 100644 index 00000000..9b924646 --- /dev/null +++ b/src/materials/SolubilityRatioMaterial.C @@ -0,0 +1,48 @@ +//* This file is part of the MOOSE framework +//* https://www.mooseframework.org +//* +//* All rights reserved, see COPYRIGHT for full restrictions +//* https://github.com/idaholab/moose/blob/master/COPYRIGHT +//* +//* Licensed under LGPL 2.1, please see LICENSE for details +//* https://www.gnu.org/licenses/lgpl-2.1.html + +#include "SolubilityRatioMaterial.h" + +registerMooseObject("TMAPApp", SolubilityRatioMaterial); + +InputParameters +SolubilityRatioMaterial::validParams() +{ + InputParameters params = InterfaceMaterial::validParams(); + params.addClassDescription("Calculates the jump in concentration across an interface."); + params.addRequiredParam( + "solubility_primary", "The material property on the primary side of the interface"); + params.addRequiredParam( + "solubility_secondary", "The material property on the secondary side of the interface"); + params.addRequiredCoupledVar("concentration_primary", + "Primary side non-linear variable for jump computation"); + params.addRequiredCoupledVar("concentration_secondary", + "Secondary side non-linear variable for jump computation"); + return params; +} + +SolubilityRatioMaterial::SolubilityRatioMaterial(const InputParameters & parameters) + : InterfaceMaterial(parameters), + _solubility_primary_name(getParam("solubility_primary")), + _solubility_secondary_name(getParam("solubility_secondary")), + _solubility_primary(getMaterialPropertyByName(_solubility_primary_name)), + _solubility_secondary(getNeighborMaterialPropertyByName(_solubility_secondary_name)), + _concentration_primary(coupledValue("concentration_primary")), + _concentration_secondary(coupledNeighborValue("concentration_secondary")), + _jump(declareProperty("solubility_ratio")) +{ +} + +void +SolubilityRatioMaterial::computeQpProperties() +{ + mooseAssert(_neighbor_elem, "Neighbor elem is NULL!"); + //* _jump[_qp] = _concentration_primary[_qp] - _concentration_secondary[_qp]; + _jump[_qp] = _concentration_primary[_qp]*(_solubility_secondary[_qp]/_solubility_primary[_qp] - 1); +} diff --git a/test/tests/val-2b/comparison_val-2b.py b/test/tests/val-2b/comparison_val-2b.py new file mode 100644 index 00000000..daf25a3c --- /dev/null +++ b/test/tests/val-2b/comparison_val-2b.py @@ -0,0 +1,43 @@ +import csv +import matplotlib.pyplot as plt +import numpy as np +from matplotlib import gridspec +import pandas as pd +from scipy import special + +#============ Comparison of 'Material Diffusion Experiment ===================== +#==================(Beryllium absorption/desorption) =========================== + +# Reference 1: +# R.G. Macaulay-Newcombe, D.A. Thompson, and W.W. Smeltzer, "Deuterium diffusion, +# trapping and release in ion-implanted beryllium," Fusion Engineering and Design +# 18 (1991) 419-424. + +# Reference 2: +# GR Longhurst, SL Harms, ES Marwil, and BG Miller. Verification and validation of +# tmap4. Technical Report, EG and G Idaho, Inc., Idaho Falls, ID (United States), 1992. + +fig = plt.figure(figsize=[6.5,3.9]) +gs = gridspec.GridSpec(1,1) +ax = fig.add_subplot(gs[0]) + +exp_data = pd.read_csv("val-2b-expdata.csv") +ax.plot(exp_data['Temperature'],exp_data['Flux'], label=r"Experimental Data", c='black', marker='+', linestyle='none') + +#tmap_sol = pd.read_csv("./gold/val-2b_out.csv") +#tmap_temp = tmap_sol['temp'] +#tmap_flux = tmap_sol['scaled_outflux'] +#ax.plot(tmap_temp,tmap_flux,label=r"TMAP8",c='black') + +ax.set_xlabel(u'Temperature (\u00B0C)') +ax.set_ylabel(u"Desorbed D ($10^{15}$ atom/m$^2\cdot$s)") +ax.legend(loc="best") +ax.set_xlim(left=400) +ax.set_xlim(right=800) +ax.set_ylim(bottom=-10) +ax.set_ylim(top=70) +plt.grid(b=True, which='major', color='0.65', linestyle='--', alpha=0.3) + +ax.minorticks_on() +plt.savefig('val-2b_comparison.png', bbox_inches='tight'); +plt.close(fig) diff --git a/test/tests/val-2b/val-2b.i b/test/tests/val-2b/val-2b.i new file mode 100644 index 00000000..3a210fdd --- /dev/null +++ b/test/tests/val-2b/val-2b.i @@ -0,0 +1,235 @@ +#cl=3.1622e18 +cl=1.0 + +[Mesh] + [./cmg] + type = CartesianMeshGenerator + dim = 1 + dx = '1e-9 1e-9 1e-9 1e-9 1e-9 1e-9 1e-9 1e-9 1e-9 1e-9 1e-9 1e-9 1e-9 1e-9 1e-9 1e-9 1e-9 1e-9 + 1e-9 1e-8 1e-7 1e-6 1e-5 1.888e-5 1.888e-5 1.888e-5 1.888e-5 1.888e-5 1.888e-5 1.888e-5 1.888e-5 1.888e-5 1.888e-5' + subdomain_id = '0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1' + [../] + [./subdomain_id] + input = cmg + type = SubdomainBoundingBoxGenerator + bottom_left = '18e-9 0 0' + top_right = '1.99929e-4 0 0' # sum of all dx's + block_id = 1 + [../] + [./interface] + type = SideSetsBetweenSubdomainsGenerator + input = subdomain_id + primary_block = '0' + paired_block = '1' + new_boundary = 'interface' + [../] +[] + +[Variables] + [./C_BeO] + block = 0 + [../] + [./C_Be] + block = 1 + [../] + [./T_BeO] + block = 0 + [../] + [./T_Be] + block = 1 + [../] +[] + +[Kernels] + [./diff_BeO] + type = FunctionDiffusion + variable = C_BeO + function = diffusivity_BeO + block = 0 + [../] + [./time_BeO] + type = TimeDerivative + variable = C_BeO + block = 0 + [../] + [./diff_Be] + type = FunctionDiffusion + variable = C_Be + function = diffusivity_Be + block = 1 + [../] + [./time_Be] + type = TimeDerivative + variable = C_Be + block = 1 + [../] + [./heat_BeO] + type = HeatConduction + variable = T_BeO + diffusion_coefficient = thermal_conductivity_BeO + block = 0 + [../] + [./time_heat_BeO] + type = SpecificHeatConductionTimeDerivative + variable = T_BeO + density = density_BeO + specific_heat = specific_heat_BeO + block = 0 + [../] + [./heat_Be] + type = HeatConduction + variable = T_Be + diffusion_coefficient = thermal_conductivity_Be + block = 1 + [../] + [./time_heat_Be] + type = SpecificHeatConductionTimeDerivative + variable = T_Be + block = 1 + [../] +[] + +[InterfaceKernels] + [tied] + type = PenaltyInterfaceDiffusion + variable = C_BeO + neighbor_var = C_Be + jump_prop_name = "solubility_ratio" + penalty = 1e6 + boundary = 'interface' + [] +[] + +[BCs] + [C_left] + type = EquilibriumBC + K = ${solubility_BeO} + boundary = left + enclosure_scalar_var = enclosure_pressure + temp = ${T_BeO} + variable = C_BeO + p = 0.5 + [] + [C_right] + type = NeumannBC + variable = C_Be + value = 0 + boundary = right + [] + [T_left] + type = FunctionDirichletBC + variable = T_BeO + boundary = left + function = temperature_history + [] + [T_right] + type = NeumannBC + variable = T_Be + value = 0 + boundary = right + [] +[] + +[Functions] + [enclosure_pressure] + type = ParsedFunction + value = 'if(t<180015.0, 13300.000001, if(t<182400.0, 1e-6, 0.001))' #========================= + [] + [temperature_history] + type = ParsedFunction + value = 'if(t<180000.0, 773.0, if(t<182400.0, 773.0-((1-exp(-(t-180000)/2700))*475), 300+0.05*(t-182400)))' #========================= + [] +[] + +[Materials] + [./BeO_d] + type = ParsedMaterial + f_name = 'diffusivity_BeO' + args = 'T_BeO' + function = 'if(t<182400, 1.40e-4*exp(-24408/T_BeO), 7e-5*exp(-24408/T_BeO))' + [../] + [./BeO_s] + type = ParsedMaterial + f_name = 'solubility_BeO' + args = 'T_BeO' + function = '5.00e20*exp(9377.7/T_BeO)' + [../] + [./BeO_HT] + type = GenericConstantMaterial + prop_names = 'density_BeO thermal_conductivity_BeO specific_heat_BeO' + prop_values = '3010.0 159.2 996.67774' + [../] + [./Be_d] + type = ParsedMaterial + f_name = 'diffusivity_Be' + args = 'T_Be' + function = '8.0e-9*exp(-4220/T_Be)' + [../] + [./Be_s] + type = ParsedMaterial + f_name = 'solubility_Be' + args = 'T_Be' + function = '7.156e27*exp(-11606/T_Be)' + [../] + [./Be_HT] + type = GenericConstantMaterial + prop_names = 'density_Be thermal_conductivity_Be specific_heat_Be' + prop_values = '1850 168.0 1821.62162' + [../] + [./interface_jump] + type = SolubilityRatioMaterial + mat_prop_primary = solubility_BeO + mat_prop_secondary = solubility_Be + boundary = interface + concentration_primary = C_BeO + concentration_secondary = C_Be + [../] +[] + +[Postprocessors] + [outflux] + type = SideDiffusiveFluxAverage + boundary = 'left' + diffusivity = diffusivity_BeO + variable = C_BeO + [] + [scaled_outflux] + type = ScalePostprocessor + value = outflux + scaling_factor = ${cl} + [] +[] + +[Preconditioning] + [smp] + type = SMP + full = true + [] +[] + +[Executioner] + type = Transient + end_time = 197860 + dt = 60.0 + dtmin = .01 + solve_type = NEWTON + petsc_options_iname = '-pc_type -pc_hypre_type' + petsc_options_value = 'hypre boomeramg' + automatic_scaling = true + verbose = true + compute_scaling_once = false +[] + +[Outputs] + csv = true + [out] + type = Exodus + execute_on = 'initial timestep_end' + [] + [dof] + type = DOFMap + execute_on = initial + [] + perf_graph = true +[]