I am trying to create a 2d axisymmetric solver for the navier stokes. My formulation is written in cylindrical polar coordinates. I am using the r,z, theta components to forumale the navier stokes. I have the formulation down and written below. I’m having trouble defining the trial and test functions for this problem. I’m hoping to solve for the velocity components of the z-direction(x-axis) and r-direction(yaxis) in separate variational forms. Once solved, I would like it to assembled to one u vector consisting of the z,r (x,y) components. It seems that I will need 3 test functions (ux,uy,p). and 2 trial functions w,q. Any guidance with this will be very helpful. Thank you.
This is what I have so far with the solver.
Define function space for velocity (ux,uy) and pressure (p)
W = VectorFunctionSpace(mesh, ‘P’, 2)
Define function space for system of velocities
P1 = FiniteElement(‘P’, triangle, 1)
element = MixedElement([P1, P1])
V = FunctionSpace(mesh, element)
Q = FunctionSpace(mesh, ‘P’, 1)Define boundaries
inflow = ‘near(x[0], 0)’
outflow = ‘near(x[0], 2.2)’
walls = ‘near(x[1], 0) || near(x[1], 0.41)’
cylinder = ‘on_boundary && x[0]>0.1 && x[0]<0.3 && x[1]>0.1 && x[1]<0.3’Define inflow profile
inflow_profile = (‘400.01.5x[1]*(0.41 - x[1]) / pow(0.41, 2)’, ‘0’)
Define boundary conditions
bcu_inflow = DirichletBC(V, Expression(inflow_profile, degree=2), inflow)
bcu_walls = DirichletBC(V, Constant((0, 0)), walls)
bcu_cylinder = DirichletBC(V, Constant((0, 0)), cylinder)
bcp_outflow = DirichletBC(Q, Constant(0), outflow)
bcu = [bcu_inflow, bcu_walls, bcu_cylinder]
bcp = [bcp_outflow]Define trial and test functions
w0, w1 = TestFunctions(V)
w = Function(W)
u = TrialFunction(V)
ux, uy = split(u)v = TestFunction(V)
p = TrialFunction(Q)
q = TestFunction(Q)Define functions for solutions at previous and current time steps
u_n = Function(V)
u_ = Function(V)
p_n = Function(Q)
p_ = Function(Q)ux_n, uy_n = split(u_n)
ux_, uy_ = split(u_)Define expressions used in variational forms
U = 0.5*(u_n + u)
n = FacetNormal(mesh)
f = Constant((0, 0))
k = Constant(dt)
mu = Constant(mu)
rho = Constant(rho)x,y = SpatialCoordinate(mesh)
uy_x = uy.dx(0)
uy_y = uy.dx(1)
ux_x = ux.dx(0)
ux_y = ux.dx(1)Define 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)]]))Define stress tensor
def sigma(u, p):
return 2muepsilon(u) - p*Identity(3)F1 = rho*((ux - ux_n) / k)w0dx + rhodot(dot(u_n, nabla_grad(u_n)), w)dx + inner(sigma(U, p_n), epsilon(w))dx + dot(p_nn, w)ds
- dot(munabla_grad(U)n, w)ds - muw0(1/y)ux_ydx
+ rho((uy - uy_n) / k)w1dx - muw1*(1/y)uy_ydx + muw1uy*(1/(y*y))*dx
a1 = lhs(F1)
L1 = rhs(F1)]
