sandbox/ghigo/artery1D/viscous-dissipation.c
Viscous dissipation
We solve the propagation of a wave using the viscous 1D blood flow equations in a straight artery. As the wave propagates, viscous dissipation modifies the amplitude of the wave.
Asymptotic viscous dissipation
For a small amplitude QIN of the flow rate signal, an asymptotic solution can be obtained, where the envelop of the viscous dissipation is described by the following expression:
\displaystyle QIN * \exp\left(-\frac{CF}{2}x \right)
Code
#include "grid/cartesian1D.h"
#include "bloodflow.h"
double R0 = 1.;
double K0 = 1.e4;
double L = 200.;
double CF = 0.1;
double T0 = 1.;
double SH = 1.e-3;
double AIN = 0., QIN = 0.;
double celerity (double a, double k) {
return sqrt(0.5*k*sqrt(a));
}
double inlet (double t, double t0) {
if (t < t0) return max(0., 0.5*(1. + cos(pi + 2.*pi/T0*t)));
else return 0.;
}
int main() {
origin (0., 0.);
L0 = L;
AIN = pi*pow(R0*(1 + SH), 2.);
QIN = SH*AIN*celerity(AIN, K0);
N = 2048;
run();
return 0;
}
q[left] = dirichlet (QIN*inlet (t, T0));
event defaults (i=0) {
gradient = order1;
riemann = hll_glu;
}
event init (i=0) {
foreach() {
k[] = K0;
a0[] = pi*pow(R0, 2.);
a[] = a0[];
q[] = 0.;
}
}
event viscosity (i++, last) {
scalar sa[], sq[];
foreach() {
sa[] = 0.;
sq[] = -CF*q[]/a[];
q[] += dt*sq[];
}
boundary({sa,sq});
}
event field (t = {0.5, 1., 1.5, 2., 2.5}) {
foreach() {
fprintf (stderr, "%g, %.6f, %.6f\n", x, a[], q[]) ;
}
fprintf (stderr, "\n\n") ;
}
event end (t = 2.5) {
printf ("#Viscous dissipation 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'
set output 'A.png'
set xlabel 'x [cm]'
set ylabel 'A [cm^2]'
set key center center
plot '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.5', \
'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=1', \
'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=1.5', \
'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=2', \
'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=2.5'
{0.5, 1, 1.5, 2, 2.5}
k0 = 1.e4
a0 = pi
l = 200.
cf = 0.1
sh = 1.e-3
ain = a0*(1 + sh)**2.;
qin = sh*ain*sqrt(0.5*k0*sqrt(ain));
amp(x) = qin*exp(-cf/2.*x/l)
set output 'Q.png'
set xlabel 'x [cm]'
set ylabel 'Q [cm^3/s]'
set key right center
plot amp(x) w l ls 1 t 'Asymptotic', \
'log' i 0 u 1:3 every 5 w p lt word(mypoints, 2) lc rgb word(mycolors, 2) ps 2 lw 1 t 't=1', \
'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=1', \
'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=1.5', \
'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=2', \
'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=2.5'
{0.5, 1, 1.5, 2, 2.5}