Fixing errors in gran/wall reff calculation and adding heat capabilities

doc/src/fix_wall_gran.rst

@@ -46,7 +46,7 @@ Syntax
          radius = cylinder radius (distance units)
 
 * zero or more keyword/value pairs may be appended to args
-* keyword = *wiggle* or *shear* or *contacts*
+* keyword = *wiggle* or *shear* or *contacts* or *temperature*
 
   .. parsed-literal::
 
@@ -59,6 +59,8 @@ Syntax
          vshear = magnitude of shear velocity (velocity units)
       *contacts* value = none
          generate contact information for each particle
+      *temperature* value = temperature
+         specify temperature of wall
 
 
 Examples
@@ -71,7 +73,7 @@ Examples
    fix 2 all wall/gran hooke 100000.0 20000.0 50.0 30.0 0.5 1 zcylinder 15.0 wiggle z 3.0 2.0
    fix 3 all wall/gran/region granular hooke 1000.0 50.0 tangential linear_nohistory 1.0 0.4 damping velocity region myBox
    fix 4 all wall/gran/region granular jkr 1e5 1500.0 0.3 10.0 tangential mindlin NULL 1.0 0.5 rolling sds 500.0 200.0 0.5 twisting marshall region myCone
-   fix 5 all wall/gran/region granular dmt 1e5 0.2 0.3 10.0 tangential mindlin NULL 1.0 0.5 rolling sds 500.0 200.0 0.5 twisting marshall damping tsuji region myCone
+   fix 5 all wall/gran/region granular dmt 1e5 0.2 0.3 10.0 tangential mindlin NULL 1.0 0.5 rolling sds 500.0 200.0 0.5 twisting marshall damping tsuji heat 10 region myCone temperature 1.0
    fix 6 all wall/gran hooke  200000.0 NULL 50.0 NULL 0.5 0 xplane -10.0 10.0 contacts
 
 Description
@@ -177,6 +179,15 @@ the clockwise direction for *vshear* > 0 or counter-clockwise for
 *vshear* < 0.  In this case, *vshear* is the tangential velocity of
 the wall at whatever *radius* has been defined.
 
+The *temperature* keyword is used to assign a temperature to the wall.
+The following value can either be a numeric value or an equal-style
+:doc:`variable <variable>`.  If the value is a variable, it should be
+specified as v_name, where name is the variable name.  In this case, the
+variable will be evaluated each timestep, and its value used to determine
+the temperature. This option must be used in conjunction with a heat
+conduction model defined in :doc:`pair_style granular <pair_granular>`
+and a compatible atom style which stores temperature and a heat flux
+as defined by the :doc:`atom_style sphere/temp <atom_style>` command.
 
 Restart, fix_modify, output, run start/stop, minimize info
 """""""""""""""""""""""""""""""""""""""""""""""""""""""""""

doc/src/fix_wall_gran_region.rst

@@ -36,12 +36,14 @@ Syntax
 
 * wallstyle = region (see :doc:`fix wall/gran <fix_wall_gran>` for options for other kinds of walls)
 * region-ID = region whose boundary will act as wall
-* keyword = *contacts*
+* keyword = *contacts* or *temperature*
 
   .. parsed-literal::
 
       *contacts* value = none
          generate contact information for each particle
+      *temperature* value = temperature
+         specify temperature of wall
 
 Examples
 """"""""
@@ -200,6 +202,16 @@ values for the 6 wall/particle coefficients than for particle/particle
 interactions.  E.g. if you wish to model the wall as a different
 material.
 
+The *temperature* keyword is used to assign a temperature to the wall.
+The following value can either be a numeric value or an equal-style
+:doc:`variable <variable>`.  If the value is a variable, it should be
+specified as v_name, where name is the variable name.  In this case, the
+variable will be evaluated each timestep, and its value used to determine
+the temperature. This option must be used in conjunction with a heat
+conduction model defined in :doc:`pair_style granular <pair_granular>`
+and a compatible atom style which stores temperature and a heat flux
+as defined by the :doc:`atom_style sphere/temp <atom_style>` command.
+
 Restart, fix_modify, output, run start/stop, minimize info
 """""""""""""""""""""""""""""""""""""""""""""""""""""""""""
 

src/GRANULAR/contact.cpp

@@ -328,20 +328,24 @@ void ContactModel::read_restart(FILE *fp)
 
 /* ---------------------------------------------------------------------- */
 
-bool ContactModel::check_contact(double rtemp)
+bool ContactModel::check_contact()
 {
   if (contact_type == WALL) {
     // Used by fix_wall_gran.cpp
+    //   radj = radius of wall
+    //   dx already provided
     rsq = lensq3(dx);
     radsum = radi;
-    if (rtemp == 0) Reff = radi;
-    else Reff = radi * rtemp/(radi + rtemp);
+    if (radj == 0) Reff = radi;
+    else Reff = radi * radj / (radi + radj);
   } else if (contact_type == WALLREGION) {
     // Used by fix_wall_gran_region.cpp
-    rsq = rtemp * rtemp;
-    radsum = radi + radi;
-    if (rtemp == 0) Reff = radi;
-    else Reff = radi * rtemp/(radi + rtemp);
+    //   radj = radius of wall
+    //   dx and r already provided
+    rsq = r * r;
+    radsum = radi;
+    if (radj == 0) Reff = radi;
+    else Reff = radi * radj / (radi + radj);
   } else {
     sub3(xi, xj, dx);
     rsq = lensq3(dx);

src/GRANULAR/contact.h

@@ -51,7 +51,7 @@ class ContactModel : protected Pointers {
   ContactModel(class LAMMPS *);
   ~ContactModel();
   void init();
-  bool check_contact(double = 0);
+  bool check_contact();
   void prep_contact();
   void calculate_forces();
   double pulloff_distance(double, double);
@@ -84,7 +84,7 @@ class ContactModel : protected Pointers {
   // Contact properties/output
   double forces[3], torquesi[3], torquesj[3], dq;
 
-  double radi, radj, rwall, meff, dt, Ti, Tj, area;
+  double radi, radj, meff, dt, Ti, Tj, area;
   double Fntot, magtortwist;
 
   double *xi, *xj, *vi, *vj, *omegai, *omegaj;

src/GRANULAR/fix_wall_gran.cpp

@@ -354,6 +354,15 @@ void FixWallGran::init()
   if (utils::strmatch(update->integrate_style,"^respa"))
     nlevels_respa = (dynamic_cast<Respa *>( update->integrate))->nlevels;
 
+  // check for compatible heat conduction atom style
+
+  if (heat_flag) {
+    if (!atom->temperature_flag)
+      error->all(FLERR,"Heat conduction in fix wall/gran requires atom style with temperature property");
+    if (!atom->heatflux_flag)
+      error->all(FLERR,"Heat conduction in fix wall/gran requires atom style with heatflux property");
+  }
+
   // check for FixRigid so can extract rigid body masses
 
   fix_rigid = nullptr;
@@ -529,9 +538,10 @@ void FixWallGran::post_force(int /*vflag*/)
     model->dx[1] = dy;
     model->dx[2] = dz;
     model->radi = radius[i];
+    model->radj = rwall;
     if (model->beyond_contact) model->touch = history_one[i][0];
 
-    touchflag = model->check_contact(rwall);
+    touchflag = model->check_contact();
 
     if (!touchflag) {
       if (use_history)

src/GRANULAR/fix_wall_gran_region.cpp

@@ -230,6 +230,8 @@ void FixWallGranRegion::post_force(int /*vflag*/)
       model->dx[1] = region->contact[ic].dely;
       model->dx[2] = region->contact[ic].delz;
       model->radi = radius[i];
+      model->radj = region->contact[ic].radius;
+      model->r = region->contact[ic].r;
       if (model->beyond_contact) model->touch = history_many[i][c2r[ic]][0];
 
       touchflag = model->check_contact(region->contact[ic].r);

src/GRANULAR/pair_granular.cpp

@@ -442,6 +442,13 @@ void PairGranular::init_style()
   if (comm->ghost_velocity == 0)
     error->all(FLERR,"Pair granular requires ghost atoms store velocity");
 
+  if (heat_flag) {
+    if (!atom->temperature_flag)
+      error->all(FLERR,"Heat conduction in pair granular requires atom style with temperature property");
+    if (!atom->heatflux_flag)
+      error->all(FLERR,"Heat conduction in pair granular requires atom style with heatflux property");
+  }
+
   // allocate history and initialize models
   int size_max[NSUBMODELS] = {0};