Stokes equations from Cartesian to Cylindrical, results mismatch

variational formulation
1

Hello,
I am solving a test example, which is Stokes equations (slow flow equations) in Cartesian coordinates and cylindrical coordinates, so that I can compare the two and make sure I am changing coordinate system appropriately. For the simple conditions I have imposed, I expect a quadratic profile for the velocity, which I successfully recover for the Cartesian case, but not the cylindrical case. I believe the code should be the same between the two, apart from the stress tensor definition and the integration measure, which (I think) I did correctly. I’ve been stuck on this for days, I would really appreciate any help/comments.

Code for Cartesian case:

 from dolfin import *

# Define mesh and geometry - We solve for half of the domain we need, and impose symmetry

mesh = RectangleMesh(Point(0, 0), Point(1, 1), 100, 100)
n = FacetNormal(mesh)

# Define Taylor--Hood function space W
V = VectorElement("CG", triangle, 2)
Q = FiniteElement("CG", triangle, 1)
W = FunctionSpace(mesh, MixedElement([V, Q]))

# Define Function and TestFunction(s)
w = Function(W)
(u, p) = split(w)
(v, q) = split(TestFunction(W))

# Define the viscosity and bcs

mu = Constant(1)
u_c = Constant(0.0)
p0 = Constant(8.0)
p1 = Constant(0.0)


inflow = 'near(x[1], 1.0)'
wall1 = 'near(x[0], 1.0)'
wall2 = 'near(x[0], 0.0)'
outflow = 'near(x[1], 0.0)'
bcu_inflow = DirichletBC(W.sub(1), p0, inflow)
bcu_wall1 = DirichletBC(W.sub(0), (0.0, u_c), wall1)
bcu_wall2 = DirichletBC(W.sub(0), (0.0, u_c), wall2)
bcu_outflow = DirichletBC(W.sub(1), p1, outflow)
bcs = [bcu_inflow, bcu_wall1, bcu_outflow, bcu_wall2]
# Define stress tensor
x = SpatialCoordinate(mesh)

def epsilon(v):
    return sym(as_tensor([[v[0].dx(0), v[0].dx(1), 0],
                          [v[1].dx(0), v[1].dx(1), 0],
                          [0, 0, 0]]))
# stress tensor
def sigma(v, p):
    return 2*mu*epsilon(v)-Id(p)

def Id(p):
    return as_tensor([[p, 0, 0],
                      [0, p, 0],
                     [0, 0, p]])

# Define the variational form
f = Constant((0, -1))



colors = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
colors.set_all(0)  # default to zero
CompiledSubDomain("near(x[0], 0.0)").mark(colors, 5)
CompiledSubDomain("near(x[1], 1.0)").mark(colors, 1)
CompiledSubDomain("near(x[0], 1.0)").mark(colors, 3)  # wall
CompiledSubDomain("near(x[1], 0.0)").mark(colors, 4)  # outflow

# Create the measure
ds = Measure("ds", subdomain_data=colors)

a1 = (inner(sigma(u, p), epsilon(v))) * dx
a2 = (- div(u) * q - dot(f, v)) * dx
F = a1 + a2

# Solve problem
solve(F == 0, w, bcs)

# Plot solutions
(u, p) = w.split()
File("StokesWithBC/velocitySym.pvd") << u
File("StokesWithBC/pressureSym.pvd") << p

Code for the 2D cylindrical case (almost exactly the same, apart from integration measure and stress tensor)

from dolfin import *
import numpy as np

# Define mesh and geometry - We solve for half of the domain we need, and impose symmetry
mesh = RectangleMesh(Point(0, 0), Point(1, 1), 100, 100)
n = FacetNormal(mesh)

# Define Taylor--Hood function space W
V = VectorElement("CG", triangle, 2)
Q = FiniteElement("CG", triangle, 1)
W = FunctionSpace(mesh, MixedElement([V, Q]))

# Define Function and TestFunction(s)
w = Function(W)
(u, p) = split(w)
(v, q) = split(TestFunction(W))

# Define the viscosity and bcs

mu = Constant(1)
u_c = Constant(0.0)
p0 = Constant(8.0)
p1 = Constant(0.0)

# Note, x[0] is r and x[1] is x, and x[1] == 0 is the bottom.
inflow = 'near(x[1], 1.0)'
wall1 = 'near(x[0], 1.0)'
wall2 = 'near(x[0], 0.0)'
outflow = 'near(x[1], 0.0)'
bcu_inflow = DirichletBC(W.sub(1), p0, inflow)
bcu_wall1 = DirichletBC(W.sub(0), (0.0, u_c), wall1)
bcu_wall2 = DirichletBC(W.sub(0), (0.0, u_c), wall2)
bcu_outflow = DirichletBC(W.sub(1), p1, outflow)
bcs = [bcu_inflow, bcu_wall1, bcu_outflow, bcu_wall2]
# Define stress tensor
# epsilon = sym(grad(u))
x = SpatialCoordinate(mesh)


# symmetric gradient
def epsilon(v):
    return sym(as_tensor([[v[0].dx(0), 0, v[0].dx(1)],
                          [0, v[0] / x[0], 0],
                          [v[1].dx(0), 0, v[1].dx(1)]]))


def Id(p):
    return as_tensor([[p, 0, 0],
                      [0, p, 0],
                     [0, 0, p]])


# stress tensor
def sigma(v, p):
    return 2*mu*epsilon(v)-Id(p)


# Define the variational form
f = Constant((0, -1))



colors = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
colors.set_all(0)
CompiledSubDomain("near(x[0], 0.0)").mark(colors, 5)
CompiledSubDomain("near(x[1], 1.0)").mark(colors, 1)
CompiledSubDomain("near(x[0], 1.0)").mark(colors, 3)  # wall
CompiledSubDomain("near(x[1], 0.0)").mark(colors, 4)  # outflow

# Create the measure
ds = Measure("ds", subdomain_data=colors)

a1 = (inner(sigma(u, p), epsilon(v))) * x[0] * dx
# I add the extra u[0] *q*dx because I am using the 'div' operator in Python, which has this extra term in cylindrical coordinates
a2 = (- div(u) * q - dot(f, v)) * x[0] * dx - u[0] * q * dx 
F = a1 + a2

# Solve problem
solve(F == 0, w, bcs)

# Plot solutions
(u, p) = w.split()
File("StokesWithBC/velocityCyl.pvd") << u
File("StokesWithBC/pressureCyl.pvd") << p

Cartesian Coordinates

Cartesian
Cartesian1486×512 28.3 KB

Cylindrical Coordinates

Cylindrical
Cylindrical1486×512 19.4 KB

1 Like
2

I should note that if I remove the ‘x[0]’ factors in variational form, then I get matching answers. But, the integration measure has to change when you change coordinate system, right? And the measure changes to r dr dz , so I need the x[0] factor. Maybe my Geometry is wrong ? Although in 2D cylindrical coordinates, a rectangle should represent the right geometry… Anyway, I’ve been going in circles for a while. Let me know if you have any ideas!

3

Why would you expect the cylindrical solution to be the same as the 2D Cartesian solution?

One is a 2D problem which could be considered a slice of an infinite slab. And assuming you’re solving in the (r, z) plane, the other is an azimuthally symmetric slice of a 3D cylinder.

See for example here.

1 Like
4

Hi,
Thanks for replying. You’re right, of course. It shouldn’t be equal, I had clearly misunderstood the problem. These past few hours I have tried to solve the Stokes problem imposing pressure at the top and bottom, as well as velocity on the walls. Since it’s unidirectional flow, it’s easy to calculate an analytical solution. My Cartesian coordinates converges to the analytical solution, whereas the cylindrical case does not (they have slightly different analytical solutions, as you indicated). This is the simplest example I can think of, and I can’t get it to work. Is it obvious why? I attach the code. Note that I have two options for the variational formulation, I can either use the Python/Dolfin inbuilt functions, or I can put each component myself. The way I have written it, I get the same answer whichever formulation I use, which suggests that there is something conceptual I don’t understand.

from dolfin import *
import numpy as np

# Define mesh and geometry - We solve for half of the domain we need, and impose symmetry
mesh = RectangleMesh(Point(0, 0), Point(1, 1), 100, 100)
n = FacetNormal(mesh)

# Define Taylor--Hood function space W
V = VectorElement("CG", triangle, 2)
Q = FiniteElement("CG", triangle, 1)
W = FunctionSpace(mesh, MixedElement([V, Q]))

# Define Function and TestFunction(s)
w = Function(W)
(u, p) = split(w)
(v, q) = split(TestFunction(W))

# Define the viscosity and bcs

mu = Constant(1)

u_c = Constant(0.0)
p0 = Constant(-9.0)
p1 = Constant(0.0)

# Note, x[0] is r and x[1] is x, and x[1] == 0 is the bottom.
inflow = 'near(x[1], 1.0)'
wall1 = 'near(x[0], 1.0)'
wall2 = 'near(x[0], 0.0)'
outflow = 'near(x[1], 0.0)'
bcu_inflow = DirichletBC(W.sub(1), p0, inflow)
bcu_wall1 = DirichletBC(W.sub(0), (0.0, 0.0), wall1)
bcu_wall2 = DirichletBC(W.sub(0), (0.0, -2.0), wall2)
bcu_outflow = DirichletBC(W.sub(1), p1, outflow)
bcs = [bcu_inflow, bcu_wall1, bcu_outflow, bcu_wall2]

x = SpatialCoordinate(mesh)


# symmetric gradient
def epsilon(v):
    return sym(as_tensor([[v[0].dx(0), 0, v[0].dx(1)],
                          [0, v[0] / x[0], 0],
                          [v[1].dx(0), 0, v[1].dx(1)]]))


# stress tensor
def sigma(v, p):
    return 2*mu*epsilon(v)-p*Identity(3)

def sigmabc(v, p):
    return mu * (nabla_grad(v) + nabla_grad(v).T) - p*Identity(2)


# Define the variational form
f = Constant((0, -1))

colors = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
colors.set_all(0)
CompiledSubDomain("near(x[0], 0.0)").mark(colors, 5)
CompiledSubDomain("near(x[1], 1.0)").mark(colors, 1)
CompiledSubDomain("near(x[0], 1.0)").mark(colors, 3)  # wall
CompiledSubDomain("near(x[1], 0.0)").mark(colors, 4)  # outflow

# Create the measure
ds = Measure("ds", subdomain_data=colors)
#ap1 = (2*u[0].dx(0)*v[0].dx(0) + (u[0].dx(1) + u[1].dx(0))*(v[0].dx(1) + v[1].dx(0)) + 2*(u[0] / x[0]) * (v[0] / x[0]) + 2*u[1].dx(1)*v[1].dx(1)) * x[0] * dx
#ap2 = (-p*v[0].dx(0) -p*(v[0] / x[0]) - p*v[1].dx(1)) * x[0] * dx
#a1 = ap1 + ap2
a1 = (inner(sigma(u, p), epsilon(v))) * x[0] * dx
a2 = (- div(u) * q - dot(f, v)) * x[0] * dx - u[0] * q * dx
#a2 = (- ((1/x[0])*(x[0]*u[0]).dx(0) + u[1].dx(1)) * q - f[1]*v[1]) * x[0] * dx
b1 = - dot(dot(sigmabc(u, p), n), v) * x[0] * ds
#b3 = - dot(dot(sigmabc(u, p), v), n) * x[0] * ds(3)
#b4 = - asp*dot(dot(sigmabc(u, p), v), n) * x[0] * ds(4)
F = a1 + a2 + b1

# Solve problem
solve(F == 0, w, bcs)

# Plot solutions
(u, p) = w.split()

# Compute error
u_e = Expression(('0', '2*(x[0]*x[0]-1)'), degree=2)
W2 = FunctionSpace(mesh, V)
u_e = interpolate(u_e, W2)
u_r = interpolate(u, W2)
error = np.abs(np.array(u_e.vector()) - np.array(u_r.vector())).max()

print('error = ', error)
print('max u:', np.array(u_r.vector()).max())

File("StokesWithBC/velocityCyl.pvd") << u
File("StokesWithBC/pressureCyl.pvd") << p
File("StokesWithBC/velocityAnalyticCyl.pvd") << u_e

The error that I get is 1.244, which is very large. Even worse, when I refine the mesh it gets larger rather than smaller, so I have clearly done it wrong. If you have any ideas, I’d really appreciate it! In any case, thanks for your earlier response.

5

Turns out my maths was wrong, which is unfortunate. I have fixed it, so if anyone has a similar issue, this is the working code:

from dolfin import *
import numpy as np

# Define mesh and geometry - We solve for half of the domain we need, and impose symmetry
mesh = RectangleMesh(Point(0, 0), Point(1, 1), 100, 100)
n = FacetNormal(mesh)

# Define Taylor--Hood function space W
V = VectorElement("CG", triangle, 2)
Q = FiniteElement("CG", triangle, 1)
W = FunctionSpace(mesh, MixedElement([V, Q]))

# Define Function and TestFunction(s)
w = Function(W)
(u, p) = split(w)
(v, q) = split(TestFunction(W))

# Define the viscosity and bcs

mu = Constant(1)

u_c = Constant(0.0)

# Analytic expression for the velocity 
u_in = Expression(('0', '-2*(1-x[0]*x[0])'), degree=2)

# Note, x[0] is r and x[1] is x, and x[1] == 0 is the bottom.
inflow = 'near(x[1], 1.0)'
wall1 = 'near(x[0], 1.0)'
wall2 = 'near(x[0], 0.0)'
outflow = 'near(x[1], 0.0)'
bcu_inflow = DirichletBC(W.sub(0), u_in,  inflow)
bcu_wall1 = DirichletBC(W.sub(0), (0.0, 0.0), wall1)
bcu_wall2 = DirichletBC(W.sub(0), (0.0, -2.0), wall2)
bcs = [bcu_inflow, bcu_wall1, bcu_wall2]

x = SpatialCoordinate(mesh)


# symmetric gradient
def epsilon(v):
    return sym(as_tensor([[v[0].dx(0), 0, v[0].dx(1)],
                          [0, v[0] / x[0], 0],
                          [v[1].dx(0), 0, v[1].dx(1)]]))

#gradient
def epsilon2(v):
    return as_tensor([[v[0].dx(0), 0, v[0].dx(1)],
                          [0, v[0] / x[0], 0],
                          [v[1].dx(0), 0, v[1].dx(1)]])


# stress tensor
def sigma(v, p):
    return 2*mu*epsilon(v)-p*Identity(3)


def div_cyl(v):
    return (1/x[0])*(x[0]*v[0]).dx(0) + v[1].dx(1)


# analytic expression of what the stress tensor should be with my analytic velocity
def sigmabc2(v):
    return 2*mu*sym(as_tensor([[0, 0],
                               [4*x[0], 0]]))


# Define the variational form
f = Constant((0, -1))

colors = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
colors.set_all(0)
CompiledSubDomain("near(x[0], 0.0)").mark(colors, 5)
CompiledSubDomain("near(x[1], 1.0)").mark(colors, 1)
CompiledSubDomain("near(x[0], 1.0)").mark(colors, 3)  # wall
CompiledSubDomain("near(x[1], 0.0)").mark(colors, 4)  # outflow
vext = as_vector([v[0], 0, v[1]])
next = as_vector([n[0], 0, n[1]])
# Create the measure
ds = Measure("ds", subdomain_data=colors)

#Note: both the commented and un-commented version of the variational form works, one is just written explicitly, and the other one uses the Dolfin functions.   

#ap1 = (2*u[0].dx(0)*v[0].dx(0) + (u[0].dx(1) + u[1].dx(0))*(v[0].dx(1) + v[1].dx(0)) + 2*(u[0] / x[0]) * (v[0] / x[0]) + 2*u[1].dx(1)*v[1].dx(1)) * x[0] * dx
#ap2 = (-p*v[0].dx(0) -p*(v[0] / x[0]) - p*v[1].dx(1)) * x[0] * dx
#a1 = ap1 + ap2

# Here I can use either the symmetric gradient or the non-symmetric gradient.
a1 = (inner(sigma(u, p), epsilon(v))) * x[0] * dx
a2 = (- div_cyl(u) * q - dot(f, v)) * x[0] * dx
#a2 = (- ((1/x[0])*(x[0]*u[0]).dx(0) + u[1].dx(1)) * q - f[1]*v[1]) * x[0] * dx
b1 = - dot(dot(sigmabc2(u), v), n) * x[0] * ds(4)


F = a1 + a2 + b1

# Solve problem
solve(F == 0, w, bcs)

# Plot solutions
(u, p) = w.split()

# Compute error
u_e = Expression(('0', '-2*(1-x[0]*x[0])'), degree=2)
W2 = FunctionSpace(mesh, V)
u_e = interpolate(u_e, W2)
u_r = interpolate(u, W2)
error = np.abs(np.array(u_e.vector()) - np.array(u_r.vector())).max()

print('error velocity = ', error)
print('max u:', np.array(u_r.vector()).max())

p_e = Expression('7*x[1]', degree=1)
W3 = FunctionSpace(mesh, Q)

vertex_values_p_e = p_e.compute_vertex_values(mesh)
vertex_values_p_r = p.compute_vertex_values(mesh)

error = np.max(np.abs(vertex_values_p_e - vertex_values_p_r))
print('error pressure = ', error)
7 Likes
6

Hi @kamilova-maths Thanks for your code I am solving the Navier-Stokes in 2D but for cylindrical coordinate. I would appreciate if you can share variational formulation for stokes equation along with the boundary conditions you used. Thank you very much!