2VertDrops.py 6.72 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
31
# TESTCASE DESCRIPTION
# This tescase presents the interaction between two interacting Stokes clouds.
# Clouds are made of particles 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"

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

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

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

Matthieu Constant's avatar
Matthieu Constant committed
82
# Enable computation of the initial velocity
Matthieu Constant's avatar
Matthieu Constant committed
83
init = True
Matthieu Constant's avatar
Matthieu Constant committed
84
85
86
87
88
89
90
91
92
93
94

# Physical parameters
g = -9.81                                        # gravity
rhop = 2400                                      # particles density
r = 25e-6                                        # particles 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
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
98
99
100
# Numerical parameters
outf = 1                                         # number of iterations between output files
dt = 5e-2                                        # time step
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
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
126
127
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())

Matthieu Constant's avatar
Matthieu Constant committed
128
129
# Initial fluid velocities are computed based on the forces applied by the particles without moving the particles
# The initial mesh is obtained by adapting for times the mesh on the converged fields
Matthieu Constant's avatar
Matthieu Constant committed
130
if init:
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
131
    for jj in range(4):
Matthieu Constant's avatar
Matthieu Constant committed
132
133
        dt1 = dt/100
        t = 0
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
134
        if jj!=0:
Matthieu Constant's avatar
Matthieu Constant committed
135
136
            # Adaptation of the mesh.
            fluid.adapt_mesh(5e-3,8e-4,50000)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
137
138
139
        fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
        
        while t < tEnd1 :
Matthieu Constant's avatar
Matthieu Constant committed
140
141
142
143
            # Fluid solver
            niter = fluid.implicit_euler(dt1)
            forces = fluid.compute_node_force(dt1)
            # Changes the velocities of the particles without moving the particles
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
144
145
            p.velocity()[:] += dt1*forces/p.mass()
            fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
Matthieu Constant's avatar
Matthieu Constant committed
146
            # Increase time step if convergence
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
147
            if niter < 5 and dt1 < 20:
Matthieu Constant's avatar
Matthieu Constant committed
148
149
150
                dt1 *= 3
            elif niter > 5 and dt1 > dt/100:
                dt1 /= 3
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
151
            print(dt1)
Matthieu Constant's avatar
Matthieu Constant committed
152
            t += dt1
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
153

Matthieu Constant's avatar
Matthieu Constant committed
154
t = 0
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
155
fluid.export_vtk(outputdir,0,0)
Matthieu Constant's avatar
Matthieu Constant committed
156
tic = time.time()
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
157
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
Matthieu Constant's avatar
Matthieu Constant committed
158
159
160
161

#
# COMPUTATION LOOP
#
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
162
while t < tEnd : 
Matthieu Constant's avatar
Matthieu Constant committed
163
    # Fluid solver
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
164
    fluid.implicit_euler(dt)
Matthieu Constant's avatar
Matthieu Constant committed
165

Matthieu Constant's avatar
Matthieu Constant committed
166
    # Adaptation of the mesh.
Matthieu Constant's avatar
Matthieu Constant committed
167
    if (ii%15==0 and ii != 0):
Matthieu Constant's avatar
Matthieu Constant committed
168
       fluid.adapt_mesh(5e-3,8e-4,50000)
Matthieu Constant's avatar
Matthieu Constant committed
169
    
Matthieu Constant's avatar
Matthieu Constant committed
170
    # Computation of the new velocities
Matthieu Constant's avatar
Matthieu Constant committed
171
    forces = fluid.compute_node_force(dt)
Matthieu Constant's avatar
Matthieu Constant committed
172
173
174
175
    vn = p.velocity() + forces * dt / p.mass()
    vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
    # Number of sub time step
    nsub = max(1, int(np.ceil((vmax * dt * 4)/min(p.r()))))
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
176
    print("NSUB", nsub,"VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(p.r()))
Matthieu Constant's avatar
Matthieu Constant committed
177
    # NLGS iterations
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
178
179
    for i in range(nsub) :
        p.iterate(dt/nsub, forces)  
Matthieu Constant's avatar
Matthieu Constant committed
180

Matthieu Constant's avatar
Matthieu Constant committed
181
    t += dt
Matthieu Constant's avatar
Matthieu Constant committed
182
183
    fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())

Matthieu Constant's avatar
Matthieu Constant committed
184
    # Output files writting
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
185
186
187
188
189
    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
190
    print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))