delestre.c

Analytic solution proposed by O. Delestre
Analytic solution
Code
Plots

Analytic solution proposed by O. Delestre

We solve the inviscid 1D blood flow equations in a straight artery. The computed blood flow is compared to an analytic solution.

Analytic solution

We search for a simple solution of the form:

$\displaystyle \left\{ \begin{aligned} & A\left(x,\, t \right) = a\left(t\right)x + b\left(t\right) \\ & U\left(x,\, t \right) = c\left(t\right)x + d\left(t\right). \end{aligned} \right.$

Injecting these expressions in the inviscid 1D blood flow equations written in non-conservative form, we obtain the following ordinary differential equations for the coefficients $a$ , $b$ , $c$ and $d$ :

$\displaystyle \left\{ \begin{aligned} & a = 0 \\ & d^{'} + bc = 0 \\ & c^{'} + c^2 = 0 \\ & d^{'} + cd = 0. \end{aligned} \right.$

We finally obtain the following analytic solution for the cross-sectional area $A$ and the flow rate $Q$ :

$\displaystyle \left\{ \begin{aligned} & A\left(x,\, t\right) = \frac{C_2}{C_1 + t}\\ & U\left(x,\, t\right) = \frac{C_3 + x}{C_1 + t} \\ & Q\left(x,\, t\right) = AU, \end{aligned} \right.$

where $C_1$ , $C_2$ and $C_3$ are constants chosen through the inlet and outlet boundary conditions.

Code

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

double R0 = 1.;
double K0 = 1.e4;
double L = 10.;

double C1 = -1.;
double C2 = -pi;
double C3 = 5.;

scalar ea[], eq[];
double ea1, ea2, eamax, eq1, eq2, eqmax;
int ne = 0;

double celerity (double a, double k) {

  return sqrt(0.5*k*sqrt(a));
}

double analyticA (double t, double x,
		  double C1, double C2) {

  return (C2)/(C1 + t);
}

double analyticU (double t, double x,
		  double C1, double C3) {

  return (C3 + x)/(C1 + t);
}

int main() {

  origin (0., 0.);
  L0 = L;

  for (N = 32; N <= 256; N *= 2) {

    ea1 = ea2 = eamax = 0.;
    eq1 = eq2 = eqmax = 0.;
    ne = 0;
    
    run();

    printf ("ea, %d, %g, %g, %g\n", N, ea1/ne,
	    sqrt(ea2/ne), eamax);
    printf ("eq, %d, %g, %g, %g\n", N, eq1/ne,
	    sqrt(eq2/ne), eqmax);
  }
  
  return 0; 
}

q[left] = dirichlet(analyticA (t, 0., C1, C2)*analyticU (t, 0., C1, C3));
a[right] = dirichlet(analyticA (t, L0, C1, C2));

event defaults (i=0) {

  gradient = order1;
  riemann = hll_glu;
}

event init (i=0) {
  
  foreach() {
    k[] = K0;
    a0[] = pi*pow(R0, 2.);
    a[] = analyticA (0., x, C1, C2);
    q[] = analyticA (0., x, C1, C2)*analyticU (0., x, C1, C3);
  }
}

event field (t = {0., 0.1, 0.2, 0.3, 0.4}) {

  if (N == 128) {
    foreach() {
      fprintf (stderr, "%g, %.6f, %.6f\n", x, a[], q[]) ;
    }
    fprintf (stderr, "\n\n") ;
  }
}

event error (i++) {

  ne++;
  
  foreach() {
    ea[] = a[] - analyticA (t, x, C1, C2);
    eq[] = q[] - analyticA (t, x, C1, C2)*
      analyticU (t, x, C1, C3);
  }
  
  norm na = normf (ea);
  ea1 += na.avg;
  ea2 += na.rms*na.rms;
  if (na.max > eamax)
    eamax = na.max;
  
  norm nq = normf (eq);
  eq1 += nq.avg;
  eq2 += nq.rms*nq.rms;
  if (nq.max > eqmax)
    eqmax = nq.max;
}

event end (t = 0.4) {
  printf ("#Inviscid Delestre test case completed\n") ;
}

Plots

Spatial evolution of the cross-sectional area A computed at $t= (script) {0, 0.1, 0.2, 0.3, 0.4}

Spatial evolution of the flow rate Q computed at $t= (script) {0, 0.1, 0.2, 0.3, 0.4}

Comparison of the evolution of the relative error norms for the cross-sectional area $A$ with the number of cells $N$ (script)

Comparison of the evolution of the relative error norms for the flow rate $Q$ with the number of cells $N$ (script)