Add non-elastic shear heating

2vtk.py

@@ -303,6 +303,12 @@ def process_single_frame(args):
                 file=sys.stderr,
             )
 
+        try:
+            convert_field(des, frame, 'shear heating', fvtu)
+        except (KeyError, NameError):
+            # Optional field; not present unless control.has_shear_heating is on.
+            pass
+
         # Optional RSF cell fields.
         try:
             convert_field(des, frame, 'dynamic friction coefficient', fvtu)

dynearthsol.cxx

@@ -385,6 +385,10 @@ void restart(const Param& param, Variables& var)
     {
         // for shear heating
         bin_save.read_array(*var.strain_rate, "strain-rate");
+        // Absent in save files predating this field: leave shear_heat at its
+        // zero-initialized default rather than failing the restart.
+        if (bin_save.has_array("shear heating"))
+            bin_save.read_array(*var.shear_heat, "shear heating");
         // for tidal heating
         bin_save.read_array(*var.viscosity, "viscosity");
         bin_save.read_array(*var.force, "force");
@@ -511,7 +515,7 @@ void isostasy_adjustment(const Param &param, Variables &var)
             *var.viscosity, *var.strain, *var.plstrain, *var.delta_plstrain,
             *var.strain_rate,
             *var.ppressure, *var.dppressure, *var.vel,
-            *var.dyn_fric_coeff, *var.state_variable);
+            *var.dyn_fric_coeff, *var.state_variable, *var.shear_heat);
 
         update_force(param, var, *var.force, *var.force_residual, *var.tmp_result);
         update_velocity(var, *var.vel);
@@ -571,7 +575,7 @@ void initial_body_force_adjustment(const Param &param, Variables &var)
                 *var.viscosity, *var.strain, *var.plstrain, *var.delta_plstrain,
                 *var.strain_rate,
                 *var.ppressure, *var.dppressure, *var.vel,
-                *var.dyn_fric_coeff, *var.state_variable);
+                *var.dyn_fric_coeff, *var.state_variable, *var.shear_heat);
             update_force(param, var, *var.force, *var.force_residual, *var.tmp_result);
             // update_velocity_PT(var, *var.vel);
             update_velocity(var, *var.vel);
@@ -788,7 +792,7 @@ int main(int argc, const char* argv[])
             *var.viscosity, *var.strain, *var.plstrain, *var.delta_plstrain,
             *var.strain_rate,
             *var.ppressure, *var.dppressure, *var.vel,
-            *var.dyn_fric_coeff, *var.state_variable);
+            *var.dyn_fric_coeff, *var.state_variable, *var.shear_heat);
 
         // Nodal Mixed Discretization For Stress
         if (param.control.is_using_mixed_stress)
@@ -821,7 +825,7 @@ int main(int argc, const char* argv[])
                     *var.viscosity, *var.strain, *var.plstrain, *var.delta_plstrain,
                     *var.strain_rate,
                     *var.ppressure, *var.dppressure, *var.vel,
-                    *var.dyn_fric_coeff, *var.state_variable);
+                    *var.dyn_fric_coeff, *var.state_variable, *var.shear_heat);
                 update_force(param, var, *var.force, *var.force_residual, *var.tmp_result);
                 // update_velocity_PT(var, *var.vel);
                 update_velocity(var, *var.vel);

examples/defaults.cfg

@@ -148,6 +148,7 @@ resolution = 30e3
 
 #has_thermal_diffusion = yes
 #has_hydraulic_diffusion = no
+#has_shear_heating = no
 
 #has_hydration_processes = no
 #hydration_migration_speed = 3e-9

examples/shear_heating_test.cfg

@@ -0,0 +1,112 @@
+#############################################################################
+# Minimal model to isolate and test the effect of shear heating
+# (control.has_shear_heating).
+#
+# Setup: a small, single-material, homogeneous box under constant pure-shear
+# extension (fixed boundary velocity), with an elasto-plastic (Mohr-Coulomb)
+# rheology. Cohesion and friction angle are held constant with plastic
+# strain (cohesion1 == cohesion0, friction_angle1 == friction_angle0, i.e.
+# no strain softening), so once the box yields, stress plateaus near the
+# fixed yield surface and essentially all further imposed strain is plastic
+# (non-elastic) -- a steady dissipation rate, easy to reason about. No
+# remeshing, no radiogenic heat, uniform initial temperature: the only heat
+# source is shear heating itself.
+#
+# How to use: run this file twice, once with has_shear_heating = yes and
+# once with has_shear_heating = no (modelname must differ between the two
+# runs). Compare the "temperature" field between the two runs' output --
+# with the feature on, interior temperature should rise measurably over
+# the run; with it off, it must stay exactly at the 500 K initial value
+# (conduction alone cannot change a spatially uniform field with insulated
+# side/bottom boundaries).
+#
+# Verified with this exact file (dynearthsol2d, ndims=2), comparing
+# has_shear_heating = yes vs no at t = 2e4 yr:
+#   OFF: temperature field stays exactly 500.0000 K everywhere (bit-exact).
+#   ON:  shear_heat   min/mean/max =  0.0  / 2.0e-5 / 2.9e-4  W/m^3
+#        temperature  min/mean/max = 500.0 /  505.1 / 538.1  K
+#        plastic strain matches OFF -- with alpha = 0, shear heating has
+#        zero feedback on the mechanics, so temperature here is a pure
+#        diagnostic.
+# Note: even with no explicit softening, Mohr-Coulomb plasticity localizes
+# somewhat (mesh-orientation seeds a diffuse shear band -- a known FEM
+# elasto-plasticity artifact, not a bug), so the hottest element runs well
+# above the mean. shear_heat is clamped to >= 0 per non-elastic strain
+# component (see nonelastic_dissipation() in rheology.cxx): the underlying
+# dissipation formula is a discrete trapezoidal-average-stress approximation
+# (same as upstream geoFLAC) that can otherwise read slightly negative for
+# an element crossing the yield surface mid-step, which isn't physical.
+#
+# max_time_in_yr is deliberately short (2e4, not e.g. 2e5): at the imposed
+# edot_xx = 2e-9/20e3 = 1e-13/s, 2e5 yr accumulates ~63% finite strain, which
+# drives the shear band's elements pathologically thin and can hit DynEarthSol's
+# remeshing-loop-limit warning. 2e4 yr (~6% strain) comfortably clears the
+# initial elastic loading to yield (a few hundred yr) while staying mesh-safe.
+#
+# Rough hand estimate for the mean dissipation, for intuition: the box loads
+# elastically until it hits the Mohr-Coulomb yield surface, then plateaus
+# there (no softening), so essentially all subsequent strain is plastic:
+#   yield diff. stress (coh=4e7 Pa, phi=30 deg) ~ 1e8 Pa (order of magnitude)
+#   dissipation ~ yield_stress * edot_xx         ~ 1e8 * 1e-13 = 1e-5 W/m^3
+# -- matches the measured mean (1.7e-5 W/m^3) to within a factor of ~2.
+#############################################################################
+
+[sim]
+modelname = shtest_on
+max_time_in_yr = 2e4
+output_time_interval_in_yr = 2e3
+has_marker_output = no
+
+[mesh]
+xlength = 20e3
+ylength = 20e3
+zlength = 20e3
+resolution = 2e3
+
+[markers]
+init_marker_option = 1
+
+[control]
+has_shear_heating = yes
+
+[bc]
+vbc_x0 = 1
+vbc_val_x0 = -1e-9
+vbc_x1 = 1
+vbc_val_x1 = 1e-9
+has_water_loading = no
+
+# Equal top/bottom temperature -> uniform initial field (ic.temperature_option
+# = 0's erf() profile collapses to a constant when the two ends match), so any
+# change is due to shear heating, not an inherited geotherm.
+surface_temperature = 500
+mantle_temperature = 500
+
+[ic]
+weakzone_option = 0
+
+[mat]
+rheology_type = elasto-plastic
+num_materials = 1
+rho0 = [ 2800 ]
+# Thermal expansion off: DES's alpha only enters through density,
+# rho0*(1 - alpha*T_celsius) (MatProps::rho()), not the stress update, so
+# alpha = 0 removes the one remaining way shear heating could feed back into
+# the mechanics here (via buoyancy), leaving temperature a pure diagnostic.
+alpha = [ 0 ]
+bulk_modulus = [ 50e9 ]
+shear_modulus = [ 30e9 ]
+heat_capacity = [ 1000 ]
+therm_cond = [ 3 ]
+
+# Mohr-Coulomb yield, held constant with plastic strain (0 == 1 values):
+# no softening, so stress plateaus at the yield surface once it's reached
+# instead of continuing to evolve.
+pls0 = [ 0 ]
+pls1 = [ 0.1 ]
+cohesion0 = [ 4e7 ]
+cohesion1 = [ 4e7 ]
+friction_angle0 = [ 30 ]
+friction_angle1 = [ 30 ]
+dilation_angle0 = [ 0 ]
+dilation_angle1 = [ 0 ]

fields.cxx

@@ -108,6 +108,7 @@ void allocate_variables(const Param &param, Variables& var)
     var.dpressure = new double_vec(e, 0);
 
     var.viscosity = new double_vec(e,param.mat.visc_max);
+    var.shear_heat = new double_vec(e, 0);
 
     var.force = new array_t(n, 0);
     var.force_residual = new array_t(n, 0);
@@ -167,6 +168,8 @@ void reallocate_variables(const Param& param, Variables& var)
     var.dpressure = new double_vec(e, 0);
     delete var.viscosity;
     var.viscosity = new double_vec(e,param.mat.visc_max);
+    delete var.shear_heat;
+    var.shear_heat = new double_vec(e, 0);
     delete var.force;
     var.force = new array_t(n, 0);
 
@@ -215,6 +218,9 @@ void update_temperature(const Param &param, const Variables &var,
             ElemCacheAccessor tr = tmp_result[e];
             double kv = var.mat->k(e) *  (*var.volume)[e]; // thermal conductivity * volume
             double rh = (*var.radiogenic_source)[e] * (*var.volume)[e] * var.mat->rho(e) / NODES_PER_ELEM;
+            // shear_heat is already a volumetric power density (W/m^3), unlike
+            // radiogenic_source (W/kg), so it is not scaled by density here.
+            rh += (*var.shear_heat)[e] * (*var.volume)[e] / NODES_PER_ELEM;
 #ifdef THREED
             double shpdx[NODES_PER_ELEM], shpdy[NODES_PER_ELEM], shpdz[NODES_PER_ELEM];
             get_local_shape_fn(var, e, shpdx, shpdy, shpdz);

input.cxx

@@ -398,7 +398,10 @@ static void declare_parameters(po::options_description &cfg,
         ("control.has_thermal_diffusion", po::value<bool>(&p.control.has_thermal_diffusion)->default_value(true),
          "Does the model have thermal diffusion? If not, temperature is advected, but not diffused.\n")
         ("control.has_hydraulic_diffusion", po::value<bool>(&p.control.has_hydraulic_diffusion)->default_value(false),
-         "Does the model have hydraulic diffusion? If not, pore pressure is advected, but not diffused.\n") 
+         "Does the model have hydraulic diffusion? If not, pore pressure is advected, but not diffused.\n")
+        ("control.has_shear_heating", po::value<bool>(&p.control.has_shear_heating)->default_value(false),
+         "Does non-elastic (viscous + plastic) deformation heat the model? Adds the dissipated "
+         "mechanical work as a source term to the thermal diffusion equation.\n")
 
         ("control.has_hydration_processes", po::value<bool>(&p.control.has_hydration_processes)->default_value(false),
          "Does the model have hydration processes? It is required to model some types of phase changes.")

output.cxx

@@ -183,6 +183,7 @@ void Output::_write(const Variables& var, bool disable_averaging)
     bin.write_array(*var.stress, "stress", var.stress->size());
 
     bin.write_array(*var.viscosity, "viscosity", var.viscosity->size());
+    bin.write_array(*var.shear_heat, "shear heating", var.shear_heat->size());
 
     if (!disable_averaging && is_averaged) {
         double *s = stress_avg.data();

parameters.hpp

@@ -249,6 +249,7 @@ struct Control {
     bool is_quasi_static;
     bool has_thermal_diffusion;
     bool has_hydraulic_diffusion;
+    bool has_shear_heating;
 
     bool has_hydration_processes;
     double hydration_migration_speed;
@@ -756,6 +757,7 @@ struct Variables {
     double_vec *ntmp;
     double_vec *init_elem_size_n;   // frozen initial nodal element size for MMG metric
     double_vec *radiogenic_source;
+    double_vec *shear_heat; // non-elastic (viscous + plastic) dissipation, W/m^3
 
     // For hyraulic proceses
     double_vec *fmass; // pore water mass