I trying to solving this PDE:
but my code output Nuselt number is "NaN"，I dont know what problem in this
And this is my Matlab code:
% CFD Project - Finite Volume Method for Temperature Field in Concentric Annuli (Using Point SOR and TDMA)
clear;
clc;
% Problem parameters
a_values = [0.1, 0.5, 0.8];
U_prime_values = [-1, -2];
gamma_values = [0, 0.5];
Nr = 51; % Radial grid points
Nx = 101; % Axial grid points
Lx = 2; % Axial length
R_inner = 0.1; % Inner radius
R_outer = 1; % Outer radius
omega = 1.5; % SOR relaxation factor
% Generate radial and axial grids
r = linspace(R_inner, R_outer, Nr); % Radial grid points
x = linspace(0, Lx, Nx); % Axial grid points
dr = r(2) - r(1); % Radial step size
dx = x(2) - x(1); % Axial step size
% Check if dx and dr are zero
if dx == 0 || dr == 0
error('dx or dr is zero. Please check the grid generation code.');
end
% Initial conditions and boundary conditions
theta = rand(Nr, Nx) * 1e-6; % Initialize theta with very small random values
% Define boundary conditions
for i = 1:Nr
theta(i, 1) = 0; % Neumann boundary at x = 0
theta(i, Nx) = 0; % Neumann boundary at x = Lx
end
for j = 1:Nx
if x(j) >= 1 && x(j) <= 1.5
theta(1, j) = cos(4 * pi * (x(j) - 1)); % Inner boundary at r = R_inner
else
theta(1, j) = 0;
end
if x(j) >= 1 && x(j) <= 1.5
theta(Nr, j) = -sin(4 * pi * (x(j) - 1)); % Outer boundary at r = R_outer
else
theta(Nr, j) = 0;
end
end
% Set parameters for finite volume discretization
scheme = 'FOU';
tol = 1e-8; % Convergence tolerance
max_iter = 10000; % Maximum iterations
% SOR Iteration
for t = 1:max_iter
theta_old = theta; % Save the previous iteration result
for i = 2:Nr-1
for j = 2:Nx-1
[conv_term, diff_term] = computeTerms(theta, i, j, dr, dx, scheme);
theta_new = (1 - omega) * theta(i, j) + omega * 0.5 * (conv_term + diff_term);
theta(i, j) = theta_new;
end
end
% Check for convergence
if max(max(abs(theta - theta_old))) < tol
fprintf('SOR iteration converged in %d steps\n', t);
break;
end
% Display intermediate results
fprintf('Iteration %d: max theta = %f, min theta = %f\n', t, max(theta(:)), min(theta(:)));
end
% Post-processing: Calculate the mean temperature and Nusselt number at the outer cylinder
theta_m = computeMeanTemperature(theta, r, dr);
Nu_a = computeNusseltNumber(theta, r, dr, R_inner, R_outer, theta_m);
% Display results
fprintf('Mean temperature: %f\n', theta_m);
fprintf('Nusselt number at the outer cylinder: %f\n', Nu_a);
% Plot results
figure;
imagesc(x, r, theta);
colorbar;
title('Temperature distribution \theta');
xlabel('x');
ylabel('r');
set(gca, 'YDir', 'normal');
figure;
plot(r, mean(theta, 2), '-o');
title('Radial mean temperature distribution');
xlabel('Radial position r');
ylabel('Mean temperature \theta_m');
% --- Function definitions ---
function [conv_term, diff_term] = computeTerms(theta, i, j, dr, dx, scheme)
% Compute convection and diffusion terms based on the selected scheme
switch scheme
case 'FOU'
conv_term = (theta(i, j) - theta(i, j-1)) / dx;
diff_term = (theta(i+1, j) - 2 * theta(i, j) + theta(i-1, j)) / dr^2;
case 'SOCD'
conv_term = (theta(i, j+1) - theta(i, j-1)) / (2 * dx);
diff_term = (theta(i+1, j) - 2 * theta(i, j) + theta(i-1, j)) / dr^2;
end
end
function theta_m = computeMeanTemperature(theta, r, dr)
% Compute the mean temperature using the integral definition
r_matrix = repmat(r', 1, size(theta, 2));
integral_num = sum(sum(r_matrix .* theta)) * dr;
integral_den = sum(r) * dr;
% Display intermediate calculation results
fprintf('Integral numerator (integral_num): %f\n', integral_num);
fprintf('Integral denominator (integral_den): %f\n', integral_den);
% Avoid division by zero
if integral_den == 0
theta_m = NaN;
else
theta_m = integral_num / integral_den;
end
end
function Nu_a = computeNusseltNumber(theta, r, dr, R_inner, R_outer, theta_m)
% Compute the Nusselt number at the outer cylinder
dtheta_dr_outer = (theta(end, :) - theta(end-1, :)) / dr;
% Display intermediate calculation results
fprintf('dtheta/dr at outer (dtheta_dr_outer): %f\n', dtheta_dr_outer);
fprintf('theta(end, :) - theta_m: %f\n', theta(end, :) - theta_m);
% Avoid division by zero
if all(theta(end, :) == theta_m)
Nu_a = NaN;
else
Nu_a = -2 * (1 - R_inner) * dtheta_dr_outer / (theta(end, :) - theta_m);
end
end
