script figures intracting drop 3d

testcases/ImagesArticle/InteractingDrop/Diag/2DiagDrops.py

+# Marblesbag - Copyright (C) <2010-2018>
+# <Universite catholique de Louvain (UCL), Belgium
+#  Universite de Montpellier, France>
+# 	
+# List of the contributors to the development of Marblesbag: see AUTHORS file.
+# Description and complete License: see LICENSE file.
+# 	
+# This program (Marblesbag) 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
+from marblesbag import fluid3 as fluid
+from marblesbag import scontact3
+
+import numpy as np
+import os
+import time
+import shutil
+import random
+
+#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'''
+def genInitialPosition(filename, r, rout, rhop, compacity) :
+    p = scontact3.ParticleProblem()
+    #Loading of the mesh.msh file specifying physical boundaries name
+    p.load_msh_boundaries("mesh.msh", ["Top", "Bottom","X","Z"])
+    #Space between the grains to obtain the expected compacity
+    N = compacity*rout**3/(r**3)
+    e = 2*rout/(6*N/np.pi)**(1./3.)
+    #Definition of the points where the grains are located
+    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
+                #First drop
+                if R2<(rout-r)**2:
+                    p.add_particle((x+random.uniform(-e/2+r,e/2-r)+rout, 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 * rhop);
+                    
+    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
+                #Second drop
+                if R2<(rout-r)**2:
+                    p.add_particle((x+random.uniform(-e/2+r,e/2-r)-rout, 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 * rhop);
+    #Vertical shift of the grains to the top of the box
+    p.position()[:, 1] += 0.22
+    print("Number of grains=%d" % len(p.position()[:, 1] ))
+    p.write_vtk(filename,0,0)
+
+
+outputdir = "outputIm"
+if not os.path.isdir(outputdir) :
+    os.makedirs(outputdir)
+
+t = 0
+ii = 0
+#Enable computation of the initial velocity
+initialize = 1
+
+#physical parameters
+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)
+
+#numerical parameters
+lcmin = 0.0004                                  # mesh size
+dt = 5e-2                                       # time step
+alpha = 5e-4                                    # stabilization coefficient
+epsilon = alpha*lcmin**2 /nu                    # stabilization parametre
+print('epsilon',epsilon)
+
+shutil.copy("mesh.msh", outputdir +"/mesh.msh")
+
+#Object particle creation
+genInitialPosition(outputdir, r, R, rhop, compacity)
+p = scontact3.ParticleProblem()
+p.read_vtk(outputdir,0)
+
+print("r = %g, m = %g\n" %  (p.r()[0], p.mass()[0]))
+outf = 10                                       # number of iterations between output files
+
+#Object fluid creation + Boundary condition of the fluid (field 0 is horizontal-width velocity; field 1 is vertical velocity; field 2 is horizontal-depth velocity; field 3 is pressure)
+#Format: strong_boundaries = [(Boundary tag, Fluid field, Value)
+strong_boundaries = [("Top",1,0.),("Bottom",1,0.),("X",0,.0),("Z",2,0.),("Top",3,0.)]
+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
+tEnd1 = 100
+
+#Initial fluid velocities are computed based on the forces applied by the grains without moving the grains
+if initialize:
+    for jj in range(4):
+        dt1 = dt
+        t = 0
+        if jj!=0:
+        #Adaptation of the mesh. Args are minimal mesh radius, maximal mesh radius and number of elements
+            fluid.adapt_mesh(5e-3,8e-4,600000)
+        fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
+        
+        while t < tEnd1 : 
+            #Fluid solver
+            niter = fluid.implicit_euler(dt1)
+            forces = fluid.compute_node_force(dt1)
+            #Changes the velocities of the grains without moving the grains
+            p.velocity()[:] += dt1*forces/p.mass()
+            fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
+            #Increase time step if convergence
+            if niter < 5 and dt1 < 20:
+                dt1 *= 3
+            elif niter > 5 and dt1 > dt:
+                dt1 /= 3
+            print(dt1)
+            t += dt1
+            ii += 1
+
+ii = 0
+t = 0
+tic = time.clock()
+fluid.export_vtk(outputdir,0,0)
+#Computation loop
+fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
+while t < tEnd : 
+#Adaptation of the mesh. Args are minimal mesh radius, maximal mesh radius and number of elements
+    if (ii%5==0 and ii != 0):
+       fluid.adapt_mesh(5e-3,8e-4,600000)
+       fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
+    #Fluid solver
+    fluid.implicit_euler(dt)
+    forces = fluid.compute_node_force(dt)
+    #Computation of the new velocities
+    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()))))
+    print("NSUB", nsub,"VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(p.r()))
+    #Contact solver
+    for i in range(nsub) :
+        p.iterate(dt/nsub, forces)  
+    fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
+    t += dt
+    #Output files writting
+    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))

testcases/ImagesArticle/InteractingDrop/Diag/test.py

+from math import *
+import matplotlib
+from matplotlib.tri import Triangulation
+from matplotlib import pyplot
+from matplotlib.mlab import griddata
+import matplotlib.cm as cm
+#!/usr/bin/env python
+from marblesbag import scontact3
+
+import numpy as np
+import os
+import time
+import shutil
+import random
+
+
+p = scontact3.ParticleProblem()
+p.read_vtk("outputIm",20)
+
+
+
+i = 70
+R = 2.7e-3
+r = []
+x = []
+y = []
+
+b1 = False
+b2 = False
+
+rhop=2400
+
+n = len(p.mass())
+print(n)
+print(p.position()[0,1])
+
+
+fig, ax = pyplot.subplots()
+
+for i in range(n) :
+    if i<n/2.:
+        if p.position()[i,2]<.00002 and p.position()[i,2]>-.00002:
+            circle = pyplot.Circle((p.position()[i,0]/R, p.position()[i,1]/R), 2*(3*p.mass()[i]/(4*np.pi*rhop))**(1./3.)/R, edgecolor='r', facecolor=(0,0,1),linewidth=0.1)
+            ax.add_artist(circle)
+    else:
+        if p.position()[i,2]<.00002 and p.position()[i,2]>-.00002:
+            circle = pyplot.Circle((p.position()[i,0]/R, p.position()[i,1]/R), 2*(3*p.mass()[i]/(4*np.pi*rhop))**(1./3.)/R, edgecolor='b', facecolor=(1,0,0),linewidth=0.1)
+            ax.add_artist(circle)
+
+
+
+ax.set_xlim((-1.5, 3.5))
+ax.set_ylim((.17/R, .18/R))
+ax.set_aspect('equal')
+ax.get_xaxis().set_visible(False)
+ax.get_yaxis().set_visible(False)
+fig.savefig('state_3d_0.png',bbox_inches='tight', dpi=900)

testcases/ImagesArticle/InteractingDrop/Vert/2VertDrops.py

+# Marblesbag - Copyright (C) <2010-2018>
+# <Universite catholique de Louvain (UCL), Belgium
+#  Universite de Montpellier, France>
+# 	
+# List of the contributors to the development of Marblesbag: see AUTHORS file.
+# Description and complete License: see LICENSE file.
+# 	
+# This program (Marblesbag) 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
+from marblesbag import fluid3 as fluid
+from marblesbag import scontact3
+
+import numpy as np
+import os
+import time
+import shutil
+import random
+#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'''
+def genInitialPosition(filename, r, rout, rhop, compacity) :
+    p = scontact3.ParticleProblem()
+    #Loading of the mesh.msh file specifying physical boundaries name
+    p.load_msh_boundaries("mesh.msh", ["Top", "Bottom","X","Z"])
+    #Space between the grains to obtain the expected compacity
+    N = compacity*rout**3/(r**3)
+    e = 2*rout/(6*N/np.pi)**(1./3.)
+    #Definition of the points where the grains are located
+    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
+                #First drop
+                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 * rhop);
+                    
+    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
+                #Second drop
+                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 * rhop);
+    #Vertical shift of the grains to the top of the box
+    p.position()[:, 1] += 0.22
+    print("Number of grains=%d" % len(p.position()[:, 1] ))
+    p.write_vtk(filename,0,0)
+
+
+outputdir = "outputIm"
+if not os.path.isdir(outputdir) :
+    os.makedirs(outputdir)
+
+t = 0
+ii = 0
+
+#Enable computation of the initial velocity
+initialize = 1
+
+#physical parameters
+g =  -9.81                                      # gravity
+rhop = 2400                                     # grains density
+nu = 1e-4                                       # kinematic viscosity
+compacity = 0.013                               # solid volume fraction in the drop
+drho = 15.6                                     # density differences rhop-rho
+rho = rhop-drho/compacity                       # fluid density
+r=25e-6                                         # grains radii
+mu = nu*rho                                     # dynamic viscosity
+R = 2e-3                                        # drop radius
+tEnd = 1000                                     # final time
+
+#numerical parameters
+lcmin = 0.0004                                  # mesh size
+dt = 5e-2                                       # time step
+alpha = 1e-3                                    # stabilization coefficient
+epsilon = alpha*lcmin**2 /nu                    # stabilization parametre
+print('epsilon',epsilon)
+
+shutil.copy("mesh.msh", outputdir +"/mesh.msh")
+
+#Object particles creation
+genInitialPosition(outputdir, r, R, rhop, compacity)
+p = scontact3.ParticleProblem()
+p.read_vtk(outputdir,0)
+
+print("r = %g, m = %g\n" %  (p.r()[0], p.mass()[0]))
+print("RHOP = %g" % rhop)
+outf = 10                                       # number of iterations between output files
+
+#Object fluid creation + Boundary condition of the fluid (field 0 is horizontal-width velocity; field 1 is vertical velocity; field 2 is horizontal-depth velocity; field 3 is pressure)
+#Format: strong_boundaries = [(Boundary tag, Fluid field, Value)
+strong_boundaries = [("Top",1,0.),("Bottom",1,0.),("X",0,.0),("Z",2,0.),("Top",3,0.)]
+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
+tEnd1 = 100
+
+#Initial fluid velocities are computed based on the forces applied by the grains without moving the grains
+if initialize:
+    for jj in range(4):
+        dt1 = dt
+        t = 0
+        if jj!=0:
+        #Adaptation of the mesh. Args are minimal mesh radius, maximal mesh radius and number of elements
+            fluid.adapt_mesh(5e-3,8e-4,600000)
+        fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
+        
+        while t < tEnd1 : 
+            #Fluid solver
+            niter = fluid.implicit_euler(dt1)
+            forces = fluid.compute_node_force(dt1)
+            #Changes the velocities of the grains without moving the grains
+            p.velocity()[:] += dt1*forces/p.mass()
+            fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
+            #Increase time step if convergence
+            if niter < 5 and dt1 < 20:
+                dt1 *= 3
+            elif niter > 5 and dt1 > dt:
+                dt1 /= 3
+            print(dt1)
+            t += dt1
+            ii += 1
+
+ii = 0
+t = 0
+tic = time.clock()
+#Computation loop
+fluid.export_vtk(outputdir,0,0)
+fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
+while t < tEnd : 
+    #Adaptation of the mesh. Args are minimal mesh radius, maximal mesh radius and number of elements
+    if (ii%5==0 and ii != 0):
+       fluid.adapt_mesh(5e-3,8e-4,600000)
+       fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
+    #Fluid solver
+    fluid.implicit_euler(dt)
+    forces = fluid.compute_node_force(dt)
+    #Computation of the new velocities
+    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()))))
+    print("NSUB", nsub,"VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(p.r()))
+    #Contact solver
+    for i in range(nsub) :
+        p.iterate(dt/nsub, forces)  
+    fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
+    t += dt
+    #Output files writting
+    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))

testcases/ImagesArticle/InteractingDrop/Vert/test.py

+from math import *
+import matplotlib
+from matplotlib.tri import Triangulation
+from matplotlib import pyplot
+from matplotlib.mlab import griddata
+import matplotlib.cm as cm
+#!/usr/bin/env python
+from marblesbag import scontact3
+
+import numpy as np
+import os
+import time
+import shutil
+import random
+
+
+p = scontact3.ParticleProblem()
+p.read_vtk("outputIm",50)
+
+
+
+i = 70
+R = 2e-3
+r = []
+x = []
+y = []
+
+b1 = False
+b2 = False
+
+rhop=2400
+
+n = len(p.mass())
+print(n)
+print(p.position()[0,1])
+
+
+fig, ax = pyplot.subplots()
+
+for i in range(n) :
+    if i<n/2.:
+        if p.position()[i,2]<.0001 and p.position()[i,2]>-.0001:
+            circle = pyplot.Circle((p.position()[i,0]/R, p.position()[i,1]/R), 2*(3*p.mass()[i]/(4*np.pi*rhop))**(1./3.)/R, edgecolor='r', facecolor=(0,0,1),linewidth=0.1)
+            ax.add_artist(circle)
+    else:
+        if p.position()[i,2]<.0001 and p.position()[i,2]>-.0001:
+            circle = pyplot.Circle((p.position()[i,0]/R, p.position()[i,1]/R), 2*(3*p.mass()[i]/(4*np.pi*rhop))**(1./3.)/R, edgecolor='b', facecolor=(1,0,0),linewidth=0.1)
+            ax.add_artist(circle)
+
+
+
+ax.set_xlim((-2.5, 2.5))
+ax.set_ylim((.181/R, .191/R))
+ax.set_aspect('equal')
+ax.get_xaxis().set_visible(False)
+ax.get_yaxis().set_visible(False)
+fig.savefig('stateVert_3d_0.png',bbox_inches='tight', dpi=900)