myviscosity-viscoplastic.h

Regularization method for a yield stress fluid
- Initialization
- Regularization

Regularization method for a yield stress fluid

We define the yield stress T and the regularization parameters.

(const) face vector T = zerof;
double eps = 1.e-4;

When dealing with multiphase flows, we can define T using an arithmetic average. The user can overload this definition.

#ifndef T
#define T(f,T1,T2) (clamp(f,0.,1.)*(T1 - T2) + T2)
#endif

We also define a scalar quantity to characterize the yielded and unyielded regions of the domain: 0 for yielded, 1 for unyielded. Note that when using a regularization method, this quantity is post-processed and is not computed directly.

scalar yuyz[];

static inline void refine_yuyz (Point point, scalar s)
{
  refine_bilinear (point, s);
  foreach_child()
    s[] = s[] > 0. ? 1. : 0.;
}

static inline void restriction_yuyz (Point point, scalar s)
{
  restriction_average (point, s);
  s[] = s[] > 0. ? 1. : 0.;
}

Initialization

event defaults (i = 0)
{
#if TREE
  yuyz.refine = yuyz.prolongation = refine_yuyz;
  yuyz.restriction = restriction_yuyz;
#endif // TREE
}

event init (i = 0)
{
  foreach()
    yuyz[] = 0.; // 0 : yielded 
}

Regularization

Helpful functions

To compute the second invariant of the stress tensor S_2 on each face of the grid, we define helpfull functions to compute gradients in the directions orthogonal to the face normal.

#if EMBED
foreach_dimension()
double yield_face_avg_gradient_t1_x (Point point, scalar a, int i) {
  double left = nodata, right = nodata;
  if (cs[i,0] && (emerged || csm1[i,0]))
    right = (fs.y[i,0] && fs.y[i,1] &&
	     cs[i,1] && cs[i,-1] &&
	     (emerged || (csm1[i,1] && csm1[i,-1])) ?  (a[i,1] - a[i,-1])/(2.*Delta) :
	     fs.y[i,1] && fs.y[i,2] &&
	     cs[i,1] && cs[i,2] &&
	     (emerged || (csm1[i,1] && csm1[i,2])) ?   (-a[i,2] + 4.*a[i,1] - 3.*a[i,0])/(2.*Delta) :
	     fs.y[i,0] && fs.y[i,-1] &&
	     cs[i,-1] && cs[i,-2] &&
	     (emerged || (csm1[i,-1] && csm1[i,-2])) ? (a[i,-2] - 4.*a[i,-1] + 3.*a[i,0])/(2.*Delta) :
	     fs.y[i,1] && cs[i,1] &&
	     (emerged || csm1[i,1]) ?                  (a[i,1] - a[i])/Delta :
	     fs.y[i,0] && cs[i,-1] &&
	     (emerged || csm1[i,-1]) ?                 (a[i] - a[i,-1])/Delta : nodata);
  if (cs[i-1,0] && (emerged || csm1[i-1,0]))
    left = (fs.y[i-1,0] && fs.y[i-1,1] &&
	     cs[i-1,1] && cs[i-1,-1] &&
	     (emerged || (csm1[i-1,1] && csm1[i-1,-1])) ?  (a[i-1,1] - a[i-1,-1])/(2.*Delta) :
	     fs.y[i-1,1] && fs.y[i-1,2] &&
	     cs[i-1,1] && cs[i-1,2] &&
	     (emerged || (csm1[i-1,1] && csm1[i-1,2])) ?   (-a[i-1,2] + 4.*a[i-1,1] - 3.*a[i-1,0])/(2.*Delta) :
	     fs.y[i-1,0] && fs.y[i-1,-1] &&
	     cs[i-1,-1] && cs[i-1,-2] &&
	     (emerged || (csm1[i-1,-1] && csm1[i-1,-2])) ? (a[i-1,-2] - 4.*a[i-1,-1] + 3.*a[i-1,0])/(2.*Delta) :
	     fs.y[i-1,1] && cs[i-1,1] &&
	     (emerged || csm1[i-1,1]) ?                    (a[i-1,1] - a[i-1])/Delta :
	     fs.y[i-1,0] && cs[i-1,-1] &&
	     (emerged || csm1[i-1,-1]) ?                   (a[i-1] - a[i-1,-1])/Delta : nodata);  
  return (right == nodata && left == nodata ? 0. :
	  right == nodata ? left :
	  left == nodata ? right :
	  fs.x[i] ? (left + right)/2. : 0.);
}

foreach_dimension()
double yield_face_avg_gradient_t2_x (Point point, scalar a, int i) {
  double left = nodata, right = nodata;
  if (cs[i] && (emerged || csm1[i]))
    right = (fs.z[i,0,0] && fs.z[i,0,1] &&
	     cs[i,0,1] && cs[i,0,-1] &&
	     (emerged || (csm1[i,0,1] && csm1[i,0,-1])) ?  (a[i,0,1] - a[i,0,-1])/(2.*Delta) :
	     fs.z[i,0,1] && fs.z[i,0,2] &&
	     cs[i,0,1] && cs[i,0,2] &&
	     (emerged || (csm1[i,0,1] && csm1[i,0,2])) ?   (-a[i,0,2] + 4.*a[i,0,1] - 3.*a[i,0,0])/(2.*Delta) :
	     fs.z[i,0,0] && fs.z[i,0,-1] &&
	     cs[i,0,-1] && cs[i,0,-2] &&
	     (emerged || (csm1[i,0,-1] && csm1[i,0,-2])) ? (a[i,0,-2] - 4.*a[i,0,-1] + 3.*a[i,0,0])/(2.*Delta) :
	     fs.z[i,0,1] && cs[i,0,1] &&
	     (emerged || csm1[i,0,1]) ?                    (a[i,0,1] - a[i])/Delta :
	     fs.z[i,0,0] && cs[i,0,-1] &&
	     (emerged || csm1[i,0,-1]) ?                   (a[i] - a[i,0,-1])/Delta : nodata);
  if (cs[i-1] && (emerged || csm1[i-1]))
    left = (fs.z[i-1,0,0] && fs.z[i-1,0,1] &&
	     cs[i-1,0,1] && cs[i-1,0,-1] &&
	     (emerged || (csm1[i-1,0,1] && csm1[i-1,0,-1])) ?  (a[i-1,0,1] - a[i-1,0,-1])/(2.*Delta) :
	     fs.z[i-1,0,1] && fs.z[i-1,0,2] &&
	     cs[i-1,0,1] && cs[i-1,0,2] &&
	     (emerged || (csm1[i-1,0,1] && csm1[i-1,0,2])) ?   (-a[i-1,0,2] + 4.*a[i-1,0,1] - 3.*a[i-1,0,0])/(2.*Delta) :
	     fs.z[i-1,0,0] && fs.z[i-1,0,-1] &&
	     cs[i-1,0,-1] && cs[i-1,0,-2] &&
	     (emerged || (csm1[i-1,0,-1] && csm1[i-1,0,-2])) ? (a[i-1,0,-2] - 4.*a[i-1,0,-1] + 3.*a[i-1,0,0])/(2.*Delta) :
	     fs.z[i-1,0,1] && cs[i-1,0,1] &&
	     (emerged || csm1[i-1,0,1]) ?                      (a[i-1,0,1] - a[i-1])/Delta :
	     fs.z[i-1,0,0] && cs[i-1,0,-1] &&
	     (emerged || csm1[i-1,0,-1]) ?                     (a[i-1] - a[i-1,0,-1])/Delta : nodata);  
  return (right == nodata && left == nodata ? 0. :
	  right == nodata ? left :
	  left == nodata ? right :
	  fs.x[i] ? (left + right)/2. : 0.);
}
#else
foreach_dimension()
double yield_face_avg_gradient_t1_x (Point point, scalar a, int i) { // avg dady, dy = 2*Delta
  double right = (a[i,1]   - a[i,-1])/(2.*Delta);
  double left  = (a[i-1,1] - a[i-1,-1])/(2.*Delta);
  return (left + right)/2.;
}

foreach_dimension()
double yield_face_avg_gradient_t2_x (Point point, scalar a, int i) { // avg dadz, dz = 2*Delta
  double right = (a[i,0,1]   - a[i,0,-1])/(2.*Delta);
  double left  = (a[i-1,0,1] - a[i-1,0,-1])/(2.*Delta);
  return (left + right)/2.;
}
#endif // EMBED

The following function updates the variable yuyz that characterizes yielded and unyielded regions.

void yielded_region ()
{

We first compute the second invariant of the stress tensor S_2.

  face vector S2[];
  foreach_face() {
    double D_11 = face_gradient_x (u.x, 0); // uxx
    double D_22 = yield_face_avg_gradient_t1_x (point, u.y, 0); // avg of uyy
    double D_12 = 0.5*(face_gradient_x (u.y,0) +
		       yield_face_avg_gradient_t1_x (point, u.x, 0)); // 0.5*(uyx + avg uxy)
    S2.x[] = (sq (D_11) + sq (D_22) + 2.*sq (D_12));
#if dimension == 3 
    double D_33 = yield_face_avg_gradient_t2_x (u.z, 0); // avg uzz
    double D_13 = 0.5*(face_gradient_x (u.z,0) +
		       yield_face_avg_gradient_t2_x (point, u.x, 0)); // 0.5*(uzx + avg uxz)
    double D_23 = 0.5*(yield_face_avg_gradient_t1_x (point, u.z, 0) +
		       yield_face_avg_gradient_t2_x (point, u.y, 0)); // 0.5*(avg uzy + avg uyz)
    S2.x[] += (sq (D_33) + 2.*sq (D_13) + 2.*sq (D_23));
#endif // dimension == 3
  }
#if AXI
  foreach_face(x)
    S2.x[] += sq ((u.y[] + u.y[-1])/(2.*y));
  foreach_face(y)
    S2.y[] += sq ((u.y[] + u.y[0,-1])/(2.*y + 1.e-20)); // y=r, so avoid r=0
#endif //AXI

#if EMBED

As S_2 is computed on each face of the domain, we do not account for the stress tensor on the embedded boundary.

#endif // EMBED

  foreach_face()
    S2.x[] = mu.x[]*(2.*sqrt (S2.x[]/2.));

  boundary ((scalar *) {S2});

Finally, we use S2 to determine if a cell is yielded or unyielded.

  foreach() {

By default, the cell is unyielded.

    yuyz[] = 1;

We compute the surface-weighted average yield stress and invariant stress tensor in the cell.

    double Ta = 0., Sa = 0., fa = 0.;
    foreach_dimension() {
      Ta += T.x[]  + T.x[1];
      Sa += S2.x[] + S2.x[1];
      fa += fm.x[] + fm.x[1];
    }
    Ta /= (fa + 1.e-20);
    Sa /= (fa + 1.e-20);
  
    if (cm[] <= 0. || (Ta && Sa <= Ta))
      yuyz[] = 0;
  }
  boundary ({yuyz});
}

Modification of the viscosity

event properties (i++)
{

We compute here the square of second invariant D2 (using the Frobenius norm) of the strain rate tensor (or deformation tensor). Note that D2 is not multiplied by the metric.

  face vector D2[];
  foreach_face() {
    double D_11 = face_gradient_x (u.x, 0); // uxx
    double D_22 = yield_face_avg_gradient_t1_x (point, u.y, 0); // avg of uyy
    double D_12 = 0.5*(face_gradient_x (u.y,0) +
		       yield_face_avg_gradient_t1_x (point, u.x, 0)); // 0.5*(uyx + avg uxy)
    D2.x[] = (sq (D_11) + sq (D_22) + 2.*sq (D_12));
#if dimension == 3 
    double D_33 = yield_face_avg_gradient_t2_x (u.z, 0); // avg uzz
    double D_13 = 0.5*(face_gradient_x (u.z,0) +
		       yield_face_avg_gradient_t2_x (point, u.x, 0)); // 0.5*(uzx + avg uxz)
    double D_23 = 0.5*(yield_face_avg_gradient_t1_x (point, u.z, 0) +
		       yield_face_avg_gradient_t2_x (point, u.y, 0)); // 0.5*(avg uzy + avg uyz)
    D2.x[] += (sq (D_33) + 2.*sq (D_13) + 2.*sq (D_23));
#endif // dimension == 3
  }

When using cylindrical coordinates, we add the term u_{\theta\theta} = u_r/r.

#if AXI
  foreach_face(x)
    D2.x[] += sq ((u.y[] + u.y[-1])/(2.*y));
  foreach_face(y)
    D2.y[] += sq ((u.y[] + u.y[0,-1])/(2.*y + 1.e-20)); // y=r, so avoid r=0
#endif //AXI

#if EMBED

As D_2 is computed on each face of the domain, we do not account for the stress tensor on the embedded boundary.

#endif // EMBED  

  boundary ((scalar *) {D2});

We now modify the viscosity to account for yielded regions using a regularization method.

  face vector muv = mu;

  foreach_face() {

We introduce here a regularized viscosity, depending on the value of ||D_2||=\sqrt{D_2/2}. Note here that the yield stress T is multiplied by the metric fm.

    double DD = sqrt (D2.x[]/2.);

    muv.x[] += T.x[]*min (1./eps, 1./(2.*DD + 1.e-20)*(1. - exp (-2.*DD/eps)));
  }
  boundary ((scalar *) {muv});
}

Finally, we determine if a cell is yielded or unyielded.

event end_timestep (i++)
{
  yielded_region ();
}