COMPUTATIONAL HEAT TRANSFER (ME 673) FINAL PROJECT LID DRIVEN CAVITY FLOW
Name: Eldwin Djajadiwinata Student ID: 434107763 Lecturer: Dr. Shereef Sadek
1
Numerical Heat Transfer ME673 Final-Project Second Semester 1435 th
Due date Wednesday June 4
Consider the two-dimensional steady state flow of air in a closed lid-driven square cavity as shown in Figure in Figure 1, Lid-driven cavity flow Figure 1,
⁄ Figure 1, Lid-driven cavity flow
defined over the domain , subject to the shown boundary boun dary conditions. Requirements: 1. Using vorticity-stream function approach write down the non-dimensional governing equations using as reference values 2. Discretize the governing PDEs using finite volume discretization 3. Solve the flow problem described above using ADI scheme for the following parameters: 4. Compare results from part 3 with results in the literature 5. Plot stream-function and vorticity contours, estimate the position and magnitude of maximum stream-function and vorticity values 6. Compute the temperature field 7. Calculate and plot number along the horizontal walls defined as
Where (*) denotes non-dimensional variables. Note: All submitted work should be entirely yours. Include source code in your report.
2
Question 1 Transforming the governing equations into non -dimensionalized forms. The governing equations are as follows.
⃗ ⃗
(1)
(2)
(3)
Let us first define the following non-dimensional variables.
(4)
a) Transforming the continuity/stream function equation
(5)
Substitute the non-dimensional variables in Eq. (4) Eq. (4) into Eq. (5) 3
Finally, we will get the non-dimensional continuity/function equ ation
(6)
b) Transforming the vorticity equation
⃗ If we open the Del operator and multiply the dot matrix we get:
(7)
Again, substitute the non-dimensional variable into the equation above.
4
Finally we obtain the non-dimensional vorticity equation as shown in Eq. (8)
(8)
c) Transforming the energy equation If we expand the Del term, the energy equation will be as follows.
(9)
Substitute the non-dimensional variables in Eq. (4) into the Eq. (9)
(10)
5
⁄ Multiplying the R.H.S with equation.
, then we will get the final form of the non-dimensional energy
(11)
Where
Question 2 The non-dimensional governing equations are as follows. The ‘star’ superscript is intentionally not included for simplicity.
⃗ ⃗
(12)
(13)
(14)
Discretizing the governing equations
The index in this discretization is based on the figure below. 6
Figure 2: The control volumes used in the discretization
a) The stream (continuity) equation
(15)
From the previous homework on finite volume, we have developed in detail the L.H.S term of Eq. (15). Thus, here we will directly show the final form which was obtained from that homework.
( )
(16)
And
7
(17)
Thus the discretized equation is
( )
(18)
b) Discretize the vorticity equation
( ⃗)
Discretizing the first term on the L.H.S,
() () Using central difference at time level
(19)
∭
:
(20)
we obtain Eq. (21). We intentionally expand the
equation at this time level in order to be used in ADI scheme which method is by splitting a single time step into two of half size time step.
()
(21)
8
Discretizing the second term on L.H.S,
∭ ⃗
:
⃗ ⃗ ⃗
(22)
Since our problem is an incompressible flow then the divergence of velocity will be zero. Thus the equation becomes:
⃗ ⃗
(23)
Using the divergence theorem we will convert the volume integral on the R.H.S to become a surface integral.
⃗ ∬ ⃗ ⃗ ⃗ ̂ ̂
(24)
⃗
(25)
Evaluate the R.H.S of Eq. (25) one by one.
∮ ∫ ∫ ∫ ∫ Evaluate
:
(26)
9
{[ ] [ ]} [ ] ∮ Evaluate the
(27)
(28)
:
∫ ∫ ∫ ∫
(29)
10
[ ] [ { ]} [ ]
(30)
(31)
Substituting Eq. (28) and Eq. (31) into Eq. (25) we obtain:
⃗ [ ] [ ] ∭
(32)
Evaluate
Using the same method as that to find the previous volume integral of Laplacian we know that:
( ) ( )
(33)
11
Substitute Eq. (21), Eq. (32), and Eq. (33) into Eq. (19) and divide it by discretized equation for vorticity equation.
we will get the
[ ] [ ] ( )
(34)
We notice on the R.H.S that the space discretization has not been indexed with the time level. We will put the index based on the ADI method that we are going to apply when solving.
Step one for ADI:
[ ] [ ]
(35)
12
Step two for ADI:
[ ] [ ]
(36)
c) Discretize the energy equation
The differential form energy equation has exactly the same form as the vorticity equation. Thus,
the discretized equation will also be the same except that the and
are substituted by
and
, respectively. Please note, as mentioned before, all of the parameters here have been non-
dimensionalized. The discretized energy equation used for ADI scheme can be seen in Eq. (37) and Eq. (38). If we solve the vorticity and stream function equation decoupled from the energy equation, i.e., the temperature gradient does not affect the velocity gradient, and if our target is the steady state condition, then the transient terms in Eq. (37) and Eq. (38) will vanish.
The first step for ADI:
[ ] [ ]
(37)
13
The second step for ADI:
[ ] [ ]
(38)
Question 3 In this question we are asked to solve the flow problem presented numerically using ADI scheme with finite volume approach. However, after several attempts on solving the flow problem using ADI scheme with finite volume were not successful, another method was searched for. Explicit upwind scheme with finite difference approach was finally chosen to solve the set of equations. The form of the discretized governing equations remains the same ex cept for the boundary. The boundary condition for the control volume (finite volume) approach is as follows. Using the “Ghost cell” method together with the Taylor Series expansion we can relate the information on the boundary to the information of the cell nodes. Also, since the given boundary conditions are not explicitly written in terms of vorticity and stream function (instead, it is in terms of velocities), we need to convert/relate the velocity terms into the vorticity and stream function terms. For the finite difference approach, in this particular problem, creating “Ghost nodes” is not required. Only Taylor series expansion is used as well as the relation between the related
14
three parameters, i.e., velocity, stream function, and vorticity as it is also done in finite volume method.
Boundary condition treatment used
At the bottom wall
To find the boundary condition for stream function we will use its definition as follows.
(39)
The u will be zero because of no slip condition while the v will be zero due to the impermeable
wall. Thus,
and it follows that
will be a constant that can be chosen to be zero.
Boundary condition for vorticity is obtained by the following method.
It is also known by definition that
To evaluate
(40)
, we will use Taylor series expansion at the bottom wall.
15
() Rearranging the equation above we get
()
(41)
Based on the definition of stream function, we also know that
(42)
Substituting Eq. (42) and Eq. (41) into Eq. (40) and also knowing that for bottom wall, the stream function and u velocity is zero, we will obtain the boundary condition for vorticity at wall bottom as follows.
(43)
Using the grid notation (instead of distance, x) shown at the beginning of this report, we can write the B.C. for bottom wall as
(44)
Using the same procedure, boundary conditions on the other walls can be obtained. These boundary conditions can be seen inside the code that will be presented in this section.
16
The upwind explicit scheme is intentionally used to reduce the stability constraint on the time step which arises when central differencing scheme is used for the convective term. Moreover, upwind scheme also takes into account the effect of flow direction since information is travelling in the same direction with the flow. Hence, a more realistic result is expected to be obtained. To do this, upwind scheme is using backward differencing when the velocity is positive and forward differencing when the velocity is negative. Therefore, in the FORTRAN code, if-conditions are used to switch from one to another. The FORTRAN code used is presented at the end of this report.
Question 4 In this section we will compare our results with those from Ghia et al. As can be seen in Figure 3 and Figure 4 that the the current results match very closely to the published-Ghia et al.s results. Moreover, it can also be seen that the grid independence has been reached since the results at 128x128 are very close to that at 64x64. However, when higher Re is attempted, i.e. 400, it is found that the discrepancy is quite large. This might be because the different scheme used or due to the current simulation has not reached steady state yet as will be shown next. From Figure 5 it can be seen that, also for Re = 400, the higher the grid and the longer the flow time, the closer the current solution to that of Ghia’s et al. This result explains that higher Re requires finer mesh and longer simulation time to reach steady state. The current result also indicates that upwind explicit scheme can also predict the flow as good as the scheme used by Ghia et al.
17
1.00 0.90 0.80
* y , n 0.70 o i t c 0.60 e r i d y 0.50 s s e l n 0.40 o i s n e 0.30 m i D
Ghia et al., 128x128, Re=100 eldwin, 128x128, 6.5s, Re=100
0.20
eldwin, 64x64, 6s, Re=100
0.10 0.00 -0.30
-0.10
0.10
0.30
0.50
0.70
0.90
1.10
Dimensionless x-velocity, u*
Figure 3: Comparison of u velocity between the current results and the Ghia et al.’s results at Re=100 at the vertical bisection of the cavity
0.30
* v , y t i c o l e v y s s e l n o i s n e m i D
0.20
0.10
Dimensionless x-direction 0.00 0.00
0.10
0.20
0.30
0.40
0.50
0.60
0.70
0.80
0.90
1.00
-0.10
Ghia et al., 128x128, Re=100 eldwin, 64x64, 6s, Re=100 -0.20
eldwin, 128x128, 6.5s, Re=100
-0.30
Figure 4: Comparison of v velocity between the current results and the Ghia et al.’s results at Re=100 at the horizontal bisection of the cavity
18
140
120
* y , 100 n o i t c e r 80 i d y s s e 60 l n i o s n e 40 m i D
Ghia et al., 128x128, Re=400 eldwin, 64x64, 15.6s, Re=400 eldwin, 64x64, 6s, Re=400
20
eldwin, 128x128, 7.8s, Re=400 eldwin, 128x128, 20.8s, Re=400
0 -0.4
-0.2
0
0.2
0.4
0.6
0.8
1
1.2
Dimensionless x-velocity, u*
Figure 5: Comparison of u velocity between the current results and the Ghia et al.’s results at Re=400 at the vertical bisection of the cavity
Question 5 & 6 Plots of various parameters are presented in this section. In Figure 6, it can be seen that the maximum vorticity happens on the top-corner of the domain for Re=100. However, it is confined on a small region.
19
Figure 6: Vorticity contour plot for grid size 64x64, Re=100, and 6 sec. flow time.
For the stream function (Figure 7), the maximum value occurs between the center of the cavity and the moving lid and it is rather shifted to the right (the direction of the lid movement). The maximum (absolute) value is predicted around -0.1. The velocities contour and temperature contour are presented in Figure 8 Figure 10 below.
20
Figure 7: Stream function contour plot for grid size 64x64, Re=100, and 6 sec. flow time.
Figure 8: velocity at u-dir. contour plot for grid size 64x64, Re=100, and 6 sec. flow time.
21
Figure 9: velocity at v-dir. contour plot for grid size 64x64, Re=100, and 6 sec. flow time.
Figure 10: Temperature contour plot for grid size 64x64, Re=100, and 6 sec. flow time.
22
Question 7 The Nusselt number was calculated using Excel program based on the simulation data. Figure 11 shows the plot of Nusselt number along the horizontal walls.
Figure 11: Plot of Nusselt number along the top and the bottom boundary for grid size 128x128, Re=400, and 7.8 sec. flow time.
23
PROGRAM lid_driven_cavity_flow
!SOLVING THE LID DRIVEN CAVITY FLOW !USING EXPLICIT UPWIND METHOD
IMPLICIT NONE
INTEGER :: i,j,m,n,k,nts
!m is number of segments in x direction !n is number of segments in y direction !nts is the number of time step
REAL, PARAMETER :: Pr=0.7, Re = 400, Gr=0
DOUBLE PRECISION :: xsize, ysize, as, bs, cs, a, b, p, omega
!xsize is the nondimensional size of the domain in x dir (0
DOUBLE PRECISION :: dx, dy, dt
!dx is the size of the segment in x direction !dy is the size of the segment in y direction !dt is the time step size
DOUBLE PRECISION, ALLOCATABLE :: z(:,:),s(:,:),u(:,:),v(:,:),& &T(:,:),T_old(:,:)
!z is the dimensionless vorticity. !s is the dimensionless stream function 24
!T is the dimensionless temperature !It is designated to be !allocatable, i.e., 2D matrix with adjustable size.
WRITE(*,10,ADVANCE='no') READ(*,*) xsize WRITE(*,11,ADVANCE='no') READ(*,*) ysize WRITE(*,12,ADVANCE='no') READ(*,*) m WRITE(*,13,ADVANCE='no') READ(*,*) n WRITE(*,*) 'Number of time steps =' READ(*,*) nts WRITE(*,*) 'Relaxation factor for temperature calculation =' READ(*,*) omega
WRITE(*,*) 'The values that you have entered are:' WRITE(*,*) WRITE(*,*) 'x =', xsize WRITE(*,*) 'y =', ysize WRITE(*,*) 'm =', m WRITE(*,*) 'n =', n WRITE(*,*) 'Number of time steps =', nts
!Calculate dx, dy, dt
dx = dble(xsize)/dble(m) dy = dble(ysize)/dble(n) dt = dble(dx)/3.d0
25
OPEN (unit = 1 , file = "result_cavity_flow_explicit_upwind")
WRITE(1,*) 'The values that you use are :' WRITE(1,*) WRITE(1,*) 'x =', xsize WRITE(1,*) 'y =', ysize WRITE(1,*) 'm =', m WRITE(1,*) 'n =', n WRITE(1,*) 'time step size =', dt WRITE(1,*) 'the actual flow time =', nts*dt
WRITE(1,*)
ALLOCATE(z(1:m+1,1:n+1)) ALLOCATE(s(1:m+1,1:n+1)) ALLOCATE(u(1:m+1,1:n+1)) ALLOCATE(v(1:m+1,1:n+1))
!Initial condition
z(:,:) = 1.d-32 s(:,:) = 1.d-32 u(:,:) = 1.d-32 v(:,:) = 1.d-32
!Boundary condition
s(1,:)=0 ; s(:,1)=0; s(m+1,:)=0;s(:,n+1)=0; u(1,:)=0 ; u(:,1)=0; v(1,:)=0 ; v(:,1)=0; 26
v(m+1,:)=0;v(:,n+1)=0; u(m+1,:)=0; u(:,n+1)=1
!Calculate the vorticity distribution
DO k=1,nts DO j=2,n DO i=2,m
z(1,:)= -(2/dx/dx)*(s(2,:)) !left wall z(m+1,:)= -(2/dx/dx)*(s(m,:)) !right wall z(:,1)= -(2/dy/dy)*(s(:,2))!bottom wall z(:,n+1)= -(2/dy/dy)*(1.d0*dy + s(:,n)) !top wall
IF (u(i,j)>=0) THEN a = (z(i,j)-z(i-1,j))/dx ELSE a = (z(i+1,j)-z(i,j))/dx ENDIF
IF (v(i,j)>=0) THEN b = (z(i,j)-z(i,j-1))/dy ELSE b = (z(i,j+1)-z(i,j))/dy ENDIF
z(i,j)=dt*(-u(i,j)*(a)-v(i,j)*(b)& & + 1/Re*((z(i-1,j)-2*z(i,j)+z(i+1,j))/dx/dx& & +(z(i,j-1)-2*z(i,j)+z(i,j+1))/dy/dy))+z(i,j)
27
ENDDO ENDDO
!CALCULATE THE STREAM FUNCTION
as = 1.d0/(dy*dy) bs = 1.d0/(dx*dx) cs = -2.d0*(as + bs) !----------------------------------------------------------------------------------DO j=2,n DO i=2,m
s(i,j)=1/cs*(-as*s(i,j-1) - bs*s(i-1,j) - bs*s(i+1,j)& &- as*s(i,j+1) - z(i,j))
ENDDO ENDDO
!------------------------------------------------------------------------------
!FINDING THE VELOCITY U AND V
DO j=2,n DO i=2,m
u(i,j)= (s(i,j+1)-s(i,j-1))/2/dy
v(i,j) = -(s(i+1,j)-s(i-1,j))/2/dx
28
ENDDO ENDDO
ENDDO !ENDDO for the time step
!-----------------------------------------------------------
!CALCULATE THE TEMPERATURE DISTRIBUTION
ALLOCATE(T(1:m+1,1:n+1)) ALLOCATE(T_old(1:m+1,1:n+1))
!INITIALIZATION TO START THE ITERATION
T(:,:) = 0 !choose any random value
!BOUNDARY CONDITION
T(:,1) = 1 !Temp at the bottom T(:,n+1)= 0 !Temp at the top !T(1,:) = T(2,:) !USING FORWARD DIFFERENCE !T(m+1,:)= T(m,:) !USING FINTE DIFFERENCE
p = 1 !first guess for p (any number > 1e-5)
DO WHILE (p >= 1e-5)
T_old(:,:) = T(:,:)
DO j=2,n DO i=2,m 29
T(1,:) = T(2,:) !USING FORWARD DIFFERENCE T(m+1,:)= T(m,:) !USING FINTE DIFFERENCE T(i,j) = T(i,j) + omega*(1/cs*(-as*T(i,j-1) - bs*T(i-1,j) bs*T_old(i+1,j)& &- as*T_old(i,j+1) + Re*Pr*(u(i,j)*(T_old(i+1,j)-T(i-1,j))/2/dx +& & v(i,j)*(T_old(i,j+1)-T(i,j-1))/2/dy))-T(i,j)) T(1,:) = T(2,:) !USING FORWARD DIFFERENCE T(m+1,:)= T(m,:) !USING FINTE DIFFERENCE
ENDDO ENDDO
p = ABS(MAXVAL(T(2:n,2:m)-T_old(2:n,2:m))/MAXVAL((T(2:n,2:m))))
ENDDO !-------------------------------------------------------------------
!PRINTING RESULTS TO DATA FILE NO.1
WRITE(1,*)'The vorticity distribution:' WRITE(1,*) WRITE(1,*)'
i=1 i=2 i=3
...(so on if available)'
DO j=1,n+1
DO i=1,m+1 WRITE(1,21,ADVANCE='no') z(i,j) ENDDO WRITE(1,*) 30
ENDDO
WRITE(1,*)
WRITE(1,*)
!------------------------------------------------------------------------------WRITE(1,*)'The stream function distribution:' WRITE(1,*) WRITE(1,*)'
i=1 i=2 i=3
...(so on if available)'
DO j=1,n+1
DO i=1,m+1 WRITE(1,21,ADVANCE='no') s(i,j) ENDDO WRITE(1,*) ENDDO
WRITE(1,*) !--------------------------------------------------------------------------------------------
WRITE(1,*)'The u distribution:' WRITE(1,*) WRITE(1,*)'
i=1 i=2 i=3
...(so on if available)'
DO j=1,n+1
DO i=1,m+1 WRITE(1,21,ADVANCE='no') u(i,j) 31
ENDDO WRITE(1,*) ENDDO
WRITE(1,*)
!--------------------------------------------------------------------------------------
WRITE(1,*)'The v distribution:' WRITE(1,*) WRITE(1,*)'
i=1 i=2 i=3
...(so on if available)'
DO j=1,n+1
DO i=1,m+1 WRITE(1,21,ADVANCE='no') v(i,j) ENDDO WRITE(1,*) ENDDO
WRITE(1,*)
!-------------------------------------------------------------------
WRITE(1,*)'The temperature distribution:' WRITE(1,*) WRITE(1,*)'
i=1 i=2 i=3
...(so on if available)'
DO j=1,n+1
DO i=1,m+1 32
