drop.py 5.13 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

Matthieu Constant's avatar
Matthieu Constant committed
24
25
26
27
28
29
30
# TESTCASE DESCRIPTION
# This tescase presents the fall of a cloud made of particles in a viscous fluid (Stokes cloud).
# Physical parameters for the drops are the ones presented by Metzger et al. (2007) 
# "Falling clouds of particles in viscous fluids"

from migflow import fluid
from migflow import scontact
31
from migflow import time_integration
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

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
39
def genInitialPosition(filename, r, rout, rhop, compacity) :
Matthieu Constant's avatar
Matthieu Constant committed
40
41
42
43
44
45
46
47
48
49
50
51
    """Set all the particles centre positions and create the particles objects to add in the computing structure
    
    Keyword arguments:    
    filename -- name of the output file
    r -- max radius of the particles
    rout -- outer radius of the cloud
    rhop -- particles density        
    compacity -- initial compacity in the cloud
    """
    # Particles structure builder
    p = scontact.ParticleProblem(2)
    # Load mesh.msh file specifying physical boundaries names
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
52
    p.load_msh_boundaries("mesh.msh", ["Top", "Bottom","Lateral"])
Matthieu Constant's avatar
Matthieu Constant committed
53

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

Matthieu Constant's avatar
Matthieu Constant committed
58
59
60
61
62
63
64
    # Definition of the points where the particles are located
    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]
Matthieu Constant's avatar
Matthieu Constant committed
65

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
66
    for i in range(x.shape[0]) :
Matthieu Constant's avatar
Matthieu Constant committed
67
        # Add particles object with properties defined above
Matthieu Constant's avatar
Matthieu Constant committed
68
        p.add_particle((x[i]+2*random.uniform(-e/2+r,e/2-r), y[i]+2*random.uniform(-e/2+r,e/2-r)), r, r**2 * np.pi * rhop)
Matthieu Constant's avatar
Matthieu Constant committed
69
    
Matthieu Constant's avatar
Matthieu Constant committed
70
    # Vertical shift of the particles to the top of the box
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
71
72
73
    p.position()[:, 1] += 0.52
    p.write_vtk(filename,0,0)

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

Matthieu Constant's avatar
Matthieu Constant committed
79
80
81
# Physical parameters
g = -9.81                                       # gravity
rhop = 2450                                     # particles density
Matthieu Constant's avatar
Matthieu Constant committed
82
r = 154e-6/2.5                                  # particles radii
Matthieu Constant's avatar
Matthieu Constant committed
83
84
85
86
87
compacity = 0.20                                # solid volume fraction in the drop
rho = 1030                                      # fluid density
nu = 1.17/rho                                   # kinematic viscosity
rout = 3.3e-3                                   # cloud radius
mu = nu*rho                                     # dynamic viscosity
Matthieu Constant's avatar
Matthieu Constant committed
88
print("RHOP = %g, r = %g, RHO = %g" % (rhop,r,rho))
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
89

Matthieu Constant's avatar
Matthieu Constant committed
90
91
# Numerical parameters
outf = 1                                        # number of iterations between output files
Jonathan Lambrechts's avatar
heun    
Jonathan Lambrechts committed
92
dt = 2e-2                                       # time step
Matthieu Constant's avatar
Matthieu Constant committed
93
tEnd = 100                                      # final time
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
94

Matthieu Constant's avatar
Matthieu Constant committed
95
96
97
98
#
# PARTICLE PROBLEM
#
# Initialise particles
Matthieu Constant's avatar
Matthieu Constant committed
99
genInitialPosition(outputdir, r, rout, rhop,compacity)
Matthieu Constant's avatar
Matthieu Constant committed
100
p = scontact.ParticleProblem(2)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
101
102
p.read_vtk(outputdir,0)

Matthieu Constant's avatar
Matthieu Constant committed
103
104
105
106
107
108
109
110
111
# Initial time and iteration
t = 0
ii = 0
jj = 0
tEnd1 = 0.2

#
# FLUID PROBLEM
#
Jonathan Lambrechts's avatar
heun    
Jonathan Lambrechts committed
112
fluid = fluid.FluidProblem(2,g,[nu*rho],[rho],drag_in_stab=1)
Matthieu Constant's avatar
Matthieu Constant committed
113
# Set the mesh geometry for the fluid computation
114
fluid.load_msh("mesh.msh")
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
115
fluid.set_wall_boundary("Top", pressure=0)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
116
117
fluid.set_wall_boundary("Bottom")
fluid.set_wall_boundary("Lateral")
Matthieu Constant's avatar
Matthieu Constant committed
118
# Set location of the particles in the mesh and compute the porosity in each computation cell
Matthieu Constant's avatar
Matthieu Constant committed
119
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.contact_forces(),init=True)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
120

Matthieu Constant's avatar
Matthieu Constant committed
121
t = 0
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
122
fluid.solution()[:,2] = (fluid.coordinates()[:,1]-0.6)*rho*g
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
123
fluid.export_vtk(outputdir,0,0)
Matthieu Constant's avatar
Matthieu Constant committed
124
125
tic = time.time()

126
127
G = np.zeros_like(p.velocity())
G[:,1] = p.mass()[:,0]*g
Matthieu Constant's avatar
Matthieu Constant committed
128
129
130
#
# COMPUTATION LOOP
#
Jonathan Lambrechts's avatar
heun    
Jonathan Lambrechts committed
131
while t < tEnd :
Matthieu Constant's avatar
Matthieu Constant committed
132

Matthieu Constant's avatar
Matthieu Constant committed
133
134
135
136
    if ((ii+1)%15==0):
        number_p = fluid.n_particles
        position_p = fluid.particle_position()
        volume_p = fluid.particle_volume()
Jonathan Lambrechts's avatar
heun    
Jonathan Lambrechts committed
137
    if (ii%15==0 and ii != 0):
Matthieu Constant's avatar
Matthieu Constant committed
138
139
        fluid.adapt_mesh(2e-2,1e-3,10000,old_n_particles=number_p,old_particle_position=position_p,old_particle_volume=volume_p)

140
    time_integration.iterate(fluid, p, dt, external_particles_forces=G)
Jonathan Lambrechts's avatar
heun    
Jonathan Lambrechts committed
141
    t += dt
Matthieu Constant's avatar
Matthieu Constant committed
142

Matthieu Constant's avatar
Matthieu Constant committed
143
    # Output files writting
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
144
145
146
147
148
    if ii %outf == 0 :
        ioutput = int(ii/outf) + 1
        p.write_vtk(outputdir, ioutput, t)
        fluid.export_vtk(outputdir, t, ioutput)
    ii += 1
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
149
    print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))