drag.py 4.41 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
import numpy as np
import os
31
32
33
import sys
index = int(sys.argv[1])
vitesse = ["36","64","95","130","303","408"]
34
# Physical parameters
35
vit = -0.001*int(vitesse[index])               #stream velocity
36
rhop = 2500                                    #particles density
Nathan Coppin's avatar
Nathan Coppin committed
37
38
rho = 1000                                     #fluid density
nu = 1e-6                                      #fluid kinematic viscosity
Nathan Coppin's avatar
Nathan Coppin committed
39
g = np.array([0,-9.81])                        #gravity
40
41
42
# Numerical parameters
dt = 1e-3                                      #time step
t = 0                                          #initial time
43
tEnd = 10                                      #final time
44
ii = 0                                         #iteration number
Nathan Coppin's avatar
Nathan Coppin committed
45
outf = 10                                      #iterations between data frames  
46
# Define output directory
47
outputdir = "notOkTrial/" + vitesse[index] 
48
49
if not os.path.isdir(outputdir) :
    os.makedirs(outputdir)
50
51
52
53
54
filename = outputdir + "/Drag.csv"
# 
# PARTICLE PROBLEM
#
#Injected particles
55
p = scontact.ParticleProblem(2,True,True)
56
p.read_vtk("depot/outputNotOk",10)
57
p.remove_particles_flag((p.position()[:,1]-p.r()[:,0]>0.5)*(p.position()[:,1]+p.r()[:,0] <0.61))
58
#Flow particles
59
p2 = scontact.ParticleProblem(2,True,True)
60
p2.read_vtk("depot/outputNotOk",10)
61
p2.remove_particles_flag(p2.position()[:,1]+p2.r()[:,0] <0.5)
62
p2.clear_boundaries()
Nathan Coppin's avatar
Nathan Coppin committed
63
64
65
p2.load_msh_boundaries("mesh.msh", ["Inner","Left","Right"],material = "Steel") 
p2.set_friction_coefficient(0.3,"Sand","Sand")
p2.set_friction_coefficient(0.3,"Sand","Steel")
66
p2.write_vtk(outputdir, 0, 0)
67
68
69
#
# FLUID PROBLEM
#
70
71
72
73
74
75
76
#fluid = fluid.FluidProblem(2,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])
#fluid.set_wall_boundary("Right", velocity = [0,0])
#fluid.set_open_boundary("Bottom", velocity = [0,vit])
#fluid.set_open_boundary("Top",pressure = 0)
77
#fluid.set_particles(p2.mass(), p2.volume(), p2.position(), p2.velocity(), p2.contact_forces(),init=True)
78
#fluid.export_vtk(outputdir, 0.,0)
79
80
# 
# COMPUTATION LOOP
Nathan Coppin's avatar
Nathan Coppin committed
81
#
82
def accumulate(F,nsub):
83
  F+=np.sum(p2.get_boundary_forces("Inner"),axis=0)/nsub
84
while t < tEnd : 
85
     p2.forced_flag()[(p2.position()[:,1] + p2.r()[:,0] <= -0.5)] = 1
86
     if t > 2:
87
       p2.velocity()[p2.forced_flag()==1,:] = [0, vit]
88
     else :
89
90
91
92
93
       p2.velocity()[p2.forced_flag()==1,:] = [0, 0]
     p2.omega()[p2.forced_flag()==0] = 1
     if  max(p2.position()[:,1]+p2.r()[:,0]) < 0.5 and t < 9 :
       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.55))
94
     F = np.zeros(2)
95
96
     #time_integration.particle_changed(fluid,p2)
     time_integration.iterate(None,p2,dt,1,contact_tol=5e-7,external_particles_forces=g*p2.mass(),after_sub_iter=lambda nsub : accumulate(F,nsub))
97
     with open(filename,"a") as file1:
98
        #F += fluid.boundary_force()[0,:]
Nathan Coppin's avatar
Nathan Coppin committed
99
        file1.write( str(F[0])+";"+str(F[1])+";"+str(t)+"\n")
100
101
102
103
     t += dt   
     # Output files writing
     if ii %outf == 0 :
         ioutput = int(ii/outf) + 1
104
         p2.write_vtk(outputdir, ioutput, t)
105
         #fluid.export_vtk(outputdir+"/Flow", t, ioutput)       
106
     ii += 1
Nathan Coppin's avatar
Nathan Coppin committed
107
     print("%i : %.2g/%.2g" % (ii, t, tEnd))