avalanch2fluids.py 6.08 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
Matthieu Constant's avatar
Matthieu Constant committed
33
34
35
36
37
38
from pylmgc90 import pre
import numpy as np
import os
import time
import shutil
import random
Matthieu Constant's avatar
Matthieu Constant committed
39

Matthieu Constant's avatar
Matthieu Constant committed
40
def genInitialPosition(filename, N, r, lx, ly, rhop) :
Matthieu Constant's avatar
Matthieu Constant committed
41
42
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
    p = scontact.ParticleProblem(2)
    # Load mesh.msh file specifying physical boundaries names
Matthieu Constant's avatar
Matthieu Constant committed
54
    p.load_msh_boundaries("mesh.msh", ["Top", "Bottom", "Right", "Left"])
Matthieu Constant's avatar
Matthieu Constant committed
55
    # Radii defining the interval of the granulo
Matthieu Constant's avatar
Matthieu Constant committed
56
57
    Ra = 0.9*r
    Rb = 1.1*r
Matthieu Constant's avatar
Matthieu Constant committed
58
    # Random list of radii in the interval
Matthieu Constant's avatar
Matthieu Constant committed
59
    radii = pre.granulo_Random(N, Ra, Rb)
Matthieu Constant's avatar
Matthieu Constant committed
60
    # Nb of particles and positions in the geometry defined by arguments: [0,lx*4]*[0,ly/10]
Matthieu Constant's avatar
Matthieu Constant committed
61
    [nb_laid_particles, coors] = pre.depositInBox2D(radii,lx*4,ly/10)
Matthieu Constant's avatar
Matthieu Constant committed
62
    # Add a grain at each centre position
Matthieu Constant's avatar
Matthieu Constant committed
63
    for i in range(nb_laid_particles):
Matthieu Constant's avatar
Matthieu Constant committed
64
       p.add_particle(coors[2*i:2*i+2], radii[i], radii[i]**2 * np.pi * rhop)
Matthieu Constant's avatar
Matthieu Constant committed
65
    # Particles horizontal translation
Matthieu Constant's avatar
Matthieu Constant committed
66
    p.position()[:, 0] += 0.1
67

Matthieu Constant's avatar
Matthieu Constant committed
68
    # Second particles area
Matthieu Constant's avatar
Matthieu Constant committed
69
70
71
    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
72
       p.add_particle(coors[2*i:2*i+2], radii[i], radii[i]**2 * np.pi * rhop)
Matthieu Constant's avatar
Matthieu Constant committed
73
74
    p.write_vtk(filename,0,0)

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

Matthieu Constant's avatar
Matthieu Constant committed
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
# Physical parameters
g = -9.81                                     # gravity
# 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
Matthieu Constant's avatar
Matthieu Constant committed
104
genInitialPosition(outputdir, N, r, lx, ly, rhop)
Matthieu Constant's avatar
Matthieu Constant committed
105
friction = 0.3                                #friction coefficient
Frédéric Dubois's avatar
Frédéric Dubois committed
106
107
lmgc90Interface.scontactTolmgc90(outputdir, 2, 0, friction)
p = lmgc90Interface.ParticleProblem(2)
Matthieu Constant's avatar
Matthieu Constant committed
108

Matthieu Constant's avatar
Matthieu Constant committed
109
110
111
112
113
114
115
116
117
118
# 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")
Nathan Coppin's avatar
Nathan Coppin committed
119
120
121
122
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
123
# Set location of the particles in the mesh and compute the porosity in each computation cell
Matthieu Constant's avatar
Matthieu Constant committed
124
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
125
126
127
128
#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
129
130
131
132
133
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
134
tic       =  time.time()
Matthieu Constant's avatar
Matthieu Constant committed
135
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
Matthieu Constant's avatar
Matthieu Constant committed
136
137
138
139

# 
# COMPUTATION LOOP
# 
Matthieu Constant's avatar
Matthieu Constant committed
140
while t < tEnd : 
Matthieu Constant's avatar
Matthieu Constant committed
141
    # Fluid solver
Matthieu Constant's avatar
Matthieu Constant committed
142
    fluid.implicit_euler(dt)
143

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

Matthieu Constant's avatar
Matthieu Constant committed
148
    # Computation of the new velocities
149
150
151
    forces = fluid.compute_node_force(dt)
    vn     = p.velocity() + forces * dt / p.mass()
    vmax   = np.max(np.hypot(vn[:, 0], vn[:, 1]))
Matthieu Constant's avatar
Matthieu Constant committed
152
    # number of contact sub time step
153
    nsub   = max(1, int(np.ceil((vmax * dt * 4)/min(p.r()))))
Matthieu Constant's avatar
Matthieu Constant committed
154
    print("NSUB", nsub,"VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(p.r()))
Matthieu Constant's avatar
Matthieu Constant committed
155
    # Contact solver
Matthieu Constant's avatar
Matthieu Constant committed
156
157
    for i in range(nsub) :
        p.iterate(dt/nsub, forces)  
158
159

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

Matthieu Constant's avatar
Matthieu Constant committed
162
    # Output files writing
Matthieu Constant's avatar
Matthieu Constant committed
163
164
165
166
    if ii %outf == 0 :
        ioutput = int(ii/outf) + 1
        p.write_vtk(outputdir, ioutput, t)
        fluid.export_vtk(outputdir, t, ioutput)
167
    ii    += 1
Matthieu Constant's avatar
Matthieu Constant committed
168
    print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))