Commit 801a8b81 authored by Margaux Boxho's avatar Margaux Boxho
Browse files

Problem with boundaries

parent ec7d6a36
Pipeline #5385 passed with stage
in 26 seconds
This diff is collapsed.
L = 0.2;
L = 0.05;
H = 0.6;
y = 0;
lc = 0.05;
lc_min = 0.01;
lc = 40*1e-3;
lc_min = 20*1e-3;
// Defining the rotation
alpha = 30*Pi/180;
a1 = -L; b1 = -H;
a2 = L; b2 = -H;
a3 = L; b3 = H;
a4 = -L; b4 = H;
a5 = -L; b5 = 0;
a6 = L; b6 = 0;
a1 = -L; b1 = -2*H;
a2 = L; b2 = -2*H;
a3 = L; b3 = 0;
a4 = -L; b4 = 0;
a5 = -L; b5 = -H;
a6 = L; b6 = -H;
s = Sin(alpha);
c = Cos(alpha);
mx = 0; my = 0.52;
mx = 0; my = -0.52;
Point(1) = {c*a4-s*b4, s*a4+c*b4, 0,lc};
Point(2) = {c*a1-s*b1, s*a1+c*b1, 0,lc};
Point(3) = {c*a2-s*b2, s*a2+c*b2, 0,lc};
Point(4) = {c*a3-s*b3, s*a3+c*b3, 0,lc};
Point(5) = {c*a5-s*b5, s*a5+c*b5, 0,lc};
Point(6) = {c*a6-s*b6, s*a6+c*b6, 0,lc};
Point(1) = {c*a4-s*b4, s*a4+c*b4, 0,lc_min};
Point(2) = {c*a1-s*b1, s*a1+c*b1, 0,lc_min}; //
Point(3) = {c*a2-s*b2, s*a2+c*b2, 0,lc_min}; //
Point(4) = {c*a3-s*b3, s*a3+c*b3, 0,lc_min};
Point(5) = {c*a5-s*b5, s*a5+c*b5, 0,lc_min};
Point(6) = {c*a6-s*b6, s*a6+c*b6, 0,lc_min};
Line(1) = {1, 5}; // Lateral left
Line(2) = {5, 2}; // Lateral left 2
Line(3) = {2, 3}; // Bottom
......@@ -41,7 +41,7 @@ Physical Line("Lateral") = {1,2, 4, 5};
Physical Line("Top") = {6};
Physical Surface("Domain") = {1};
Physical Point("PtFix") = {1};
Mesh.Algorithm = 5;
Mesh.Algorithm = 8;
Merge "lc.pos";
Field[1] = PostView;
......
......@@ -66,14 +66,14 @@ def genInitialPosition(filename, r, H, ly, lx, rhop, compacity) :
#x = x.flat
#y = y.flat
N = compacity*(lx*ly)/(np.pi*r**2)
N = 4*compacity*(lx*ly)/(np.pi*r**2)
e = 2*ly/(N)**(1./2.)
# Definition of the points where the particles are located
shift = 0.1; alpha = 30*np.pi/180;
lx = lx*np.cos(alpha)
x = np.arange(lx/2+50*r, -lx/2-50*r, -e)
x = np.arange(lx/2+25*r, -lx/2-25*r, -e)
y = np.arange(-3*r-shift, -ly+3*r-shift, -e)
#shift = 0.26;
......@@ -104,20 +104,20 @@ if not os.path.isdir(outputdir) :
# Physical parameters
g = -9.81 # gravity
r = 1e-3 # particles radius
rhop = 1500 # particles density
rhop = 2200 # particles density as aluminium
rho = 1000 # fluid density
nu = 1e-6 # kinematic viscosity
mu = rho*nu;
compacity = 0.1
compacity = 0.7
# Numerical parameters
outf = 5 # number of iterations between output files
dt = 1e-3 # time step
tEnd = 2 # final time
outf = 5 # number of iterations between output files
dt = 0.1e-2 # time step
tEnd = 10 # final time
# Geometrical parameters
ly = 0.15 # particles area height
lx = 4e-1 # particles area widht
ly = 0.8 # particles area height
lx = 0.05 # particles area widht
H = 0.6 # domain height
#
......@@ -135,8 +135,8 @@ U_0 = 2*9.81*(rhop-rho)*r**2/(9*mu);
phi = compacity
f_phi = ((1-phi)**2)/((1+phi**(1/3))*np.exp(5/3*phi/(1-phi)));
U_0 = U_0*f_phi;
P_0 = (U_0**2)*0.5*rho;
U_0 = 0.065;
P_0 = 130.755;
Re = rho*U_0*r/mu;
N_0 = len(p.position());
e_target = 1e-7
......@@ -171,14 +171,39 @@ fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(), p.contact_
# COMPUTATION LOOP
#
while t < tEnd :
if (ii==jj*delta):
print("Save state at %i" %(ii))
fluid.export_vtk(outputdir,0,0)
p.write_vtk(outputdir,0,0)
# Fluid solver
fluid.implicit_euler(dt)
# Adaptation of the mesh.
# No adapt for this reference case !!!
#if (ii%15==0 and ii != 0):
# fluid.adapt_mesh(500*r,8*r,6602*2) # for the moment we cannot set the number of elements in the new mesh
if (ii==0):
fluid.adapt_mesh(60*r, 5*r, 6602*2, e_target, 1, U_0, P_0)
if (ii%(delta*(jj+1))!=0):
fluid.adapt_mesh(60*r, 5*r, 6602*2, e_target, 0, U_0, P_0)
# Adaptation of the mesh and loading the past save state
if (ii%(delta*(jj+1))==0 and ii != 0):
print("Adapt mesh at %i" %(ii))
# Adapt the mesh by creating the lc.pos file
fluid.adapt_gen_mesh(60*r, 5*r, 6602*2, e_target, U_0, P_0)
ii = ii-delta
t = t-dt*delta
jj += 1
print("ii=%i, t=%i, jj=%i" %(ii,t,jj))
# load the save past state
fluid.import_vtk("output/fluid_00000.vtu")
p.read_vtk(outputdir,0)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(), p.contact_forces())
# create the .msh on the adapt/lc.pos
fluid.adapt_load_mesh("adapt/mesh.msh")
#fluid.implicit_euler(dt)
print("Restart computation %i iterations back at %i" %(delta, ii))
# Initialize the min_h_target attribut of the Mesh structure
fluid.adapt_mesh(60*r, 5*r, 6602*2, e_target, 1, U_0, P_0)
# Computation of the new velocities
forces = fluid.compute_node_force(dt)
......
L = 0.2;
L = 0.05;
H = 0.6;
y = 0;
lc = 60*1e-3;
lc = 40*1e-3;
lc_min = 20*1e-3;
// Defining the rotation
......@@ -41,6 +41,6 @@ Physical Line("Lateral") = {1,2, 4, 5};
Physical Line("Top") = {6};
Physical Surface("Domain") = {1};
Physical Point("PtFix") = {1};
Mesh.Algorithm = 5;
Mesh.Algorithm = 8;
Mesh.MshFileVersion = 2;
......@@ -34,7 +34,7 @@ import time
import shutil
import random
def genInitialPosition(filename, r, H, ly, lx, rhop) :
def genInitialPosition(filename, r, H, ly, lx, rhop, compacity) :
"""Set all the particles centre positions and create the particles objects to add in the computing structure
Keyword arguments:
......@@ -51,11 +51,29 @@ def genInitialPosition(filename, r, H, ly, lx, rhop) :
p.load_msh_boundaries("mesh.msh", ["Top", "Lateral","Bottom"])
#Definition of the points where the particles are located
x = np.arange(-lx/2+r, lx/2-r, 2*r)
y = np.arange(H/2-r, H/2-ly+r, -2*r)
#x = np.arange(-lx/2+r, lx/2-r, 2*r)
#y = np.arange(0-r, -ly+r, -2*r)
#x, y = np.meshgrid(x, y)
#x = x.flat
#y = y.flat
#N = compacity*(lx*ly)/(np.pi*r**2)
#e = 2*ly/(N)**(1./2.)
#print(N)
# Definition of the points where the particles are located
#x = np.arange(lx/2-r, -lx/2+r, -e)
#y = np.arange(-r, -ly+r, -e)
#x, y = np.meshgrid(x, y)
x = np.arange(lx/2-r, -lx/2+r, -2*r)
y = np.arange(-r, -ly+r, -2*r)
x, y = np.meshgrid(x, y)
x = x.flat
y = y.flat
# Add a grain at each centre position
for i in range(len(x)) :
p.add_particle((x[i], y[i]), r, r**2 * np.pi * rhop)
......@@ -72,11 +90,11 @@ r = 1e-3 # particles radius
rhop = 1500 # particles density
rho = 1000 # fluid density
nu = 1e-6 # kinematic viscosity
compacity = 0.25
# Numerical parameters
outf = 5 # number of iterations between output files
dt = 0.25e-2 # time step
tEnd = 100 # final time
dt = 0.1e-2 # time step
tEnd = 2 # final time
# Geometrical parameters
ly = 5e-2 # particles area height
......@@ -87,7 +105,7 @@ H = 1.2 # domain height
# PARTICLE PROBLEM
#
# Initialise particles
genInitialPosition(outputdir, r, H, ly, lx, rhop)
genInitialPosition(outputdir, r, H, ly, lx, rhop, compacity)
p = scontact.ParticleProblem(2)
p.read_vtk(outputdir,0)
......@@ -104,9 +122,9 @@ ii = 0
fluid = fluid.FluidProblem(2,g,[nu*rho],[rho])
# Set the mesh geometry for the fluid computation
fluid.load_msh("mesh.msh")
fluid.set_wall_boundary("Bottom")
fluid.set_wall_boundary("Lateral")
fluid.set_wall_boundary("Top",pressure=0)
fluid.set_wall_boundary("Bottom", velocity=[0,0])
fluid.set_wall_boundary("Lateral", velocity=[0,0])
fluid.set_wall_boundary("Top",pressure=0,velocity=[0,0])
# Set location of the particles in the mesh and compute the porosity in each computation cell
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.contact_forces())
fluid.export_vtk(outputdir,0,0)
......
L = 0.2;
H = 0.6;
y = 0;
lc = 0.01;
lc = 40*1e-3;
lc_min = 8*1e-3;
Point(1) = {-L, H, 0,lc};
Point(2) = {-L, -H, 0,lc};
Point(3) = {L, -H, 0,lc};
Point(4) = {L, H, 0,lc};
Point(1) = {-L, 0, 0,lc_min};
Point(2) = {-L, -2*H, 0,lc_min};
Point(3) = {L, -2*H, 0,lc_min};
Point(4) = {L, 0, 0,lc_min};
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 1};
Point(5) = {0,0.52,0};
Point(5) = {0,-0.52,0};
Line Loop(1) = {1:4};
......@@ -21,5 +22,6 @@ Physical Line("Lateral") = {1,3};
Physical Line("Top") = {4};
Physical Surface("Domain") = {1};
Physical Point("PtFix") = {1};
Mesh.Algorithm = 8;
Mesh.MshFileVersion = 2;
Supports Markdown
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