driven_cavity_fds_rk2.m

Matlab script for solution to the driven cavity problem using the FDS algorithm, which discretizes t

% McDermott
% 1-15-2007
% driven_cavity_fds_RK2.m
%
% Simple staggered primitive variable treatment for the driven
% cavity problem using an RK2 time integrator.

close all
clear all

Lx = 1; % cavity dimension in x
Ly = 0.25; % cavity dimension in y
nx = 64; % number of pcells in x
ny = 16; % number of pcells in y
dx = Lx/nx; % uniform grid spacing in x
dy = Ly/ny; % uniform grid spacing in y
dxdx = dx^2;
dydy = dy^2;
U0 = 1; % lid velocity
nu = 1/10000; % kinematic viscosity (1/Re)
Fo = 0.5; % Fourier number
CFL = 0.5; % Courant number
dt_dif = Fo/(nu*(1/dxdx+1/dydy)) % time step based on Fo
dt_adv = CFL*min(dx,dy)/U0 % time step based on CFL
dt = min([dt_dif,dt_adv]) % time step used in Forward Euler integrator
T = 1000*dt; % total simulation time

onethird = 1/3;
eightthirds = 8/3;

% STAGGERED GRID ARRANGEMENT:
%
% Let this represent the bottom left pcell control volume.
% Velocities are stored on their respective faces.  Normal stresses
% and normal advective fluxes are stored at the pcell center.
% Off-diagonal stresses and advective fluxes are stored at vertices
% marked by the Xs.  The bottom left most element is prescribed the
% indices (1,1).  Because of this, some of the differencing
% below looks confusing, but this seems somewhat unavoidable.
%
%       omega3(1,2)                              omega3(2,2)
%       tau12(1,2)          ^ v(1,2)             tau12(2,2)
%       X-------------------|--------------------X
%       |                                        |
%       |                                        |
%       |                                        |
%       |                                        |
%       |                                        |
%       |                                        |
%       |                                        |
%       |                                        |
%      ---> u(1,1)          O                   ---> u(2,1)
%       |                  H(1,1)                |
%       |                  tau11(1,1)            |
%       |                  tau22(1,1)            |
%       |                                        |
%       |                                        |
%       |                                        |
%       |                                        |
%       |                   ^ v(1,1)             |
%       X-------------------|--------------------X
%       omega3(1,1)                              omega3(2,1)
%       tau12(1,1)                               tau12(2,1)
%

% initialization
u = zeros(nx+1,ny);
v = zeros(nx,ny+1);
u_star = zeros(nx+1,ny);
v_star = zeros(nx,ny+1);
tau11 = zeros(nx,ny);
tau22 = zeros(nx,ny);
tau12 = zeros(nx+1,ny+1);
omega3 = zeros(nx+1,ny+1);
H = zeros(nx,ny);
b_vec = zeros(nx*ny,1);

% grid
x = linspace(dx/2,Lx-dx/2,nx);
y = linspace(dy/2,Ly-dy/2,ny);
for j = 1:ny
    X(:,j) = x;
end
for i = 1:nx
    Y(i,:) = y;
end

% Build A
A = build_sparse_matrix_neumann(nx,ny,dx,dy);


for t = dt:dt:T
    
    % STAGE 1
    u0 = u;
    v0 = v;

    % compute normal stresses
    for j = 1:ny
        for i = 1:nx     
            tau11(i,j) = -2*nu*( u(i+1,j)-u(i,j) )/dx;
            tau22(i,j) = -2*nu*( v(i,j+1)-v(i,j) )/dy;
        end
    end
    
    % compute interior vorticity and off-diagonal stresses
    for j = 2:ny
        for i = 2:nx
            omega3(i,j) = (v(i,j)-v(i-1,j))/dx - (u(i,j)-u(i,j-1))/dy;
            tau12(i,j) = -nu*( (u(i,j)-u(i,j-1))/dy + (v(i,j)-v(i-1,j))/dx );
        end
    end
    
    % bottom wall stress and vorticity
    j = 1;
    for i = 2:nx
        dudy_wall = ( 0 + 3*u(i,j) - onethird*u(i,j+1) )/dy;
        tau12(i,j) = -nu*dudy_wall;
        omega3(i,j) = -dudy_wall;
    end
    % top wall stress and vorticity
    j = ny+1;
    for i = 2:nx
        dudy_wall = ( eightthirds*U0 - 3*u(i,ny) + onethird*u(i,ny-1) )/dy;
        tau12(i,j) = -nu*dudy_wall;
        omega3(i,j) = -dudy_wall;
    end
    % left wall stress and vorticity
    i = 1;
    for j = 2:ny
        dvdx_wall = ( 0 + 3*v(i,j) - onethird*v(i+1,j) )/dx;
        tau12(i,j) = -nu*dvdx_wall;
        omega3(i,j) = dvdx_wall;
    end
    % right wall stress and vorticity
    i = nx+1;
    for j = 2:ny
        dvdx_wall = ( 0 - 3*v(nx,j) + onethird*v(nx-1,j) )/dx;
        tau12(i,j) = -nu*dvdx_wall;
        omega3(i,j) = dvdx_wall;
    end

    
    % u momentum -- predictor
    for j = 1:ny
        for i = 2:nx-1
            
            Fu_visc = -( (tau11(i,j)-tau11(i-1,j))/dx + (tau12(i,j+1)-tau12(i,j))/dy );
            
            vn = 0.5*( v(i-1,j+1) + v(i,j+1) );
            vs = 0.5*( v(i-1,j)   + v(i,j)   );
            
            eps_vs = 0; %vs*dt/dy;
            eps_vn = 0; %vn*dt/dy;
            Fx = 0.5*( (1+eps_vs)*vs*omega3(i,j) + (1-eps_vn)*vn*omega3(i,j+1) );
             
            u_star(i,j) = u(i,j) + dt*(Fx + Fu_visc);
        end
    end
    
    % v momentum -- predictor
    for j = 2:ny-1
        for i = 1:nx
            
            Fv_visc = -( (tau12(i+1,j)-tau12(i,j))/dx + (tau22(i,j)-tau22(i,j-1))/dy );
            
            ue = 0.5*( u(i+1,j-1) + u(i+1,j) );
            uw = 0.5*( u(i,j-1)   + u(i,j)   );
            
            eps_uw = 0; %uw*dt/dx;
            eps_ue = 0; %ue*dt/dx;
            Fy = -0.5*( (1+eps_uw)*uw*omega3(i,j) + (1-eps_ue)*ue*omega3(i+1,j) );
            
            v_star(i,j) = v(i,j) + dt*(Fy + Fv_visc);
        end
    end


    % build source
    for j = 1:ny
        for i = 1:nx
            b(i,j) = ( u_star(i+1,j)-u_star(i,j) )/dx + ( v_star(i,j+1)-v_star(i,j) )/dy;
        end
    end
    b = b/dt;

    % map b to source vector
    for j = 1:ny
        for i = 1:nx
            p = (j-1)*nx+i;
            b_vec(p) = b(i,j);
        end
    end
    % subtract mean for discrete compatibility condition
    b_vec = b_vec - mean(b_vec);

    % solve for H vector
    H_vec = A\b_vec;
    % subtract mean for arbitrary solution
    H_vec = H_vec - mean(H_vec);
    H_min = min(H_vec);
    H_max = max(H_vec);

    % map H_vec to grid
    for j = 1:ny
        for i = 1:nx
            p = (j-1)*nx+i;
            H(i,j) = H_vec(p);
        end
    end

    % project velocities
    for j = 1:ny
        for i = 2:nx
            u(i,j) = u_star(i,j) - dt*( H(i,j)-H(i-1,j) )/dx;
        end
    end
    for j = 2:ny
        for i = 1:nx
            v(i,j) = v_star(i,j) - dt*( H(i,j)-H(i,j-1) )/dy;
        end
    end

    % check divergence
    for j = 1:ny
        for i = 1:nx
            b(i,j) = ( u(i+1,j)-u(i,j) )/dx + ( v(i,j+1)-v(i,j) )/dy;
        end
    end
    display(['max stage 1 velocity divergence = ',num2str( max(max(abs(b))) )])
    
    % STAGE 2
    
    % compute normal stresses
    for j = 1:ny
        for i = 1:nx     
            tau11(i,j) = -2*nu*( u(i+1,j)-u(i,j) )/dx;
            tau22(i,j) = -2*nu*( v(i,j+1)-v(i,j) )/dy;
        end
    end
    
    % compute interior vorticity and off-diagonal stresses
    for j = 2:ny
        for i = 2:nx
            omega3(i,j) = (v(i,j)-v(i-1,j))/dx - (u(i,j)-u(i,j-1))/dy;
            tau12(i,j) = -nu*( (u(i,j)-u(i,j-1))/dy + (v(i,j)-v(i-1,j))/dx );
        end
    end
    
    % bottom wall stress and vorticity
    j = 1;
    for i = 2:nx
        dudy_wall = ( 0 + 3*u(i,j) - onethird*u(i,j+1) )/dy;
        tau12(i,j) = -nu*dudy_wall;
        omega3(i,j) = -dudy_wall;
    end
    % top wall stress and vorticity
    j = ny+1;
    for i = 2:nx
        dudy_wall = ( eightthirds*U0 - 3*u(i,ny) + onethird*u(i,ny-1) )/dy;
        tau12(i,j) = -nu*dudy_wall;
        omega3(i,j) = -dudy_wall;
    end
    % left wall stress and vorticity
    i = 1;
    for j = 2:ny
        dvdx_wall = ( 0 + 3*v(i,j) - onethird*v(i+1,j) )/dx;
        tau12(i,j) = -nu*dvdx_wall;
        omega3(i,j) = dvdx_wall;
    end
    % right wall stress and vorticity
    i = nx+1;
    for j = 2:ny
        dvdx_wall = ( 0 - 3*v(nx,j) + onethird*v(nx-1,j) )/dx;
        tau12(i,j) = -nu*dvdx_wall;
        omega3(i,j) = dvdx_wall;
    end
    
    
    
    % u momentum -- corrector
    for j = 1:ny
        for i = 2:nx-1
            
            Fu_visc = -( (tau11(i,j)-tau11(i-1,j))/dx + (tau12(i,j+1)-tau12(i,j))/dy );
            
            vn = 0.5*( v(i-1,j+1) + v(i,j+1) );
            vs = 0.5*( v(i-1,j)   + v(i,j)   );
            
            eps_vs = 0; %vs*dt/dy;
            eps_vn = 0; %vn*dt/dy;
            Fx = 0.5*( (1-eps_vs)*vs*omega3(i,j) + (1+eps_vn)*vn*omega3(i,j+1) );
            
            u_star(i,j) = 0.5*( u0(i,j) + u(i,j) + dt*(Fx + Fu_visc) );
        end
    end
    
    % v momentum
    for j = 2:ny-1
        for i = 1:nx
            
            Fv_visc = -( (tau12(i+1,j)-tau12(i,j))/dx + (tau22(i,j)-tau22(i,j-1))/dy );
            
            ue = 0.5*( u(i+1,j-1) + u(i+1,j) );
            uw = 0.5*( u(i,j-1)   + u(i,j)   );
            
            eps_uw = 0; %uw*dt/dx;
            eps_ue = 0; %ue*dt/dx;
            Fy = -0.5*( (1-eps_uw)*uw*omega3(i,j) + (1+eps_ue)*ue*omega3(i+1,j) );
             
            v_star(i,j) = 0.5*( v0(i,j) + v(i,j) + dt*(Fy + Fv_visc) );
        end
    end

    % build source
    for j = 1:ny
        for i = 1:nx
            b(i,j) = ( u_star(i+1,j)-u_star(i,j) )/dx + ( v_star(i,j+1)-v_star(i,j) )/dy;
        end
    end
    b = b/dt;

    % map b to source vector
    for j = 1:ny
        for i = 1:nx
            p = (j-1)*nx+i;
            b_vec(p) = b(i,j);
        end
    end
    % subtract mean for discrete compatibility condition
    b_vec = b_vec - mean(b_vec);

    % solve for H vector
    H_vec = A\b_vec;
    % subtract mean for arbitrary solution
    H_vec = H_vec - mean(H_vec);
    H_min = min(H_vec);
    H_max = max(H_vec);

    % map H_vec to grid
    for j = 1:ny
        for i = 1:nx
            p = (j-1)*nx+i;
            H(i,j) = H_vec(p);
        end
    end

    % project velocities
    for j = 1:ny
        for i = 2:nx
            u(i,j) = u_star(i,j) - dt*( H(i,j)-H(i-1,j) )/dx;
        end
    end
    for j = 2:ny
        for i = 1:nx
            v(i,j) = v_star(i,j) - dt*( H(i,j)-H(i,j-1) )/dy;
        end
    end

    % check divergence
    for j = 1:ny
        for i = 1:nx
            b(i,j) = ( u(i+1,j)-u(i,j) )/dx + ( v(i,j+1)-v(i,j) )/dy;
        end
    end
    display(['max stage 2 velocity divergence = ',num2str( max(max(abs(b))) )])
    
    
    % interpolate velocities to cell centers
    for j = 1:ny
        for i = 1:nx
            up(i,j) = 0.5*( u(i+1,j)+u(i,j) );
            vp(i,j) = 0.5*( v(i,j+1)+v(i,j) );
        end
    end
    
    SUBPLOT(2,1,1)
    velmag = sqrt( up.*up + vp.*vp );
    V = linspace(0,1,20);
    %contourf(X,Y,velmag,V); hold on
    %colorbar
    quiver(X,Y,up,vp); hold off
    title('velocity vectors')
    axis([0 Lx 0 Ly])
    set(gca,'PlotBoxAspectRatio',[Lx Ly 1])
    
    SUBPLOT(2,1,2)
    V = linspace(H_min,H_max,20);
    contourf(X,Y,H,V)
    title('pressure contours')
    colorbar
    set(gca,'PlotBoxAspectRatio',[Lx Ly 1])
    
    pause(0.001)

end