1D arterial flow

    A 1D model for arterial flows can be derived from the Navier-Stokes equations, in terms of the cross sectional area A and flow rate Q, we have \displaystyle \partial_t A +\partial_x Q = 0 \displaystyle \partial_t Q +\partial_x (Q^2/A) = - A \partial_x p/\rho - f_r where p(A) models the wall properties of the arteries, \rho is the blood density and f_r stands for the wall shear stress. For a simple linear wall relation, p = K A with K a constant, we can write the flux as F = (Q,Q^2/A + 2 e_1 A) and the source term as S = (0,-e_2 Q/A) using two parameters e_1 and e_2.

    Before including the conservation solver, we need to overload the default update function of the predictor-corrector scheme in order to add our source term.

    #include "grid/cartesian1D.h"
    #include "predictor-corrector.h"
    static double momentum_source (scalar * current, scalar * updates, double dtmax);
    event defaults (i = 0)
      update = momentum_source;
    #include "conservation.h"


    We define the conserved scalar fields a and q which are passed to the generic solver through the scalars list. We don’t have any conserved vector field.

    scalar a[], q[];
    scalar * scalars = {a,q};
    vector * vectors = NULL;

    The other parameters are specific to the example.

    double e1, e2, omega, Amp;


    We define the flux function required by the generic solver.

    void flux (const double * s, double * f, double e[2])
      double a = s[0], q = s[1], u = q/a;
      f[0] = q;
      f[1] = q*q/a + e1*a*a;
      // min/max eigenvalues
      double c = sqrt(2.*e1*a);
      e[0] = u - c; // min
      e[1] = u + c; // max

    We need to add the source term of the momentum equation. We define a function which, given the current states, fills the updates with the source terms for each conserved quantity.

    static double momentum_source (scalar * current, scalar * updates, double dtmax)

    We first compute the updates from the system of conservation laws.

      double dt = update_conservation (current, updates, dtmax);

    We recover the current fields and their variations from the lists…

      scalar a = current[0], q = current[1], dq = updates[1];

    We add the source term for q.

        dq[] += - e2*q[]/a[];
      return dt;

    Boundary conditions

    We impose a sinusoidal flux Q(t) at the left of the domain.

    q[left] = dirichlet(Amp*sin(2.*pi*omega*t));


    For small amplitudes Amp = 0.01 at the input boundary condition the system has analytical solutions for e1 < e2, in this case the spatial envelope of the flux rate behaves like Q=Amp\times e^{-e2/2x} [Wang at al., 2013].

    int main() {
      init_grid (400);
      e1 = 0.5 ;
      e2 = 0.1 ;
      omega = 1.;
      Amp = 0.01 ;

    Initial conditions

    The initial conditions are A=1 and Q=0.

    event init (i = 0) {
      theta = 1.3; // tune limiting from the default minmod
        a[] = 1.;


    We print to standard error the spatial profile of the flow rate Q.

    event printdata (t += 0.1; t <= 1.) {
        fprintf (stderr, "%g %.6f \n", x, q[]);
      fprintf (stderr, "\n\n");

    We get the following comparison between the numerical solution and the linear theory for the flow rate Q.

    Amp = 0.01
    e2 = 0.1
    set yrange [0.008:]
    set ylabel 'Q'
    set xlabel 'x'
    plot 'log' w l t 'numerical', Amp*exp(-e2/2.*x) t 'linear theory'