Commit 42f30f84 authored by Matthieu Constant's avatar Matthieu Constant
Browse files

correction g

parent a3eac645
......@@ -253,7 +253,7 @@ output = open("%s/f_%05i.vtu" %(odir, oiter), "w")
}
#if DIMENSION == 2
static void particle_drag(double u[2], double mu, double rho, double d, double c, double drag[2], double dt, double mass, double rhop, double gradp[2])
static void particle_drag(double u[2], double mu, double rho, double d, double c, double drag[2], double dt, double mass, double rhop, double gradp[2], double g)
{
double normu = hypot(u[0],u[1]);
//Reynolds/|u_p-u_f|
......@@ -269,10 +269,10 @@ static void particle_drag(double u[2], double mu, double rho, double d, double c
double F = (mass*GoU)/(mass+h*dt*GoU);
drag[0] = (u[0]-gradp[0]*dt/rhop)*F;
drag[1] = (u[1]-9.81*dt*(1-rho/rhop)-gradp[1]*dt/rhop)*F;
drag[1] = (u[1]-g*dt*(1-rho/rhop)-gradp[1]*dt/rhop)*F;
}
#else
static void particle_drag(double u[3], double mu, double rho, double d, double c, double drag[3], double dt, double mass, double rhop, double gradp[3])
static void particle_drag(double u[3], double mu, double rho, double d, double c, double drag[3], double dt, double mass, double rhop, double gradp[3], double g)
{
double normu = pow(u[0]*u[0]+u[1]*u[1]+u[2]*u[2],0.5);
//Reynolds/|u_p-u_f|
......@@ -288,7 +288,7 @@ static void particle_drag(double u[3], double mu, double rho, double d, double c
double F = (mass*GoU)/(mass+h*dt*GoU);
drag[0] = (u[0]-gradp[0]*dt/rhop)*F;
drag[1] = (u[1]-9.81*dt*(1-rho/rhop)-gradp[1]*dt/rhop)*F;
drag[1] = (u[1]-g*dt*(1-rho/rhop)-gradp[1]*dt/rhop)*F;
drag[2] = (u[2]-gradp[2]*dt/rhop)*F;
}
......@@ -352,7 +352,7 @@ void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt,
double rhop = mass/vol;
double d = 2*sqrt(vol/M_PI);
double drag[2];
particle_drag(du,problem->mu,problem->rho,d,c,drag, dt, mass, rhop, gradp);
particle_drag(du,problem->mu,problem->rho,d,c,drag, dt, mass, rhop, gradp,problem->g);
/*FILE *log = fopen("DragForce.txt", "at");
if (!log) log = fopen("DragForce.txt", "wt");
fprintf(log, "%d, %e, %e, %e, %e, %e, %e\n", i, du[0], du[1], drag[0], drag[1], mass, vol);
......@@ -444,7 +444,7 @@ void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt,
double rhop = mass/vol;
double d = 2*pow(3./4.*vol/M_PI,1./3.);
double drag[3];
particle_drag(du,problem->mu,problem->rho,d,c,drag, dt, mass, rhop,gradp);
particle_drag(du,problem->mu,problem->rho,d,c,drag, dt, mass, rhop,gradp,problem->g);
particle_force[i*3+0] += -drag[0]-vol*gradp[0];
particle_force[i*3+1] += -drag[1]-vol*gradp[1];
particle_force[i*3+2] += -drag[2]-vol*gradp[2];
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment