Second order PDE solving CFD problem

29 Ansichten (letzte 30 Tage)
BoYun
BoYun am 7 Nov. 2024 um 16:06
Kommentiert: BoYun am 11 Nov. 2024 um 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 Kommentare
Torsten
Torsten am 7 Nov. 2024 um 20:46
Bearbeitet: Torsten am 7 Nov. 2024 um 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 am 11 Nov. 2024 um 7:19
Ok，thank for yoyr help!

Melden Sie sich an, um zu kommentieren.

Melden Sie sich an, um diese Frage zu beantworten.

Antworten (1)

Walter Roberson
Walter Roberson am 7 Nov. 2024 um 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.

Kategorien

Mehr zu Computational Fluid Dynamics (CFD) finden Sie in Help Center und File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!

Translated by