Commit 20f6c880 authored by Nathan Coppin's avatar Nathan Coppin
Browse files

is this the real life?

parent 6e4da62f
Pipeline #5719 passed with stage
in 2 minutes and 11 seconds
This diff is collapsed.
......@@ -44,24 +44,27 @@ def genInitialPosition(filename, rhop) :
# Particles structure builder
p = scontact.ParticleProblem(3,friction_enabled=True, rotation_enabled=True)
# Load the mesh.msh file specifying physical boundaries names
p.load_msh_boundaries("mesh.msh", ["Top", "Bottom", "X", "Z"])
r = 2.1*0.005
p.load_msh_boundaries("mesh2.msh", ["Top", "Bottom", "X", "Z"],material = "small")
r = 0.1
i = 0
# Positioning the particles on a regular grid
#for y in np.arange(.0-r, -0.1+r, -2*r):
# for z in np.arange(.0-r, -.1+r, -2*r):
# for x in np.arange(.0-r, -.1+r, -2*r):
# r1 = random.choice([r,r/2,r/2,r/2,r/2,r/2])
# p.add_particle((x, y, z), r1, 4./3.* r1**3 * np.pi * rhop,"small");
# i += 1;
p.add_particle((-0.05+r, -0.025, -0.050), r, 4./3.* r**3 * np.pi * rhop,"small");
p.add_particle((-0.05, -0.07, -0.050), r, 4./3.* r**3 * np.pi * rhop,"small");
for y in np.arange(-4.3-r, -10, -2*r*1.05):
for z in np.arange(.0-1.1*r, -1.+1*r, -2*r*1.05):
for x in np.arange(.0-1.1*r, -2.+r, -2*r*1.05):
r1 = np.random.uniform(.95*r, 1.05*r)
p.add_particle((x, y, z), r1, 4./3.* r1**3 * np.pi * rhop,"small");
i += 1;
#p.position()[:,2] -= r
p.position()[:,1] += 0.4
p.velocity()[:,2] += 0.5
#p.add_particle((-0.05+r, -0.025, -0.050), r, 4./3.* r**3 * np.pi * rhop,"small");
#p.add_particle((-0.05, -0.07, -0.050), r, 4./3.* r**3 * np.pi * rhop,"small");
print("Number of particles=%d" % len(p.position()[:, 1]))
p.write_vtk(filename, 0, 0)
# Define output directory
outputdir = "Test"
outputdir = "aval3"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
......@@ -71,8 +74,8 @@ rhop = 1500 # particles density
# Numerical parameters
outf = 2 # number of iterations between output files
dt = 1e-3 # time step
tEnd = 4 # final time
dt = 5e-4 # time step
tEnd = 10 # final time
#
# PARTICLE PROBLEM
......@@ -81,36 +84,40 @@ tEnd = 4 # final time
genInitialPosition(outputdir, rhop)
p = scontact.ParticleProblem(3,friction_enabled=True, rotation_enabled=True)
p.read_vtk("Test",0)
p.set_friction_coefficient(0.25,"small","small")
# Initial time and iteration
p.read_vtk(outputdir,0)
p.set_friction_coefficient(0.,"small","small")
# Initial time and iteration
t = 0
ii = 0
tic = time.time()
#
# COMPUTATION LOOP
#
keep = p.particle_material() == 1
p.velocity()[0,1] = -0.1
opened = 0
forces = np.zeros_like(p.velocity())
print(keep)
while t < tEnd :
#forces[keep,1] = g*p.mass()[keep,0]
if t>1.5 and opened == 0:
p.clear_boundaries()
p.load_msh_boundaries("mesh.msh", ["Top", "Bottom", "X", "Z"],material = "small")
p.velocity()[:,:] = 0
p.omega()[:,:] = 0
p.contact_forces()[:,:] = 0
p.set_friction_coefficient(0.41,"small","small")
opened=1
forces[:,1] = g*p.mass()[:,0]
vn = p.velocity() + forces*dt/p.mass()
vmax = np.max(np.linalg.norm(vn))
# number of contact sub time step
nsub = max(1, int(np.ceil((vmax * dt * 4)/min(p.r()))))
nsub = max(150, int(np.ceil((vmax * dt * 4)/min(p.r()))))
print("NSUB", nsub,"VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(p.r()))
# NLGS iterations
for i in range(nsub) :
p.iterate(dt/nsub, forces)
p.iterate(dt/nsub, forces, tol = 5e-7)
t += dt
# Output files writing
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
ioutput = int(ii/outf) + 1 + 2000
p.write_vtk(outputdir, ioutput, t)
ii += 1
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))
SetFactory("OpenCASCADE");
L = 6;
H = 10;
P = 1;
y = 0;
lc= 2e-3;
SetFactory("OpenCASCADE");
Box(1) = {-L,-H,-P,L,H,P};
Physical Volume("Domain") = {1};
Physical Surface("Top") = {4};
Physical Surface("Bottom") = {3};
Physical Surface("X") = {1,2};
Physical Surface("Z") = {5,6};
Physical Point("PtFix") = {5};
Characteristic Length{ PointsOf{ Volume{:}; } } = lc;
Mesh.MshFileVersion = 2;
SetFactory("OpenCASCADE");
L = 2;
H = 10;
P = 1;
y = 0;
lc= 2e-3;
SetFactory("OpenCASCADE");
Box(1) = {-L,-H,-P,L,H,P};
Physical Volume("Domain") = {1};
Physical Surface("Top") = {4};
Physical Surface("Bottom") = {3};
Physical Surface("X") = {1,2};
Physical Surface("Z") = {5,6};
Physical Point("PtFix") = {5};
Characteristic Length{ PointsOf{ Volume{:}; } } = lc;
Mesh.MshFileVersion = 2;
# Migflow - Copyright (C) <2010-2018>
# <Universite catholique de Louvain (UCL), Belgium
# Universite de Montpellier, France>
#
# List of the contributors to the development of Marblesbag: see AUTHORS file.
# Description and complete License: see LICENSE file.
#
# This program (Migflow) is free software:
# you can redistribute it and/or modify it under the terms of the GNU Lesser General
# Public License as published by the Free Software Foundation, either version
# 3 of the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with this program (see COPYING and COPYING.LESSER files). If not,
# see <http://www.gnu.org/licenses/>.
#!/usr/bin/env python
from migflow import scontact
import numpy as np
import os
import time
import shutil
import random
from math import sin,cos,radians
outputdir = "output"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
#grains creation and placement + physical boundaries definition
'''r is the radius of the grains
rout is the outer radius of the geometry
rin is the inner radius of the geometry
rhop is the particles density'''
def genInitialPosition(filename, r, rout, rhop) :
p = scontact.ParticleProblem(3,friction_enabled=True,rotation_enabled=True)
p.add_boundary_triangle((0,0,0),(rout,0,0),(rout,rout,0),"Bottom",material="gold")
p.add_boundary_triangle((0,0,0),(0,rout,0),(rout,rout,0),"Bottom",material="gold")
p.add_boundary_triangle((0,0,0),(rout,0,0),(rout,0,rout),"Front",material="gold")
p.add_boundary_triangle((0,0,0),(0,0,rout),(rout,0,rout),"Front",material="gold")
p.add_boundary_triangle((0,rout,0),(0,rout,rout),(rout,rout,rout),"Rear",material="gold")
p.add_boundary_triangle((0,rout,0),(rout,rout,0),(rout,rout,rout),"Rear",material="gold")
p.add_boundary_triangle((0,0,0),(0,rout,0),(0,rout,rout),"Left",material="gold")
p.add_boundary_triangle((0,0,0),(0,0,rout),(0,rout,rout),"Left",material="gold")
p.add_boundary_triangle((rout,0,0),(rout,rout,0),(rout,rout,rout),"Right",material="gold")
p.add_boundary_triangle((rout,0,0),(rout,0,rout),(rout,rout,rout),"Right",material="gold")
p.add_boundary_disk((0,0,0),0,"Front",material="gold")
p.add_boundary_segment((0,0,0), (0,rout,rout), "Left", material = "gold")
p.add_particle((0+1.5*r,0.5,0+r), r, r**2 * np.pi * rhop,"wut");
p.write_vtk(filename,0,0)
t = 0
ii = 0
#physical parameters
g = -9.81 # gravity
rhop = 2450 # grains density
tEnd = 0.20 # final time
#numerical parameters
dt = 1e-3 # time step
#geometry parameters
rout = 1 # outer radius
r = 10e-3 # grains radius
shutil.copy("mesh.msh", outputdir +"/mesh.msh")
#Object particles creation
genInitialPosition(outputdir, r, rout, rhop)
p = scontact.ParticleProblem(3,friction_enabled=True,rotation_enabled=True)
p.read_vtk(outputdir,0)
outf = 1 # number of iterations between output files
#Taking friction into account
p.set_friction_coefficient(5.8,"wut","gold")
#Computation loop
forces = np.zeros_like(p.velocity())
k = 0
while t < tEnd :
forces[:,2] = g*p.mass()[0]*cos(radians(25))
forces[:,0] = -g*p.mass()[0]*sin(radians(25))*(2**0.5/2)
forces[:,1] = -g*p.mass()[0]*sin(radians(25))*(2**0.5/2)
#Computation of the new velocities
vn = p.velocity() + forces * dt / p.mass()
vmax = np.max(np.sqrt(vn[:, 0]**2+ vn[:, 1]**2 + vn[:,2]**2))
#number of sub time step
nsub = max(1, int(np.ceil((vmax * dt * 4)/min(p.r()))))
print("i", k, "NSUB", nsub,"VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(p.r()))
k +=1
#Contact solver
for i in range(nsub) :
tol = 1e-12
p.iterate(dt/nsub, forces, tol)
t += dt
#Output files writting
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
p.write_vtk(outputdir, ioutput, t)
ii += 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