中心線上の値の出力

Question

Seiya Arakawa 2019-10-23

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/486930-

评论： michio 2019-10-24

在 MATLAB Online 中打开

MATLABで2次元正方キャビティ内の流れを解析したいのですが

中心線上の流速

u(0,y),v(x,0)の出力をどうやれば出力できるのかわからないのでご教授いただければと思います。

解析にはSMAC法を用いました。

clear all
%% 変数の定義
Re = 100; %レイノルズ数
nt = 500; %maximum number of time steps
dt = 1.e-2; %Time step;
lx = 1; ly = 1; %size of domain
nx = 20; ny = 20; %Number of computational cells
alpha = 1.9; %Acceleration parameter for
             %  SOR method
%% Initialization
             
dx = lx/nx; %Size fo computational cell;
dy = ly/ny;
x = ((1:nx)-0.5)*dx; %Coordinates(cell-center);
y = ((1:ny)-0.5)*dy;
%% Initialize arrays
u = zeros(nx+1,ny+2); %Velocity in x direction
v = zeros(nx+2,ny+1); %Velocity in y direction
p = zeros(nx+2,ny+2); %Pressure
phi = zeros(nx+2,ny+2); %Presssure correction
b = zeros(nx+2,ny+2); %RHS of poisson eq.
uc = zeros(nx,ny); %u at cell center
                   % (for plotting)
vc = zeros(nx,ny); %v at cell center
                   % (for plotting)
                   
%% Prepare figure
[X,Y] = meshgrid(x,y); 
h = quiver(X',Y',uc,vc);
%% Main loop
for it = 1:nt
    %% 1.「仮の速度」の計算
    
    for j = 2:ny+1
        for i = 2:nx
            u(i,j) = u(i,j) + dt *( ...
                -((u(i+1,j)+u(i,j))*(u(i+1,j)+u(i,j)) ...
                -(u(i,j)+u(i-1,j))*(u(i,j)+u(i-1,j)))/dx/4. ...
                -((u(i,j+1)+u(i,j))*(v(i+1,j)+v(i,j)) ...
                -(u(i,j)+u(i,j-1))*(v(i+1,j-1)+v(i,j-1)))/dy/4. ...
                -(p(i+1,j)-p(i,j))/dx ...
                +((u(i+1,j)-2.*u(i,j)+u(i-1,j))/dx^2 ...
                +(u(i,j+1)-2.*u(i,j)+u(i,j-1))/dy^2)/Re ...
                );                
        end
    end
    
    for j = 2:ny
        for i = 2:nx+1
            v(i,j) = v(i,j) + dt *( ...
                -((u(i,j+1)+u(i,j))*(v(i+1,j)+v(i,j)) ...
                -(u(i-1,j+1)+u(i-1,j))*(v(i,j)+v(i-1,j)))/dx/4. ...
                -((v(i,j+1)+v(i,j))*(v(i,j+1)+v(i,j)) ...
                -(v(i,j)+v(i,j-1))*(v(i,j)+v(i,j-1)))/dy/4. ...
                -(p(i,j+1)-p(i,j))/dy ...
                +((v(i+1,j)-2.*v(i,j)+v(i-1,j))/dx^2 ...
                +(v(i,j+1)-2.*v(i,j)+v(i,j-1))/dy^2)/Re ...
                );
        end
    end
     
    %% 2.境界条件の適用
    % Apply boundary conditions
    
    u(:,1) = -u(:,2);  v(:,1) = 0.;             %bottom
    u(:,ny+2) = 2.-u(:,ny+1);  v(:,ny+1) = 0.;  %top
    u(1,:) = 0.;  v(1,:) = -v(2,:);             %left
    u(nx+1,:) = 0.;  v(nx+2,:) = -v(nx+1,:);    %left
    
    %% 3.Poisson方程式の求解+速度+圧力修正
    %Correction of velocity and pressure
    %(2nd step of SMAC method)
    
    %RHS of pressure Poisson eq.
    
    for j = 2:ny+1
        for i = 2:nx+1
            b(i,j) = ((u(i,j)-u(i-1,j))/dx ...
                    + (v(i,j)-v(i,j-1))/dy)/dt;
        end
    end
    
    %Solve matrix eq. using SOR method
    
    err = 0.; maxerr = 1.;
    while maxerr > 1.e-8
        
        %Boundary condition for phi(すべてNeumann条件)
        
        phi(:,1) = phi(:,2); %bottom
        phi(:,nx+2) = phi(:,ny+1); %top
        phi(1,:) = phi(2,:); %left
        phi(nx+2,:) = phi(nx+1,:); %right
        
        for j = 2:ny+1
            for i = 2:nx+1
                err = ((phi(i+1,j)-2.*phi(i,j)+phi(i-1,j))/dx^2 ....
                     +(phi(i,j+1)-2.*phi(i,j)+phi(i,j-1))/dy^2 ...
                 -b(i,j));
                 phi(i,j) = phi(i,j) ...
                           + alpha/2./(2./dx^2+2./dy^2)*err;
                       if(maxerr > abs(err))
                          maxerr = abs(err); %残差の最大値を更新
                       end
                       
            end
        end
    end
    
    %Velocity and pressure correction using phi
    
    p = p + phi;
    
    for j = 2:ny+1
        for i = 2:nx
            u(i,j) = u(i,j) - dt * (phi(i+1,j)-phi(i,j))/dx;
        end
    end
    
    for j = 2:ny
        for i = 2:nx+1
            v(i,j) = v(i,j) - dt * (phi(i,j+1)-phi(i,j))/dy;
        end
    end
    
    %% 4.境界条件の適用
    % Apply boundary conditions
    
    u(:,1) = -u(:,2);  v(:,1) = 0.;  %bottom
    u(:,ny+2) = 2.-u(:,ny+1);  v(:,ny+1) = 0.;  %top
    u(1,:) = 0.;  v(1,:) = -v(2,:);  %left
    u(nx+1,:) = 0.;  v(nx+2,:) = -v(nx+1,:);  %left
    
    %% 5.結果の出力
    %Monitor a velocity around the center to check convergence
    u(round(nx/2),round(ny/2))
 
    %Draw vector field
    for j = 1:ny
        for i = 1:nx
            uc(i,j) = (u(i+1,j+1)+u(i,j+1))/2.;
            vc(i,j) = (v(i+1,j+1)+v(i+1,j))/2.;
        end
    end
    xlim([0 lx]); ylim([0 ly]);
    set(h,'Udata',uc,'Vdata',vc);
    drawnow
end