drop.py 3.76 KB
 Matthieu Constant committed Jun 23, 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 30 31 32 33 34 35 36 37 38 39 40 41 42 43 ``````#!/usr/bin/env python from marblesbag import fluid as fluid from marblesbag import scontact2 import numpy as np import os import time import shutil import random def genInitialPosition(filename, r1, r2, rout, rhop1, rhop2) : p = scontact2.ParticleProblem() p.load_msh_boundaries("mesh.msh", ["Top", "Bottom","Left","Right"]) e = .7e-4 x = np.arange(rout, -rout, -e) y = np.arange(rout, -rout, -e) x, y = np.meshgrid(x, y) R2 = x**2 + y**2 keep = R2 < (rout - r2/2)**2 x = x[keep] y = y[keep] N1 = 0 N2 = 0 for i in range(x.shape[0]) : a = random.choice([1,2]) if a == 1 and N1 <= 128: N1 += 1 p.add_particle((x[i]+random.uniform(-e/3,e/3), y[i]+random.uniform(-e/3,e/3)), r1, r1**2 * np.pi * rhop1); elif a == 2 and N2 <= 128: N2 += 1 p.add_particle((x[i]+random.uniform(-e/3,e/3), y[i]+random.uniform(-e/3,e/3)), r2, r2**2 * np.pi * rhop2); else: if N1 > 128 and a == 1: p.add_particle((x[i]+random.uniform(-e/3,e/3), y[i]+random.uniform(-e/3,e/3)), r2, r2**2 * np.pi * rhop2); if N2 > 128 and a == 2: p.add_particle((x[i]+random.uniform(-e/3,e/3), y[i]+random.uniform(-e/3,e/3)), r1, r1**2 * np.pi * rhop1); print('N1',N1,'N2',N2,'N',N1+N2) p.position()[:, 1] += 0.45 p.write(filename) `````` Matthieu Constant committed Jul 10, 2017 44 ``````outputdir = "outputGoutteStokesGoodNew" `````` Matthieu Constant committed Jun 23, 2017 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 ``````if not os.path.isdir(outputdir) : os.makedirs(outputdir) t = 0 ii = 0 r1=25e-6/2 r2 = 25e-6/4 R = 1e-3 p = scontact2.ParticleProblem() #R = np.random.uniform(45e-06, 90e-06, len(x)) #physical parameters g = -9.81 rho = 1200 rhop1 = 2400 rhop2 = 2400 `````` Matthieu Constant committed Jul 10, 2017 62 ``````nu = 1e-4#/(100**0.5) `````` Matthieu Constant committed Jun 23, 2017 63 64 ``````V = 0.5 # todo : estimate V base on limit velocity print('V',V) `````` Matthieu Constant committed Jul 10, 2017 65 ``````tEnd = 10 `````` Matthieu Constant committed Jun 23, 2017 66 67 68 `````` #numerical parameters lcmin = 0.00005 # approx r*100 but should match the mesh size `````` Matthieu Constant committed Jul 10, 2017 69 ``````dt = 5e-4 `````` Matthieu Constant committed Jun 23, 2017 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 ``````alpha = 2.5e-4 epsilon = alpha*lcmin**2 /nu#2e-2*c*h*2#2e-2*c*h # ?? not sure ??1e-5 autoEpsilon = False print('epsilon',epsilon) shutil.copy("mesh.msh", outputdir +"/mesh.msh") #scontact2Interface.MeshLoader(p, "funnel.msh", ("Funnel", "Top", "Bottom", "Lateral")) p.write(outputdir+"/part-00000") genInitialPosition(outputdir + "/part-00000", r1, r2, R, rhop1, rhop2) p = scontact2.ParticleProblem(outputdir+"/part-00000") print("r = %g, m = %g\n" % (p.r()[0], p.mass()[0])) print("RHOP = %g" % rhop1) `````` Matthieu Constant committed Jul 10, 2017 85 ``````outf = 50 `````` Matthieu Constant committed Jun 23, 2017 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 ``````outf1 = 100000 strong_boundaries = [("Top",2,0.),("Top",1,0.),("Top",0,0.),("Bottom",1,0.),("Bottom",0,0.),("Left",0,0.),("Left",1,0.),("Right",0,0.),("Right",1,0.)] fluid = fluid.fluid_problem("mesh.msh",g,nu*rho,rho,epsilon,autoEpsilon,strong_boundaries) fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity()) fluid.export(outputdir,0,0) ii = 0 tic = time.clock() forces = g*p.mass() while t < tEnd : `````` Matthieu Constant committed Jul 10, 2017 101 102 103 `````` if (ii%100==0 and ii != 0) or ii==10: fluid.adapt_mesh(0.01,100,5e-5,autoEpsilon) forces = fluid.compute_node_force(dt,10*r1) `````` Matthieu Constant committed Jun 23, 2017 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 `````` 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) fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity()) fluid.implicit_euler(dt) t += dt if ii %outf == 0 : ioutput = int(ii/outf) + 1 p.write_vtk(outputdir, ioutput, t) p.write(outputdir+"/part-%05d" % ioutput) p.write_boundary_vtk(outputdir, ioutput, t) fluid.export(outputdir, t, ioutput) ii += 1 print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.clock() - tic))``````