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

7 visualizzazioni (ultimi 30 giorni)
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 Commenti
Torsten
Torsten il 3 Mar 2023
Modificato: Torsten il 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.

Accedi per commentare.

Risposte (0)

Community Treasure Hunt

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

Start Hunting!

Translated by