Added SolubilityRatioMaterial to handle dissimilar solid material int…

…erfaces ('linked segments' in legacy TMAP parlance). Also added problem val-2d (Refs idaholab#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

@@ -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<Real> & _solubility_primary;
+  const MaterialProperty<Real> & _solubility_secondary;
+  const VariableValue & _concentration_primary;
+  const VariableValue & _concentration_secondary;
+  MaterialProperty<Real> & _jump;
+};

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<std::string>(
+      "solubility_primary", "The material property on the primary side of the interface");
+  params.addRequiredParam<std::string>(
+      "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<std::string>("solubility_primary")),
+    _solubility_secondary_name(getParam<std::string>("solubility_secondary")),
+    _solubility_primary(getMaterialPropertyByName<Real>(_solubility_primary_name)),
+    _solubility_secondary(getNeighborMaterialPropertyByName<Real>(_solubility_secondary_name)),
+    _concentration_primary(coupledValue("concentration_primary")),
+    _concentration_secondary(coupledNeighborValue("concentration_secondary")),
+    _jump(declareProperty<Real>("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);
+}

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)

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
+[]