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

조회 수: 8 (최근 30일)
Kumaresh Kumaresh
Kumaresh Kumaresh 2023년 3월 3일
댓글: Kumaresh Kumaresh 2023년 3월 6일
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
Torsten
Torsten 2023년 3월 3일
편집: Torsten 2023년 3월 3일
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 2023년 3월 6일
Thank you for your comments

댓글을 달려면 로그인하십시오.

답변 (0개)

Community Treasure Hunt

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

Start Hunting!

Translated by