inject friction

testcases/injectionWGrains/Q0.05Lyon/adapt/mesh.geo

+L = .136;
+H = .25;
+l = 2e-3;
+lc = 0.01;
+f = 10;
+
+Point(1) = {0, H, 0};
+Point(2) = {0, 0, 0};
+Point(3) = {L/2-l/2,0,0};
+Point(4) = {L/2+l/2,0,0};
+Point(5) = {L,0,0};
+Point(6) = {L,H,0};
+
+Point(9) = {L/2,0,0};
+Point(10) = {L/2,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 Loop(1) = {1:6};
+
+Plane Surface(1) = {1};
+Physical Line("Bottom") = {2,4};
+Physical Line("Lateral") = {1,5};
+Physical Line("Top") = {6};
+Physical Line("Injection") = {3};
+Physical Surface("Domain") = {1};
+Physical Point("PtFix") = {1};
+
+Mesh.Algorithm = 5;
+
+Merge "lc.pos";
+Field[1] = PostView;
+Field[1].IView = 0;
+Background Field = 1;
+
+Mesh.CharacteristicLengthExtendFromBoundary = 0;
+Mesh.CharacteristicLengthFromPoints = 0;
+Mesh.MshFileVersion = 2;

testcases/injectionWGrains/Q0.05Lyon/inject2fGrains.py

+# 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
+
+
+
+
+outputdir = "outputVid1"
+outputdir1 = "outputVid"
+if not os.path.isdir(outputdir) :
+    os.makedirs(outputdir)
+
+t = 0
+ii = 0
+
+def genInitialPosition(filename, N, r, lx, ly, rhop) :
+
+    p = scontact.ParticleProblem(2)
+    #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], 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
+print(g)
+rho0 = 1.117                                      # fluid density
+rho1 = 785.92
+nu0 = 1.57e-5                                   # kinematic viscosity
+nu1 = 1.2e-3/rho1
+tEnd = 50                                     # final time
+r = 5e-4/2
+N = 100000
+lx = .136
+ly = .09
+rhop = 1059
+#numerical parameters
+dt = .001                                  # time step
+
+shutil.copy("mesh.msh", outputdir +"/mesh.msh")
+genInitialPosition(outputdir, N, r, lx, ly, rhop)
+friction=0.3
+lmgc90Interface.scontactTolmgc90(outputdir, 2, 0, friction)
+p = lmgc90Interface.ParticleProblem(2)
+p.write_vtk(outputdir,0,0)
+#p = scontact.ParticleProblem(2)
+#p.read_vtk(outputdir,0)
+# number of iterations between output files
+outf = 25
+ii = 0
+t = 0
+def outerBndV(x) :
+    #print(0.265258*min((6*t**5-15*t**4+10*t**3),1))
+    return 0.265258*min((6*t**5-15*t**4+10*t**3),1)
+
+fluid = fluid.FluidProblem(2,g,[nu0*rho0,nu1*rho1],[rho0,rho1],coeff_stab=0.001)
+fluid.load_msh("mesh.msh")
+fluid.set_weak_boundary("Bottom","Wall")
+fluid.set_weak_boundary("Lateral","Wall")
+fluid.set_weak_boundary("Top","pressure",[0,1])
+fluid.set_weak_boundary("Injection","velocity", [0,outerBndV,1])
+
+fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
+
+#set initial_condition
+fluid.solution()[fluid.coordinates()[:,1]>0.18,3] = .5
+fluid.solution()[fluid.coordinates()[:,1]>0.19,3] = 1
+
+fluid.export_vtk(outputdir,0,0)
+
+ii = 0
+tic = time.time() 
+
+fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
+while t < tEnd : 
+#Adaptation of the mesh. Args are minimal mesh radius, maximal mesh radius and number of elements
+    #Fluid solver
+    fluid.implicit_euler(dt, newton_max_it=20)
+    if (ii%50==0):
+       fluid.adapt_mesh(1e-2,1e-2/5,5000)
+    forces = fluid.compute_node_force(dt)
+    #Computation of the new velocities
+    vn = p.velocity() + forces * dt / p.mass()
+    vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
+    #number of sub time step
+    nsub = max(1, int(np.ceil((vmax * dt * 4)/min(p.r()))))
+    print("NSUB", nsub,"VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(p.r()))
+    #Contact solver
+    for i in range(nsub) :
+        p.iterate(dt/nsub, forces)  
+    t += dt
+    fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
+    #Output files writting
+    if ii %outf == 0 :
+        ioutput = int(ii/outf) + 1
+        p.write_vtk(outputdir, ioutput, t)
+        fluid.export_vtk(outputdir, t, ioutput)
+    ii += 1
+    print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))

testcases/injectionWGrains/Q0.05Lyon/mesh.geo

+L = .136;
+H = .25;
+l = 2e-3;
+lc = 0.01;
+f = 4;
+
+Point(1) = {0, H, 0};
+Point(2) = {0, 0, 0};
+Point(3) = {L/2-l/2,0,0};
+Point(4) = {L/2+l/2,0,0};
+Point(5) = {L,0,0};
+Point(6) = {L,H,0};
+
+Point(9) = {L/2,0,0};
+Point(10) = {L/2,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 Loop(1) = {1:6};
+
+Plane Surface(1) = {1};
+Physical Line("Bottom") = {2,4};
+Physical Line("Lateral") = {1,5};
+Physical Line("Top") = {6};
+Physical Line("Injection") = {3};
+Physical Surface("Domain") = {1};
+Physical Point("PtFix") = {1};
+
+Field[1] = Attractor;
+Field[1].EdgesList = {9};
+Field[1].NNodesByEdge = 200;
+
+Field[2] = Threshold;
+Field[2].DistMax = 10*l;
+Field[2].DistMin = l;
+Field[2].LcMax = lc/f;
+Field[2].LcMin = lc/f;
+Field[2].IField = 1;
+
+Background Field = 2;
+Mesh.MshFileVersion = 2;