melangeur.py 5.47 KB
Newer Older
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
1
# MigFlow - Copyright (C) <2010-2018>
2
3
# <Universite catholique de Louvain (UCL), Belgium
#  Universite de Montpellier, France>
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
4
# 	
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
5
# List of the contributors to the development of MigFlow: see AUTHORS file.
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
6
7
# Description and complete License: see LICENSE file.
# 	
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
8
# This program (MigFlow) is free software: 
9
# you can redistribute it and/or modify it under the terms of the GNU Lesser General 
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
10
11
12
13
14
15
# 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
16
# GNU Lesser General Public License for more details.
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
17
# 
18
19
# 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, 
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
20
21
# see <http://www.gnu.org/licenses/>.

22
#!/usr/bin/env python
Matthieu Constant's avatar
Matthieu Constant committed
23
24

# TESTCASE DESCRIPTION
Frédéric Dubois's avatar
Frédéric Dubois committed
25
26
27
28
# Grains in circular shear flow generated by a rotating drum.
# This example shows how to set boundary conditions as a function of some parameters.

# A boolean parameter gives the possibility to use lmgc90 to solve contacts instead of scontact.
Nathan Coppin's avatar
Nathan Coppin committed
29
use_lmgc90 = True
Matthieu Constant's avatar
Matthieu Constant committed
30
31
32

from migflow import fluid
from migflow import scontact
Frédéric Dubois's avatar
Frédéric Dubois committed
33
34
if use_lmgc90 :
    from migflow import lmgc90Interface
35
36
37
38
39

import numpy as np
import os
import time
import shutil
Matthieu Constant's avatar
Matthieu Constant committed
40
import random
41

Matthieu Constant's avatar
Matthieu Constant committed
42
def genInitialPosition(filename, r, rout, rin, rhop) :
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
    r -- max particle radius
    rout -- outer radius of the drum
    rhop -- particle density
    """
    # Grains structure builder
    p = scontact.ParticleProblem(2)
    # Load mesh.msh file specifying physical boundaries names
54
    p.load_msh_boundaries("mesh.msh", ["Outer", "Inner"])
Nathan Coppin's avatar
Nathan Coppin committed
55
    
Matthieu Constant's avatar
Matthieu Constant committed
56
57
    # Definition of the points where the grains are located
    x = np.arange(rout, -rout, 2 * -r)
58
59
    x, y = np.meshgrid(x, x)
    R2 = x**2 + y**2
Matthieu Constant's avatar
Matthieu Constant committed
60
61
    # Condition to be inside the outer boundary
    keep = R2 < (rout-r)**2
62
63
    x = x[keep]
    y = y[keep]
Matthieu Constant's avatar
Matthieu Constant committed
64
    R2 = x**2 + y**2
Matthieu Constant's avatar
Matthieu Constant committed
65
66
    # Condition to be outside the inner boundary
    keep = R2 > (rin+r)**2
Matthieu Constant's avatar
Matthieu Constant committed
67
68
    x = x[keep]
    y = y[keep]
Matthieu Constant's avatar
Matthieu Constant committed
69
    # Add a grain at each centre position
70
    for i in range(x.shape[0]) :
Matthieu Constant's avatar
Matthieu Constant committed
71
72
        if y[i]<0:
            rhop1 = random.choice([rhop*.9,1.1*rhop,rhop])
Nathan Coppin's avatar
Nathan Coppin committed
73
            p.add_particle((x[i], y[i]), r, r**2 * np.pi * rhop1)
Matthieu Constant's avatar
Matthieu Constant committed
74
    p.write_vtk(filename,0,0)
Matthieu Constant's avatar
Matthieu Constant committed
75
 
76
77


Frédéric Dubois's avatar
Frédéric Dubois committed
78
79
80
81
82
# Define output directory
outputdir = "output"
if not os.path.isdir(outputdir) :
    os.makedirs(outputdir)

Matthieu Constant's avatar
Matthieu Constant committed
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
# Physical parameters
g = 0                                          # gravity
rho = 1.253e3                                  # fluid density
rhop = 1000                                    # grains density
nu = 1e-4                                      # kinematic viscosity

# Numerical parameters
dt = 1e-2                                      # time step
tEnd = 50                                      # final time

# Geometry parameters
outf = 5                                       # number of iterations between output files
rout = 0.0254                                  # outer radius
rin = 0.0064                                   # inner radius
r = 397e-6/2                                   # grains radius

#
# PARTICLE PROBLEM
#
# Object particles creation
Matthieu Constant's avatar
Matthieu Constant committed
103
genInitialPosition(outputdir, r, rout, rin, rhop)
Frédéric Dubois's avatar
Frédéric Dubois committed
104
105
106
107
108
109
110
if use_lmgc90 :
    friction = 0.1                                    # friction coefficient
    lmgc90Interface.scontactTolmgc90(outputdir,2,0,friction)
    p = lmgc90Interface.ParticleProblem(2)
else :
  p = scontact.ParticleProblem(2)
  p.read_vtk(outputdir,0)
111

Matthieu Constant's avatar
Matthieu Constant committed
112
113
114
# Initial time and iteration
t = 0
ii = 0
115

Matthieu Constant's avatar
Matthieu Constant committed
116
117
118
119
120
#
# FLUID PROBLEM
#
fluid = fluid.FluidProblem(2,g,[nu*rho],[rho])
# Set the mesh geometry for the fluid computation
Matthieu Constant's avatar
Matthieu Constant committed
121
122
123
124
125
126
fluid.load_msh("mesh.msh")
fluid.set_strong_boundary("Outer",0,lambda x : -x[:, 1])
fluid.set_strong_boundary("Outer",1,lambda x : x[:, 0])
fluid.set_strong_boundary("Inner",0,0)
fluid.set_strong_boundary("Inner",1,0)
fluid.set_strong_boundary("PtFix",2,0)
Nathan Coppin's avatar
Nathan Coppin committed
127
128
fluid.set_wall_boundary("Inner")
fluid.set_wall_boundary("Outer")
Matthieu Constant's avatar
Matthieu Constant committed
129
# Set location of the grains in the mesh and compute the porosity in each computation cell
130
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
Matthieu Constant's avatar
Matthieu Constant committed
131
fluid.export_vtk(outputdir,0,0)
132

Matthieu Constant's avatar
Matthieu Constant committed
133
tic = time.time()
Matthieu Constant's avatar
Matthieu Constant committed
134
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
Matthieu Constant's avatar
Matthieu Constant committed
135
136
137
138

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

Matthieu Constant's avatar
Matthieu Constant committed
143
    # Computation of the new velocities
144
    forces = fluid.compute_node_force(dt)
Matthieu Constant's avatar
Matthieu Constant committed
145
    vn = p.velocity() + forces*dt/p.mass()
146
    vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
Matthieu Constant's avatar
Matthieu Constant committed
147
148
    # Number of contact sub time step
    nsub = max(1, int(np.ceil((vmax*dt*4)/min(p.r()))))
149
    print("NSUB", nsub,"VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(p.r()))
Matthieu Constant's avatar
Matthieu Constant committed
150
    # NLGS iterations
151
    for i in range(nsub) :
Matthieu Constant's avatar
Matthieu Constant committed
152
        tol = 1e-6                                # Numerical tolerance for grains intersection
Matthieu Constant's avatar
Matthieu Constant committed
153
        p.iterate(dt/nsub, forces, tol)
Matthieu Constant's avatar
Matthieu Constant committed
154

Matthieu Constant's avatar
Matthieu Constant committed
155
    # Localisation of the grains in the fluid
156
157
    fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
    t += dt
Matthieu Constant's avatar
Matthieu Constant committed
158

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