Hi guys,
I am working on solving a 2D energy, mass balance equation with a arrenius equation in a cylindrical coordiates.
I believe when r = 0, z is the centreline and r = radius, z is the wall. However, in my case, the temperature is fixed at wall and inlet, and the centreline is neumann BC. The equation has a sink term so the temperature should be lower closing to the centreline.
My results gives strange result as shown below, the inlet closing to the centreline shows higher temperature compare with r = radius/2, which does not make sense. I have checked all the code and BC, there is no problem. If the centreline temperature is fixed, the result makes sense.
It would be appreciated if anyone can help or discuss!
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.0005 # 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
# 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(4000) # Catalyst density, kg/m^3
u_s = Constant(0.00049) # Superficial velocity, m/s
DeltaH = Constant(65000) # 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(704)
Cpf = Constant(3037)
# 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] # r coordinate
# 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
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")
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")
# 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()