# Help Varational Formulation of 2D PDE in Cylindrical Coordinates

Hi there,
I am working on the varational formulation for my advection-diffusion equations with a sink term as shwon.
My code is attached below. I have checked many times the varatrional formulation seems correct but it is not converging at all…
Thank you!

``````from fenics import *
import matplotlib.pyplot as plt
import numpy as np

# Mesh and function space
z_range, r_range = 0.25, 0.01  # domain dimensions
grid_size = 0.001 # m
z_mesh = z_range/grid_size
r_mesh = r_range/grid_size
mesh = RectangleMesh(Point(0, 0), Point(z_range, r_range), int(z_mesh), int(r_mesh))
V = FunctionSpace(mesh, 'P', 1)

def inlet(x, on_boundary):
return on_boundary and near(x[0], 0)    # x[0] = z axis
def wall(x, on_boundary):
return on_boundary and near(x[1], r_range)  # x[1] = r axis
# def central(x, on_boundary):
#    return on_boundary and near(x[1], 0)    # x[0] = z axis

def on_axis(r, on_boundary):
return on_boundary and near(r[1], 0)

# Boundary conditions
bc_T_inlet = DirichletBC(V, Constant(570), inlet)
bc_T_wall = DirichletBC(V, Constant(570), wall)
#bc_T_central = DirichletBC(V, Constant(550), central)

bc_C_inlet = DirichletBC(V, Constant(4456), inlet)

bcs_T = [bc_T_inlet, bc_T_wall]  # Temperature boundary conditions
bcs_C = [bc_C_inlet]  # Concentration boundary condition at inlet

# Constants and parameters
lambda_er = Constant(1.923)  # Effective radial conductivity, W/mK
D_er = Constant(0.0000083)  # Effective radial diffusion coefficient, m^2/s
rho_catalyst = Constant(1000)  # Catalyst density, kg/m^3
u_s = Constant(0.00049)  # Superficial velocity, m/s
DeltaH = Constant(63500)  # Reaction heat, J/kg, negative for exothermic
k0 = Constant(50000)  # Pre-exponential factor, m^3/kg s
R = Constant(8.314)  # Gas constant, J/(mol K)
E = Constant(117000)  # Activation energy, J/mol
# rho_f = Constant(8.687)
rho_f = Constant(853)
Cpf = Constant(2407)

# Solution variables
T = Function(V)
C = Function(V)
v = TestFunction(V)

T.interpolate(Constant(570))  # Intial temperature
C.interpolate(Constant(4456))  # Initial concentration

# Define the reaction term for both equations
def reaction_term(T, C):
return -k0 * exp(-E / (R * (T))) * C

# Set up equations
def setup_equations(T, C, v):
r = SpatialCoordinate(mesh)[1]

# Temperature equation
Q_MCH_in = reaction_term(T, C)

#F_T = (u_s * rho_f * Cpf * T.dx(1) * v - lambda_er * (T.dx(0) * v.dx(0) + T.dx(1) * v.dx(1)) - lambda_er * (1/r) * T.dx(0) * v - DeltaH * rho_catalyst * Q_MCH_in * v) * dx
# Concentration equation
#F_C = (u_s * C.dx(1) * v - D_er * (C.dx(0) * v.dx(0) + C.dx(1) * v.dx(1)) - D_er * (1/r) * C.dx(0) * v - rho_catalyst * Q_MCH_in * v) * dx

F_T = (r * u_s * rho_f * Cpf * T.dx(0) * v - lambda_er * (T.dx(0) * v.dx(0) * r + T.dx(1) * v.dx(1) * r + 1/r * T.dx(1) * v * r) - DeltaH * rho_catalyst * Q_MCH_in * v * r) * dx

F_C = (r * u_s * C.dx(0) * v - D_er * (C.dx(0) * v.dx(0) * r + C.dx(1) * v.dx(1) * r + 1/r * C.dx(1) * v * r) - r * rho_catalyst * Q_MCH_in * v) * dx

return F_T, F_C

# Solve equations iteratively
tol = 1E-10  # Convergence tolerance
error = 1.0
iter_count = 0
max_iter = 100

while error > tol and iter_count < max_iter:
T_prev = T.copy(deepcopy=True)
C_prev = C.copy(deepcopy=True)

# Solve temperature equation
F_T, F_C = setup_equations(T, C, v)
solve(F_T == 0, T, bcs_T)

# Solve concentration equation
solve(F_C == 0, C, bcs_C)

# Error and convergence check
error = errornorm(T, T_prev, 'L2') + errornorm(C, C_prev, 'L2')
iter_count += 1
print(f"Iteration {iter_count}: Error = {error}")

boundary_coords = mesh.coordinates()  # mesh coordinates
boundary_values = T.compute_vertex_values(mesh)
for i, x in enumerate(boundary_coords):
if near(x[0], 0):  # check if z = 0
print(f"Temperature at inlet ({x[0]}, {x[1]}) = {boundary_values[i]} K")
if near(x[1], r_range):  # check if r = r_range
print(f"Temperature at wall ({x[0]}, {x[1]}) = {boundary_values[i]} K")
if near(x[1], 0):  # check if r = 0
print(f"Temperature at centreline ({x[0]}, {x[1]}) = {boundary_values[i]} K")

# Plot solution
plt.figure(figsize=(10, 5))
p = plot(T)
plt.colorbar(p)
plt.title('Temperature Distribution')
plt.show()

# Plot solution
plt.figure(figsize=(10, 5))
p = plot(C)
plt.colorbar(p)
plt.title('Concentration Distribution')
plt.show()

# Define the line for plotting (r = 0 and r = r_range / 2)
r_0 = 0
r_mid = r_range / 2
r_r = r_range
r_range_points = [Point(z, r_0) for z in np.linspace(0, z_range, 100)]
r_mid_points = [Point(z, r_mid) for z in np.linspace(0, z_range, 100)]
r_r_points = [Point(z, r_r) for z in np.linspace(0, z_range, 100)]

# Temperatures and concentrations at r = 0 and r = r_range / 2
temperatures_r_0 = [T(point) for point in r_range_points]
concentrations_r_0 = [C(point) for point in r_range_points]
temperatures_r_mid = [T(point) for point in r_mid_points]
concentrations_r_mid = [C(point) for point in r_mid_points]
temperatures_r_r = [T(point) for point in r_r_points]
concentrations_r_r = [C(point) for point in r_r_points]

# Plotting temperature
plt.figure(figsize=(10, 5))
plt.plot(np.linspace(0, z_range, 100), temperatures_r_0, label='Temperature at r=0')
plt.plot(np.linspace(0, z_range, 100), temperatures_r_mid, label='Temperature at r=r/2')
plt.plot(np.linspace(0, z_range, 100), temperatures_r_r, label='Temperature at r=r')
plt.xlabel('Axial position z (m)')
plt.ylabel('Temperature (K)')
plt.title('Temperature Profile')
plt.legend()
plt.grid(True)
plt.show()

# Plotting concentration
plt.figure(figsize=(10, 5))
plt.plot(np.linspace(0, z_range, 100), concentrations_r_0, label='Concentration at r=0')
plt.plot(np.linspace(0, z_range, 100), concentrations_r_mid, label='Concentration at r=r/2')
plt.plot(np.linspace(0, z_range, 100), concentrations_r_r, label='Concentration at r=r')
plt.xlabel('Axial position z (m)')
plt.ylabel('Concentration (mol/m^3)')
plt.title('Concentration Profile')
plt.legend()
plt.grid(True)
plt.show()

``````