Lid Driven Cavity - Post 2

(You can find the PDF version of the below blogpost here. )
(Citation : R, AAKASH. "Lid Driven Cavity." Classroom fundaes. Aakash R, 2 May 2017. Web. http://instafundae.blogspot.in/2017/05/lid-driven-cavity-post-2.html.)

Abstract

The lid driven cavity is a standard benchmark problem for any CFD software/Code, Lid driven cavity flows are important in many industrial processing applications such as short-dwell and flexible blade coaters.This problem consists of a 2-D cavity , whose top side moves with  a uniform velocity $u_0$.There are numerous results form the literature for the above problem , which have been used to benchmark the code,its respective results presented in this paper.The benchmarking references used in this paper are A. AbdelMigid [1] ,O. BOTELLA [3],Ghia $et$ $al.$[4]
This paper presents steady state solutions to Modified Lid Driven Cavity(Same as lid driven cavity problem but with inflow and outflow at the left and right hand side boundary) introduced in the book by Hoffman [5]  for low to Middle  $Re$s which can be extended to produce results for the Lid Driven Cavity . It is to be duly noted that the results presented in this paper are just a small part of what the Code (The Matlab code written to solve the Modified LDC) has been tested for,the same code can be used for middle to high level $Re$s too with appropriate grid sizes.
The approach used in solving the above mentioned problem is the standard stream vorticity approach . The numerical schemes used for producing results in this include central ,forward ,backward differencing schemes,Guass seidel PSOR(Point Successive Over Relaxation) methods in solving equation,etc.

The schematic diagram of the Modified Lid Driven Cavity is as follows :

Governing Equations

The governing equations of the stream vorticity approach are Vorticity transport equation,stream function  equation , continuity equations which are coupled .The below equations can be easily derived .
\begin{eqnarray}
\frac{\delta \omega}{\delta t} + u \frac{\delta \omega}{\delta x} +v\frac{\delta \omega}{\delta y} &=& \gamma \left(\frac{\delta ^2 \omega}{\delta x ^2} + \frac{\delta ^2 \omega}{\delta y ^2}\right)\\
  \nabla ^2\psi&=& - \omega\\
u &=&  \frac{\delta \psi}{\delta y}\\
v &=& -\frac{\delta \psi}{\delta x}
\end{eqnarray}
Note the continuity equation is implicitly used in deriving the stream vorticity formulations.

Boundary Conditions
Modified Lid Driven Cavity

Boundary conditions for $u,v,\psi$

Bottom side :
\begin{eqnarray}
u &=& 0\\
v &=& 0\\
\psi &=& c_1
\end{eqnarray}
Top side :
\begin{eqnarray}
u &=& 1\\
v &=& 0\\
\psi &=& c_2
\end{eqnarray}
Right side:
For $1:(j_3 -1)$
\begin{eqnarray}
u &=& 0\\
v &=& 0\\
\psi &=& c_1
\end{eqnarray}
 For $j_3:j_4$
\begin{eqnarray}
\frac{\delta \psi}{\delta x} &=& 0
\end{eqnarray}
 For $(j4+1):end$
\begin{eqnarray}
u &=& 0\\
v &=& 0\\
\psi &=& c_2
\end{eqnarray}

Left side :
 For $1:(j_1 -1)$
\begin{eqnarray}
u &=& 0\\
v &=& 0\\
\psi &=& c_1
\end{eqnarray}
 For $j_1:j_2$
\begin{eqnarray}
\frac{\delta \psi}{\delta x} &=& 0
\end{eqnarray}
 For $(j2+1):end$
\begin{eqnarray}
u &=& 0\\
v &=& 0\\
\psi &=& c_2
\end{eqnarray}

Boundary conditions for $\omega$

\begin{equation}
  \omega  = -\left(  \frac{\delta ^2 \psi}{\delta x^2} + \frac{\delta ^2 \psi}{\delta y^2}\right)
\end{equation}

Discretized Equations

Vorticity Transport Equation

Inner Grid Formulation
Central Differencing Scheme
\begin{equation}
u(i,j)\frac{\omega(i+1,j)-\omega(i-1,j) }{2\Delta x}+v(i,j)\frac{\omega(i,j+1)-\omega(i,j-1) }{2\Delta y} =  \gamma \left( \frac{\omega(i+1,j)-2\omega(i,j)+\omega(i-1,j)}{(\Delta x)^2} + \frac{\omega(i,j+1)-2\omega(i,j)+\omega(i,j-1)}{(\Delta y)^2}\right)
\end{equation}
Rearranging , equating $\Delta x = \Delta y$,we get:
\begin{equation}
\omega(i,j) = \frac{1}{4} \left( \omega(i+1,j)\left(1 - \frac{u(i,j)\Delta x}{2\gamma}\right) + \omega(i-1,j)\left(1 + \frac{u(i,j)\Delta x }{2\gamma}\right) \\ + \omega(i,j+1)\left(1 - \frac{v(i,j)\Delta x }{2\gamma}\right) +  \omega(i,j-1)\left(1 + \frac{v(i,j)\Delta x}{2\gamma}\right)\right)
\end{equation}
Finally applying the PSOR
\begin{equation}
\omega(i,j) =(1-\alpha)\omega(i,j) +  \frac{\alpha}{4} \left( \omega(i+1,j)\left(1 - \frac{u(i,j)\Delta x}{2\gamma}\right) + \omega(i-1,j)\left(1 + \frac{u(i,j)\Delta x }{2\gamma}\right) \\ + \omega(i,j+1)\left(1 - \frac{v(i,j)\Delta x }{2\gamma}\right) +  \omega(i,j-1)\left(1 + \frac{v(i,j)\Delta x}{2\gamma}\right)\right)
\end{equation}
Boundary Formulation

Left Boundary 
\begin{eqnarray}
\omega(1,j )&=& - \left. \frac{\delta ^2 \psi}{\delta x^2} \right \vert _{(1,j)} - \left.\frac{\delta u}{\delta y}\right \vert _{(1,j)}\\
\implies \omega(1,j) &=& 2\frac{\psi(1,j)-\psi(2,j)}{\Delta x^2} - 2\frac{v(1,j)}{\Delta x} -\frac{u(1,j+1)-u(1,j-1)}{(2\Delta y)}
\end{eqnarray}
Right Boundary 
\begin{eqnarray}
\omega(im,j )&=& - \left. \frac{\delta ^2 \psi}{\delta x^2} \right \vert _{(1,j)} - \left.\frac{\delta u}{\delta y}\right \vert _{(1,j)}\\
\implies \omega(1,j) &=& 2\frac{\psi(1,j)-\psi(2,j)}{\Delta x^2} - 2\frac{v(1,j)}{\Delta x} -\frac{u(1,j+1)-u(1,j-1)}{(2\Delta y)}
\end{eqnarray}

Stream Vorticity Equation

Inner Grid Formulation
Central Differencing Scheme
\begin{eqnarray}
\omega = - \left( \frac{\psi(i+1,j)-2\omega(i,j)+\psi(i-1,j)}{(\Delta x)^2} + \frac{\psi(i,j+1)-2\psi(i,j)+\psi(i,j-1)}{(\Delta y)^2}\right)\\
\implies \psi(i,j) = \frac{1}{4}\left(\psi(i,j-1)+\psi(i,j+1)+\psi(i-1,j)+\psi(i+1,j)+ \omega (i,j)\Delta x ^2 \right)
\end{eqnarray}
Finally applying PSOR we get
\begin{equation}
 \psi(i,j) = (1-\alpha)\psi(i,j)+\frac{\alpha}{4}(\psi(i,j-1)+\psi(i,j+1)+\psi(i-1,j)+\psi(i+1,j)+ \omega (i,j)\Delta x ^2 )
\end{equation}

Boundary Formulation

Left Boundary ($j_1 : j_2$):
\begin{eqnarray}
\left. \frac{\delta \psi}{\delta x} \right \vert _{(1,j)} = v
\end{eqnarray}
Using 2nd order accurate scheme we get
\begin{eqnarray}
\psi(1,j) = \frac{4\psi(2,j)-\psi(3,j) + 2\Delta x v(1,j)}{3}
\end{eqnarray}
 Right Boundary ($j_3 : j_4$):
Using 2nd order accurate scheme we get
\begin{eqnarray}
\psi(im,j) = \frac{4\psi(im-1,j)-\psi(im-2,j) - 2\Delta x v(1,j)}{3}
\end{eqnarray}

Velocity Equations

Inner Grid Formulation

Central Differencing Scheme
\begin{eqnarray}
u = \frac{\psi(i,j+1)-\psi(i,j-1)}{(2\Delta  y)}\\
v = \frac{\psi(i+1,j)-\psi(i-1,j)}{(2\Delta  x)}
\end{eqnarray}

Boundary Formulation

Forward Differencing Scheme
\begin{eqnarray}
u &=& \frac{(\psi(j+1,i) - \psi(j,i))}{\Delta y}\\
v &=& 0
\end{eqnarray}

Parameters used for convergence measure
Before we see the algorithm used in solving the problem ,I would introduce few parameters used in checking the convergence of the solution.
The parameters are as follows :

1.  $e_1$: RMSE $\left ( \frac{\psi  ^{k+1} _{i,j} -\psi  ^k _{i,j}}{\psi  ^k _{i,j}} \right )$ (where k is the iteration number)
   Variable name used in the program ("err").
2.  $e_2$ :RMSE$\left ( \frac{\omega  ^{k+1} _{i,j} -\omega  ^k _{i,j}}{\omega  ^k _{i,j}}\right )$ .
   Variable name used in the program ("err1")
3.  $e_3$ : Numerical solution accuracy of Vorticity equation .
   Variable name used in the program ("ERR1")
4. $e_4$ : Numerical solution accuracy of Navier stokes equation .
   Variable name used in the program ("ERR2")

It is to be noted that $e_1<e_2<e_3<e_4$ for the algorithm used in this paper , hence convergence is measured using $e_4$ , divergence measured using $e_1$.

Algorithm

1. Assuming a specific grid size $N_x,N_y $  .
2.  Construct matrices for  $\psi,\omega,u,v$  using the above mentioned size using random values .  
3. Assume appropriate boundary conditions for  $\psi,\omega,u,v$ .
4. Solve for  $\psi$  using the PSOR method using Stream Vorticity function .
5.  Compute  $u,v$  using central difference scheme .
6. Iterate the PSOR method for  $\omega$   using vorticity transport Equation once (using the pre solved  $\psi,u,v $ .
   (Note : Don't solve for  $\omega$ run the PSOR over the entire grid once.)
7. Compute necessary Error Functions, check for convergence.
8. If converged , STOP. Else go to Step 4.

Results

Presented below are few results which can be used for used for benchmarking purposes . 

129x129 Grid 

S.No	U Velocity  (m/s)					V velocity  (m/s)
Distance (m)	Re 100	Re 400	Re 800	Re 1000	Re 3200	Re 100	Re 400	Re 800	Re 1000	Re 3200
1.00	1.0000	1.0000	1.0000	1.0000	1.0000	0.0000	0.0000	0.0000	0.0000	0.0000
0.9760	0.8431	0.7568	0.6799	0.6537	0.5189	-0.0465	-0.0885	-0.1349	-0.1655	-0.3331
0.9687	0.7911	0.6825	0.5955	0.569	0.4747	-0.0621	-0.1224	-0.1885	-0.2311	-0.4441
0.9609	0.7402	0.6155	0.5283	0.5068	0.4609	-0.0777	-0.1577	-0.2434	-0.2967	-0.5218
0.9531	0.6906	0.5567	0.4775	0.4632	0.4583	-0.0930	-0.1936	-0.2974	-0.3586	-0.5537
0.8515	0.2357	0.2883	0.3211	0.3350	0.3465	-0.2388	-0.4452	-0.4241	-0.4061	-0.3679
0.7344	0.0037	0.1603	0.1806	0.1892	0.2000	-0.2134	-0.2522	-0.2295	-0.2356	-0.2300
0.6172	-0.1388	0.0198	0.0483	0.0551	0.0775	-0.0738	-0.0900	-0.0981	-0.0999	-0.1057
0.5000	-0.2087	-0.1159	-0.0705	-0.0620	-0.0367	0.0574	0.0526	0.0295	0.0258	0.0143
0.4531	-0.2134	-0.1722	-0.1161	-0.1089	-0.0800	0.0970	0.1110	0.0798	0.0771	0.0605
0.2813	-0.1570	-0.3244	-0.2982	-0.2808	-0.2405	0.1752	0.2855	0.2758	0.2690	0.2361
0.1719	-0.1013	-0.2376	-0.3564	-0.3852	-0.3434	0.1696	0.2862	0.3515	0.3715	0.3556
0.1016	-0.0641	-0.1418	-0.2456	-0.2923	-0.4297	0.1323	0.2357	0.3070	0.3370	0.4291
0.0703	-0.0464	-0.1001	-0.1783	-0.2236	-0.4036	0.1031	0.1947	0.2625	0.2958	0.4071
0.0625	-0.0418	-0.0897	-0.1610	-0.2029	-0.3798	0.0944	0.1812	0.2475	0.2804	0.3914
0.0547	-0.0371	-0.0792	-0.1435	-0.1817	-0.3499	0.0850	0.1661	0.2302	0.2625	0.3719
0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000

Comparison of my code with Ghia et al (129x129 Grid) 

Grid Points	U-Velocity		Grid Points	V-Velocity
	Present	Ghia et al		Present	Ghia et al
	Re 100			Re 100
129	1	1	129	0	0
126	0.84311	0.84123	125	-0.0621	-0.05906
125	0.79114	0.78871	124	-0.07766	-0.07391
124	0.74015	0.73722	123	-0.09305	-0.08864
123	0.69055	0.68717	122	-0.10817	-0.10313
110	0.23571	0.23151	117	-0.17657	-0.16914
95	0.00375	0.00332	111	-0.23289	-0.22445
80	-0.1388	-0.13641	104	-0.2526	-0.24533
65	-0.2087	-0.20581	65	0.057388	0.05454
59	-0.2134	-0.2109	31	0.178879	0.17527
37	-0.157	-0.15662	30	0.178671	0.17507
23	-0.1013	-0.1015	21	0.164132	0.16077
14	-0.0641	-0.06434	13	0.125833	0.12317
10	-0.0464	-0.04775	11	0.11128	0.1089
9	-0.0418	-0.04192	10	0.103134	0.10091
8	-0.0371	-0.03717	9	0.094379	0.09233
1	0	0	1	0	0

251x251 Grid

S.No	U Velocity  (m/s)				V velocity  (m/s)
Distance (m)	Re 100	Re 400	Re 1000	Re 3200	Re 100	Re 400	Re 1000	Re 3200
1.00	1.0000	1.0000	1.0000	1.0000	0.0000	0.0000	0.0000	0.0000
0.9760	0.8431	0.7568	0.6537	0.5189	-0.0465	-0.0885	-0.1655	-0.3331
0.9687	0.7911	0.6825	0.569	0.4747	-0.0621	-0.1224	-0.2311	-0.4441
0.9609	0.7402	0.6155	0.5068	0.4609	-0.0777	-0.1577	-0.2967	-0.5218
0.9531	0.6906	0.5567	0.4632	0.4583	-0.0930	-0.1936	-0.3586	-0.5537
0.8515	0.2357	0.2883	0.3350	0.3465	-0.2388	-0.4452	-0.4061	-0.3679
0.7344	0.0037	0.1603	0.1892	0.2000	-0.2134	-0.2522	-0.2356	-0.2300
0.6172	-0.1388	0.0198	0.0551	0.0775	-0.0738	-0.0900	-0.0999	-0.1057
0.5000	-0.2087	-0.1159	-0.0620	-0.0367	0.0574	0.0526	0.0258	0.0143
0.4531	-0.2134	-0.1722	-0.1089	-0.0800	0.0970	0.1110	0.0771	0.0605
0.2813	-0.1570	-0.3244	-0.2808	-0.2405	0.1752	0.2855	0.2690	0.2361
0.1719	-0.1013	-0.2376	-0.3852	-0.3434	0.1696	0.2862	0.3715	0.3556
0.1016	-0.0641	-0.1418	-0.2923	-0.4297	0.1323	0.2357	0.3370	0.4291
0.0703	-0.0464	-0.1001	-0.2236	-0.4036	0.1031	0.1947	0.2958	0.4071
0.0625	-0.0418	-0.0897	-0.2029	-0.3798	0.0944	0.1812	0.2804	0.3914
0.0547	-0.0371	-0.0792	-0.1817	-0.3499	0.0850	0.1661	0.2625	0.3719
0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000

Comparison with Ghia etal ,and O .BOTELLA:

Grid Points	U-Velocity					Grid Points	V-Velocity
	Present	Ghia et al	Present	Ghia et al	O.BOTELLA		Present	Ghia et al	Present	Ghia et al	O.BOTELLA
	Re 400		Re 1000				Re 400		Re 1000
251	1	1	1	1	1	251	0	0	0	0	0
245	0.75389422	0.75837	0.653659036	0.65928	0.6644227	243	-0.127997501	-0.12146	-0.2311	-0.21388	-0.22792
243	0.67893746	0.68439	0.569419488	0.57492	0.5808359	241	-0.16480199	-0.15663	-0.2967	-0.27669	-0.29369
241	0.61153269	0.61756	0.506788795	0.51117	0.5169277	239	-0.202194664	-0.19254	-0.3586	-0.33714	-0.35532
239	0.5527114	0.55892	0.463245781	0.46604	0.4723329	237	-0.239373072	-0.22847	-0.4134	-0.39188	-0.41038
214	0.29147429	0.29093	0.334964485	0.33304	0.3372212	228	-0.382661909	-0.23827	-0.5226	-0.5155	-0.52644
185	0.16391904	0.16256	0.189150105	0.18719	0.1886747	216	-0.452317228	-0.44993	-0.426	-0.42665	-0.42645
155	0.01930732	0.02135	0.055104972	0.05702	0.0570178	202	-0.383530806	-0.38598	-0.3172	-0.31966	-0.32021
126	-0.1152734	-0.11477	-0.061988317	-0.0608	-0.062056	126	0.052196894	0.05188	0.02584	0.02526	0.0258
114	-0.1730285	-0.17119	-0.108910482	-0.10648	-0.1082	60	0.302110802	0.30174	0.32157	0.32235	0.32536
71	-0.3276167	-0.32726	-0.280799838	-0.27805	-0.28037	58	0.30257526	0.30203	0.33037	0.33075	0.33399
44	-0.2423208	-0.24299	-0.385150406	-0.38289	-0.388569	40	0.282024218	0.28124	0.37335	0.37095	0.37692
26	-0.1428726	-0.14612	-0.292346467	-0.2973	-0.300456	24	0.228366656	0.22965	0.32684	0.32627	0.33304
19	-0.1047669	-0.10338	-0.223618183	-0.2222	-0.222896	21	0.212590596	0.2092	0.30931	0.30353	0.30991
17	-0.0938794	-0.09266	-0.202882065	-0.20196	-0.20233	19	0.200473862	0.19713	0.29582	0.29012	0.29627
15	-0.0829239	-0.08186	-0.181688508	-0.18109	-0.181288	17	0.1868535	0.1836	0.28043	0.27485	0.28071
1	0	0	0	0	0	1	0	0	0	0	0

501x501 Grid:

Note : The Underlined values show some discrepancies with the results of Ghia etal.

Grid Points	U-Velocity			Grid Points	V-Velocity
	Present	Ghia et al	AbdelMigid		Present	Ghia et al	AbdelMigid
	Re 3200	Re 3200	Re 3200		Re 3200	Re 3200	Re 3200
501	1	1	1	501	0	0	0
489	0.518868	0.53236	0.536	485	-0.44412	-0.39017	−0.4106
485	0.474687	0.48296	0.4886	481	-0.52184	-0.47425	−0.5049
481	0.460908	0.46547	0.4661	478	-0.55368	-0.52357	−0.5560
478	0.458339	0.46101	0.461	474	-0.56399	-0.54053	−0.5671
427	0.346469	0.34682	0.3491	454	-0.44061	-0.44307	−0.4462
368	0.200026	0.19791	0.2035	431	-0.37699	-0.37401	−0.3799
310	0.077457	0.07156	0.078	403	-0.31086	-0.31184	−0.3148
251	-0.03674	-0.04272	−0.0369	251	0.014275	0.00999	0.0142
228	-0.07997	-0.86636	−0.0808	118	0.286802	0.28188	0.2864
142	-0.24047	-0.24427	−0.2414	114	0.295331	0.2903	0.2954
87	-0.34339	-0.34323	−0.3445	79	0.37531	0.37119	0.3753
52	-0.42972	-0.41933	−0.4314	48	0.4294	0.42768	0.4318
36	-0.40356	-0.37827	−0.4081	40	0.418656	0.41906	0.4224
32	-0.37977	-0.35344	−0.3839	36	0.407058	0.40917	0.4109
28	-0.34989	-0.32407	−0.3530	32	0.391374	0.3956	0.3949
1	0	0	0	1	0	0	0

Matlab Code URL

Link : https://github.com/RAAKASH/Intro-to-CFD-/tree/master/LDC
Check LDC - Post 3 for description of the Matlab code .

References

[1] Tamer A. AbdelMigid, Khalid M. Saqr, Mohamed A. Kotb, and Ahmed A. Aboelfarag. Revisiting the lid-driven cavity flow problem: Review and new steady state benchmarking results using {GPU} accelerated code. Alexandria Engineering Journal, 56(1):123 – 135, 2017.

[2] C. K. Aidun, N. G. Triantafillopoulos, and J. D. Benson. Global stability of a liddriven cavity with throughflow: Flow visualization studies. Physics of Fluids A: Fluid Dynamics, 3(9):2081–2091, 1991.

[3] O Botella and R Peyret. Benchmark spectral results on the lid-driven cavity flow. Computers & Fluids, 27(4):421–433, 1998.

[4] UKNG Ghia, Kirti N Ghia, and CT Shin. High-re solutions for incompressible flow using the navier-stokes equations and a multigrid method. Journal of computational physics, 48(3):387–411, 1982.

[5] K.A. Hoffman. Computational Fluid Dynamics for Engineers. Number v. 2. Engineering Education System, 1989.