hourglass 2d

pyfefp/Brinkman.py

@@ -82,17 +82,18 @@ class Brinkman :
         F = solgrad[..., _p_, :]
         #drag
         GoU = f * Cd_u * self.rhof / 2 * Ap
-        reaction = (F * self.pf.volume + (GoU*U) / (1 + self.stab/self.pf.mass * GoU))*0.3
+        reaction = (F * self.pf.volume + (GoU*U) / (1 + self.stab/self.pf.mass * GoU))*self.dragFactor
         if fonpart is not None:
             fonpart[...] = -reaction - self.g * self.pf.volume * self.rhof
         val[..., _U_] =  reaction / self.rhof 
         val[..., _p_] = 0
         return val
 
-    def __init__(self, mu, pf, g, rho, rhop, epsilon, ns = False, stab = 0, D = 2, c=1) :
+    def __init__(self, mu, pf, g, rho, rhop, epsilon, ns = False, stab = 0, D = 2, c=1, dragFactor=0.3) :
         self.D = D
         self.c = c
         self.ns = ns
+        self.dragFactor = dragFactor
         self.rhof = rho
         self.rhop = rhop
         self.stab = stab
@@ -126,7 +127,6 @@ class Brinkman :
                 jac.add_to_field(f, g, pp[:, None, :] * ff[:, :, None], self.pf.eid)
 
     def fluidTerm(self, meshJ, solution, res, jac) :
-        print("FLUIDTERM")
         dofm = solution.dof_manager()
         qp, qw = quadrature.get_triangle(3) if self.D == 2 else quadrature.get_tetrahedron(3)
         self.vs = self.pf.momentum.at_qp(qp)

pyfefp/fluidSolver.py

@@ -9,7 +9,7 @@ from .legendre import *
 import time
 
 class FluidSolver :
-    def __init__(self, mesh, dt, rhop, rho, mu, c, g, epsilon, TypeEMu, TypeEMp, timeScheme="Implicit",imex=True) :
+    def __init__(self, mesh, dt, rhop, rho, mu, c, g, epsilon, TypeEMu, TypeEMp, dragFactor=0.3,timeScheme="Implicit", imex=True) :
         if isinstance(mesh, str):
             mesh = Mesh(mesh)
         self._mesh = mesh
@@ -28,6 +28,7 @@ class FluidSolver :
             rho = rho,
             rhop = rhop,
             epsilon = epsilon,
+            dragFactor=dragFactor,
             D = self._jac.dim,
             stab = dt)
         d = DofManager(mesh)
@@ -173,7 +174,6 @@ class FluidSolver :
     def solution(self) :
         return self._sol
 
-
     def set_mesh_position(self, x) :
         self._jac.update(x)
 

testcases/sablier-2d/depot-fluid.py

@@ -9,12 +9,11 @@ import os
 import time
 import random
 
-def genInitialPosition(filename, r, ox, oy, nx, ny, rhop) :
+def genInitialPosition(filename, r, ox, oy, lx, ly, rhop) :
     p = scontact2.ParticleProblem()
     scontact2Interface.MeshLoader(p, "mesh.msh", ["Top", "Box"])
-    x = np.arange(ox+2*r, ox+(2*nx*r), 2*r)
-    y = np.arange(oy+2*r, oy+2*ny*r, 2*r)
-
+    x = np.arange(ox+r, ox+lx, 2*r)
+    y = np.arange(oy+r, oy+ly, 2*r)
     for i in range(x.shape[0]) :
        for j in range(y.shape[0]) :
           rr=r*(0.8+random.random()*0.2)
@@ -22,47 +21,46 @@ def genInitialPosition(filename, r, ox, oy, nx, ny, rhop) :
     p.write(filename)
 
 
-outputdir = "output8"
+outputdir = "output-initpressure"
 if not os.path.isdir(outputdir) :
     os.makedirs(outputdir)
 
 t = 0
 ii = 0
 
-d=1e-3 # la valeur de ref dans le .geo
-r=0.1*d
+r=1e-4*3
 p = scontact2.ParticleProblem()
 
 #physical parameters
 g =  [0, -9.81]
 rho = 1200
 rhop = 1500 #2400
-nu = 1e-4
-V = 1e-3  #0.5 # todo : estimate V base on limit velocity
-print('V',V)
+nu = 1e-6
 tEnd = 20
 
 #numerical parameters
-h = 2e-4 # approx r*100 but should match the mesh size
-
-dt = 5e-4 #0.5*r/V
-
-c =  1e-2 #h/dt* 7./50.
+#h = 2e-4 # approx r*100 but should match the mesh size
 
+c =  20 #h/dt* 7./50.
+dt = 1e-4 #0.5*r/V
 epsilon = 1e-5 #2e-2*c*h # ?? not sure ??
 
 print("c = %g, eps = %g" % (c, epsilon))
 
-genInitialPosition(outputdir + "/part-00000", r, 0., 51*d, 75, 200, rhop)
-
+y0part = 0.051
+lypart = 0.025
+lxpart = 0.015
+genInitialPosition(outputdir + "/part-00000", r, 0., y0part, lxpart, lypart, rhop)
 p = scontact2Interface.ScontactInterface(outputdir+"/part-00000")
 
+deltamass = p.volume().sum()*(rhop-rho)
+deltapinit = deltamass*-g[1]/lxpart
+
 print("r = %g, m = %g\n" %  (p.r()[0], p.mass()[0]))
 print("RHOP = %g" % rhop)
-outf = 10 #000
-outf1 = 10000
+outf = 100 #000
 
-fluid = pyfefp.FluidSolver("mesh.msh", dt, rhop=rhop, rho=rho, g = g, mu=nu*rho, c=c*c, epsilon = epsilon,timeScheme="Explicit", TypeEMu=TriangleP1, TypeEMp=TriangleP1)
+fluid = pyfefp.FluidSolver("mesh.msh", dt, rhop=rhop, rho=rho, g = g, mu=nu*rho, c=c*c, epsilon = epsilon,timeScheme="Explicit", dragFactor=0.3, TypeEMu=TriangleP1, TypeEMp=TriangleP1)
 
 # impermeability
 fluid.add_boundary_condition("Top", 1, 0.)
@@ -71,47 +69,28 @@ fluid.add_boundary_condition("Box", 1, 0.)
 fluid.add_boundary_condition("Box", 0, 0.)
 
 fluid.complete()
-
+y = fluid.mesh_position()[...,1].transpose()
+pinit = (y0part+lypart-y)*deltapinit/lypart
+pinit[pinit<0] = 0
+pinit[pinit>deltapinit] = deltapinit
+fluid._dof.setField(2, pinit, fluid.solution())
 fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
 
-filenergy =open(outputdir+'/filenergy.txt','w')
-filenergy.write('Energie cinetique solide, Energie potentielle solide, Energie cinetique fluide, Energie totale\n')
-
-p.write_vtk(outputdir, 1, 0.)
-p.write_boundary_vtk(outputdir, 1, 0.)
-
 tic = time.clock()
 while t < tEnd :
     forces = fluid.solve() + g * p.mass()
-    vn = p.velocity() + forces * dt / p.mass()
-
-    vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
-    Ecs = sum(np.hypot(vn[:, 0], vn[:, 1])**2 *0.5)
-    Eps = sum(rhop*p.position()[:,1]*9.81)
-    if ii % outf1 == 0:
-        filenergy.write('%.12f, %.12f, %.12f, %.12f\n'% (Ecs,Eps,fluid._law.Ecf,Ecs+Eps+fluid._law.Ecf))
-    if (vmax > 10*V) :
-        print("CRASH !!")
-        #exit()
-    if ii==10000:
-        fluid._law.epsilon=2e-1*c*h
-    #print("NSUB", nsub,
-    print("VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(p.r()))        
-    nsub = max(1, int(np.ceil((vmax * dt * 4)/min(p.r()))))
-
-    for i in range(nsub) :
-        p.iterate(dt/nsub, forces)
+    NS = 10
+    if ii%NS == 0 :
+        p.iterate(dt*NS, forces)
     fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
     t += dt
-    print("dt",dt)
     if ii %outf == 0 :
         ioutput = int(ii/outf) + 1
-        p.write_vtk(outputdir, ioutput, t)
+        #p.write_vtk(outputdir, ioutput, t)
         p.write(outputdir, ioutput, t)
-        p.write_boundary_vtk(outputdir, ioutput, t)
-        #fluid.write_solution(outputdir, ioutput, t, "msh")
-        fluid.write_solution(outputdir, ioutput, t, "vtk")
+        #p.write_boundary_vtk(outputdir, ioutput, t)
+        fluid.write_solution(outputdir, ioutput, t, "msh")
+        #fluid.write_solution(outputdir, ioutput, t, "vtk")
     ii += 1
     print("%i : %.2g/%.2g (cpu %.2g)" % (ii, t, tEnd, time.clock() - tic))
 
-filenergy.close()

testcases/sablier-2d/mesh.geo

 d=1e-3;
 L = 15*d;
 H = 50*d;
-lc = 2*d;
-R=3*d;
+lc = 3*d;
+R=2*d;
 
 Point(1) = {0., 0., 0, lc};
 Point(2) = {L, 0., 0, lc};
 Point(3) = {L, H, 0, lc};
-Point(4) = {0.5*(L+R), H, 0, lc};
-Point(5) = {0.5*(L+R), H+d, 0, lc};
+Point(4) = {0.5*(L+R), H, 0, lc/3};
+Point(5) = {0.5*(L+R), H+d, 0, lc/3};
 Point(6) = {L, H+d, 0, lc};
 Point(7) = {L, H+d+H, 0, lc};
 Point(8) = {0., H+d+H, 0, lc};
 Point(9) = {0., H+d, 0, lc};
-Point(10) = {0.5*(L-R), H+d, 0, lc};
-Point(11) = {0.5*(L-R), H, 0, lc};
+Point(10) = {0.5*(L-R), H+d, 0, lc/3};
+Point(11) = {0.5*(L-R), H, 0, lc/3};
 Point(12) = {0., H, 0, lc};
 Line(1) = {1, 2};
 Line(2) = {2, 3};