Integral constraint for a Poiseuille flow

1

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

2

I suggest you to add error messages as well.

3

Thank you for the suggestion. I got the following error

Can't add expressions with different shapes.
Traceback (most recent call last):
    Fc = (u - Constant(C0)) * ksi * dx + q*phi*ds(inlet)

In this topic they also added the following command but I do not know why

u.assign(Constant([C0,C0]))