drag.py 4.69 KB
Newer Older
Nathan Coppin's avatar
Nathan Coppin committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
# MigFlow - Copyright (C) <2010-2018>
# <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)
Nathan Coppin's avatar
Nathan Coppin committed
63
64
65
state = p.state()
p.remove_particles_flag((state.x[:,1]-p.r()[:,0]>0.5)*(state.x[:,1]+p.r()[:,0] <0.52))
state.v[:,1] = vit
Nathan Coppin's avatar
Nathan Coppin committed
66
67
#Flow particles
p2 = scontact.ParticleProblem(3,True,True)
Nathan Coppin's avatar
Nathan Coppin committed
68
p.read_vtk("depot/"+ str(numDep),2)
Nathan Coppin's avatar
Nathan Coppin committed
69
70
state2 = p2.state()
p2.remove_particles_flag(state2.x[:,1]+p2.r()[:,0] <0.5)
Nathan Coppin's avatar
Nathan Coppin committed
71
p2.clear_boundaries()
72
73
74
75
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
76
77
78
#
# FLUID PROBLEM
#
Nathan Coppin's avatar
Nathan Coppin committed
79
80
81
82
83
84
85
86
87
88
89
#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(),init=True)
#fluid.export_vtk(outputdir, 0.,0)
Nathan Coppin's avatar
Nathan Coppin committed
90
91
92
# 
# COMPUTATION LOOP
#
93
def accumulate(F):
Nathan Coppin's avatar
Nathan Coppin committed
94
  F+=np.sum(p2.get_boundary_impulses("Inner"),axis=0)
Nathan Coppin's avatar
Nathan Coppin committed
95
while t < tEnd : 
Nathan Coppin's avatar
Nathan Coppin committed
96
97
98
99
100
101
     p2.forced_flag()[(state2.x[:,1] + p2.r()[:,0] <= -0.5)] = 1
     state2.v[p2.forced_flag()==1,:] = [0, vit, 0]
     state2.omega[p2.forced_flag()==1,:] = [0,0,0]
     if  max(state2.x[:,1]+p2.r()[:,0]) < 0.5 and t < 21 :
       p2.add_particles(state.x,p.r(),p.mass(),v=state.v,tag="Sand",forced=p.forced_flag(),contact_forces=p.contact_forces())
     p2.remove_particles_flag( (state2.x[:,1] + p2.r()[:,0] >-0.52))
102
     F = np.zeros(3)
Nathan Coppin's avatar
Nathan Coppin committed
103
104
     #time_integration.particle_changed(fluid,p2)
     time_integration.iterate(None,p2,dt,1,contact_tol=1e-6,external_particles_forces=g*p2.mass(),after_sub_iter=lambda subdt : accumulate(F))
Nathan Coppin's avatar
Nathan Coppin committed
105
     with open(filename,"a") as file1:
Nathan Coppin's avatar
Nathan Coppin committed
106
        F /= dt
Nathan Coppin's avatar
Nathan Coppin committed
107
        #F += fluid.boundary_force()[0,:]
Nathan Coppin's avatar
Nathan Coppin committed
108
        file1.write( str(F[0])+";"+str(F[1])+";"+str(F[2])+";"+str(t)+"\n")
Nathan Coppin's avatar
Nathan Coppin committed
109
110
111
112
     t += dt   
     # Output files writing
     if ii %outf == 0 :
         ioutput = int(ii/outf) + 1
113
         p2.write_vtk(outputdir, ioutput, t)
Nathan Coppin's avatar
Nathan Coppin committed
114
         #fluid.export_vtk(outputdir, t, ioutput)       
Nathan Coppin's avatar
Nathan Coppin committed
115
116
     ii += 1
     print("%i : %.2g/%.2g" % (ii, t, tEnd))