drag.py 7.14 KB
Newer Older
1
# MigFlow - Copyright (C) <2010-2018>
2
3
4
# <Universite catholique de Louvain (UCL), Belgium
#  Universite de Montpellier, France>
# 	
5
# List of the contributors to the development of MigFlow: see AUTHORS file.
6
7
# Description and complete License: see LICENSE file.
# 	
8
# This program (MigFlow) is free software: 
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
# 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

24
25
26
27
28
# TESTCASE DESCRIPTION
# Study of the drag on a fixed disk in a 2D immersed granular flow
from migflow import scontact
from migflow import fluid
from migflow import time_integration
29
30
31
32
33
34
import numpy as np
import os
import time
import shutil
import random

35
# Physical parameters
Nathan Coppin's avatar
Nathan Coppin committed
36
vit = -0.48                                    #stream velocity
37
rhop = 2500                                    #particles density
Nathan Coppin's avatar
Nathan Coppin committed
38
39
rho = 500                                     #fluid density
nu = 1e-8                                      #fluid kinematic viscosity
Nathan Coppin's avatar
Nathan Coppin committed
40
g = np.array([0,-21.81])
41
42
43
# Numerical parameters
dt = 1e-3                                      #time step
t = 0                                          #initial time
Nathan Coppin's avatar
Nathan Coppin committed
44
tEnd = 13                                      #final time
45
ii = 0                                         #iteration number
Nathan Coppin's avatar
Nathan Coppin committed
46
outf = 1                                      #iterations between data frames  
47
48

# Define output directory
49
outputdir = "ellipseForced"
50
51
if not os.path.isdir(outputdir) :
    os.makedirs(outputdir)
52
53
54
55
56
57
58
59
60
61
filename = outputdir + "/Drag.csv"
if os.path.exists(filename):
  print("Remove old results? (y/n)\n")
  decision = input()
  if decision.lower() == 'y':
    os.remove(filename)

# 
# PARTICLE PROBLEM
#
62

63
#Injected particles
64
p = scontact.ParticleProblem(2,True,True)
Nathan Coppin's avatar
Nathan Coppin committed
65
66
67
68
69
70
71
72
p.read_vtk("depot/depotEllipse",450) 
p.remove_particles_flag( (p.position()[:,1]  < -0.2339))                   
p.position()[:,:] = p.position()[:,(1,0)]
p.position()[:,:] += [0.38005,0.647]
p.remove_particles_flag( (p.position()[:,1]  > 0.507))
p.velocity()[:,:] = [0,0] 
flipped = np.copy(p.position())
flipped[:,0] = -flipped[:,0]
73
p.contact_forces()[:,:] = 0
Nathan Coppin's avatar
Nathan Coppin committed
74
p.omega()[:] = 0
Nathan Coppin's avatar
Nathan Coppin committed
75
p.forced_flag()[:] = 1
Nathan Coppin's avatar
Nathan Coppin committed
76
p.mass()[:] *= 0.05
77
78
p.write_vtk(outputdir + "/In_particles", 0, 0)
#Injected particles
79
80


81
#Flow particles
82
p2 = scontact.ParticleProblem(2,True,True)
83
p2.read_vtk(outputdir + "/In_particles",0)
84
85
86
p2.clear_boundaries()
p2.load_msh_boundaries("mesh.msh", ["Inner"],material = "Steel") 
p2.load_msh_boundaries("mesh.msh", ["Left", "Right"],material = "Plexi")
87
#p2.remove_particles_flag( (p2.position()[:,1]  < 0.5-2*p2.r()[:,0]))
Nathan Coppin's avatar
Nathan Coppin committed
88
#p2.add_particles(p.position(),p.r(),p.mass(),v=p.velocity(),tag="Sand",forced=p.forced_flag(),contact_forces=p.contact_forces())
Nathan Coppin's avatar
Nathan Coppin committed
89
90
91
92
93
94
95
96
97
98
99
100
101
102
p2.add_particles(p.position()[:,:] - [0,0.291],p.r(),p.mass(),v=np.zeros_like(p.velocity()),tag="Sand",forced=np.zeros_like(p.forced_flag()),contact_forces=p.contact_forces())
p2.add_particles(flipped - [0,0.291*2],p.r(),p.mass(),v=np.zeros_like(p.velocity()),tag="Sand",forced=np.zeros_like(p.forced_flag()),contact_forces=p.contact_forces())
p2.add_particles(p.position()[:,:] - [0,0.291*3],p.r(),p.mass(),v=np.zeros_like(p.velocity()),tag="Sand",forced=np.zeros_like(p.forced_flag()),contact_forces=p.contact_forces())
p2.add_particles(flipped - [0,0.291*4],p.r(),p.mass(),v=np.zeros_like(p.velocity()),tag="Sand",forced=np.zeros_like(p.forced_flag()),contact_forces=p.contact_forces())
p2.remove_particles_flag( (p2.position()[:,1]**2 + p2.position()[:,0]**2  > (0.01905 + max(p2.r()))**2))
p2.r()[:] *= 0.99


p2.set_friction_coefficient(0.28,"Sand","Sand")
p2.set_friction_coefficient(0.28,"Sand","Steel")
p2.set_friction_coefficient(0.14,"Sand","Plexi")
p2.set_tangent_boundary_velocity("Left",0)#vit)
p2.set_tangent_boundary_velocity("Right",0)#-vit)
p2.forced_flag()[p2.position()[:,1] + 2*p2.r()[:,0] >= 0.5] = 1
103
p2.contact_forces()[:,:] = 0
Nathan Coppin's avatar
Nathan Coppin committed
104
p2.velocity()[:,0] = 0
Nathan Coppin's avatar
switch    
Nathan Coppin committed
105
p2.velocity()[p2.position()[:,1] < 0,1] = 0
Nathan Coppin's avatar
Nathan Coppin committed
106
p2.omega()[:] = 0
Nathan Coppin's avatar
Nathan Coppin committed
107
p2.mass()[:,0] *= 0.05
108
109
p2.set_get_momentum(1) 
p2.write_vtk(outputdir+"/Flow", 0, 0)
110

111
112
113
#
# FLUID PROBLEM
#
114

115
#fluid = fluid.FluidProblem(2,0,[nu*rho],[rho])
116
#Set the mesh geometry for the fluid computation
Nathan Coppin's avatar
Nathan Coppin committed
117
#fluid.load_msh("mesh.msh")
118
119
120
121
122
123
124
#fluid.set_wall_boundary("Inner")
#fluid.set_wall_boundary("Left", velocity = [0,vit])
#fluid.set_wall_boundary("Right", velocity = [0,vit])
#fluid.set_open_boundary("Bottom", velocity = [0,vit])
#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+"/Flow", 0.,0)
125
126
127

# 
# COMPUTATION LOOP
128
#++	
Nathan Coppin's avatar
Nathan Coppin committed
129
print(p2.n_particles())
130
tic = time.time()
Nathan Coppin's avatar
Nathan Coppin committed
131
vitC = -0.01
132
while t < tEnd : 
Nathan Coppin's avatar
Nathan Coppin committed
133
134
135
136
137
138
139
140
141
142
143
     momentum = np.zeros_like(p2.momentum())
     if vitC >= vit and t>=500*dt:
       vitC += 2*dt*vit
     vitB = vitC#*0.72 if int((t+dt)/dt)%2==0 else vitC#-0.002 if int((t/dt))%3==0 else vitC+0.002
     vitC2 = vitC#*72.0 if int(t/dt)%2==0 else vitC
     print("vitC", vitC, "vitB", vitB)
     p2.forced_flag()[p2.position()[:,1] + 2*p2.r()[:,0] < 0.5] = 0
     p2.forced_flag()[(p2.position()[:,1] + 2*p2.r()[:,0] <= -0.5)+(p2.position()[:,1] + 2*p2.r()[:,0] >= 0.5)] = 1
     p2.velocity()[p2.forced_flag()==1,:] = [0, vitC2]
     p2.velocity()[(p2.position()[:,1] + 2*p2.r()[:,0] <= -0.5),1] = vitB
     p2.omega()[p2.forced_flag()==1] = 0
144
     if  max(p2.position()[:,1]) < 0.5 :
Nathan Coppin's avatar
Nathan Coppin committed
145
       p2.add_particles(p.position() if random.choice([True,False]) else flipped,p.r(),p.mass(),v=p.velocity(),tag="Sand",forced=p.forced_flag(),contact_forces=p.contact_forces())
Nathan Coppin's avatar
Nathan Coppin committed
146
     p2.remove_particles_flag( (p2.position()[:,1] + p2.r()[:,0] >-0.551))
147
     #time_integration.particle_changed(fluid,p2)
Nathan Coppin's avatar
Nathan Coppin committed
148
149
150
     G = np.zeros_like(p2.velocity())
     G[:,0] = p2.mass()[:,0]*g[1]
     time_integration.iterate(None,p2,dt,20,contact_tol=3e-8,momentum=momentum)
Nathan Coppin's avatar
Nathan Coppin committed
151
     #momentum += fluid.boundary_force()[0:3,:]
152
     print(momentum)
Nathan Coppin's avatar
Nathan Coppin committed
153
154
     speedAbove = np.mean(p.velocity()[(p.position()[:,1] < 0.45)*(p.position()[:,1] > 0.3),1])
     speedBelow = np.mean(p.velocity()[(p.position()[:,1] > -0.45)*(p.position()[:,1] < -0.3),1])
155
     with open(filename,"a") as file1:
Nathan Coppin's avatar
Nathan Coppin committed
156
157
158
        file1.write(str(momentum[0,0]) + ";" + str(momentum[0,1]) + ";" + str(momentum[1,0]) + ";" 
                    + str(momentum[1,1]) + ";" + str(momentum[2,0]) + ";" + str(momentum[2,1]) + ";" 
                    + str(t) + ";" + str(0 if np.isnan(speedAbove) else speedAbove) + ";" + str(0 if np.isnan(speedBelow) else speedBelow) + "\n")
159
160
161
162
163
     t += dt   
     # Output files writing
     if ii %outf == 0 :
         ioutput = int(ii/outf) + 1
         p2.write_vtk(outputdir+"/Flow", ioutput, t)
164
         #fluid.export_vtk(outputdir+"/Flow", t, ioutput)       
165
166
     ii += 1
     print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))
167