2DDDrops.py 8.1 KB
Newer Older
Matthieu Constant's avatar
Matthieu Constant committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
# MigFlow - Copyright (C) <2010-2018>
# <Universite catholique de Louvain (UCL), Belgium
#  Universite de Montpellier, France>
# 	
# List of the contributors to the development of MigFlow: see AUTHORS file.
# Description and complete License: see LICENSE file.
# 	
# This program (MigFlow) is free software: 
# you can redistribute it and/or modify it under the terms of the GNU Lesser General 
# 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
# GNU Lesser General Public License for more details.
# 
# 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, 
# see <http://www.gnu.org/licenses/>.

#!/usr/bin/env python
Matthieu Constant's avatar
Matthieu Constant committed
23
24
25
26
27

# TESTCASE DESCRIPTION
# This tescase presents the fall of interacting clouds made of particles in a viscous fluid (Stokes cloud).
# Physical parameters for the drops are the ones presented in Constant et al. (2018)

Matthieu Constant's avatar
Matthieu Constant committed
28
#!/usr/bin/env python
Matthieu Constant's avatar
Matthieu Constant committed
29
30
from migflow import fluid
from migflow import scontact
Matthieu Constant's avatar
Matthieu Constant committed
31
32
33
34
35
36
37
38

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

def genInitialPosition(filename, r, rout, rhop1, rhop2, compacity) :
Matthieu Constant's avatar
Matthieu Constant committed
39
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
    rhop1 -- heavy particles density        
    rhop2 -- light particles density        
    compacity -- initial compacity in the cloud
    """
    # Particles structure builder
    p = scontact.ParticleProblem(3)
    # Load mesh.msh file specifying physical boundaries name
Matthieu Constant's avatar
Matthieu Constant committed
52
53
    p.load_msh_boundaries("mesh.msh", ["Top", "Bottom","X","Z"])
    
Matthieu Constant's avatar
Matthieu Constant committed
54
    # Space between the particles to obtain the expected compacity
Matthieu Constant's avatar
Matthieu Constant committed
55
56
57
    N = compacity*rout**3/(r**3)
    e = 2*rout/(6*N/np.pi)**(1./3.)

Matthieu Constant's avatar
Matthieu Constant committed
58
    # Definition of the points where the particles are located
Matthieu Constant's avatar
Matthieu Constant committed
59
60
61
62
    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
Matthieu Constant's avatar
Matthieu Constant committed
63
                # First drop
Matthieu Constant's avatar
Matthieu Constant committed
64
65
66
67
68
69
70
                if R2<(rout-r)**2:
                    p.add_particle((x+random.uniform(-e/2+r,e/2-r), y+random.uniform(-e/2+r,e/2-r)-6*rout, z+random.uniform(-e/2+r,e/2-r)), r, 4./3.* r**3 * np.pi * rhop2)
                    
    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
Matthieu Constant's avatar
Matthieu Constant committed
71
                # Second drop
Matthieu Constant's avatar
Matthieu Constant committed
72
73
74
75
76
77
78
                if R2<(rout-r)**2:
                    p.add_particle((x+random.uniform(-e/2+r,e/2-r), y+random.uniform(-e/2+r,e/2-r)+6*rout, z+random.uniform(-e/2+r,e/2-r)), r, 4./3.* r**3 * np.pi * rhop1)
    
    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
Matthieu Constant's avatar
Matthieu Constant committed
79
                # Third drop
Matthieu Constant's avatar
Matthieu Constant committed
80
81
82
83
84
85
86
                if R2<(rout-r)**2:
                    p.add_particle((x+random.uniform(-e/2+r,e/2-r), y+random.uniform(-e/2+r,e/2-r)+2*rout, z+random.uniform(-e/2+r,e/2-r)), r, 4./3.* r**3 * np.pi * rhop1)
                    
    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
Matthieu Constant's avatar
Matthieu Constant committed
87
                # Fourth drop
Matthieu Constant's avatar
Matthieu Constant committed
88
89
90
91
                if R2<(rout-r)**2:
                    p.add_particle((x+random.uniform(-e/2+r,e/2-r), y+random.uniform(-e/2+r,e/2-r)-2*rout, z+random.uniform(-e/2+r,e/2-r)), r, 4./3.* r**3 * np.pi * rhop1)
    p.write_vtk(filename,0,0)

Matthieu Constant's avatar
Matthieu Constant committed
92
# Define output directory
Matthieu Constant's avatar
Matthieu Constant committed
93
94
95
96
outputdir = "output"
if not os.path.isdir(outputdir) :
    os.makedirs(outputdir)

Matthieu Constant's avatar
Matthieu Constant committed
97
# Enable computation of the initial velocity
Matthieu Constant's avatar
Matthieu Constant committed
98
init = True
Matthieu Constant's avatar
Matthieu Constant committed
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114

# Physical parameters
g = -9.81                                        # gravity
rhop1 = 1600                                     # heavy particles density
rhop2 = 400                                      # light particles density
r = 25e-6                                        # particles radii
compacity = 0.013/4                              # solid volume fraction in the drop
rho = 1000                                       # fluid density
nu = 1e-6                                        # kinematic viscosity
R = 2e-3                                         # cloud radius
mu = nu*rho                                      # dynamic viscosity

# Numerical parameters
outf = 1                                         # number of iterations between output files
dt = 5e-3                                        # time step
tEnd = 100                                       # final time
Matthieu Constant's avatar
Matthieu Constant committed
115

Matthieu Constant's avatar
Matthieu Constant committed
116
117
118
119
# 
# PARTICLE PROBLEM
#
# Initialise particles
Matthieu Constant's avatar
Matthieu Constant committed
120
genInitialPosition(outputdir, r, R, rhop1, rhop2, compacity)
Matthieu Constant's avatar
Matthieu Constant committed
121
p = scontact.ParticleProblem(3)
Matthieu Constant's avatar
Matthieu Constant committed
122
123
p.read_vtk(outputdir,0)

Matthieu Constant's avatar
Matthieu Constant committed
124
125
126
127
# Initial time and iteration
t = 0
ii = 0
jj = 0
128
tEnd1 = 100
Matthieu Constant's avatar
Matthieu Constant committed
129
130
131
132

#
# FLUID PROBLEM
#
Matthieu Constant's avatar
Matthieu Constant committed
133
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
134
# Set the mesh geometry for the fluid computation
Matthieu Constant's avatar
Matthieu Constant committed
135
fluid.load_msh("mesh.msh")
Matthieu Constant's avatar
Matthieu Constant committed
136
137
138
139
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
140

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

Matthieu Constant's avatar
Matthieu Constant committed
144
145
# 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
146
147
if init:
    for jj in range(4):
148
        dt1 = dt
Matthieu Constant's avatar
Matthieu Constant committed
149
        t = 0
Matthieu Constant's avatar
Matthieu Constant committed
150
        if jj!=0:
Matthieu Constant's avatar
Matthieu Constant committed
151
            # Adaptation of the mesh.
Matthieu Constant's avatar
Matthieu Constant committed
152
153
154
155
            fluid.adapt_mesh(8e-3,8e-4,250000,cmax=.1,cmin=20)
        fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
        
        while t < tEnd1 : 
Matthieu Constant's avatar
Matthieu Constant committed
156
157
158
159
            # Fluid solver
            niter = fluid.implicit_euler(dt1)
            forces = fluid.compute_node_force(dt1)
            # Changes the velocities of the particles without moving the particles
Matthieu Constant's avatar
Matthieu Constant committed
160
161
            p.velocity()[:] += dt1*forces/p.mass()
            fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
Matthieu Constant's avatar
Matthieu Constant committed
162
            # Increase time step if convergence
Matthieu Constant's avatar
Matthieu Constant committed
163
            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:
Matthieu Constant's avatar
Matthieu Constant committed
166
                dt1 /= 3
Matthieu Constant's avatar
Matthieu Constant committed
167
            print(dt1)
Matthieu Constant's avatar
Matthieu Constant committed
168
            t += dt1
Matthieu Constant's avatar
Matthieu Constant committed
169

Matthieu Constant's avatar
Matthieu Constant committed
170
t = 0
Matthieu Constant's avatar
Matthieu Constant committed
171
172
173
fluid.export_vtk(outputdir,0,0)
tic       =  time.time()
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
Matthieu Constant's avatar
Matthieu Constant committed
174
175
176
177

#
# COMPUTATION LOOP
#
Matthieu Constant's avatar
Matthieu Constant committed
178
while t < tEnd : 
Matthieu Constant's avatar
Matthieu Constant committed
179
    # Fluid solver
Matthieu Constant's avatar
Matthieu Constant committed
180
181
    fluid.implicit_euler(dt)

Matthieu Constant's avatar
Matthieu Constant committed
182
    # 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
183
184
185
    if (ii%5==0 and ii != 0):
       fluid.adapt_mesh(8e-3,8e-4,250000,cmin=20,cmax=0.1)

Matthieu Constant's avatar
Matthieu Constant committed
186
    # Computation of the new velocities
Matthieu Constant's avatar
Matthieu Constant committed
187
    forces = fluid.compute_node_force(dt)
Matthieu Constant's avatar
Matthieu Constant committed
188
    vn = p.velocity() + forces * dt / p.mass()
189
    vmax = np.max(np.linalg.norm(vn))
Matthieu Constant's avatar
Matthieu Constant committed
190
191
    # Number of subtime step
    nsub = max(1, int(np.ceil((vmax * dt * 4)/min(p.r()))))
Matthieu Constant's avatar
Matthieu Constant committed
192
    print("NSUB", nsub,"VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(p.r()))
Matthieu Constant's avatar
Matthieu Constant committed
193
    # NLGS iterations
Matthieu Constant's avatar
Matthieu Constant committed
194
195
196
197
    for i in range(nsub) :
        p.iterate(dt/nsub, forces)  
    
    fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
Matthieu Constant's avatar
Matthieu Constant committed
198
    t += dt
Matthieu Constant's avatar
Matthieu Constant committed
199
    
Matthieu Constant's avatar
Matthieu Constant committed
200
    # Output files writting
Matthieu Constant's avatar
Matthieu Constant committed
201
202
203
204
    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
205
    ii += 1
Matthieu Constant's avatar
Matthieu Constant committed
206
    print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))