oneGrain.py 4.71 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
# Marblesbag - Copyright (C) <2010-2018>
# <Universite catholique de Louvain (UCL), Belgium
#  Universite de Montpellier, France>
# 	
# List of the contributors to the development of Marblesbag: see AUTHORS file.
# Description and complete License: see LICENSE file.
# 	
# This program (Marblesbag) 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
from marblesbag import fluid as mbfluid
from marblesbag import scontact2

import numpy as np
import os
import subprocess
import time
import shutil
import random
import unittest

class Poiseuille(unittest.TestCase) :
    def runTest(self) :
        dir_path = os.path.dirname(os.path.realpath(__file__))
        os.chdir(dir_path)


        outputdir = "output"
        if not os.path.isdir(outputdir) :
            os.makedirs(outputdir)
        subprocess.call(["gmsh", "-2", "mesh.geo","-clscale","1"])
        t = 0
        ii = 0

        R = [1 ,2, 3, 4, 6, 8, 10]
        R = np.array(R)*.5e-3
        V = []
        Vref = [0.15754716981132066,0.27735849056603756,0.3745283018867921,0.45188679245282976,0.5698113207547165,0.6490566037735843,0.7283018867924522]

        for r in R:
            print("Start with r=%g"%(r/2))
            #physical parameters
            g =  -9.81                                      # gravity
            rho = 1000                                      # fluid density
            rhop = 2640
            nu = 1e-6                                   # kinematic viscosity
            mu = nu*rho                                     # dynamic viscosity
            tEnd = 10                                     # final time

            #numerical parameters
            lcmin = .1                                  # mesh size
            dt = .01                                       # time step

            shutil.copy("mesh.msh", outputdir +"/mesh.msh")
            p = scontact2.ParticleProblem()
            p.load_msh_boundaries("mesh.msh", ["Top", "Bottom","Lateral"])
            p.add_particle((0,4.5),r,r**2*np.pi*rhop)
            p.write_vtk(outputdir,0,0)
            p = scontact2.ParticleProblem()
            p.read_vtk("output",0)
            p.write_vtk(outputdir, 0, 0)
            print("r = %g, m = %g\n" %  (p.r()[0], p.mass()[0]))
            print("RHOP = %g" % rhop)

            outf = 1                                       # number of iterations between output files

79
            fluid = mbfluid.fluid_problem(g,nu*rho,rho)
80
            fluid.load_msh("mesh.msh")
81
            fluid.set_strong_boundary("Top",2,0.)
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
            fluid.set_weak_boundary("Bottom","Wall")
            fluid.set_weak_boundary("Lateral","Wall")
            fluid.set_weak_boundary("Top","Wall")
            ii = 0
            t = 0


            fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())

            fluid.export_vtk(outputdir,0,0)

            ii = 0
            tic = time.clock()
            while ii < 100 : 
                #Fluid solver
                fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
                fluid.implicit_euler(dt)
                forces = fluid.compute_node_force(dt)
                vn = p.velocity() + forces * dt / p.mass()
                vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
                nsub = max(1, int(np.ceil((vmax * dt * 4)/min(10*p.r()))))
                print("NSUB", nsub,"VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(p.r()))
                for i in range(nsub) :
                    p.iterate(dt/nsub, forces) 
                t += dt
                #Output files writting
                if ii %outf == 0 :
                    ioutput = int(ii/outf) + 1
                    fluid.export_vtk(outputdir, t, ioutput)
                    p.write_vtk(outputdir, ioutput, t)
                ii += 1
                print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.clock() - tic))
            V.append(abs(p.velocity()[0,1]))

        V = np.array(V)
        Vref = np.array(Vref)
        R = 0
        for i in range(len(V)):
            R += abs(V[i]-Vref[i])/Vref[i]
        R /= len(V)

        self.assertLess(R,5./100., "error is too large in oneGrain")