Commit 3193e0f6 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

fix poiseuille-2fluids.py

parent b40fadaa
Pipeline #9736 passed with stages
in 3 minutes and 50 seconds
......@@ -21,17 +21,14 @@
#!/usr/bin/env python
from migflow import fluid as mbfluid
from migflow import scontact2
from migflow import time_integration
import numpy as np
import os
import subprocess
import time
import shutil
import random
import subprocess
#Physical parameters for the drops are the ones presented by Metzger et al. (2007) "Falling clouds of particles in viscous fluids"
import unittest
class Poiseuille(unittest.TestCase) :
......@@ -41,14 +38,14 @@ class Poiseuille(unittest.TestCase) :
outputdir = "output"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
subprocess.call(["gmsh", "-2", "mesh.geo","-clscale", "2"])
subprocess.call(["gmsh", "-2", "mesh.geo","-clscale", "2"])
t = 0
ii = 0
#physical parameters
g = 0 # gravity
g = np.array([0,0]) # gravity
rho = 1000 # fluid density
nu0 = 1e-3 # kinematic viscosity
nu1 = 1e-3 # kinematic viscosity
......@@ -58,62 +55,34 @@ class Poiseuille(unittest.TestCase) :
lcmin = .1 # mesh size
dt = 10 # time step
alpha = 1e-4 # stabilization coefficient
epsilon = alpha*lcmin**2 /min(nu0,nu1) # stabilization parametre
print('epsilon',epsilon)
shutil.copy("mesh.msh", outputdir +"/mesh.msh")
outf = 1 # number of iterations between output files
#Object fluid creation + Boundary condition of the fluid (field 0 is horizontal velocity; field 1 is vertical velocity; field 2 is pressure)
#Format: strong_boundaries = [(Boundary tag, Fluid field, Value)
u0,v0,q0 = 0,1,2
u1,v1,q1 = 3,4,5
p = 6
fluid = mbfluid.fluid_problem(g,[nu0*rho,nu1*rho],[rho,rho])
fluid.set_strong_boundary("Bottom",u0,0)
fluid.set_strong_boundary("Bottom",v0,0.)
fluid.set_strong_boundary("Bottom",u1,0)
fluid.set_strong_boundary("Bottom",v1,0.)
fluid.set_strong_boundary("Top",u0,0)
fluid.set_strong_boundary("Top",v0,0.)
fluid.set_strong_boundary("Top",u1,0)
fluid.set_strong_boundary("Top",v1,0.)
fluid.set_strong_boundary("LeftUp",u0,lambda x : x[:,1]/(20*nu0*rho)*x[:,1]*(1-x[:, 1]))
fluid.set_strong_boundary("LeftUp",v0,0)
fluid.set_strong_boundary("LeftUp",q0,lambda x : x[:,1])
fluid.set_strong_boundary("LeftUp",u1,lambda x : (1-x[:,1])/(20*nu1*rho)*x[:,1]*(1-x[:, 1]))
fluid.set_strong_boundary("LeftUp",v1,0)
fluid.set_strong_boundary("LeftDown",u0,lambda x : (x[:,1])/(20*nu0*rho)*x[:,1]*(1-x[:, 1]))
fluid.set_strong_boundary("LeftDown",v0,0)
fluid.set_strong_boundary("LeftDown",q0,lambda x : x[:,1])
fluid.set_strong_boundary("LeftDown",u1,lambda x : (1-x[:,1])/(20*nu1*rho)*x[:,1]*(1-x[:, 1]))
fluid.set_strong_boundary("LeftDown",v1,0)
fluid.set_strong_boundary("RightDown",v0,0)
fluid.set_strong_boundary("RightDown",v1,0)
fluid.set_strong_boundary("RightUp",v0,0)
fluid.set_strong_boundary("RightUp",v1,0)
fluid.load_mesh("mesh.msh")
fluid = mbfluid.FluidProblem(2,g,[nu0*rho,nu1*rho],[rho,rho])
fluid.load_msh("mesh.msh")
fluid.set_open_boundary("LeftUp",velocity=[lambda x : 1/(20*nu0*rho)*x[:,1]*(1-x[:, 1]),0],concentration=1)
fluid.set_open_boundary("LeftDown",velocity=[lambda x : 1/(20*nu1*rho)*x[:,1]*(1-x[:, 1]),0],concentration=0)
fluid.set_wall_boundary("Bottom",velocity=[0,0])
fluid.set_wall_boundary("Top",velocity=[0,0])
fluid.set_open_boundary("RightUp",pressure=0)
fluid.set_open_boundary("RightDown",pressure=0)
ii = 0
t = 0
#set initial_condition
x = fluid.coordinates()
s = fluid.solution()
s[:,2] = x[:,1]
s[:,5] = 1-s[:,2]
mu = nu0*rho
vI = 1/(120*mu)
fluid.write_vtk(outputdir,0,0)
ii = 0
tic = time.time()
while ii < 100:
while ii < 1000:
#Fluid solver
fluid.implicit_euler(dt)
time_integration.iterate(fluid,None,dt)
t += dt
#Output files writting
if ii %outf == 0 :
......@@ -124,9 +93,9 @@ class Poiseuille(unittest.TestCase) :
s = fluid.solution()
x = fluid.coordinates()
vel = (s[:,0]+s[:,3]-1/(20*nu0*rho)*x[:,1]*(1-x[:, 1]))**2
vel = (s[:,0]-1/(20*nu0*rho)*x[:,1]*(1-x[:, 1]))**2
vS = np.sum((1/(20*nu0*rho)*x[:,1]*(1-x[:, 1]))**2)
print('Error', (vel.sum())**.5)
self.assertLess(vel.sum()**.5,(vS**0.5)/100, "error is too large in Poiseuille")
self.assertLess(vel.sum()**.5,(vS**0.5)/50, "error is too large in Poiseuille")
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment