# Integral constraint for a Poiseuille flow

Hello everyone,

I want to solve a 2D Poiseuille flow with a flow rate constraint. I follow the steps proposed in this topic but I am getting errors. I want to solve the Navier-Stokes equation in a 2D rectangular domain. I would like to have a fully-developed flow at the inlet so I want to impose a constraint regarding the total volumetric flow rate. The unknowns of the problem are the velocity U(Ux,Uy), pressure P and I also defined a “Real” element for flow rate Q. At the wall I impose no-slip condition and at the inlet/outlet I include the term from the integration by parts. The code is converged for Dirichlet BCs at the inlet. Here is a sample of my code:

``````*# Define function spaces *
V           = VectorElement('CG'  ,mesh.ufl_cell(),2)
P           = FiniteElement('CG'  ,mesh.ufl_cell(),1)
Q           = FiniteElement('Real',mesh.ufl_cell(),0)
VPQ         = FunctionSpace(mesh,MixedElement([V,P,Q]))
phi,eta,ksi = TestFunctions(VPQ)

*# Set datum pressure*
def PinPoint(x, on_boundary ):
tol = 1.e-4
pointPres = near( x[0], 0.995833    , tol) and near(x[1],  0.005  , tol)

if (pointPres is True):
return (pointPres)

*#BCs*
Pref = Constant(0.0)
ns2  = Constant((0.0,0.0))

vFluidWall  = DirichletBC(VPQ.sub(0), ns2 ,boundaries, wall  )

bcp = DirichletBC(VPQ.sub(1),Pref,PinPoint,'pointwise')

bcs     = [bcp , vFluidWall ]
bcs_    = [bcp , vFluidWall ]

*# Create functions*
W         = Function(VPQ)
[bc.apply(W.vector()) for bc in bcs]
u,p,q     = split(W)

*#	stresses*
def tau(X):
nu = Constant(0.1)

def cauchy(X,P):
return -P*Identity(len(X)) + tau(X)

#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
# 						*weak forms*
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-

dx = Measure("dx")(subdomain_data=subdomains)
ds = Measure("ds")(subdomain_data=boundaries) # outer boundary
dS = Measure("dS")(subdomain_data=boundaries) # inner boundary

n       =   FacetNormal(mesh)
Diff    =   inner( cauchy(u,p) , grad(phi))*dx()
outFBC  = - inner( dot( cauchy(u,p),n ),phi )*ds(outlet)
inFBC   = - inner( dot( cauchy(u,p),n ),phi )*ds(inlet)

Cont  = eta*div(u)*dx()

*# Conservation term*
Co = 0.01
Fc = (u - Constant(C0)) * ksi * dx + q*phi*ds(inlet)

Res    = Diff + Cont + Conv + outFBC + inFBC + Fc
Jac    = derivative(Res,W)
Corr   = Function(VPQ)

*# Solve*
NewtonRaphson(Res,Jac,bcs_,W,Corr,rank)
``````

``````Can't add expressions with different shapes.
`u.assign(Constant([C0,C0]))`