2 fluids

src/fluid_problem.c

@@ -196,7 +196,9 @@ static void node_force_f(FluidProblem *problem, double *f0, double *sol, double
   double rhop = mass/vol;
   double drho = rhop-rho;
   double g = problem->g;
-  double gamma = particle_drag_coeff(Due,mu,rho,vol,c);
+  double gamma0 = particle_drag_coeff(Due,problem->mu[0],problem->rho[0],vol,c);
+  double gamma1 = particle_drag_coeff(Due,problem->mu[1],problem->rho[1],vol,c);
+  double gamma = gamma0*a + gamma1*(1-a);
   double fcoeff = mass/(mass+dt*gamma);
 
   for (int j = 0; j < D; ++j){
@@ -472,10 +474,10 @@ static void compute_stab_parameters(FluidProblem *problem, double dt) {
       mu = problem->mu[0]*a + problem->mu[1]*(1-a);
     }
     double nu = mu/rho;
-    problem->taup[iel] = 1/sqrt(pow2(2/dt)+pow2(2*normu/h)+pow2(4*nu/pow2(h)));
+    problem->taup[iel] = 1/sqrt(pow2(6/dt)+pow2(2*normu/h)+pow2(4*nu/pow2(h)));
     problem->tauc[iel] = h*normu*fmin(h*normu/(6*nu),0.5);
-    problem->taua[iel] = 1/sqrt(pow2(2/dt)+pow2(2*normu/h));
-    problem->extraa[iel] = fmax(0.,fmax(-amin,amax-1))*0.1;
+    problem->taua[iel] = 1/sqrt(pow2(6/dt)+pow2(2*normu/h));
+    problem->extraa[iel] = fmax(0.,fmax(-amin,amax-1))*0.5;
     maxtaup = fmax(maxtaup, problem->taup[iel]);
     mintaup = fmin(mintaup, problem->taup[iel]);
   }

testcases/injectionWGrains/inject2fGrains.py

@@ -58,7 +58,7 @@ def genInitialPosition(filename, N, r, lx, ly, rhop) :
 alpha = 0*np.pi/4.
 g =  -9.81*np.cos(alpha)                                      # gravity
 print(g)
-rho0 = 1.117*10                                      # fluid density
+rho0 = 1.117                                      # fluid density
 rho1 = 1000
 nu1 = 1.57e-5                                   # kinematic viscosity
 nu0 = 1e-6
@@ -70,7 +70,7 @@ ly = .14
 rhop = 2500
 #numerical parameters
 lcmin = .1                                  # mesh size
-dt = .001                                       # time step
+dt = .001                                      # time step
 
 shutil.copy("mesh.msh", outputdir +"/mesh.msh")
 genInitialPosition(outputdir, N, r, lx, ly, rhop)
@@ -81,7 +81,7 @@ outf = 1                                       # number of iterations between ou
 ii = 0
 t = 0
 def outerBndV(x) :
-    return max(0.1*t,0.1)
+    return max(1*t,0.1)
 
 strong_boundaries = [("Top",2,2,0.),("Injection",1,1,outerBndV),("Injection",3,3,1),("Injection",0,0,0)
 #        ("Top",3,3,1),("Lateral",3,3,1),("Bottom",3,3,1)

testcases/injectionWGrains/inject2fGrainsReload.py

@@ -44,7 +44,7 @@ ii = 0
 alpha = 0*np.pi/4.
 g =  -9.81*np.cos(alpha)                                      # gravity
 print(g)
-rho0 = 1.117*10                                      # fluid density
+rho0 = 1.117                                      # fluid density
 rho1 = 1000
 nu1 = 1.57e-5                                   # kinematic viscosity
 nu0 = 1e-6
@@ -57,7 +57,7 @@ rhop = 2500
 #numerical parameters
 lcmin = .1                                  # mesh size
 dt = .001                                       # time step
-iReload = 30
+iReload = 1270
 
 
 shutil.copy("mesh.msh", outputdir +"/mesh.msh")
@@ -70,7 +70,7 @@ t = 0
 
 
 def outerBndV(x) :
-    return max(0.1*t,0.1)
+    return max(t,0.1)
 
 strong_boundaries = [("Top",2,2,0.),("Injection",1,1,outerBndV),("Injection",3,3,1),("Injection",0,0,0)
 #        ("Top",3,3,1),("Lateral",3,3,1),("Bottom",3,3,1)