2DiagDrops.py 5.97 KB
Newer Older
Matthieu Constant's avatar
Matthieu Constant committed
1
#!/usr/bin/env python
Matthieu Constant's avatar
Matthieu Constant committed
2
3
from marblesbag import fluid as fluid
from marblesbag import scontact2
Matthieu Constant's avatar
Matthieu Constant committed
4
5
6
7
8
9
10

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

Matthieu Constant's avatar
Matthieu Constant committed
11
12
13
14
15
16
17
18
#Physical parameters for the drops are the ones presented by Machu et al. (2001) "Coalescence, torus formation and breakup of sedimenting drops: experiments and computer simulations"

#grains creation and placement + physical boundaries definition
'''filename is the name of the output file
r is the radius of the grains
rout is the radius of the drop
rhop is the particles density
compacity is the solid volume fraction inside the drop'''
Matthieu Constant's avatar
Matthieu Constant committed
19
def genInitialPosition(filename, r, rout, rhop, compacity) :
Matthieu Constant's avatar
Matthieu Constant committed
20
    p = scontact2.ParticleProblem()
Matthieu Constant's avatar
Matthieu Constant committed
21
    #Loading of the mesh.msh file specifying physical boundaries name
Matthieu Constant's avatar
Matthieu Constant committed
22
    p.load_msh_boundaries("mesh.msh", ["Top", "Bottom","Lateral"])
Matthieu Constant's avatar
Matthieu Constant committed
23
    #Space between the grains to obtain the expected compacity
Matthieu Constant's avatar
Matthieu Constant committed
24
25
    N = compacity*4*rout**2/(np.pi*r**2)
    e = 2*rout/(N)**(1./2.)
Matthieu Constant's avatar
Matthieu Constant committed
26
    #Definition of the points where the grains are located
Matthieu Constant's avatar
Matthieu Constant committed
27
28
29
30
31
32
    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
33
    #First drop
Matthieu Constant's avatar
Matthieu Constant committed
34
35
    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
36
    #Second drop
Matthieu Constant's avatar
Matthieu Constant committed
37
38
    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
39
    #Vertical shift of the grains to the top of the box
Matthieu Constant's avatar
Matthieu Constant committed
40
    p.position()[:, 1] += 0.22
Matthieu Constant's avatar
Matthieu Constant committed
41
42
43
    p.write_vtk(filename,0,0)


Matthieu Constant's avatar
Matthieu Constant committed
44
outputdir = "output"
Matthieu Constant's avatar
Matthieu Constant committed
45
46
47
48
49
50
if not os.path.isdir(outputdir) :
    os.makedirs(outputdir)

t = 0
ii = 0

Matthieu Constant's avatar
Matthieu Constant committed
51
#Enable computation of the initial velocity
Matthieu Constant's avatar
Matthieu Constant committed
52
initialize = 1
Matthieu Constant's avatar
Matthieu Constant committed
53
54

#physical parameters
Matthieu Constant's avatar
Matthieu Constant committed
55
56
57
58
59
60
61
62
63
64
65
g =  -9.81                                      # gravity
rhop = 2400                                     # grains density
nu = 1e-4                                       # kinematic viscosity
drho = 35.4                                     # density differences rhop-rho
compacity = 0.03                                # solid volume fraction in the drop
rho = rhop-drho/compacity                       # fluid density
r=25e-6                                         # grains radii
R = 2.7e-3                                      # drop radius
mu = nu*rho                                     # dynamic viscosity
tEnd = 1000                                     # final time
print("RHOP = %g" % rhop)
Matthieu Constant's avatar
Matthieu Constant committed
66
67

#numerical parameters
Matthieu Constant's avatar
Matthieu Constant committed
68
69
70
71
lcmin = 0.0008                                  # mesh size
dt = 5e-2                                       # time step
alpha = 1e-3                                    # stabilization coefficient
epsilon = alpha*lcmin**2 /nu                    # stabilization parametre
Matthieu Constant's avatar
Matthieu Constant committed
72
73
74
75
print('epsilon',epsilon)

shutil.copy("mesh.msh", outputdir +"/mesh.msh")

Matthieu Constant's avatar
Matthieu Constant committed
76
#Object particles creation
Matthieu Constant's avatar
Matthieu Constant committed
77
78
genInitialPosition(outputdir, r, R, rhop,compacity)
p = scontact2.ParticleProblem()
Matthieu Constant's avatar
Matthieu Constant committed
79
80
p.read_vtk(outputdir,0)
print("r = %g, m = %g\n" %  (p.r()[0], p.mass()[0]))
Matthieu Constant's avatar
Matthieu Constant committed
81
outf = 1                                        # number of iterations between output files
Matthieu Constant's avatar
Matthieu Constant committed
82

Matthieu Constant's avatar
Matthieu Constant committed
83
84
#Object fluid creation + Boundary condition of the fluid (field 0 is horizontal velocity; field 1 is vertical velocity; field 2 is pressure)
#Format: strong_boundaries = [(Boundary tag, Fluid field, Value)
Matthieu Constant's avatar
Matthieu Constant committed
85
strong_boundaries = [("Top",2,0.),("Top",1,0.),("Bottom",1,0.),("Lateral",0,0.)]
Matthieu Constant's avatar
Matthieu Constant committed
86
87
88
89
90
fluid = fluid.fluid_problem("mesh.msh",g,nu*rho,rho,epsilon,strong_boundaries)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())

ii = 0
jj = 0
Matthieu Constant's avatar
Matthieu Constant committed
91
92
tEnd1 = .2

Matthieu Constant's avatar
Matthieu Constant committed
93
#Initial fluid velocities are computed based on the forces applied by the grains without moving the grains
Matthieu Constant's avatar
Matthieu Constant committed
94
95
if initialize:
    for jj in range(4):
Matthieu Constant's avatar
Matthieu Constant committed
96
        dt1 = dt/100
Matthieu Constant's avatar
Matthieu Constant committed
97
98
        t = 0
        if jj!=0:
Matthieu Constant's avatar
Matthieu Constant committed
99
        #Adaptation of the mesh. Args are minimal mesh radius, maximal mesh radius and number of elements
Matthieu Constant's avatar
Matthieu Constant committed
100
            fluid.adapt_mesh(5e-3,8e-4,50000)
Matthieu Constant's avatar
Matthieu Constant committed
101
102
        fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
        
Matthieu Constant's avatar
Matthieu Constant committed
103
104
        while t < tEnd1 :
            #Fluid solver
Matthieu Constant's avatar
Matthieu Constant committed
105
106
            niter = fluid.implicit_euler(dt1)
            forces = fluid.compute_node_force(dt1)
Matthieu Constant's avatar
Matthieu Constant committed
107
            #Changes the velocities of the grains without moving the grains
Matthieu Constant's avatar
Matthieu Constant committed
108
            p.velocity()[:] += dt1*forces/p.mass()
Matthieu Constant's avatar
Matthieu Constant committed
109
            fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
Matthieu Constant's avatar
Matthieu Constant committed
110
            #Increase time step if convergence
Matthieu Constant's avatar
Matthieu Constant committed
111
112
            if niter < 5 and dt1 < 20:
                dt1 *= 3
Matthieu Constant's avatar
Matthieu Constant committed
113
            elif niter > 5 and dt1 > dt/100:
Matthieu Constant's avatar
Matthieu Constant committed
114
115
116
117
118
119
120
                dt1 /= 3
            print(dt1)
            t += dt1
            ii += 1

ii = 0
t = 0
Matthieu Constant's avatar
Matthieu Constant committed
121
122


Matthieu Constant's avatar
Matthieu Constant committed
123
fluid.export_vtk(outputdir,0,0)
Matthieu Constant's avatar
Matthieu Constant committed
124
125
126

ii = 0
tic = time.clock()
Matthieu Constant's avatar
Matthieu Constant committed
127
#Computation loop
Matthieu Constant's avatar
Matthieu Constant committed
128
129
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
while t < tEnd : 
Matthieu Constant's avatar
Matthieu Constant committed
130
    #Adaptation of the mesh. Args are minimal mesh radius, maximal mesh radius and number of elements
Matthieu Constant's avatar
Matthieu Constant committed
131
    if (ii%5==0 and ii != 0):
Matthieu Constant's avatar
Matthieu Constant committed
132
133
134
       fluid.adapt_mesh(5e-3,8e-4,50000)
       fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
    else :
Matthieu Constant's avatar
Matthieu Constant committed
135
       fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
Matthieu Constant's avatar
Matthieu Constant committed
136
    #Fluid solver
Matthieu Constant's avatar
Matthieu Constant committed
137
138
    fluid.implicit_euler(dt)
    forces = fluid.compute_node_force(dt)
Matthieu Constant's avatar
Matthieu Constant committed
139
    #Computation of the new velocities
Matthieu Constant's avatar
Matthieu Constant committed
140
141
    vn = p.velocity() + forces * dt / p.mass()
    vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
Matthieu Constant's avatar
Matthieu Constant committed
142
143
    #number of sub time step
    nsub = max(1, int(np.ceil((vmax * dt * 4)/min(p.r()))))
Matthieu Constant's avatar
Matthieu Constant committed
144
    print("NSUB", nsub,"VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(p.r()))
Matthieu Constant's avatar
Matthieu Constant committed
145
    #Contact solver
Matthieu Constant's avatar
Matthieu Constant committed
146
147
148
    for i in range(nsub) :
        p.iterate(dt/nsub, forces)  
    t += dt
Matthieu Constant's avatar
Matthieu Constant committed
149
    #Output files writting
Matthieu Constant's avatar
Matthieu Constant committed
150
151
152
153
154
155
    if ii %outf == 0 :
        ioutput = int(ii/outf) + 1
        p.write_vtk(outputdir, ioutput, t)
        fluid.export_vtk(outputdir, t, ioutput)
    ii += 1
    print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.clock() - tic))