Computational Fluid Dynamics: Homework #2
Due on March 16, 2020 at 5:00 pm
1
Computational Fluid Dynamics : Homework #2
Problem 1
Problem 1
[5 Points] Solve the initial value problem for the model hyperbolic equation (the one-dimensional linear
advection equation), described by the partial differential equation
∂u
∂u
+c
= 0,
∂t
∂x
(1)
using the five explicit methods provided below. Bonus points will be given for correct solutions obtained
from the MacCormack and Warming-Beam methods. For each numerical method, use N = 41 points to
span the interval 0 ≤ x ≤ Lx with a Lx = 2. Assume the mesh is uniform in the x-direction, where the
location at each point is given by
xi = (i − 1)∆x
with ∆x = Lx /(N − 1)
For each method use the stability condition,
c∆t
= 0.9
∆x
to compute the time step, ∆t. Use the value c = 1. Run each method for 10 time steps and compare the
results to the analytical solution, u(x, t) = u0 (x − ct) For reference, the derivation for the analytical solution
to the linear convection problem using the method of characteristics is provided below. Use the following
initial condition,
(
1.0 x ≤ 0.5
u0 (x, 0) =
0.5 x > 0.5
Since the wave is moving to the right (c > 0), the boundary condition at x = 0 is fixed by the initial
condition, u(0, t) = u0 (0, t). At the right boundary x = 2, always use the backward difference formula,
(D− /∆x)·, for the spatial derivative to avoid needing to specify a ”ghost” node at that boundary.
Explicit Numerical Methods
• Explicit Backward [1 Point]
D− · n
u
∆x i
∆x
Stability Criteria : c > 0, ∆t ≤
|c|
Order of Accuracy : O (∆t, ∆x)
Method : un+1
= uni − c∆t
i
(2)
• Explicit Forward [1 Point]
D+ · n
u
∆x i
∆x
Stability Criteria : c > 0, ∆t ≤
|c|
Order of Accuracy : O (∆t, ∆x)
Method :
un+1
= uni − c∆t
i
(3)
• Explicit Central [1 Point]
Method : un+1
= uni − c∆t
i
Stability Criteria :
Always unstable
Order of Accuracy : O ∆t, ∆x2
Problem 1 continued on next page. . .
2
D0 · n
u
∆x i
(4)
Computational Fluid Dynamics : Homework #2
Problem 1 (continued)
• Lax [1 Point]
uni−1 + uni+1
D0 · n
− c∆t
u
2
∆x i
∆x
Stability Criteria : ∆t ≤
|c|
Order of Accuracy : O ∆t, ∆x2
Method : un+1
=
i
(5)
• Lax-Wendroff [1 Point]
Method :
Stability Criteria :
un+1
= uni − c∆t
i
∆t ≤
D0 · n 1 2 2 D+ · D− · n
u + c ∆t
u
∆x i
2
∆x ∆x i
(6)
∆x
|c|
Order of Accuracy : O ∆t2 , ∆x2
• MacCormack [Bonus: 1 Point]
D+ · n
Method : un+1
u
= uni − c∆t
i
∆x i
1 n
D− · n+1
n+1
n+1
ui
=
u + ui − c∆t
u
2 i
∆x i
∆x
Stability Criteria : ∆t ≤
|c|
2
Order of Accuracy : O ∆t , ∆x2
(7)
• Warming-Beam [Bonus: 1 Point]
c∆t D− · n
ui
2 ∆x
D− · n+1/2 ∆x D− · n
n
un+1
=
u
−
c∆t
u
u
+
i
i
i
∆x
2 ∆x i
2∆x
Stability Criteria : ∆t ≤
|c|
Order of Accuracy : O ∆t2 , ∆x2
Method :
n+1/2
ui
= uni −
(8)
Difference Formulas
Assuming a uniform mesh size of ∆x, the difference operators, D, used in the difference equations for the
explicit numerical methods given in the previous section are defined as
• Forward difference, first-order accurate
φi+1 − φi
D+ ·
φi =
∆x
∆x
(9)
• Backward difference, first-order accurate
D− ·
φi − φi−1
φi =
∆x
∆x
(10)
• Central difference, second-order accurate
D0 ·
φi+1 − φi−1
φi =
∆x
2∆x
Problem 1 continued on next page. . .
3
(11)
Computational Fluid Dynamics : Homework #2
Problem 1 (continued)
Analytical Solution to a Hyperbolic PDE
The solution to Eq. 1 can be determined from the initial condition u(x, 0) = u0 (x) and the boundary
condition, u(0, t) = u0 (0, t) = ul . If c ≥ 0, then we only need the boundary condition at x = 0. Also, note
that the boundary condition is informed by the initial condition; hence we are providing solutions to the
initial value problem. The method of characteristics can be used to find an analytical solution to the above
equation. Along a characteristic line, we have
dx
= c.
(12)
dt
As a result, along a particular curve x = x(t), we have
i ∂u dx ∂u
dh
∂u
∂u
u(x(t), t) =
+
=
+c .
(13)
dt
∂t
dt ∂x
∂t
∂x
The value u is therefore constant along the characteristics (the material derivative du/dt is zero, thus u is
constant along lines x − t lines). We also know that these characteristic lines are straight, since the curvature
(second-derivative) of the characteristics is zero,
d dx
dc
d2 x
=
=0
(14)
=
2
dt
dt dt
dt
Knowing that the value of u is constant along a characteristic, we can define a series of parametric lines that
only depend on the initial value of u at a point emanating from (η, 0). Direct integration of
dx
= c.
dt
(15)
x = ct + ξ
(16)
ξ = x − ct
(17)
results in the characteristic line
or
where the solution of u(x, t) of the initial value problem is given by the
u(x, t) = u(ξ, 0) = u0 (ξ) = u0 (x − ct)
(18)
Problem 2
[2 Points] Solve the initial value problem for the model hyperbolic equation,
∂u
∂u
+c
=0
∂t
∂x
using the implicit central difference method, given as
Method : un+1
= uni − c∆t
i
Stability Criteria :
Always stable
Order of Accuracy : O ∆t, ∆x2
(19)
D0 · n+1
u
∆x i
(20)
(21)
(22)
Use the same mesh as used in Problem 1, but for c∆t
∆x use the values 0.9, 1.8, 2.0. As in Problem 1, compare
the computed solutions to the analytical solution after 10 time steps.
Note: The value c∆t
∆x is often referred to as the CFL number, in reference to Richard Courant, Kurt
Friedrichs, and Hans Lewy, who described the stability of hyperbolic PDEs in a 1928 paper. While the
criteria is a direct result of von Neumann stability analysis (deduced many years later), Courant, Friedrichs,
and Lewy derived the result more intuitively. They reasoned that during a single discrete time step, a wave
(characteristic) moving across a mesh must not move further than the adjacent grid point during that time
step.
4
Computational Fluid Dynamics : Homework #2
Problem 3
Problem 3
[3 Points] Solve the boundary value problem for the model parabolic equation, described by the partial
differential equation
∂2u
∂u
(23)
=ν 2
∂t
∂x
using both the full explicit (α = 0) and full implicit (α = 1), central difference methods described as
Dxx ·
(1 − α)uni + αun+1
i
∆x
Stability Criteria : α = 1 : Always stable
Method : un+1
= uni − ν∆t
i
α = 0 : ∆t ≤ ∆x2 /(2k)
Order of Accuracy : O ∆t, ∆x2
(24)
(25)
where the second-order spatial difference equation is
• Central difference, second-order accurate
– Uniform mesh
Dxx ·
φi+1 − 2φi + φi−1
φi =
∆x
∆x2
(26)
φi+1 −φi
φi −φi−1
xi+1 −xi − xi −xi−1
1
2 (xi+1 − xi−1 )
(27)
– Non-uniform mesh
Dxx ·
φi =
∆x
Use N = 41 points to span the interval 0 ≤ x ≤ Lx with Lx = 1. Assume the mesh is uniform in the
x-direction, where the location at each point is given by
xi = (i − 1)∆x
with ∆x = Lx /(N − 1)
(28)
The parabolic PDE is a boundary value problem and requires us to specify a boundary condition at each
boundary in the domain. Assume ν = 1 and use a Dirichlet boundary condition, i.e., fixed value of u at each
end of the domain, x = 0 and x = Lx = 1. The initial condition and boundary conditions are given as
Left boundary condition x = 0 : u(0, t) = 0
Right boundary condition x = Lx : u(Lx , t) = 1
Initial condition t = 0 : u(x, 0) = 0
(29)
(30)
(31)
For the full explicit (α = 0) method, use the stability constraint to compute the time step, ∆t, and run
for 100 time steps. For the full implicit (α = 1) method, use two values of ∆t. For one simulation use
∆t = 1/8 for 10 time steps, and for a second simulation use ∆t = 109 for 1 time step. Plot the results. The
steady-state solution of u(x) should be a linear function bounded by 0 at x = 0 and 1 at x = Lx .
Problem 4
[5 Points] Use the implicit central difference method as described Eqs. 24 to solve heat equation,
∂T
k ∂2T
∂2T
=
=ν 2,
2
∂t
ρcp ∂x
∂x
(32)
where ν = k/(ρcp ) is the thermal diffusivity. In the previous problem, we essentially solved a homogenous,
non-dimensional form of the heat equation. In this problem, let us find the solution to a more practical
Problem 4 continued on next page. . .
5
Computational Fluid Dynamics : Homework #2
Problem 4 (continued)
problem. Use the following values, ρ = 8000 kg/m3 , cp = 500 J/kg · K, and k = 10 W/m · K. At the left
boundary, assume we have been given a known constant heat flux, q0 , and at theright boundary, we will
assume the surface is adiabatic, meaning the heat flux across that boundary is zero. The initial and boundary
conditions as
Left BC - Specified Heat Flux, x = 0 : qw (0) = q0 = −7.5 × 105 W/m2
Right BC - Adiabatic, x = Lx : qw (Lx ) = 0.0
Initial condition, t = 0 : T (x, 0) = T0 = 300K
(33)
(34)
(35)
Note that the heat flux is defined as negative. Following sign convention, this is the heat flux into the
material. We can find an analytical solution to the heat condition, which is given by the equation,
1
νt
x
1
T (x, t) − T0
= 2 + −
+
Lx
Lx
3 Lx
2
|q0 | k
x
Lx
2
−2
∞
X
nπx
1
2 2 νt
exp
−n
π
cos
n2
L2x
Lx
n=1
(36)
Compare your numerical solutions to the analytical solution at t = 40s. To evaluate the Fourier series in the
analytical solution, you can truncate the series at n = 1000.
Note: For long times the scenario described by the boundary conditions above is unphysical, since as
time continues (t → ∞) more and more energy is added to the system, q0 is constant. In practice, q0 is
not constant. At high temperatures, black-body radiation at this surface must be taken into account. This
would imply that the boundary condition is a function of temperature. A mixed boundary condition like
this is called a Robin boundary condition.
Bonus: 3 points
Find a solution to the heat conduction equation as described above, but instead of using a uniform mesh,
use a non-uniform mesh that is hyperbolically spaced. Compare your results to those obtained from using
a uniform mesh. The mesh sizing can be computed using the following procedure. Given the length of the
domain as Lx , the number of points in the mesh as N , the location of the first node index as x1 , and the
minimum spacing as ∆xmin = x2 −x1 , then the x-location of each node in the mesh is given by the hyperbolic
formula
i−1
xi = x1 + Lx
e N −1 − 1
eκ − 1
(37)
where the value of ∆xmin can be used to determine stretching factor κ. Using Newton’s method (a rootfinding algorithm), we can iteratively solve for the value of κ using the equation,
1
x2 − x1 = ∆xmin = Lx
e N −1 − 1
eκ − 1
(38)
Once κ is determined, the mesh can be defined using Eq. 37. Use the same implicit central difference
equation, however, the difference formula for the second derivative should account for the non-uniform mesh
size.
6
...

Purchase answer to see full
attachment