Coding the Thomas algorithm for a heat and mass transfer problem
8 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
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
0 Commenti
Risposta accettata
Torsten
il 26 Giu 2024
Modificato: Torsten
il 26 Giu 2024
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.
3 Commenti
Torsten
il 27 Giu 2024
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.
Più risposte (0)
Vedere anche
Categorie
Scopri di più su Startup and Shutdown in Help Center e File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!