change Patankar

marblesbag/src/fluid_problem.c

@@ -40,7 +40,7 @@ int fluid_problem_export(const FluidProblem *problem, const char *output_dir, do
   return 0;
 }
 
-static void particle_drag(double u[2], double mu, double rho, double d, double c, double drag[2], double dt, double mass, double rhop)
+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 normu = hypot(u[0],u[1]);
   //Reynolds/|u_p-u_f|
@@ -49,8 +49,11 @@ static void particle_drag(double u[2], double mu, double rho, double d, double c
   Cd_u = Cd_u*Cd_u;
   double f = pow(c,-1.8);
   double GoU = f*Cd_u*rho/2*d;
-  drag[0] = GoU*u[0]/(1+(rhop+c/(1-c)*rho)*dt/(c/(1-c)*rho*mass)*GoU);
-  drag[1] = GoU*u[1]/(1+(rhop+c/(1-c)*rho)*dt/(c/(1-c)*rho*mass)*GoU);
+  double h = 1;//1+(1-c)*rhop/(c*rho);
+  drag[0] = GoU*(rhop*u[0])/(1+(rhop+c/(1-c)*rho)*dt/(c/(1-c)*rho*mass)*GoU);
+  drag[1] = GoU*(rhop*u[1]-9.81*(rhop-rho)*dt)/(rhop+(rhop*(1-c)+c*rho)*dt*GoU/(c*rho*mass/rhop));
+  drag[1] = GoU*(rhop*u[1]-9.81*dt*(rhop-rho)-gradp[1]*dt)/(rhop+dt*GoU*h/(mass/rhop));
+  drag[0] = GoU*(rhop*u[0]-gradp[0]*dt)/(rhop+dt*GoU*h/(mass/rhop));
 }
 
 void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt, double *particle_force) {
@@ -110,9 +113,15 @@ 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);
+    particle_drag(du,problem->mu,problem->rho,d,c,drag, dt, mass, rhop, gradp);
+    /*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);
+    fclose(log);*/
+
     particle_force[i*2+0] += -drag[0]-vol*gradp[0];
     particle_force[i*2+1] += -drag[1]-vol*gradp[1];
+    //printf("%e,%e,%e,%e,%e,%e\n",particle_force[i*2+0],particle_force[i*2+1],-drag[0],-drag[1],0.,problem->g*(mass-problem->rho*vol));
     for (int iphi = 0; iphi < 3; ++iphi) {
       problem->node_force[tri[iphi]*2+0] += drag[0]*phi[iphi];
       problem->node_force[tri[iphi]*2+1] += drag[1]*phi[iphi];