Commit a7b88f33 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

reduce diffusivity 2 fluids

parent d9fce575
Pipeline #4676 passed with stage
in 22 seconds
......@@ -506,7 +506,7 @@ static void compute_stab_parameters(FluidProblem *problem, double dt) {
problem->taup[iel] = 1/sqrt(pow2(2/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+0.01,amax-1+0.01))*problem->coeffStab;
problem->extraa[iel] = fmax(0.,fmax(-amin,amax-1))*problem->coeffStab;
}
}
......@@ -582,13 +582,12 @@ static void fluid_problem_f(FluidProblem *problem, double *f0, double *f1, doubl
f1[(U+i)*D+i] += (stab_param==0 ? 1 : 0)*divu*tauc*rho;
f1[P*D+i] = -u[i] + (stab_param==0 ? taup*R : stab_param*dp[i]);
if(problem->n_fluids == 2){
f1[P*D+i] += taua*(divu+(c-cold)/dt)*u[i];
f1[A*D+i] = -u[i]*a + taua*Ra*u[i] + taup*R*a + da[i]*(extraa);
f1[A*D+i] = -u[i]*a + 2*taua*Ra*u[i] +da[i]*extraa;
}
}
if(problem->n_fluids == 2){
/*if(problem->n_fluids == 2){
f1[A*D+1] += -fmax(a,0)*fmax((1-a),0)*1e-6;
}
}*/
}
static void compute_weak_boundary_conditions(FluidProblem *problem, double dt, double *all_local_vector, double *all_local_matrix)
......
......@@ -65,7 +65,7 @@ def _is_string(s) :
class FluidProblem :
"""Create numerical representation of the fluid."""
def __init__(self, dim, g, mu, rho, coeff_stab=0.01, petsc_solver_type="-pc-type lu"):
def __init__(self, dim, g, mu, rho, coeff_stab=0.01, petsc_solver_type="-pc_type lu -ksp_monitor"):
"""Build the fluid structure.
Keyword arguments:
......
......@@ -25,10 +25,10 @@ Line(9) = {9,10};
Line Loop(1) = {1:6};
Plane Surface(1) = {1};
Physical Line("Bottom") = {2,4};
Physical Line("Bottom") = {2,3,4};
Physical Line("Lateral") = {1,5};
Physical Line("Top") = {6};
Physical Line("Injection") = {3};
//Physical Line("Injection") = {3};
Physical Surface("Domain") = {1};
Physical Point("PtFix") = {1};
......
......@@ -51,31 +51,15 @@ outf = 1 # number of iterations between ou
lx = .111
def outerBndV(x) :
v = np.zeros_like(x)
if t < 1:
for i in range(len(x[:,1])):
if (x[:,0]-lx/2)**2<(8e-3)**2:
v[i,1] = 0.265258*(3*t**2-2*t**3)#(6*t**5-15*t**4+10*t**3)
else:
v[i,1] = 0
else:
for i in range(len(x[:,1])):
if ((x[:,0]-lx/2)**2<(8e-3)**2):
v[i,1] = 0.265258
else:
v[i,1] = 0
return v
x = 1-(x[:,0]-lx)**2/(2*8e-3)**2
return np.where(x>0, 0.265258*x*x, 0.)
fluid = fluid.FluidProblem(2,g,[nu0*rho0,nu1*rho1],[rho0,rho1])
fluid = fluid.FluidProblem(2,g,[nu0*rho0,nu1*rho1],[rho0,rho1],1e-2)
fluid.load_msh("mesh.msh")
fluid.set_strong_boundary("Top",2,0)
fluid.set_strong_boundary("Injection",1,outerBndV)
fluid.set_strong_boundary("Injection",3,1)
fluid.set_open_boundary("Bottom",velocity=outerBndV,concentration=1)
fluid.set_wall_boundary("Lateral","Wall")
fluid.set_open_boundary("Top",pressure=0,concentration=1)
fluid.set_open_boundary("Injection",velocity=outerBndV,concentration=1)
fluid.set_open_boundary("Bottom",velocity=[0,outerBndV],concentration=1)
fluid.set_open_boundary("Injection",velocity=[0,outerBndV],concentration=1)
fluid.set_wall_boundary("Lateral")
fluid.set_open_boundary("Top",pressure=0,concentration=0)
ii = 0
t = 0
......@@ -87,8 +71,10 @@ while t < tEnd :
fluid.implicit_euler(dt)
t += dt
#Output files writting
if ii%2 == 0:
fluid.adapt_mesh(4e-2,0.2e-2,7000)
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
fluid.export_vtk(outputdir, t, ioutput)
ii += 1
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))
\ No newline at end of file
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))
......@@ -25,10 +25,10 @@ Line(9) = {9,10};
Line Loop(1) = {1:6};
Plane Surface(1) = {1};
Physical Line("Bottom") = {2,4};
Physical Line("Bottom") = {2,3,4};
Physical Line("Lateral") = {1,5};
Physical Line("Top") = {6};
Physical Line("Injection") = {3};
//Physical Line("Injection") = {3};
Physical Surface("Domain") = {1};
Physical Point("PtFix") = {1};
......
Markdown is supported
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