I am developing a 2d axisymmetric solver using the following formulation.
\begin{align*}
& Continuity \xrightarrow{} \int_\Omega m\frac{\partial u}{\partial x}d\Omega + \int_\Omega m\frac{\partial v}{\partial y}d\Omega + \int_\Omega m\frac{v}{y}~d\Omega = 0 \ \
& z-component(x) \xrightarrow{} \int_\Omega w\rho\frac{\partial u}{\partial t}~d\Omega + \int_\Omega w\rho(u \frac{\partial u}{\partial x})~d\Omega + \int_\Omega w\rho(v \frac{\partial u}{\partial y})~d\Omega - \int_\Omega p(\frac{\partial w}{\partial x})~d\Omega\
&~~~~~~~~~~~~~~~~~~~~~~~~ - \int_\Gamma wt~d\Gamma - \mu \int_\Omega w(\frac{1}{y} \frac{\partial u}{\partial y})~d\Omega = 0 \ \
& r-component(y) \xrightarrow{} \int_\Omega w\rho\frac{\partial v}{\partial t}~d\Omega + \int_\Omega w\rho(u \frac{\partial v}{\partial x})~d\Omega + \int_\Omega w(v \frac{\partial v}{\partial y})~d\Omega - \int_\Omega p(\frac{\partial w}{\partial y})~d\Omega - \int_\Gamma wt~d\Gamma \
&~~~~~~~~~~~~~~~~~~~~~~~~ - \mu\int_\Omega w(\frac{1}{y} \frac{\partial v}{\partial y})~d\Omega + \int_\Omega w\mu (\frac{v}{y^2})~d\Omega = 0 \
\end{align*}
My first attempt was to modify an existing Navier Stokes solver to solve for u_x,u_y separetely. This is what I have so far
from future import print_function
from fenics import *
from mshr import *
import numpy as np
T = 3.0 # final time
num_steps = 5000 # number of time steps
dt = T / num_steps # time step size
mu = 0.001 # dynamic viscosity
rho = 1 # density
Create mesh
channel = Rectangle(Point(0, 0), Point(2.2, 0.41))
cylinder = Circle(Point(0.2, 0.2), 0.05)
domain = channel - cylinder
mesh = generate_mesh(domain, 64)
V = VectorFunctionSpace(mesh, ‘P’, 2)
Q = FunctionSpace(mesh, ‘P’, 1)
Define function space for velocity
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)
# U2 = 0.5*(uy_n + uy)
n = FacetNormal(mesh)
f = Constant((0, 0))
k = Constant(dt)
mu = Constant(mu)
rho = Constant(rho)
x,y = SpatialCoordinate(mesh)
ux = u[0]
uy = u[1]
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)
Define variational problem for step 1
F1 = rho*dot((u - u_n) / k, v)*dx \
+ rho*dot(dot(u_n, nabla_grad(u_n)), v)*dx \
+ inner(sigma(U, p_n), epsilon(v))*dx \
+ dot(p_n*n, v)ds - dot(munabla_grad(U)*n, v)*ds \
- dot(f, v)*dx
a1 = lhs(F1)
L1 = rhs(F1)
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)
Define variational problem for step 2
a2 = dot(nabla_grad§, nabla_grad(q))*dx
L2 = dot(nabla_grad(p_n), nabla_grad(q))*dx - (1/k)div(u_)qdx - (1/k)muuy_(1/y)qdx
# # Define variational problem for step 3
a3 = dot(u, w)dx + (1/k)muu(1/y)wdx
L3 = dot(u_, w)dx - kdot(nabla_grad(p_ - p_n),w)*dx
Assemble matrices
A1 = assemble(a1)
A2 = assemble(a2)
A3 = assemble(a3)
Apply boundary conditions to matrices
[bc.apply(A1) for bc in bcu]
[bc.apply(A2) for bc in bcp]
Create XDMF files for visualization output
xdmffile_u = XDMFFile(‘navier_stokes_cylinder/velocity.xdmf’)
xdmffile_p = XDMFFile(‘navier_stokes_cylinder/pressure.xdmf’)
Create time series (for use in reaction_system.py)
timeseries_u = TimeSeries(‘navier_stokes_cylinder/velocity_series’)
timeseries_p = TimeSeries(‘navier_stokes_cylinder/pressure_series’)
Save mesh to file (for use in reaction_system.py)
File(‘navier_stokes_cylinder/cylinder.xml.gz’) << mesh
Create progress bar
progress = Progress(‘Time-stepping’)
progress = Progress(‘Time-stepping’, num_steps)
set_log_level(PROGRESS)
Time-stepping
t = 0
count = 0
for n in range(num_steps):
# Update current time
t += dt
# Step 1: Tentative velocity step
b1 = assemble(L1)
[bc.apply(b1) for bc in bcu]
solve(A1, u_.vector(), b1, 'bicgstab', 'hypre_amg')
# Step 2: Pressure correction step
b2 = assemble(L2)
[bc.apply(b2) for bc in bcp]
solve(A2, p_.vector(), b2, 'bicgstab', 'hypre_amg')
# # Step 3: Velocity correction step
# b3 = assemble(L3)
# solve(A3, u_.vector(), b3, 'cg', 'sor')
# Plot solution
# plot(u_, title='Velocity')
# plot(p_, title='Pressure')
# Save solution to file (XDMF/HDF5)
N = 50 #save every N time steps
if (count/N).is_integer():
xdmffile_u.write(u_, t)
xdmffile_p.write(p_, t)
# Save nodal values to file
timeseries_u.store(u_.vector(), t)
timeseries_p.store(p_.vector(), t)
# Update previous solution
u_n.assign(u_)
p_n.assign(p_)
# Update progress bar
# progress.update(t / T)
set_log_level(LogLevel.PROGRESS)
progress += 1
set_log_level(LogLevel.ERROR)
print('u max:', u_.vector().get_local().max())
count = count + 1