New Reynolds

marblesbag/src/fluid_problem.c

@@ -44,16 +44,20 @@ static void particle_drag(double u[2], double mu, double rho, double d, double c
 {
   double normu = hypot(u[0],u[1]);
   //Reynolds/|u_p-u_f|
-  double Re_O_u = d*c/mu*rho;
+  double Re_O_u = sqrt(d)*c/mu*rho;
   double Cd_u = 0.63*sqrt(normu)+4.8/sqrt(Re_O_u);
   Cd_u = Cd_u*Cd_u;
   double f = pow(c,-1.8);
   double GoU = f*Cd_u*rho/2*d;
-  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));
+
+  double a = 0;
+  double h = (1+(1-c)*rhop/(c*rho)+a)/(a+1);
+  //printf("h = %.6g\n", h);
+  double F = (mass*GoU)/(mass+h*dt*GoU);
+  //double F = mass/(h*dt);
+
+  drag[0] = (u[0]-gradp[0]*dt/rhop)*F;
+  drag[1] = (u[1]-9.81*dt*(1-rho/rhop)-gradp[1]*dt/rhop)*F; 
 }
 
 static double mesh_get_min_edge_length(Mesh *mesh) {
@@ -71,21 +75,23 @@ static double mesh_get_min_edge_length(Mesh *mesh) {
   return sqrt(emin);
 }
 
-static void fluid_problem_smooth_velocity(FluidProblem *problem, double lsmooth, double *smooth_velocity) {
+static void fluid_problem_smooth_velocity(FluidProblem *problem, double lsmooth, double *smooth_velocity, double *smooth_c) {
   Mesh *mesh = problem->mesh;
   double emin = mesh_get_min_edge_length(mesh);
-  int n = ceil(lsmooth*lsmooth/(emin*emin));
+  int n = ceil(10*lsmooth*lsmooth/(emin*emin));
   double alpha = lsmooth*lsmooth/n;
-  double *rhs = malloc(mesh->n_nodes*2*mesh->n_nodes);
+  double *rhs = malloc(mesh->n_nodes*3*mesh->n_nodes);
   double *mass = malloc(mesh->n_nodes*mesh->n_nodes);
   for (int i = 0; i < mesh->n_nodes; ++i){
     smooth_velocity[i*2+0] = problem->solution[i*3+0];
     smooth_velocity[i*2+1] = problem->solution[i*3+1];
+    smooth_c[i] = problem->porosity[i];
   }
   for (int iter = 0; iter < n; ++iter) {
     for (int i = 0; i < mesh->n_nodes; ++i){
-      rhs[i*2+0] = 0;
-      rhs[i*2+1] = 0;
+      rhs[i*3+0] = 0;
+      rhs[i*3+1] = 0;
+      rhs[i*3+2] = 0;
       mass[i] = 0;
     }
     for (int iel=0; iel < mesh->n_triangles; ++iel) {
@@ -109,6 +115,7 @@ static void fluid_problem_smooth_velocity(FluidProblem *problem, double lsmooth,
       };
       double U[] = {smooth_velocity[tri[0]*2+0],smooth_velocity[tri[1]*2+0],smooth_velocity[tri[2]*2+0]};
       double V[] = {smooth_velocity[tri[0]*2+1],smooth_velocity[tri[1]*2+1],smooth_velocity[tri[2]*2+1]};
+      double C[] = {smooth_c[tri[0]],smooth_c[tri[1]],smooth_c[tri[2]]};
       double du[2] = {
         dphi[0][0]*U[0]+dphi[1][0]*U[1]+dphi[2][0]*U[2],
         dphi[0][1]*U[0]+dphi[1][1]*U[1]+dphi[2][1]*U[2]
@@ -117,15 +124,21 @@ static void fluid_problem_smooth_velocity(FluidProblem *problem, double lsmooth,
         dphi[0][0]*V[0]+dphi[1][0]*V[1]+dphi[2][0]*V[2],
         dphi[0][1]*V[0]+dphi[1][1]*V[1]+dphi[2][1]*V[2]
       };
+      double dc[2] = {
+        dphi[0][0]*C[0]+dphi[1][0]*C[1]+dphi[2][0]*C[2],
+        dphi[0][1]*C[0]+dphi[1][1]*C[1]+dphi[2][1]*C[2]
+      };
       for (int i = 0; i < 3; ++i) {
-        rhs[tri[i]*2+0] -= (dphi[i][0]*du[0] + dphi[i][1]*du[1])*detj/2;
-        rhs[tri[i]*2+1] -= (dphi[i][0]*dv[0] + dphi[i][1]*dv[1])*detj/2;
+        rhs[tri[i]*3+0] -= (dphi[i][0]*du[0] + dphi[i][1]*du[1])*detj/2;
+        rhs[tri[i]*3+1] -= (dphi[i][0]*dv[0] + dphi[i][1]*dv[1])*detj/2;
+        rhs[tri[i]*3+2] -= 0*(dphi[i][0]*dc[0] + dphi[i][1]*dc[1])*detj/2;
         mass[tri[i]] += detj/6;
       }
     }
     for (int i = 0; i < mesh->n_nodes; ++i) {
-      smooth_velocity[i*2+0] += rhs[i*2+0]/mass[i]*alpha;
-      smooth_velocity[i*2+1] += rhs[i*2+1]/mass[i]*alpha;
+      smooth_velocity[i*2+0] += rhs[i*3+0]/mass[i]*alpha;
+      smooth_velocity[i*2+1] += rhs[i*3+1]/mass[i]*alpha;
+      smooth_c[i] += rhs[i*3+2]/mass[i]*alpha;
     }
   }
   free(mass);
@@ -136,7 +149,8 @@ void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt,
   double *porosity = problem->porosity;
   Mesh *mesh = problem->mesh;
   double *smoothv = malloc(mesh->n_nodes*2*sizeof(double));
-  fluid_problem_smooth_velocity(problem, lsmooth, smoothv);
+  double *smoothc = malloc(mesh->n_nodes*sizeof(double));
+  fluid_problem_smooth_velocity(problem, lsmooth, smoothv,smoothc);
   const double *solution = problem->solution;
   for (int i = 0; i < mesh->n_nodes*2; ++i) {
     problem->node_force[i] = 0.0;
@@ -153,7 +167,7 @@ void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt,
       continue;
     }
     int *tri = problem->mesh->triangles+iel*3;
-    double C[] = {porosity[tri[0]],porosity[tri[1]],porosity[tri[2]]};
+    double C[] = {smoothc[tri[0]],smoothc[tri[1]],smoothc[tri[2]]};
     double U[] = {smoothv[tri[0]*2+0],smoothv[tri[1]*2+0],smoothv[tri[2]*2+0]};
     double V[] = {smoothv[tri[0]*2+1],smoothv[tri[1]*2+1],smoothv[tri[2]*2+1]};
     double P[] = {solution[tri[0]*3+2],solution[tri[1]*3+2],solution[tri[2]*3+2]};