2D unsteady heat advection diffusion equation with Crank-Nicolson method scheme

how can i solve a 2D unsteady heat advection diffusion equation with Crank-Nicolson method scheme using Matlab? the convective flows are given by Taylor-Green vortex solution.

Answer 1

[SAI SRUJAN]

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.

https://www.mathworks.com/help/matlab/ref/circshift.html

https://www.mathworks.com/help/matlab/ref/speye.html

https://www.mathworks.com/help/matlab/ref/spdiags.html

I hope this helps!