---

tension.c

Go to the documentation of this file.

#include "utilf.h"
#include "markers.h"
#include "surfaces.h"
#include "tension.h"
#include "extra.h"
#include "inout.h"


#define  PHI_R_LINEAR(in, i, t1) \
         (((t1 - (in).t[i])*(in).p[i+1] + ((in).t[i+1] - t1)*(in).p[i])/ \
          ((in).t[i+1] - (in).t[i]))


/* Compute r\Phi given two polynomials describing the interface and the
   parameter t */
static real phi_r(polynom3 px, polynom3 py, real t, real sigma, real mu,
                 real2D u, real2D v, real2D ap, int nx, int ny)
{
  real n1, n2, xp, yp, dl, x, y, kappa;
  real dudx, dvdy, dudy, dvdx, eu, ev;

  xp = splint1(px, t);
  yp = splint1(py, t);
  dl = sqrt(sq(xp) + sq(yp));
  n1 = - yp/dl; n2 = xp/dl;
  y = splint(py, t) - 2.0;
  kappa = y*(yp*splint2(px, t) - xp*splint2(py, t))/sq(dl) + xp/dl;
/*cld  kappa = y*(yp*splint2(px, t) - xp*splint2(py, t))/(sq(dl)*dl) + xp/dl; cld*/

  x = splint(px, t); y = splint(py, t);
#ifdef EXTRA_LINEAR
  extra_velocity_normals(x, y, n1, n2, x, y, u, v, ap, 
                         &eu, &ev, 
                         &dudx, &dvdy, &dudy, &dvdx, nx, ny);
#else
  extra_velocity_normals2(x, y, n1, n2, x, y, u, v, ap, 
                          &eu, &ev, 
                          &dudx, &dvdy, &dudy, &dvdx, nx, ny);
#endif
  y -= 2.;

  return - sigma*kappa 
    + 2.*y*mu*(dudx*sq(n1) + dvdy*sq(n2) + (dudy + dvdx)*n1*n2);
}


/* Compute \sum_{t1}^{t2} \Phi r dr */
static real phi_r_dr(interface in, int j, real t1, real t2)
{
  /* coefficients of Newton-Cotes formula */
  static real coef[7] = {1., 3., 3., 2., 3., 3., 1.}; 
  int i;
  real yp, h = (t2 - t1)/6., t, sum = 0.0;
  real a = (in.p[j+1] - in.p[j])/(in.t[j+1] - in.t[j]);
  real b = (in.t[j+1]*in.p[j] - in.t[j]*in.p[j+1])/(in.t[j+1] - in.t[j]);
  polynom3 spy = in.spy[j];
  
  for (i = 0; i < 7; i++) {
    t = t1 + h*i;
    yp = splint1(spy, t);
    sum += coef[i]*yp*(a*t + b);
  }
  return 3.*h*sum/8.;
}


/* Compute \sum_{t1}^{t2} \Phi r dz */
static real phi_r_dz(interface in, int j, real t1, real t2)
{
  /* coefficients of Newton-Cotes formula */
  static real coef[7] = {1., 3., 3., 2., 3., 3., 1.}; 
  int i;
  real xp, h = (t2 - t1)/6., t, sum = 0.0;
  real a = (in.p[j+1] - in.p[j])/(in.t[j+1] - in.t[j]);
  real b = (in.t[j+1]*in.p[j] - in.t[j]*in.p[j+1])/(in.t[j+1] - in.t[j]);
  polynom3 spx = in.spx[j];
  
  for (i = 0; i < 7; i++) {
    t = t1 + h*i;
    xp = splint1(spx, t);
    sum += coef[i]*xp*(a*t + b);
  }
  return 3.*h*sum/8.;
}


void interface_tensionu(interface in, interface_cut cut,
                        real2D u, real2D v, 
                        real2D cu, real2D ap, real2D syu, real2D p,
                        real2D pi, real2D ni,
                        int nx, int ny)
{
  int i, j, l;
  real t, t1, t2;
  polynom3 px, py;

  fill0(pi, nx, ny);
  fill0(ni, nx, ny);

  /* For each intersection */
  for (l = 0; l < cut.ncut - 1; l++) {
    px = in.spx[cut.list[l].ip];
    py = in.spy[cut.list[l].ip];

    t1 = cut.list[l].t; t2 = cut.list[l+1].t;
    t = 0.5*(t1 + t2);
    i = splint(px, t) - cut.xo;
    j = splint(py, t) - cut.yo;

    /* Intersection with the verticals */
    if (cut.list[l].type == VER_TYPE) {
      if (cut.list[l].coord == i) {
        if (!INP(i-1,j)) {
          pi[i-1][j] += PHI_R_LINEAR(in, cut.list[l].ip, t1)/
            (splint(py, t1) - 2.);
          ni[i-1][j] += 1.0;
        }
      }
      else if (!INP(i,j)) {
        pi[i][j] += PHI_R_LINEAR(in, cut.list[l].ip, t1)/
          (splint(py, t1) - 2.);
        ni[i][j] += 1.0;
      }
    }
    
    if (INU(i,j)) u[i][j] -= phi_r_dr(in, cut.list[l].ip, t1, t2)/cu[i][j];
  }
  for (i = 2; i <= nx; i++)
    for (j = 2; j < ny; j++) 
      if (INU(i,j)) {
        if (pi[i-1][j] != 0.0)
          u[i][j] += syu[i][j]*pi[i-1][j]/(cu[i][j]*ni[i-1][j]);
        if (pi[i][j] != 0.0)
          u[i][j] -= syu[i+1][j]*pi[i][j]/(cu[i][j]*ni[i][j]);
      }
}


void interface_tensionv(interface in, interface_cut cut,
                        real2D u, real2D v, real2D cv, real2D ap, real2D sxv, 
                        real2D av, real2D p,
                        real2D pi, real2D ni,
                        int nx, int ny)
{
  int i, j, l;
  real t, t1, t2;
  polynom3 px, py;

  fill0(pi, nx, ny);
  fill0(ni, nx, ny);

  /* For each intersection */
  for (l = 0; l < cut.ncut - 1; l++) {
    px = in.spx[cut.list[l].ip];
    py = in.spy[cut.list[l].ip];

    t1 = cut.list[l].t; t2 = cut.list[l+1].t;
    t = 0.5*(t1 + t2);
    i = splint(px, t) - cut.xo;
    j = splint(py, t) - cut.yo;

    /* Intersection with the horizontals */
    if (cut.list[l].type == HOR_TYPE) {
      if (cut.list[l].coord == j) {
        if (!INP(i,j-1)) {
          pi[i][j-1] += PHI_R_LINEAR(in, cut.list[l].ip, t1)
            /(splint(py, t1) - 2.);
          ni[i][j-1] += 1.0;
        }
      }
      else if (!INP(i,j)) {
        pi[i][j] += PHI_R_LINEAR(in, cut.list[l].ip, t1)
          /(splint(py, t1) - 2.);
        ni[i][j] += 1.0;
      }
    }

    if (INV(i,j)) v[i][j] += phi_r_dz(in, cut.list[l].ip, t1, t2)/cv[i][j];
  }
  for (i = 2; i < nx; i++)
    for (j = 3; j <= ny; j++) 
      if (INV(i,j)) {
        if (pi[i][j-1] != 0.0)
          v[i][j] += (sxv[i][j] + av[i][j]/2.)*pi[i][j-1]
            /(cv[i][j]*ni[i][j-1]);
        else if (!INP(i,j-1))
          v[i][j] += av[i][j]/2.*p[i][j]/cv[i][j];
        if (pi[i][j] != 0.0)
          v[i][j] -= (sxv[i][j+1] - av[i][j]/2.)*pi[i][j]
            /(cv[i][j]*ni[i][j]);
        else if (!INP(i,j))
          v[i][j] += av[i][j]/2.*p[i][j-1]/cv[i][j];
      }
}


void interfacial_pressure(interface in, real2D u, real2D v, real2D ap,
                          real sigma, real mu, real pg, real tau,
                          int nx, int ny)
{
  real h = 1.0/(nx - 2.);
  int i;

  sigma *= sq(tau)/(h*h*h);
  mu *= tau/sq(h);
  pg *= sq(tau/h);

  for (i = 0; i < in.n - 1; i++)
    in.p[i] = pg*(in.y[i] - 2.) + 
      phi_r(in.spx[i], in.spy[i], in.t[i], sigma, mu, 
                         u, v, ap, nx, ny);
  in.p[in.n-1] = pg*(in.y[in.n-1] - 2.) +
    phi_r(in.spx[in.n-2], in.spy[in.n-2], in.t[in.n-1], 
          sigma, mu, u, v, ap, nx, ny);
}

---