bloodflow-hr.h

A hydrostatic reconstruction solver for the 1D blood flow equations
- Variables and parameters
- Time-integration

A hydrostatic reconstruction solver for the 1D blood flow equations

One-dimensional (1D) blood flow models have been successfully used to describe the flow of blood in large arteries of the systemic circulation. They are valuable and efficient tools to capture pulse wave propagation in large arterial networks and obtain satisfactory evaluations of average quantities such as the flow rate q, the cross-sectional area a or the pressure p.

The 1D blood equations are: \displaystyle \left\{ \begin{aligned} & \partial_t a + \partial_x q = 0 \\ & \partial_t q +\partial_x (q^2/a) = - a/\rho\partial_x p - f_r, \end{aligned} \right. where a is the cross-sectional area of the artery, q is the flow rate, p is the pressure, \rho is the blood density and f_r is the viscous friction term, defined as: \displaystyle f_r = \phi \nu q/a, with \nu the dynamic viscosity of blood. We typically set \phi = 8\pi.

To close the 1D system of equations, we introduce the following pressure law, describing the elastic behavior of the arterial wall: \displaystyle p = p_0 + k [\sqrt{a} - \sqrt{a_0}], where a_0 is the cross-sectional area at rest of the artery, k is the arterial wall rigidity and p_0 is the pressure applied on the exterior of the artery.

Variables and parameters

The primary fields are the arterial cross-sectional area a and the flow rate q. The secondary fields are the vessel properties and correspond to the cross-sectional area at rest a_0 and the wall rigidity k. We also introduce the terciary fields zb=k\sqrt{a_0}, h=k\sqrt{a} and eta=h-zb, which are equilibrium related quantities. The order of the declaration of these fields is important as it conditions the order in which the variables will be refined (k before zb, zb before a and a before eta).

scalar k[], zb[], a[], eta[];
scalar q[];

The following parameters are used to define the zero and to characterize algorithm convergence. Their values should not be modified.

double dry = 1.e-30;

The SEPS constant is used to avoid division by zero.

#undef SEPS
#define SEPS 1.e-30

Time-integration

The time-integration is performed using the generic second-order predictor-corrector scheme.

#include "predictor-corrector.h"

Setup

This generic time-integration scheme needs to know which fields are updated. This list will be updated in the defaults event below.

scalar * evolving = NULL;

We need to overload the default advance function of the predictor-corrector scheme.

trace
static void advance_bloodflow (scalar * output, scalar * input, 
			       scalar * updates, double dt)
{
  // recover scalar and vector fields from lists
  scalar ai = input[0], ao = output[0], da = updates[0];
  scalar qi = input[1], qo = output[1], dq = updates[1];

  foreach() {
    ao[] = ai[] + dt*da[];
    qo[] = qi[] + dt*dq[];

    eta[] = k[]*sqrt (ao[]) - zb[];
  }
    
  // fixme: on trees eta is defined as eta = h - zb and not ho -
  // zb in the refine_eta() and restriction_eta() functions below
  boundary ((scalar *) {ao, qo, eta});
}

When using an adaptive discretisation (i.e. a tree)., we need to make sure that eta is maintained as h-zb whenever cells are refined or restricted.

#if TREE
static void refine_eta (Point point, scalar eta)
{
  foreach_child()
    eta[] = k[]*sqrt (a[]) - zb[];
}

static void restriction_eta (Point point, scalar eta)
{
  eta[] = k[]*sqrt (a[]) - zb[];
}
#endif

HLL flux

void hll(double k, double al, double ar, double ql, double qr,
	 double Delta, double * fa, double * fq, double * dtmax)
{
  double ul = ql/al;
  double ur = qr/ar;
  double cl = sqrt (0.5*k*sqrt (al));
  double cr = sqrt (0.5*k*sqrt (ar));
  double SL = min (ul - cl, ur - cr);
  double SR = max (ul + cl, ur + cr);

  if (0. <= SL) {
    *fa = al*ul;
    *fq = al*(ul*ul + k*sqrt (al)/3.);
  }
  else if (0. >= SR) {
    *fa = ar*ur;
    *fq = ar*(ur*ur + k*sqrt (ar)/3.);
  }
  else {
    double fal = al*ul;
    double fql = al*(ul*ul + k*sqrt (al)/3.);
    double far = ar*ur;
    double fqr = ar*(ur*ur + k*sqrt (ar)/3.);
    *fa = (SR*fal - SL*far + SL*SR*(ar - al))/(SR - SL);
    *fq = (SR*fql - SL*fqr + SL*SR*(qr - ql))/(SR - SL);
  }

  double a = max (fabs (SL), fabs (SR));
  if (a > dry) {
    double dt = CFL*Delta/a;
    if (dt < *dtmax)
      *dtmax = dt;
  }
}

Update function used in predictor-corrector

trace
double update_bloodflow (scalar * evolving, scalar * updates, double dtmax) {

We first recover the currently evolving cross-sectional area and flow rate.

  scalar q = evolving[1];

The variables Fa and Fq will contain the fluxes for a and q respectively. The field S is necessary to store the topography source term.

  face vector Fa[], Fq[], S[];

The gradients, used for high-order spatial reconstruction of the conserved fields, are stored in locally-allocated fields. First-order reconstruction is used for the gradient fields. Note that to preserve the equilibrium states, the spatial reconstruction accounts for the equilibrium states the method is designed to preserve.

vector gk[], gzb[], geta[], gq[];
  for (scalar s in {gk, gzb, geta, gq}) {
    s.gradient = zero;
#if TREE
    s.prolongation = refine_linear;
#endif
  }
  gradients ({k, zb, eta, q}, {gk, gzb, geta, gq});

Left/right state reconstruction

We are ready to compute the fluxes through each face of the domain.

  foreach_face (reduction (min:dtmax)) {

We use the central values of each scalar/vector quantity and the gradients to compute the left and right states of each field at each cell face. The index l and r correspond respectively to the left and right of each cell face.

    double dx = Delta/2.;

    double kr = k[] - dx*gk.x[];
    double kl = k[-1] + dx*gk.x[-1];

    double zbr = zb[] - dx*gzb.x[];
    double zbl = zb[-1] + dx*gzb.x[-1];

    double etar = eta[] - dx*geta.x[];
    double etal = eta[-1] + dx*geta.x[-1];

    double qr = q[] - dx*gq.x[];
    double ql = q[-1] + dx*gq.x[-1];

We then compute the left and right states of a.

    double ar = sq ((zbr + etar)/kr);
    double al = sq ((zbl + etal)/kl);

As part of the hydrostatic reconstruction procedure, we define left and right states (m,p) at each cell face of the conserved variables a and q.

    double klr = max (kl, kr);
    double zblr = min (zbl, zbr);
    
    double ap = sq (max (0., zblr + etar)/klr);
    double am = sq (max (0., zblr + etal)/klr);

    double qp = qr; // qr/ar*ap
    double qm = ql; // ql/al*am

We then call the generic HLL Riemann solver and store the computed fluxes in temporary fields.

    double fa, fq;  
    hll (klr, am, ap, qm, qp, Delta*cm[]/fm.x[], &fa, &fq, &dtmax);

Topography source term

    double sr = klr*ap*sqrt (ap)/3. - kr*ar*sqrt (ar)/3.;
    double sl = klr*am*sqrt (am)/3. - kl*al*sqrt (al)/3.;
    
    foreach_dimension() {
      Fa.x[] = fm.x[]*fa;
      Fq.x[] = fm.x[]*(fq - sr); // Used to update the face 0 of the cell
      S.x[]  = fm.x[]*(fq - sl); // Used to update the face 1 of the cell
    }
  }
  boundary_flux ({Fa, S, Fq});

Updates for evolving quantities

The update for each scalar quantity is the divergence of the fluxes.

  scalar da = updates[0], dq = updates[1];
  foreach() {
    da[] = 0.;
    dq[] = 0.;
    foreach_dimension() {
      da[] += (Fa.x[] - Fa.x[1])/(cm[]*Delta);
      dq[] += (Fq.x[] - S.x[1,0])/(cm[]*Delta);
    }

An extra cell centered term must be added to the momentum update to preserve consistency when second-order reconstruction is used. This term is 0 for first-order reconstruction. In the following, to match the indices above, the indices (l,r) now refer respectively to the left of the right face and to the right of the left face of each cell.

    double dx = Delta/2.;

    double kr = k[] - dx*gk.x[];
    double kl = k[] + dx*gk.x[];

    double zbr = zb[] - dx*gzb.x[];
    double zbl = zb[] + dx*gzb.x[];

    double etar = eta[] - dx*geta.x[];
    double etal = eta[] + dx*geta.x[];

    double ar = sq ((zbr + etar)/kr);
    double al = sq ((zbl + etal)/kl);
    
    double kc = sqrt (kl*kr);
    double zbc = 0.5*(zbl + zbr);
    double acr = sq ((zbc + etar)/(kc + SEPS));
    double acl = sq ((zbc + etal)/(kc + SEPS));

    dq[] -= (
    	     (kc*acl*sqrt(acl)/3. - kl*al*sqrt (al)/3.) -
    	     (kc*acr*sqrt(acr)/3. - kr*ar*sqrt (ar)/3.)
    	     )/(cm[]*Delta);
  }
  return dtmax;
}

Initialisation and cleanup

We use the defaults event defined in predictor-corrector to setup the initial conditions.

event defaults (i = 0)
{
  evolving = list_concat ({a}, {q});

  foreach()
    for (scalar s in evolving)
      s[] = 0.;
  boundary (evolving);

  advance = advance_bloodflow;
  update = update_bloodflow;
  
#if TREE
  for (scalar s in {k, zb, a, eta, q}) {
    s.refine = s.prolongation = refine_linear;
    s.restriction = restriction_volume_average;
    }

We also define special reconstructions for eta.

  eta.refine  = refine_eta;
  eta.restriction = restriction_eta;
#endif
}

The event below will be executed after all the other initial events to take into account user-defined field initialisations.

event init (i = 0)
{
  foreach() {
    eta[] = k[]*sqrt (a[]) - zb[];
  }
  boundary (all);
}

At the end of the simulation we need to free any memory allocated in the defaults event.

event cleanup (i = end, last)
{
  free (evolving);
}

Well-balanced reconstruction of a

#include "elevation.h"