CEVE 455/555 CMOR 423/523
Homework 4
Due: 11:59 pm CST, Apr. 26th, 2024
Please give complete, well-written solutions to the following exercises and submit via Canvas. Collaboration is OK, but not encouraged. Indicate in your submission if you have collaborated with others and with whom.
1. 50 points (20 + 10 + 20). Consider 2D wave equation utt = c2 (uxx + uyy) with periodic boundary conditions on a square [ −L, L] × [ −L, L].
(a) Implement the leapfrog scheme (central difference in both temperal and spatial directions). Use 4th order RK to initiate the computation.
(b) u(x,y, t) = cos cos cos is a solution to the equation. Use this smooth ex- act solution to check the second order accuracy of the scheme. Let L = 2π and c = 1. The initial conditions are u(x,y,0) = cosxcosy and ut (x,y,0) = 0. Boundary conditions are pe- riodic. Fix ∆y = ∆x and ∆t = 0.9∆x/√2. Define the L∞ error at the final time T as maxi,j |Ui,j − u (xi
, yj , T)| where Ui,j is the numerical solution at the final time at the grid point (xi
, yj ). Run the code till T = 1 with Nx = Ny = 8, 16, 32, 64, 128, 256. Show both L∞ errors and the corresponding convergence order.
(c) Let L = 1.1 and T = 1. Use initial condition u(x,y,0) = e − (x2 +y2 )/a2 and ut (x,y,0) = 0 with
a = 0.02. Fix ∆t = 0.99∆x/√2. Run with N = 800 till T = 1. Plot three figures: i. the solution Ui,j as a 2D image with color bars;
ii. the solution along the line x = 0, i.e., a 1D slice of the 2D solution along x = 0;
iii. the diagonal of the solution, i.e., the 1D cut along the line x = y, label the horizontal axis with the distance to the origin
Do the same thing to another initial condition u(x,y,0) = 0 and ut (x,y,0) = e − (x2 +y2 )/a2 with a = 0.02.
2. 50 points (5 + 15 + 15 + 15). The vorticity stream-function formulation of the 2D incompressible NavierStokes is given by:
Here ψ is the stream function and ω is the vorticity, which is the curl of the velocity field u = ⟨u,v⟩ . Given ω, to find the velocity, first find ψ by solving ∆ψ = ω, then we have the velocity by computing u = −ψy , v = ψx. For simplicity, we only consider periodic cases to bypass the boundary condition.
(a) (2D Poisson Solver) We need to solve a Poisson equation every time step/stage. You can modify your code in problem 1 of HW2.
(b) (Linear Convection Diffusion) To solve the nonlinear convection diffusion, test the discretization on the linear one first. Consider the following equation with periodic boundary condition on [0, 2π] × [0, 2π] :
ut + ux + uy = d(uxx + uyy ), d > 0.
We can use centered difference for all spatial derivatives to achieve second order accuracy in space.
Let h = ∆x = ∆y denote mesh size of a uniform. mesh and ∆h denote the discrete Laplacian:
If we use explicit methods, the time step constaint for the Laplacian part will be ∆t ∼ d/h2 which is acceptable if d is very small (for instance d ∼ h ). However, using forward Euler (2) with d = 0 is never stable because its stability regions do not contain any imaginary axis. The stability regions of RK3 and RK4 contain part of the imaginary axis thus we could have a reasonable stable time step for small d. Use the Strong Stability Preserving (SSP) RK3 given in problem 3 of HW3 to solve (2) with the time step ∆t = min {0.3d/h2, 0.3h} (h = ∆x = ∆y). Test your code using the exact solution u(x, t) = exp(−2dt) sin(x+ y − 2t) with d = 0.01 on [0, 2π] × [0, 2π]. Run it till T = 0.2 with N = 16, 32, 64, 128, 256, 512. Show the table of L2
-errors and the corresponding convergence.
(c) (2D Incompressible Flow) Plugging the equations in the second line of (1) into the equation in the first line of (1), we get
Let Dx and Dy denote the central difference for the first order partial derivatives, we get
Use the RK3 above for the time derivative. At each time step tn , first compute the maximum of ve- locity by Un = max{|un | , |vn |}, then time step can betaken as ∆t = min {0.3Re2, 0.3Un/1 = ∆x = ∆y). For each time stage, you need to compute/update the velocity by solving the Poisson equation ψ. Test the accuracy with the exact solution: ω(x, y, t) = −2 exp(−2t/Re) sin x sin y on domain [0, 2π] × [0, 2π] with Re = 100. Notice that the 2D Poisson equation needs to be solved for each time stage in one RK step. Run your code till T = 0.2 with N = 16, 32, 64, 128, 256. List the L2-errors and the corresponding convergence order as a table.
(d) (Implicit diffusion treatment) For small Re, the time step ∆t = 0.3Reh2 is too small to use. Instead, we can consider using the IMEX method, i.e., consider the following scheme (first order accurate in time):
Take Re = 70. Implement the scheme (3) with the time step size ∆t = ∆x, and the initial data
on domain [0, 2π] × [0, 2π], where we take ρ = π/15 and δ = 0.05. Plot your vorticity at T = 7 with N = 512.