avalanch1fluidnofriction.py 4.2 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
# MigFlow - Copyright (C) <2010-2018>
# <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
26
27
28

# TESTCASE DESCRIPTION:
# Collapsing of a grain column in fluid. 
# This test case cannot be runned on its own, it requires the output files of the depot test case (../depot/dep.py) to create a compact column as initial situation
# Positions, masses and radii of grains are setted in the depot.py file.
from migflow import fluid
Matthieu Constant's avatar
Matthieu Constant committed
29
from migflow import scontact
30
31
32
33
34
35
36
37
38
39
40

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

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

41
42
43
44
45
# Physical parameters
g = -9.81                                 # gravity
rhop = 1500                               # grains density
rhof = 1000                               # fluid density
nuf = 1e-6                                # kinematic fluid viscosity
46

47
48
# Numerical parameters
outf = 10                                 # number of iterations between output files
Matthieu Constant's avatar
Matthieu Constant committed
49
dt = 1e-4                                 # time step
50
tEnd = 10                                 # final time
51

Matthieu Constant's avatar
Matthieu Constant committed
52
53
54
55
# 
# PARTICLES PROBLEM
#
# Particles container structure builder to load file created with depot.py.
Matthieu Constant's avatar
Matthieu Constant committed
56
p1 = scontact.ParticleProblem(2)
Matthieu Constant's avatar
Matthieu Constant committed
57
# Use particles deposit computed by the depot.py file of "depot" folder.
58
p1.read_vtk("../depot/output",25)
Matthieu Constant's avatar
Matthieu Constant committed
59
# New particles container structure builder to set other boundary conditions than in depot.py.
Matthieu Constant's avatar
Matthieu Constant committed
60
p = scontact.ParticleProblem(2)
61
62
63
64
p.load_msh_boundaries("mesh.msh", ["Top", "Bottom", "Left", "Right"])
for i in range(len(p1.position()[:,0])):
    p.add_particle((p1.position()[i,0], p1.position()[i,1]), p1.r()[i], p1.r()[i]**2 * np.pi * rhop)
p.write_vtk(outputdir, 0, 0)
Matthieu Constant's avatar
Matthieu Constant committed
65
66

# Initial time and iteration.
67
68
t = 0
ii = 0
69

Matthieu Constant's avatar
Matthieu Constant committed
70
71
72
73
#
# FLUID PROBLEM
#
# Fluid structure builder.
Matthieu Constant's avatar
Matthieu Constant committed
74
fluid = fluid.FluidProblem(2,g,[nuf*rhof],[rhof])
Matthieu Constant's avatar
Matthieu Constant committed
75
# Set the mesh geometry for the fluid computation.
76
fluid.load_msh("mesh.msh")
Nathan Coppin's avatar
Nathan Coppin committed
77
78
79
80
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
81
# Set locations of the grains in the mesh and compute the porosity in each computation cell.
82
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
Matthieu Constant's avatar
Matthieu Constant committed
83
# Write output files for post-visualization.
84
fluid.export_vtk(outputdir,0,0)
Matthieu Constant's avatar
Matthieu Constant committed
85

86
tic = time.time()
87
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
Matthieu Constant's avatar
Matthieu Constant committed
88
89
90
91

# 
# COMPUTATION LOOP
# 
92
while t < tEnd : 
Matthieu Constant's avatar
Matthieu Constant committed
93
    # Solve averaged Navier-Stokes equations.
94
95
    fluid.implicit_euler(dt)

Matthieu Constant's avatar
Matthieu Constant committed
96
    # Adaptation of the mesh.
97
    if (ii%15==0 and ii != 0):
98
       fluid.adapt_mesh(1e-2, 3e-3, 10000)
99

Matthieu Constant's avatar
Matthieu Constant committed
100
    # Computation of the new velocities.
101
    forces = fluid.compute_node_force(dt)
102
103
    vn = p.velocity() + forces*dt/p.mass()
    vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
Matthieu Constant's avatar
Matthieu Constant committed
104
    # Number of contact sub time step.
105
    nsub = max(1, int(np.ceil((vmax * dt * 4)/min(p.r()))))
106
    print("NSUB", nsub,"VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(p.r()))
Matthieu Constant's avatar
Matthieu Constant committed
107
    # NLGS iterations.
108
109
110
    for i in range(nsub) :
        p.iterate(dt/nsub, forces)  

111
    t += dt
112
113
    fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())

Matthieu Constant's avatar
Matthieu Constant committed
114
    # Output files writing.
115
    if ii%outf == 0 :
116
117
118
        ioutput = int(ii/outf) + 1
        p.write_vtk(outputdir, ioutput, t)
        fluid.export_vtk(outputdir, t, ioutput)
119
    ii += 1
Matthieu Constant's avatar
Matthieu Constant committed
120
    print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))