---

make_bc_rad.c

Go to the documentation of this file.

/*---------------------------------------------------------------------------*/
/* SCCS Information: @(#)make_bc.c      1.1    3/7/96 */
/*---------------------------------------------------------------------------*/
/* Implements the following boundary conditions
  
   Imposed pressure on x=0,1 and y=1 verifying the mass conservation and
   Bernouilli's equation
   Neuman conditions for velocity on x=0,1 and y=1
   axis of symmetry on y=0
 
*/   
#include "utilf.h"
#include "stypes.h"
#include "markers.h"
#include "extra.h"
#include "make_bc.h"
#include "make_bc_periodic.h"


static real rp_pressure(real pe, real pR,
                        real R, real Rp, real x, real y, real xc)
{
  real r = sqrt(sq(x - xc) + sq(y - 2.));
  /*  return pe*(1. - R/r) + 
    (R*pR + 0.5*R*sq(Rp)*(1. - cube(R/r)))/r;
    */
  return pe;
}


void bc_vector_bound(real2D u, real2D v, 
                     real R, real Rp, real xc,
                     int nx, int ny)
{
  int i, j;
  real r;

  R *= R;
  /* top and bottom walls */
  for (i = 1; i <= nx; i++) {
    u[i][1] = u[i][2];
    v[i][2] = 0.0;
    r = sqrt(sq((real)i - xc) + sq((real)ny + 0.5 - 2.));
    u[i][ny] = R/cube(r)*Rp*((real)i - xc);
    r = sqrt(sq((real)i + 0.5 - xc) + sq((real)ny - 2.));
    v[i][ny] = R/cube(r)*Rp*((real)ny - 2.);
  }
  /* left and right walls */
  for (j = 1; j <= ny; j++) {
    r = sqrt(sq(2. - xc) + sq((real)j + 0.5 - 2.));
    u[2][j] = R/cube(r)*Rp*(2. - xc);
    r = sqrt(sq((real)nx - xc) + sq((real)j + 0.5 - 2.));
    u[nx][j] = R/cube(r)*Rp*((real)nx - xc);
    r = sqrt(sq(1.5 - xc) + sq((real)j - 2.));
    v[1][j] = R/cube(r)*Rp*((real)j - 2.);
    r = sqrt(sq((real)nx + 0.5 - xc) + sq((real)j - 2.));
    v[nx][j] = R/cube(r)*Rp*((real)j - 2.);
  }
}


void bc_vector(real2D u, real2D a, real2D v, real2D c,
               interface in,
               real R, real Rp,
               int nx, int ny)
{
  int i, j;

  /* extrapolate the interfacial cells */
  for (i = 2; i < nx; i++)
    for (j = 2; j < ny; j++) {
      if (a[i-1][j] > 0.0 && a[i-1][j] < 1.0)
        u[i][j] = interplate_u(u, a, NULL,
                               (real)i, 0.5 + (real)j, nx, ny, INTERPLATE_FULL);
      if (i != 2 && a[i-1][j] == 0.0)
        u[i][j] = 0.0;
      if (c[i][j-1] > 0.0 && c[i][j-1] < 1.0)
        v[i][j] = interplate_v(v, c, NULL, 
                               0.5 + (real)i, (real)j, nx, ny, INTERPLATE_FULL);
      if (c[i][j-1] == 0.0)
        v[i][j] = 0.0;
    }

  bc_vector_bound(u, v, R, Rp, 0.5*(nx - 2.) + 2., nx, ny);
}


void bc_vector_div(real2D u, real2D a, real2D v, real2D c,
                   interface in,
                   real R, real Rp,
                   int nx, int ny)
{
  int i, j;
  real eu, ev, a11, a22, a12, a21;

  /* extrapolate the interfacial cells */
  for (i = 2; i < nx; i++)
    for (j = 2; j < ny; j++) {
      if (a[i-1][j] > 0.0 && a[i-1][j] < 1.0) {
        extra_velocity((real)i, 0.5 + (real)j, u, a, v, c, in, 
                       &eu, &ev, &a11, &a22, &a12, &a21, nx, ny);
        u[i][j] = eu;
      }
      if (i != 2 && a[i-1][j] == 0.0)
        u[i][j] = 0.0;
      if (c[i][j-1] > 0.0 && c[i][j-1] < 1.0) {
        extra_velocity(0.5+(real)i, (real)j, u, a, v, c, in, 
                       &eu, &ev, &a11, &a22, &a12, &a21, nx, ny);
        v[i][j] = ev;
      }
      if (c[i][j-1] == 0.0)
        v[i][j] = 0.0;
    }

  bc_vector_bound(u, v, R, Rp, 0.5*(nx - 2.) + 2., nx, ny);
}


void bc_pressure(real2D p, real2D cc, interface in,
                 real pe, real pg, real pR, real R, real Rp,
                 int nx, int ny)
{
  int i, j;

  /* extrapolate the interfacial cells */
  for (i = 2; i < nx; i++) {
    for (j = 2; j < ny; j++) {
      if (cc[i][j] > 0.0 && cc[i][j] < 1.0)
        p[i][j] = interplate_p(p, cc, in,
                               0.5 + (real)i, 0.5 + (real)j, nx, ny);
      if (cc[i][j] == 0.0)
        p[i][j] = pg;
    }
  }

  for (i = 1; i <= nx; i++)
    p[i][ny - 1] = rp_pressure(pe, pR, R, Rp, 
                               0.5 + (real)i, 0.5 + (real)(ny - 1), 
                               0.5*(nx - 2.) + 2.);
  for (j = 1; j <= ny; j++) {
    p[2][j] = rp_pressure(pe, pR, R, Rp, 
                          2.5, 0.5 + (real)j, 
                          0.5*(nx - 2.) + 2.);
    p[nx - 1][j] = rp_pressure(pe, pR, R, Rp, 
                               0.5 + (real)(nx - 1), 0.5 + (real)j,
                               0.5*(nx - 2.) + 2.);
  }
  bc_scalar(p, nx, ny, NULGRAD);
}


void bc_tensor(real2D S11, real2D S22, real2D S12, int nx, int ny)
{
  int i, j;
  
  /*  slip conditions on the walls S12 = 0  */
  for (i = 2; i <= nx - 1; i++) {
    S12[i][2] = 0.0;
    S12[i][ny] = S12[i][ny-1];
  }
  for (j = 2; j <= ny - 1; j++) {
    S12[2][j] = S12[3][j];
    S12[nx][j] = S12[nx-1][j];
  }
}


void bc_scalar(real2D scal, int nx, int ny, char sw)
{
  int i, j;
  
  /*  div: gradient equal to zero on y=0,1 planes  */
  switch(sw) {
    case NULGRAD:
      for (i = 1; i <= nx; i++)
        {
          scal[i][1] = scal[i][2];
          scal[i][ny] = scal[i][ny - 1];
        }
      for (j = 1; j <= ny; j++) {
        scal[1][j] = scal[2][j];
        scal[nx][j] = scal[nx-1][j];
      }
      break;
    case NULGRAD2:
      for (i = 1; i <= nx; i++)
        scal[i][1] = scal[i][2];
      break;
    case NUL:
      /*  residue equal to zero on the walls  */ 
      for (i = 1; i <= nx; i++)
        scal[i][1] = scal[i][ny] = 0.0;
      for (j = 1; j <= ny; j++)
        scal[1][j] = scal[nx][j] = 0.0;
      break;
    default:
      printf("illegal switch value for a scalar field\n");
      exit(1);
  }
}


int bc_xcoord(int i, int nx)
{
  /* mirror symmetric boundary conditions */
  if (i <= 1)
    return 3 - i;
  if (i >= nx)
    return 2*nx - i - 1;
  return i;
}


int bc_ycoord(int j, int ny)
{
  /* mirror symmetric boundary conditions */
  if (j <= 1)
    return 3 - j;
  if (j >= ny)
    return 2*ny - j - 1;
  return j;
}

---