sandbox/easystab/potential.m

    Potential flow in 2D

    This is a code to solve the equation for potential flow in a square.

    clear all; clf

% parameters
Nx=10; % gridpoints in x
Ny=10; % gridpoints in y 
Lx=2; % domain size in x
Ly=2; % domain size in y
pts=5; % number of points in finite difference stencils
method='fd'; % discretization

% differentiation 
[d.x,d.xx,d.wx,x]=dif1D(method,-1,Lx,Nx,5);
[d.y,d.yy,d.wy,y]=dif1D(method,-1,Ly,Ny,5);
[D,l,X,Y,Z,I,NN]=dif2D(d,x,y);

    The equation for the velocity potential \phi(x,y) is \displaystyle \phi_{xx}+\phi_{yy}=0. 

    % system matrices
A=D.xx+D.yy;
b=zeros(NN,1);

    We start with the known analystical solution of a potential flow \displaystyle \phi=xy thus u(x,y)=y and v(x,y)=x.

    % initial guess
phi0=X.*Y;

    The boundary conditions are the velocity normal to the walls. The components of the velocity are \displaystyle u=\phi_x and \displaystyle v=\phi_y.

    We use the analytical expression of the flow to impose the value of the nonhomogeneous boundary conditions.

    % boundary conditions
loc=[l.top;l.left;l.bot;l.right];

C=[D.x([l.left;l.right],:); ...
   D.y([l.bot;l.top],:)];

A(loc,:)=C;
b(loc,:)=C*phi0(:);

% solve system
phi=A\b;

% validation
u=reshape(D.x*phi,Nx,Ny);
v=reshape(D.y*phi,Nx,Ny);

utheo=reshape(D.x*phi0(:),Nx,Ny);
vtheo=reshape(D.y*phi0(:),Nx,Ny);

% show velocity field
subplot(1,2,1);
quiver(X,Y,u,v,'b');
hold on
quiver(X,Y,utheo,vtheo,'r');
axis equal; axis([-1 -1+Lx -1 -1+Ly])
grid on; 
xlabel('x'); ylabel('y'); title('velocity field')

% show error
subplot(1,2,2);
phi=phi-phi(1)+phi0(1);
mesh(X,Y,abs(reshape(phi,Ny,Nx)-phi0));
xlabel('x'); ylabel('y'); title('phi error')

set(gcf,'paperpositionmode','auto');
print('-dpng','-r80','potential.png');

    The error is close to machine accuracy even thought the resolution is low since the solution \phi is linear. The two velocity fields (numerical and computed) are on top of each other. The figure