Commit bb625247 authored by Matthieu Constant's avatar Matthieu Constant
Browse files

cas test

parent df93a42a
Pipeline #5534 passed with stage
in 2 minutes and 20 seconds
# MigFlow - Copyright (C) <2010-2018>
# <Universite catholique de Louvain (UCL), Belgium
# Universite de Montpellier, France>
#
# List of the contributors to the development of MigFlow: 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 fluid
from migflow import scontact
from migflow import lmgc90Interface
from pylmgc90 import pre
import numpy as np
import os
import time
import shutil
import random
import sys
outputdir = "outputv5e-5grainshg4.5-air-volumedrag15-deb.013l-forceinstabtry-testsigma-07smallolddragverysmall"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
t = 0
ii = 0
use_lmgc = False
no_particle = False
def genInitialPosition(filename, N, r, lx, ly, rhop) :
p = scontact.ParticleProblem(2,friction_enabled=True)
#Loading of the mesh.msh file specifying physical boundaries name
p.load_msh_boundaries("mesh.msh", ["Top", "Bottom", "Lateral","Injection"])
#Definition of the points where the grains are located
Ra = 0.9*r
Rb = 1.1*r
radii = pre.granulo_Random(N, Ra, Rb)
[nb_laid_particles, coors] = pre.depositInBox2D(radii,lx,ly)
for i in range(nb_laid_particles):
p.add_particle(coors[2*i:2*i+2], radii[i], 4./3.*radii[i]**2 * np.pi * rhop);
p.write_vtk(filename,0,0)
#physical parameters
alpha = 0*np.pi/4.
g = -9.81*np.cos(alpha) # gravity
rho0 = 1.117 # fluid density
rho1 = 785.92
nu0 = 1.57e-5#1.57e-5 # kinematic viscosity
nu1 = 1.2e-3/rho1
tEnd = 50 # final time
r = 5e-4/2
N = 100000
lx = .02
ly = .025
rhop = 1059
#numerical parameters
dt = .001/4 # time step
H = .05
allp = []
class Injector :
def __init__(self, V, M, P) :
self._V = V
self._M0 = M
self._M = M
self._P = P
self._P0 = P
print("intialize injection : M = %g, P = %g" %(self._M, self._P))
@property
def P(self) :
return self._P
def iterate(self, QM) :
self._M += QM
self._P = self._M/self._M0 * self._P0
print("updating injection : QM = %g M = %g, P = %g" %(QM, self._M,self._P))
return self._P
Q=0.013e-3/60
Vin = 5e-5#todo callibrate
P0in = -rho1*(H)*g+1e5
injector = Injector(Vin,Vin*rho0,P0in)
w = 2e-3
drag = 15/(w*w)
shutil.copy("mesh.msh", outputdir +"/mesh.msh")
if not no_particle :
genInitialPosition(outputdir, N, r, lx, ly, rhop)
if use_lmgc:
friction=0.3
lmgc90Interface.scontactTolmgc90(outputdir, 2, 0, friction)
p = lmgc90Interface.ParticleProblem(2)
p.write_vtk(outputdir,0,0)
else:
p = scontact.ParticleProblem(2,friction_enabled=True)
p.read_vtk(outputdir,0)
#p = scontact.ParticleProblem(2,friction_enabled=True)
#for i in range(len(p1.mass())):
# p.add_particle(p1.position()[i,:],p1.r()[i],p1.mass()[i])
#p.load_msh_boundaries("mesh.msh", ["Top", "Bottom", "Lateral","Injection"])
p.set_friction_coefficient(.4)
p.write_vtk(outputdir,0,0)
print(p.r())
print(p.mass())
fluid = fluid.FluidProblem(2,g,[nu0*rho0,nu1*rho1],[rho0,rho1],volume_drag=0,sigma=0.007,drag_in_stab=1);
fluid.load_msh("mesh.msh")
l = 2e-3
L = lx
sc_v = Q/(l*w)
# number of iterations between output files
outf = 200
ii = 0
t = 0
def outerBndV(x) :
v = np.ndarray((len(x[:,0])))
v[:] = 0
v[np.abs(x[:,0]-L/2)<l/2] = sc_v
return v
def inBndP(x) :
v = np.ones_like(x[:,0])*(injector.P-1e5)
return v
fluid.set_wall_boundary("Lateral",pressure=None,velocity=[0,0])
fluid.set_wall_boundary("Bottom",pressure=None,velocity=[0,0])
fluid.set_open_boundary("Top",pressure=0,porosity=1,concentration=0)
#fluid.set_wall_boundary("Injection",pressure=None,velocity=[0,0])
fluid.set_open_boundary("Injection",pressure=inBndP,porosity=1,concentration=1)
fluid.set_strong_boundary("Bottom",1,0)
fluid.set_strong_boundary("Lateral",0,0)
#fluid.set_strong_boundary("Bottom",0,0)
#fluid.set_strong_boundary("Lateral",1,0)
if not no_particle :
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(), p.contact_forces(),init=True)
#set initial_condition
#c = np.zeros_like(fluid.coordinates()[:,1])
x = fluid.coordinates()
d = np.hypot(x[:,0]-0.034,x[:,1]-0.034)
#c[d<0.01] = 1
#fluid.set_concentration_cg(c)
c = np.zeros_like(fluid.coordinates()[:,1])
fluid.set_concentration_cg(c)
fluid.solution()[:,2] = (H-fluid.coordinates()[:,1])*(-g)*rho1
ii = 0
tic = time.time()
fluid.export_vtk(outputdir,0,0)
ioutput = 0
#for i in range(100) :
# t += dt
# if i%10 == 0 :
# ioutput += 1
# fluid.export_vtk(outputdir, t, ioutput)
# fluid.implicit_euler(dt*10, newton_max_it=20,skip_concentration=True)
alpha = 1./2.
while t < tEnd :
#Adaptation of the mesh. Args are minimal mesh radius, maximal mesh radius and number of elements
#Fluid solver
sf0 = np.copy(fluid.solution())
vp0 = np.copy(p.velocity())
pp0 = np.copy(p.position())
af0 = np.copy(fluid.concentration_dg())
cp0 = np.copy(p.contact_forces())
fluid.implicit_euler(dt, newton_max_it=20)
sf1 = np.copy(fluid.solution())
if not no_particle :
#if (ii%50==0):
# fluid.adapt_mesh(1e-2,1e-2/5,5000)
f0 = np.copy(fluid.compute_node_force(dt))
#Computation of the new velocities
vn = p.velocity() + f0 * dt / p.mass()
vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
#number of sub time step
nsub = max(5, int(np.ceil((vmax * dt * 8)/min(p.r()))))
print("NSUB", nsub,"VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(p.r()))
#Contact solver
for i in range(nsub) :
#tol = 5e-9
p.iterate(dt/nsub, f0)
p.position()[:,:] = pp0
p.contact_forces()[:,:] = cp0
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.contact_forces())
fluid.solution()[:,:] = sf0
#fluid.concentration_dg()[:] = af0
fluid.implicit_euler(dt, newton_max_it=20)
fluid.solution()[:,:] = (alpha*fluid.solution()+(1-alpha)*sf1)
fluid.advance_concentration(dt)
if not no_particle :
#if (ii%50==0):
# fluid.adapt_mesh(1e-2,1e-2/5,5000)
f1 = np.copy(fluid.compute_node_force(dt))
p.velocity()[:,:] = vp0
#Computation of the new velocities
vn = vp0 + (alpha*f0+(1-alpha)*f1)*dt/p.mass()
vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
#number of sub time step
nsub = max(5, int(np.ceil((vmax * dt * 8)/min(p.r()))))
print("NSUB", nsub,"VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(p.r()))
#Contact solver
for i in range(nsub) :
#tol = 5e-9
p.iterate(dt/nsub, (alpha*f0+(1-alpha)*f1))
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.contact_forces())
t += dt
flux = fluid.volume_flux("Injection")
injector.iterate((Q+flux*w)*rho0*dt)
print("fluxes : Q %g Domain : %g\n" % (Q, flux*w))
allp.append(injector.P)
#Output files writting
if ii %outf == 0 :
ioutput += 1
if not no_particle :
p.write_vtk(outputdir, ioutput, t)
fluid.export_vtk(outputdir, t, ioutput)
np.savetxt("pressurev5e-5sigma07", allp)
ii += 1
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))
L = .02;
H = .05;
l = 2e-3;
h = 4e-3;
lc = 0.02;
f = 8;
Point(1) = {0 , H, 0, lc/f};
Point(2) = {0 , 0, 0, lc/f};
Point(3) = {L , 0, 0, lc/f};
Point(4) = {L , H, 0, lc/f};
Point(5) = {L/2-l, 0, 0, lc/f};
Point(6) = {L/2+l, 0, 0, lc/f};
/*Point(7) = {L/2-l, h/4, 0, lc/f};
Point(8) = {L/2+l, h/4, 0, lc/f};
Point(9) = {L/2-l/2, h/4, 0, lc/f};
Point(10) = {L/2+l/2, h/4, 0, lc/f};
Point(11) = {L/2-l/2, h, 0, lc/f};
Point(12) = {L/2+l/2, h, 0, lc/f};*/
Line(1) = {1, 2};
Line(2) = {2, 5};
Line(3) = {5, 6};
Line(4) = {6, 3};
Line(5) = {3, 4};
Line(6) = {4, 1};
/*Line(7) = {12, 10};
Line(8) = {10, 8};
Line(9) = {8, 6};
Line(10) = {6, 3};
Line(11) = {3, 4};
Line(12) = {4, 1};*/
Line Loop(1) = {1:6};
Plane Surface(1) = {1};
Physical Line("Bottom") = {2,4};
Physical Line("Lateral") = {1,5};
Physical Line("Injection") = {3};
Physical Line("Top") = {6};
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