12-2. 格子生成の基礎：2D格子 - 凪またりの垂れ流し

はじめに
2D格子とは何か？
2Dで平板流れの初期値を設定してみる
おわりに：2Dの世界へ

はじめに

こんにちは、凪またりです。

今回は、前回の1D格子を2D（2次元）に拡張します。基本的には1Dで実施した同じことを繰り返すだけなので簡単です。3D（3次元）への拡張も同じことをするだけになります。

他のトピックについては『こちら』でまとめていますので、初めてご覧いただいた方はご確認いただければと思います。

2D格子とは何か？

2D格子は、平面（ $x-y$ 座標）を小さなセルに分割したものです。
例えば、長さ $L=1$ m、高さ $H=0.2$ mの領域があれば、それを $n_x×n_y$ の格子に分割して、各点で物理量（ $ρ, u_x, u_y, p, T$ ）を定義します。1Dでは $x$ 方向だけでしたが、2Dでは $y$ 方向が加わり、流れの広がりや境界層を表現できます。
平板流れでイメージするならば、平板が $x$ 軸に沿って置かれ、空気が $x$ 方向に流れつつ、 $y$ 方向では壁面から離れるにつれて速度が変化する感じです。

2Dで平板流れの初期値を設定してみる

均一2D格子の生成

1Dと同様、まずは「均一格子」を作ってみます。
$L=1$ m、 $H=0.2$ m を $n_x=10$ 、 $n_y=5$ で分割すると、

$Δx=L/n_x=0.1$ m
$Δy=H/n_y=0.04$ m
格子点座標： $x_i=(i-1)Δx$ , 　 $y_j = (j-1)Δy$

です。これを基に、Fortranで格子点の座標を設定してみます。

program grid_2d
implicit none
integer :: i, j, nx, ny
real :: L, H, dx, dy
real, dimension(:,:), allocatable :: x, y

! 格子条件を設定
nx = 10
ny = 5
L = 1.0 ! x方向長さ (m)
H = 0.2 ! y方向高さ (m)
dx = L / real(nx)
dy = H / real(ny)

allocate(x(nx+1, ny+1))
allocate(y(nx+1, ny+1))

! 格子の座標を設定
do i = 1, nx+1
do j = 1, ny+1
x(i, j) = (i-1) * dx
y(i, j) = (j-1) * dy
end do
end do

print *, "dx:", dx, "dy:", dy
do j = 1, ny+1
print *, "y:", y(1, j), "x:", x(:, j)
end do

deallocate(x, y)
end program grid_2d

実行した結果は⇩です。

これで、2D格子が完成しました。

格子に物理量を割り当てる

流れを計算するには、まず各格子点に物理量を割り当てることが必要なので、実際に割り当ててみます。
今回は、⇩のイメージ図のような平板上の流れの初期条件を模擬して設定してみます。
なお、板の壁面付近には境界層が存在するので、簡単にそれっぽい式で表現しています。

program grid_with_data
implicit none
integer :: i, j, nx, ny
real :: L, H, dx, dy, u_in, rho_in
real, dimension(:,:), allocatable :: x, y, u_x, rho

! 計算領域の作成
nx = 10
ny = 5
L = 1.0
H = 0.2
dx = L / real(nx)
dy = H / real(ny)
u_in = 170.0 ! 流入速度 (m/s, Ma = 0.5)
rho_in = 1.2 ! 流入密度 (kg/m^3)

allocate(x(nx+1, ny+1))
allocate(y(nx+1, ny+1))
allocate(u_x(nx+1, ny+1))
allocate(rho(nx+1, ny+1))

! 格子生成、速度場（u_x）と密度場（ρ）の定義
do i = 1, nx+1
do j = 1, ny+1
x(i, j) = (i-1) * dx
y(i, j) = (j-1) * dy
u_x(i, j) = u_in * (1.0 - exp(-10.0 * y(i, j))) ! 壁付近で遅くて、上に行くほど速くなる
rho(i, j) = rho_in * (1.0 - 0.05 * y(i, j) / H) ! 上に行くほど少し密度が減る
end do
end do

do j = 1, ny+1
print *, "y:", y(1, j), "u_x:", u_x(1, j), "rho:", rho(1, j)
end do

deallocate(x, y, u_x, rho)
end program grid_with_data

$u_x$ ：壁面（ $y=0$ ）で0、遠くで $u_{in}$ に近づく指数関数で境界層を模擬
$ρ$ ： $y$ 方向で5%減少。圧縮性で圧力や温度が影響

実行した結果は⇩です。

平板流れの初期場ができました。

グラフで見るとこんな感じです。( $y$ 座標が縦軸です)

データの出力

2Dデータをファイルに保存します。ParaViewで読めるように.vtk形式で出力します。
先ほどのコードにオレンジ色部分を付け足せばOKです。
grid.vtkが生成されるので、メモ帳などで中身を確認してみてください。

program grid_output
implicit none
integer :: i, j, nx, ny
real :: L, H, dx, dy, u_in, rho_in
real, dimension(:,:), allocatable :: x, y, u_x, rho
integer :: unit_xyz

! 計算領域の作成
nx = 10
ny = 5
L = 1.0
H = 0.2
dx = L / real(nx)
dy = H / real(ny)
u_in = 170.0 ! 流入速度 (m/s, Ma = 0.5)
rho_in = 1.2 ! 流入密度 (kg/m^3)

allocate(x(nx+1, ny+1))
allocate(y(nx+1, ny+1))
allocate(u_x(nx+1, ny+1))
allocate(rho(nx+1, ny+1))

! 格子生成、速度場（u_x）と密度場（ρ）の定義
do i = 1, nx+1
do j = 1, ny+1
x(i, j) = (i-1) * dx
y(i, j) = (j-1) * dy
u_x(i, j) = u_in * (1.0 - exp(-10.0 * y(i, j))) ! 壁付近で遅くて、上に行くほど速くなる
rho(i, j) = rho_in * (1.0 - 0.05 * y(i, j) / H) ! 上に行くほど少し密度が減る
end do
end do

! ===== ParaView用出力 ===== !

! == XYZファイル（座標） == !
open(newunit=unit_xyz, file="grid.xyz", status="replace")

! 格子サイズ
write(unit_xyz,*) nx+1, ny+1, 1

! x 座標
do j = 1, ny+1
do i = 1, nx+1
write(unit_xyz,'(E15.7)') x(i,j)
end do
end do

! y 座標
do j = 1, ny+1
do i = 1, nx+1
write(unit_xyz,'(E15.7)') y(i,j)
end do
end do

! z 座標
do j = 1, ny+1
do i = 1, nx+1
write(unit_xyz,'(E15.7)') 0.0
end do
end do

close(unit_xyz)

print *, "grid.xyz written"

deallocate(x, y, u_x, rho)

end program grid_output

ParaViewで見ると下記のような格子が出来ています。

※ParaViewを開いて、左上にあるFileからOpenでunit.xyzを開いて「Apply」を押すと格子ファイルが読み込まれます。最初はSurfaceでの表示になっているかもですが、上図はWireframeで表示しています。
※操作方法ですが、「Ctrl+右クリック」で回転できたり、移動は「左クリック」でできたり、時には「Shift+左クリック」で回転だったり、何基準で操作が決まっているかよく知りません...。