# 中心線上の値の出力

2 ビュー (過去 30 日間)
Seiya Arakawa 2019 年 10 月 23 日
コメント済み: michio 2019 年 10 月 24 日
MATLABで2次元正方キャビティ内の流れを解析したいのですが

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

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 件のコメント

サインイン to comment.

### 回答 (2 件)

michio 2019 年 10 月 23 日
キャビティ流れ、懐かしいですね。私も学生時代に授業でやりました。
さて、
「中心線上の流速 u(0,y),v(x,0) をどうやれば出力できるのかわからない」
とのことですが、もう少し状況を説明いただけますか？

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

#### 1 件のコメント

michio 2019 年 10 月 24 日
コメントを追いやすいようにこちらにて回答します。
＊＊＊ Seiya Arakawa さんのコメント引用＊＊＊＊

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

u, v それぞれに対して griddedinterpolant が使えるんじゃないかと思います。
ページ内の
「フル グリッドとグリッド ベクトルを使用した 3 次元内挿」
の例題を確認してみてください。特に後半の"グリッド ベクトル" を使用した例がそのまま使えると思います。

サインイン to comment.

Seiya Arakawa 2019 年 10 月 23 日

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

#### 0 件のコメント

サインイン to comment.

サインイン してこの質問に回答します。