Problem with velocity

1

Hi there
i have a problem with the velocity . i dont undersatand why the velocity is 0.0
the full code is:

import matplotlib.pyplot as plt
from fenics import *
from dolfin import *
import numpy as np
from mshr import*

p1 = Point(np.array([1.0,  0.0]))
p2 =  Point(np.array([2.0,  -0.2]))




p3 = Point(np.array([6.0,  0.2]))


c0 = Circle(p1,1)
R = Rectangle (p2,p3)
vessie = R+c0

mesh = generate_mesh (vessie , 64)


plt.figure(figsize=(20, 8))
plot(mesh, title="Mesh")


mark = {"generic": 0,
"wall1": 1,
 "wall2": 2,      
"right": 3 }

#Fonction pour le maillage
subdomains = MeshFunction("size_t", mesh, 1)
subdomains.set_all(mark["generic"])

#Partie droite des conditions aux bords
class Right(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], 6.0)

#Parois

class Wall1(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near (x[1], -0.2) or near (x[1],0.2 ) 
             

class Wall2(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near (x[1]*x[1]+x[0]*x[0] ,1)

#application parois
wall1 = Wall1() 
wall2 = Wall2()
wall1.mark(subdomains, mark["wall1"])
wall2.mark(subdomains, mark["wall2"])



#application a droite
right = Right()
right.mark(subdomains, mark["right"])




#définitioon de l'espace des fonctions 
P2 = VectorElement('P', triangle, 2)
P1 = FiniteElement('P', 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_wall2 = 2000.00

p_right = 0.0

# pression des deux coté  codées pour DOLFIN:

pressure_wall2 = Constant(p_wall2)
pressure_right = Constant(p_right)

#force
force = Constant((0.0, 0.0))

#Utilisation de DOLFIN
a = inner(grad(u), grad(v))*dx - p*div(v)*dx + q*div(u)*dx
L = inner(force, v) * dx \
- p_right * inner(n, v) * ds(mark["right"]) - p_wall2 * inner(n, v) * ds(mark["wall2"])

#conditions de vitesse aux bords de la parois
noslip = Constant((0.0, 0.0))
bcu_wall = DirichletBC(W.sub(0), noslip, subdomains, mark["wall1"])





#Résolution
w = Function(W)
#solve(a == L, w,  bcu)
solve(a == L, w,  bcu_wall)

# Split using deep copy
(u, p) = w.split(True)

u

u.vector().get_local()

# Plot solution
plt.figure(figsize=(20, 8))
plot(u, title="Velocity")

plot§

2

There are two problems with your code.

  1. The two domains are not connected to eachother.
  2. The marker function for “wall2” is not marking anything.

You can verify this by saving the meshfunction to pvd and visualize it in paraview.

File("subdomains.pvd") << subdomains

subdomains
subdomains2954×1266 500 KB

3

ok i have conected the two domains
i have used the ''wall2" function to describe the circle. maybe i have done a mistake in the définition

4

i have change wall 2 to:
class Wall2(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near (x[1]x[1] + ((x[0]-1)(x[0]-1)) - 1, 0 )

5

And? Did that solve your problem ?

6

no but i know where is the problem.
i have put a pressure auround all the circle.
now i want to put it in a portion of the circle.
i want to remove the part,where de circle meet the rectangle. between y= -0.2 and y = 0.2