betonlmgc.py 4.37 KB
Newer Older
Matthieu Constant's avatar
Matthieu Constant committed
1
#!/usr/bin/env python
2

Matthieu Constant's avatar
Matthieu Constant committed
3
#variable to use lmgc or not
Matthieu Constant's avatar
Matthieu Constant committed
4
use_lmgc = 1
5
6


Matthieu Constant's avatar
Matthieu Constant committed
7
8
from marblesbag import fluid as fluid
from marblesbag import scontact2
9
10
if use_lmgc :
    from marblesbag import lmgc2Interface
Matthieu Constant's avatar
Matthieu Constant committed
11
12
13
14
15
16
17
18
19
20

import numpy as np
import os
import time
import shutil

outputdir = "output"
if not os.path.isdir(outputdir) :
    os.makedirs(outputdir)

Matthieu Constant's avatar
Matthieu Constant committed
21
22
23
24
#grains creation and placement + physical boundaries definition
'''r is the radius of the grains
rout is the outer radius of the geometry
rhop is the particles density'''
Matthieu Constant's avatar
Matthieu Constant committed
25
26
def genInitialPosition(filename, r, rout, rhop) :
    p = scontact2.ParticleProblem()
Matthieu Constant's avatar
Matthieu Constant committed
27
    #Loading of the mesh.msh file specifying physical boundaries name
Matthieu Constant's avatar
Matthieu Constant committed
28
    p.load_msh_boundaries("mesh.msh", ["Outer", "Inner"])
Matthieu Constant's avatar
Matthieu Constant committed
29
    #Definition of the points where the grains are located
Matthieu Constant's avatar
Matthieu Constant committed
30
31
32
33
    x = np.arange(rout, -rout, 2 * -r)
    y = np.arange(-rout, -2*rout/3., 2 * r)
    x, y = np.meshgrid(x, y)
    R2 = x**2 + y**2
Matthieu Constant's avatar
Matthieu Constant committed
34
    #condition to be inside the outer boundary
Matthieu Constant's avatar
Matthieu Constant committed
35
36
37
    keep = R2 < (rout - r)**2
    x = x[keep]
    y = y[keep]
Matthieu Constant's avatar
Matthieu Constant committed
38
    #First layer of grains
Matthieu Constant's avatar
Matthieu Constant committed
39
40
41
    for i in range(x.shape[0]) :
        p.add_particle((x[i], y[i]), r, r**2 * np.pi * rhop);
    
Matthieu Constant's avatar
Matthieu Constant committed
42
43
    x = np.arange(rout, -rout, -r)
    y = np.arange(-2*rout/3.+r, -rout/2.,  r)
Matthieu Constant's avatar
Matthieu Constant committed
44
45
    x, y = np.meshgrid(x, y)
    R2 = x**2 + y**2
46
    keep = R2 < (rout - r)**2
Matthieu Constant's avatar
Matthieu Constant committed
47
48
    x = x[keep]
    y = y[keep]
Matthieu Constant's avatar
Matthieu Constant committed
49
    #Second layer of grains
Matthieu Constant's avatar
Matthieu Constant committed
50
    for i in range(x.shape[0]) :
Matthieu Constant's avatar
Matthieu Constant committed
51
        p.add_particle((x[i], y[i]), r/2., r**2/4. * np.pi * rhop);
Matthieu Constant's avatar
Matthieu Constant committed
52
        
Matthieu Constant's avatar
Matthieu Constant committed
53
54
55
56
57
58
59
    x = np.arange(rout, -rout, 2 * -r/3)
    y = np.arange(-rout/2.+r, -5*rout/12., 2 * r/3)
    x, y = np.meshgrid(x, y)
    R2 = x**2 + y**2
    keep = R2 < (rout - r/3)**2
    x = x[keep]
    y = y[keep]
Matthieu Constant's avatar
Matthieu Constant committed
60
    #Third layer of grains
Matthieu Constant's avatar
Matthieu Constant committed
61
62
63
    for i in range(x.shape[0]) :
        p.add_particle((x[i], y[i]), r/3, r**2 /9 * np.pi * rhop);
    p.write_vtk(filename,0,0)
Matthieu Constant's avatar
Matthieu Constant committed
64
65
66
67
68
69


t = 0
ii = 0

#physical parameters
Matthieu Constant's avatar
Matthieu Constant committed
70
71
72
73
74
g =  -9.81                              # gravity
rho = 1e3                               # fluid density
rhop = 1600                             # grains density
nu = 1e-4                               # kinematic viscosity
tEnd = 100                              # final time
Matthieu Constant's avatar
Matthieu Constant committed
75
76

#numerical parameters
Matthieu Constant's avatar
Matthieu Constant committed
77
78
79
h = 0.002                               # mesh size
dt = 5e-3                               # time step
epsilon = 1e-4                          # stabilization parametre
Frédéric Dubois's avatar
Frédéric Dubois committed
80

Matthieu Constant's avatar
Matthieu Constant committed
81
82
83
#Geometrical parametres
rout = 0.0254                           # outer boundary
r = 397e-6/2                            # grains radius
Matthieu Constant's avatar
Matthieu Constant committed
84
85
shutil.copy("mesh.msh", outputdir +"/mesh.msh")

Matthieu Constant's avatar
Matthieu Constant committed
86
#Object particles creation
Matthieu Constant's avatar
Matthieu Constant committed
87
genInitialPosition(outputdir, r, rout, rhop)
Matthieu Constant's avatar
Matthieu Constant committed
88
if use_lmgc :
Matthieu Constant's avatar
Matthieu Constant committed
89
    friction=0.1                        # friction coefficient
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
90
    lmgc2Interface.scontactToLmgc90(outputdir,0,friction)
Matthieu Constant's avatar
Matthieu Constant committed
91
92
    p = lmgc2Interface.LmgcInterface()
else :
Matthieu Constant's avatar
Matthieu Constant committed
93
94
    p = scontact2.ParticleProblem()
    p.read_vtk(outputdir,0)
Matthieu Constant's avatar
Matthieu Constant committed
95

Matthieu Constant's avatar
Matthieu Constant committed
96
97
98
outf = 5                                # number of iterations between output files
#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
99
strong_boundaries = [("Outer",0,lambda x : -4*x[:, 1]),("Outer",1,lambda x : 4*x[:, 0]),("Inner",0,0.),("Inner",1,0.),("PtFix",2,0.)]
Matthieu Constant's avatar
Matthieu Constant committed
100
fluid = fluid.fluid_problem("mesh.msh",g,nu*rho,rho,epsilon,strong_boundaries)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
101
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
Matthieu Constant's avatar
Matthieu Constant committed
102
103
104
fluid.export(outputdir,0,0)

tic = time.clock()
Matthieu Constant's avatar
Matthieu Constant committed
105
106
#Computation loop
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
Matthieu Constant's avatar
Matthieu Constant committed
107
while t < tEnd :
Matthieu Constant's avatar
Matthieu Constant committed
108
109
    #Fluid solver
    fluid.implicit_euler(dt)
Matthieu Constant's avatar
Matthieu Constant committed
110
    forces = fluid.compute_node_force(dt)
Matthieu Constant's avatar
Matthieu Constant committed
111
    #Computation of the new velocities
Matthieu Constant's avatar
Matthieu Constant committed
112
113
    vn = p.velocity() + forces * dt / p.mass()
    vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
Matthieu Constant's avatar
Matthieu Constant committed
114
    #number of sub time step
Matthieu Constant's avatar
Matthieu Constant committed
115
116
    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()))
Matthieu Constant's avatar
Matthieu Constant committed
117
    #Contact solver
Matthieu Constant's avatar
Matthieu Constant committed
118
    for i in range(nsub) :
Matthieu Constant's avatar
Matthieu Constant committed
119
120
121
        tol = 1e-6
        p.iterate(dt/nsub, forces, tol)
    #Localisation of the grains in the fluid
Matthieu Constant's avatar
Matthieu Constant committed
122
123
    fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
    t += dt
Matthieu Constant's avatar
Matthieu Constant committed
124
    #Output files writting
Matthieu Constant's avatar
Matthieu Constant committed
125
126
    if ii %outf == 0 :
        ioutput = int(ii/outf) + 1
Matthieu Constant's avatar
Matthieu Constant committed
127
        p.write_vtk(outputdir, ioutput, t)
128
        fluid.export_vtk(outputdir, t, ioutput)
Matthieu Constant's avatar
Matthieu Constant committed
129
130
    ii += 1
    print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.clock() - tic))