Coding the Thomas algorithm for a heat and mass transfer problem

I am trying to solve, vary t and y, and plot for temperature, concentration and velocity in a heat and mass transfer problem. I have two codes that I am trying to merge into one. When I run the function, I get:

Not enough input arguments.

Error in solvePDEs (line 27)

y(n) = linspace(y_span(1), y_span(2), num_y);

The function is:

function [u, theta, phi] = solvePDEs(Gr, Re, Pr, Sc, k, N_c, N_t, R, n, t_span, y_span, num_y, num_t)

% Solves the system of partial differential equations (PDEs) describing

% natural convection heat and mass transfer in a porous medium.

%

% Args:

% Gr: Grashof number.

% Re: Reynolds number.

% Pr: Prandtl number.

% Sc: Schmidt number.

% k: Thermophoretic parameter.

% N_c: Concentration buoyancy ratio.

% N_t: Thermal buoyancy ratio.

% R: Reaction rate constant.

% n: Order of reaction.

% t_span: Time span for simulation.

% y_span: Spatial domain for simulation.

% num_y: Number of grid points in y-direction.

% num_t: Number of time steps.

%

% Returns:

% u: Velocity field.

% theta: Temperature field.

% phi: Concentration field.

i=1:10;

for n=i

y(n) = linspace(y_span(1), y_span(2), num_y);

dy = y(2) - y(1);

end

% Discretize the time domain

t = linspace(t_span(1), t_span(2), num_t);

dt = t(2) - t(1);

% Discretize the spatial domain

% Initialize the solution matrices

u = zeros(num_y, num_t);

theta = zeros(num_y, num_t);

phi = zeros(num_y, num_t);

% Apply initial conditions

u(:, 1) = 0;

theta(:, 1) = 0;

phi(:, 1) = 0;

% Apply boundary conditions

u(1, :) = 1;

theta(1, :) = 1/2;

phi(1, :) = 1/2;

u(end, :) = 0;

theta(end, :) = -1/2;

phi(end, :) = -1/2;

% Implement the finite difference method

for i = 2:num_t

for j = 2:num_y-1

% Calculate the second-order central difference approximations

d2u_dy2 = (u(j+1, i-1) - 2*u(j, i-1) + u(j-1, i-1)) / dy^2;

d2theta_dy2 = (theta(j+1, i-1) - 2*theta(j, i-1) + theta(j-1, i-1)) / dy^2;

d2phi_dy2 = (phi(j+1, i-1) - 2*phi(j, i-1) + phi(j-1, i-1)) / dy^2;

% Calculate the first-order forward difference approximation

dtheta_dy = (theta(j+1, i-1) - theta(j, i-1)) / dy;

% Update the solution using the finite difference equations

u(j, i) = u(j, i-1) + dt * (d2u_dy2 + (Gr/Re) * (theta(j, i-1) + b*phi(j, i-1)) - alpha);

theta(j, i) = theta(j, i-1) + dt/Pr * d2theta_dy2;

phi(j, i) = phi(j, i-1) + dt/Sc * (d2phi_dy2 + k * (phi(j, i-1) + N_c)/(theta(j, i-1) + N_t) * dtheta_dy - R*phi(j, i-1)^n);

end

end

% Return the solution matrices

return;

end

AND

function solvePDEs()

% Parameters (example values, these need to be defined or adjusted)

Gr = 1; Re = 100; alpha = 0.1; b = a0.5;

Pr = 0.7; Sc = 1; k = 0.1; N_c = 0.1; N_t = 0.1; R = 0.01; n = 2;

% Spatial and temporal discretization

Ny = 50; % Number of spatial points

y = linspace(0, 1, Ny);

dy = y(2) - y(1);

dt = 0.01; % Time step

T = 1; % Total time

Nt = floor(T/dt);

% Initial conditions

u = zeros(Ny, 1);

theta = zeros(Ny, 1);

phi = zeros(Ny, 1);

% Preallocate for speed

u_new = u;

theta_new = theta;

phi_new = phi;

% Time-stepping loop

for t = 1:Nt

phi_term = zeros(1,100/Ny+1)% preallocated to the correct size

% Solve for u

u_new = thomasAlgorithm(Ny, dy, dt, u, (Gr/Re) * (theta + b*phi) - alpha);

% Solve for theta

theta_new = thomasAlgorithm(Ny, dy, dt/Pr, theta, zeros(Ny, 1));

% Solve for phi

dtheta_dy = diff(theta) / dy;

dtheta_dy = [dtheta_dy; dtheta_dy(end)]; % Approximate derivative at the boundary

phi_term = k * diff((phi + N_c) ./ (theta + N_t) .* dtheta_dy) / dy;

phi_term = [0; phi_term; 0]; % Boundary conditions

phi_new = thomasAlgorithm(Ny, dy, dt/Sc, phi, phi_term - R * phi.^n);

% Update solutions

u = u_new;

theta = theta_new;

phi = phi_new;

% Plotting at selected time steps

if mod(t, floor(Nt/4)) == 0 || t == 1

figure;

subplot(3, 1, 1);

plot(y, u, '-o');

title(['Velocity u at t = ', num2str(t*dt)]);

subplot(3, 1, 2);

plot(y, theta, '-o');

title(['Temperature \theta at t = ', num2str(t*dt)]);

subplot(3, 1, 3);

plot(y, phi, '-o');

title(['Concentration \phi at t = ', num2str(t*dt)]);

end

end

end

function x_new = thomasAlgorithm(N, dy, dt, x, source)

% Coefficients for the tridiagonal matrix

a = repmat(-dt/dy^2, N, 1);

b = repmat(1 + 2*dt/dy^2, N, 1);

c = repmat(-dt/dy^2, N, 1);

d = x + dt * source;

% Boundary conditions

b(1) = 1; b(N) = 1;

a(1) = 0; c(N) = 0;

d(1) = 1; d(N) = 0; % Example boundary values

% Thomas algorithm

for i = 2:N

m = a(i) / b(i-1);

b(i) = b(i) - m * c(i-1);

d(i) = d(i) - m * d(i-1);

end

x_new = zeros(N, 1);

x_new(N) = d(N) / b(N);

for i = N-1:-1:1

x_new(i) = (d(i) - c(i) * x_new(i+1)) / b(i);

end

end

Answer 1

[Torsten]

When I run the function, I get:

Not enough input arguments.

Error in solvePDEs (line 27)

y(n) = linspace(y_span(1), y_span(2), num_y);

My guess is that you don't supply any numerical values for the input parameters Gr, Re, Pr, Sc, k, N_c, N_t, R, n, t_span, y_span, num_y, num_t when you "run" the function.

And why do both functions have the same name ? This will give conflicts.

[Sanee]

Thank you, I changed the second function name to ThomasA. Now taking just the second function, what do you think is wrong with this line?

Error in ThomasA (line 38)

phi_new = thomasAlgorithm(Ny, dy, dt/Sc, phi, phi_term - R * phi.^n);

[Torsten]

phi_term - R * phi.^n

phi_term and phi must be vectors of the same length for that the above expression is defined, but they aren't.

[Sanee]

Get it! Thanks.