Merge pull request #307 from ludwig-cf/feature-305-lubrication

Feature 305 lubrication

README.md

@@ -3,7 +3,7 @@
 
 A lattice Boltzmann code for complex fluids
 
-![Build Status](https://github.com/ludwig-cf/ludwig/actions/workflows/regression.yml/badge.svg?branch=develop)
+![Build Status](https://github.com/ludwig-cf/ludwig/actions/workflows/regression.yml/badge.svg)
 [![CII Best Practices](https://bestpractices.coreinfrastructure.org/projects/1998/badge)](https://bestpractices.coreinfrastructure.org/projects/1998)
 
 

src/bbl.c

@@ -49,6 +49,7 @@ struct bbl_s {
   int active;           /* Global flag for active particles. */
   int ellipsoid_didt;   /* Switch for unsteady term method */
   int ndist;            /* Number of LB distributions active */
+  double eta;           /* Dynamic viscosity for lubrication correction */
   double deltag;        /* Excess or deficit of phi between steps */
   double stress[3][3];  /* Surface stress diagnostic */
 };
@@ -77,6 +78,9 @@ void record_force_torque(colloid_t * pc);
 void bbl_ellipsoid_unsteady_mI(const double q[4], const double mi[3],
 			       const double omega[3], double didt[3][3]);
 
+int bbl_wall_lubr_correction_ellipsoid(bbl_t * bbl, wall_t * wall,
+				       colloid_t * pc, double wdrag[3]);
+
 /*****************************************************************************
  *
  *  bbl_create
@@ -103,6 +107,15 @@ int bbl_create(pe_t * pe, cs_t * cs, lb_t * lb, bbl_t ** pobj) {
   bbl->ellipsoid_didt = BBL_ELLIPSOID_UPDATE_QUATERNION;
   lb_ndist(lb, &bbl->ndist);
 
+  /* I would like to obtain the viscosity from the lb data structure;
+   * but it's not present at initialisation, so ... */
+
+  {
+    physics_t * phys = NULL;
+    physics_ref(&phys);
+    physics_eta_shear(phys, &bbl->eta);
+  }
+
   *pobj = bbl;
 
   return 0;
@@ -898,7 +911,6 @@ int bbl_update_colloids(bbl_t * bbl, wall_t * wall, colloids_info_t * cinfo) {
   /* All colloids, including halo */
 
   colloids_info_all_head(cinfo, &pc);
-
   for ( ; pc; pc = pc->nextall) {
 
     if (pc->s.bc != COLLOID_BC_BBL) continue;
@@ -1186,7 +1198,9 @@ void bbl_ladd_ellipsoid(bbl_t * bbl, colloid_t * pc, wall_t * wall,
 
   util_q4_inertia_tensor(pc->s.quat, mI_P, mI);
 
-  wall_lubr_sphere(wall, pc->s.ah, pc->s.r, dwall);
+
+  /* Lubrication correction at wall ... */
+  bbl_wall_lubr_correction_ellipsoid(bbl, wall, pc, dwall);
 
   /* Add inertial terms to diagonal elements */
 
@@ -1370,9 +1384,8 @@ int bbl_6x6_gaussian_elimination(double a[6][6], double xb[6]) {
 static int bbl_wall_lubrication_account(bbl_t * bbl, wall_t * wall,
 					colloids_info_t * cinfo) {
 
-  double f[3] = {0.0, 0.0, 0.0};
-  double dwall[3];
   colloid_t * pc = NULL;
+  double f[3]    = {0.0, 0.0, 0.0};
 
   assert(bbl);
   assert(cinfo);
@@ -1382,8 +1395,19 @@ static int bbl_wall_lubrication_account(bbl_t * bbl, wall_t * wall,
   colloids_info_local_head(cinfo, &pc);
 
   for (; pc; pc = pc->nextlocal) {
+    double dwall[3] = {0};
     if (pc->s.bc != COLLOID_BC_BBL) continue;
-    wall_lubr_sphere(wall, pc->s.ah, pc->s.r, dwall);
+
+    if (pc->s.shape == COLLOID_SHAPE_SPHERE) {
+      wall_lubr_sphere(wall, pc->s.ah, pc->s.r, dwall);
+    }
+    else if (pc->s.shape == COLLOID_SHAPE_ELLIPSOID) {
+      bbl_wall_lubr_correction_ellipsoid(bbl, wall, pc, dwall);
+    }
+    else {
+      /* Specifically, no COLLOID_SHAPE_DISK */
+      assert(0);
+    }
     f[X] -= pc->s.v[X]*dwall[X];
     f[Y] -= pc->s.v[Y]*dwall[Y];
     f[Z] -= pc->s.v[Z]*dwall[Z];
@@ -1496,3 +1520,117 @@ void bbl_ellipsoid_unsteady_mI(const double q[4],
 
   return;
 }
+
+/*****************************************************************************
+ *
+ *  bbl_wall_lubr_correction_ellipsoid
+ *
+ *  We will compute the distance between the centre of the ellipsoid and
+ *  the tangent plane normal to each wall. This will provide the
+ *  normal separation. One computation per co-ordinate direction is
+ *  required.
+ *
+ *  This is fed into a fudge for the lubrication correction based on the
+ *  analytical expression for a sphere with a hydrodynamic radius of the
+ *  semi-major axis. A more representative effective radius could be
+ *  identified with some geometry.
+ *
+ *  A velocity-independent drag is returned.
+ *
+ *  Returns zero on success. If an overlap between the surface of the
+ *  ellipsoid and a wall position is detected, a non-zero value will
+ *  be returned.
+ *
+ *****************************************************************************/
+
+int bbl_wall_lubr_correction_ellipsoid(bbl_t * bbl, wall_t * wall,
+				       colloid_t * pc, double wdrag[3]) {
+  int ifail = 0;
+  double lmin[3] = {0};
+  double ltot[3] = {0};
+  double ah      = pc->s.elabc[X];  /* semi-major axis length */
+
+  cs_lmin(bbl->cs, lmin);
+  cs_ltot(bbl->cs, ltot);
+
+  if (wall->param->isboundary[X]) {
+    double xtan    = -1.0;     /* colloid centre to colloid tangent plane */
+    double nhat[3] = {1.0, 0.0, 0.0};                  /* +/- wall normal */
+    double hc      = wall->param->lubr_rc[X];                  /* cut off */
+    double dh      = wall->param->lubr_dh[X];    /* offset distance  wall */
+
+    double xc      = pc->s.r[X];
+    double drag    = 0.0;
+
+    xtan = util_q4_distance_to_tangent_plane(pc->s.elabc, pc->s.quat, nhat);
+
+    /* Compute a correction at (potentially) each end. */
+    {
+      double hbot = xc - (lmin[X] + dh) - xtan;    /* surface to bottom wall */
+      drag += wall_lubr_drag(bbl->eta, ah, hbot, hc);
+      if (hbot < 0.0) ifail = -1;
+    }
+    {
+      double htop = lmin[X] + (ltot[X] - dh) - xc - xtan;  /* surface to top */
+      drag += wall_lubr_drag(bbl->eta, ah, htop, hc);
+      if (htop < 0.0) ifail = +1;
+    }
+    wdrag[X] = drag;
+  }
+
+  /* Repeat for y-direction */
+
+  if (wall->param->isboundary[Y]) {
+    double ytan    = -1.0;     /* colloid centre to colloid tangent plane */
+    double nhat[3] = {0.0, 1.0, 0.0};                  /* +/- wall normal */
+    double hc      = wall->param->lubr_rc[Y];                  /* cut off */
+    double dh      = wall->param->lubr_dh[Y];    /* offset distance  wall */
+
+    double yc      = pc->s.r[Y];
+    double drag    = 0.0;
+
+    ytan = util_q4_distance_to_tangent_plane(pc->s.elabc, pc->s.quat, nhat);
+
+    /* Compute a correction at (potentially) each end. */
+    {
+      double hbot = yc - (lmin[Y] + dh) - ytan;    /* surface to bottom wall */
+      drag += wall_lubr_drag(bbl->eta, ah, hbot, hc);
+      if (hbot < 0.0) ifail = -1;
+    }
+    {
+      double htop = lmin[Y] + (ltot[Y] - dh) - yc - ytan;  /* surface to top */
+      drag += wall_lubr_drag(bbl->eta, ah, htop, hc);
+      if (htop < 0.0) ifail = +1;
+    }
+    wdrag[Y] = drag;
+  }
+
+  /* Repeat for z-direction */
+
+  if (wall->param->isboundary[Z]) {
+    double ztan    = -1.0;     /* colloid centre to colloid tangent plane */
+    double nhat[3] = {0.0, 0.0, 1.0};                  /* +/- wall normal */
+    double hc      = wall->param->lubr_rc[Z];                  /* cut off */
+    double dh      = wall->param->lubr_dh[Z];    /* offset distance  wall */
+
+    double zc      = pc->s.r[Z];
+    double drag    = 0.0;
+
+    ztan = util_q4_distance_to_tangent_plane(pc->s.elabc, pc->s.quat, nhat);
+
+    /* Compute a correction at (potentially) each end. */
+    {
+      double hbot = zc - (lmin[Z] + dh) - ztan;    /* surface to bottom wall */
+      drag += wall_lubr_drag(bbl->eta, ah, hbot, hc);
+      if (hbot < 0.0) ifail = -1;
+    }
+    {
+      double htop = lmin[Z] + (ltot[Z] - dh) - zc - ztan;  /* surface to top */
+      drag += wall_lubr_drag(bbl->eta, ah, htop, hc);
+      if (htop < 0.0) ifail = +1;
+    }
+    wdrag[Z] = drag;
+  }
+
+  return ifail;
+}

src/colloids.h

@@ -7,7 +7,7 @@
  *  Edinburgh Soft Matter and Statistical Physics Group and
  *  Edinburgh Parallel Computing Centre
  *
- *  (c) 2010-2023 The University of Edinburgh
+ *  (c) 2010-2024 The University of Edinburgh
  *
  *  Contributing authors:
  *  Kevin Stratford ([email protected])
@@ -184,7 +184,6 @@ __host__ int colloid_rb(colloids_info_t * info, colloid_t * pc, int index,
 __host__ int colloid_rb_ub(colloids_info_t * info, colloid_t * pc, int index,
 			   double rb[3], double ub[3]);
 
-__host__ int is_site_inside_colloid(colloid_t * pc, double rsep[3]);
 __host__ int colloids_type_check(colloids_info_t * info);
 __host__ int colloids_ellipsoid_abc_check(colloids_info_t * info);
 

src/util_ellipsoid.c

@@ -589,7 +589,7 @@ int util_ellipsoid_euler_from_vectors(const double a0[3], const double b0[3],
  *  cf2 = (16/3) e^3 [+2e +(3e^2 - 1) log(1+e/1-e)]^-1
  *
  *  The Stokes' relationship is:
- *  F = 6 pi mu a (U_1 cf1 xhat + U_2 cf2 yhat)
+ *  F = 6 pi eta a (U_1 cf1 xhat + U_2 cf2 yhat)
  *
  *  a    is the semi-major axis
  *  b    is the semi-minor axis (b < a)
@@ -634,14 +634,15 @@ int util_ellipsoid_prolate_settling_velocity(double a,
  *  util_discrete_volume_ellipsoid
  *
  *  A utility to return the discrete volume of an ellipsoid with axes
- *  abs[3] and orientation described by quaternion q placed on the
+ *  abc[3] and orientation described by quaternion q placed on the
  *  unit lattice at position r.
  *
  *  We just draw a large box around the central position based on the
  *  semi-principal radius.
  *  The inside/outside determination is via util_q4_is_inside_ellipsoid().
  *
  *  The return value is the volume as an integer.
+ *  The volume is also returned as a double in the argument list.
  *
  *****************************************************************************/
 
@@ -670,3 +671,113 @@ int util_discrete_volume_ellipsoid(const double abc[3], const double r0[3],
 
   return inside;
 }
+
+/****************************************************************************
+ *
+ *  util_ellipsoid_eqn_lhs
+ *
+ *  Calculate the lhs of the equation for an ellipsoid using the equation
+ *  of the ellipsoid
+ *
+ ****************************************************************************/
+
+double ellipsoid_eqn_lhs(const double elA[3][3], const double rsep[3]) {
+
+  double lhs =
+    + elA[0][0]*rsep[X]*rsep[X]
+    + elA[1][1]*rsep[Y]*rsep[Y]
+    + elA[2][2]*rsep[Z]*rsep[Z]
+    + (elA[0][1]+elA[1][0])*rsep[X]*rsep[Y]
+    + (elA[0][2]+elA[2][0])*rsep[X]*rsep[Z]
+    + (elA[1][2]+elA[2][1])*rsep[Y]*rsep[Z];
+
+  return lhs;
+}
+
+/*****************************************************************************
+ *
+ *  util_ellipsoid_gaussian_rad_curv
+ *
+ *  Calculate the Gaussian curvature of an ellipsoid in the body fixed
+ *  coordinate system
+ *
+ ****************************************************************************/
+
+double Gaussian_RadCurv_ellipsoid(const double *elabc, const double x[3]) {
+
+  double a2 = elabc[0]*elabc[0];
+  double b2 = elabc[1]*elabc[1];
+  double c2 = elabc[2]*elabc[2];
+  double s1= ((x[0]*x[0])/(a2*a2)) + ((x[1]*x[1])/(b2*b2)) + ((x[2]*x[2])/(c2*c2));
+
+  return elabc[0]*elabc[1]*elabc[2]*s1;
+}
+
+/*****************************************************************************
+ *
+ *  util_q4_to_r
+ *
+ *  Return the rotation matrix which corresponds to q.
+ *
+ *  The rotation matrix may be constructed directly from
+ *
+ *     b = (2q_0 - 1) a + 2(q.a)q + 2q_0 q x a
+ *
+ *  See util_q4_rotate_vector() above.
+ *
+ *****************************************************************************/
+
+void util_q4_to_r(const double q[4], double r[3][3]) {
+
+  r[0][0] = 2.0*(q[0]*q[0] - 0.5 + q[1]*q[1]            );
+  r[0][1] = 2.0*(                  q[2]*q[1] - q[0]*q[3]);
+  r[0][2] = 2.0*(                  q[3]*q[1] + q[0]*q[2]);
+
+  r[1][0] = 2.0*(                  q[1]*q[2] + q[0]*q[3]);
+  r[1][1] = 2.0*(q[0]*q[0] - 0.5 + q[2]*q[2]            );
+  r[1][2] = 2.0*(                  q[3]*q[2] - q[0]*q[1]);
+
+  r[2][0] = 2.0*(                  q[1]*q[3] - q[0]*q[2]);
+  r[2][1] = 2.0*(                  q[2]*q[3] + q[0]*q[1]);
+  r[2][2] = 2.0*(q[0]*q[0] - 0.5 + q[3]*q[3]            );
+
+  return;
+}
+
+/*****************************************************************************
+ *
+ *  util_q4_distance_to_tangent_plane
+ *
+ *  For ellipsoid with principal semi-axis lengths abc and orientation
+ *  described by q, what is the normal distance from the centre to the
+ *  tangent plane with plane normal nhat (nhat in the lab frame).
+ *
+ *  See e.g., Ferguson, Mathematical Geology, 11, 329--336 (1979).
+ *  https://link.springer.com/content/pdf/10.1007/BF01034997.pdf
+ *
+ *  In brief:
+ *  If the normal is n_i then beta_j = n_i R_ij where R_ij is the
+ *  rotation to move from the lab frame to the principal frame.
+ *
+ *  The distance is then d = sqrt(beta_i^2 abc_i^2)
+ *
+ *  Here, the rotation matrix is replaced by the rotation q^*.
+ *
+ *****************************************************************************/
+
+double util_q4_distance_to_tangent_plane(const double abc[3],
+					 const double q[4],
+					 const double nhat[3]) {
+  double d       = 0.0;
+  double qbar[4] = {q[0], -q[1], -q[2], -q[3]};
+  double beta[3] = {0};
+
+  util_q4_rotate_vector(qbar, nhat, beta);
+
+  for (int i = 0; i < 3; i++) {
+    d += beta[i]*beta[i]*abc[i]*abc[i];
+  }
+  assert(d >= 0.0);
+
+  return sqrt(d);
+}

src/util_ellipsoid.h

@@ -27,6 +27,7 @@ int  util_q4_is_inside_ellipsoid(const double q[4], const double elabc[3],
 				 const double r[3]);
 void util_q4_inertia_tensor(const double q[4], const double moment[3],
 			    double mI[3][3]);
+
 int util_ellipsoid_euler_from_vectors(const double a0[3], const double b0[3],
 				      double euler[3]);
 int util_ellipsoid_prolate_settling_velocity(double a, double b, double eta,
@@ -38,6 +39,11 @@ void matrix_product(const double a[3][3], const double b[3][3],
 		    double result[3][3]);
 void matrix_transpose(const double a[3][3], double result[3][3]);
 
+void util_q4_to_r(const double q[4], double r[3][3]);
+double util_q4_distance_to_tangent_plane(const double abc[3],
+					 const double q[4],
+					 const double nhat[3]);
+
 /*****************************************************************************
  *
  *  __host__ __device__ static inline