wave-propagation.c

Inviscid wave propagation in a straight elastic artery
Results for first order

Inviscid wave propagation in a straight elastic artery

We solve the 1D blood flow equations in a straight artery initially deformed in its center. At t>0, the vessel relaxes towards its steady-state at rest. Consequently, two waves are created that propagate towards both extremeties of the artery. The solution of the flow of blood generated by this elastic relaxation is obtained numerically and compared to the anaytic solution obtained using linear wave theory.

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

We define the artery’s geometrical and mechanical properties.

#define R0 (1.)
#define DR (1.e-3)
#define XS (2./5.*(L0))
#define XE (3./5.*(L0))
#define shape(x) (((XS) <= (x) && (x) <= (XE)) ? (DR)/2.*(1. + cos(pi + 2.*pi*((x) - (XS))/((XE) - (XS)))) : 0.)

#define K0 1.e4
#define C0 (sqrt (0.5*(K0)*sqrt(pi)*(R0)))

We define the linear analytic solution for the cross-sectional area a and the flow rate q, which depends on the shape of initial vessel deformation.

#define analytic_a(t,x) (pi*sq (R0)*sq (1. + 0.5*((shape ((x) - (C0)*(t))) + (shape ((x) + (C0)*(t))))))
#define analytic_q(t,x) (-(C0)*(-(shape ((x) - (C0)*(t))) + (shape ((x) + (C0)*(t))))*(analytic_a((t),(x))))

int main() {

The domain is 10.

  L0 = 10.;
  size (L0);
  origin (0.);
  
  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[]*(1. + (shape (x))));
    q[] = 0.;
  }
}

Post-processing

We output the computed fields.

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

  if (N == 256) {
    
    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))/(pi*sq ((R0))), a[]/(pi*sq ((R0))),
	       (analytic_q (t,x)), q[]
	       );
  }
}

Next, we compute the spatial error for the flow rate.

event error (t = 0.03) {

  scalar err_q[];
  foreach()
    err_q[] = fabs (q[] - (analytic_q (t, x)));
  boundary ((scalar *) {err_q});

  norm nq = normf (err_q);
  
  fprintf (ferr, "%d %g %g %g\n",
	   N,
	   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=256.

reset
set xlabel 'x'
set ylabel 'a/a_0'
set yrange [1:1 + 2.5e-3]
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=256. (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=256. (script)

Convergence

Finally, we plot the evolution of the error for the flow rate q with the number of cells N.

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($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 '|q|_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 '|q|_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 '|q|_{max}, '.ftitle(a3, b3), \
     exp (f3(log(x))) ls 1 lc rgb "red" notitle

Spatial convergence for q (script)