sandbox/yonghui/biomeca/2drrc.c

    An poiseuille flow With output resistance

    Here, we study the effect of the output resistance in the axisymmetric tube and compare our results with the theoretical values. Mesh adaptation is used by default.

    Here P_ {out} is similar to a paraller RC circuit (Details in KVmodel ) instead of 0, (R_1, R_2, C) are resistors and capacitors ,and Q is the value of the mass flow. This method is commonly used to simulate the effects of external vascular models.

    #include "axi.h"
#include "navier-stokes/centered.h"
#include "view.h"
#define ADAPT 1 
#define alpha_w 10. 
#define MIN_LEVEL 3
#define LEVEL 5
#define MAX_LEVEL 6
#define tmax 1*2.*M_PI
#define  ORR 0.1 //R1
#define  ORC 0.1 //R2
#define  OCC 10  //C

FILE * fp1; FILE * fp2; FILE * fp3;
double	Pout,Pold;
//double	Pold;
double	DEBIT, Qold;
//double	Qold;
double	PPPP;

    main fonction

    int main(){
  N = 2 << LEVEL;
  fp1 = fopen ("testpoint", "w");
  fp2 = fopen ("pressure.dat", "w");
  run();
}

    we set the no slip Boundary conditions for lateral wall

    u.n[top]  = dirichlet(0.);
u.t[top]  = dirichlet(0.); 
//We apply a poiseuille flow in the x direction
u.n[left] = dirichlet(0.25 * (1. - sq(y)));
//u.n[left] = dirichlet(0.25*cos(t) * (1. - sq(y)));
u.t[left] = dirichlet(0.);
u.n[right] = neumann(0.);
u.t[right] = dirichlet(0.);
p[left]  = neumann(0.);
pf[left]  = neumann(0.);

    initial event

    here we calculate the \mu based on womersley number \alpha.

    event init (t = 0) {
  double viscosity = 1. / sq(alpha_w);
  const face vector muc[] = {viscosity , viscosity};
  mu = muc;
  fprintf(stderr,"Corresponding viscosity: %g\n", viscosity);   
  Pout =0.;
  Pold =0.;
  Qold =0.;
}

    Time integration

    In each iteration, we calculate the value of the mass flow and apply a new value to the output pressure based on it. We record the pressure value of a choosen point in order to check if the result is converged after the simulation (periodic form).

    event midpressure(t <= tmax; i += 1) {
  DEBIT = 0.; 
  // we calculate the Q
  foreach_boundary(right){
    DEBIT += u.x[]*Delta;
  }

    DETAILS OF THE THEORY CALCULATION

    Here R_e corresponding the resistance equivelent of our model,R_1 is the first resistance outside, $R_2 & C $ are the paraller RC after, We can define 2 variable:

    \alpha = \frac{R_e+R_1}{R_2}+1 and \beta = 1 - \frac{R_e}{R_{tot}}

    We can see in the video that

    \mathbf{i}_{tot} = \frac{V}{R_{tot}}(1-exp(\frac{-\alpha t}{(R_e+R_1)C})) +\frac{V}{R_e+R_1}exp(\frac{-\alpha t}{(R_e+R_1)C})

    With simple caculation we can get our theory outlet pressure

    V_{out} = \beta V (1.-exp(\frac{-\alpha t}{(R_e+R_1)C}))

      double RRe =0.8;
  double alpha = (RRe + ORR) / ORC + 1.;
  double beta =  1. - RRe/(RRe + ORR +ORC);
  double QQQ = beta * (1.-exp( -alpha*t/((RRe+ORR)*OCC)))* DEBIT;

    DETAILS OF THE OUTPUT PRESSURE CALCULATION

    The equation od the KV model is

    $P + R_2 C _t P = (R_1 + R_2) Q + R_1 R_2 C _t Q $

    which can be write in the form

    P (1 + \frac{R_2 C}{dt} )= \frac{R_2 C}{dt} P_{old} + (R_1 +R_2)Q + R_1 R_2 C \frac{Q - Q_{old}}{dt}

      double term1 = ORC * OCC * Pold / dt;
  double term2 = (ORR + ORC) * DEBIT ;
  double term3 = ORR * ORC * OCC * (DEBIT - Qold) / dt;
  double deno = 1. + ORC * OCC / dt ;
  //Pout = t < tmax / 10.? (ORR * DEBIT) : ((term1 + term2 + term3) / deno);
  Pout =  (term1 + term2 + term3) / deno;
  fprintf(fp1,"%g %g %g %g %g %g\n" , t, DEBIT, Qold, Pout, Pold,QQQ);
  PPPP = (ORR +ORC) * DEBIT;

    We set the new value of outlet pressure to BC, and give the P Q to the Pold Qold

      p[right]  = dirichlet(Pout);
  pf[right]  = dirichlet(Pout); 
  Pold = Pout;
  Qold = DEBIT;
}

    We output the pressure in position y=0.5.

    event tracer (t = end) {
  for (double xx = 0.; xx < 1.; xx += L0/100. ){
    double pyy = L0/2.;
    fprintf(fp2,"%d %g %g %g %g %g\n", N, t, xx, interpolate(p , xx, pyy), Pout,PPPP);
  }
}

    Mesh adaptation

    We adapt the mesh according to the error on the volume fraction field and the velocity.

    #if ADAPT
event adapt (i++) {
  adapt_wavelet((scalar *){u}, (double[]){5e-4,1e-3}, MAX_LEVEL, MIN_LEVEL) ;
}
#endif

    Results

    plot 'testpoint' u 1:2 w l t'debit' ,\
'testpoint' u 1:3 w l t'Q' ,\
'testpoint' u 1:4 w l t'Pressure' ,\
'testpoint' u 1:5 w l t'Pold'
    testpoint (script)

    testpoint (script)

    plot 'testpoint' u 1:5 w l t'Pressure-time' ,\
'testpoint' u 1:6 w l t'theory'
    convergence (script)

    convergence (script)

    plot 'pressure.dat' u 3:4 w l t'pressure' ,\
'pressure.dat' u 3:5 w l t'pressure out'
    pressure at y = L0/2 (script)

    pressure at y = L0/2 (script)