drop.py 6.68 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
# 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"

28
from migflow import fluid
Matthieu Constant's avatar
Matthieu Constant committed
29
from migflow import scontact
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
30
31
32
33
34
35

import numpy as np
import os
import time
import shutil
import random
Matthieu Constant's avatar
Matthieu Constant committed
36

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
37
def genInitialPosition(filename, r, rout, rhop, compacity) :
Matthieu Constant's avatar
Matthieu Constant committed
38
    """Set all the particles centre positions and create the particles objects to add in the computing structure
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
39
    
Matthieu Constant's avatar
Matthieu Constant committed
40
41
42
43
44
45
46
47
48
49
50
51
52
    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(3)
    # Load mesh.msh file specifying physical boundaries name
    p.load_msh_boundaries("mesh.msh", ["Top", "Bottom","X","Z"])

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

Matthieu Constant's avatar
Matthieu Constant committed
56
    # Definition of the points where the particles are located
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
57
58
59
60
61
    for x in np.arange(rout, -rout, -e):
        for y in np.arange(rout, -rout, -e):
            for z in np.arange(rout, -rout, -e):
                R2 = x**2 + y**2 + z**2
                if R2<(rout-r)**2:
Matthieu Constant's avatar
Matthieu Constant committed
62
                    # Add particles object with properties defined above
Matthieu Constant's avatar
Matthieu Constant committed
63
64
                    p.add_particle((x, y, z), r, 4./3.* r**3 * np.pi * rhop)
    
Matthieu Constant's avatar
Matthieu Constant committed
65
    # Vertical shift of the particles to the top of the box
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
66
67
68
    p.position()[:, 1] += 0.52
    p.write_vtk(filename,0,0)

Matthieu Constant's avatar
Matthieu Constant committed
69
# Define output directory
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
70
71
72
73
outputdir = "output"
if not os.path.isdir(outputdir) :
    os.makedirs(outputdir)

Matthieu Constant's avatar
Matthieu Constant committed
74
# Enable computation of the initial velocity
Matthieu Constant's avatar
Matthieu Constant committed
75
init = True
Matthieu Constant's avatar
Matthieu Constant committed
76
77
78
79
80
81
82
83
84
85

# Physical parameters
g = -9.81                                       # gravity
rhop = 2450                                     # particles density
r = 154e-6                                      # particles radii
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
86
print("RHOP = %g, r = %g, RHO = %g" % (rhop,r,rho))
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
87

Matthieu Constant's avatar
Matthieu Constant committed
88
# Numerical parameters
89
outf = 10                                       # number of iterations between output files
Matthieu Constant's avatar
Matthieu Constant committed
90
dt = 5e-2                                       # time step
Matthieu Constant's avatar
Matthieu Constant committed
91
tEnd = 1000                                      # final time
Matthieu Constant's avatar
Matthieu Constant committed
92

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

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

#
# FLUID PROBLEM
#
Matthieu Constant's avatar
Matthieu Constant committed
110
fluid = fluid.FluidProblem(3,g,[nu*rho],[rho],petsc_solver_type="-pc_type ilu -ksp_max_it 30")
Matthieu Constant's avatar
Matthieu Constant committed
111
# Set the mesh geometry for the fluid computation
Matthieu Constant's avatar
Matthieu Constant committed
112
fluid.load_msh("mesh.msh")
Matthieu Constant's avatar
Matthieu Constant committed
113
114
115
116
fluid.set_wall_boundary("Top", pressure=0)
fluid.set_wall_boundary("Bottom")
fluid.set_wall_boundary("X")
fluid.set_wall_boundary("Z")
Matthieu Constant's avatar
Matthieu Constant committed
117

Matthieu Constant's avatar
Matthieu Constant committed
118
# Set location of the particles in the mesh and compute the porosity in each computation cell
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
119
120
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())

Matthieu Constant's avatar
Matthieu Constant committed
121
122
# 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
123
if init:
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
124
    for jj in range(4):
125
        dt1 = dt
Matthieu Constant's avatar
Matthieu Constant committed
126
        t = 0
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
127
        if jj!=0:
Matthieu Constant's avatar
Matthieu Constant committed
128
            # Adaptation of the mesh.
129
            fluid.adapt_mesh(8e-3,8e-4,150000,cmax=.1,cmin=20)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
130
131
132
        fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
        
        while t < tEnd1 : 
Matthieu Constant's avatar
Matthieu Constant committed
133
            # Fluid solver
134
            niter = fluid.implicit_euler(dt1,newton_atol=1e-7,newton_max_it=100)
Matthieu Constant's avatar
Matthieu Constant committed
135
136
            forces = fluid.compute_node_force(dt1)
            # Changes the velocities of the particles without moving the particles
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
137
138
            p.velocity()[:] += dt1*forces/p.mass()
            fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
Matthieu Constant's avatar
Matthieu Constant committed
139
            # Increase time step if convergence
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
140
            if niter < 5 and dt1 < 20:
Matthieu Constant's avatar
Matthieu Constant committed
141
                dt1 *= 3
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
142
            elif niter > 5 and dt1 > dt:
Matthieu Constant's avatar
Matthieu Constant committed
143
                dt1 /= 3
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
144
            print(dt1)
Matthieu Constant's avatar
Matthieu Constant committed
145
            t += dt1
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
146

Matthieu Constant's avatar
Matthieu Constant committed
147
t = 0
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
148
fluid.export_vtk(outputdir,0,0)
Matthieu Constant's avatar
Matthieu Constant committed
149
tic       =  time.time()
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
150
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
Matthieu Constant's avatar
Matthieu Constant committed
151
152
153
154

#
# COMPUTATION LOOP
#
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
155
while t < tEnd : 
Matthieu Constant's avatar
Matthieu Constant committed
156
    # Fluid solver
157
    fluid.implicit_euler(dt,newton_atol=1e-7,newton_max_it=100)
Matthieu Constant's avatar
Matthieu Constant committed
158

Matthieu Constant's avatar
Matthieu Constant committed
159
    # Adaptation of the mesh. Args are maximal mesh radius, minimal mesh radius and number of elements. Optional arguments can be given to weigh the max and min gradient used in the adapt formula
Matthieu Constant's avatar
Matthieu Constant committed
160
    if (ii%5==0 and ii != 0):
161
       fluid.adapt_mesh(8e-3,8e-4,150000,cmin=20,cmax=0.1)
Matthieu Constant's avatar
Matthieu Constant committed
162

Matthieu Constant's avatar
Matthieu Constant committed
163
    # Computation of the new velocities
Matthieu Constant's avatar
Matthieu Constant committed
164
    forces = fluid.compute_node_force(dt)
Matthieu Constant's avatar
Matthieu Constant committed
165
    vn = p.velocity() + forces * dt / p.mass()
166
167
    vmax = np.max(np.linalg.norm(vn))
    
Matthieu Constant's avatar
Matthieu Constant committed
168
169
    # Number of subtime step
    nsub = max(1, int(np.ceil((vmax * dt * 4)/min(p.r()))))
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
170
    print("NSUB", nsub,"VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(p.r()))
Matthieu Constant's avatar
Matthieu Constant committed
171
    # NLGS iterations
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
172
173
    for i in range(nsub) :
        p.iterate(dt/nsub, forces)  
Matthieu Constant's avatar
Matthieu Constant committed
174
    
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
175
    fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
Matthieu Constant's avatar
Matthieu Constant committed
176
    t += dt
Matthieu Constant's avatar
Matthieu Constant committed
177
    
Matthieu Constant's avatar
Matthieu Constant committed
178
    # Output files writting
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
179
180
181
182
    if ii %outf == 0 :
        ioutput = int(ii/outf) + 1
        p.write_vtk(outputdir, ioutput, t)
        fluid.export_vtk(outputdir, t, ioutput)
Matthieu Constant's avatar
Matthieu Constant committed
183
    ii += 1
184
    print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))