%% Below is a Matlab code to integrate in time the vorticity equation
% First change to $mu=0.1$, $nx=11$, $dt=0.035$, $dpic=0.025$ and $tfinal=0.525$.
% The plot of the streamfunction should show the numerical instability.
% Decreasing $dt$ to 0.3 nearly stabilises the result at this $tfinal$, while $dt=0.25$ works to $tfinal=3.0$.
% This value of $dt$ is the marginal value, but that is for large grids, and small grids are slightly more stable.
% Best to work at $dt=0.2*h*h/mu$
% The code first calculates the boundary conditions to first-order, and then extrapolates these results to apply the conditions to second-order.
% Use this second-order code.
% Change to $mu=0.1$, $nx=11$, $dpic=0.1$ and $tfinal=1.0$ and $dt=0.02$.
% Find the value of $\omega(x=0.5,y=0.5,t=1.0)$.
% Now decrease $dt$ to 0.01, 0.005, 0.0025 and 0.001.
% Then with $nx=15$ try $dt$ = 0.01 (largest stable value), 0.005, 0.0025 and 0.001.
% Finally with $nx=21$ try $dt=$ 0.005 (largest stable value), 0.0025 and 0.001.
% Plot these results for $\psi(x=0.5,y=0.5,t=1.0)$ as a function of $dt$.
% You could experiment by deleting the lines of the code that apply the boundary conditions at second-order.
% Now set $tfinal=3.0$, $dpic=0.1$ and $dt=0.2*h*h/mu$ and obtain $\omega(x=0.5,y=0.5,t)$ to find how long it takes the vorticity to attain a steady value within 4 significant figures.
% Compare your plots for the steady state of the streamfunction and the vorticity.
% Gather results for different spatial resolutions $nx$ for the steady horizontal velocity $u$ at the mid-section $x=0.5 $, and plot on top of one another.
% The code calculates the force on the top plate to second-order accuracy. Find the steady force for different spatial resolutions, $nx=11$, 15, 21, 29 and 41. Plot the force as a function of the grid size $h$.
% Change the top slip boundary condition from $u=\sin^2\pi x$ to $u=1$, i.e change
%% w(i,nx)= -(sin(pi*(i-1)*h)*sin(pi*(i-1)*h)...
%% to w(i,nx) = -1...
% Find the force on the top plate for various spatial resolutions, say $nx= 10$, 14, 20, 28 and 40. Show that the force diverges as the resolution increases as $F = 4.32\ln(1/h) - 3.75$.
%% Streamfunction vorticity code
clear all %clear all variables
clc % clear output screen
%Constants
dpic = 0.5;
tfinal = 3.0;
mu = 0.1; %i.e. Re=10
nx = 10+1;
ny = nx; %different nx & ny gives non-square
itmax = 2*nx;
w = zeros(nx,ny); % form a matrix of zeros (1..nx, 1..ny)
p = w; % w = vorticty, p = streamfunction
wt = w; % wt = dw/dt
%Initialization
h = 1/(nx-1); %h = dx (=dy); so x-size is 1, y-size is ny/nx
t = 0;
r = 2/(1+ pi/(nx-1)); %SOR parameter
dt = 0.2*h*h/mu;
for i= 1:nx
for j=1:ny
p(i,j) = 0;
w(i,j) = 0;
end
end
%% Main program starts here
tpic = t + dpic - 0.5*t;
while(t