2DiagDrops.py 5.41 KB
Newer Older
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
1
# MigFlow - Copyright (C) <2010-2020>
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
31
# TESTCASE DESCRIPTION
# This tescase presents the interaction between two interacting Stokes clouds.
# Clouds are made of particles and placed one behind another with an horizontal offset
# 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"

from migflow import fluid
from migflow import scontact
32
from migflow import time_integration
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
33
34
35
36
37
38
39
40

import numpy as np
import os
import time
import shutil
import random

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

Matthieu Constant's avatar
Matthieu Constant committed
55
56
57
    # 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
58

Matthieu Constant's avatar
Matthieu Constant committed
59
60
61
62
63
64
65
    # 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
66

Matthieu Constant's avatar
Matthieu Constant committed
67
    # First cloud
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
68
69
    for i in range(x.shape[0]) :
        p.add_particle((x[i]+random.uniform(-e/2+r,e/2-r)+rout, 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
70
    # Second cloud
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
71
    for i in range(x.shape[0]) :
Matthieu Constant's avatar
Matthieu Constant committed
72
        p.add_particle((x[i]+random.uniform(-e/2+r,e/2-r)-rout, 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
73
    
Matthieu Constant's avatar
Matthieu Constant committed
74
    # Vertical shift of the particles to the top of the box
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
75
76
77
    p.position()[:, 1] += 0.22
    p.write_vtk(filename,0,0)

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

Matthieu Constant's avatar
Matthieu Constant committed
83
84

# Physical parameters
Matthieu Constant's avatar
Matthieu Constant committed
85
g = np.array([0,-9.81])                         # gravity
Matthieu Constant's avatar
Matthieu Constant committed
86
87
88
89
90
91
92
93
rhop = 2400                                     # particles density
r = 25e-6                                       # particles radii
compacity = 0.03                                # solid volume fraction in the drop
nu = 1e-4                                       # kinematic viscosity
drho = 35.4                                     # density differences rhop-rho
rho = rhop-drho/compacity                       # fluid density
rout = 2.7e-3                                   # cloud radius
mu = nu*rho                                     # dynamic viscosity
Matthieu Constant's avatar
Matthieu Constant committed
94
H = 0.5
Matthieu Constant's avatar
Matthieu Constant committed
95
print("RHOP = %g, r = %g, RHO = %g" % (rhop,r,rho))
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
96

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

Matthieu Constant's avatar
Matthieu Constant committed
102
103
104
105
#
# PARTICLE PROBLEM
#
# Initialise particles
Matthieu Constant's avatar
Matthieu Constant committed
106
genInitialPosition(outputdir, r, rout, rhop,compacity)
Matthieu Constant's avatar
Matthieu Constant committed
107
p = scontact.ParticleProblem(2)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
108
109
p.read_vtk(outputdir,0)

Matthieu Constant's avatar
Matthieu Constant committed
110
111
112
113
114
115
116
117
118
# Initial time and iteration
t = 0
ii = 0
jj = 0
tEnd1 = 0.2

#
# FLUID PROBLEM
#
Matthieu Constant's avatar
Matthieu Constant committed
119
fluid = fluid.FluidProblem(2,g,[nu*rho],[rho],petsc_solver_type="-pc_type ilu -ksp_max_it 30")
Matthieu Constant's avatar
Matthieu Constant committed
120
# Set the mesh geometry for the fluid computation
Matthieu Constant's avatar
Matthieu Constant committed
121
fluid.load_msh("mesh.msh")
Matthieu Constant's avatar
Matthieu Constant committed
122
123
124
fluid.set_wall_boundary("Top", pressure=0)
fluid.set_wall_boundary("Bottom")
fluid.set_wall_boundary("Lateral")
Matthieu Constant's avatar
Matthieu Constant committed
125
# Set location of the particles in the mesh and compute the porosity in each computation cell
Nathan Coppin's avatar
Nathan Coppin committed
126
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.contact_forces())
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
127

Matthieu Constant's avatar
Matthieu Constant committed
128
t = 0
Nathan Coppin's avatar
Nathan Coppin committed
129
fluid.write_vtk(outputdir,0,0)
Matthieu Constant's avatar
Matthieu Constant committed
130
131
132
133
134
tic = time.time()

#
# COMPUTATION LOOP
#
135
while t < tEnd :
Matthieu Constant's avatar
Matthieu Constant committed
136
137
138
139
    if ((ii+1)%15==0):
        number_p = fluid.n_particles
        position_p = fluid.particle_position()
        volume_p = fluid.particle_volume()
Matthieu Constant's avatar
Matthieu Constant committed
140
    # Adaptation of the mesh.
Matthieu Constant's avatar
Matthieu Constant committed
141
    if (ii%15==0 and ii != 0):
Matthieu Constant's avatar
Matthieu Constant committed
142
       fluid.adapt_mesh(5e-3,8e-4,50000,old_n_particles=number_p,old_particle_position=position_p,old_particle_volume=volume_p)
Nathan Coppin's avatar
Nathan Coppin committed
143
    time_integration.iterate(fluid, p, dt, external_particles_forces=g*p.mass())
144
    t+=dt
Matthieu Constant's avatar
Matthieu Constant committed
145
    # Output files writting
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
146
147
148
    if ii %outf == 0 :
        ioutput = int(ii/outf) + 1
        p.write_vtk(outputdir, ioutput, t)
Nathan Coppin's avatar
Nathan Coppin committed
149
        fluid.write_vtk(outputdir, ioutput, t)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
150
    ii += 1
151
    print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))