2VertDrops.py 8.03 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/>.

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
22
#!/usr/bin/env python
Matthieu Constant's avatar
Matthieu Constant committed
23
24
25
26
27
28
29

'''
This tescase presents the interaction between two interacting Stokes clouds.
Clouds are made of grains and placed one behind another 
Physical parameters for the clouds are the ones presented by Machu et al. (2001) 
"Coalescence, torus formation and breakup of sedimenting drops: experiments and computer simulations"
'''
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
30
31
from migflow import fluid as fluid
from migflow import scontact2
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
32
33
34
35
36
37

import numpy as np
import os
import time
import shutil
import random
Matthieu Constant's avatar
Matthieu Constant committed
38
39
40
41
42
43
44
45
46

'''
genInitialPosition is a function that sets all the grains centre positions and create the grains objects to add in the computing structure
Args:    -filename    : name of the output file
         -r           : max radius of the grains
         -rout        : outer radius of the cloud
         -rhop        : grains density        
         -compacity   : initial compacity in the cloud
'''
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
47
def genInitialPosition(filename, r, rout, rhop, compacity) :
Matthieu Constant's avatar
Matthieu Constant committed
48
49
50
    #Grains structure builder
    p     = scontact2.ParticleProblem()
    #Load mesh.msh file specifying physical boundaries names
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
51
    p.load_msh_boundaries("mesh.msh", ["Top", "Bottom","Lateral"])
Matthieu Constant's avatar
Matthieu Constant committed
52

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
53
    #Space between the grains to obtain the expected compacity
Matthieu Constant's avatar
Matthieu Constant committed
54
55
56
    N     = compacity*4*rout**2/(np.pi*r**2)
    e     = 2*rout/(N)**(1./2.)

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
57
    #Definition of the points where the grains are located
Matthieu Constant's avatar
Matthieu Constant committed
58
59
60
61
62
63
64
65
    x     = np.arange(rout, -rout, -e)
    x, y  = np.meshgrid(x, x)
    R2    = x**2 + y**2
    keep  = R2 < (rout - r)**2
    x     = x[keep]
    y     = y[keep]

    #First cloud
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
66
67
    for i in range(x.shape[0]) :
        p.add_particle((x[i]+random.uniform(-e/2+r,e/2-r), y[i]+random.uniform(-e/2+r,e/2-r)-2*rout), r, r**2 * np.pi * rhop);
Matthieu Constant's avatar
Matthieu Constant committed
68
    #Second cloud
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
69
70
    for i in range(x.shape[0]) :
        p.add_particle((x[i]+random.uniform(-e/2+r,e/2-r), y[i]+random.uniform(-e/2+r,e/2-r)+2*rout), r, r**2 * np.pi * rhop);
Matthieu Constant's avatar
Matthieu Constant committed
71
    
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
72
73
74
75
    #Vertical shift of the grains to the top of the box
    p.position()[:, 1] += 0.22
    p.write_vtk(filename,0,0)

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

#Enable computation of the initial velocity
Matthieu Constant's avatar
Matthieu Constant committed
82
init      = 1
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
83
84

#physical parameters
Matthieu Constant's avatar
Matthieu Constant committed
85
86
87
88
89
90
91
92
93
94
g         = -9.81                                     # gravity
rhop      =  2400                                     # grains density
r         =  25e-6                                    # grains radii
compacity =  0.013                                    # solid volume fraction in the drop
nu        =  1e-4                                     # kinematic viscosity
drho      =  15.6                                     # density differences rhop-rho
rho       =  rhop-drho/compacity                      # fluid density
rout      =  2e-3                                     # cloud radius
mu        =  nu*rho                                   # dynamic viscosity
print("RHOP = %g, r = %g, RHO = %g" % (rhop,r,rho))
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
95
96

#numerical parameters
Matthieu Constant's avatar
Matthieu Constant committed
97
98
99
outf      =  1                                        # number of iterations between output files
dt        =  5e-2                                     # time step
tEnd      =  100                                      # final time
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
100

Matthieu Constant's avatar
Matthieu Constant committed
101
102
#Initialise grains
genInitialPosition(outputdir, r, rout, rhop,compacity)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
103
104
105
p = scontact2.ParticleProblem()
p.read_vtk(outputdir,0)

Matthieu Constant's avatar
Matthieu Constant committed
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
#Initial time and iteration
t         =  0
ii        =  0
jj        =  0
tEnd1     =  0.2

'''
fluid_problem is the function builder of the fluid class
Args:     -g             : gravity
          -viscosity     : list of fluid phases viscosities
          -density       : list of fluid phases densities
'''
fluid = fluid.fluid_problem(g,[nu*rho],[rho])
#Set the mesh geometry for the fluid computation
fluid.load_msh("mesh.msh")
'''
set_strong_boundary is the function setting the strong boundaries (=constraints) for the fluid problem
Args:     -Tag           : physical tag (set in the mesh.geo file) of the boundary on which you want to add a constraint
          -Field         : O: x-velocity; 1: y-velocity; 2: pressure
          -Value         : value assigned to the specified field on the specified boundary
'''
fluid.set_strong_boundary("Top",2,0)
fluid.set_strong_boundary("Top",1,0)
fluid.set_strong_boundary("Bottom",1,0)
fluid.set_strong_boundary("Lateral",0,0)
'''
set_weak_boundary is the function setting the weak boundaries (=normal fluxes) for the fluid problem
Args:     -Tag           : physical tag (set in the mesh.geo file) of the boundary on which you want to add a constraint
          -Bnd Type      : boundaries can be of type: wall, inflow, outflow
'''
fluid.set_weak_boundary("Top","Wall")
fluid.set_weak_boundary("Bottom","Wall")
fluid.set_weak_boundary("Lateral","Wall")

#Set location of the grains in the mesh and compute the porosity in each computation cell
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
141
142
143
144
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())


#Initial fluid velocities are computed based on the forces applied by the grains without moving the grains
Matthieu Constant's avatar
Matthieu Constant committed
145
146
#The initial mesh is obtained by adapting for times the mesh on the converged fields
if init:
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
147
    for jj in range(4):
Matthieu Constant's avatar
Matthieu Constant committed
148
149
        dt1                  = dt/100
        t                    = 0
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
150
151
        if jj!=0:
        #Adaptation of the mesh. Args are minimal mesh radius, maximal mesh radius and number of elements
Matthieu Constant's avatar
Matthieu Constant committed
152
            fluid.adapt_mesh(5e-3,5e-4,50000)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
153
154
155
156
        fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
        
        while t < tEnd1 :
            #Fluid solver
Matthieu Constant's avatar
Matthieu Constant committed
157
158
            niter            = fluid.implicit_euler(dt1)
            forces           = fluid.compute_node_force(dt1)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
159
160
161
162
163
            #Changes the velocities of the grains without moving the grains
            p.velocity()[:] += dt1*forces/p.mass()
            fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
            #Increase time step if convergence
            if niter < 5 and dt1 < 20:
Matthieu Constant's avatar
Matthieu Constant committed
164
                dt1         *= 3
Matthieu Constant's avatar
Matthieu Constant committed
165
            elif niter > 5 and dt1 > dt/1000:
Matthieu Constant's avatar
Matthieu Constant committed
166
                dt1         /= 3
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
167
            print(dt1)
Matthieu Constant's avatar
Matthieu Constant committed
168
            t               += dt1
Matthieu Constant's avatar
Matthieu Constant committed
169
            print(t)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
170

Matthieu Constant's avatar
Matthieu Constant committed
171
t         =  0
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
172
173
fluid.export_vtk(outputdir,0,0)

Matthieu Constant's avatar
Matthieu Constant committed
174
tic       =  time.time()
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
175
176
177
178
179
#Computation loop
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
while t < tEnd : 
    #Fluid solver
    fluid.implicit_euler(dt)
Matthieu Constant's avatar
Matthieu Constant committed
180
181
182

    #Adaptation of the mesh. Args are minimal mesh radius, maximal mesh radius and number of elements
    if (ii%15==0 and ii != 0):
Matthieu Constant's avatar
Matthieu Constant committed
183
       fluid.adapt_mesh(5e-3,5e-4,50000)
Matthieu Constant's avatar
Matthieu Constant committed
184
    
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
185
    #Computation of the new velocities
Matthieu Constant's avatar
Matthieu Constant committed
186
187
188
    forces = fluid.compute_node_force(dt)
    vn     = p.velocity() + forces * dt / p.mass()
    vmax   = np.max(np.hypot(vn[:, 0], vn[:, 1]))
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
189
    #number of sub time step
Matthieu Constant's avatar
Matthieu Constant committed
190
    nsub   = max(1, int(np.ceil((vmax * dt * 4)/min(p.r()))))
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
191
192
193
194
    print("NSUB", nsub,"VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(p.r()))
    #Contact solver
    for i in range(nsub) :
        p.iterate(dt/nsub, forces)  
Matthieu Constant's avatar
Matthieu Constant committed
195
196
197
198

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

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
199
200
201
202
203
204
    #Output files writting
    if ii %outf == 0 :
        ioutput = int(ii/outf) + 1
        p.write_vtk(outputdir, ioutput, t)
        fluid.export_vtk(outputdir, t, ioutput)
    ii += 1
Matthieu Constant's avatar
Matthieu Constant committed
205
    print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))