Hi Hesam,
I understand that you are trying to solve a 2D unsteady heat advection diffusion equation with Crank-Nicolson method.
The Taylor-Green vortex is a classic solution in fluid dynamics that describes a particular type of incompressible flow. The Taylor-Green vortex solution is characterized by sinusoidal velocity fields that create a pattern of rotating vortices.
In two dimensions, the velocity components ( 'u' ) and ( 'v' ) of the Taylor-Green vortex are defined as:
u(x, y) = - sin(x) * cos(y)
v(x, y) = cos(x) * sin(y)
Refer to the following code sample to proceed further,
% Parameters
Lx = 2*pi; % Length of the domain in x
Ly = 2*pi; % Length of the domain in y
Nx = 50;
Ny = 50;
alpha = 0.01;
dt = 0.01;
t_final = 1.0;
dx = Lx / (Nx - 1);
dy = Ly / (Ny - 1);
x = linspace(0, Lx, Nx);
y = linspace(0, Ly, Ny);
[X, Y] = meshgrid(x, y);
% Initial condition
T = sin(X) .* sin(Y);
% Precompute coefficients
rx = alpha * dt / (2 * dx^2);
ry = alpha * dt / (2 * dy^2);
% sparse matrices for Crank-Nicolson
Ix = speye(Nx);
Iy = speye(Ny);
ex = ones(Nx, 1);
ey = ones(Ny, 1);
% 1D Laplacian matrices
Lx = spdiags([ex -2*ex ex], [-1 0 1], Nx, Nx);
Ly = spdiags([ey -2*ey ey], [-1 0 1], Ny, Ny);
% Adjust for periodic boundary conditions
Lx(1, end) = 1; Lx(end, 1) = 1;
Ly(1, end) = 1; Ly(end, 1) = 1;
% 2D Laplacian operator
L = kron(Iy, Lx) + kron(Ly, Ix);
% Identity matrix for 2D problem
I = speye(Nx * Ny);
% Crank-Nicolson matrices
A = I - rx * L;
B = I + rx * L;
% Time-stepping loop
T = T(:); % Flatten T for vectorized operations
for t = 0:dt:t_final
% Update velocity field
u = -sin(X) .* cos(Y);
v = cos(X) .* sin(Y);
u = u(:);
v = v(:);
Tx = (circshift(T, -1) - circshift(T, 1)) / (2 * dx);
Ty = (circshift(T, -Nx) - circshift(T, Nx)) / (2 * dy);
RHS = B * T - dt * (u .* Tx + v .* Ty);
T_new = A \ RHS;
T = T_new;
% Reshape and plot
T_plot = reshape(T, [Nx, Ny]);
surf(X, Y, T_plot);
title(['Time = ', num2str(t)]);
xlabel('x');
ylabel('y');
zlabel('Temperature');
drawnow;
end
For a comprehensive understanding of the "circshift", "spdiags" and "speye" functions in MATLAB, please refer to the following documentation.
I hope this helps!
