depot.py 7.29 KB
Newer Older
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
1
# MigFlow - Copyright (C) <2010-2020>
Matthieu Constant's avatar
Matthieu Constant committed
2
3
# <Universite catholique de Louvain (UCL), Belgium
#  Universite de Montpellier, France>
Michel Henry's avatar
Michel Henry committed
4
#
Matthieu Constant's avatar
Matthieu Constant committed
5
6
# List of the contributors to the development of MigFlow: see AUTHORS file.
# Description and complete License: see LICENSE file.
Michel Henry's avatar
Michel Henry committed
7
8
9
#
# This program (MigFlow) is free software:
# you can redistribute it and/or modify it under the terms of the GNU Lesser General
Matthieu Constant's avatar
Matthieu Constant committed
10
11
# Public License as published by the Free Software Foundation, either version
# 3 of the License, or (at your option) any later version.
Michel Henry's avatar
Michel Henry committed
12
#
Matthieu Constant's avatar
Matthieu Constant committed
13
14
15
16
# 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.
Michel Henry's avatar
Michel Henry committed
17
#
Matthieu Constant's avatar
Matthieu Constant committed
18
# You should have received a copy of the GNU Lesser General Public License
Michel Henry's avatar
Michel Henry committed
19
# along with this program (see COPYING and COPYING.LESSER files).  If not,
Nathan Coppin's avatar
Nathan Coppin committed
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
# MigFlow - Copyright (C) <2010-2020>
# <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,
Matthieu Constant's avatar
Matthieu Constant committed
39
40
41
42
43
44
45
46
47
48
# see <http://www.gnu.org/licenses/>.

#!/usr/bin/env python

# TESTCASE DESCRIPTION
# Bidimensional particles sedimentation in fluid


from migflow import fluid
from migflow import scontact
49
from migflow import time_integration
Matthieu Constant's avatar
Matthieu Constant committed
50
51
52
53

import numpy as np
import os
import shutil
Nathan Coppin's avatar
Nathan Coppin committed
54
55
import matplotlib.pyplot as plt
import time
Matthieu Constant's avatar
Matthieu Constant committed
56
57
import random

Michel Henry's avatar
Michel Henry committed
58
def genInitialPosition(filename, r, H, ly, lx, rhop) :
Matthieu Constant's avatar
Matthieu Constant committed
59
    """Set all the particles centre positions and create the particles objects to add in the computing structure
Michel Henry's avatar
Michel Henry committed
60
61

    Keyword arguments:
Matthieu Constant's avatar
Matthieu Constant committed
62
63
64
65
66
    filename -- name of the output file
    r -- max radius of the particles
    H -- domain height
    ly - particles area height
    lx -- particles area width
Michel Henry's avatar
Michel Henry committed
67
    rhop -- particles density
Matthieu Constant's avatar
Matthieu Constant committed
68
69
    """
    # Particles structure builder
Michel Henry's avatar
Michel Henry committed
70
    p = scontact.ParticleProblem(2)
Nathan Coppin's avatar
Nathan Coppin committed
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97

    p.add_boundary_segment((-0.07,-0.21),(-0.07,0.23/5),0,"Steel")
    p.add_boundary_segment((-0.07,0.23/5),(0.07,0.23/5),0,"Steel")
    p.add_boundary_segment((0.07,0.23/5),(0.07,-0.21),0,"Steel")
    p.add_boundary_segment((0.07,-0.21),(-0.07,-0.21),0,"Steel")
    R = np.ones((8))*r
    angle = 5
    Dx = 2*r*np.cos(np.deg2rad(angle))
    Dy = 2*r*np.sin(np.deg2rad(angle))
    pos = np.array([[s*Dx,s*Dy] for s in range(8)])
    material = np.array(["Sand" for i in range(8)],dtype=object)
    pos2 = np.copy(pos)
    R2 = np.copy(R)
    mat2 = np.copy(material)
    for i in range (1):
      pos = np.concatenate((pos,pos2+np.array([-(i+1)*Dy,(i+1)*Dx])))
      R = np.concatenate((R,R2))
      material = np.concatenate((material,mat2))
    #R = np.concatenate((R,np.array([r,r,r,r,r,r])))
    #material = np.concatenate((material,np.array(["Sand","Sand","Sand","Sand","Sand","Sand"],dtype=object)))
    #pos3 = np.array([[s*Dx,s*Dy+0.013] for s in [0,15]])
    #for i in range (3): 
    #  pos = np.concatenate((pos,pos3+np.array([-(i+1)*Dy,(i+1)*Dx])))
    #for i in range(1,4):
    #  for j in range(1,15):
    #    material[16*i + j] = "ignore"
    p.add_particle_body(pos,R,np.pi*rhop*R**2,material)
Matthieu Constant's avatar
Matthieu Constant committed
98
99
100
    p.write_vtk(filename,0,0)

# Define output directory
Nathan Coppin's avatar
Nathan Coppin committed
101
outputdir = "outputLowDragHighI*"
Matthieu Constant's avatar
Matthieu Constant committed
102
103
104
105
if not os.path.isdir(outputdir) :
    os.makedirs(outputdir)

# Physical parameters
Matthieu Constant's avatar
Matthieu Constant committed
106
g = np.array([0,-9.81])                         # gravity
Nathan Coppin's avatar
Nathan Coppin committed
107
108
r = 2.15e-3/4                                   # particles radius
rhop = 1840                                     # particles density
Matthieu Constant's avatar
Matthieu Constant committed
109
110
111
112
rho = 1000                                      # fluid density
nu = 1e-6                                       # kinematic viscosity

# Numerical parameters
Nathan Coppin's avatar
Nathan Coppin committed
113
114
115
outf = 1                                        # number of iterations between output files
dt = 1e-3                                       # time step
tEnd = 2.5                                        # final time
Matthieu Constant's avatar
Matthieu Constant committed
116
117

# Geometrical parameters
Matthieu Constant's avatar
Matthieu Constant committed
118
119
ly = 1e-1                                       # particles area height
lx = 4e-1                                       # particles area widht
120
H = 0.6                                         # domain height
Matthieu Constant's avatar
Matthieu Constant committed
121
122
123
124
125

#
# PARTICLE PROBLEM
#
# Initialise particles
Michel Henry's avatar
Michel Henry committed
126
genInitialPosition(outputdir, r, H, ly, lx, rhop)
Matthieu Constant's avatar
Matthieu Constant committed
127
p = scontact.ParticleProblem(2)
Michel Henry's avatar
Michel Henry committed
128
p.read_vtk(outputdir,0)
Matthieu Constant's avatar
Matthieu Constant committed
129

Nathan Coppin's avatar
Nathan Coppin committed
130
p.set_friction_coefficient(0.,"Sand","Steel")
Matthieu Constant's avatar
Matthieu Constant committed
131
132
133
# Initial time and iteration
t         =  0
ii        =  0
Nathan Coppin's avatar
Nathan Coppin committed
134
135
p.body_invert_mass()[4] = 1/(rhop*6.92e-3*1.25*2.15e-3)
p.body_invert_inertia()[4] = 12*p.body_invert_mass()[4]/(6.92e-3**2+(1.25*2.15e-3)**2)
Matthieu Constant's avatar
Matthieu Constant committed
136
137
138
139
140
#
# FLUID PROBLEM
#
fluid = fluid.FluidProblem(2,g,[nu*rho],[rho])
# Set the mesh geometry for the fluid computation
Nathan Coppin's avatar
Nathan Coppin committed
141
fluid.load_msh("rect.msh")
Matthieu Constant's avatar
Matthieu Constant committed
142
143
144
fluid.set_wall_boundary("Bottom")
fluid.set_wall_boundary("Lateral")
fluid.set_wall_boundary("Top",pressure=0)
Matthieu Constant's avatar
Matthieu Constant committed
145
# Set location of the particles in the mesh and compute the porosity in each computation cell
Nathan Coppin's avatar
Nathan Coppin committed
146
fluid.set_particles(p.delassus(), p.volume(), p.position(), p.velocity(), np.zeros_like(p.velocity()))
Nathan Coppin's avatar
Nathan Coppin committed
147
fluid.write_vtk(outputdir,0,0)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
148
p.write_vtk(outputdir, 0, 0)
Matthieu Constant's avatar
Matthieu Constant committed
149
150
151
152
153

tic = time.time()
#
# COMPUTATION LOOP
#
Nathan Coppin's avatar
Nathan Coppin committed
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186

xdata = []
ydata1 = []
ydata2 = []
ydata3 = []
ydata4 = []
fig, ((ax1,ax2),(ax3,ax4)) = plt.subplots(2,2)


ax1.plot(np.linspace(0,6,1000),np.zeros(1000))
ax1.set_xlim(0, 6)
ax1.set_ylim(-np.pi, np.pi)
ax1.set_title('Theta')
line1, = ax1.plot(xdata, ydata1, 'r-')

ax2.plot(np.linspace(0,6,1000),np.zeros(1000))
ax2.set_xlim(0, 6)
ax2.set_ylim(-30, 30)
ax2.set_title('Omega')
line2, = ax2.plot(xdata, ydata2, 'c-')

ax3.plot(np.linspace(0,6,1000),np.zeros(1000))
ax3.set_xlim(0, 6)
ax3.set_ylim(-0.1, 0.1)
ax3.set_title('Vx')
line3, = ax3.plot(xdata, ydata3, 'k-')

ax4.plot(np.linspace(0,6,1000),np.zeros(1000))
ax4.set_xlim(0, 6)
ax4.set_ylim(-0.2, 0.1)
ax4.set_title('Vy')
line4, = ax4.plot(xdata, ydata4, 'g-')

Matthieu Constant's avatar
Matthieu Constant committed
187
while t < tEnd :
Nathan Coppin's avatar
Nathan Coppin committed
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
    #xdata.append(t)
    #ydata1.append(p.body_theta()[4][0])
    #line1.set_xdata(xdata)
    #line1.set_ydata(ydata1)

    #ydata2.append(p.body_omega()[4][0])
    #line2.set_xdata(xdata)
    #line2.set_ydata(ydata2)

    #ydata3.append(p.body_velocity()[4][0])
    #line3.set_xdata(xdata)
    #line3.set_ydata(ydata3)

    #ydata4.append(p.body_velocity()[4][1])
    #line4.set_xdata(xdata)
    #line4.set_ydata(ydata4)

    #plt.draw()
    #plt.pause(1e-21)
    time_integration.iterate(fluid, p, dt, min_nsub=1, external_particles_forces=g*p.r()**2*rhop*np.pi)
Matthieu Constant's avatar
Matthieu Constant committed
208
209
210
211
212
    t += dt

    # Output files writting
    if ii %outf == 0 :
        ioutput = int(ii/outf) + 1
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
213
        p.write_vtk(outputdir, ioutput, t)
Nathan Coppin's avatar
Nathan Coppin committed
214
        fluid.write_vtk(outputdir, ioutput, t)
Matthieu Constant's avatar
Matthieu Constant committed
215
216
    ii += 1
    print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))
Nathan Coppin's avatar
Nathan Coppin committed
217
218
219
    print("\n\nCluster velocity : %f\n\n"%(np.linalg.norm(p.body_velocity()[4,:])))

plt.show()