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
