diff --git a/docs_part1/vanilaCavityFlow_JP.md b/docs_part1/vanilaCavityFlow_JP.md index fd6be6e..8bab354 100644 --- a/docs_part1/vanilaCavityFlow_JP.md +++ b/docs_part1/vanilaCavityFlow_JP.md @@ -115,7 +115,7 @@ CFD に関わる方にとっては基本的なことかと思いますが、自 この手法自体は Harlow and Welch (1965) [2] や Chorin (1968) [3] などで提案されたと紹介されることが多く、MAC (Marker and Cell method) 法とも呼ばれます。圧力項を残したバージョンは SMAC (Simplified MAC) とか。この辺の歴史や呼び名は正直はっきり理解できていないのですが、詳しい方コメントください!時間発展の高精度化・安定化を図った Kim and Moin (1985) [4] の手法も実用的です。 -## 圧力の境界条件 +## 仮の速度場の境界条件 ちなみにこの手法は @@ -131,7 +131,7 @@ CFD に関わる方にとっては基本的なことかと思いますが、自 -個人的には Perot (1993) [5] や Chang (2002) [7] で紹介されている、離散化した Navier-Stokes 方程式の LU 分解をベースにした議論が気に入っているので、これをここでは採用します。**圧力の境界条件は不要!**と鮮やかな見解を提示しています。 +個人的には Perot (1993) [5] や Chang (2002) [7] で紹介されている、離散化した Navier-Stokes 方程式の LU 分解をベースにした議論が気に入っているので、これをここでは採用します。**仮の速度場の境界条件なんて考える必要もない!**と鮮やかな見解を提示しています。 # 離散化手法 diff --git a/docs_part1/vanilaCavityFlow_JP.mlx b/docs_part1/vanilaCavityFlow_JP.mlx index 8a25f3e..d3f9a42 100644 Binary files a/docs_part1/vanilaCavityFlow_JP.mlx and b/docs_part1/vanilaCavityFlow_JP.mlx differ diff --git a/functions/getIntermediateU_xyRK3_dst.m b/functions/getIntermediateU_xyRK3_dst.m new file mode 100644 index 0000000..3cde381 --- /dev/null +++ b/functions/getIntermediateU_xyRK3_dst.m @@ -0,0 +1,30 @@ +% Copyright 2020 The MathWorks, Inc. +function xu = getIntermediateU_xyRK3_dst(u, b, dt, Re, nx, ny, dx, dy, id) + +kx = [1:nx-1]'; +ax = pi*kx/(nx); +mwx = 2*(cos(ax)-1)/dx^2;% DST-I + +ky = [1:ny]'; +ay = pi*ky/(ny); +mwy = 2*(cos(ay)-1)/dy^2; % DST-II + +mw = mwx+mwy'; % Modified Wavenumber + +% 各 column 毎に変換されるので、 +% column 方向(x 方向に)dst1, 転置して +% y 方向に dst2、転置し返して元に戻す処理 +rhshat = mydst2(mydst1(b)')'; + +switch id + case 1 + uhat = rhshat./(1-(4/15)*dt/(Re)*mw); + case 2 + uhat = rhshat./(1-(1/15)*dt/(Re)*mw); + case 3 + uhat = rhshat./(1-(1/6)*dt/(Re)*mw); +end + +xu = mydst1(mydst3(uhat')'); + +end \ No newline at end of file diff --git a/functions/getIntermediateV_xyRK3_dst.m b/functions/getIntermediateV_xyRK3_dst.m new file mode 100644 index 0000000..2fbbcfe --- /dev/null +++ b/functions/getIntermediateV_xyRK3_dst.m @@ -0,0 +1,32 @@ +% Copyright 2020 The MathWorks, Inc. +function xv = getIntermediateV_xyRK3_dst(v, b, dt, Re, nx, ny, dx, dy,id) + + +kx = [1:nx]'; +ax = pi*kx/(nx); +mwx = 2*(cos(ax)-1)/dx^2;% DST-I + +ky = [1:ny-1]'; +ay = pi*ky/(ny); +mwy = 2*(cos(ay)-1)/dy^2; % DST-II + +mw = mwx+mwy'; % Modified Wavenumber + +% 各 column 毎に変換されるので、 +% column 方向(x 方向に)dst2, 転置して +% y 方向に dst1、転置し返して元に戻す処理 +rhshat = mydst1(mydst2(b)')'; + +switch id + case 1 + vhat = rhshat./(1-(4/15)*dt/(Re)*mw); + case 2 + vhat = rhshat./(1-(1/15)*dt/(Re)*mw); + case 3 + vhat = rhshat./(1-(1/6)*dt/(Re)*mw); +end + +xv = mydst3(mydst1(vhat')'); + +end + diff --git a/functions/updateVelocityField_RK3.m b/functions/updateVelocityField_RK3.m index b7afc4f..6b96d19 100644 --- a/functions/updateVelocityField_RK3.m +++ b/functions/updateVelocityField_RK3.m @@ -81,10 +81,12 @@ % Implicit treatment for xy direction b = u(2:end-1,2:end-1) - dt*(gamma(id)*Nu + zeta(id)*Nu_old - alpha(id)/(Re)*(Lux+Luy+Lubc)); -xu = getIntermediateU_xyRK3(u, b, dt, Re, nx, ny, dx, dy, id); +% xu = getIntermediateU_xyRK3(u, b, dt, Re, nx, ny, dx, dy, id); +xu = getIntermediateU_xyRK3_dst(u, b, dt, Re, nx, ny, dx, dy, id); b = v(2:end-1,2:end-1) - dt*(gamma(id)*Nv + zeta(id)*Nv_old - alpha(id)/(Re)*(Lvx+Lvy+Lvbc)); -xv = getIntermediateV_xyRK3(v, b, dt, Re, nx, ny, dx, dy, id); +% xv = getIntermediateV_xyRK3(v, b, dt, Re, nx, ny, dx, dy, id); +xv = getIntermediateV_xyRK3_dst(v, b, dt, Re, nx, ny, dx, dy, id); u(2:end-1,2:end-1) = xu; v(2:end-1,2:end-1) = xv;