%% Below is a Matlab code for the primitive variable formulation with a staggered grid.
%% First try some tests to check the stability of time-stepping, the second-order accuracy, and the time to reach the steady state.
%% Gather the results for different spatial resolutions $nx$ for the steady horizontal velocity at the mid-section $x=0.5$, and plot these results on top of one another. Include a result from the streamfunction vorticity formulation.
%% The code calculates the force to second-order accuracy on the top plate, on the sides and on the bottom.
%% Experiment with different geometries of the cavity, say $ny=2*nx$.
% The optimal parameter for the SOR relaxation will change with the geometry, but this should not affect the calculation of the steady state.
% The force on the bottom plate should decrease rapidly with deeper cavities, while the force on the top plate should decrease only slightly. On the other hand for shallower cavities, the force on the top plate increases as the height decreases.
%% Experiment with different Reynolds numbers, $1/mu$, say $mu=0.2$, 0.5, and 0.05. The code reduces the time-step as $mu$ is increased. That will increase the run-time, although the steady state should be reached earlier. At smaller values of $mu$ the time-step may be limited by the CFL condition if the spatial resolution is not increased. The forces change very little in $mu>0.1$, i.e.~are effectively in the low Reynolds number limit. For higher Reynolds numbers, $mu<0.02$, the forces increase as the top boundary layer begins to thin. Note at these higher Reynolds number the steady state is achieved later, say by $t=10$.
%% primative variable on staggered grid, algorithm 3
clear all %clear all variables
clc % clear output screen
%% Constants
mu = 0.1; %i.e. Re=10
nx = 11+1; % +1 because index 1..nx not 0..nx
ny = nx; %different nx & ny gives non-square
tfinal = 3.0;
dtpic = 0.5;
w = zeros(nx,ny); % form a matrix of zeros (1..nx, 1..ny)
p = w; % w = vorticty, p = streamfunction
u = w; % u velocty, sim v
v = w;
ut = w; % ut = du/dt, sim vt
vt = w;
bcb = zeros(nx); % for BC on Bottom. Top, Left and Right
bct = bcb;
bcl = zeros(ny);
bcr = bcl;
%% Initialization
h = 1/(nx-2); %h = dx (=dy); so x-size is 1, y-size is ny/nx
t = 0;
r = 2/(1+ 0.8/(nx-1)); %SOR parameter
dt = 0.2*h*h/mu;
tpic = t + dtpic - 0.5*dt;
for i= 1:nx
for j=1:ny
p(i,j) = 0;
u(i,j) = 0;
v(i,j) = 0;
end
end
%% Main program starts here
while(ttpic
% Force calculation
% top, uses u = 0 at ends to save correction to int
q = 0;
for i = 1:nx-1
q = q + (u(i,ny)-u(i,ny-1));
end
% bottom
qb = 0;
for i = 1:nx-1
qb = qb - (u(i,2)-u(i,1));
end
% sides
qs = 0;
for j = 2:ny-1
qs = qs + h*0.5*(p(1,j)+p(1,j) - p(nx-1,j)-p(nx,j))/mu;
end
formatSpecA='t= %8.4f; top= %8.4f;sides=%8.4f;bottom=%8.4f; -qs-qb=%8.4f\n ';
fprintf(formatSpecA, t, q, qs, qb, -qs-qb);
% End of Force calculation
tpic = t + dtpic - 0.5*dt;
end % endif t>tpic
t = t + dt;
end %% end of time-step, endwhile of t