drag.py 4.74 KB
Newer Older
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
1
# MigFlow - Copyright (C) <2010-2020>
Nathan Coppin's avatar
Nathan Coppin committed
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
# <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

# TESTCASE DESCRIPTION
Nathan Coppin's avatar
Nathan Coppin committed
25
# Study of the drag on a fixed cylinder in a 3D immersed granular flow
Nathan Coppin's avatar
Nathan Coppin committed
26
27
28
29
30
31
32
from migflow import scontact
from migflow import fluid
from migflow import time_integration
import numpy as np
import os
import time
import shutil
33
import random
Nathan Coppin's avatar
Nathan Coppin committed
34
import sys
Nathan Coppin's avatar
Nathan Coppin committed
35
listeVit = [14,36,64,95,130,169,211,255,303,354,408,467]
Nathan Coppin's avatar
Nathan Coppin committed
36
37
38
inter = int(sys.argv[1])
index = int(inter/5)
numDep = inter%5
Nathan Coppin's avatar
Nathan Coppin committed
39
# Physical parameters
Nathan Coppin's avatar
Nathan Coppin committed
40
vit = -0.001*listeVit[index]                   #stream velocity
Nathan Coppin's avatar
Nathan Coppin committed
41
muwall = 0.3                                   #friction coefficient between the walls and the particles
Nathan Coppin's avatar
Nathan Coppin committed
42
mupart = 0.15                                  #friction coefficient between particles
Nathan Coppin's avatar
Nathan Coppin committed
43
rho = 1000                                     #fluid density
Nathan Coppin's avatar
Nathan Coppin committed
44
nu = 1e-6                                      #fluid kinematic viscosity
45
g = np.array([0,-9.81,0])                      #gravity
Nathan Coppin's avatar
Nathan Coppin committed
46
47
48
# Numerical parameters
dt = 1e-3                                      #time step
t = 0                                          #initial time
Nathan Coppin's avatar
Nathan Coppin committed
49
tEnd = 30                                      #final time
Nathan Coppin's avatar
Nathan Coppin committed
50
ii = 0                                         #iteration number
Nathan Coppin's avatar
Nathan Coppin committed
51
outf = 2500                                    #iterations between data frames  
Nathan Coppin's avatar
Nathan Coppin committed
52
# Define output directory
Nathan Coppin's avatar
Nathan Coppin committed
53
outputdir = "output"
Nathan Coppin's avatar
Nathan Coppin committed
54
55
56
57
58
59
60
61
if not os.path.isdir(outputdir) :
    os.makedirs(outputdir)
filename = outputdir + "/Drag.csv"
# 
# PARTICLE PROBLEM
#
#Injected particles
p = scontact.ParticleProblem(3,True,True)
Nathan Coppin's avatar
Nathan Coppin committed
62
p.read_vtk("depot/"+ str(numDep),2)
63
64
p.remove_particles_flag((p.position()[:,1]-p.r()[:,0]>0.5)*(p.position()[:,1]+p.r()[:,0] <0.52))
p.velocity()[:,1] = vit
Nathan Coppin's avatar
Nathan Coppin committed
65
66
#Flow particles
p2 = scontact.ParticleProblem(3,True,True)
Nathan Coppin's avatar
Nathan Coppin committed
67
p.read_vtk("depot/"+ str(numDep),2)
68
p2.remove_particles_flag(p2.position()[:,1]+p2.r()[:,0] <0.5)
Nathan Coppin's avatar
Nathan Coppin committed
69
p2.clear_boundaries()
70
71
72
73
p2.load_msh_boundaries("mesh.msh", ["Inner","Left", "Right","Front","Rear"],material = "Steel")
p2.set_friction_coefficient(mupart,"Sand","Sand")
p2.set_friction_coefficient(muwall,"Sand","Steel")
p2.write_vtk(outputdir, 0, 0)
Nathan Coppin's avatar
Nathan Coppin committed
74
75
76
#
# FLUID PROBLEM
#
Nathan Coppin's avatar
Nathan Coppin committed
77
78
79
80
81
82
83
84
85
86
87
fluid = fluid.FluidProblem(3,g,[nu*rho],[rho],drag_in_stab=0)
fluid.load_msh("mesh.msh")
fluid.set_wall_boundary("Inner")
fluid.set_wall_boundary("Left", velocity = [0,0,0])
fluid.set_wall_boundary("Right", velocity = [0,0,0])
fluid.set_open_boundary("Bottom", velocity = [0,vit,0])
fluid.set_wall_boundary("Front", velocity = [0,0,0])
fluid.set_wall_boundary("Rear", velocity = [0,0,0])
fluid.set_open_boundary("Top",pressure = 0)
fluid.set_particles(p2.mass(), p2.volume(), p2.position(), p2.velocity(), p2.contact_forces())
fluid.export_vtk(outputdir, 0,0)
Nathan Coppin's avatar
Nathan Coppin committed
88
89
90
# 
# COMPUTATION LOOP
#
91
92
def accumulate(F,n_divide):
  F+=np.sum(p2.get_boundary_forces("Inner"),axis=0)/n_divide
Nathan Coppin's avatar
Nathan Coppin committed
93
while t < tEnd : 
94
95
96
97
98
99
     p2.forced_flag()[(p2.position()[:,1] + p2.r()[:,0] <= -0.5)] = 1
     p2.velocity()[p2.forced_flag()==1,:] = [0, vit, 0]
     p2.omega()[p2.forced_flag()==1,:] = [0,0,0]
     if  max(p2.position()[:,1]+p2.r()[:,0]) < 0.5 and t < 21 :
       p2.add_particles(p.position(),p.r(),p.mass(),v=p.velocity(),tag="Sand",forced=p.forced_flag(),contact_forces=p.contact_forces())
     p2.remove_particles_flag( (p2.position()[:,1] + p2.r()[:,0] >-0.52))
100
     F = np.zeros(3)
Nathan Coppin's avatar
Nathan Coppin committed
101
     fluid.set_particles(p2.mass(), p2.volume(), p2.position(), p2.velocity(), p2.contact_forces())
102
     time_integration.iterate(None,p2,dt,1,contact_tol=1e-6,external_particles_forces=g*p2.mass(),after_sub_iter=lambda n_divide : accumulate(F,n_divide))
Nathan Coppin's avatar
Nathan Coppin committed
103
     with open(filename,"a") as file1:
Nathan Coppin's avatar
Nathan Coppin committed
104
        F += fluid.boundary_force()[0,:]
Nathan Coppin's avatar
Nathan Coppin committed
105
        file1.write( str(F[0])+";"+str(F[1])+";"+str(F[2])+";"+str(t)+"\n")
Nathan Coppin's avatar
Nathan Coppin committed
106
107
108
109
     t += dt   
     # Output files writing
     if ii %outf == 0 :
         ioutput = int(ii/outf) + 1
110
         p2.write_vtk(outputdir, ioutput, t)
Nathan Coppin's avatar
Nathan Coppin committed
111
         fluid.write_vtk(outputdir, ioutput, t)       
Nathan Coppin's avatar
Nathan Coppin committed
112
113
     ii += 1
     print("%i : %.2g/%.2g" % (ii, t, tEnd))