Commit ea88fc48 authored by matthieu's avatar matthieu
Browse files

correction Brinkman.py + ajout des entrées pour la stabilisation

parent ba6b4425
......@@ -41,7 +41,7 @@ class Brinkman :
val = np.ndarray(sol.shape + (self.D,))
val[..., _U_, :] = -solgrad[..., _U_, :] * self.mu
val[..., _p_, :] = 0#self.vs# + sol[..., _U_]
val[..., _p_, :] = -1e-8 * solgrad[..., _p_, :]
val[..., _p_, :] = -self.epsilon * solgrad[..., _p_, :]
if self.ns :
val[..., _u_, :] += - sol[..., [_u_]] * sol[..., _U_] / self.ca * self.rhof
val[..., _v_, :] += - sol[..., [_v_]] * sol[..., _U_] / self.ca * self.rhof
......@@ -55,27 +55,35 @@ class Brinkman :
_U_ = range(self.D)
val = np.zeros_like(sol)
if self.D == 2 :
#diametre
d = 2 * (self.pf.volume[0, 0] / np.pi)**0.5
#aire projetee sur le plan perpendiculaire a l'ecoulement
Ap=d
else :
d = 2 * (3./(4 * np.pi) * self.pf.volume[0, 0])**(1./3)
#aire projetee sur le plan perpendiculaire a l'ecoulement
Ap=np.pi*(d/2)**2
#difference des vitesses solide et fluide u_p-u_f
U = self.vs - sol[..., _U_]/self.ca
normU = np.sqrt(self.dot(U, U))[...,np.newaxis]
#Reynolds/|u_p-u_f|
Re_O_u = d * self.ca/self.mu * self.rhof
Cd_u = (0.63 * normU + 4.8/(Re_O_u**0.5))**2
#coefficient du drag
Cd_u = (0.63 * normU**0.5 + 4.8/(Re_O_u**0.5))**2
f = self.ca**(-1.8)
F = solgrad[..., _p_, :]
GoU = f * Cd_u * self.rhof / 2 * d
#mass = self.pf.volume * 3/2 * self.rhop
#drag
GoU = f * Cd_u * self.rhof / 2 * Ap
aext = self.g * (self.rhop - self.rhof) /self.rhop
UImp = U + aext * self.stab
if fonpart is not None:
fonpart[...] = (-F*self.pf.volume - GoU*UImp) / (1 + self.stab/self.pf.mass * GoU) + aext * self.pf.mass
val[..., _U_] = (F*self.pf.volume + GoU* UImp) / (1 + self.stab/self.pf.mass * GoU)
fonpart[...] = (-F * self.pf.volume - GoU*UImp) / (1 + self.stab/self.pf.mass * GoU) + aext * self.pf.mass
val[..., _U_] = (F *self.pf.volume + GoU* UImp*self.pf.volume) / (1 + self.stab/self.pf.mass * GoU)
val[..., _p_] = 0
return val
def __init__(self, mu, pf, g, rho, rhop, ns = False, stab = 0, D = 2) :
def __init__(self, mu, pf, g, rho, rhop, epsilon, ns = False, stab = 0, D = 2) :
self.D = D
self.ns = ns
self.rhof = rho
......@@ -84,6 +92,7 @@ class Brinkman :
self.mu = mu
self.pf = pf
self.g = np.array(g)
self.epsilon = epsilon
def particleInteractionTerm(self, meshJ, solution, res, jac) :
_p_ = self.D
......
......@@ -8,7 +8,7 @@ from .legendre import *
import time
class FluidSolver :
def __init__(self, mesh, dt, rhop, rho, mu, g, imex=True) :
def __init__(self, mesh, dt, rhop, rho, mu, g, epsilon, TypeEMu, TypeEMp, imex=True) :
if isinstance(mesh, str):
mesh = Mesh(mesh)
self._mesh = mesh
......@@ -24,18 +24,19 @@ class FluidSolver :
g = np.array(g),
rho = rho,
rhop = rhop,
epsilon = epsilon,
D = self._jac.dim,
stab = dt)
d = DofManager(mesh)
if self._jac.dim == 2 :
d.add_field(TriangleP1)
d.add_field(TriangleP1)
d.add_field(TriangleP1)
d.add_field(TypeEMu)
d.add_field(TypeEMu)
d.add_field(TypeEMp)
elif self._jac.dim == 3:
d.add_field(TetrahedronP1Mini)
d.add_field(TetrahedronP1Mini)
d.add_field(TetrahedronP1Mini)
d.add_field(TetrahedronP1)
d.add_field(TypeEMu)#TetrahedronP1Mini
d.add_field(TypeEMu)
d.add_field(TypeEMu)
d.add_field(TypeEMp)#TetrahedronP1
else :
raise ValueError("invalid mesh dimension : %i" % self._jac.dim)
self._dt = dt
......
......@@ -10,7 +10,7 @@ class ScontactInterface :
p.velocity()[:, :] = 0;
self._r = p.particles()[:, 0]
self._volume = np.pi * p.particles()[:, [0]]**2 * 2./3 ## 2./3 = 2d-3d correction
self._volume = np.pi * p.particles()[:, [0]]**2
self._mass = p.particles()[:, [1]]
self._position = p.position()
self._velocity = p.velocity()
......@@ -38,7 +38,7 @@ class ScontactInterface :
def iterate(self, dt, forces) :
self._p.externalForces()[:] = forces
self._p.iterate(np.min(self._r), dt, 1e-5)
self._p.iterate(np.min(self._r), dt, 1e-7)
def write(self, odir, i, t) :
outputname = "%s/part-%05i" % (odir, i)
......
......@@ -4,6 +4,7 @@ from pyfefp import scontact2Interface
import numpy as np
import time
import os
from pyfefp.legendre import *
def genInitialPosition(filename, r, lx, ly, rhop) :
p = scontact2.ParticleProblem()
......@@ -18,14 +19,14 @@ def genInitialPosition(filename, r, lx, ly, rhop) :
p.position()[:, 1] += 4
p.write(filename)
dt = 0.0125
tEnd = 15000
dt = 0.005
tEnd = 2
outputdir = "output"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
rhop=1100
d = 0.04
rhop=2500
d = 0.005
mu = 0.001
rho = 1000
ii = 0
......@@ -47,7 +48,7 @@ print("r = %g, m = %g\n" % (p.r()[0], p.mass()[0]))
print("dtmax = %g Re_p/u = %g" % (p.mass()[0] / (150 * mu * d), d * rho * 0.3/mu))
print("RHOP = %g" % rhop)
outf = 5
fluid = pyfefp.FluidSolver("mesh.msh", dt, rhop=rhop, rho=rho, g = g, mu=mu, imex=False)
fluid = pyfefp.FluidSolver("mesh.msh", dt, rhop=rhop, rho=rho, mu=mu, g = g, epsilon=1e-02, imex=False, TypeEMu=TriangleP1, TypeEMp=TriangleP1)
fluid.add_boundary_condition("Box", 0, 0.)
fluid.add_boundary_condition("Box", 1, 0.)
fluid.add_boundary_condition("Top", 0, 0.)
......@@ -67,7 +68,7 @@ while t < tEnd :
vn = p.velocity() + forces * dt / p.mass()
vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
nsub = max(1, int(np.ceil((vmax * dt * 2)/min(p.r()))))
print("NSUB", nsub, "VMAX * dt", vmax * dt, "r", 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())
......
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