melangeur.py 5.17 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
25
26
27
28
29
30

# TESTCASE DESCRIPTION
# Grains in circular shear flow.
# This example shows how to set boundary conditions as a function of some parametres.
# A boolean parametre give the possibility to use lmgc90 to solve contacts or not.

from migflow import fluid
from migflow import scontact
31
32
33
34
35

import numpy as np
import os
import time
import shutil
Matthieu Constant's avatar
Matthieu Constant committed
36
import random
37

Matthieu Constant's avatar
Matthieu Constant committed
38
# Define output directory
39
40
41
42
outputdir = "output"
if not os.path.isdir(outputdir) :
    os.makedirs(outputdir)

Matthieu Constant's avatar
Matthieu Constant committed
43
def genInitialPosition(filename, r, rout, rin, rhop) :
Matthieu Constant's avatar
Matthieu Constant committed
44
45
46
47
48
49
50
51
52
53
54
    """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
55
    p.load_msh_boundaries("mesh.msh", ["Outer", "Inner"])
Matthieu Constant's avatar
Matthieu Constant committed
56
    
Matthieu Constant's avatar
Matthieu Constant committed
57
58
59
60
61
62
63
64
65
66
67
68
69
70
    # Definition of the points where the grains are located
    x = np.arange(rout, -rout, 2 * -r)
    x, y = np.meshgrid(x, x)
    R2 = x**2 + y**2
    # Condition to be inside the outer boundary
    keep = R2 < (rout-r)**2
    x = x[keep]
    y = y[keep]
    R2 = x**2 + y**2
    # Condition to be outside the inner boundary
    keep = R2 > (rin+r)**2
    x = x[keep]
    y = y[keep]
    # Add a grain at each centre position
71
    for i in range(x.shape[0]) :
Matthieu Constant's avatar
Matthieu Constant committed
72
73
        if y[i]<0:
            rhop1 = random.choice([rhop*.9,1.1*rhop,rhop])
Matthieu Constant's avatar
Matthieu Constant committed
74
            p.add_particle((x[i], y[i]), r, r**2 * np.pi * rhop1)
Matthieu Constant's avatar
Matthieu Constant committed
75
    p.write_vtk(filename,0,0)
Matthieu Constant's avatar
Matthieu Constant committed
76
 
Matthieu Constant's avatar
Matthieu Constant committed
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
# 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
97
genInitialPosition(outputdir, r, rout, rin, rhop)
Matthieu Constant's avatar
Matthieu Constant committed
98
p = scontact.ParticleProblem(2)
Matthieu Constant's avatar
Matthieu Constant committed
99
p.read_vtk(outputdir,0)
100

Matthieu Constant's avatar
Matthieu Constant committed
101
102
103
104
105
106
107
108
109
# Initial time and iteration
t = 0
ii = 0

#
# 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
110
111
112
113
114
115
116
117
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)
fluid.set_weak_boundary("Inner","Wall")
fluid.set_weak_boundary("Outer","Wall")
Matthieu Constant's avatar
Matthieu Constant committed
118
# Set location of the grains in the mesh and compute the porosity in each computation cell
119
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
Matthieu Constant's avatar
Matthieu Constant committed
120
fluid.export_vtk(outputdir,0,0)
121

Matthieu Constant's avatar
Matthieu Constant committed
122
tic = time.time()
Matthieu Constant's avatar
Matthieu Constant committed
123
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
Matthieu Constant's avatar
Matthieu Constant committed
124
125
126
127

# 
# COMPUTATION LOOP
#
Matthieu Constant's avatar
Matthieu Constant committed
128
while t < tEnd :
Matthieu Constant's avatar
Matthieu Constant committed
129
    # Fluid solver
Matthieu Constant's avatar
Matthieu Constant committed
130
    fluid.implicit_euler(dt)
Matthieu Constant's avatar
Matthieu Constant committed
131

Matthieu Constant's avatar
Matthieu Constant committed
132
133
134
135
136
137
    # Computation of the new velocities
    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()))))
138
    print("NSUB", nsub,"VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(p.r()))
Matthieu Constant's avatar
Matthieu Constant committed
139
    # NLGS iterations
140
    for i in range(nsub) :
Matthieu Constant's avatar
Matthieu Constant committed
141
        tol = 1e-6                                # Numerical tolerance for grains intersection
Matthieu Constant's avatar
Matthieu Constant committed
142
        p.iterate(dt/nsub, forces, tol)
Matthieu Constant's avatar
Matthieu Constant committed
143

Matthieu Constant's avatar
Matthieu Constant committed
144
    # Localisation of the grains in the fluid
145
    fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
Matthieu Constant's avatar
Matthieu Constant committed
146
    t += dt
Matthieu Constant's avatar
Matthieu Constant committed
147

Matthieu Constant's avatar
Matthieu Constant committed
148
    # Output files writting
149
150
151
    if ii %outf == 0 :
        ioutput = int(ii/outf) + 1
        p.write_vtk(outputdir, ioutput, t)
Matthieu Constant's avatar
Matthieu Constant committed
152
        fluid.export_vtk(outputdir, t, ioutput)
Matthieu Constant's avatar
Matthieu Constant committed
153
154
    ii += 1
    print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time()-tic))