Second order PDE solving CFD problem

Question

BoYun 2024-11-7，16:06

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2164680-second-order-pde-solving-cfd-problem

评论： BoYun 2024-11-11，7:19

在 MATLAB Online 中打开

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

% 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

% 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

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

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

I run my code and Nuselt and Mean temperature and theta is "NaN"，I dont know what problem in my code.

2 个评论
显示无隐藏无

Torsten 2024-11-7，20:46

编辑：Torsten 2024-11-7，21:18

在 MATLAB Online 中打开

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 2024-11-11，7:19

Ok，thank for yoyr help!

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Walter Roberson 2024-11-7，19:35

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2164680-second-order-pde-solving-cfd-problem#answer_1542555

在 MATLAB Online 中打开

            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 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

Second order PDE solving CFD problem

2 个评论
显示无隐藏无

回答（1 个）

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

另请参阅

类别

标签

Community Treasure Hunt

Second order PDE solving CFD problem

2 个评论 显示 无隐藏 无

回答（1 个）

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

另请参阅

类别

标签

Community Treasure Hunt

2 个评论
显示无隐藏无

0 个评论
显示 -2更早的评论隐藏 -2更早的评论