avalanch1fluidnofriction.py 4.8 KB
Newer Older
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
1
# MigFlow - Copyright (C) <2010-2020>
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
# <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
23
24
25

# TESTCASE DESCRIPTION:
# Collapsing of a grain column in fluid. 
Matthieu Constant's avatar
Matthieu Constant committed
26
27
# This test case can be runned on its own, but it is possible to use the output files of the depot test case (../depot/depot.py) to create a compact column as initial situation
# Otherwise an hexagonal packing is created
28
from migflow import fluid
Matthieu Constant's avatar
Matthieu Constant committed
29
from migflow import scontact
30
from migflow import time_integration
31
32
33
34
35
36
37
38
39
40
41

import numpy as np
import os
import time
import shutil
import random

outputdir = "output"
if not os.path.isdir(outputdir) :
    os.makedirs(outputdir)

Matthieu Constant's avatar
Matthieu Constant committed
42

43
# Physical parameters
Matthieu Constant's avatar
Matthieu Constant committed
44
g = np.array([0,-9.81])                   # gravity
45
46
47
rhop = 1500                               # grains density
rhof = 1000                               # fluid density
nuf = 1e-6                                # kinematic fluid viscosity
48
H = 0.2                                   # domain height
49

50
51
# Numerical parameters
outf = 10                                 # number of iterations between output files
Matthieu Constant's avatar
Matthieu Constant committed
52
dt = 1e-3                                 # time step
53
tEnd = 10                                 # final time
Nathan Coppin's avatar
Nathan Coppin committed
54
use_pre_deposit = True                   # set if the deposit in ../depot/depot.py has been achieved. If false an hexagonal packing is used for the inital column
55

Matthieu Constant's avatar
Matthieu Constant committed
56
57
58
# 
# PARTICLES PROBLEM
#
Matthieu Constant's avatar
Matthieu Constant committed
59
60

# Particles container for the current problem
Matthieu Constant's avatar
Matthieu Constant committed
61
p = scontact.ParticleProblem(2)
62
p.load_msh_boundaries("mesh.msh", ["Top", "Bottom", "Left", "Right"])
Matthieu Constant's avatar
Matthieu Constant committed
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
def hexagonal_packing(x0,y0,lx,ly,rmin,rmax) :
    """Set all the particles centre positions using an hexagonal arrangement
    
    Keyword arguments:
    x0 -- X-axis bottom left corner
    y0 -- Y-axis bottom left corner 
    ly -- domain height
    lx -- domain width 
    rmin -- minimum radius
    rmax -- maximum radius
    """
    for x in np.arange(x0+rmax,x0+lx-rmax,2*rmax) :
        for y in np.arange(y0+rmax,y0+ly-rmax,2*(3**0.5)*rmax):
            z = np.random.random((1))
            r = rmin+(rmax-rmin)*z
            p.add_particle((x, y), r, r**2*np.pi*rhop)
    for x in np.arange(2*rmax+x0,x0+lx-rmax,2*rmax) :
        for y in np.arange(y0+3**0.5*rmax+rmax,y0+ly-rmax,2*(3**0.5)*rmax):
            z = np.random.random((1))
            r = rmin+(rmax-rmin)*z
            p.add_particle((x, y), r, r**2*np.pi*rhop)
if use_pre_deposit:
    # Particles container structure builder to load file created with depot.py.
    p1 = scontact.ParticleProblem(2)
    # Use particles deposit computed by the depot.py file of "depot" folder.
    p1.read_vtk("../depot/output",25)
Nathan Coppin's avatar
Nathan Coppin committed
89
    p.add_particles(p1.position(),p1.r(),p1.mass(),v=p1.velocity())
Matthieu Constant's avatar
Matthieu Constant committed
90
91
92
93
else:
    r = 5e-4                              # particles radius
    hexagonal_packing(0,0,0.08,0.165,0.9*r,1.1*r)
    
94
p.write_vtk(outputdir, 0, 0)
Matthieu Constant's avatar
Matthieu Constant committed
95
96

# Initial time and iteration.
97
98
t = 0
ii = 0
99

Matthieu Constant's avatar
Matthieu Constant committed
100
101
102
103
#
# FLUID PROBLEM
#
# Fluid structure builder.
Matthieu Constant's avatar
Matthieu Constant committed
104
fluid = fluid.FluidProblem(2,g,[nuf*rhof],[rhof])
Matthieu Constant's avatar
Matthieu Constant committed
105
# Set the mesh geometry for the fluid computation.
106
fluid.load_msh("mesh.msh")
Matthieu Constant's avatar
Matthieu Constant committed
107
108
109
110
fluid.set_wall_boundary("Bottom")
fluid.set_wall_boundary("Left")
fluid.set_wall_boundary("Top",pressure=0)
fluid.set_wall_boundary("Right")
Matthieu Constant's avatar
Matthieu Constant committed
111
112
113
114
fluid.set_strong_boundary("Bottom",1,0)
fluid.set_strong_boundary("Top",1,0)
fluid.set_strong_boundary("Left",0,0)
fluid.set_strong_boundary("Right",0,0)
Matthieu Constant's avatar
Matthieu Constant committed
115
# Set locations of the grains in the mesh and compute the porosity in each computation cell.
116
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(), p.contact_forces())
Matthieu Constant's avatar
Matthieu Constant committed
117
# Write output files for post-visualization.
Nathan Coppin's avatar
Nathan Coppin committed
118
fluid.write_vtk(outputdir,0,0)
Matthieu Constant's avatar
Matthieu Constant committed
119

120
tic = time.time()
Matthieu Constant's avatar
Matthieu Constant committed
121

122
#
Matthieu Constant's avatar
Matthieu Constant committed
123
# COMPUTATION LOOP
124
125
#
while t < tEnd :
Nathan Coppin's avatar
Nathan Coppin committed
126
    time_integration.iterate(fluid, p, dt, min_nsub=10, external_particles_forces=g*p.mass())
127
    t += dt
128

Matthieu Constant's avatar
Matthieu Constant committed
129
    # Output files writing.
130
    if ii%outf == 0 :
131
132
        ioutput = int(ii/outf) + 1
        p.write_vtk(outputdir, ioutput, t)
Nathan Coppin's avatar
Nathan Coppin committed
133
        fluid.write_vtk(outputdir, ioutput, t)
134
    ii += 1
Matthieu Constant's avatar
Matthieu Constant committed
135
    print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))