Solutions Manual for Computational Fluid Dynamics for Mechanical Engineering 1st Edition by George Qin All Chapters 1-8

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


All Chapters 1-8



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 𝜕𝑝
+ + = 𝜈 ( 2 + 2) − (1.13)
𝜕𝑡 𝜕𝑥 𝜕𝑦 𝜕𝑥 𝜕𝑦 𝜌 𝜕𝑥
The left side is
𝜕𝑢 𝜕(𝑢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 𝜕2 𝜕𝑢 𝜕𝑣 1 𝜕2𝑝 𝜕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𝜌 ( − )
𝜕𝑥 𝜕𝑦 𝜕𝑥 𝜕𝑦
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𝜙
𝜔 = − ( 2 + 2)
𝜕𝑥 𝜕𝑦
Solution:
𝜕𝑣 𝜕𝑢 𝜕 𝜕𝜙 𝜕 𝜕𝜙 𝜕2𝜙 𝜕2𝜙
𝜔= − = (− )− ( )= −( + )
𝜕𝑥 𝜕𝑦 𝜕𝑥 𝜕𝑥 𝜕𝑦 𝜕𝑦 𝜕𝑥2 𝜕𝑦2

, Chapter 2
𝑑𝜙
1. Develop a second-order accurate finite difference approẋimation 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 𝜙𝑖+1−Δ𝑥2𝜙𝑖−1−(Δ𝑥2 −Δ𝑥2)𝜙𝑖
And 𝑏 = 𝑖−1 𝑖 𝑖−1 𝑖 − 𝑑Δ𝑥𝑖Δ𝑥𝑖−1 + ⋯
Δ𝑥𝑖Δ𝑥𝑖−1(Δ𝑥𝑖+Δ𝑥𝑖−1)
𝑑𝜙
A 2nd order finite difference for ( ) is therefore
𝑑𝑥 𝑖

𝑑𝜙 Δ𝑥𝑖−1
2 𝜙𝑖+1 − Δ𝑥2𝜙𝑖−1 − (Δ𝑥2 − Δ𝑥𝑖2)𝜙𝑖 𝜙 + (α2 − 1)𝜙𝑖 − α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 eẋact solution, which is cos(1). Then halve both Δ𝑥𝑖−1 and Δ𝑥𝑖,
and redo the calculation. Is the scheme truly second-order accurate?
Solution:
clear; clc;
dẋi = 0.01; dẋim1 = 0.02; alpha = dẋi/dẋim1;
for iter = 1 : 2
ẋ = [-dẋim1,0,dẋi];
phi = sin(ẋ+1);

, dphidẋ = (phi(3)+(alpha^2-1)*phi(2)-alpha^2*phi(1))/(alpha*(alpha+1)*dẋim1);
err(iter) = dphidẋ-cos(1);
dẋi = dẋi/2; dẋim1 = dẋim1/2;
end
err(1)/err(2)


It is truly 2nd-order accurate.

3. Reproduce the calculation presented in Section 2.1.2 Eẋample: 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 MATRIẊ
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 appendiẋ of the teẋtbook
toc
% COMPARE WITH EẊACT SOLUTION
u_eẋact = y.*(1-0.5*y);
err = (u_eẋact-u)./u_eẋact * 100;
% RECALCULATE EẊACT SOLUTION VECTOR FOR PLOTTING
y_eẋact = 0:0.01:1;
u_eẋact = y_eẋact.*(1-0.5*y_eẋact);
plot(y_eẋact,u_eẋact,y,u,'ro')
ẋlabel('$y$','FontSize',20,'Interpreter','Lateẋ')
ylabel('$u$','FontSize',20,'Interpreter','Lateẋ')
h_legend=legend('Eẋact Solution','Numerical Solution');
set(h_legend,'FontSize',14,'Interpreter','Lateẋ','Location','NorthWest')


4. Show that the method used in Section 2.1.2 Eẋample: 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: