Description
Problem 1. Let us consider non-dimensional Navier-Stokes equations for incompressible flow in nondimensional form
| ∂tu + (u · ∇)u = −∇p + Re−1∇2u, | (1) |
| ∇ · u = 0, | (2) |
in a 2D periodic domain 0 ≤ x < 2π and 0 ≤ y < 2π. In order to illustrate the fractional step method, let us integrate the non-linear terms with an explicit Euler method and the linear terms including the pressure with an implicit Euler method, to arrive at
 ,                                         (3)
where A = I − Re−1∆tL, G is the discretized gradient operator, D is the discretized divergence operator, L = D · G is the discretized Laplacian operator and H are the discretized non-linear terms. To keep the illustration simple, let us consider a Fourier spatial discretization in both the x and y directions, with Nx = Ny = 32, Re = 1000 and ∆t = 0.01.
- Write a code to determine block matrices A, G, and D, in wavenumber space (α, β) and assemble matrix M. Plot each of these 4 matrices vs (α, β).
Hints: Plot matrices in 2D using the MATLAB command imagesc or a method that yields a similar output. Note that matrix A has a size 2n × 2n where n = NxNy, D has a size n × 2n and G = DT. You can define these matrices by means of for or do loops but you may also find the sparse, diag and repmat commands useful if you are working with MATLAB. For instance, the identity matrix term in the A block of M is generated by
%
II_2n = diag(ones(2*n,1));
% eye_in_A = sparse(1:2*n,1:2*n,II_2n,3*n,3*n);
% and the ∂x term in G can be generated by
% alp = [0:Nx/2-1 -Nx/2:-1];
% ialp = repmat(i*alp,[1,Ny]);
%
Galp = sparse(1:n     ,2*n+1:3*n,ialp,3*n,3*n); %
- The resulting matrix M is singular. In this Fourier representation, the singularity is easy to grasp since M has a row of zeros. Which one? (i.e. which value of (α, β)?) What does this mean physically in terms of the pressure?
- The artificial compressibility method removes the singularity of M by solving the approximate problem
 .                                                  (4)
Write a code to assemble the new M matrix of the artificial compressibility method for λ = 0.01. Plot the matrix and calculate the eigenvalue closest to zero.
- The fractional step method is equivalent to an approximate LU factorization of the original problem
(5)
| with B ≈ inv(A), leading to the scheme: | |
| Au∗ = un − Hn∆t | (6) |
,                                                                    (7)
un+1 = u∗ + BGpn+1∆t.                                                                         (8)
It turns out that working with a Fourier spatial discretization makes computing inv(A) trivial since A is a diagonal matrix. Thus, we will use an exact fractional step replacing B with inv(A) in the scheme above, and then use the solution as reference to compute the errors of the artificial compressibility scheme.
- Assuming that un = sinxsiny and vn = cosxcosy, which satisfies continuity, show that the non-linear terms are
| u∂xu + v∂yu = sin(2x)/2, | (9) |
| u∂xv + v∂yv = −sin(2y)/2. | (10) |
so that un −Hn∆t does not satisfy continuity and therefore u∗ does not satisfy continuity either. Fourier transform un and the non-linear terms and perform an exact iteration of the fractional step method. Assume pn = 0 everywhere. There is no aliasing in this case given that we chose Nx and Ny to be large enough. Plot the intermediate solution in the physical domain u∗(x,y), v∗(x,y) and p∗(x,y). Plot the solution in the physical domain un+1(x,y), vn+1(x,y) and pn+1(x,y). How does it differ from the intermediate solution? (Hint: Matrix DBG in 7 has one row of zeros for the same reason that M has one. But this indetermination is not important because the product BGpn+1 in 8 cancels it. So, it is ok to choose a bogus value of pbα=0,β=0.)
- Perform one iteration with the artificial compressibility scheme. Plot the solution in the physicaldomain un+1(x,y), vn+1(x,y) and pn+1(x,y).
- Plot the L2-norm of the errors in u, v and p as a function of ∆t and λ. Vary ∆t in the range [10−6,10−1] and λ in the range [10−4,0.1].
- Extra credit (kind of difficult): Incidentally, it also turns out that A being diagonal in Fourier discretizations also makes any matrix B produce an exact fractional step as long as B is diagonal, non-singular and B(n + j,n + j) = B(j,j) for j = 1,n. Prove this property for extra credit and try it numerically with B(j,j), j = 1,n defined as random matrix. The lesson to learn here is that depending on your spatial discretization scheme, there may be many simple ways to make the projection operation
un+1 = hI − BG(DBG)−1 Diun+1
accurate or even exact, and that choosing B to be the Taylor expansion of inv(A) might not always be the most efficient method.
