Stokes in axisimmetrical cylindrical coordinates - DOLFINX

dolfinx
2

First of all, you should use the collapsed spaces for your boundary conditions:

# Function spaces
P2 = ufl.VectorElement("Lagrange", msh.ufl_cell(), 2)
P1 = ufl.FiniteElement("Lagrange", msh.ufl_cell(), 1)
TH = P2 * P1
W = FunctionSpace(msh, TH)
V, _ = W.sub(0).collapse()
Q, _ = W.sub(1).collapse()
# Functions to geometrically locate subsets
def null_v_boundary(x):
    return np.isclose(np.abs(x[1]), 0.5)

def inlet_boundary(x):
    return np.isclose(x[0], 0.0)

def outlet_boundary(x):
    return np.isclose(x[0], 2.0)

# Inflow velocity
def in_velocity_expression(x):
    return np.vstack((v_max*(np.ones(x[1].shape[0])-4*x[1]**2), np.zeros(x.shape[1])))

# Boundary conditions
# Null velocity
null_v = Function(V)
facets = locate_entities_boundary(msh, 1, null_v_boundary)
dofs = locate_dofs_topological((W.sub(0), V), 1, facets)
bc0 = dirichletbc(null_v, dofs, W.sub(0))

# Null pressure
null_p = Function(Q)
facets = locate_entities_boundary(msh, 1, outlet_boundary)
dofs = locate_dofs_topological((W.sub(1), Q), 1, facets)
bc1 = dirichletbc(null_p, dofs, W.sub(1))

# Inflow velocity condition
in_velocity = Function(V)
in_velocity.interpolate(in_velocity_expression)
facets = locate_entities_boundary(msh, 1, inlet_boundary)
dofs = locate_dofs_topological((W.sub(0), V), 1, facets)
bc2 = dirichletbc(in_velocity, dofs, W.sub(0))

Secondly,as your mesh contains y=0, which should be treated carefully as you both multiply and divide by x[1]. If you say that there are existing implementations in older versions of FEniCS that work, please refer us to them and point out differences in your implementation and theirs.

1 Like