中心線上の値の出力

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

出力された図です。

Answer 1

[michio]

キャビティ流れ、懐かしいですね。私も学生時代に授業でやりました。

さて、

「中心線上の流速 u(0,y),v(x,0) をどうやれば出力できるのかわからない」

とのことですが、もう少し状況を説明いただけますか？

中心線がグリッドと一致していればそこに該当する値をみれば良いと思いますし、もし一致していないのであれば何らかの内挿を実施する必要があると思います。コードがしっかり書かれていることも考えると、その辺でトラブっている様にも思えません。

「ここまで考えてこうやってみたが、この辺でうまくいかない」という形ですと答えも得やすいかと思いますので、よろしくお願いいたします。

[michio]

コメントを追いやすいようにこちらにて回答します。

＊＊＊ Seiya Arakawa さんのコメント引用＊＊＊＊

丁寧な返信恐縮です。

中心線がグリッドと一致してないと考えたので、内挿を実施しようと考えました。

その場合、おそらく「%%5.結果の出力」の部分に中心流速をプロットするプログラミングを内挿すると考えたのですが、その内挿するプログラムがうまく書けずに困っていました。

（抽象的な質問になってしまって大変心苦しいのですが）

＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊

具体的に何か試したコードがあればと思ったんですが、構想段階ということでしょうか。

u, v それぞれに対して griddedinterpolant が使えるんじゃないかと思います。

https://jp.mathworks.com/help/matlab/ref/griddedinterpolant.html

ページ内の

「フル グリッドとグリッド ベクトルを使用した 3 次元内挿」

の例題を確認してみてください。特に後半の"グリッド ベクトル" を使用した例がそのまま使えると思います。

Answer 2

[Seiya Arakawa]

丁寧な返信恐縮です。

中心線がグリッドと一致してないと考えたので、内挿を実施しようと考えました。

その場合、おそらく「%%5.結果の出力」の部分に中心流速をプロットするプログラミングを内挿すると考えたのですが、その内挿するプログラムがうまく書けずに困っていました。

（抽象的な質問になってしまって大変心苦しいのですが）