betonlmgc.py 6.18 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/>.

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

Matthieu Constant's avatar
Matthieu Constant committed
24
25
26
27
28
# TESTCASE DESCRIPTION
# Mixing of grains in a rotating drum.
# 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.

Matthieu Constant's avatar
Matthieu Constant committed
29
#variable to use lmgc or not
Matthieu Constant's avatar
Matthieu Constant committed
30
use_lmgc = 1
Matthieu Constant's avatar
Matthieu Constant committed
31
32
from migflow import fluid
from migflow import scontact
33
if use_lmgc :
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
34
    from migflow import lmgc2Interface
Matthieu Constant's avatar
Matthieu Constant committed
35
36
37
38
39
40
import numpy as np
import os
import time
import shutil

def genInitialPosition(filename, r, rout, rhop) :
Matthieu Constant's avatar
Matthieu Constant committed
41
42
43
44
45
46
47
48
49
50
51
    """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
Matthieu Constant's avatar
Matthieu Constant committed
52
    p.load_msh_boundaries("mesh.msh", ["Outer", "Inner"])
53

Matthieu Constant's avatar
Matthieu Constant committed
54
55
56
57
58
59
60
61
62
63
64
    # First layer of grains
    # Definition of the points where the grains are located
    x = np.arange(rout, -rout, 2 * -r)
    y = np.arange(-rout, -2*rout/3., 2 * r)
    x, y = np.meshgrid(x, y)
    R2 = x**2 + y**2
    # Condition to be inside the outer boundary
    keep = R2 < (rout-r)**2
    x = x[keep]
    y = y[keep]
    # Add a grain at each centre position
Matthieu Constant's avatar
Matthieu Constant committed
65
66
67
    for i in range(x.shape[0]) :
        p.add_particle((x[i], y[i]), r, r**2 * np.pi * rhop);
    
Matthieu Constant's avatar
Matthieu Constant committed
68
69
70
71
72
73
74
75
76
77
78
    # Second layer of grains
    # Definition of the points where the grains are located
    x = np.arange(rout, -rout, -r)
    y = np.arange(-2*rout/3.+r, -rout/2.,  r)
    x, y = np.meshgrid(x, y)
    # Condition to be inside the outer boundary
    R2 = x**2 + y**2
    keep = R2 < (rout-r)**2
    x = x[keep]
    y = y[keep]
    # Add a grain at each centre position
Matthieu Constant's avatar
Matthieu Constant committed
79
    for i in range(x.shape[0]) :
Matthieu Constant's avatar
Matthieu Constant committed
80
        p.add_particle((x[i], y[i]), r/2., r**2/4. * np.pi * rhop);
81
    
Matthieu Constant's avatar
Matthieu Constant committed
82
83
84
85
86
87
88
89
90
91
92
    # Third layer of grains
    # Definition of the points where the grains are located
    x = np.arange(rout, -rout, 2 * -r/3)
    y = np.arange(-rout/2.+r, -5*rout/12., 2 * r/3)
    x, y = np.meshgrid(x, y)
    R2 = x**2 + y**2
    # Condition to be inside the outer boundary
    keep = R2 < (rout-r/3)**2
    x = x[keep]
    y = y[keep]
    # Add a grain at each centre position
Matthieu Constant's avatar
Matthieu Constant committed
93
    for i in range(x.shape[0]) :
Matthieu Constant's avatar
Matthieu Constant committed
94
        p.add_particle((x[i], y[i]), r/3, r**2 /9*np.pi*rhop);
Matthieu Constant's avatar
Matthieu Constant committed
95
    p.write_vtk(filename,0,0)
Matthieu Constant's avatar
Matthieu Constant committed
96

Matthieu Constant's avatar
Matthieu Constant committed
97
# Define output directory
98
99
100
outputdir = "output"
if not os.path.isdir(outputdir) :
    os.makedirs(outputdir)
Matthieu Constant's avatar
Matthieu Constant committed
101

Matthieu Constant's avatar
Matthieu Constant committed
102
103
104
105
106
107
# Physical parameters
g = -9.81                                       # gravity
rhop = 1600                                     # grains density
r = 397e-6/2                                    # grains radius
rho = 1e3                                       # fluid density
nu = 1e-4                                       # kinematic viscosity
Matthieu Constant's avatar
Matthieu Constant committed
108

Matthieu Constant's avatar
Matthieu Constant committed
109
110
111
112
# Numerical parameters
outf = 1                                        # number of iterations between output files
dt = 5e-3                                       # time step
tEnd = 100                                      # final time
Frédéric Dubois's avatar
Frédéric Dubois committed
113

Matthieu Constant's avatar
Matthieu Constant committed
114
115
# Geometrical parametres
rout = 0.0254                                   # outer boundary
Matthieu Constant's avatar
Matthieu Constant committed
116

Matthieu Constant's avatar
Matthieu Constant committed
117
118
119
120
# 
# PARTICLE PROBLEM
#
# Initialise grains
Matthieu Constant's avatar
Matthieu Constant committed
121
genInitialPosition(outputdir, r, rout, rhop)
Matthieu Constant's avatar
Matthieu Constant committed
122
if use_lmgc :
123
    friction = 0.1                                    # friction coefficient
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
124
    lmgc2Interface.scontactToLmgc90(outputdir,0,friction)
Matthieu Constant's avatar
Matthieu Constant committed
125
    p = lmgc2Interface.LmgcInterface()
Matthieu Constant's avatar
Matthieu Constant committed
126
else :
Matthieu Constant's avatar
Matthieu Constant committed
127
    p = scontact.ParticleProblem(2)
Matthieu Constant's avatar
Matthieu Constant committed
128
    p.read_vtk(outputdir,0)
Matthieu Constant's avatar
Matthieu Constant committed
129

Matthieu Constant's avatar
Matthieu Constant committed
130
131
132
133
134
135
136
137
# Initial time and iteration
t = 0
ii = 0

#
# FLUID PROBLEM
#
fluid = fluid.FluidProblem(2,g,[nu*rho],[rho])
138
139
140
141
142
143
144
145
146
147
#Set the mesh geometry for the fluid computation
fluid.load_msh("mesh.msh")
fluid.set_strong_boundary("Outer",0,lambda x : -4*x[:, 1])
fluid.set_strong_boundary("Outer",1,lambda x : 4*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")
#Set location of the grains in the mesh and compute the porosity in each computation cell
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
148
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
149
fluid.export_vtk(outputdir,0,0)
Matthieu Constant's avatar
Matthieu Constant committed
150

Matthieu Constant's avatar
Matthieu Constant committed
151
tic = time.time()
Matthieu Constant's avatar
Matthieu Constant committed
152
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
Matthieu Constant's avatar
Matthieu Constant committed
153
154
155
156

#
# COMPUTATION LOOP
#
Matthieu Constant's avatar
Matthieu Constant committed
157
while t < tEnd :
Matthieu Constant's avatar
Matthieu Constant committed
158
    # Fluid solver
Matthieu Constant's avatar
Matthieu Constant committed
159
    fluid.implicit_euler(dt)
160

Matthieu Constant's avatar
Matthieu Constant committed
161
162
163
164
165
166
    # 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()))))
Matthieu Constant's avatar
Matthieu Constant committed
167
    print("NSUB", nsub,"VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(p.r()))
Matthieu Constant's avatar
Matthieu Constant committed
168
    # NLGS iterations
Matthieu Constant's avatar
Matthieu Constant committed
169
    for i in range(nsub) :
Matthieu Constant's avatar
Matthieu Constant committed
170
        tol = 1e-6                                #numerical tolerance for grains intersection
Matthieu Constant's avatar
Matthieu Constant committed
171
        p.iterate(dt/nsub, forces, tol)
172

Matthieu Constant's avatar
Matthieu Constant committed
173
    # Localisation of the grains in the fluid
Matthieu Constant's avatar
Matthieu Constant committed
174
    fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
Matthieu Constant's avatar
Matthieu Constant committed
175
    t += dt
176

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