HOMEWORK # 3: FINITE VOLUMES FOR TWO DIMENSIONAL PROBLEMS, V1.2
DUE: WEDNESDAY APRIL 18TH 2018
1. Two-Dimensional, Time Dependent Advection Equation
In this homework problem, you will solve the two-dimensional convection equation using Finite Volume Methods. The governing partial differential equation we will solve is:
∂ω + ∇ · V⃗ (x, y) · ω = 0 ∂t
Where, ω is the vorticity in the domain. We are interested in how this ω moves through the domain due to the velocity V⃗ (x, y) in the domain. This velocity, will in general, be different at different locations of the domain. The equation will be solved on a unit-square domain with cell-centered ω-values while the fluxes (ω · V⃗ ) will be defined at the cell faces/interfaces (see figure 1).
Figure 1. The finite volume method domain and boundary conditions for this homework.
1
2 DUE: WEDNESDAY APRIL 18TH 2018
The FVM method will be used to advect/convect the distribution of ω(x,y) in time. A structured grid in the x− and y− coordinate systems will be used. The ω value will be assumed given at time t = 0 and have a single value in any cell of the domain (each cell would have a different value, however, the value in a given cell is single-valued or constant). Periodic boundary conditions will be applied on all boundaries to allow long time evolution of the convection/advection equation. A simple forward Euler explicit time integration method will be implemented.
Figure 2. The spatial discretization assumes a constant valued ω within the cell (with the value stored at the ”cell center” and the flux through each of the N,S,E and W boundaries is uniform and evaluated at the face center.
The finite volume method will assume that each 2D finite volume (finite area in the 2D setting) will contain a constant amount of vorticity ω. The flux through the North, South, East and West faces of the cell will be assumed constant across the length of the face. The boundary conditions for the convection equation are periodic on all sides. This means that any flow or flux exiting the boundary will re-enter on the opposite boundary. The domain in Figure 1 shows ghost cells around the boundary that could help with the implementation of the periodic boundary condition. The first order finite volume spatial discretization for a cell inside of the domain is shown in Figure 2.
The governing equation can be discretized using Finite Volume Methods (with k repre- senting the current time-step) as:
(1) ∂ω + ∇ · (V⃗ (x, y)ω) = 0 ∂t
Which, when discretized in 2D on a regular grid becomes: (2)
ωk+1=ωk −∆tf(ωk,Vk) +∆tf(ωk,Vk) −∆tf(ωk,Vk) +∆tf(ωk,Vk) i,j i,j ∆x x E ∆x x W ∆y y N ∆y y S
(3)
(4)
(5)
(6)
f(ωk,Vk) = x i,j E
f(ωk,Vk) = x i,j W
f(ωk,Vk) = y i,j N
f(ωk,Vk) = y i,j S
1V ωk +ωk − 2xi+1,ji,j
1|V |ωk −ωk
HOMEWORK # 3: FINITE VOLUMES FOR TWO DIMENSIONAL PROBLEMS, V1.2 3
Where, we use an up-winding scheme to calculate the fluxes between cells:
2 x i+1,j 1|V |ωk −ωk
i,j
1V ωk +ωk 2xi,ji−1,j
− 1V ωk +ωk −
2 x i,j i−1,j
2yi,j+1i,j
1V ωk +ωk −
1|V |ωk
2 y i,j+1
−ωk i,j
2yi,ji,j−1
1|V |ωk −ωk 2 y i,j i,j−1
This equation, once correctly implemented in Matlab, will be able to describe the convection of ω based on an arbitrary velocity field. In the questions below, different velocities are used to determine how the ω evolves over time.
2. Introducing the bonus problem – A Vorticity-Velocity-Streamfunciton CFD solver
Figure 3. The grid used for solving the velocity from the vorticity via the streamfunction. The vorticity ω is projected to the corner vertices (red squares) where the finite difference method can be used to determine the streamfunction Ψ. The derivatives in x− and y− directions can be taken to recover velocity at the cell face centers. This is a relatively approximate determination of the velocity.
Later in the problem set, you may explore a bonus question in which ω is responsible for creating the velocity field V⃗ . This effectively becomes a velocity-vorticity method in which
4 DUE: WEDNESDAY APRIL 18TH 2018
distributed vorticity drives the velocity that advects that same vorticity distribution. The vorticity is related to the velocity by introducing a 2D stream-function:
With the boundary conditions:
∇2Ψ(x, y) = −ω(x, y) Ψ=0
On the outside of the domain. This effectively sets up a no-normal flow through the outer boundary.
This equation can be solved on a grid using the methods developed at the beginning of the class for elliptic problems (e.g., the membrane problem).
The velocity at the face centers can be calculated from the streamfunction solution as:
∂Ψ ∂y ∂Ψ ∂x
3. Homework Questions
Please attempt all questions. The bonus is optional. Please hand in code on turn-it-in. Please provide paper copies of all pertinent graphs and the answers to all questions unless otherwise explicitly indicated.
(1) Fill in the steps of the derivation between the continuous PDE (1) and the finite volume discretization of that PDE (equation 2).
(2) Show/describe that your final equation is a conservation equation.
(3) Show/describe how the up-winding flux equations select the correct ω at the cell
faces to apply to the convection.
(4) Describe how the East and West fluxes (or North and South) fluxes are related
from adjacent cells.
(5) Describe how you will implement your flux calculations (equations 3-6) and the
ω update equation (equation 2). Based on your previous answer, indicate how you might (you don’t have to) develop an efficient solver/implementation that only calculates the fluxes once rather than twice.
(6) Write a pseudo code for the FVM computer code.
(7) Implement the Matlab code with the following t = 0 initial condition for ω(x, y):
(x−0.5)2+(y−0.5)2 ω(x,y) = e− l2
l = 0.1
(8) The study of the method: For the following constant convection speeds and timesteps across the domain, plot how the maximum value of the wave being convected changes w.r.t time. In addititon, on a separate plot, show how the total ω (i.e.,
= Vx(x,y) = Vy(x,y)
HOMEWORK # 3: FINITE VOLUMES FOR TWO DIMENSIONAL PROBLEMS, V1.2 5
ωdA) varies with respect to time on the entire domain when dt = dx and domain
dt = dx/2. Do this for dx = dy = 0.1,0.05,0.025. Comment on/explain your findings:
(a) Configuration 1: Vx = 1, Vy = 0, dt = dx/2 (b) Configuration2: Vx =1,Vy =0,dt=dx
(c) Configuration 3: Vx = 1, Vy = 1, dt = dx/2 (d) Configuration4: Vx =1,Vy =1,dt=dx
(9) Using the method for the good of humankind: In this problem, you will take your code developed in the prior sections and make some small modifications to solve an actual problem – the evolution of a concentration of ω due to a Stommel Gyre 1. This could be thought of as assessing where a pollutant (ω) flows over time due to the background velocity of the Gyre.
To do this, you will need to modify the code so that the following changes occur: • The velocity is represented by the expressions of a Strommel Gyre – you may
wish to quiver plot this velocity to make sure it looks reasonable:
u(x, y) = L v(x,y) = −πC1
−1 −C2x L +C2e
πy cos( L )
πy sin(L)
π2C1 (L − x)
−C2x L − e
1 20
C2 = 10 L=1
C1 =
Once again, the distribution of the pollutant will be a 2D-Gaussian defined by:
(x−0.5)2+(y−0.5)2 ω(x,y) = e− l2
l = 0.1
• You will need to setup boundary conditions that impose no flux through the outer boundary of the domain.
(10) Describe (but do not implement) the changes you would need to the equations/to your code to:
• Be able to use an unstructured grid of triangular cells with a finite volume solver?
• Improve the accuracy of your solver?
• Implement timesteps dt > dx
(11) Bonus 1: Vectorize your code for improved performance. Determine how much
faster the vectorized version is than the loop version.
(12) Bonus 2: Setup a general vorticity-velocity solver by convecting the vorticity based
on the velocity determined from a finite difference solver of ∇2Ψ = −ω. For this
1http://www.meteo.psu.edu/holocene/public html/Mann/courses/METEO470SPR17/Readings/StommelGyre Haldvogel Trenberth.pdf
6
DUE: WEDNESDAY APRIL 18TH 2018
problem, use the code from question 9 where you have no-flux boundary conditions in the domain. This means you may also set Ψ = 0 for the outside nodes of the finite difference solver.
For this problem, implement the following ’dancing vortex pair’ problem:
(x−0.45)2+(y−0.45)2 (x−0.65)2+(y−0.65)2 ω(x,y) = e− l2 +e− l2
l = 0.1