Commit 375d3b07 authored by Maud Schauss's avatar Maud Schauss
Browse files

depot sable

parent d53ab791
Pipeline #7563 failed with stage
in 18 seconds
from migflow import fluid
from migflow import scontact
from migflow import time_integration
import subprocess
import numpy as np
import os
import time
import shutil
import random
def genInitialPosition(filename, r, H, ly, lx, rhop) :
#Particles structure builder
p = scontact.ParticleProblem(2,True,True)
#Load mesh.msh file specifying physical boundaries names
p.load_msh_boundaries("mesh.msh", ["Top_Left","Top_Right","Lateral","Bottom","Pile_Left","Pile_Right","Pile_Pointe"], material="Pieu")
#Definition of the points where the particles are located
x = np.arange(-lx/2+r-1e-8, lx/2-r+1e-8, 2*r)
y = np.arange(-H/2+r, -H/2+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)) :
if abs(x[i]) > 0.04 :
r2 = (1-random.random()*0.3)*r
p.add_particle((x[i], y[i]), r2, r2**2*np.pi*rhop,"Sable")
continue
if y[i]<-0.0375 and abs(x[i]) < 0.04 :
r2 = (1-random.random()*0.3)*r
p.add_particle((x[i], y[i]), r2, r2**2*np.pi*rhop,"Sable")
p.write_vtk(filename,0,0)
#Define output directory
outputdir = "output"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
#Physical parameters
g = -9.81 # gravity
rho = 1000 # fluid density
nu = 1e-6 # kinematic viscosity
mu= 1 # dynamic viscosity
#Numerical parameters
outf = 50 # number of iterations between output files
dt = 5e-4 # time step
tEnd = 4 # final time
#Geometrical parameters
H=0.4 # domain height
L=0.8 # domain lenght
L_pile=0.05 # pile width
#Parameters of pile movement
A=H/4 # amplitude
w=2*np.pi # frequency
#Particle parameters
r = 4e-3 # particles radius
ly = H # particles area height
lx = L # particles area width
rhop = 1500 # particles density
#Initial time and iteration
t = 0
ii = 0
iter = 0
#
# PARTICLE PROBLEM
#
# Initialise particles
genInitialPosition(outputdir, r, H, ly, lx, rhop)
p = scontact.ParticleProblem(2,True,True)
p.read_vtk(outputdir,0)
p.set_friction_coefficient(0.1, "Pieu", "Sable")
p.set_friction_coefficient(0.3, "Sable", "Sable")
G = np.zeros_like(p.velocity())
G[:,1] = p.mass()[:,0]*g
tic = time.time()
#
# COMPUTATION LOOP
#
while t < tEnd:
time_integration.iterate(None, p, dt, min_nsub=5, external_particles_forces=G)
#time_integration.iterate(fluid, None, dt)
t += dt
# Output files writting
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
p.write_vtk(outputdir, ioutput, t)
iter += 1
ii += 1
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))
// Gmsh project created on Mon Sep 23 22:59:42 2019
SetFactory("OpenCASCADE");
L = 0.4;
H = 0.2;
l=0.025;
h=0.15;
lc = 0.05;
Point(1) = {-L, H, 0,lc/2};
Point(2) = {-L, -H, 0,lc};
Point(3) = {L, -H, 0, lc};
Point(4) = {L, H, 0, lc/2};
Point(11)= {-l,H, 0, lc/2};
Point(12)= {l, H, 0, lc/2};
//Point(5) = {-l, 2*h, 0,lc/2};
Point(6) = {-l, 1/4*h, 0,lc/2};
Point(7) = {l, 1/4*h, 0,lc/2};
//Point(8) = {l, 2*h, 0,lc/2};
Point(9) = {0, 0, 0,lc/2};
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 12};
Line(8) = {11,1};
Line(9) = {11, 6};
Line(10) = {6, 9};
Line(11) = {9, 7};
Line(12) = {7, 12};
Line Loop(1) = {1,2,3,4,12,11,10,9,8};
Plane Surface(1) = {1};
Physical Line("Bottom") = {2};
Physical Line("Lateral") = {3,1};
Physical Line("Pile_Left") = {9};
Physical Line("Pile_Right") = {12};
Physical Line("Pile_Pointe") = {10,11};
Physical Line("Top_Right") = {4};
Physical Line("Top_Left") = {8};
Physical Surface("Domain") = {1};
Physical Point("PtFix") = {1};
Mesh.MshFileVersion = 2;
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