Commit 999580da authored by Matthieu Constant's avatar Matthieu Constant
Browse files

liquid drop with adapt mesh

parent d1a787fa
......@@ -40,31 +40,42 @@ ii = 0
#physical parameters
g = -9.81 # gravity
rho0 = 2000 # fluid density
rho1 = 1000
nu = 1e-4 # kinematic viscosity
mu = nu*rho0 # dynamic viscosity
tEnd = 10 # final time
g = -9.81
rhof = 1030 # fluid density
rhop = 2450 # grains density
r=154e-6 # grains radii
R = 3.3e-3 # drop radius
compacity = .20 # solid volume fraction in the drop
phi = 1 - compacity
nuf = 1.17/rhof # kinematic viscosity
muf = nuf*rhof # dynamic viscosity
##mixture properties
rhom = (1-phi)*rhop+phi*rhof #mixture density
num = nuf*(1+5/2*phi) #Einstein viscosity
print("rhom = %g, num = %g" % (rhom,num))
tEnd = 200 # final time
#numerical parameters
lcmin = .1 # mesh size
dt = .01 # time step
dt = 1e-2 # time step
shutil.copy("mesh.msh", outputdir +"/mesh.msh")
outf = 1 # number of iterations between output files
outf = 10 # number of iterations between output files
def outerBndV(x) :
print(.4*max(np.sin(t*np.pi*2./1),0))
return 0.4*max(np.sin(t*np.pi*2./1),0)
strong_boundaries = [("Top",2,2,0.)]#,("Lateral",0,0,0),("Bottom",1,1,0)]
fluid = fluid.fluid_problem("mesh.msh",g,[nu*rho0,nu*rho1],[rho0,rho1],strong_boundaries,1)
strong_boundaries = [("Top",2,2,0.),("Top",1,1,0.),("Bottom",1,1,0.),("Lateral",0,0,0.)]#,("Lateral",0,0,0),("Bottom",1,1,0)]
fluid = fluid.fluid_problem("mesh.msh",g,[num*rhom,nuf*rhof],[rhom,rhof],strong_boundaries,1)
fluid.set_weak_boundary("Bottom","Wall")
fluid.set_weak_boundary("Lateral","Wall")
fluid.set_weak_boundary("Top","Wall")
fluid.set_weak_boundary("Injection","Wall")
ii = 0
t = 0
......@@ -72,7 +83,7 @@ s = fluid.solution()
x = fluid.coordinates()
for i in range(len(x[:,0])):
if (x[i,0])**2+(x[i,1]-.14)**2<.01**2:
if (x[i,0])**2+(x[i,1]-.52)**2<R**2:
s[i,3] = 1
else:
s[i,3] = 0
......@@ -86,9 +97,9 @@ ii = 0
tic = time.clock()
while t < tEnd :
#Fluid solver
dt = min(dt+.1*dt,.1)
print(dt)
fluid.implicit_euler(dt)
if (ii%11==0 and ii != 0):
fluid.adapt_mesh(1e-3,1e-4,10000)
t += dt
#Output files writting
if ii %outf == 0 :
......
L = .1;
H = .2;
L = 0.05;
H = 0.6;
y = 0;
lc = 0.01;
f = 4;
Point(1) = {-L, H, 0,lc};
Point(2) = {-L, -H, 0,lc};
Point(3) = {-0.004,-H,0,lc/f};
Point(4) = {0.004,-H,0,lc/f};
Point(5) = {L, -H, 0,lc};
Point(6) = {L, H, 0,lc};
Point(9) = {0,-H,0,lc/f};
Point(10) = {0,H,0,lc/f};
lc = 0.001;
Point(1) = {-L, H, 0};
Point(2) = {-L, -H, 0};
Point(3) = {L, -H, 0};
Point(4) = {L, H, 0};
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 5};
Line(5) = {5, 6};
Line(6) = {6, 1};
Line(9) = {9,10};
Line(4) = {4, 1};
Point(5) = {0,0.52,0};
Line Loop(1) = {1:6};
Line Loop(1) = {1:4};
Plane Surface(1) = {1};
Physical Line("Bottom") = {2,4};
Physical Line("Lateral") = {1,5};
Physical Line("Top") = {6};
Physical Line("Injection") = {3};
Physical Line("Bottom") = {2};
Physical Line("Lateral") = {1,3};
Physical Line("Top") = {4};
Physical Surface("Domain") = {1};
Physical Point("PtFix") = {1};
Field[1] = Attractor;
Field[1].EdgesList = {9};
Field[1].NNodesByEdge = 200;
Field[1].NodesList = {5};
Field[2] = Threshold;
Field[2].DistMax = .05;
Field[2].DistMax = 0.08;
Field[2].DistMin = 0.02;
Field[2].LcMax = lc/f;
Field[2].LcMin = lc/f;
Field[2].LcMax = 0.008;
Field[2].LcMin = 0.001;
Field[2].IField = 1;
Background Field = 2;
Mesh.CharacteristicLengthFromPoints = 0;
Mesh.CharacteristicLengthExtendFromBoundary = 0;
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