wave-propagation2.c

Inviscid wave propagation in a straight elastic artery
Results for second 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 = minmod;
}

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 second 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$ .

$a/a_0$ for $N=256$ . (script)

$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$ . We compare the results with those obtained with a first-order scheme.

Spatial convergence for $q$ (script)