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

3 visualizaciones (últimos 30 días)
Kumaresh Kumaresh
Kumaresh Kumaresh el 3 de Mzo. de 2023
Comentada: Kumaresh Kumaresh el 6 de Mzo. de 2023
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 comentarios
Torsten
Torsten el 3 de Mzo. de 2023
Editada: Torsten el 3 de Mzo. de 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.
Kumaresh Kumaresh
Kumaresh Kumaresh el 6 de Mzo. de 2023
Thank you for your comments

Iniciar sesión para comentar.

Respuestas (0)

Community Treasure Hunt

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

Start Hunting!

Translated by