# finite difference method for continuity equation in polar form

13 views (last 30 days)
aktham mansi on 9 Apr 2022
Commented: aktham mansi on 11 Apr 2022
the continuty equation for incompressible and 2 d irrotational flow in polar form is written as:
Using direct method, Lu decomposition
the domain is a circle with inner radius 1 and outer radius 2 and theta with a domain [0; 2*pi].
boundary conditions are u(1; theta) = 1 and u(2; theta) = 0
discretization with uniform (theta*r) (81 * 41), Implement Jacobi,
the discretization equation is : (j index for r and i index for theta)
use a direct solver (LU factorization) to solve. The domain can be imagined as rectangle.
%%%%%LU factorization
clc
clear all
%%%%%%%%%%%%%%% Inputs
r_in=1; % Inside Radius of polar coordinates, r_in, say 1 m
r_out = 2; % Outside Radins of polar coordinates, r_out, say 2 m
j_max = 40; % no. of sections divided between r_in and r_out, choose: 40, 80, 160
dr = (r_out - r_in)/j_max; % section length, m
% total angle = 2*pi
i_max= 80; % no. of angle steps, choose: 80, 160,320
dtheta= 2*pi/i_max; % angle step, rad
r = 1:dr:2;
theta=0:dtheta:2*pi;
[r,theta]=meshgrid(r,theta);
%%%%%initialize solution array
u=zeros(i_max+1,j_max+1);%%%%81*41 matrix
u_0=zeros(i_max+1,j_max+1);%%%% old
u(:,1)=u(:,1)+1; %%%%B.C
u(:,2)=u(:,2)+0;%%%%% B.C
beta=dr^2/dtheta^2;
aktham mansi on 11 Apr 2022
can you help me to Plot Error function versus delta r and delta theta in a log-log scale.
%%%%%jacobi
clc
clear all
%%%%%%%%%%%%%%% Inputs
r_in=1; % Inside Radius of polar coordinates, r_in, say 1 m
r_out = 2; % Outside Radins of polar coordinates, r_out, say 2 m
j_max = 40; % no. of sections divided between r_in and r_out, choose: 40, 80, 160
dr = (r_out - r_in)/j_max; % section length, m
% total angle = 2*pi
i_max= 80; % no. of angle steps, choose: 80, 160,320
dtheta= 2*pi/i_max; % angle step, rad
r = 1:dr:2;
theta=0:dtheta:2*pi;
[r,theta]=meshgrid(r,theta);
%%%%%initialize solution array
u=zeros(i_max+1,j_max+1);%%%%81*41 matrix
u_0=zeros(i_max+1,j_max+1);%%%% old
u(:,1)=u(:,1)+1; %%%%B.C
u(:,2)=u(:,2)+0;%%%%% B.C
u(1,:)=1:-1/j_max:0;
u(i_max+1,:)=1:-1/j_max:0;
beta=dr^2/dtheta^2;
n=1;
k=0;
%%% j index for radius r and i index for theta%%%%
while k==0
u_0=u;
k=1;
for i=2:i_max
for j=2:j_max
r(j)=1+(j-1)*dr;
theta(i)=1+(i-1)*dtheta;
u(i,j)=(((r(j)+r(j+1))/2)*u_0(i,j+1)+((r(j)+r(j-1))/2)*u_0(i,j-1)+beta*u_0(i+1,j)+beta*u_0(i-1,j))/(((r(j)+r(j+1))/2)+((r(j)+r(j-1))/2)+2*beta);
error(i,j)=abs(u(i,j)-u_0(i,j))/sqrt(i_max*j_max);
if abs(u(i,j)-u_0(i,j))>(10^-5)
k=0;
end
end
end
n=n+1;
end
n
r = linspace(1,2,j_max+1);
theta = linspace(0,2*pi,i_max+1);
[r,theta] = meshgrid(r,theta);
uana = -log(r)/log(2) + 1; %%%%%analytical solution
[x,y]=pol2cart(theta,r);
figure(1);
loglog(error);
figure(2);
surface(x,y,uana);%%%%%analytical solution
title('analytical solution');
colorbar;
figure(3);
surface(x,y,u);%%%%%%% FDM jacobi solution
title('FDM JACOBI solution');
colorbar;

### Categories

Find more on Matrix Indexing in Help Center and File Exchange

R2018a

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by