sandbox/yonghui/biomeca/torusrrc.c
blood flow in curved pipe
Idea from Antoonvh tube
Some exemples of the functions used in this code. RRC & Womersley
We want to test if the new “embed” function can better treat with the non-slip wall conditions (and the sharp velocity/pressure differences at the inlet outlet surface).
We want to simulate the blood circulation in a torus shape pipe (a simplified model for the aortic arche) with major radius R_t and minor radius r_t.
We propose a sinusoidal shape input flux and a KV model as output condition to simulate the effect of heart beat and following artery.
#include "grid/octree.h"
#include "embed.h"
#include "navier-stokes/centered.h"
#include "view.h"
face vector av[], muc[];
The torus is defined using these variables and macros:
double minor_rt = 1., major_rt = 3.;
#define R (sqrt(sq(x) + sq(y)))
#define SINPSI (y/R)
#define COSPSI (x/R)
#define TUBE (R <= 0.01 ? major_rt \
: (sqrt(sq(x - major_rt*COSPSI) + \
sq(y - major_rt*SINPSI) + sq(z))))
Macros for KV model
#define alpha_w 10.
#define ORR 0.1
#define ORC 0.1
#define OCC 10
#define tend 50.
double Pout,Pold,Qold,PPPP,DEBIT;
We set the inlet and outlet conditions on the bottum.
Z - front (z =L0) back (z=0)
Y - top bottom
X - left right
u.n[embed] = dirichlet (0);
u.t[embed] = dirichlet (0);
#if (dimension == 3)
u.r[embed] = dirichlet (0);
#endif
u.n[bottom] = ( x >= 0)? dirichlet(5.*cs[]*cos(t)):neumann(0);
u.t[bottom] = ( x >= 0)? dirichlet(0):neumann(0);
#if (dimension == 3)
u.r[bottom] = ( x >= 0)? dirichlet(0):neumann(0);
#endif
//p[bottom] = ( y >= 0)? dirichlet(0):neumann(0) ;
//pf[bottom] = ( y >= 0)? dirichlet(0):neumann(0) ;
p[top] = neumann(0) ;
pf[top] = neumann(0) ;
int j;
int main() {
L0 = 10;
X0 = Z0 = -L0/2.;
N = 64;
mu = muc;
a = av;
run();
}
event init (t = 0) {
vertex scalar phi[];
foreach_vertex()
phi[] = minor_rt - TUBE;
fractions (phi, cs, fs);
view (phi = 0.1, psi = 0., theta = 0.8);
box();
draw_vof ("cs", edges = true, lw = 0.5);
save ("vof.png");// f0.refine = f0.prolongation = fraction_refine;
Pout = DEBIT = Pold = Qold = 0.;
j=0;
}
KV model bc
event mybc (i++){
double RRe =0.8;
double alpha = (RRe + ORR) / ORC + 1.;
double beta = 1. - RRe/(RRe + ORR +ORC);
double QQQ = beta * (1.-exp( -alpha*t/((RRe+ORR)*OCC)))* DEBIT;
double term1 = ORC * OCC * Pold / dt;
double term2 = (ORR + ORC) * DEBIT ;
double term3 = ORR * ORC * OCC * (DEBIT - Qold) / dt;
double deno = 1. + ORC * OCC / dt ;
Pout = (term1 + term2 + term3) / deno;
p[bottom] = ( y >= 0)? dirichlet(Pout):neumann(0) ;
pf[bottom] = ( y >= 0)? dirichlet(Pout):neumann(0) ;
Pold = Pout;
Qold = DEBIT;
}
movie for fun
event movies (t += 0.05) {
scalar omega[];
vorticity(u,omega);
view(width=1280,height=640);
box();
cells();
squares ("u.y",max=5.,min=-5.);
//squares ("omega");
save ("cs.mp4");
}
We note the velocity of axial direction at the middle of the model (x,z=0)
event veloprofil (t += tend*0.04)
{
if (t > tend*0.8){
j += 1;
char name2[80];
sprintf (name2, "test-%d", j);
FILE *fp2;
fp2 = fopen (name2, "w");
for (double uyy = 1.5; uyy <= 4.5; uyy += 0.01)
fprintf(fp2,"%g %g %g\n", t, uyy,interpolate(u.x , 0., uyy,0.));
}
}
Dynamic adaption, with maximum resolution corresponds to a 64^3.
event adapt (i++) {
fprintf(stderr, "%d %g %g\n", i, t,interpolate(u.x, 0., 3., 0.));
adapt_wavelet ((scalar*) {cs,u},(double[]){1e-4,0.01, 0.01, 0.01}, 7);
}
event stop (t = tend);
Results
plot 'test-1' u 2:3 w l t'time1',\
'test-2' u 2:3 w l t'time2',\
'test-3' u 2:3 w l t'time3',\
'test-4' u 2:3 w l t'time4',\
'test-5' u 2:3 w l t'time5'