Commit bf6c4877 authored by Matthieu Constant's avatar Matthieu Constant
Browse files

poiseuille avec 2 viscosités diff

parent 1b280b55
......@@ -48,28 +48,34 @@ tEnd = 100000 # final time
#numerical parameters
lcmin = .1 # mesh size
dt = 0.5 # time step
dt = 10 # time step
alpha = 1e-4 # stabilization coefficient
epsilon = alpha*lcmin**2 /nu # stabilization parametre
print('epsilon',epsilon)
shutil.copy("mesh.msh", outputdir +"/mesh.msh")
outf = 10 # number of iterations between output files
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)
strong_boundaries = [("Top",1,1,0.),("Bottom",1,1,0.),("Left",1,1,0.),("Right",1,1,0.),("Top",0,0,0.),("Bottom",0,0,0.),("Left",0,0,lambda x : 1/(20*mu)*x[:,1]*(1-x[:, 1]))]
fluid = fluid.fluid_problem("mesh.msh",g,[nu*rho],[rho],epsilon,strong_boundaries,1,1)
strong_boundaries = [("Top",1,1,0.),("Bottom",1,1,0.),("LeftUp",1,1,0.),("RightUp",1,1,0.),("LeftDown",1,1,0.),("RightDown",1,1,0.),("Top",0,0,0.),("Bottom",0,0,0.),("LeftUp",0,0,lambda x : 1/(20*mu)*x[:,1]*(1-x[:, 1])),("LeftDown",0,0,lambda x : 1/(20*mu)*x[:,1]*(1-x[:, 1]))]
fluid = fluid.fluid_problem("mesh.msh",g,[nu*rho],[rho],epsilon,strong_boundaries,1)
ii = 0
t = 0
#set initial_condition
x = fluid.coordinates()
s = fluid.solution()
a = (x[:,0]+5)/10
s[:,2] = 1.
fluid.export_vtk(outputdir,0,0)
ii = 0
tic = time.clock()
while t < tEnd :
while ii < 100 :
#Fluid solver
fluid.implicit_euler(dt)
t += dt
......@@ -80,3 +86,12 @@ while t < tEnd :
ii += 1
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.clock() - tic))
s = fluid.solution()
x = fluid.coordinates()
vel = (s[:,0]-1/(20*nu*rho)*x[:,1]*(1-x[:, 1]))**2
print('Error', (vel.sum())**.5)
if (vel.sum())**.5<1.5e-3:
exit(1)
else:
exit(0)
\ No newline at end of file
......@@ -41,14 +41,15 @@ ii = 0
#physical parameters
g = -9.81 # gravity
rho = 1000 # fluid density
nu0 = 1e-3 # kinematic viscosity
rho0 = 1000 # fluid density
rho1 = 1000
nu0 = 1e-2 # kinematic viscosity
nu1 = 1e-3 # kinematic viscosity
tEnd = 100000 # final time
#numerical parameters
lcmin = .1 # mesh size
dt = 10 # time step
dt = 1 # time step
alpha = 1e-4 # stabilization coefficient
epsilon = alpha*lcmin**2 /min(nu0,nu1) # stabilization parametre
print('epsilon',epsilon)
......@@ -71,15 +72,15 @@ strong_boundaries = [
("Top",v0,v0,0.),
("Top",u1,u1,0),
("Top",v1,v1,0.),
("LeftUp",u0,u0,lambda x : 1./(20*nu0*rho)*x[:,1]*(1-x[:, 1])),
("LeftUp",u0,u0,lambda x : (6*x[:,1]**5-15*x[:,1]**4+10*x[:,1]**3)/(20*nu0*rho0)*x[:,1]*(1-x[:, 1])),
("LeftUp",v0,v0,0),
("LeftUp",q0,q0,1),
("LeftUp",u1,u1,lambda x : 0/(20*nu1*rho)*x[:,1]*(1-x[:, 1])),
("LeftUp",q0,q0,lambda x : 6*x[:,1]**5-15*x[:,1]**4+10*x[:,1]**3),
("LeftUp",u1,u1,lambda x : (1-(6*x[:,1]**5-15*x[:,1]**4+10*x[:,1]**3))/(20*nu1*rho1)*x[:,1]*(1-x[:, 1])),
("LeftUp",v1,v1,0),
("LeftDown",u0,u0,lambda x : 1./(20*nu0*rho)*x[:,1]*(1-x[:, 1])),
("LeftDown",u0,u0,lambda x : (6*x[:,1]**5-15*x[:,1]**4+10*x[:,1]**3)/(20*nu0*rho0)*x[:,1]*(1-x[:, 1])),
("LeftDown",v0,v0,0),
("LeftDown",q0,q0,1),
("LeftDown",u1,u1,lambda x : 0/(20*nu1*rho)*x[:,1]*(1-x[:, 1])),
("LeftDown",q0,q0,lambda x : 6*x[:,1]**5-15*x[:,1]**4+10*x[:,1]**3),
("LeftDown",u1,u1,lambda x : (1-(6*x[:,1]**5-15*x[:,1]**4+10*x[:,1]**3))/(20*nu1*rho1)*x[:,1]*(1-x[:, 1])),
("LeftDown",v1,v1,0),
("RightDown",v0,v0,0),
("RightDown",v1,v1,0),
......@@ -87,7 +88,7 @@ strong_boundaries = [
("RightUp",v1,v1,0),
]
fluid = fluid.fluid_problem("mesh.msh",g,[nu0*rho,nu1*rho],[rho,rho],epsilon,strong_boundaries,2)
fluid = fluid.fluid_problem("mesh.msh",g,[nu0*rho0,nu1*rho1],[rho0,rho1],epsilon,strong_boundaries,2)
ii = 0
t = 0
......@@ -95,9 +96,12 @@ t = 0
x = fluid.coordinates()
s = fluid.solution()
a = (x[:,0]+5)/10
s[:,2] = 1-a*0.8
s[:,2] = 6*x[:,1]**5-15*x[:,1]**4+10*x[:,1]**3
s[:,5] = 1-s[:,2]
fluid.export_vtk(outputdir,0,0)
ii = 0
......@@ -106,6 +110,10 @@ while t < tEnd :
#Fluid solver
fluid.implicit_euler(dt)
t += dt
s = fluid.solution()
x = fluid.coordinates()
vel = (s[:,0]+s[:,3]-1/(20*nu0*rho)*x[:,1]*(1-x[:, 1]))**2
print('Error', (vel.sum())**.5)
#Output files writting
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
......
# 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 fluid
from marblesbag import scontact2
import numpy as np
import os
import time
import shutil
import random
#Physical parameters for the drops are the ones presented by Metzger et al. (2007) "Falling clouds of particles in viscous fluids"
outputdir = "output"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
t = 0
ii = 0
#physical parameters
g = -9.81 # gravity
rho0 = 1000 # fluid density
rho1 = 1000
nu0 = 1e-2 # kinematic viscosity
nu1 = 1e-3 # kinematic viscosity
tEnd = 100000 # final time
#numerical parameters
lcmin = .1 # mesh size
dt = 1 # 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
strong_boundaries = [
("Bottom",u0,u0,0),
("Bottom",v0,v0,0.),
("Bottom",u1,u1,0),
("Bottom",v1,v1,0.),
("Top",u0,u0,0),
("Top",v0,v0,0.),
("Top",u1,u1,0),
("Top",v1,v1,0.),
("LeftUp",u0,u0,lambda x : (6*x[:,1]**5-15*x[:,1]**4+10*x[:,1]**3)/(20*nu0*rho0)*x[:,1]*(1-x[:, 1])),
("LeftUp",v0,v0,0),
("LeftUp",q0,q0,lambda x : 6*x[:,1]**5-15*x[:,1]**4+10*x[:,1]**3),
("LeftUp",u1,u1,lambda x : (1-(6*x[:,1]**5-15*x[:,1]**4+10*x[:,1]**3))/(20*nu1*rho1)*x[:,1]*(1-x[:, 1])),
("LeftUp",v1,v1,0),
("LeftDown",u0,u0,lambda x : (6*x[:,1]**5-15*x[:,1]**4+10*x[:,1]**3)/(20*nu0*rho0)*x[:,1]*(1-x[:, 1])),
("LeftDown",v0,v0,0),
("LeftDown",q0,q0,lambda x : 6*x[:,1]**5-15*x[:,1]**4+10*x[:,1]**3),
("LeftDown",u1,u1,lambda x : (1-(6*x[:,1]**5-15*x[:,1]**4+10*x[:,1]**3))/(20*nu1*rho1)*x[:,1]*(1-x[:, 1])),
("LeftDown",v1,v1,0),
("RightDown",v0,v0,0),
("RightDown",v1,v1,0),
("RightUp",v0,v0,0),
("RightUp",v1,v1,0),
]
fluid = fluid.fluid_problem("mesh.msh",g,[nu0*rho0,nu1*rho1],[rho0,rho1],epsilon,strong_boundaries,2)
ii = 0
t = 0
#set initial_condition
x = fluid.coordinates()
s = fluid.solution()
a = (x[:,0]+5)/10
s[:,2] = 6*x[:,1]**5-15*x[:,1]**4+10*x[:,1]**3
s[:,5] = 1-s[:,2]
fluid.export_vtk(outputdir,0,0)
ii = 0
tic = time.clock()
while t < tEnd :
#Fluid solver
fluid.implicit_euler(dt)
t += dt
s = fluid.solution()
x = fluid.coordinates()
vel = (s[:,0]+s[:,3]-1/(20*nu0*rho)*x[:,1]*(1-x[:, 1]))**2
print('Error', (vel.sum())**.5)
#Output files writting
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
fluid.export_vtk(outputdir, t, ioutput)
ii += 1
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.clock() - tic))
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