second order center differencing scheme soluting to PDE

I am trying to find a solution to this PDE:

u* (∂θ/∂x) = (1-a^2)[(1/r) (∂/∂r) (r *(∂θ/∂r))+(∂^2θ/∂x^2]+γ*(∂u/∂r)^2 ,0≤ x ≤ 2, a≤ r ≤ 1

u=2(1-EU*)(1-r^2+B*ln r)/M, B=(a^2-1)/ln r ,E=(a^2-0.5B)/(a^2-1) ,B*=B(1-0.5U*)/(1-EU*), M=1+a^2-B

for r=a

∂θ/∂r=0 ,0≤ x <1

∂θ/∂r=cos[4π(x-1)] ,1≤ x ≤ 1.5

∂θ/∂r=0 ,1.5< x ≤ 2

for r=1

∂θ/∂r=0 ,0≤ x <1

θ=sin[4π(x-1)] ,1≤ x ≤ 1.5

∂θ/∂r=0 ,1.5≤ x ≤<2

for x=0 ,or 2

∂θ/∂r=0

a=0.1,U*=-1,γ=0.5

grid sizes:51*101

I wrote the following code to solve this PDE, but the result is incorrect. I feel like my problem-solving approach has stalled. Could someone please help?

Thank you!

code:

% MATLAB code for solving a partial differential equation using the finite volume method,

% and solving with SOR iteration and TDMA.

% Basic parameters and grid setup

a = 0.1;

gamma = 0.5;

U_star = -1;

Nx = 101; % Number of grid points in x-direction

Nr = 51; % Number of grid points in r-direction

dx = 2 / (Nx - 1); % Grid spacing in x-direction

dr = (1 - a) / (Nr - 1); % Grid spacing in r-direction

x = linspace(0, 2, Nx);

r = linspace(a, 1, Nr);

% Initialize variables

theta = zeros(Nr, Nx);

B = (a^2 - 1) / log(r(end));

E = (a^2 - 0.5 * B) / (a^2 - 1);

B_star = B * (1 - 0.5 * U_star) / (1 - E * U_star);

M = 1 + a^2 - B;

% Set SOR parameters

omega = 1.5; % Over-relaxation factor

tol = 1e-6; % Convergence criterion

max_iter = 10000; % Maximum number of iterations

% Initialize residual and iteration count

residual = inf;

iter = 0;

% Iterative solution

display('Starting SOR iteration...');

while residual > tol && iter < max_iter

theta_old = theta;

% Update in horizontal direction using SOR

for i = 2:Nr-1

for j = 2:Nx-1

% Calculate partial derivatives in r and x directions

dtheta_dr_r = (theta(i+1, j) - theta(i-1, j)) / (2 * dr);

dtheta_dx2 = (theta(i, j+1) - 2 * theta(i, j) + theta(i, j-1)) / (dx^2);

% Discretize the governing equation

theta_new = (1-a^2) * (1/r(i) * dtheta_dr_r + dtheta_dx2) + gamma * (U_star / dr)^2;

% Update using SOR formula

theta(i, j) = (1 - omega) * theta(i, j) + omega * theta_new;

end

end

% Calculate residual

residual = max(max(abs(theta - theta_old)));

iter = iter + 1;

end

display(['SOR iteration completed, total iterations: ', num2str(iter), '.']);

% Use TDMA to solve the matrix problem

for j = 2:Nx-1

% TDMA preparation

A = zeros(Nr, Nr);

B = zeros(Nr, 1);

for i = 2:Nr-1

A(i, i-1) = -1 / dr^2;

A(i, i) = 2 / dr^2 + 2 / dx^2;

A(i, i+1) = -1 / dr^2;

B(i) = theta(i, j);

end

% Set boundary conditions

if x(j) >= 1 && x(j) <= 1.5

B(1) = cos(4 * pi * (x(j) - 1)); % Boundary condition dr=cos[4*pi*(x-1)]

else

B(1) = 0; % Boundary condition dr=0

end

A(1, 1) = 1;

A(Nr, Nr) = 1;

B(Nr) = 0; % Boundary condition dr=0

% Solve using TDMA

theta(:, j) = A\B;

end

% Plot the results

figure;

surf(x, r, theta);

xlabel('x');

ylabel('r');

zlabel('theta');

title('Numerical solution of theta');

% Handling boundary conditions

for j = 1:Nx

if x(j) >= 1 && x(j) <= 1.5

theta(1, j) = -sin(4 * pi * (x(j) - 1));

end

end

% Final result

figure;

contourf(x, r, theta);

colorbar;

xlabel('x');

ylabel('r');

title('Contour plot of theta');

[Torsten]

Your boundary conditions for x = 0 or x = 2 cannot be prescribed as dtheta/dr = 0.

[Xu]

I see now that my initial input was incorrect. The correct boundary condition should be for x=0 or 2, with dθ/dx=0

[Torsten]

Since convection is in x-direction, I find it strange that you don't prescribe theta at x = 0 (or x = 2, depending on the sign of u*).

What do you model with your PDE ?

Answer 1

[Abhimenyu]

Hi Xu,

I understand that the goal of the MATLAB code shared by you is to solve a partial differential equation (PDE) using the finite volume method (FVM) with the Successive Over-Relaxation (SOR) technique.

Upon reviewing the original code, I have identified a few areas that require correction and improvement.

Please find the revised MATLAB code below with explanation on the changes made:

% MATLAB code for solving a partial differential equation using the finite volume method,
% and solving with SOR iteration and TDMA.
% Basic parameters and grid setup
a = 0.1;
gamma = 0.5;
U_star = -1;
Nx = 101; % Number of grid points in x-direction
Nr = 51; % Number of grid points in r-direction
dx = 2 / (Nx - 1); % Grid spacing in x-direction
dr = (1 - a) / (Nr - 1); % Grid spacing in r-direction
x = linspace(0, 2, Nx);
r = linspace(a, 1, Nr);

% Initialize variables
theta = zeros(Nr, Nx);
[X, R] = meshgrid(x, r);
B = (a^2 - 1) / log(a);
E = (a^2 - 0.5 * B) / (a^2 - 1);
B_star = B * (1 - 0.5 * U_star) / (1 - E * U_star);
M = 1 + a^2 - B;

• Velocity Profile Issue and Initialization: The velocity field 'u' was missing, which is crucial for the convective term in the PDE. Also the initialisation at r=1 was required.

% Calculate velocity profile
u = 2 * (1 - E * U_star) * (1 - R.^2 + B_star * log(R)) / M;

% Set SOR parameters
omega = 1.5; % Over-relaxation factor
tol = 1e-6; % Convergence criterion
max_iter = 10000; % Maximum number of iterations

% Initialize residual and iteration count
residual = inf;
iter = 0;

% Set initial boundary conditions at r = 1
for j = 1:Nx
    if x(j) >= 1 && x(j) <= 1.5
        theta(end, j) = sin(4 * pi * (x(j) - 1));
    end
end

% Iterative solution

• Discretization Error: In the code shared, there was incomplete treatment of cylindrical coordinate and missing convective term (u*∂θ/∂x).
• Boundary condition: The code needed implementation of gradient boundary conditions, and inclusion of missing proper treatment at x = 0 and x = 2.

display('Starting SOR iteration...');

Starting SOR iteration...

while residual > tol && iter < max_iter
    theta_old = theta;
    
    % Update interior points using SOR
    for i = 2:Nr-1
        for j = 2:Nx-1
            % Calculate coefficients
            ae = (1-a^2)/(dx^2);
            aw = (1-a^2)/(dx^2);
            an = (1-a^2)/(dr^2) + (1-a^2)/(2*r(i)*dr);
            as = (1-a^2)/(dr^2) - (1-a^2)/(2*r(i)*dr);
            ap = -(2*(1-a^2)/(dx^2) + 2*(1-a^2)/(dr^2));
            
            % Source term with convection
            Su = gamma * (u(i,j)/dr)^2 - u(i,j)*(theta(i,j+1) - theta(i,j-1))/(2*dx);
            
            % Calculate new value
            theta_new = -(ae*theta(i,j+1) + aw*theta(i,j-1) + an*theta(i+1,j) + as*theta(i-1,j) + Su)/ap;
            
            % Apply SOR update
            theta(i,j) = (1 - omega) * theta(i,j) + omega * theta_new;
        end
    end
    
    % Apply boundary conditions at r = a
    for j = 1:Nx
        if x(j) >= 1 && x(j) <= 1.5
            theta(1,j) = theta(2,j) - dr*cos(4*pi*(x(j)-1));
        else
            theta(1,j) = theta(2,j);  % Zero gradient
        end
    end
    
    % Apply boundary conditions at x = 0 and x = 2
    theta(:,1) = theta(:,2);
    theta(:,end) = theta(:,end-1);
    
    % Calculate residual
    residual = max(max(abs(theta - theta_old)));
    iter = iter + 1;
end
display(['SOR iteration completed, total iterations: ', num2str(iter), '.']);

SOR iteration completed, total iterations: 4656.

• TDMA improvements: The source term (b vector) only contained the previous iteration value and the boundary conditions were not properly incorporated into the system. Now the code switches between different boundary types based on x position.

% Use TDMA for final solution refinement
for j = 2:Nx-1
    % TDMA preparation
    A = zeros(Nr, Nr);
    b = zeros(Nr, 1);
    
    % Set up tridiagonal system
    for i = 2:Nr-1
        A(i,i-1) = as;
        A(i,i) = ap;
        A(i,i+1) = an;
        b(i) = -(ae*theta(i,j+1) + aw*theta(i,j-1) + Su);
    end
    
    % Boundary conditions
    A(1,1) = 1; A(1,2) = -1;  % Zero gradient or specified at r = a
    A(Nr,Nr-1) = -1; A(Nr,Nr) = 1;  % Zero gradient at r = 1
    
    if x(j) >= 1 && x(j) <= 1.5
        b(1) = -dr*cos(4*pi*(x(j)-1));  % Specified gradient
        b(Nr) = sin(4*pi*(x(j)-1));     % Specified value
    else
        b(1) = 0;  % Zero gradient
        b(Nr) = 0; % Zero gradient
    end
    
    % Solve system
    theta(:,j) = A\b;
end

% Plot results
figure;
surf(x, r, theta);
xlabel('x');
ylabel('r');
zlabel('theta');
title('Numerical solution of theta');

% Final result
figure;
contourf(x, r, theta);
colorbar;
xlabel('x');
ylabel('r');
title('Contour plot of theta');

I hope this helps!

[Xu]

Thank you very much for reviewing and adjusting this code! With your help, it runs more smoothly, and I’ve gained a better understanding of the details. I really appreciate your time and assistance!