Updates to concrete moisture transport models

doc/content/source/materials/ConcreteThermalMoisture.md

@@ -358,8 +358,8 @@ $W$ = total water content (g/g)\\
 $H$ = pore relative humidity, and $H=P_v$/$P_{vs}$ \\
 $P_{vs}$ = saturate vapor pressure $= P_{atm}\exp\left(4871.3\frac{T-100}{373.15T}\right)$ [!cite](Bary2012) (where $T$ is the temperature in K))\\
 $P_{atm}$ = standard atmospheric pressure $= 101.325 Pa$ \\
-$D_h$ = moisture diffusivity (also referred as humidity diffusivity), $m^2/s$\\
-$D_{ht}$= coupled moisture diffusivity under the influence of a temperature gradient, $m^2/s$\\
+$D_h$ = moisture diffusivity (also referred as humidity diffusivity), $kg/(m^2\cdot s)$\\
+$D_{ht}$= coupled moisture diffusivity under the influence of a temperature gradient, $kg/(m^2\cdot s\cdot K)$\\
 $W_d$= total mass of free evaporable water released into the pores by dehydration of the cement paste\\
 $t$  = time, $s$
 

src/bcs/SpecifiedVaporPressureBC.C

@@ -22,8 +22,8 @@ SpecifiedVaporPressureBC::validParams()
   InputParameters params = NodalBC::validParams();
   params.set<Real>("duration") = 0.0;
   params.addRequiredParam<Real>("vapor_pressure", "in Pa");
-  params.addParam<Real>("T_ref", 20.0, "initial temperature");
-  params.addParam<Real>("rh_ref", 0.96, "initial humidity");
+  params.addParam<Real>("T_ref", 20.0, "Initial temperature in C");
+  params.addParam<Real>("rh_ref", 0.96, "initial relative humidity");
   params.addCoupledVar("temperature", "nonlinear variable name holding temperature field");
   params.addClassDescription(
       "Prescribed vapor pressure boundary condition for moisture transport in concrete.");
@@ -45,9 +45,9 @@ SpecifiedVaporPressureBC::SpecifiedVaporPressureBC(const InputParameters & param
 Real
 SpecifiedVaporPressureBC::computeQpResidual()
 {
-  Real T = _T_ref;
+  Real T_ref_K = _T_ref + 273.15;
 
-  Real p_sat0 = 101325.0 * std::exp(4871.3 * (_T_ref - 100.0) / 373.15 / (T + 273.15));
+  Real p_sat0 = 101325.0 * std::exp(4871.3 * (T_ref_K - 100.0) / 373.15 / T_ref_K);
   _initial = _rh_ref * p_sat0;
   _final = _vapor_pressure;
 
@@ -58,10 +58,11 @@ SpecifiedVaporPressureBC::computeQpResidual()
   else
     value = _final;
 
+  Real T = T_ref_K;
   if (_has_temperature)
-    T = _temp[_qp];
+    T = _temp[_qp] + 273.15;
 
-  Real p_sat = 101325.0 * std::exp(4871.3 * (T - 100.0) / 373.15 / (T + 273.15));
+  Real p_sat = 101325.0 * std::exp(4871.3 * (T - 100.0) / 373.15 / T);
   Real rh = value / p_sat;
 
   return _u[_qp] - rh;

test/tests/concrete_moisture_heat_transfer/maqbeth.i → test/tests/concrete_moisture_heat_transfer/maqbeth_1d.i

@@ -1,107 +1,109 @@
 [Problem]
-   coord_type = RZ
+  type = ReferenceResidualProblem
+  reference_vector = 'ref'
+  extra_tag_vectors = 'ref'
+  coord_type = RZ
 []
 
 [Mesh]
-  file = maqbethRZ.e
-  #uniform_refine = 1
+  file = maqbeth_1d.e
+  construct_side_list_from_node_list = true
 []
 
 [Functions]
-  [./ramp1]
+  [temp_hist]
     type = PiecewiseLinear
     x = '0.0  36000.0 114120.0 138600.0  213120'
     y = '20.0  20.0      150.0     80.0   200.0'
-    #scale_factor = 0.2
-  [../]
+  []
 []
 
 [Variables]
-  [./T]
+  [T]
     order = FIRST
     family = LAGRANGE
     initial_condition = 20.0
-  [../]
-  [./rh]
+  []
+  [rh]
     order = FIRST
     family = LAGRANGE
     initial_condition = 0.96
-  [../]
+  []
 []
 
 [AuxVariables]
-  [./thermal_conductivity]
+  [thermal_conductivity]
     order = CONSTANT
     family = Monomial
-  [../]
+  []
 
-  [./thermal_capacity]
+  [thermal_capacity]
     order = CONSTANT
     family = Monomial
-  [../]
+  []
 
-  [./moisture_capacity]
+  [moisture_capacity]
     order = CONSTANT
     family = Monomial
-  [../]
+  []
 
-  [./humidity_diffusivity]
+  [humidity_diffusivity]
     order = CONSTANT
     family = Monomial
-  [../]
+  []
 
-  [./water_content]
+  [water_content]
     order = CONSTANT
     family = Monomial
-  [../]
-  [./water_hydrated]
+  []
+  [water_hydrated]
     order = CONSTANT
     family = Monomial
-  [../]
+  []
 []
 
 [Kernels]
-  [./T_td]
+  [T_td]
     type     = ConcreteThermalTimeIntegration
     variable = T
-  [../]
-  [./T_diff]
+    extra_vector_tags = 'ref'
+  []
+  [T_diff]
     type     = ConcreteThermalConduction
     variable = T
-  [../]
+    extra_vector_tags = 'ref'
+  []
 
-  [./T_conv]
+  [T_conv]
     type     = ConcreteThermalConvection
     variable = T
     relative_humidity = rh
-  [../]
+    extra_vector_tags = 'ref'
+  []
 
-  [./T_adsorption]
+  [T_adsorption]
     type     = ConcreteLatentHeat
     variable = T
     H = rh
-  [../]
+    extra_vector_tags = 'ref'
+  []
 
-  [./rh_td]
+  [rh_td]
     type     = ConcreteMoistureTimeIntegration
     variable = rh
-  [../]
+    extra_vector_tags = 'ref'
+  []
 
-  [./rh_diff]
+  [h_diff]
     type     = ConcreteMoistureDiffusion
     variable = rh
     temperature = T
-  [../]
-
-  #[./rh_dehydration]
-  #  type     = ConcreteMoistureDehydration
-  #  variable = rh
-  #  temperature = T
-  #[../]
+    extra_vector_tags = 'ref'
+  []
 []
 
 [AuxKernels]
-  [./k]
+  [k]
     type = MaterialRealAux
     variable = thermal_conductivity
     property = thermal_conductivity
@@ -120,13 +122,13 @@
     property = moisture_capacity
     execute_on = 'timestep_end'
   [../]
-  [./rh_duff]
+  [./rh_diff]
     type = MaterialRealAux
     variable = humidity_diffusivity
     property = humidity_diffusivity
     execute_on = 'timestep_end'
   [../]
-  [./wc_duff]
+  [./wc_diff]
     type = MaterialRealAux
     variable = water_content
     property = moisture_content
@@ -144,7 +146,6 @@
   [./concrete]
     type = ConcreteThermalMoisture
     block = 1
-    # setup thermal property models and parameters
     # options available: CONSTANT ASCE-1992 KODUR-2004 EUROCODE-2004 KIM-2003
     thermal_conductivity_model =  KODUR-2004
     thermal_capacity_model     =  KODUR-2004
@@ -164,7 +165,7 @@
     concrete_cure_time          = 23.0       #curing time in (days)
 
     # options available for humidity diffusivity:
-    moisture_diffusivity_model = Bazant      #options: Bazant Xi Mensi
+    moisture_diffusivity_model =  Bazant      #options: Bazant Xi Mensi
     D1 = 3.0e-10
     aggregate_vol_fraction = 0.7             #used in Xi's moisture diffusivity model
 
@@ -181,7 +182,7 @@
     type = FunctionDirichletBC
     variable = T
     boundary = '1'
-    function = ramp1
+    function = temp_hist
   [../]
 
   [./T_right]
@@ -196,47 +197,70 @@
     type = SpecifiedVaporPressureBC
     variable = rh
     boundary = '1'
-    duration = 36000
+    duration = 3600
     vapor_pressure = 2500.0
     temperature = T
   [../]
   [./rh_right]
     type = SpecifiedVaporPressureBC
     variable = rh
     boundary = '2'
-    duration = 36000
+    duration = 3600
     vapor_pressure = 2500.0
     temperature = T
   [../]
 []
 
+[VectorPostprocessors]
+  [profiles]
+    type = LineValueSampler
+    start_point = '0.5 0 0'
+    end_point = '1.1 0 0'
+    num_points = 100
+    variable = 'T rh'
+    sort_by = id
+    outputs = csv
+  []
+[]
+
+[Preconditioning]
+  [smp]
+    type = SMP
+    full = true
+  []
+[]
+
 [Executioner]
   type       = Transient
   solve_type = 'PJFNK'
 
-  petsc_options_iname = '-pc_type -pc_hypre_type -ksp_gmres_restart -snes_ls -pc_hypre_boomeramg_strong_threshold'
-  petsc_options_value = 'hypre boomeramg 201 cubic 0.7'
+  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
+  petsc_options_value = 'lu superlu_dist'
+
+  automatic_scaling = true
 
-  dt = 100
-  num_steps = 5
+  dt = 10000
 
-  #dtmax = 600.0
-  #end_time = 900000.0
-  #
-  #[./TimeStepper]
-  #  type = SolutionTimeAdaptiveDT
-  #  dt = 100.0
-  #[../]
+  end_time = 900000.0
 
   l_max_its  = 50
   l_tol      = 1e-6
   nl_max_its = 10
   nl_rel_tol = 1e-6
   nl_abs_tol = 1e-10
+
+  [Predictor]
+    type = SimplePredictor
+    scale = 1.0
+    skip_times_old = '0.0  36000.0 114120.0 138600.0  213120'
+  []
 []
 
 [Outputs]
-  file_base  = maqbeth_out
   exodus     = true
-  sync_times = '69599.88 110400.12 121200.12 151200.00 211200.12 221400.00 240599.88 281400.12 454200.12 550800.00 701399.88'
+  [csv]
+    type = CSV
+    file_base = 'csv/out'
+    time_data = true
+  []
 []

test/tests/concrete_moisture_heat_transfer/maqbeth_1d.jou

@@ -0,0 +1,19 @@
+## /Applications/Cubit-15.2/Cubit.app/Contents/MacOS/cubitclx
+## Cubit Version 15.2
+## Cubit Build 405468
+## Revised 2016-09-23 10:47:14 -0600 (Fri, 23 Sep 2016)
+## Running 05/03/2021 03:26:23 PM
+## Command Options:
+create vertex 0.5 0 0
+create vertex 1.1 0 0
+create curve vertex 1 2
+curve 1 scheme dualbias fine size 0.015 coarse size 0.03
+mesh curve 1
+
+block 1 curve 1
+block 1 element type BAR2
+
+nodeset 1 vertex 1
+nodeset 2 vertex 2
+
+export mesh 'maqbeth_1d.e' overwrite

test/tests/concrete_moisture_heat_transfer/temp_rh_plots.py

@@ -0,0 +1,44 @@
+#!/usr/bin/env python
+
+import pandas
+import mooseutils
+import matplotlib.pyplot as plt
+
+f1 = plt.figure(1, figsize=(6,4.5))
+ax = plt.gca()
+data = mooseutils.VectorPostprocessorReader('csv/out_profiles_*.csv')
+
+t_times = [19.5, 30.5, 41.8, 61.3, 125.0, 194.6]
+rh_times = [33.4, 58.4, 66.7, 77.9, 152.8, 194.5]
+
+#T = data('T', time = 0)
+#print(T)
+#print(data.data.keys().values)
+
+#Plot a line for every output time:
+#for time in data.data.data.keys().values:
+#  data.data.data[time].plot(ax=ax,x='x',y='T')
+
+#Plot lines for selected interpolated output times:
+for t_hr in t_times:
+  t = t_hr * 3600
+  data.update(time=t)
+  data.data.data.plot(ax=ax,x='id',y='T',label=str(t_hr)+' hr')
+
+plt.xlabel('Radial Position (m)')
+plt.ylabel('Temperature ($\degree$C)')
+
+plt.savefig('temp_plot.pdf')
+
+f2 = plt.figure(2, figsize=(6,4.5))
+ax = plt.gca()
+
+for t_hr in rh_times:
+  t = t_hr * 3600
+  data.update(time=t)
+  data.data.data.plot(ax=ax,x='id',y='rh',label=str(t_hr)+' hr')
+
+plt.xlabel('Radial Position (m)')
+plt.ylabel('$h_r$')
+
+plt.savefig('rh_plot.pdf')