Second order PDE solving CFD problem
31 views (last 30 days)
Show older comments
I try to solving this PDE and this is my project
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
I run my code and Nuselt and Mean temperature and theta is "NaN",I dont know what problem in my code.
2 Comments
Torsten
on 7 Nov 2024 at 20:46
Edited: Torsten
on 7 Nov 2024 at 21:18
As already noted here:
, it does not make sense to set dtheta/dx = 0 at x = 0 and x = 2 if your equation is based on convection in x-direction.
Further note that
% 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
sets theta = 0, not dtheta/dx = 0 resp. dtheta/dr = 0.
Answers (1)
Walter Roberson
on 7 Nov 2024 at 19:35
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;
Those values creep upwards. By the time of i = 39, you are getting infinite results.
0 Comments
See Also
Categories
Find more on Computational Fluid Dynamics (CFD) in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!