How can I modify my code from explicit to implicit method ?

3 views (last 30 days)
Hello all,
I have written my code explicitly. I would like to see how my code works implicitly. Can someone help me share their thoughts about it ?
Thank you
clear; clc; close all;
W =5; L=0.25;
Nx=10; Ny=10;
x = linspace(0,L,Nx);
y = linspace(0,W,Ny);
[X, Y] =meshgrid(x,y);
dx = L/Nx; dy = W/Ny;
k1 = dx/dy;
k2 = dx^2/dy^2;
rhoo = 0.8; mu =0.8E-5;
per = 1e-11; %1e-09
u = zeros(Ny, Nx);
v = zeros(Ny, Nx);
P = 101325 * ones(Ny, Nx);
% Boundary conditions
u(1,:) = u(2,:); %0 % left
v(1,:) = v(2,:); % left
P(1,:) = 101325;
u(:,1) = u(:,2); %0 % top - open to atmosphere
v(:,1) = v(:,2); % top - open to atmosphere
P(:,1) = 101325;
u(Ny,:) = u(Ny-1,:); % right - symmetry
v(Ny,:) = v(Ny-1,:); % right - symmetry
P(Ny,:) = P(Ny-1,:); % right - symmetry
u(:,Nx) = 0; % bottom
v(:,Nx) = 0; % bottom
P(:,Nx) = P(:,Nx-1);
S = [
0.1458 0.1929 0.2123 0 0 0 0 0 0 0
0.1458 0.1929 0.2123 0 0 0 0 0 0 0
0.1458 0.1929 0.2123 0 0 0 0 0 0 0
0.1458 0.1929 0.2123 0 0 0 0 0 0 0
0.1458 0.1929 0.2123 0 0 0 0 0 0 0
0.1458 0.1929 0.2123 0 0 0 0 0 0 0
0.1458 0.1929 0.2123 0 0 0 0 0 0 0
0.1458 0.1929 0.2123 0 0 0 0 0 0 0
0.1458 0.1929 0.2123 0 0 0 0 0 0 0
0.1458 0.1929 0.2123 0 0 0 0 0 0 0
];
uold = u; vold = v; Pold = P;
err_uvelocity = 1; err_vvelocity = 1; err_pressure = 1;
while err_uvelocity > 1e-10 || err_vvelocity > 1e-10 || err_pressure > 1e-10
%for test = 1:1
for i=2:Ny-1
for j=2:Nx-1
%%%%%%%%%%%%%%%%%%%%%
u(1,:) = u(2,:); %0 % left
v(1,:) = v(2,:); % left
u(:,1) = u(:,2); %0 % top - open to atmosphere
v(:,1) = v(:,2); % top - open to atmosphere
u(Ny,:) = u(Ny-1,:); % right - symmetry
v(Ny,:) = v(Ny-1,:); % right - symmetry
P(Ny,:) = P(Ny-1,:); % right - symmetry
P(:,Nx) = P(:,Nx-1); % bottom
% P(:,1) = 101325;
%%%%%%%%%%%%%%%%%%%%%
%%Calculating pressure
XX = (S(i,j)*mu*(dx^2))/ (per*rhoo) ;
YY = -k2*(P(i,j-1) + P(i,j+1)) ;
ZZ = -P(i-1,j) - P(i+1,j);
denom = -2 -2*k2;
P(i,j) = (XX + YY + ZZ) / denom;
%%%%%%%%%%%%%%%%%%%%%
%%Calculating velocities
%AlongX = (P(i+1,j) - P(i-1,j))/(2*dx);
%AlongY = (P(i,j+1) - P(i,j-1))/(2*dy);
AlongX = (P(i+1,j) - P(i,j))/(dx);
AlongY = (P(i,j+1) - P(i,j))/(dy);
Const = -per/mu;
u(i,j) = Const*AlongX;
v(i,j) = Const*AlongY;
%%%%%%%%%%%%%%%%%%%%%
err_uvelocity = norm(uold(i,j) - u(i,j));
err_vvelocity = norm(vold(i,j) - v(i,j));
err_pressure = norm(Pold(i,j) - P(i,j));
end
end
uold=u;
vold=v;
Pold = P;
%PPrime = P.^2;
end
% contourf(x, y, v);
% colorbar;
## quiver(x,y,u,v);
  2 Comments
Torsten
Torsten on 3 Mar 2023
Edited: Torsten on 3 Mar 2023
I don't see anything "explicit" in the code above that could be made "implicit".
Seems this is an iteration between pressure and velocity, maybe a pressure-velocity coupling method like SIMPLE or something similar. It might be true that there are explicit and implicit methods to handle this coupling, but I think it's unrealistic that you will get help for this in a MATLAB forum. Maybe you can ask CFD experts.

Sign in to comment.

Answers (0)

Community Treasure Hunt

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

Start Hunting!