sandbox/yonghui/biomeca/3darchrrc.c
real 3d arch with KV model
This is a rough attempt to merge all the solvers together, plz understande the code by reading those simple exemples. RRC & Womersley & ReadSTL & MassFlow &
#include "grid/octree.h"
#include "navier-stokes/centered.h"
#include "utils.h"
#include "distance.h"
#include "fractions.h"
#include "view.h"
#include "lambda2.h"
#define MIN_LEVEL 5
#define LEVEL 7
#define MAX_LEVEL 10
#define eq_r 1.
#define ADAPT 0 //1 = use
#define tmax 1*2.*M_PI
#define alpha_w 10.
#define ORR 0.1
#define ORC 0.1
#define OCC 10
FILE * fp1;
FILE * fp2;
FILE * fp3;
double Pout;
double Pold;
double DEBIT;
double Qold;
double PPPP;
Importing the geometry
void fraction_from_stl (scalar f, FILE * fp)
{
We read the STL file and compute the bounding box of the model.
coord * p = input_stl (fp);
coord min, max;
bounding_box (p, &min, &max);
double maxl = -HUGE;
foreach_dimension()
if (max.x - min.x > maxl)
maxl = max.x - min.x;
scalar d[];
distance (d, p);
vertex scalar phi[];
foreach_vertex()
phi[] = (d[] + d[-1] + d[0,-1] + d[-1,-1] +
d[0,0,-1] + d[-1,0,-1] + d[0,-1,-1] + d[-1,-1,-1])/8.;
fractions (phi, f);
}
volume fraction field declaire
scalar f0[];
int main()
{
init_grid (8);
size (8.);
origin (-1*L0/4 ,-2*L0/5, 0.);
fp1 = fopen ("testpoint", "w");
fp2 = fopen ("pressure.dat", "w");
fp3 = fopen ("uuu", "w");
run();
}
Boundary Conditions
u.n[back] = ( y <= 0)? dirichlet(10.*f0[]*cos(t)):neumann(0);
u.t[back] = ( y <= 0)? dirichlet(0):neumann(0);
u.r[back] = ( y <= 0)? dirichlet(0):neumann(0);
p[front] = neumann(0) ;
pf[front] = neumann(0) ;
Initial Event
event init (t = 0) {
if (!restore (file = "restart")) {
FILE * fp = fopen ("aortic_arch.stl", "r");
char name2[10];
for (int ii=MIN_LEVEL; ii<=LEVEL; ii++){
fraction_from_stl (f0, fp);
adapt_wavelet ({f0}, (double[]){0.00001}, ii);
//dump file
sprintf (name2, "dump_%d",ii);
dump(file = name2);
}
fraction_from_stl (f0, fp);
view (fov = 19.9677, quat = {0.614306,0.390179,0.155159,0.66807},
tx = -0.120804, ty = -0.340923, width = 1920, height = 1080);
draw_vof ("f0");
squares ("f0", min = -10, max =10);
save ("init.png");
double viscosity = sq(eq_r) / sq(alpha_w);
const face vector muc[] = {viscosity, viscosity, viscosity};
mu = muc;
}
DEBIT = 0.;
Pout =0.;
Pold =0.;
Qold =0.;
}
event bc (t <= tmax; i += 1) {
//BC : lateral velocity = 0
foreach()
foreach_dimension()
u.x[] = f0[]*u.x[];
boundary ((scalar *){u , p});
inlet surface position [-1:3],[-3:0.5] in xy plan, cut z=0.4 seems good, you can also choose a plan by yourself (Recommande try it first in bview3D).
DEBIT = 0.;
double zz = 0.4;
for (double xx = -1.; xx < 3.; xx += 0.05 ){
for (double yy = -3.; yy < 0.5; yy += 0.05 ){
DEBIT += interpolate(u.z, xx, yy, zz)*sq(0.05);
}
}
RRC model , in this case, it seems still has some problem
double Reee = 0.5;
double alpha = (Reee + ORR) / ORC + 1.;
double beta = 1. - Reee/(Reee + ORR +ORC);
double QQQ = beta * (1.-exp( -alpha*t/((Reee+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 ;
//out puressure calculate //RRC model
Pout = (term1 + term2 + term3) / deno;
//Pout = ORR * DEBIT;
//Pout=0. ;
fprintf(fp1,"%g %g %g %g %g\n" , t, DEBIT, Qold, Pout, Pold);
//position du test point
double px = 3.5;
double py = 0.8;
double pz = 2.;
fprintf(stderr, "%d %g %g %g \n", i, t, dt,interpolate(u.x , px, py, pz));
p[back] = ( y >= 0)? dirichlet(Pout):neumann(0) ;
pf[back] = ( y >= 0)? dirichlet(Pout):neumann(0) ;
//we give the P Q to the Pold Qold
Pold = Pout;
Qold = DEBIT;
}
Output
you can delete some commande, it woundn’t change the simulation
event snapshot (i += 1000)
{
char name[80];
sprintf (name, "dump-%d", i);
scalar l2[];
lambda2 (u,l2);
dump (file = name);
}
event video (t += tmax/300.){
scalar l2[];
lambda2 (u,l2);
input and output view
view (fov = 20, quat = {0,0,0,1}, tx = -0.19684, ty = 0.00981146, bg = {0.3,0.4,0.6}, width = 640, height = 640);
// draw_vof("f0");
box();
squares ("u.z", alpha = 0.4);
cells( alpha = 0.4);
save ("input_uz.mp4");
main chanel view
clear();
view (fov = 18.0787, quat = {0.126228,-0.561414,-0.620277,0.533052}, tx = 0.00499754, ty = -0.373772, bg = {0.3,0.4,0.6}, width = 640, height = 640, samples = 1);
//draw_vof("f0");
squares("u.z",n = {1,-1,0},alpha =2.4 );
squares("u.z",n = {-1,-1,1.5},alpha = -3.5 );
cells(n = {1,-1,0},alpha =2.4 );
cells(n = {-1,-1,1.5},alpha = -3.5 );
save ("mainchannel_uz.mp4");
clear();
view (fov = 18.0787, quat = {0.126228,-0.561414,-0.620277,0.533052}, tx = 0.00499754, ty = -0.373772, bg = {0.3,0.4,0.6}, width = 640, height = 640, samples = 1);
//draw_vof("f0");
squares("l2",n = {1,-1,0},alpha =2.4 );
squares("l2",n = {-1,-1,1.5},alpha = -3.5 );
cells(n = {1,-1,0},alpha =2.4 );
cells(n = {-1,-1,1.5},alpha = -3.5 );
save ("mainchannel_l2.mp4");
clear();
view (fov = 18.0787, quat = {0.126228,-0.561414,-0.620277,0.533052}, tx = 0.00499754, ty = -0.373772, bg = {0.3,0.4,0.6}, width = 640, height = 640, samples = 1);
//draw_vof("f0");
squares("u.x",n = {1,-1,0},alpha =2.4 );
squares("u.x",n = {-1,-1,1.5},alpha = -3.5 );
cells(n = {1,-1,0},alpha =2.4 );
cells(n = {-1,-1,1.5},alpha = -3.5 );
save ("mainchannel_ux.mp4");
cross section view
view (fov = 11.1239, quat = {0.658889,-0.274942,-0.26525,0.648006}, tx = -0.207846, ty = -0.306235, bg = {0.3,0.4,0.6}, width = 640, height = 640, samples = 1);
// draw_vof("f0");
squares("u.x",n = {-1,-1,0},alpha =-2.4);
cells(n = {-1,-1,0},alpha =-2.4);
save ("section_ux.mp4");
squares("l2",n = {-1,-1,0},alpha =-2.4);
cells(n = {-1,-1,0},alpha =-2.4);
save ("section_l2.mp4");
}
Adaption
#if ADAPT
event adapt (i++) {
double uemax = 0.01;
adapt_wavelet ({u}, (double[]){uemax,uemax,uemax}, MAX_LEVEL, 4);
}
#endif
Results
plot 'testpoint' u 1:2 w l t'debit' ,\
'testpoint' u 1:3 w l t'Q' ,\
'testpoint' u 1:4 w l t'Pressure' ,\
'testpoint' u 1:5 w l t'Pold'
plot 'pressure.dat' u 3:4 w l t'pressure' ,\
'pressure.dat' u 3:5 w l t'pressure out'