sandbox/ghigo/artery1D/tourniquet.c
Inviscid tourniquet
We solve the inviscid 1D blood flow equations in a straight artery. The artery is initially inflated on its left half, mimicking the action of a tourniquet. At t>0, the vessel relaxes towards its steady-state at rest. The solution of the flow of blood generated by this elastic relaxation is obtained numerically and compared to the solution of the corresponding Riemann problem.
Analytic solution of the corresponding Riemann problem
The initial conditions are:
\displaystyle A(t=0,\, x) = \left\{ \begin{aligned} & A_L \quad \mathrm{if} \: x < x_m \\ & A_R \quad \mathrm{if} \: x >= x_m \end{aligned} \right. , \qquad Q(t=0,\, x) = 0,
with A_L > A_R and x_m = 0. Using the characteristic method, we obtain the following analytic solution:
\displaystyle \begin{aligned} & \mathrm{if} \: x < x_A = - c_L t \quad \left\{ \begin{aligned} & A(t,\, x) = A_L \\ & U(t,\, x) = 0 \end{aligned} \right. \\ & \mathrm{if} \: x_A <= x < x_B = (u_M - c_M)t \quad \left\{ \begin{aligned} & U(t,\, x) = \frac{4}{5}\frac{x}{t} + \frac{4}{5}c_L \\ & c(t,\, x) = -\frac{1}{5}\frac{x}{t} + \frac{4}{5}c_L \end{aligned} \right. \\ & \mathrm{if} \: x_B <= x < x_C = \frac{A_M U_M}{A_M - A_R} t \quad \left\{ \begin{aligned} & A(t,\, x) = A_M \\ & U(t,\, x) = U_M \end{aligned} \right. \\ & \mathrm{if} \: x_C <= x \quad \left\{ \begin{aligned} & A(t,\, x) = A_R \\ & U(t,\, x) = 0 \end{aligned} \right. . \end{aligned}
The variables A_M and U_M are obtain by solving the following system iteratively:
\displaystyle \left\{ \begin{aligned} & U_M + 4 c_M = 4 c_L\\ & Q_M = \frac{A_M U_M}{A_M - A_R} \left[A_M - A_R\right]\\ & \left[ \frac{Q_M^2}{A_M} + \frac{K}{3\rho} A_M^{\frac{3}{2}} \right] - \frac{K}{3 \rho} A_R^{\frac{3}{2}} = \frac{A_M U_M}{A_M - A_R} Q_M . \end{aligned} \right.
Code
#include "grid/cartesian1D.h"
#include "bloodflow.h"
double R0 = 1.;
double K0 = 1.e4;
double L = 10.;
double XM = 0.;
double DR = 1.e-1;
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 dr,
double a0, double k0) {
double al = a0*pow(1 + dr, 2.);
double cl = celerity(al, k0);
double ar = a0;
double am = 3.459578046858399;
double cm = celerity(am, k0);
double um = 9.192473939896399;
double s = am*um/(am-ar);
double xa = -cl*t;
double xb = (um - cm)*t;
double xc = s*t;
if (x <= xa) return al;
else if (xa < x && x <= xb) {
return pow(2./k0, 2.)*pow( -1./5.*x/t + 4./5.*cl, 4.);
}
else if (xb < x && x <= xc) return am;
else return ar;
}
double analyticU (double t, double x, double dr,
double a0, double k0) {
double al = a0*pow(1 + dr, 2.);
double cl = celerity(al, k0);
double ar = a0;
double am = 3.459578046858399;
double cm = celerity(am, k0);
double um = 9.192473939896399;
double s = am*um/(am-ar);
double xa = -cl*t;
double xb = (um - cm)*t;
double xc = s*t;
if (x <= xa) return 0.;
else if (xa < x && x <= xb) return 4./5.*x/t + 4./5.*cl;
else if (xb < x && x <= xc) return um;
else return 0.;
}
int main() {
origin (-L/2., 0.);
L0 = L;
for (N = 128; N <= 1024; 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] = neumann(0.);
q[right] = neumann(0.);
event defaults (i=0) {
gradient = order1;
riemann = hll_glu;
}
event init (i=0) {
foreach() {
k[] = K0;
a0[] = pi*pow(R0, 2.);
a[] = a0[]*(x < XM ? pow(1 + DR, 2.) : 1);
q[] = 0.;
}
}
event field (t = {0., 0.01, 0.02, 0.03, 0.04}) {
if (N == 1024) {
foreach() {
fprintf (stderr, "%g, %.6f, %.6f\n", x, a[], q[]) ;
}
fprintf (stderr, "\n\n") ;
}
}
event error (i++) {
ne++;
foreach() {
ea[] = a[] - analyticA (t, x, DR, a0[], k[]);
eq[] = q[] - analyticU (t, x, DR, a0[], k[])*
analyticA (t, x, DR, a0[], k[]);
}
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 = 4.e-2) {
printf ("#Tourniquet test case completed\n") ;
}Plots
reset
mycolors = "dark-blue red sea-green dark-violet orange"
mypoints = '1 2 3 4 6'
set datafile separator ','
set style line 1 lt -1 lw 2 lc 'black'
k0 = 1.e4
a0 = pi
dr = 0.1
al = a0*(1 + dr)**2.
cl = sqrt(0.5*k0*sqrt(al))
ar = a0
am = 3.459578046858399;
cm = sqrt(0.5*k0*sqrt(am))
um = 9.192473939896399;
s = 100.01113797047884;
xa(t) = -cl*t;
xb(t) = (um - cm)*t;
xc(t) = s*t;
a(t, x) = x <= xa(t) ? al : \
xa(t) < x && x <= xb(t) ? ((2./k0)**2.)*(-1./5.*x/t + 4./5.*cl)**4. : \
xb(t) < x && x <= xc(t) ? am : ar
u(t, x) = x <= xa(t) ? 0. : \
xa(t) < x && x <= xb(t) ? 4./5.*x/t + 4./5.*cl : \
xb(t) < x && x <= xc(t) ? um : 0.
q(t, x) = a(t, x)*u(t, x)
set output 'A.png'
set xlabel 'x [cm]'
set ylabel 'A [cm^2]'
set key top right
plot a(0., x) w l ls 1 t 'Analytic', \
a(0.01, x) w l ls 1 notitle, \
a(0.02, x) w l ls 1 notitle, \
a(0.03, x) w l ls 1 notitle, \
a(0.04, x) w l ls 1 notitle, \
'log' i 0 u 1:2 every 5 w p lt word(mypoints, 1) lc rgb word(mycolors, 1) ps 2 lw 1 t 't=0', \
'log' i 1 u 1:2 every 5 w p lt word(mypoints, 2) lc rgb word(mycolors, 2) ps 2 lw 1 t 't=0.01', \
'log' i 2 u 1:2 every 5 w p lt word(mypoints, 3) lc rgb word(mycolors, 3) ps 2 lw 1 t 't=0.02', \
'log' i 3 u 1:2 every 5 w p lt word(mypoints, 4) lc rgb word(mycolors, 4) ps 2 lw 1 t 't=0.03', \
'log' i 4 u 1:2 every 5 w p lt word(mypoints, 5) lc rgb word(mycolors, 5) ps 2 lw 1 t 't=0.04'
{0, 0.01, 0.02, 0.03, 0.04}
set output 'Q.png'
set xlabel 'x [cm]'
set ylabel 'Q [cm^3/s]'
set key center center
plot q(0., x) w l ls 1 t 'Analytic', \
q(0.01, x) w l ls 1 notitle, \
q(0.02, x) w l ls 1 notitle, \
q(0.03, x) w l ls 1 notitle, \
q(0.04, x) w l ls 1 notitle, \
'log' i 0 u 1:3 every 5 w p lt word(mypoints, 1) lc rgb word(mycolors, 1) ps 2 lw 1 t 't=0', \
'log' i 1 u 1:3 every 5 w p lt word(mypoints, 2) lc rgb word(mycolors, 2) ps 2 lw 1 t 't=0.01', \
'log' i 2 u 1:3 every 5 w p lt word(mypoints, 3) lc rgb word(mycolors, 3) ps 2 lw 1 t 't=0.02', \
'log' i 3 u 1:3 every 5 w p lt word(mypoints, 4) lc rgb word(mycolors, 4) ps 2 lw 1 t 't=0.03', \
'log' i 4 u 1:3 every 5 w p lt word(mypoints, 5) lc rgb word(mycolors, 5) ps 2 lw 1 t 't=0.04'
{0, 0.01, 0.02, 0.03, 0.04}
set xtics 128,2,1024
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) '< grep ^ea out' u (log($2)):(log($3)) via a1, b1
fit f2(x) '< grep ^ea out' u (log($2)):(log($4)) via a2, b2
fit f3(x) '< grep ^ea out' u (log($2)):(log($5)) via a3, b3
set output 'errA.png'
set xlabel 'N'
set ylabel 'Error norms'
set key top right
set cbrange[1e-10:1e10]
plot '< grep ^ea out' u 2:3 w p lt word(mypoints, 1) lc rgb word(mycolors, 1) ps 3 lw 2 t '|A|_1, '.ftitle(a1, b1), \
exp (f1(log(x))) ls 1 notitle, \
'< grep ^ea out' u 2:4 w p lt word(mypoints, 2) lc rgb word(mycolors, 2) ps 3 lw 2 t '|A|_2, '.ftitle(a2, b2), \
exp (f2(log(x))) ls 1 notitle, \
'< grep ^ea out' u 2:5 w p lt word(mypoints, 3) lc rgb word(mycolors, 3) ps 3 lw 2 t '|A|_{max}, '.ftitle(a3, b3), \
exp (f3(log(x))) ls 1 notitle
Comparison of the evolution of the error norms for the cross-sectional area A with the number of cells N (script)
fit f1(x) '< grep ^eq out' u (log($2)):(log($3)) via a1, b1
fit f2(x) '< grep ^eq out' u (log($2)):(log($4)) via a2, b2
fit f3(x) '< grep ^eq out' u (log($2)):(log($5)) via a3, b3
set output 'errQ.png'
set xlabel 'N'
set ylabel 'Error norms'
set key top right
plot '< grep ^eq out' u 2:3 w p lt word(mypoints, 1) lc rgb word(mycolors, 1) ps 3 lw 2 t '|Q|_1, '.ftitle(a1, b1), \
exp (f1(log(x))) ls 1 notitle, \
'< grep ^eq out' u 2:4 w p lt word(mypoints, 2) lc rgb word(mycolors, 2) ps 3 lw 2 t '|Q|_2, '.ftitle(a2, b2), \
exp (f2(log(x))) ls 1 notitle, \
'< grep ^eq out' u 2:5 w p lt word(mypoints, 3) lc rgb word(mycolors, 3) ps 3 lw 2 t '|Q|_{max}, '.ftitle(a3, b3), \
exp (f3(log(x))) ls 1 notitle
Comparison of the evolution of the error norms for the flow rate Q with the number of cells N (script)
