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

Question

Kumaresh Kumaresh 2023-3-3

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1922205-how-can-i-modify-my-code-from-explicit-to-implicit-method

评论： 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 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 2023-3-6