Copyright (c) 2020, The MathWorks, Inc.
In this document, we describe the implicit treatment of the viscous terms in the Navier-Stokes equation for better stability at low Reynolds numbers. At the same time, two steps Adams-Bashforth and three steps Runge-Kutta are described and implemented for better time-integration stability and accuracy.
This document is part 2 of the CFD 101 series.
- Understand how to implement Crank-Nicolson for viscous terms
- Understand how to appy Adams-Bashforth for convective terms.
- Simulate at high Reynolds number (Re = 5,000 below)
The discussion below is based on Adams-Bashforth for simplicity, but more robust (but expensive) Runge-Kutta is also implemented. See Appendix.
- MATLAB R2019b
- Signal Processing Toolbox*
*: It's recommended. The code below uses dct function from Signal Processing Toolbox to solve the Poisson equation for pressure, but this part can be replaced by iterative solver or directly inverting the matrices, which does not require any additional toolboxes. Examples are also provided.
Popular Neumann's stability analysis and Courant number discussions. The point is the following:
The stability condition of the diffusion equation with a central difference and explicit Euler integration is
and we need 1/4 of the time step if we 1/2 the grid size. Also, the condition is severe for lower Reynolds numbers, which explains why the solution in part 1 diverges at a low Reynolds number. It is known that implicit treatment provides unconditional stability.
For the convection equations, the Courant number (c is the velocity）
comes into play and we need to keep the Courant number not exceeding which depends on the method used. Based on the observation we use
for Adams-Bashforth.
The 1st order explicit Euler integration gives
and
in a matrix form.
The implicit Crank-Nicolson scheme is used for linear (viscous) terms and second-order Adams-Bashforth scheme for non-linear terms. The discrete form of the above equation can be written as
Or in the matrix form, as
where
where is the implicit operator for the advection-diffusion part of the momentum equation,
is the convective operator,
is the gradient operator,
is the divergence operator,
is the explicit right-hand side of the momentum equation,
is the boundary condition vector from the viscous term and
is the boundary condition vector for the incompressibility constraint. Only the viscous term is treated implicitly here and has the velocity field of the next time step (
).
The fractional step method is related to the block LU factorization of the equation above (Perot 1993 [1])
This is exact but computationally very expensive since calculating the inverse of is not practical. Hence it is usually solved approximating
. Different approximations to the inverse result in different classes of fractional step method.
The classic fractional step method corresponds to using , which results in a first-order error term. By choosing
, the resulting error is 2nd order. In the case of
,
The single time step reduces to the following sequence of operations,
The difference from part 1 is the first equation, namely
and we need to solver the system of equation for both
and
.
It should be noted that the operator acts only on the velocity inside the domain excluing the boundary conditions, thus we have
on the right hand side.
Enough equations already. Let's put everything in MATLAB. Please refer to each function on the repository for details, but we would like to highlight the following three points.
- The matrices that we need to solve for viscous term
- How to deal with the boundary condition
- How to solve the equations
Again, the sole difference from the previous Euler method is the following equation to solve for the intermediate velocity.
In the MATLAB code (see updateVelocityField_CNAB.m for details)
% Implicit treatment for xy direction
b = u(2:end-1,2:end-1) + dt*(-3*Nu+Nu_old)/2 + dt/(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 + dt/(2*Re)*(Lvx+Lvy+Lvbc);
xv = getIntermediateV_xyCNAB(v, b, dt, Re, Nx, Ny, dx, dy);
u(2:end-1,2:end-1) = xu;
v(2:end-1,2:end-1) = xv;
Euler method (see updateVelocityField_Euler.m for details)
u(2:end-1,2:end-1) = u(2:end-1,2:end-1) - dt*Nu + dt/Re*(Lux+Luy);
v(2:end-1,2:end-1) = v(2:end-1,2:end-1) - dt*Nv + dt/Re*(Lvx+Lvy);
So simple!
Anyways, you can find the essense in getIntermediateU_xyCNAB.m and getIntermediateV_xyCNAB.m. CNAB denotes Crank-Nicolson, Adams-Bashforth
Let us revisit the that you can find in
. The difference equation is
and is only tri-diagonal matrix if you decides to apply implicit method for x-direction only. With the implicit treatment for both x,y direction,
becomes a little wider, but sparse matrix anyways. The size of the matrix is
for
for
and the number of element is but the number of non-zero element is
. The inconsistency in size for
and
stem from the use of the staggered-grid system.
For the boundary layer flow, for example, it was shown in Simens (2008), that the viscous time-step requirement could be more severe than the convective restriction only in the wall-normal direction in the case of a stretched grid that is usually employed in wall-bounded turbulence to resolve small-scale phenomena close to the wall. In such a case, it makes sense to apply an implicit method for one direction, and the equation can be solved more efficiently.
Here we review where the term comes from. The staggered-grid cell configuration and its indexes are sketched below.
The velocity in orange is the "ghost velocity" and calculated as (for example)
The idea is that the mean (2nd order interpolation) of the two velocities across the boundary becomes the boundary value. Since the ghost velocity has the value of the inside the domain, how we treat the boundary affects the definition of .
For example, at the bottom boundary, the y-component of the viscous term for becomes
And the coefficients for is now -3 and the second term of the right hand side goes to
. On the contrary, the x-component of the viscous term for
at the left boudnary becomes
The coeffients for is kept being -2. The ghost velocity is not required and
directly refers to the boundary value.
The situation is oppisite for the , and
requires ghost velocity for the x-component of the viscous terms.
For small number of grid sizes, let getL4u.m calculate the for
.
clear, close all
addpath('../functions/')
Lx = 1; Ly = 1;
nx = 5; ny = 5;
dx = Lx/nx; dy = Ly/ny;
maskU = false(nx+1,ny+2);
maskU(2:end-1,2:end-1) = true;
L4u = getL4u(nx,ny,dx,dy,maskU);
spy(L4u);
Sparse it is. The actual values are..
full(L4u(1:nx,1:nx))
ans = 5x5
-125.0000 25.0000 0 0 25.0000
25.0000 -125.0000 25.0000 0 0
0 25.0000 -125.0000 25.0000 0
0 0 25.0000 -125.0000 0
25.0000 0 0 0 -100.0000
Sidenote: maskU is the logical array with true where is defined.
getL4u.m calculates for more complex geometry as long as we defined
maskU properly. For the current case, we have trues simply surrounded by falses (boundary)
maskU
maskU = 6x7 の logical 配列
0 0 0 0 0 0 0
0 1 1 1 1 1 0
0 1 1 1 1 1 0
0 1 1 1 1 1 0
0 1 1 1 1 1 0
0 0 0 0 0 0 0
MATLAB is indeed good at solving and there exist so many build-in functions. We have implemented three. For the actual implementation, please refer to
getIntermediateU_xyCNAB.m.
As observed above is a sparse matrix.
dt = 0.01; Re = 100; % a random test case
A4u = eye(size(L4u),'like',L4u)-dt/(2*Re)*L4u; % A matrix for u
spy(A4u)
MATLAB backslash works nicely for sparse matrix too.
r = rand(nx-1,ny);
xu = A4u\r(:);
If we have fixed time-step size dt and Reynolds number Re, we can pre-decompose the and the backslash operation becomes much faster by reuse the matrix decompositions, though it requires more memory to store them.
decomposition Object uses the decomposition type is automatically chosen based on the properties of the input matrix.
A4u = decomposition(A4u);
xu1 = A4u\r(:);
For this method to work, we need to make use of the persistent variable.
MATLAB has 11 algorithms ready.
Here we choose cgs since it does not require symmetry nor positive definite.
A4u = eye(size(L4u),'like',L4u)-dt/(2*Re)*L4u;
xu2 = cgs(A4u,r(:));
cgs は、相対残差 6.5e-11 をもつ解に 反復 2 で収束しました。
The use of preconditioner helps convergence, thus the computational time, but still, it requires the constant time step size dt and Reynolds number.
Sidenote: The useful nature of the iterative method is that we do not have to explicitly prepare the matrices.
xu4 = cgs(@(x) operatorAu_CNAB(x,dt,Re,nx,ny,dx,dy),r(:));
cgs は、相対残差 6.5e-11 をもつ解に 反復 2 で収束しました。
You can specify the first input argument as a function handle such that the function returns (here
operatorAu_CNAB.m)
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(:);
The difference equation can be directly written in code.
Let's animate the flow with Re = 5,000. The numerical integration process discussed above is now in the function updateVelocityField_CNAB.m.
NOTE: Somehow, the code gives an error with recordGIF = true if you are running in R2019b or earlier. If you want to record the animation to GIF, please run script_AnimateVelocityField_part2.m instead.
Re = 5000; % Reynolds number
nt = 3000; % max time steps
Lx = 1; Ly = 1; % domain size
Nx = 80; Ny = 80; % Number of grids
dt = 0.003; % time step;
% Grid size (Equispaced)
dx = Lx/Nx;
dy = Ly/Ny;
% Coordinate of each grid (cell center)
xce = ((1:Nx)-0.5)*dx;
yce = ((1:Ny)-0.5)*dy;
% Coordinate of each grid (cell corner)
xco = (0:Nx)*dx;
yco = (0:Ny)*dy;
Initialize the flow field.
u = zeros(Nx+1,Ny+2); % velocity in x direction (u)
v = zeros(Nx+2,Ny+1); % velocity in y direction (v)
uce = (u(1:end-1,2:end-1)+u(2:end,2:end-1))/2; % u at cell center
vce = (v(2:end-1,1:end-1)+v(2:end-1,2:end))/2; % v at cell center
Downsample the velocity field at the interval of visRate so that the plot is not going to be filled with arrows. Also GIF is updated (appended) at every recordRate time steps.
visRate = 4; % downsample rate of the data for quiver
recordGIF = false; % set to true if you wanna make GIF
recordRate = 20;
filename = 'animation_sample.gif'; % Specify the output file name
Contour plot
figure
[Xce,Yce] = meshgrid(xce,yce); % cell center grid
[~,h_abs] = contourf(Xce',Yce',sqrt(uce.^2+vce.^2)); % contour
警告: 等高線図は ZData が定数の場合はレンダリングされません
hold on
Quiver plot
% Downsample the data（d = downsampled）
xced = xce(1:visRate:end);
yced = yce(1:visRate:end);
[Xced,Yced] = meshgrid(xced, yced);
uced = uce(1:visRate:end,1:visRate:end);
vced = vce(1:visRate:end,1:visRate:end);
h_quiver = quiver(Xced',Yced',uced,vced,3,'Color',[1,1,1]);
hold off
xlim([0 Lx]); ylim([0 Ly]);
Additional arrow representing the bounary condition at the top lid.
harrow = annotation('textarrow',[0.3 0.7],[0.96 0.96],"LineWidth",2);
Delete the ticks/tick labels.
haxes = gca;
haxes.XTick = [];
haxes.YTick = [];
Just to make it a little fun, the velocity of the top lid reverses at the 1000 steps (out of 3000 steps).
initialFrame = true;
for ii = 1:nt
bctop = 1; % top velocity
if ii > 1000
bctop = -1;
harrow.X = [0.7, 0.3]; % inverse the arrow
end
% Update the velocity field (uses dct)
[u,v] = updateVelocityField_CNAB_bctop(u,v,Nx,Ny,dx,dy,Re,dt,bctop,'dct');
% [u,v] = updateVelocityField_RK3_bctop(u,v,Nx,Ny,dx,dy,Re,dt,bctop,'dct');
% Update the plot at every recordRate steps
if mod(ii,recordRate) == 0
% get velocity at the cell center (for visualization)
uce = (u(1:end-1,2:end-1)+u(2:end,2:end-1))/2; % u at cell center
vce = (v(2:end-1,1:end-1)+v(2:end-1,2:end))/2; % v at cell center
% update plot (downsample)
h_quiver.UData = uce(1:visRate:end,1:visRate:end);
h_quiver.VData = vce(1:visRate:end,1:visRate:end);
h_abs.ZData = sqrt(uce.^2+vce.^2);
drawnow
if recordGIF
frame = getframe(gcf); %#ok<UNRCH> % Figure 画面をムービーフレーム（構造体）としてキャプチャ
tmp = frame2im(frame); % 画像に変更
[A,map] = rgb2ind(tmp,256); % RGB -> インデックス画像に
if initialFrame
imwrite(A,map,filename,'gif','LoopCount',Inf,'DelayTime',0.1);
initialFrame = false;
else
imwrite(A,map,filename,'gif','WriteMode','append','DelayTime',0.1);% 画像をアペンド
end
end
end
end
We haven't yet validated the result or anything, but it seems to work fine with high Reynolds number.
In this document, we have discussed the Crank-Nicolson and Adams-Bashfoth. It's observed that the simulation still diverges with high Renolds number if the time step size is not small enough. If you find the simualtion is not working please try Runge-Kutta instead by replacing
[u,v] = updateVelocityField_CNAB_bctop(u,v,Nx,Ny,dx,dy,Re,dt,bctop,'dct');
by
[u,v] = updateVelocityField_RK3_bctop(u,v,Nx,Ny,dx,dy,Re,dt,bctop,'dct');
The low-storage third-order semi-implicit Runge–Kutta method of Spalart et al. (1991)[3] is also implemented for the temporal discretization to achieve a higher order of temporal accuracy. See updateVelocityField_RK3.m. The equations to be solved at each step then become
where
It required three times as much as a computation for a one-time step, but it allows a bigger time step size.
[1] B. Perot, An analysis of the fractional step method. J. Comp. Physics, 108: 51-58, 1993
[2] Simens, M.P., Jim´enez, J., Hoyas, S. and Mizuno, Y. 2009 A high-resolution code for turbulent boundary layers. Journal of Computational Physics 228 (11), 4218-4231.
[3] Spalart, P.R., Moser, R.D. and Rogers, M.M. 1991 Spectral methods for the Navier-Stokes equations with one infinite and two periodic directions. Journal of Computational Physics 96 (2), 297-324.