depot-fluid.py 4.7 KB
Newer Older
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
1
# Marblesbag - Copyright (C) <2010-2018>
2
3
# <Universite catholique de Louvain (UCL), Belgium
#  Universite de Montpellier, France>
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
4
# 	
5
# List of the contributors to the development of Marblesbag: see AUTHORS file.
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
6
7
8
# Description and complete License: see LICENSE file.
# 	
# This program (Marblesbag) is free software: 
9
# you can redistribute it and/or modify it under the terms of the GNU Lesser General 
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
10
11
12
13
14
15
# 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
16
# GNU Lesser General Public License for more details.
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
17
# 
18
19
# 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, 
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
20
21
# see <http://www.gnu.org/licenses/>.

dubois's avatar
maj  
dubois committed
22
#!/usr/bin/env python
23
24
from marblesbag import fluid as fluid
from marblesbag import scontact2
Matthieu Constant's avatar
Matthieu Constant committed
25
from marblesbag import lmgc2Interface
26

dubois's avatar
maj  
dubois committed
27
28
29
import numpy as np
import os
import time
30
import shutil
dubois's avatar
maj  
dubois committed
31
import random
Matthieu Constant's avatar
Matthieu Constant committed
32
#grains creation and placement + physical boundaries definition
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
33
def genInitialPosition(filename, r, ox, oy, lx, ly, rhop) :
dubois's avatar
maj  
dubois committed
34
    p = scontact2.ParticleProblem()
Matthieu Constant's avatar
Matthieu Constant committed
35
    #Loading of the mesh.msh file specifying physical boundaries name
36
    p.load_msh_boundaries("mesh.msh", ["Top", "Box"])
Matthieu Constant's avatar
Matthieu Constant committed
37
    #Definition of the points where the grains are located
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
38
39
    x = np.arange(ox+r, ox+lx, 2*r)
    y = np.arange(oy+r, oy+ly, 2*r)
dubois's avatar
maj  
dubois committed
40
41
    for i in range(x.shape[0]) :
       for j in range(y.shape[0]) :
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
42
          rr=r*(0.95+random.random()*0.05)
Matthieu Constant's avatar
Matthieu Constant committed
43
          #Addition of an particle object at each point
44
          p.add_particle((x[i], y[j]), rr, rr**2 * np.pi * rhop);
Matthieu Constant's avatar
Matthieu Constant committed
45
    p.write_vtk(filename,0,0)
dubois's avatar
maj  
dubois committed
46
47


48
outputdir = "output"
dubois's avatar
maj  
dubois committed
49
50
51
52
53
54
55
if not os.path.isdir(outputdir) :
    os.makedirs(outputdir)

t = 0
ii = 0

#physical parameters
Matthieu Constant's avatar
Matthieu Constant committed
56
57
58
59
60
61
62
63
64
g =  -9.81                                      # gravity
rho = 1200                                      # fluid density
rhop = 1800                                     # grains density
nu = 1e-6                                       # kinematic viscosity
tEnd = 7                                        # final time
r=1e-4                                          # grains radius
y0part = 40e-3                                  # vertical center of the grains
lypart = 20e-3                                  # semi-height of the grains
lxpart = 0.015                                  # semi-width of the grains
dubois's avatar
maj  
dubois committed
65
66

#numerical parameters
Matthieu Constant's avatar
Matthieu Constant committed
67
68
dt = 1e-3                                       # time step
epsilon = 1e-4                                  # stabilization parameter
dubois's avatar
maj  
dubois committed
69

Matthieu Constant's avatar
Matthieu Constant committed
70
print("eps = %g" % ( epsilon))
Matthieu Constant's avatar
Matthieu Constant committed
71
#Enable the use of lmgc90
Matthieu Constant's avatar
Matthieu Constant committed
72
use_lmgc = True
Matthieu Constant's avatar
Matthieu Constant committed
73
74
75


#Object particles creation
Matthieu Constant's avatar
Matthieu Constant committed
76
genInitialPosition(outputdir, r, 0., y0part, lxpart, lypart, rhop)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
77
if use_lmgc :
Matthieu Constant's avatar
Matthieu Constant committed
78
79
    friction=0.1                                # friction coefficient
    lmgc2Interface.scontactToLmgc90(outputdir, 0, friction)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
80
81
    p = lmgc2Interface.LmgcInterface()
else :
Matthieu Constant's avatar
Matthieu Constant committed
82
83
    p = scontact2.ParticleProblem()
    p.read_vtk(outputdir,0)
dubois's avatar
maj  
dubois committed
84
85

print("r = %g, m = %g\n" %  (p.r()[0], p.mass()[0]))
Matthieu Constant's avatar
Matthieu Constant committed
86
outf = 20                                       # number of iterations between output files
dubois's avatar
maj  
dubois committed
87
88


Matthieu Constant's avatar
Matthieu Constant committed
89
90
#Object fluid creation + Boundary condition of the fluid (field 0 is horizontal velocity; field 1 is vertical velocity; field 2 is pressure)
#Format: strong_boundaries = [(Boundary tag, Fluid field, Value)
91
strong_boundaries = [("Box",0,0),("Box",1,0),("Top",0,0.),("Top",1,0.),("Top",2,0.)]
Matthieu Constant's avatar
Matthieu Constant committed
92
fluid = fluid.fluid_problem("mesh.msh",g,nu*rho,rho,epsilon,strong_boundaries)
dubois's avatar
maj  
dubois committed
93
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
94
fluid.export(outputdir,0,0)
dubois's avatar
maj  
dubois committed
95
96

tic = time.clock()
Matthieu Constant's avatar
Matthieu Constant committed
97
#Computation loop
Matthieu Constant's avatar
Matthieu Constant committed
98
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
Matthieu Constant's avatar
Matthieu Constant committed
99
100
while t < tEnd :
    #Fluid solver
Matthieu Constant's avatar
Matthieu Constant committed
101
    fluid.implicit_euler(dt)
102
    forces = fluid.compute_node_force(dt)
Matthieu Constant's avatar
Matthieu Constant committed
103
    #Computation of the new velocities
104
105
    vn = p.velocity() + forces * dt / p.mass()
    vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
Matthieu Constant's avatar
Matthieu Constant committed
106
    #number of sub time step
Matthieu Constant's avatar
Matthieu Constant committed
107
    nsub = max(1, int(np.ceil((vmax * dt * 4)/min(p.r()))))
108
    print("NSUB", nsub,"VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(p.r()))
Matthieu Constant's avatar
Matthieu Constant committed
109
    #Contact solver
110
    for i in range(nsub) :
Matthieu Constant's avatar
Matthieu Constant committed
111
112
113
        tol = 1e-6
        p.iterate(dt/nsub, forces, tol)
    #Localisation of the grains in the fluid
dubois's avatar
maj  
dubois committed
114
115
    fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
    t += dt
Matthieu Constant's avatar
Matthieu Constant committed
116
    #Output files writting
dubois's avatar
maj  
dubois committed
117
118
    if ii %outf == 0 :
        ioutput = int(ii/outf) + 1
Matthieu Constant's avatar
Matthieu Constant committed
119
120
        p.write_vtk(outputdir, ioutput, t)
        fluid.export_vtk(outputdir, t, ioutput)
dubois's avatar
maj  
dubois committed
121
    ii += 1
122
    print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.clock() - tic))
dubois's avatar
maj  
dubois committed
123