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)
return nu*(grad(X)+grad(X).T)
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()
Conv = inner(grad(u)*u, phi)*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)
```

Thank you in advance