Need help implementing BCs for Streamfunction-vorticity formulation of steady Stokes between parallel plates

[rkailash]

Dear Community

I have a little more than a month’s experience in working with FEniCS and finite element techniques. I am trying to learn how to solve the steady Stokes flow problem between two parallel plates using the streamfunction (\psi)-vorticity(\omega) formulation, and compare against the solution obtained from the conventional velocity-pressure formulation. The problem geometry is as given below:

[Figure taken from G Linga’s tutorial handout].

The governing equations are as follows

\nabla^2\psi+\omega=0
\nabla^2\omega=0

where

u_x=\dfrac{\partial \psi}{\partial y}, u_y=-\dfrac{\partial \psi}{\partial x}, \omega=\left(\dfrac{\partial u_y}{\partial x}-\dfrac{\partial u_x}{\partial y}\right)

Defining v and q as the scalar test functions, I naively wrote the variational formulation as follows

\int_{x}\int_{y}\left[\left(\nabla \psi\right)\cdot\left(\nabla v\right)-\omega v\right]=0
\int_{x}\int_{y}\left(\nabla\omega\right)\cdot\left(\nabla q\right)=0

but I believe I might have missed something.

Once I have the BCs defined, my plan is to set up and solve the system of coupled PDEs on FEniCS.

The boundary conditions in terms of the velocity and pressure are quite easy to understand, but I am having trouble prescribing the appropriate BCs in terms of \psi and \omega. Could you please guide me/direct me to useful resources?

Thank You
Warm Regards

[Kei_Yamamoto]

Could you supply minimal working example of your code ?
What I would imagine is to use u and p to define steam function and vorticity, and thus you can directly use u to define BC, so I’m curious how you built your variational formulation in FEniCS.

[rkailash]

Dear @Kei_Yamamoto

I have not yet fully written the FEniCS code. I will write it once I have the following aspects cleared about the boundary conditions:

In terms of the velocity and pressure, the following conditions apply:

(i) \mathbf{u}=0 at y=0 [No slip]
(ii) \mathbf{u}=0 at y=d [No slip]
(iii) p_{\text{left}} at x=0 specified by user
(iv) p_{\text{right}} at x=l specified by user.

In other words, these are all Dirichlet BCs.

Now, in terms of the stream-function, I assume the no-slip BC would become

\mathbf{u}=0\implies u_x=\dfrac{\partial \psi}{\partial y}=0 and u_y=-\dfrac{\partial \psi}{\partial x}=0 at y=0 and y=d. Does this mean I could simply set \psi=c , a constant, at the top and bottom walls? Is it necessary to specify \psi at the left and right boundaries as well?

I also do not know how to specify the boundary conditions in the vorticity (\omega) at the boundaries based on the velocity and pressure information.

Once I get these doubts cleared, I think I will be in a better position to correctly set up the variational formulation and write the code.

Thank you for your time
Warm Regards

[rkailash]

Please find below the code I have written for this purpose. The weak variational formulation is provided above, in the original post. Since I do not know how to specify the boundary conditions for vorticity, I am obtaining a zero value for the solution everywhere. Any help is greatly appreciated.

## Solving Stokes flow in a rectangular pipe
## using streamfunction-vorticity formulation

from dolfin import *
import numpy as np

# Define domain and discretization
delta = 0.5
length = 10.0
diameter = 4.0
# Create mesh
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)

# Making a mark dictionary
# Note: the values should be UNIQUE identifiers. 
mark = {"generic": 0,
         "left": 1,
         "right": 2,
          "top" : 3,
      "bottom" : 4 }

subdomains = MeshFunction("size_t", mesh, 1)
subdomains.set_all(mark["generic"])

class Left(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], 0)

class Right(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], length)

class Top(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1],diameter)

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

top = Top()
top.mark(subdomains, mark["top"])
bottom = Bottom()
bottom.mark(subdomains, mark["bottom"])


geom_name="geometry.pvd"
File(geom_name) << subdomains

# Define function spaces
PSI = FiniteElement("Lagrange", triangle, 1)
OMEGA = FiniteElement("Lagrange", triangle, 1)
W = FunctionSpace(mesh, PSI*OMEGA)

# Define variational problem
(psi, omega) = TrialFunctions(W)
(v, q) = TestFunctions(W)

# Surface normal
n = FacetNormal(mesh)
# Pressures. First define the numbers (for later use):
omega_top = (diameter/4.)
omega_bottom = -(diameter/4.)
# ...and then the DOLFIN constants:
om_top = Constant(omega_top)
om_bottom = Constant(omega_bottom)
# Body force:
force = Constant((0.0))

a = inner(grad(psi), grad(v))*dx - (omega*v)*dx + inner(grad(omega), grad(q))*dx
L = inner(force, v) * dx

##Applying boundary conditions
noslip = Constant((0.0))
bc_psi_top = DirichletBC(W.sub(0), noslip, subdomains, mark["top"])
bc_psi_bottom = DirichletBC(W.sub(0), noslip, subdomains, mark["bottom"])
#bc_omega_top = DirichletBC(W.sub(1), om_top, subdomains, mark["top"])
#bc_omega_bottom = DirichletBC(W.sub(1), om_bottom, subdomains, mark["bottom"])

bcs = [bc_psi_top, bc_psi_bottom]
#
## Solve variational problem
w = Function(W)
solve(a == L, w, bcs)

# Split using deep copy
(psi, omega) = w.split(True)

## Save solution in VTK format
psi_file = File("psi.pvd")
psi_file << psi
omg_file = File("omega.pvd")
omg_file << omega

vel_x=project(psi.dx(1),W.sub(0).collapse())
vel_y=-project(psi.dx(0),W.sub(0).collapse())

vel_profile=as_vector((vel_x,vel_y))

vel_x_file = File("velocity_x.pvd")
vel_x_file << vel_x

#vel_y_file = File("velocity_y.pvd")
#vel_y_file << vel_y

Solution for \psi as obtained from code: