Second order PDE solving CFD problem

29 visualizaciones (últimos 30 días)
BoYun
BoYun el 7 de Nov. de 2024 a las 16:06
Comentada: BoYun el 11 de Nov. de 2024 a las 7:19
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
Iteration 1: max theta = Inf, min theta = -Inf Iteration 2: max theta = Inf, min theta = -Inf Iteration 3: max theta = Inf, min theta = -Inf Iteration 4: max theta = Inf, min theta = -Inf Iteration 5: max theta = Inf, min theta = -Inf Iteration 6: max theta = Inf, min theta = -Inf Iteration 7: max theta = Inf, min theta = -Inf Iteration 8: max theta = Inf, min theta = -Inf Iteration 9: max theta = Inf, min theta = -Inf Iteration 10: max theta = Inf, min theta = -Inf Iteration 11: max theta = Inf, min theta = -Inf Iteration 12: max theta = Inf, min theta = -Inf Iteration 13: max theta = Inf, min theta = -Inf Iteration 14: max theta = Inf, min theta = -Inf Iteration 15: max theta = Inf, min theta = -Inf Iteration 16: max theta = Inf, min theta = -Inf Iteration 17: max theta = Inf, min theta = -Inf Iteration 18: max theta = Inf, min theta = -Inf Iteration 19: max theta = Inf, min theta = -Inf Iteration 20: max theta = Inf, min theta = -Inf Iteration 21: max theta = Inf, min theta = -Inf Iteration 22: max theta = Inf, min theta = -Inf Iteration 23: max theta = Inf, min theta = -Inf Iteration 24: max theta = Inf, min theta = -Inf Iteration 25: max theta = Inf, min theta = -Inf Iteration 26: max theta = Inf, min theta = -Inf Iteration 27: max theta = Inf, min theta = -Inf Iteration 28: max theta = Inf, min theta = -Inf Iteration 29: max theta = Inf, min theta = -Inf Iteration 30: max theta = Inf, min theta = -Inf Iteration 31: max theta = Inf, min theta = -Inf Iteration 32: max theta = Inf, min theta = -Inf Iteration 33: max theta = Inf, min theta = -Inf Iteration 34: max theta = Inf, min theta = -Inf Iteration 35: max theta = Inf, min theta = -Inf Iteration 36: max theta = Inf, min theta = -Inf Iteration 37: max theta = Inf, min theta = -Inf Iteration 38: max theta = Inf, min theta = -Inf Iteration 39: max theta = Inf, min theta = -Inf Iteration 40: max theta = Inf, min theta = -Inf Iteration 41: max theta = Inf, min theta = -Inf Iteration 42: max theta = Inf, min theta = -Inf Iteration 43: max theta = Inf, min theta = -Inf Iteration 44: max theta = Inf, min theta = -Inf
SOR iteration converged in 45 steps
% Post-processing: Calculate the mean temperature and Nusselt number at the outer cylinder
theta_m = computeMeanTemperature(theta, r, dr);
Integral numerator (integral_num): NaN Integral denominator (integral_den): 0.504900
Nu_a = computeNusseltNumber(theta, r, dr, R_inner, R_outer, theta_m);
dtheta/dr at outer (dtheta_dr_outer): 0.000000 dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): NaN dtheta/dr at outer (dtheta_dr_outer): 0.000000 theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN theta(end, :) - theta_m: NaN
% Display results
fprintf('Mean temperature: %f\n', theta_m);
Mean temperature: NaN
fprintf('Nusselt number at the outer cylinder: %f\n', Nu_a);
Nusselt number at the outer cylinder: NaN
% 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 comentarios
Torsten
Torsten el 7 de Nov. de 2024 a las 20:46
Editada: Torsten el 7 de Nov. de 2024 a las 21:18
As already noted here:
https://uk.mathworks.com/matlabcentral/answers/2164340-second-order-center-differencing-scheme-soluting-to-pde?s_tid=srchtitle
, 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.
BoYun
BoYun el 11 de Nov. de 2024 a las 7:19
Ok，thank for yoyr help!

Iniciar sesión para comentar.

Iniciar sesión para responder a esta pregunta.

Respuestas (1)

Walter Roberson
Walter Roberson el 7 de Nov. de 2024 a las 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.

Categorías

Más información sobre Computational Fluid Dynamics (CFD) en Help Center y File Exchange.

Etiquetas

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by