avalanch2fluids.py 5.72 KB
Newer Older
1
# MigFlow - Copyright (C) <2010-2018>
Matthieu Constant's avatar
Matthieu Constant committed
2
3
4
# <Universite catholique de Louvain (UCL), Belgium
#  Universite de Montpellier, France>
# 	
5
# List of the contributors to the development of MigFlow: see AUTHORS file.
Matthieu Constant's avatar
Matthieu Constant committed
6
7
# Description and complete License: see LICENSE file.
# 	
8
# This program (MigFlow) is free software: 
Matthieu Constant's avatar
Matthieu Constant committed
9
10
11
12
13
14
15
16
17
18
19
20
21
22
# 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
Matthieu Constant's avatar
Matthieu Constant committed
23
24
25
26
27
28
29
30
31

# TESTCASE DESCRIPTION
# Collapsing of a grain column in fluid. 
# The grain column is initially mixed in sea water while the rest of the domain contain soft water. 
# particles are initially placed using a random loose packing method of the lmgc90 software.
# The lmgc90 software is also used to solve contacts in this test cases.

from migflow import fluid
from migflow import scontact
Frédéric Dubois's avatar
Frédéric Dubois committed
32
from migflow import lmgc90Interface
33
from migflow import time_integration
Matthieu Constant's avatar
Matthieu Constant committed
34
35
36
37
38
39
from pylmgc90 import pre
import numpy as np
import os
import time
import shutil
import random
40
use_lmgc = False
Matthieu Constant's avatar
Matthieu Constant committed
41

42
def genInitialPosition(filename, N, r, lx, ly, rhop,use_lmgc) :
Matthieu Constant's avatar
Matthieu Constant committed
43
44
45
46
47
48
49
50
51
52
53
    """Create a container to initialize particles and write particles in an initial output file.

    Keyword arguments:
    filename -- directory to write the output file
    N -- maximum number of grains
    r -- mean particle radius
    lx -- particle area width
    ly -- particle aera height
    rhop -- particle density
    """
    # Particles structure builder
54
55
56
57
    if use_lmgc :
      p = scontact.ParticleProblem(2)
    else:
      p = scontact.ParticleProblem(2,True,True)
Matthieu Constant's avatar
Matthieu Constant committed
58
    # Load mesh.msh file specifying physical boundaries names
Matthieu Constant's avatar
Matthieu Constant committed
59
    p.load_msh_boundaries("mesh.msh", ["Top", "Bottom", "Right", "Left"])
Matthieu Constant's avatar
Matthieu Constant committed
60
    # Radii defining the interval of the granulo
Matthieu Constant's avatar
Matthieu Constant committed
61
62
    Ra = 0.9*r
    Rb = 1.1*r
Matthieu Constant's avatar
Matthieu Constant committed
63
    # Random list of radii in the interval
Matthieu Constant's avatar
Matthieu Constant committed
64
    radii = pre.granulo_Random(N, Ra, Rb)
Matthieu Constant's avatar
Matthieu Constant committed
65
    # Nb of particles and positions in the geometry defined by arguments: [0,lx*4]*[0,ly/10]
Matthieu Constant's avatar
Matthieu Constant committed
66
    [nb_laid_particles, coors] = pre.depositInBox2D(radii,lx*4,ly/10)
Matthieu Constant's avatar
Matthieu Constant committed
67
    # Add a grain at each centre position
Matthieu Constant's avatar
Matthieu Constant committed
68
    for i in range(nb_laid_particles):
Matthieu Constant's avatar
Matthieu Constant committed
69
       p.add_particle(coors[2*i:2*i+2], radii[i], radii[i]**2 * np.pi * rhop)
Matthieu Constant's avatar
Matthieu Constant committed
70
    # Particles horizontal translation
Matthieu Constant's avatar
Matthieu Constant committed
71
    p.position()[:, 0] += 0.1
72

Matthieu Constant's avatar
Matthieu Constant committed
73
    # Second particles area
Matthieu Constant's avatar
Matthieu Constant committed
74
75
76
    radii = pre.granulo_Random(N, Ra, Rb)
    [nb_laid_particles, coors] = pre.depositInBox2D(radii,lx,ly)
    for i in range(nb_laid_particles):
Matthieu Constant's avatar
Matthieu Constant committed
77
       p.add_particle(coors[2*i:2*i+2], radii[i], radii[i]**2 * np.pi * rhop)
Matthieu Constant's avatar
Matthieu Constant committed
78
79
    p.write_vtk(filename,0,0)

Matthieu Constant's avatar
Matthieu Constant committed
80
# Define output directory
Matthieu Constant's avatar
Matthieu Constant committed
81
82
83
84
outputdir = "output2fluids"
if not os.path.isdir(outputdir) :
    os.makedirs(outputdir)

Matthieu Constant's avatar
Matthieu Constant committed
85
# Physical parameters
86
g = np.array([0,-9.81])                       # gravity
Matthieu Constant's avatar
Matthieu Constant committed
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
# Particles
r = 0.0005                                    # radius
lx = 0.1                                      # particles area width
ly = 0.18                                     # particles area height
N = 15000                                     # number of particles
rhop = 1100                                   # solid density
# Soft water
rhof = 1000                                   # soft water density
nuf = 1e-6                                    # soft water viscosity
# Mixture properties
rhom = 1050                                   # sea water density
num = 1e-6                                    # sea water kinematic viscosity

# Numerical parameters
outf = 10                                     #number of iterations between output files
dt = 1e-4                                     #time step
tEnd = 10                                     #final time

# 
# PARTICLE PROBLEM
# 
# Initialise particles
109
110
111
112
113
114
115
116
genInitialPosition(outputdir, N, r, lx, ly, rhop, use_lmgc)
if use_lmgc :
  friction = 0.3                                #friction coefficient
  lmgc90Interface.scontactTolmgc90(outputdir, 2, 0, friction)
  p = lmgc90Interface.ParticleProblem(2)
else:
  p  = scontact.ParticleProblem(2,True,True)
  p.read_vtk("output2fluids",0)
Matthieu Constant's avatar
Matthieu Constant committed
117

Matthieu Constant's avatar
Matthieu Constant committed
118
119
120
121
122
123
124
125
126
127
# Initial time and iteration
t = 0
ii = 0

# 
# FLUID PROBLEM
# 
fluid = fluid.FluidProblem(2,g,[num*rhom,nuf*rhof],[rhom,rhof])
# Set the mesh geometry for the fluid computation
fluid.load_msh("mesh.msh")
Matthieu Constant's avatar
Matthieu Constant committed
128
129
130
131
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
132
# Set location of the particles in the mesh and compute the porosity in each computation cell
133
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(), p.contact_forces(), init=True)
134
135
136
137
#Access to fluid fields and node coordiantes
s         =  fluid.solution()
x         =  fluid.coordinates()
#Set initial condition for the concentration
Matthieu Constant's avatar
Matthieu Constant committed
138
139
140
141
c = np.ndarray((fluid.n_nodes()))
c[:] = 0
c[x[:,0]<lx] = 1
fluid.set_concentration_cg(c)
Matthieu Constant's avatar
Matthieu Constant committed
142
143
fluid.export_vtk(outputdir,0,0)

Matthieu Constant's avatar
Matthieu Constant committed
144
tic       =  time.time()
Jonathan Lambrechts's avatar
fix    
Jonathan Lambrechts committed
145
G = np.zeros((p.n_particles(),p.dim())
146
G[:,1] = g[1]*p.mass()[:,0]
Matthieu Constant's avatar
Matthieu Constant committed
147
148
149
# 
# COMPUTATION LOOP
# 
Matthieu Constant's avatar
Matthieu Constant committed
150
while t < tEnd : 
151
152
    time_integration.iterate(fluid, p, dt, min_nsub=10, external_particles_forces=G)
    t += dt
153

Matthieu Constant's avatar
Matthieu Constant committed
154
    # Output files writing
Matthieu Constant's avatar
Matthieu Constant committed
155
156
157
158
    if ii %outf == 0 :
        ioutput = int(ii/outf) + 1
        p.write_vtk(outputdir, ioutput, t)
        fluid.export_vtk(outputdir, t, ioutput)
159
    ii    += 1
Matthieu Constant's avatar
Matthieu Constant committed
160
    print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))