MATLAB Answers

0

中心線上の値の出力

Asked by Seiya Arakawa on 23 Oct 2019
Latest activity Commented on by michio
on 24 Oct 2019
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
出力された図です。

  0 Comments

Sign in to comment.

2 Answers

michio
Answer by michio
on 23 Oct 2019

キャビティ流れ、懐かしいですね。私も学生時代に授業でやりました。
さて、
「中心線上の流速 u(0,y),v(x,0) をどうやれば出力できるのかわからない」
とのことですが、もう少し状況を説明いただけますか?
中心線がグリッドと一致していればそこに該当する値をみれば良いと思いますし、もし一致していないのであれば何らかの内挿を実施する必要があると思います。コードがしっかり書かれていることも考えると、その辺でトラブっている様にも思えません。
「ここまで考えてこうやってみたが、この辺でうまくいかない」という形ですと答えも得やすいかと思いますので、よろしくお願いいたします。

  1 Comment

michio
on 24 Oct 2019
コメントを追いやすいようにこちらにて回答します。
*** Seiya Arakawa さんのコメント引用****
丁寧な返信恐縮です。
中心線がグリッドと一致してないと考えたので、内挿を実施しようと考えました。
その場合、おそらく「%%5.結果の出力」の部分に中心流速をプロットするプログラミングを内挿すると考えたのですが、その内挿するプログラムがうまく書けずに困っていました。
(抽象的な質問になってしまって大変心苦しいのですが)
***********************
具体的に何か試したコードがあればと思ったんですが、構想段階ということでしょうか。
u, v それぞれに対して griddedinterpolant が使えるんじゃないかと思います。
ページ内の
「フル グリッドとグリッド ベクトルを使用した 3 次元内挿」
の例題を確認してみてください。特に後半の"グリッド ベクトル" を使用した例がそのまま使えると思います。

Sign in to comment.


Answer by Seiya Arakawa on 23 Oct 2019

丁寧な返信恐縮です。
中心線がグリッドと一致してないと考えたので、内挿を実施しようと考えました。
その場合、おそらく「%%5.結果の出力」の部分に中心流速をプロットするプログラミングを内挿すると考えたのですが、その内挿するプログラムがうまく書けずに困っていました。
(抽象的な質問になってしまって大変心苦しいのですが)

  0 Comments

Sign in to comment.