finite difference method for continuity equation in polar form

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]

@Torsten

@VBBV

@Kaid Noureddine

[Torsten]

And how do you plan to discretize the periodic boundary condition at theta=0 and theta = 2*pi and the condition

du/dr = 0 at r = 0 ?

[aktham mansi]

@Torsten

yes, i solved it using iteratition, the results are like that

[Torsten]

Little problem with the periodicity ?

[aktham mansi]

@Torsten

yes, i faced problem ar theta=0 or 2*pi

i used this code for jacobi

https://www.mathworks.com/matlabcentral/answers/1691530-finite-difference-method-scheme

the analytical results:

[aktham mansi]

@Torsten

[aktham mansi]

@Torsten

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;