avalanch2fluids.py 7.25 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
23
24
25
26
27
28
29
30
31
# 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 as fluid
from migflow import scontact2
from migflow import lmgc2Interface
from pylmgc90 import pre
import numpy as np
import os
import time
import shutil
import random
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
'''
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. 
Grains 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.
'''

'''
genInitialPosition is a function that sets all the grains centre positions and create the grains objects to add in the computing structure
Args:    -filename    : name of the output file
         -N           : max number of grains
         -r           : mean radius of the grains
         -lx          : domain width 
         -ly          : domain height
         -rhop        : grains density        
'''
Matthieu Constant's avatar
Matthieu Constant committed
47
def genInitialPosition(filename, N, r, lx, ly, rhop) :
48
    #Grains structure builder
Matthieu Constant's avatar
Matthieu Constant committed
49
    p = scontact2.ParticleProblem()
50
    #Load mesh.msh file specifying physical boundaries names
Matthieu Constant's avatar
Matthieu Constant committed
51
    p.load_msh_boundaries("mesh.msh", ["Top", "Bottom", "Right", "Left"])
52
    #Radii defining the interval of the granulo
Matthieu Constant's avatar
Matthieu Constant committed
53
54
    Ra = 0.9*r
    Rb = 1.1*r
55
    #Random list of radii in the interval
Matthieu Constant's avatar
Matthieu Constant committed
56
    radii = pre.granulo_Random(N, Ra, Rb)
57
    #Nb of grains and positions in the geometry defined by arguments: [0,lx*4]*[0,ly/10]
Matthieu Constant's avatar
Matthieu Constant committed
58
    [nb_laid_particles, coors] = pre.depositInBox2D(radii,lx*4,ly/10)
59
    #Add a grain at each centre position
Matthieu Constant's avatar
Matthieu Constant committed
60
    for i in range(nb_laid_particles):
Matthieu Constant's avatar
Matthieu Constant committed
61
       p.add_particle(coors[2*i:2*i+2], radii[i], radii[i]**2 * np.pi * rhop)
62
    #Grains horizontal translation
Matthieu Constant's avatar
Matthieu Constant committed
63
    p.position()[:, 0] += 0.1
64
65

    #Second grains area
Matthieu Constant's avatar
Matthieu Constant committed
66
67
68
    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
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
71
    p.write_vtk(filename,0,0)

72
#Define output directory
Matthieu Constant's avatar
Matthieu Constant committed
73
74
75
76
77
outputdir = "output2fluids"
if not os.path.isdir(outputdir) :
    os.makedirs(outputdir)

#physical parameters
78
79
80
81
82
83
84
85
86
87
g         = -9.81                                     #gravity
##grains
r         =  0.0005                                   #radius
lx        =  0.1                                      #grains area width
ly        =  0.18                                     #grains area height
N         =  15000                                    #number of grains
rhop      =  1100                                     #solid density
##soft water
rhof      =  1000                                     #soft water density
nuf       =  1e-6                                     #soft water viscosity
Matthieu Constant's avatar
Matthieu Constant committed
88
##mixture properties
89
90
rhom      =  1050                                     #sea water density
num       =  1e-6                                     #sea water kinematic viscosity
Matthieu Constant's avatar
Matthieu Constant committed
91
92

#numerical parameters
93
94
95
outf      =  10                                       #number of iterations between output files
dt        =  1e-3                                     #time step
tEnd      =  10                                       #final time
Matthieu Constant's avatar
Matthieu Constant committed
96

97
#Initialise grains
Matthieu Constant's avatar
Matthieu Constant committed
98
genInitialPosition(outputdir, N, r, lx, ly, rhop)
99
friction  =  0.3                                      #friction coefficient
Matthieu Constant's avatar
Matthieu Constant committed
100
lmgc2Interface.scontactToLmgc90(outputdir, 0, friction)
101
102
103
104
105
106
107
108
109
110
111
p         =  lmgc2Interface.LmgcInterface()
#Initial time and iteration
t         =  0
ii        =  0

'''
fluid_problem is the function builder of the fluid class
Args:     -g             : gravity
          -viscosity     : list of fluid phases viscosities
          -density       : list of fluid phases densities
'''
Matthieu Constant's avatar
Matthieu Constant committed
112
fluid = fluid.fluid_problem(g,[num*rhom,nuf*rhof],[rhom,rhof])
113
#Set the mesh geometry for the fluid computation
Matthieu Constant's avatar
Matthieu Constant committed
114
115
fluid.load_msh("mesh.msh")

116
117
118
119
120
121
'''
set_strong_boundary is the function setting the strong boundaries (=constraints) for the fluid problem
Args:     -Tag           : physical tag (set in the mesh.geo file) of the boundary on which you want to add a constraint
          -Field         : O: x-velocity; 1: y-velocity; 2: pressure
          -Value         : value assigned to the specified field on the specified boundary
'''
Matthieu Constant's avatar
Matthieu Constant committed
122
123
124
125
126
fluid.set_strong_boundary("Top",2,0)
fluid.set_strong_boundary("Top",1,0)
fluid.set_strong_boundary("Bottom",1,0)
fluid.set_strong_boundary("Right",0,0)
fluid.set_strong_boundary("Left",0,0)
127
128
129
130
131
'''
set_weak_boundary is the function setting the weak boundaries (=normal fluxes) for the fluid problem
Args:     -Tag           : physical tag (set in the mesh.geo file) of the boundary on which you want to add a constraint
          -Bnd Type      : boundaries can be of type: wall, inflow, outflow
'''
Matthieu Constant's avatar
Matthieu Constant committed
132
133
134
135
136
fluid.set_weak_boundary("Bottom","Wall")
fluid.set_weak_boundary("Left","Wall")
fluid.set_weak_boundary("Top","Wall")
fluid.set_weak_boundary("Right","Wall")

137
#Set location of the grains in the mesh and compute the porosity in each computation cell
Matthieu Constant's avatar
Matthieu Constant committed
138
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
139
140
141
142
#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
143
144
145
146
147
for i in range(len(x[:,0])):
    if (x[i,0])<lx:
        s[i,3] = 1
fluid.export_vtk(outputdir,0,0)

Matthieu Constant's avatar
Matthieu Constant committed
148
tic       =  time.time()
Matthieu Constant's avatar
Matthieu Constant committed
149
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
150
#Computation loop
Matthieu Constant's avatar
Matthieu Constant committed
151
152
153
while t < tEnd : 
    #Fluid solver
    fluid.implicit_euler(dt)
154
155

    #Adaptation of the mesh. Args are minimal mesh radius, maximal mesh radius and number of elements
Matthieu Constant's avatar
Matthieu Constant committed
156
157
    if (ii%15==0 and ii != 0):
       fluid.adapt_mesh(1e-2,3e-3,10000)
158

Matthieu Constant's avatar
Matthieu Constant committed
159
    #Computation of the new velocities
160
161
162
163
164
    forces = fluid.compute_node_force(dt)
    vn     = p.velocity() + forces * dt / p.mass()
    vmax   = np.max(np.hypot(vn[:, 0], vn[:, 1]))
    #number of contact sub time step
    nsub   = max(1, int(np.ceil((vmax * dt * 4)/min(p.r()))))
Matthieu Constant's avatar
Matthieu Constant committed
165
166
167
168
    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)  
169
170

    t     += dt
Matthieu Constant's avatar
Matthieu Constant committed
171
    fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
172
173

    #Output files writing
Matthieu Constant's avatar
Matthieu Constant committed
174
175
176
177
    if ii %outf == 0 :
        ioutput = int(ii/outf) + 1
        p.write_vtk(outputdir, ioutput, t)
        fluid.export_vtk(outputdir, t, ioutput)
178
    ii    += 1
Matthieu Constant's avatar
Matthieu Constant committed
179
    print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))