diff --git a/README.md b/README.md index c32df78..b633214 100644 --- a/README.md +++ b/README.md @@ -4,23 +4,34 @@ This repo provides a MATLAB example code for the lid-driven cavity flow where incompressible Navier Stokes equation is numerically solved using a simple 2nd order finite difference scheme on a staggered grid system. +![sample](./gif/animation_sample1e2and1e4.gif) +(Left: Re = 100, Right: Re = 10,000) -![sample](./gif/animation_sample.gif) -The top lid changes its direction. The arrow denotes the velocity field and the contour denotes its magnitude. +The top lid changes its direction. The arrow denotes the velocity field, and the contour denotes its magnitude. ## Part 1 + - Click [here](./docs_part1/vanilaCavityFlow_EN.md) for detailed documentation in English. - 日本語のドキュメントは[こちら](./docs_part1/vanilaCavityFlow_JP.md) から -In order to focus on understanding basic idea of the numerical schemes, the method is kept premitive; explicit treatment of viscous term (the solution will diverge at low Reynolds number) and -the time integration is Euler. +The numerical scheme is kept primitive; the explicit treatment of viscous term (the solution diverges at low Reynolds number), and the time integration is Euler. まずは単純な手法でキャビティ流れのシミュレーションを実施します。 +## Part 2 + +- Click [here](./docs_part2/vanilaCavityFlowImplicit_EN.md) for detailed documentation in English. +- 日本語のドキュメントは[こちら](./docs_part2/vanilaCavityFlowImplicit_JP.md) から + +The implicit treatments for viscous terms are implemented at low Reynolds number, namely the Crank-Nicolson method. For better stability for non-linear terms, Adams-Bashforth, and 3 steps-Runge-Kutta is also implemented. + +拡散項に対して陰解法を実装しました。対流項へアダムス・バッシュフォースを使用したもの、3段階のルンゲクッタ法の２つの時間発展を実装しています。 + + ## Next to come -Plan is to implement the implicit treatment for viscous terms to obtain better stability. +The plan is to allow arbitrary boundary conditions for more fun simulations. @@ -33,7 +44,7 @@ Plan is to implement the implicit treatment for viscous terms to obtain better s # ToDo 1. Implement implicit treatment of viscous terms -2. Imlement crank-Nicolson for the non-linear terms +2. Implement crank-Nicolson for the non-linear terms 3. Allow obstacles within the domain 4. Allow inflow from the wall 5. Make it to 3D diff --git a/docs_part2/vanilaCavityFlowImplicit_EN.mlx b/docs_part2/vanilaCavityFlowImplicit_EN.mlx new file mode 100644 index 0000000..a187a2c Binary files /dev/null and b/docs_part2/vanilaCavityFlowImplicit_EN.mlx differ diff --git a/functions/getIntermediateU_xyCNAB.m b/functions/getIntermediateU_xyCNAB.m new file mode 100644 index 0000000..242dbd9 --- /dev/null +++ b/functions/getIntermediateU_xyCNAB.m @@ -0,0 +1,24 @@ +% Copyright 2020 The MathWorks, Inc. +function xu = getIntermediateU_xyCNAB(u, b, dt, Re, nx, ny, dx, dy) + +persistent A4u dtRe_old L4u + +if isempty(dtRe_old) + dtRe_old = dt/Re; +end + +if isempty(A4u) || any(size(L4u) ~= [(nx-1)*ny,(nx-1)*ny]) || dtRe_old ~= dt/Re + dtRe_old = dt/Re; + maskU = false(nx+1,ny+2); + maskU(2:end-1,2:end-1) = true; + L4u = getL4u(nx,ny,dx,dy,maskU); + A4u = eye(size(L4u),'like',L4u)-dt/(2*Re)*L4u; + A4u = decomposition(A4u); +end + +% x0 = u(2:end-1,2:end-1); +% [xu,flag] = cgs(A4u,b(:),[],[],[],[],x0(:)); +xu = A4u\b(:); +xu = reshape(xu,[nx-1,ny]); + +end \ No newline at end of file diff --git a/functions/getIntermediateU_xyRK3.m b/functions/getIntermediateU_xyRK3.m new file mode 100644 index 0000000..bce3e88 --- /dev/null +++ b/functions/getIntermediateU_xyRK3.m @@ -0,0 +1,37 @@ +% Copyright 2020 The MathWorks, Inc. +function xu = getIntermediateU_xyRK3(u, b, dt, Re, nx, ny, dx, dy, id) + +persistent A4u1 A4u2 A4u3 dtRe_old L4u + +if isempty(dtRe_old) + dtRe_old = dt/Re; +end + +if isempty(A4u1) || any(size(L4u) ~= [(nx-1)*ny,(nx-1)*ny]) || dtRe_old ~= dt/Re + dtRe_old = dt/Re; + maskU = false(nx+1,ny+2); + maskU(2:end-1,2:end-1) = true; + L4u = getL4u(nx,ny,dx,dy,maskU); + A4u1 = eye(size(L4u),'like',L4u)-(8/15)*dt/(Re)*L4u; + A4u2 = eye(size(L4u),'like',L4u)-(2/15)*dt/(Re)*L4u; + A4u3 = eye(size(L4u),'like',L4u)-(2/6)*dt/(Re)*L4u; + A4u1 = decomposition(A4u1); + A4u2 = decomposition(A4u2); + A4u3 = decomposition(A4u3); + +end + +% x0 = u(2:end-1,2:end-1); +% [xu,flag] = cgs(A4u,b(:),[],[],[],[],x0(:)); +switch id + case 1 + xu = A4u1\b(:); + case 2 + xu = A4u2\b(:); + case 3 + xu = A4u3\b(:); +end + +xu = reshape(xu,[nx-1,ny]); + +end \ No newline at end of file diff --git a/functions/getIntermediateV_xyCNAB.m b/functions/getIntermediateV_xyCNAB.m new file mode 100644 index 0000000..f9babe5 --- /dev/null +++ b/functions/getIntermediateV_xyCNAB.m @@ -0,0 +1,25 @@ +% Copyright 2020 The MathWorks, Inc. +function xv = getIntermediateV_xyCNAB(v, b, dt, Re, nx, ny, dx, dy) + +persistent A4v dtRe_old L4v + +if isempty(dtRe_old) + dtRe_old = dt/Re; +end + +if isempty(A4v) || any(size(L4v) ~= [nx*(ny-1),nx*(ny-1)]) || dtRe_old ~= dt/Re + dtRe_old = dt/Re; + maskV = false(nx+2,ny+1); + maskV(2:end-1,2:end-1) = true; + L4v = getL4v(nx,ny,dx,dy,maskV); + A4v = eye(size(L4v),'like',L4v)-dt/(2*Re)*L4v; + A4v = decomposition(A4v); +end + +% x0 = v(2:end-1,2:end-1); +% [xv,flag] = cgs(A4v,b(:),[],[],[],[],x0(:)); +xv = A4v\b(:); +xv = reshape(xv,[nx,ny-1]); + +end + diff --git a/functions/getIntermediateV_xyRK3.m b/functions/getIntermediateV_xyRK3.m new file mode 100644 index 0000000..ee89d79 --- /dev/null +++ b/functions/getIntermediateV_xyRK3.m @@ -0,0 +1,35 @@ +% Copyright 2020 The MathWorks, Inc. +function xv = getIntermediateV_xyRK3(v, b, dt, Re, nx, ny, dx, dy,id) + +persistent A4v1 A4v2 A4v3 dtRe_old L4v + +if isempty(dtRe_old) + dtRe_old = dt/Re; +end + +if isempty(A4v1) || any(size(L4v) ~= [nx*(ny-1),nx*(ny-1)]) || dtRe_old ~= dt/Re + dtRe_old = dt/Re; + maskV = false(nx+2,ny+1); + maskV(2:end-1,2:end-1) = true; + L4v = getL4v(nx,ny,dx,dy,maskV); + A4v1 = eye(size(L4v),'like',L4v)-(8/15)*dt/(Re)*L4v; + A4v2 = eye(size(L4v),'like',L4v)-(2/15)*dt/(Re)*L4v; + A4v3 = eye(size(L4v),'like',L4v)-(2/6)*dt/(Re)*L4v; + A4v1 = decomposition(A4v1); + A4v2 = decomposition(A4v2); + A4v3 = decomposition(A4v3); + +end + +switch id + case 1 + xv = A4v1\b(:); + case 2 + xv = A4v2\b(:); + case 3 + xv = A4v3\b(:); +end +xv = reshape(xv,[nx,ny-1]); + +end + diff --git a/functions/getL4u.m b/functions/getL4u.m new file mode 100644 index 0000000..021a357 --- /dev/null +++ b/functions/getL4u.m @@ -0,0 +1,64 @@ +% Copyright 2020 The MathWorks, Inc. +function L4u = getL4u(nx,ny,dx,dy,mask) + +% Note: Boundary condition effects y derivative only +% y-derivatives +% ub = (u0+u1)/2 => u0 = 2*ub - u1; +% Lu_1 = (u2 - 2*u1 + u0)/dy^2 +% = (u2 - 2*u1 + 2*ub - u1)/dy^2 +% = (u2 - 3*u1)/dx^2 + 2*ub/dy^2 +% +% x-derivatives +% Lu_1 = (u2 - 2*u1 + ub)/dx^2 +% = (u2 - 2*u1)/dx^2 + ub/dx^2 + +matrixSize = [nx-1,ny]; + +% 3 +%152 +% 4 +coefx = 1/dx^2; +coefy = 1/dy^2; + +i = zeros(nx*ny*5,1); +j = zeros(nx*ny*5,1); +v = zeros(nx*ny*5,1); +index1 = 1; +index = 1; +for jj=1:ny + for ii=1:nx-1 + idx = [mask(ii+2,jj+1),mask(ii,jj+1),mask(ii+1,jj+2),mask(ii+1,jj),mask(ii+1,jj+1)]; + stencils = [ii+1,jj + ii-1,jj + ii,jj+1 + ii,jj-1 + ii,jj]; + stencils = stencils(idx,:); + + tmpx = 2; + tmpy = 2 + sum(~idx(3:4)); + + coeffs = [coefx + coefx + coefy + coefy + -tmpx*coefx-tmpy*coefy]; + coeffs = coeffs(idx); + linearIdx = sub2ind(matrixSize, stencils(:,1), stencils(:,2)); + + i(index1:index1+length(linearIdx)-1) = index; + j(index1:index1+length(linearIdx)-1) = linearIdx; + v(index1:index1+length(linearIdx)-1) = coeffs; + index1 = index1 + length(linearIdx); + index = index + 1; + end +end + +i(i==0) = []; +j(j==0) = []; +v(v==0) = []; +L4u = sparse(i,j,v,(nx-1)*ny,(nx-1)*ny); +L4u(:,~any(L4u,1)) = []; +L4u(~any(L4u,2),:) = []; + +end diff --git a/functions/getL4v.m b/functions/getL4v.m new file mode 100644 index 0000000..e67b264 --- /dev/null +++ b/functions/getL4v.m @@ -0,0 +1,64 @@ +% Copyright 2020 The MathWorks, Inc. +function L4v = getL4v(nx,ny,dx,dy,mask) + +% Note: Boundary condition effects x derivative only +% x-derivatives +% vb = (v0+v1)/2 => v0 = 2*vb - v1; +% Lv_1 = (v2 - 2*v1 + v0)/dx^2 +% = (v2 - 2*v1 + 2*vb - u1)/dx^2 +% = (v2 - 3*v1)/dx^2 + 2*vb/dx^2 +% +% y-derivatives +% Lv_1 = (v2 - 2*v1 + vb)/dy^2 +% = (v2 - 2*v1)/dy^2 + vb/dy^2 + +matrixSize = [nx,ny-1]; + +% 3 +%152 +% 4 +coefx = 1/dx^2; +coefy = 1/dy^2; + +i = zeros(nx*ny*5,1); +j = zeros(nx*ny*5,1); +v = zeros(nx*ny*5,1); +index1 = 1; +index = 1; +for jj=1:ny-1 + for ii=1:nx + idx = [mask(ii+2,jj+1),mask(ii,jj+1),mask(ii+1,jj+2),mask(ii+1,jj),mask(ii+1,jj+1)]; + stencils = [ii+1,jj + ii-1,jj + ii,jj+1 + ii,jj-1 + ii,jj]; + stencils = stencils(idx,:); + + tmpx = 2 + sum(~idx(1:2)); + tmpy = 2; + + coeffs = [coefx + coefx + coefy + coefy + -tmpx*coefx-tmpy*coefy]; + coeffs = coeffs(idx); + linearIdx = sub2ind(matrixSize, stencils(:,1), stencils(:,2)); + + i(index1:index1+length(linearIdx)-1) = index; + j(index1:index1+length(linearIdx)-1) = linearIdx; + v(index1:index1+length(linearIdx)-1) = coeffs; + index1 = index1 + length(linearIdx); + index = index + 1; + end +end + +i(i==0) = []; +j(j==0) = []; +v(v==0) = []; +L4v = sparse(i,j,v,nx*(ny-1),nx*(ny-1)); +L4v(:,~any(L4v,1)) = []; +L4v(~any(L4v,2),:) = []; + +end \ No newline at end of file diff --git a/functions/operatorAu_CNAB.m b/functions/operatorAu_CNAB.m new file mode 100644 index 0000000..d4aee75 --- /dev/null +++ b/functions/operatorAu_CNAB.m @@ -0,0 +1,19 @@ +% Copyright 2020 The MathWorks, Inc. +function Ax = operatorAu_CNAB(x,dt,Re,nx,ny,dx,dy) +% Getting A for u velocity (Crank-Nicolson) +% A = [I - dt/(2*Re)*Lxy]; + +x = reshape(x,nx-1,ny); +xbig = zeros(nx+1,ny+2); +xbig(2:end-1,2:end-1) = x; + +% Addting boudanry value +% ub = (u0+u1)/2 => u0 = 2*ub - u1; +xbig(:,1) = - xbig(:,2); +xbig(:,end) = - xbig(:,end-1); + +Ax = x - dt/(2*Re)*(xbig(1:end-2,2:end-1)-2*xbig(2:end-1,2:end-1)+xbig(3:end,2:end-1))/dx^2; % nx-1 * ny +Ax = Ax - dt/(2*Re)*(xbig(2:end-1,1:end-2)-2*xbig(2:end-1,2:end-1)+xbig(2:end-1,3:end))/dy^2; % nx-1 * ny +Ax = Ax(:); + +end diff --git a/functions/operatorAv_CNAB.m b/functions/operatorAv_CNAB.m new file mode 100644 index 0000000..eab8583 --- /dev/null +++ b/functions/operatorAv_CNAB.m @@ -0,0 +1,19 @@ +% Copyright 2020 The MathWorks, Inc. +function Ax = operatorAv_CNAB(x,dt,Re,nx,ny,dx,dy) +% Getting A for v velocity (Crank-Nicolson) +% A = [I - dt/(2*Re)*Lxy]; + +x = reshape(x,nx,ny-1); +xbig = zeros(nx+2,ny+1); +xbig(2:end-1,2:end-1) = x; + +% Addting boudanry value +% ub = (u0+u1)/2 => u0 = 2*ub - u1; +xbig(1,:) = - xbig(2,:); +xbig(end,:) = - xbig(end-1,:); + +Ax = x - dt/(2*Re)*(xbig(1:end-2,2:end-1)-2*xbig(2:end-1,2:end-1)+xbig(3:end,2:end-1))/dx^2; % nx-1 * ny +Ax = Ax - dt/(2*Re)*(xbig(2:end-1,1:end-2)-2*xbig(2:end-1,2:end-1)+xbig(2:end-1,3:end))/dy^2; % nx-1 * ny +Ax = Ax(:); +end + diff --git a/functions/updateVelocityField_CNAB.m b/functions/updateVelocityField_CNAB.m new file mode 100644 index 0000000..1892368 --- /dev/null +++ b/functions/updateVelocityField_CNAB.m @@ -0,0 +1,95 @@ +function [u,v] = updateVelocityField_CNAB(u,v,nx,ny,dx,dy,Re,dt,bctop,method) +% Copyright 2020 The MathWorks, Inc. + +persistent Nu_old Nv_old + +if isempty(Nu_old) || isempty(Nv_old) || any(size(Nu_old) ~= [(nx-1)*ny,(nx-1)*ny]) + Nu_old = zeros(nx-1,ny); + Nv_old = zeros(nx,ny-1); +end + +% Apply boundary conditions: +% represented as the values on ghost cells +u(:,1) = 0-u(:,2); v(:,1) = 0.; %bottom +u(:,end) = 2*bctop-u(:,end-1); v(:,end) = 0; %top +u(1,:) = 0; v(1,:) = 0-v(2,:); %left +u(end,:) = 0; v(end,:) = 0-v(end-1,:); %right + +% gU(u) Get viscous terms for u +Lux = (u(1:end-2,2:end-1)-2*u(2:end-1,2:end-1)+u(3:end,2:end-1))/dx^2; % nx-1 * ny +Luy = (u(2:end-1,1:end-2)-2*u(2:end-1,2:end-1)+u(2:end-1,3:end))/dy^2; % nx-1 * ny + +% gU(v) Get viscous terms for v +Lvx = (v(1:end-2,2:end-1)-2*v(2:end-1,2:end-1)+v(3:end,2:end-1))/dx^2; % nx * ny-1 +Lvy = (v(2:end-1,1:end-2)-2*v(2:end-1,2:end-1)+v(2:end-1,3:end))/dy^2; % nx * ny-1 + +% Η̌vZ +% Get nonlinear terms +% 1. interpolate velocity at cell center/cell cornder +uce = (u(1:end-1,2:end-1)+u(2:end,2:end-1))/2; +uco = (u(:,1:end-1)+u(:,2:end))/2; +vco = (v(1:end-1,:)+v(2:end,:))/2; +vce = (v(2:end-1,1:end-1)+v(2:end-1,2:end))/2; + +% 2. multiply +uuce = uce.*uce; +uvco = uco.*vco; +vvce = vce.*vce; + +% 3-1. get derivative for u +Nu = (uuce(2:end,:) - uuce(1:end-1,:))/dx; +Nu = Nu + (uvco(2:end-1,2:end) - uvco(2:end-1,1:end-1))/dy; + +% 3-2. get derivative for v +Nv = (vvce(:,2:end) - vvce(:,1:end-1))/dy; +Nv = Nv + (uvco(2:end,2:end-1) - uvco(1:end-1,2:end-1))/dx; + +% ̑xZo +% Implicit treatment for xy direction +Lubc = zeros(size(Luy)); +Lubc(:,end) = 2*bctop/dy^2; % effect of the top BC on Lu +Lvbc = zeros(size(Lvy)); + +b = u(2:end-1,2:end-1) - dt*((3*Nu-Nu_old)/2 - 1/(2*Re)*(Lux+Luy+Lubc)); +xu = getIntermediateU_xyCNAB(u, b, dt, Re, nx, ny, dx, dy); +b = v(2:end-1,2:end-1) - dt*((3*Nv-Nv_old)/2 - 1/(2*Re)*(Lvx+Lvy+Lvbc)); +xv = getIntermediateV_xyCNAB(v, b, dt, Re, nx, ny, dx, dy); + +Nu_old = Nu; +Nv_old = Nv; + +u(2:end-1,2:end-1) = xu; +v(2:end-1,2:end-1) = xv; + +% Vx +% ͂̎i|\jđxʕۑ𖞂ɎʑB +% velocity correction +% RHS of pressure Poisson eq. +b = ((u(2:end,2:end-1)-u(1:end-1,2:end-1))/dx ... + + (v(2:end-1,2:end)-v(2:end-1,1:end-1))/dy); + +% Solve for p +switch method + case 'direct' + % by directly inverting the matrixiږ@j + dp = solvePoissonEquation_direct(b,nx,ny,dx,dy); + case 'dct' + % by using the discrete cosine transformiRTCϊgpj + % Note: Signal Processing Toolbox required + dp = solvePoissonEquation_2dDCT(b,nx,ny,dx,dy); + + otherwise + error("Specified method: " + method + " is not supported." + ... + "It should be either direct or dct"); +end +% correction to get the final velocity +p = dp; +u(2:end-1,2:end-1) = u(2:end-1,2:end-1) - (p(2:end,:)-p(1:end-1,:))/dx; +v(2:end-1,2:end-1) = v(2:end-1,2:end-1) - (p(:,2:end)-p(:,1:end-1))/dy; + +% check the divergence +% b = ((u(2:end,2:end-1)-u(1:end-1,2:end-1))/dx ... +% + (v(2:end-1,2:end)-v(2:end-1,1:end-1))/dy); +% norm(b) + +end \ No newline at end of file diff --git a/functions/updateVelocityField_RK3.m b/functions/updateVelocityField_RK3.m new file mode 100644 index 0000000..ab3e857 --- /dev/null +++ b/functions/updateVelocityField_RK3.m @@ -0,0 +1,117 @@ +% Copyright 2020 The MathWorks, Inc. +function [u,v] = updateVelocityField_RK3(u,v,nx,ny,dx,dy,Re,dt,bctop,method) + +[u,v] = updateVelocityField_RK3substep(u,v,nx,ny,dx,dy,Re,dt,bctop,method,1); +[u,v] = updateVelocityField_RK3substep(u,v,nx,ny,dx,dy,Re,dt,bctop,method,2); +[u,v] = updateVelocityField_RK3substep(u,v,nx,ny,dx,dy,Re,dt,bctop,method,3); + +% check the divergence +% b = ((u(2:end,2:end-1)-u(1:end-1,2:end-1))/dx ... +% + (v(2:end-1,2:end)-v(2:end-1,1:end-1))/dy); +% norm(b) + +end + +function [u,v] = updateVelocityField_RK3substep(u,v,nx,ny,dx,dy,Re,dt,bctop,method,id) + +persistent Nu_old Nv_old + +if isempty(Nu_old) || isempty(Nv_old) || any(size(Nu_old) ~= [(nx-1)*ny,(nx-1)*ny]) + Nu_old = zeros(nx-1,ny); + Nv_old = zeros(nx,ny-1); +end + +% Apply boundary conditions: +% represented as the values on ghost cells +u(:,1) = 0-u(:,2); v(:,1) = 0.; %bottom +u(:,end) = 2*bctop-u(:,end-1); v(:,end) = 0.; %top +u(1,:) = 0; v(1,:) = 0-v(2,:); %left +u(end,:) = 0; v(end,:) = 0-v(end-1,:); %right + +% gU(u) Get viscous terms for u +Lux = (u(1:end-2,2:end-1)-2*u(2:end-1,2:end-1)+u(3:end,2:end-1))/dx^2; % nx-1 * ny +Luy = (u(2:end-1,1:end-2)-2*u(2:end-1,2:end-1)+u(2:end-1,3:end))/dy^2; % nx-1 * ny + +% gU(v) Get viscous terms for v +Lvx = (v(1:end-2,2:end-1)-2*v(2:end-1,2:end-1)+v(3:end,2:end-1))/dx^2; % nx * ny-1 +Lvy = (v(2:end-1,1:end-2)-2*v(2:end-1,2:end-1)+v(2:end-1,3:end))/dy^2; % nx * ny-1 + +% Η̌vZ +% Get nonlinear terms +% 1. interpolate velocity at cell center/cell cornder +uce = (u(1:end-1,2:end-1)+u(2:end,2:end-1))/2; +uco = (u(:,1:end-1)+u(:,2:end))/2; +vco = (v(1:end-1,:)+v(2:end,:))/2; +vce = (v(2:end-1,1:end-1)+v(2:end-1,2:end))/2; + +% 2. multiply +uuce = uce.*uce; +uvco = uco.*vco; +vvce = vce.*vce; + +% 3-1. get derivative for u +Nu = (uuce(2:end,:) - uuce(1:end-1,:))/dx; +Nu = Nu + (uvco(2:end-1,2:end) - uvco(2:end-1,1:end-1))/dy; + +% 3-2. get derivative for v +Nv = (vvce(:,2:end) - vvce(:,1:end-1))/dy; +Nv = Nv + (uvco(2:end,2:end-1) - uvco(1:end-1,2:end-1))/dx; + +% ̑xZo +Lubc = zeros(size(Luy)); +% Lbc(1,:) = u(1,2:end-1)/dx^2; +% Lbc(end,:) = u(end,2:end-1)/dx^2; +% Lbc(:,1) = bcbottom; % 0 in this case +Lubc(:,end) = 2*bctop/dy^2; + +Lvbc = zeros(size(Lvy)); +% Lbc(1,:) = 2*v(1,2:end-1)/dx^2; +% Lbc(end,:) = 2*v(end,2:end-1)/dx^2; +% Lbc(:,1) = v(2:end-1,1)/dy^2; +% Lbc(:,end) = v(2:end-1,end)/dy^2; + +alpha = [4/15,1/15,1/6]; +beta = alpha; +gamma = [8/15,5/12,3/4]; +zeta = [0,-17/60,-5/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); + +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); + +Nu_old = Nu; +Nv_old = Nv; + + +u(2:end-1,2:end-1) = xu; +v(2:end-1,2:end-1) = xv; + +% Vx +% ͂̎i|\jđxʕۑ𖞂ɎʑB +% velocity correction +% RHS of pressure Poisson eq. +b = ((u(2:end,2:end-1)-u(1:end-1,2:end-1))/dx ... + + (v(2:end-1,2:end)-v(2:end-1,1:end-1))/dy); + +% Solve for p +switch method + case 'direct' + % by directly inverting the matrixiږ@j + dp = solvePoissonEquation_direct(b,nx,ny,dx,dy); + case 'dct' + % by using the discrete cosine transformiRTCϊgpj + % Note: Signal Processing Toolbox required + dp = solvePoissonEquation_2dDCT(b,nx,ny,dx,dy); + otherwise + error("Specified method: " + method + " is not supported." + ... + "It should be either direct or dct"); +end +% correction to get the final velocity +p = dp; +u(2:end-1,2:end-1) = u(2:end-1,2:end-1) - (p(2:end,:)-p(1:end-1,:))/dx; +v(2:end-1,2:end-1) = v(2:end-1,2:end-1) - (p(:,2:end)-p(:,1:end-1))/dy; + +end \ No newline at end of file diff --git a/gif/animation_sample1e2and1e4.gif b/gif/animation_sample1e2and1e4.gif new file mode 100644 index 0000000..8c4fffd Binary files /dev/null and b/gif/animation_sample1e2and1e4.gif differ