Solution Manual for Computational Fluid Dynamics for Mechanical Engineering 1st Edition by George Qin

Solutions Manual for Computational Fluid
Dynamics for Mechanical Engineering, 1e by
George Qin (All Chapters)


Chapter 1
1. Show that Equation (1.14) can also be written as
𝜕𝑢 𝜕𝑢 𝜕𝑢 𝜕2𝑢 𝜕2𝑢 1 𝜕𝑝
+𝑢 +𝑣 = 𝜈 ( 2 + 2) −
𝜕𝑡 𝜕𝑥 𝜕𝑦 𝜕𝑥 𝜕𝑦 𝜌 𝜕𝑥
Solution
Equation (1.14) is
𝜕𝑢
𝜕(𝑢2) 𝜕(𝑣𝑢) 𝜕2𝑢 𝜕2𝑢 1 𝜕𝑝
+ +
𝜕𝑡 𝜕𝑥 = 𝜈( + )− (1.13)
The left side is 𝜕𝑦 𝜕𝑥2 𝜕𝑦2 𝜌 𝜕𝑥
𝜕𝑢
𝜕(𝑢2) 𝜕(𝑣𝑢) 𝜕𝑢 𝜕𝑢 𝜕𝑢 𝜕𝑣
𝜕𝑡 + + = + 2𝑢 +𝑣 +𝑢
𝜕𝑥 𝜕𝑦 𝜕𝑡 𝜕𝑥 𝜕𝑦 𝜕𝑦
𝜕𝑢 𝜕𝑢 𝜕𝑢 𝜕𝑢 𝜕𝑣 𝜕𝑢 𝜕𝑢 𝜕𝑢
= +𝑢 +𝑣 +𝑢( + ) = +𝑢 +𝑣
𝜕𝑡 𝜕𝑥 𝜕𝑦 𝜕𝑥 𝜕𝑦 𝜕𝑡 𝜕𝑥 𝜕𝑦
since
𝜕𝑢 𝜕𝑣
+ =0
𝜕𝑥 𝜕𝑦
due to the continuity equation.
2. Derive Equation (1.17).
Solution:
From Equation (1.14)
𝜕𝑢 𝜕(𝑢2) 𝜕(𝑣𝑢) 𝜕2𝑢 𝜕2𝑢
1 𝜕𝑝
+ + = 𝜈 ( + ) −
𝜕𝑡 𝜕𝑥 𝜕𝑦 𝜕𝑥2 𝜕𝑦2 𝜌 𝜕𝑥
Define 𝑥𝑖 𝑡𝑈 𝑝
𝑢̃ 𝑢 𝑣 = , 𝑡̃ = , 𝑝̃ =
𝑥̃ = , 𝑣̃ = ,
𝑈 𝑈 𝑖
𝐿 𝐿 𝜌𝑈2
Equation (1.14) becomes
𝑈𝜕𝑢̃ 𝑈 𝜕(𝑢̃ ) 𝑈2𝜕(𝑣̃𝑢 ̃)
2 2
𝜈𝑈 𝜕2𝑢̃ 𝜕2𝑢̃ 𝜌𝑈2 𝜕𝑝̃
+ + = ( + ) −
𝐿 𝐿𝜕𝑥 𝐿𝜕𝑦 𝐿2 𝜕𝑥2 𝜕𝑦2 𝜌𝐿 𝜕𝑥
𝜕𝑡̃
𝑈
Dividing both sides by 𝑈2/𝐿, Equation (1.17) follows.

3. Derive a pressure Poisson equation from Equations (1.13) through (1.15):

, 𝜕2𝑝 𝜕2𝑝 𝜕𝑢 𝜕𝑣 𝜕𝑣 𝜕𝑢
+ = 2𝜌 ( − )
𝜕𝑥2 𝜕𝑦2 𝜕𝑥 𝜕𝑦 𝜕𝑥 𝜕𝑦
Solution:
𝜕𝑢 𝜕𝑣
𝜕𝑥 + =0 (1.13)
𝜕𝑦 2
𝜕𝑢 𝜕(𝑢 )2
𝜕(𝑣𝑢) 𝜕𝑢 2
𝜕𝑢 1 𝜕𝑝
+ +
𝜕𝑡 𝜕𝑥 = 𝜈( 2
+ 2
)− (1.14)
𝜕𝑣 𝜕(𝑢𝑣) 𝜕(𝜕𝑦
𝑣2) 𝜕𝑥
𝜕2𝑣 𝜕𝑦
𝜕2𝑣 𝜌1𝜕𝑥
𝜕𝑝
+ +
𝜕𝑡 𝜕𝑥 = 𝜈 ( 2 + 2) − (1.15)
𝜕𝑦 𝜕𝑥 𝜕𝑦 𝜌 𝜕𝑦
Taking 𝑥-derivative of each term of Equation (1.14) and 𝑦-derivative of each term of Equation (1.15),
then adding them up, we have

𝜕 𝜕𝑢 𝜕𝑣 𝜕2(𝑢2) 𝜕2(𝑣𝑢) + 𝜕 (𝑣 )
2 2

𝜕𝑡 ( + ) +
𝜕𝑥 𝜕𝑦 2
+ 2 𝜕𝑥𝜕𝑦 𝜕𝑦2
𝜕𝑥
𝜕2 𝜕𝑢 𝜕𝑣 1 𝜕2𝑝 𝜕2𝑝
𝜕2
= 𝜈( + )( + )− ( 𝜕𝑥 2 + 𝜕𝑦2)
Due to continuity, we have 𝜕𝑥2 𝜕𝑦2 𝜕𝑥 𝜕𝑦 𝜌
𝜕2𝑝
𝜕2𝑝 𝜕2(𝑢2) 𝜕2(𝑣𝑢) 𝜕2(𝑣2)
+ ]
2
+ 2
= −𝜌 [ 2
+ 2 𝜕𝑥𝜕𝑦 𝜕𝑦2
𝜕𝑥 𝜕𝑦 𝜕𝑥
= −2𝜌(𝑢𝑥𝑢𝑥 + 𝑢𝑢𝑥𝑥 + 𝑢𝑥𝑣𝑦 + 𝑢𝑣𝑥𝑦 + 𝑢𝑥𝑦𝑣 + 𝑢𝑦𝑣𝑥 + 𝑣𝑦𝑣𝑦 + 𝑣𝑣𝑦𝑦)
𝜕 𝜕 𝜕𝑢 𝜕𝑣
= −2𝜌 [(𝑢𝑥 + 𝑢 +𝑣 ) ( + ) + 𝑢𝑦𝑣𝑥 + 𝑣𝑦𝑣𝑦]
𝜕𝑥 𝜕𝑦 𝜕𝑥 𝜕𝑦
𝜕𝑣 𝜕𝑢 𝜕𝑢 𝜕𝑣
= −2𝜌(𝑢𝑦𝑣𝑥 + 𝑣𝑦𝑣𝑦) = −2𝜌(𝑢𝑦𝑣𝑥 − 𝑢𝑥𝑣𝑦) = 2𝜌 ( − )
𝜕𝑥 𝜕𝑦 𝜕𝑥 𝜕𝑦
4. For a 2-D incompressible flow we can define the stream function 𝜙 by requiring
𝑢 = 𝜕𝜙 𝜕𝜙
; 𝑣=−
𝜕𝑦 𝜕𝑥
We also can define a flow variable called vorticity
𝜕𝑣 𝜕𝑢
𝜔= −
𝜕𝑥 𝜕𝑦
Show that
𝜕2𝜙 𝜕2𝜙

𝜔 = −( + )
Solution: 𝜕𝑥2 𝜕𝑦2
𝜕𝑣
𝜕𝑢 𝜕 𝜕𝜙 𝜕 𝜕𝜙 𝜕2𝜙 𝜕2𝜙
𝜔= − = (− )− (
𝜕𝑥 𝜕𝑦 𝜕𝑥 𝜕𝑥 𝜕𝑦 𝜕𝑦 ) = − ( 2 + 2 )
𝜕𝑥 𝜕𝑦

, Chapter 2 𝑑𝜙
1. Develop a second-order accurate finite difference approximation for ( ) on a non-uniform
𝑑𝑥 𝑖
mesh using information (𝜙 and 𝑥 values) from mesh points 𝑥𝑖−1, 𝑥𝑖 and 𝑥𝑖+1. Suppose 𝛿𝑥𝑖 =
𝑥𝑖+1 − 𝑥𝑖 = 𝛼𝛿𝑥𝑖−1 = 𝛼(𝑥𝑖 − 𝑥𝑖−1).

Solution:

Assume close to the 𝑖𝑡ℎ point, 𝜙(𝑥) = 𝜙𝑖 + 𝑏(𝑥 − 𝑥𝑖) + 𝑐(𝑥 − 𝑥𝑖)2 + 𝑑(𝑥 − 𝑥𝑖)3 …
Then 𝑑𝜙 = 𝑏 + 2𝑐(𝑥 − 𝑥 ) + ⋯ and 𝑑𝜙) = 𝑏.
𝑖 (
𝑑𝑥 𝑑𝑥 𝑖
Now 𝜙𝑖+1 = 𝜙(𝑥𝑖+1) = 𝜙𝑖 + 𝑏(𝑥𝑖+1 − 𝑥𝑖) + 𝑐(𝑥𝑖+1 − 𝑥𝑖)2 + ⋯ = 𝜙𝑖 + 𝑏Δ𝑥𝑖 + 𝑐Δ𝑥2 + 𝑑Δ𝑥3 …
𝑖 𝑖
And 𝜙𝑖−1 = 𝜙(𝑥𝑖−1) = 𝜙𝑖 + 𝑏(𝑥𝑖−1 − 𝑥𝑖) + 𝑐(𝑥𝑖−1 − 𝑥𝑖)2 + ⋯ = 𝜙𝑖 − 𝑏Δ𝑥𝑖−1 + 𝑐Δ𝑥 2
− 𝑑Δ𝑥3 …
𝑖−1 𝑖−1
So Δ𝑥2 𝜙𝑖+1 − Δ𝑥2𝜙𝑖−1 = (Δ𝑥2 − Δ𝑥2)𝜙𝑖 + 𝑏Δ𝑥𝑖Δ𝑥𝑖−1(Δ𝑥𝑖 + Δ𝑥𝑖−1) + 𝑑Δ𝑥2Δ𝑥2 (Δ𝑥𝑖 +
𝑖−1 𝑖 𝑖−1 𝑖 𝑖 𝑖−1
Δ𝑥𝑖−1) + ⋯
2 2 2 2
Δ𝑥 𝜙𝑖+1−Δ𝑥 𝜙𝑖−1−(Δ𝑥 −Δ𝑥 )𝜙𝑖
And 𝑏 = 𝑖−1 𝑖 𝑖−1 𝑖 − 𝑑Δ𝑥𝑖Δ𝑥𝑖−1 + ⋯
Δ𝑥𝑖Δ𝑥𝑖−1(Δ𝑥𝑖+Δ𝑥𝑖−1)
𝑑𝜙
A 2nd order finite difference for ( ) is therefore
𝑑𝑥 𝑖

𝑑𝜙 2
( ) Δ𝑥𝑖−1 𝜙𝑖+1 − Δ𝑥𝑖2𝜙𝑖−1 − (Δ𝑥2 𝑖−1 − Δ𝑥𝑖 2)𝜙𝑖 ( 2 ) 2
𝑑𝑥 𝑖 = 𝑏 ≈ = 𝜙𝑖+1 + α − 1 𝜙𝑖 − α 𝜙𝑖−1
Δ𝑥𝑖Δ𝑥𝑖−1(Δ𝑥𝑖 + Δ𝑥𝑖−1) α(α + 1)Δ𝑥𝑖−1

2. Use the scheme you developed for problem 1 to evaluate the derivative of 𝜙(𝑥) =
sin(𝑥 − 𝑥𝑖 + 1) at point 𝑖. Suppose Δ𝑥𝑖−1 = 0.02 and Δ𝑥𝑖 = 0.01. Compare your
numerical result with the exact solution, which is cos(1). Then halve both Δ𝑥𝑖−1 and Δ𝑥𝑖,
and redo the calculation. Is the scheme truly second-order accurate?
Solution:
clear; clc;
dxi = 0.01; dxim1 = 0.02; alpha = dxi/dxim1;
for iter = 1 : 2
x = [-dxim1,0,dxi];
phi = sin(x+1);

, dphidx = (phi(3)+(alpha^2-1)*phi(2)-alpha^2*phi(1))/(alpha*(alpha+1)*dxim1);
err(iter) = dphidx-cos(1);
dxi = dxi/2; dxim1 = dxim1/2;
end
err(1)/err(2)


It is truly 2nd-order accurate.

3. Reproduce the calculation presented in Section 2.1.2 Example: Laminar Channel Flow.
Solution:
% %
% FULLY-DEVELOPED CHANNEL FLOW %
% By George Qin for "A Course of Computational Fluid Dynamics" %
% %
tic
clear; clc;
% PARAMETERS
H = 1; N = 5;
% FACE LOCATIONS
yf = linspace(0,H,N+1);
% NODE LOCATIONS
y = 0.5*(yf(1:end-1)+yf(2:end));
% DELTA Y
dy = yf(2)-yf(1);
% THE THREE DIAGONAL VECTORS IN THE COEFFICIENT MATRIX
as = -(1/dy^2)*ones(1,N);
ap = (2/dy^2)*ones(1,N);
an = -(1/dy^2)*ones(1,N);
b = ones(1,N);
% SPECIAL VALUES AT THE BOUNDARIES (BOUNDARY CONDITIONS)
%ap(1) = 3/dy^2; as(1) = 0;
%ap(1) = 4/dy^2; an(1) = - (4/3)/dy^2; as(1) = 0;
ap(1) = 3/dy^2; as(1) = 0; b(1) = 3/4;
ap(end) = 1/dy^2; an(end) = 0;
% SOLVE THE SYSTEM OF LINEAR EQUATIONS WITH TDMA ALGORITHM
u = TDMA(as,ap,an,b); % this is in the appendix of the textbook
toc
% COMPARE WITH EXACT SOLUTION
u_exact = y.*(1-0.5*y);
err = (u_exact-u)./u_exact * 100;
% RECALCULATE EXACT SOLUTION VECTOR FOR PLOTTING
y_exact = 0:0.01:1;
u_exact = y_exact.*(1-0.5*y_exact);
plot(y_exact,u_exact,y,u,'ro')
xlabel('$y$','FontSize',20,'Interpreter','Latex')
ylabel('$u$','FontSize',20,'Interpreter','Latex')
h_legend=legend('Exact Solution','Numerical Solution');
set(h_legend,'FontSize',14,'Interpreter','Latex','Location','NorthWest')


4. Show that the method used in Section 2.1.2 Example: Laminar Channel Flow is first-
order accurate by using the global error estimate technique.
Solution:

These finite difference equations and their truncation errors are reproduced here: