tourniquet.c

Inviscid tourniquet
Results for first order

Inviscid tourniquet

We solve here the inviscid 1D blood flow equations in a straight artery. The artery is initially inflated on its left half, mimicking the action of a tourniquet. At t>0, the vessel relaxes towards its steady-state at rest. The solution of the flow of blood generated by this elastic relaxation is obtained numerically and compared to the solution of the corresponding Riemann problem.

Analytic solution of the corresponding Riemann problem

The initial conditions are: \displaystyle A(t=0,\, x) = \left\{ \begin{aligned} & A_L \quad \mathrm{if} \: x < x_m \\ & A_R \quad \mathrm{if} \: x >= x_m \end{aligned} \right. , \qquad Q(t=0,\, x) = 0, with A_L > A_R and x_m = 0. Using the method of characteristic, we obtain the following analytic solution: \displaystyle \begin{aligned} & \mathrm{if} \: x < x_A = - c_L t \quad \left\{ \begin{aligned} & A(t,\, x) = A_L \\ & U(t,\, x) = 0 \end{aligned} \right. \\ & \mathrm{if} \: x_A <= x < x_B = (u_M - c_M)t \quad \left\{ \begin{aligned} & U(t,\, x) = \frac{4}{5}\frac{x}{t} + \frac{4}{5}c_L \\ & c(t,\, x) = -\frac{1}{5}\frac{x}{t} + \frac{4}{5}c_L \end{aligned} \right. \\ & \mathrm{if} \: x_B <= x < x_C = \frac{A_M U_M}{A_M - A_R} t \quad \left\{ \begin{aligned} & A(t,\, x) = A_M \\ & U(t,\, x) = U_M \end{aligned} \right. \\ & \mathrm{if} \: x_C <= x \quad \left\{ \begin{aligned} & A(t,\, x) = A_R \\ & U(t,\, x) = 0 \end{aligned} \right. . \end{aligned} The variables A_M and U_M are obtain by solving the following system iteratively: \displaystyle \left\{ \begin{aligned} & U_M + 4 c_M = 4 c_L\\ & Q_M = \frac{A_M U_M}{A_M - A_R} \left[A_M - A_R\right]\\ & \left[ \frac{Q_M^2}{A_M} + \frac{K}{3\rho} A_M^{\frac{3}{2}} \right] - \frac{K}{3 \rho} A_R^{\frac{3}{2}} = \frac{A_M U_M}{A_M - A_R} Q_M . \end{aligned} \right.

double celerity (double a, double k)
{
  return sqrt(0.5*k*sqrt(a));
}

double analytic_a (double t, double x, double dr,
		   double a0, double k0)
{
  double al = a0*pow(1 + dr, 2.);
  double cl = celerity(al, k0);

  double ar = a0;

  double am = 3.459578046858399;
  double cm = celerity(am, k0);
  double um = 9.192473939896399;
  double s = am*um/(am-ar);

  double xa = -cl*t;
  double xb = (um - cm)*t;
  double xc = s*t;

  if (x <= xa) return al;
  else if (xa < x && x <= xb) {
    return pow(2./k0, 2.)*pow( -1./5.*x/t + 4./5.*cl, 4.);
  }
  else if (xb < x && x <= xc) return am;
  else return ar;	     
}

double analytic_u (double t, double x, double dr,
		  double a0, double k0)
{
  double al = a0*pow(1 + dr, 2.);
  double cl = celerity(al, k0);

  double ar = a0;
  
  double am = 3.459578046858399;
  double cm = celerity(am, k0);
  double um = 9.192473939896399;
  double s = am*um/(am-ar);

  double xa = -cl*t;
  double xb = (um - cm)*t;
  double xc = s*t;

  if (x <= xa) return 0.;
  else if (xa < x && x <= xb) return  4./5.*x/t + 4./5.*cl;
  else if (xb < x && x <= xc) return um;
  else return 0.;	     
}

double analytic_q (double t, double x, double dr,
		  double a0, double k0)
{
  return analytic_a(t, x, dr, a0, k0)*analytic_u(t, x, dr, a0, k0);
}

#include "grid/cartesian1D.h"
#include "../bloodflow-hr.h"

We define the artery’s geometrical and mechanical properties.

#define R0 (1.)
#define XM (0.)
#define DR (1.e-1)
#define shape(x) ((x) < (XM) ? (1. + (DR)) : 1.)

#define K0 (1.e4)

int main() {

The domain is 10.

  L0 = 10.;
  size (L0);
  origin (-(L0)/2.);
  
  DT = 1.e-5;

We run the computation for different grid sizes.

  for (N = 128; N <= 1024; N *= 2) {
    init_grid (N);
    run();
  }
}

Boundary conditions

We impose homogeneous Neumann boundary conditions on all variables.

a[left] = neumann (0.);
q[left] = neumann (0.);

a[right] = neumann (0.);
q[right] = neumann (0.);

Defaults conditions

event defaults (i = 0)
{
  gradient = zero;
}

Initial conditions

event init (i = 0) {

We initialize the variables k, zb, a and q.

  foreach() {
    k[] = (K0);
    zb[] = k[]*sqrt (pi)*(R0);
    a[] = sq (zb[]/k[]*(shape (x)));
    q[] = 0.;
  }
}

Post-processing

We output the computed fields.

event field (t = {0., 0.01, 0.02, 0.03, 0.04}) {

  if (N == 1024) {

    char name[80];
    sprintf (name, "fields-%.2f-pid-%d.dat", t, pid());
    FILE * ff = fopen (name, "w");
    
    foreach()
      fprintf (ff, "%g %g %g %g %g %g %g\n",
	       x, k[], sq (zb[]/k[]),
	       (analytic_a (t,x,(DR),(pi*sq ((R0))),(K0)))/(pi*sq ((R0))),
	       a[]/(pi*sq ((R0))),
	       (analytic_q (t,x,(DR),(pi*sq ((R0))),(K0))),
	       q[]
	       );
  }
}

Next, we compute the spatial error for the cross-sectional area a and the flow rate q.

event error (t = 0.04) {

  scalar err_a[], err_q[];
  foreach() {
    err_a[] = fabs (a[] - analytic_a (t, x, (DR), (pi*sq ((R0))), (K0)));
    err_q[] = fabs (q[] - analytic_q (t, x, (DR), (pi*sq ((R0))), (K0)));
  }
  boundary ((scalar *) {err_a, err_q});
  
  norm na = normf (err_a);
  norm nq = normf (err_q);
  
  fprintf (ferr, "%d %g %g %g %g %g %g\n",
	   N,
	   na.avg, na.rms, na.max,
	   nq.avg, nq.rms, nq.max);
}

End of simulation

event stop_run (t = 0.05)
{
  return 0;
}

Results for first order

Cross-sectional area and flow rate

We first plot the spatial evolution of the cross-sectional area a at t={0, 0.01, 0.02, 0.03, 0.04} for N=1024.

reset
set xlabel 'x'
set ylabel 'a/a_0'
plot '< cat fields-0.00-pid-*' u 1:4 w l lw 3 lc rgb "black" t 'analytic', \
     '< cat fields-0.01-pid-*' u 1:4 w l lw 3 lc rgb "black" notitle, \
     '< cat fields-0.02-pid-*' u 1:4 w l lw 3 lc rgb "black" notitle, \
     '< cat fields-0.03-pid-*' u 1:4 w l lw 3 lc rgb "black" notitle, \
     '< cat fields-0.04-pid-*' u 1:4 w l lw 3 lc rgb "black" notitle, \
     '< cat fields-0.00-pid-*' u 1:5 w l lw 2 lc rgb "blue" t 't=0', \
     '< cat fields-0.01-pid-*' u 1:5 w l lw 2 lc rgb "red" t 't=0.01', \
     '< cat fields-0.02-pid-*' u 1:5 w l lw 2 lc rgb "sea-green" t 't=0.02', \
     '< cat fields-0.03-pid-*' u 1:5 w l lw 2 lc rgb "coral" t 't=0.03', \
     '< cat fields-0.04-pid-*' u 1:5 w l lw 2 lc rgb "dark-violet" t 't=0.04'

a/a_0 for N=1024. (script)

reset
set key bottom right
set xlabel 'x'
set ylabel 'q'
plot '< cat fields-0.00-pid-*' u 1:6 w l lw 3 lc rgb "black" t 'analytic', \
     '< cat fields-0.01-pid-*' u 1:6 w l lw 3 lc rgb "black" notitle, \
     '< cat fields-0.02-pid-*' u 1:6 w l lw 3 lc rgb "black" notitle, \
     '< cat fields-0.03-pid-*' u 1:6 w l lw 3 lc rgb "black" notitle, \
     '< cat fields-0.04-pid-*' u 1:6 w l lw 3 lc rgb "black" notitle, \
     '< cat fields-0.00-pid-*' u 1:7 w l lw 2 lc rgb "blue" t 't=0', \
     '< cat fields-0.01-pid-*' u 1:7 w l lw 2 lc rgb "red" t 't=0.01', \
     '< cat fields-0.02-pid-*' u 1:7 w l lw 2 lc rgb "sea-green" t 't=0.02', \
     '< cat fields-0.03-pid-*' u 1:7 w l lw 2 lc rgb "coral" t 't=0.03', \
     '< cat fields-0.04-pid-*' u 1:7 w l lw 2 lc rgb "dark-violet" t 't=0.04'

q for N=1024. (script)

Convergence

Finally, we plot the evolution of the error for the cross-sectional area a and the flow rate q with the number of cells N.

reset
set xlabel 'N'
set ylabel 'L_1(a),L_2(a),L_{max}(a)'
set format y '%.1e'
set logscale

ftitle(a,b) = sprintf('order %4.2f', -b)

f1(x) = a1 + b1*x
f2(x) = a2 + b2*x
f3(x) = a3 + b3*x
fit f1(x) 'log' u (log($1)):(log($2)) via a1, b1
fit f2(x) 'log' u (log($1)):(log($3)) via a2, b2
fit f3(x) 'log' u (log($1)):(log($4)) via a3, b3

plot 'log' u 1:2 w p pt 6 ps 1.5 lc rgb "blue" t '|a|_1, '.ftitle(a1, b1), \
     exp (f1(log(x))) ls 1 lc rgb "red" notitle, \
     'log' u 1:3 w p pt 7 ps 1.5 lc rgb "navy" t '|a|_2, '.ftitle(a2, b2), \
     exp (f2(log(x))) ls 1 lc rgb "red" notitle, \
     'log' u 1:4 w p pt 5 ps 1.5 lc rgb "skyblue" t '|a|_{max}, '.ftitle(a3, b3), \
     exp (f3(log(x))) ls 1 lc rgb "red" notitle

Spatial convergence for a (script)

reset
set xlabel 'N'
set ylabel 'L_1(q),L_2(q),L_{max}(q)'
set format y '%.1e'
set logscale

ftitle(a,b) = sprintf('order %4.2f', -b)

f1(x) = a1 + b1*x
f2(x) = a2 + b2*x
f3(x) = a3 + b3*x
fit f1(x) 'log' u (log($1)):(log($5)) via a1, b1
fit f2(x) 'log' u (log($1)):(log($6)) via a2, b2
fit f3(x) 'log' u (log($1)):(log($7)) via a3, b3

plot 'log' u 1:5 w p pt 6 ps 1.5 lc rgb "blue" t '|q|_1, '.ftitle(a1, b1), \
     exp (f1(log(x))) ls 1 lc rgb "red" notitle, \
     'log' u 1:6 w p pt 7 ps 1.5 lc rgb "navy" t '|q|_2, '.ftitle(a2, b2), \
     exp (f2(log(x))) ls 1 lc rgb "red" notitle, \
     'log' u 1:7 w p pt 5 ps 1.5 lc rgb "skyblue" t '|q|_{max}, '.ftitle(a3, b3), \
     exp (f3(log(x))) ls 1 lc rgb "red" notitle

Spatial convergence for q (script)