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…
Would you please give some advice?
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()