I want to simulate a laminar flow, but i have not a parabolic profile of the vélocity.
The full code is bellow
import matplotlib.pyplot as plt
from dolfin import *
import numpy as np
#Définition du domaine rectangulaire
delta = 0.5
length = 12
diameter = 2
#discretisation du domaine
p0 = Point(np.array([0.0, 0.0]))
p1 = Point(np.array([length, diameter]))
nx = int(length/delta)
ny = int(diameter/delta)
mesh = RectangleMesh(p0, p1, nx, ny)
#visualisation du maillage
plot(mesh, title=“Mesh”)
#Numéro identification des bords
mark = {“generic”: 0,
“wall”: 1,
“left”: 2,
“right”: 3 }
#Fonction pour le maillage
subdomains = MeshFunction(“size_t”, mesh, 1)
subdomains.set_all(mark[“generic”])
#Partie gauche des conditions aux bords
class Left(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], 0)
#Partie droite des conditions aux bords
class Right(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[1], 0)
#Parois
class Wall(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], 1) and near (x[1],1)
#application a gauche
left = Left()
left.mark(subdomains, mark[“left”])
#application a droite
right = Right()
right.mark(subdomains, mark[“right”])
#application parois
wall = Wall()
wall.mark(subdomains, mark[“wall”])
#définitioon de l’espace des fonctions
P2 = VectorElement(“CG”, triangle, 2)
P1 = FiniteElement(“CG”, triangle, 1)
TH = P2 * P1
W = FunctionSpace(mesh, TH)
Definition du problème variationel
(u, p) = TrialFunctions(W)
(v, q) = TestFunctions(W)
#integral de volume et de surface
dx = Measure(“dx”, domain=mesh, subdomain_data=subdomains)
ds = Measure(“ds”, domain=mesh, subdomain_data=subdomains)
Surface normale
n = FacetNormal(mesh)
Pression des deux cotés
p_left = 1.0
p_right = 0.0
pression des deux coté codées pour DOLFIN:
pressure_left = Constant(p_left)
pressure_right = Constant(p_right)
#force
force = Constant((0.0, 0.0))
#Utilisation de DOLFIN
a = inner(grad(u), grad(v))dx - pdiv(v)dx + qdiv(u)*dx
L = inner(force, v) * dx \
- pressure_left * inner(n, v) * ds(mark[“left”]) \
- pressure_right * inner(n, v) * ds(mark[“right”])
#conditions de vitesse aux bords de la parois
noslip = Constant((0.0, 0.0))
bcu_wall = DirichletBC(W.sub(0), noslip, subdomains, mark[“wall”])
#conditions de pression a gauche
bcp_left= DirichletBC(W.sub(1), p_left, subdomains, mark[“left”])
#conditions de presson a droite
bcp_right= DirichletBC(W.sub(1), p_right, subdomains, mark[“right”])
#condition au bord vitesse en un seul vecteur
bcu = [bcu_wall,bcp_left ,bcp_right]
#Résolution
w = Function(W)
solve(a == L, w, bcu)
Split using deep copy
(u, p) = w.split(True)
Plot solution
plot(u, title=“Velocity”)
plot(p, title=“Pressure”)
#interactive()
Magnitude function pour vidulaiser le module de la vitesse
def magnitude(vec):
return sqrt(vec**2)
plot(magnitude(u), “Speed”)
Save solution in VTK format
ufile_pvd = File(“velocity.pvd”)
ufile_pvd << u
pfile_pvd = File(“pressure.pvd”)
pfile_pvd << p
Plot solution
plt.figure()
plot(u, title=“velocity”)
plt.figure()
plot(p, title=“pressure”)
Display plots
plt.show()