drop.py 4.37 KB
 Matthieu Constant committed Nov 22, 2017 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 ``````#!/usr/bin/env python from marblesbag import fluid3 as fluid from marblesbag import scontact3 import numpy as np import os import time import shutil import random def genInitialPosition(filename, r, rout, rhop, compacity) : p = scontact3.ParticleProblem() p.load_msh_boundaries("mesh.msh", ["Top", "Bottom","X","Z"]) N = compacity*rout**3/(r**3) e = 2*rout/(6*N/np.pi)**(1./3.) for x in np.arange(rout, -rout, -e): for y in np.arange(rout, -rout, -e): for z in np.arange(rout, -rout, -e): R2 = x**2 + y**2 + z**2 if R2<(rout-r)**2: p.add_particle((x+random.uniform(-e/2+r,e/2-r), y+random.uniform(-e/2+r,e/2-r), z+random.uniform(-e/2+r,e/2-r)), r, 4./3.* r**3 * np.pi * rhop); p.position()[:, 1] += 0.52 print(len(p.position()[:, 1] )) p.write(filename) `````` Matthieu Constant committed Nov 30, 2017 30 ``````outputdir = "MetzgerB" `````` Matthieu Constant committed Nov 22, 2017 31 32 33 34 35 36 37 38 ``````if not os.path.isdir(outputdir) : os.makedirs(outputdir) t = 0 ii = 0 r=154e-6 R = 3.3e-3 `````` Matthieu Constant committed Nov 30, 2017 39 ``````compacity = .20 `````` Matthieu Constant committed Nov 22, 2017 40 41 42 43 44 45 46 47 48 49 50 ``````drho = 35.4 p = scontact3.ParticleProblem() #R = np.random.uniform(45e-06, 90e-06, len(x)) #physical parameters g = -9.81 rhop = 2450 nu = 1.17/1030. rho = 1030#rhop-drho/compacity V = 0.5 # todo : estimate V base on limit velocity print('V',V) `````` Matthieu Constant committed Nov 30, 2017 51 ``````tEnd = 1000 `````` Matthieu Constant committed Nov 22, 2017 52 53 54 `````` #numerical parameters lcmin = 0.0004 # approx r*100 but should match the mesh size `````` Jonathan Lambrechts committed Nov 22, 2017 55 ``````dt = 5e-5*1000 `````` Matthieu Constant committed Nov 22, 2017 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 ``````alpha = 1e-3 epsilon = alpha*lcmin**2 /nu print('epsilon',epsilon) shutil.copy("mesh.msh", outputdir +"/mesh.msh") #scontact2Interface.MeshLoader(p, "funnel.msh", ("Funnel", "Top", "Bottom", "Lateral")) p.write(outputdir+"/part-00000") mu = nu*rho genInitialPosition(outputdir + "/part-00000", r, R, rhop, compacity) p = scontact3.ParticleProblem(outputdir+"/part-00000") print("r = %g, m = %g\n" % (p.r()[0], p.mass()[0])) print("RHOP = %g" % rhop) `````` Matthieu Constant committed Nov 30, 2017 71 ``````outf = 10 `````` Matthieu Constant committed Nov 22, 2017 72 73 74 ``````outf1 = 100000 #strong_boundaries = [("Top",0,.0),("Top",1,0.),("Top",2,0.),("Box",1,0.),("Box",0,0.),("Box",2,0.),("Top",3,0.)] `````` Matthieu Constant committed Nov 30, 2017 75 76 77 ``````#strong_boundaries = [("Top",0,.0),("Top",1,0.),("Top",2,0.),("Bottom",0,.0),("Bottom",1,0.),("Bottom",2,0.),("X",0,.0),("X",1,0.),("X",2,0.),("Z",0,.0),("Z",1,0.),("Z",2,0.),("Top",3,0.)] strong_boundaries = [("Top",1,0.),("Bottom",1,0.),("X",0,.0),("Z",2,0.),("Top",3,0.)] fluid = fluid.fluid_problem("mesh.msh",g,nu*rho,rho,epsilon,strong_boundaries, 1) `````` Matthieu Constant committed Nov 22, 2017 78 79 `````` fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity()) `````` Matthieu Constant committed Nov 30, 2017 80 `````` `````` Matthieu Constant committed Nov 22, 2017 81 82 83 84 85 86 87 88 89 `````` ii = 0 jj = 0 tEnd1 = 100 #p.velocity()[:,1] = 2./9. *(compacity)* R**2*g*(rhop-rho)/mu print(p.velocity()[1,1]) print(1-compacity) `````` Matthieu Constant committed Nov 30, 2017 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 ``````for jj in range(4): dt1 = dt t = 0 if jj!=0: fluid.adapt_mesh(1e-3,1e-4,300000, epsilon, 1) fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity()) while t < tEnd1 : niter = fluid.implicit_euler(dt1) forces = fluid.compute_node_force(dt1) vn = p.velocity() + forces * dt1 / p.mass() fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity()) if niter < 5 and dt1 < 20: dt1 *= 3 elif niter > 5 and dt1 > dt: dt1 /= 3 print(dt1) t += dt1 `````` Matthieu Constant committed Nov 22, 2017 108 ``````# if ii %outf == 0 : `````` Matthieu Constant committed Nov 30, 2017 109 110 111 112 ``````# ioutput = int(ii/outf) + 1 # p.write_vtk(outputdir, ioutput, t) # fluid.export_vtk(outputdir, t, ioutput) ii += 1 `````` Matthieu Constant committed Nov 22, 2017 113 114 115 116 `````` ii = 0 t = 0 tic = time.clock() `````` Matthieu Constant committed Nov 30, 2017 117 ``````fluid.export_vtk(outputdir,0,0) `````` Matthieu Constant committed Nov 22, 2017 118 ``````forces = g*p.mass() `````` Jonathan Lambrechts committed Nov 22, 2017 119 ``````fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity()) `````` Matthieu Constant committed Nov 22, 2017 120 121 ``````while t < tEnd : if (ii%5==0 and ii != 0): `````` Matthieu Constant committed Nov 30, 2017 122 `````` fluid.adapt_mesh(1e-3,1e-4,300000, epsilon, 1) `````` Jonathan Lambrechts committed Nov 22, 2017 123 124 `````` fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity()) fluid.implicit_euler(dt) `````` Matthieu Constant committed Nov 22, 2017 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 `````` 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(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) fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity()) t += dt if ii %outf == 0 : ioutput = int(ii/outf) + 1 p.write_vtk(outputdir, ioutput, t) fluid.export_vtk(outputdir, t, ioutput) ii += 1 print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.clock() - tic))``````