drag.py 4.21 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
import numpy as np
import os
import time
import shutil
import random
34
# Physical parameters
Nathan Coppin's avatar
Nathan Coppin committed
35
vit = -0.1                                     #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 = 13                                      #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
Nathan Coppin's avatar
Nathan Coppin committed
47
outputdir = "output" 
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)
Nathan Coppin's avatar
Nathan Coppin committed
56
p.read_vtk("depot/output",10)
Nathan Coppin's avatar
Nathan Coppin committed
57
58
state = p.state()
p.remove_particles_flag((state.x[:,1]-p.r()[:,0]>0.5)*(state.x[:,1]+p.r()[:,0] <0.53))
59
#Flow particles
60
p2 = scontact.ParticleProblem(2,True,True)
Nathan Coppin's avatar
Nathan Coppin committed
61
p2.read_vtk("depot/output",10)
Nathan Coppin's avatar
Nathan Coppin committed
62
63
state2 = p2.state()
p2.remove_particles_flag(state2.x()[:,1]+p2.r()[:,0] <0.5)
64
p2.clear_boundaries()
Nathan Coppin's avatar
Nathan Coppin committed
65
66
67
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")
68
p2.write_vtk(outputdir, 0, 0)
69
70
71
#
# FLUID PROBLEM
#
Nathan Coppin's avatar
Nathan Coppin committed
72
73
74
75
76
77
78
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)
Nathan Coppin's avatar
Nathan Coppin committed
79
fluid.set_particles(p2.mass(), p2.volume(), state2, p2.contact_forces(),init=True)
80
fluid.export_vtk(outputdir, 0.,0)
81
82
# 
# COMPUTATION LOOP
Nathan Coppin's avatar
Nathan Coppin committed
83
#
84
def accumulate(F):
Nathan Coppin's avatar
Nathan Coppin committed
85
  F+=np.sum(p2.get_boundary_impulses("Inner"),axis=0)
86
while t < tEnd : 
Nathan Coppin's avatar
Nathan Coppin committed
87
88
89
90
91
92
     p2.forced_flag()[(state2.x[:,1] + p2.r()[:,0] <= -0.5)] = 1
     state2.v[p2.forced_flag()==1,:] = [0, vit]
     state2.omega[p2.forced_flag()==1] = 0
     if  max(state2.x[:,1]+p2.r()[:,0]) < 0.5 and t < 25 :
       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.55))
93
     F = np.zeros(2)
Nathan Coppin's avatar
Nathan Coppin committed
94
95
     time_integration.particle_changed(fluid,p2)
     time_integration.iterate(fluid,p2,dt,1,contact_tol=5e-7,external_particles_forces=g*p2.mass(),after_sub_iter=lambda subdt : accumulate(F))
96
     with open(filename,"a") as file1:
Nathan Coppin's avatar
Nathan Coppin committed
97
        F += fluid.boundary_force()[0,:]
Nathan Coppin's avatar
Nathan Coppin committed
98
        file1.write( str(F[0])+";"+str(F[1])+";"+str(t)+"\n")
99
100
101
102
     t += dt   
     # Output files writing
     if ii %outf == 0 :
         ioutput = int(ii/outf) + 1
103
         p2.write_vtk(outputdir, ioutput, t)
Nathan Coppin's avatar
Nathan Coppin committed
104
         fluid.export_vtk(outputdir+"/Flow", t, ioutput)       
105
     ii += 1
Nathan Coppin's avatar
Nathan Coppin committed
106
     print("%i : %.2g/%.2g" % (ii, t, tEnd))