Solution of 2-D Incompressible Navier Stokes Equations with Artificial Compressibility Method Using Ftcs Scheme

SOLUTION OF 2-D INCOMPRESSIBLE NAVIER STOKES EQUATIONS WITH ARTIFICIAL COMRESSIBILITY METHOD USING FTCS SCHEME IMRAN AZIZ Department of Mechanical Engineering College of EME National University of Science and Technology Islamabad, Pakistan Imran_9697@hotmail. com Abstract— The paper deals with the 2-D lid-driven cavity flow governed by the non dimensional incompressible Navier-Stokes theorem in the rectangular domain. Specific boundary conditions for this case study have been defined and the flow characteristics pertaining to the scenario have been coded in MATLAB using artificial compressibility method and FTCS scheme.

The results are compared successfully with an authentic research paper by Ghia, Ghia & Shin. Keywords: Navier stokes equations, Artificial Compressibility, FTCS scheme. Introduction The Navier-Stokes equations describe the motion of fluid substances and are used to solve wide range of problems in Fluid Dynamics. These equations include conservations of mass, momentum and energy. In this documentation we present the solution of 2-d navier stokes equations in a flow driven lid cavity used as a model for subsonic bombers weapons bay.

The weapons bay is a compartment on fighter aircrafts to carry bombs [1]. It is a highly critical area for any bomber design in the aeronautical industry. The design optimization of weapon bay involves various considerations namely weapon load capacity, extent of vibration and severe aerodynamic drag. In order to approximate the problem, a MATLAB code is developed to solve two dimensional incompressible navier stokes equations. Incompressible flows are those in which density variation is not linked to the pressure. The mass conservation is a constraint on velocity field.

The continuity equation can be combined with momentum equation to derive the equation for pressure or it can be solved independently by applying artificial compressibility. We have used the later approach by employing FTCS scheme. The results have been validated on a 81 ? 81 uniform grid and at Re = 1000 by comparing the u-velocity distribution along vertical centerline and v-velocity distribution along hoizontal centerline with the results obtained in various numerical studies (Ghia, Ghia & Shin in particular). A tolerance criteria of 1e-6 is used in this regard.

The results are then verified using Grid independent studies for the above mentioned parameters. After the validation and the consequent verification of the data obtained through the coded program, the horizontal and the vertical velocity contours span wise velocity contours and the velocity vector plots have been created to graphically ascertain the results. In the ensuing passage, the behavior of convergence with time and the analysis of the convergence on various grid sizes have been performed. The analysis has been complemented with a CPU time vs. grid sizes plot.

The tolerance criteria for convergence have then further been expedited and by using a semilog plot, the CPU time vs. tolerance graph is shown. Finally, the variation in the qualitative change in the flow structure; especially on the location of primary and secondary vortices, has been determined by giving a Re variation from 100-1000. Numerical mehodology The vector q = (u (x, y), v(x, y)) and pressure p(x,y) are used to describe completely the two dimensional flow field of an incompressible fluid. These functions are the solutions of the following mass and momentum conservation equations [2]. * Mass conservation ?. q=0

Or it can be written in explicit form as ?u? x+? v? y=0 In order to solve the system, an artificial compressibility is introduced in the above continuity equation as follows: ?p?? +1? 2? uj? xj=0 Where ? =1/? 2 is taken as the “artificial compressibility” of the fluid. The compressibility is a pseudo term and can be related to the pseudo-speed of sound and to an artificial density by using the following relations [3] ?=1/? 2 ?2=p/? All the quantities are taken in non dimensional form. When the time t approaches the infinity, the time dependent term will approach to zero and the numerical scheme will move toward steady state solution.

This means that on reaching the steady state solution the original continuity equation will be recovered [4]. Thus the steady state incompressible Navier stokes equations can be expressed in pseudo transient form as ? p? t +1/? 2[ ? u? x+? v? y ]=0 u: ? u? t+? u2? x+?? yuv=-? p? x+1Re? 2u? x2+? 2u? y2 v: ? v? t+?? xuv+? v2? y=-? p? y+1Re? 2v? x2+? 2v? y2 It is to be noted that all the quantities are in non dimensional form. The solution to above equations can be categorized based on grid system employed. Two grid systems are usually used, regular grid and staggered grid.

We have applied a uniform regular grid with same grid size in the x and y directions with the following boundary conditions on u and v velocities. U0,y=0 and V0,y=0 UxL,y=0 and VxL,y=0 Ux,0=0 and Vx,0=0 Ux,yL=U? and Vx,yL=0 Where Re=XLU? v, XL=1, YL=1 and U? =1 The application of correct boundary condition is extremely important and a slight mistake can change the whole nature of the problem. This is due to the reason that for any 2-d incompressible flow simulation, the equations remains the same, only the difference in boundary conditions change the nature of the problem. Numerical Scheme:

Numerous numerical schemes can be used to solve the above equations. We have employed explicit, forward in time central in space scheme (FTCS). We begin with the discretization of the modified mass conservation equation using a first order finite difference in time and a central scheme in space. The pressure component is solved first using the following discretized equation. Discretized P equation : pi,jn+1-pi,jn? t+1/? 2[ ui+1,jn-ui-1,jn2? x+vi,j+1n-vi,j-1n2? y ]= 0 Once the pressure component is computed for the time level n+1, the velocity equations can be solved by using the pressure value calculated above as

Discretized U equation: ui,jn+1-ui,jn? t+ui+1,j2n-(ui-1,j2)n2? x+pi+1,jn-pi-1,jn2? x + ui,j+1n*vi,j+1n-ui,j-1n*vi,j-1n2? y =(1/Re)[ ui+1,jn-2ui,j+ui-1,jn? x2+ui+1,jn-2ui,j+ui-1,jn? y2 ] Discretized V equation: vi,jn+1-vi,jn? t+ vi,j+12n-vi,j-12n2? y+pi,j+1n-pi,j-1n2? y +ui+1,jn*vi+1,jn-ui-1,jn*vi-1,jn2? x =(1/Re)[ vi+1,jn-2vi,j+vi-1,jn? x2+vi+1,jn-2vi,j+vi-1,jn? y2 ] Boundary conditions for Pressure: The pressure at the boundary conditions will come from the pressure at the internal points by using the equation approximated as Neuman boundary condition on all four walls.

For pressure, Neuman boundary conditions are applied at top, bottom and left, right walls. The normal derivative ? p? n is given by two points having boundary in the middle [5]. Pi,j+1-Pi,jdy=0 => Pi,j+1=Pi,j As the point Pi,j+1 is out of domain so it is approximated as Pi,j-1. Vorticity: The vorticity at a fluid point is defined as twice the angular velocity and is mathematically expressed as Vorticity(? )=2? = ?? V For two-dimensional flow, it reduced to ?z=? u? y-? v? x In descritized form ?z=Ui,j+1n-Ui,j-1 n2? y-Vi+1,jn-Vi-1,jn2? x…….. [6]

The vorticity equation is solved in Matlab in order to see the rotational effect of flow particles. Flow is rotational and show turning effect for non zero value of vorticity. When vorticity is zero, the flow becomes non-rotational and expressed in equation form as ?=?? V = 0. Primary vortices are the major vortices and occur along the streamlines of fluid flow. Secondary vortices are generated at a place from which the fluid departs already and vacuum is developed, that vacuum causes the flow to circulate back and hence generate secondary vortices. RESULTS 1.

Validation of Results for 81? 81 grid: Results for 81? 81grid with Re=1000 are shown in Figure 1. It contains contours and graphs at Re =1000, with convergence criteria =1e-6, No of iterations = 15000, and No of time steps =50000. Numerical values of velocities for 81×81 grid weren’t found so we made the tables using 129×129 grid at Re =1000. Fig 1a: U velocity distribution along Vertical Center line. Fig 1b: Ghia, Ghia &Shin’s U velocity plots for comparison [7] Fig 1c: V velocity distribution along Horizontal Center line Fig 1d: Results of V velocities from Ghia & Ghia[8]

Table-1a: RESULTS OF U VELOCITY ALONG VERTICAL LINE| Grid point number| Corodinate number (y)| U velocity| Ghia, Ghia n Shin[9]| 1| 0| 0| 0| 8| . 0547| -0. 023| -0. 42735| 9| 0. 0625| -0. 0243| -0. 42537| 10| 0. 0703| -0. 0256| -. 41657| 14| 0. 1016| -0. 032| -0. 38000| 23| 0. 1719| -0. 0377| -0. 32709| 37| 0. 2813| -0. 0492| -. 23186| 59| 0. 4531| -0. 0752| -0. 07540| 65| 0. 5000| -0. 0815| 0. 03111| 80| 0. 6172| -0. 0735| 0. 08344| 95| 0. 7344| -0. 0176| 0. 20673| 110| 0. 8516| 0. 0279| 0. 34635| 123| 0. 9531| 0. 2483| 0. 47804| 124| 0. 9609| 0. 3227| 0. 48070| 25| 0. 9688| 0. 4164| 0. 47783| 126| 0. 9766| 0. 5136| 0. 47221| 129| 1. 000| 1| 1. 0000| | | | | Table-1b: RESULTS OF V VELOCITY ALONG HORIZONTAL CENTER LINE| Grid point number| Corodinate number (x)| V velocity| Ghia, Ghia n Shin| 1| 0| 0| 0| 9| . 0547| 0. 0311| 0. 43983| 10| 0. 0625| 0. 0331| 0. 43733| 13| 0. 0703| 0. 0374| 0. 43124| 21| 0. 1016| 0. 0408| 0. 44187| 30| 0. 1719| 0. 0402| 0. 335070| 65| 0. 2813| 0. 0346| 0. 28003| 104| 0. 8047| -0. 1202| 0. 00831| 111| 0. 8594| -0. 1017| -0. 30719| 117| 0. 8516| -0. 0552| -0. 41496| 122| 0. 9453| -0. 0204| -0. 5863| 123| 0. 9531| -0. 0155| -0. 49099| 124| 0. 9609| -0. 0114| -0. 52987| 125| 0. 9688| -0. 0083| -0. 54302| 129| 1. 000| 0| 0| 2. Results for Grid Independent Studies Grid Independent studies are carried by using different grid sizes for Re=1000, No of Iterations =50000, time step =40000. The graphs show that as the grid size increases, the results also varies and after a certain grid size they come close to each other and further change in grid size has a negligible effect on behavior of the system at particular specifications. 61×61| | | 91×91| | | 121×121| | | Fig2.

Grid Independence 3. Contours and plots: Fig 3a: U velocity Contour Fig 3b: V velocity Contour Figure 3c: Vorticity Contour Figure 3d: Vector plot (quiver(U’,V’)) Fig 3e: Vector Plot zoomed 4: Behavior of Convergence with time: Figure 4: Error plot 5. Convergence behavior against Grid size: The time results obtained for various grid sizes from matlab are tabulated in Microsoft excel and then plotted as follow: Fig 5: Graph: Grid size(x-axis) vs CPU time (y axis) The behavior of convergence against grid sizes for fixed no. of Iterations at Re =1000 is tabulated as follows:

Grid Size| CPU Time| Convergence| 21| 53. 35204| 5. 41E-06| 41| 70. 28891| 7. 92E-06| 51| 91. 9678| 9. 55E-06| 61| 112. 4502| 9. 24E-06| 71| 137. 1834| 1. 07E-05| 81| 177. 0462| 1. 08E-05| 91| 219. 8227| 1. 15E-05| 101| 269. 3492| 1. 19E-05| 121| 352. 4425| 1. 27E-05| Table 2: Convergence behavior with CPU time vs Grid size 6. Tolerance criteria for convergence: Different tolerance criteria’s are analyzed by 81×81 grid size and a semi log plot is made vs CPU time in the Microsoft excel. CPU time is plotted on semi log scale on y axis. Fig. 4: Semi log plot.

Cpu Time| Tolerance Limit| 7114. 875| 1. 00E-06| 1104. 532| 5. 00E-06| 512. 435| 8. 00E-06| 330. 2863| 1. 00E-05| 151. 7562| 2. 00E-05| 37. 341| 5. 00E-05| Table 3: Cpu time vs Tolerance limit 7. Effect of variation of Reonald’s number: Below is given the table for Reonold’s number variation from 100-10000 and its effect on minimum horizontal velocity with corresponding coordinates. Re| Y cord| Min uu’| 100| 0. 75| -0. 1016| 200| 0. 8| -0. 1004| 300| 0. 825| -0. 0928| 400| 0. 85| -0. 0855| 500| 0. 85| -0. 0794| 600| 0. 8625| -0. 0743| 700| 0. 875| -0. 0699| 800| 0. 75| -0. 0659| 900| 0. 8875| -0. 0626| 1000| 0. 8875| -0. 0596| Table. 4: Re number variation vs min horizontal velocity DISCUSSION OF RESULTS Now we start to investigate the results one by one. 1. Validation of results for 81x 81 grids. Firstly we analyze Figure 1a which shows the behavior of Horizontal Velocity along the vertical center line. As the boundary condition for U at the bottom wall is 0 and is 1 at the top wall, the graph shows that U started from 0 at the bottom, becomes negative afterwards and then continuously increases to 1 due to the effect of top wall.

This shows that behavior of our graph is consistent with boundary conditions applied. The V velocity plot along the horizontal line shows approximately a sinusoidal curve . The velocity line starts from zero at the left wall, increases in the top left portion of the graph. As the graph move towards right, V decreases to zero again and then become negative in the downright portion. It goes to zero when it reaches the wall the right wall due to the influence of boundary condition. The values of U and V at x and y coordinates are tabulated in table 1a and 1b along with Ghia ; Shin results .

The comparison of the graphs and tabular values with those of Ghia;Shin, show that although the tabular values are very different, the overall net shape of the graphs is comparable, which help us assume that our results are validated. The difference in tabular values is due to the difference in use of Numerical scheme, convergence criteria, no of time steps and the processer. Ghia; Shin used vorticity and stream function formulation and implemented it with Full approximation scheme-full multi grid procedure where as we employed artificial compressibility using FTCS scheme. 2. Grid Independence.

Grid Independent studies are performed using a range of grids to observe that where the results become independent of grid spacing. There is a slight difference in V results in 61×61 and 91×91 grid where as the later is quite similar to 121×121 grid. The behavior of U is almost the same for all these grids. 3. Contours and Plots As boundary value at the top wall is fixed for U velocity =1, it makes the fluid particle streamlines near the wall resembles a boundary layer. Since U = 0 on the other three boundaries, effect of boundary condition of the top wall is filled in the interior region.

The circular portion in the center due to the turning operation of the streamlines shows that U is negative, as conventionally u is taken as positive in the right x axis direction and negative in the left axis direction. The V velocity contour shows that it is positive in the left half (especially in the top left half) and negative in the right half resulting in the strong circulation near the top wall. The U, V velocity counters and vector plot show that the curves are dragged towards the right hand side. This shows that most of the circulation occurs on the right side.

The shifting from center to left is because of the use of high Reonold number i. e. Re = 1000 for plotting. When we look at the quiver plot it shows that at the top V=0 , when the flow vectors go down from top right, the velocity of V increases in the negative direction, when they move to left and go upward the velocity V increases in the positive direction. This is in accordance with sign convention. As initially boundary condition for v is zero at all sides, because of the motion at the top region (due to U=1 boundary condition), the magnitude and direction of V has changed to represent the physical circulation of the fluid.

The physical circulation is clearly visible at the top, and fades as we go down. This asymmetry in the flow is due to the change of boundary condition as we go from top to bottom (U changes from U=1 to U =0). In Vorticity we are only plotting the component normal to the flow ie along the z axis as: ?z=? u? y-? v? xk Where, k represents the unit vector in the z direction [10]. As our flow is planner, it gives us vector normal to our flow plane. Contour plot shows that vortices are maximum in the top right corner due to the strong turning effects in the region. There are secondary vortices in the mid right and top left region.

In between these two vortices, the lines are close to a straight, so vorticity is very low that region. Quiver plot also shows that turning is much larger at the top right region and not much turn at the bottom portion. 4. Behavior of Convergence with time: The error graph in the fig 4 is not strictly monotone as strictly monotone is the one which always go down and not oscillate. Initially the error was large and till thousand iterations it was strictly monotone but after wards its behavior becomes oscillatory. When we reduce the delta T, it results in smooth convergence and also the error become monotonic.

More reduction gives finer results and everything become realistic i. e. streamlines go near the top wall. When the time approaches infinity, steady state will be reached and the disturbance will only occur near the top wall and behavior of U will become very smooth. Not much change in V results, it become more symmetric around the center. Error become completely monotone at t=20000 steps and fast and smooth convergence occurs. 5. Behavior of Convergence against Grid size The graph shows that when we decrease grid spacing by increasing grid points, CPU time increases.

But the advantage of increasing grid size is that it increases the convergence rate against the fix number of iterations. It can be seen from table 2 which shows convergence values against the grid size. In 41×41 grid elapsed time is less and it will converge faster as compared to 61×61 grid. We can also use small number of time steps for small grid size . The disadvantage of decreasing grid size is that results will become unclear. 6. Tolerance criteria for convergence Different tolerance criteria’s are studied and results are plotted vs CPU time.

It can be seen from the semi log plot that as we increase the tolerance criteria i. e. from 1e-6 to 5e-6 and so on the CPU time for convergence decreases. So increasing the tolerance means faster convergence at the cost of accuracy. 7. Effect of variation of Re number: Increase of Re implies decrease in viscosity. Increase of Re number, cause the streamlines to shift away from the walls and also move to the top region near the top wall. when Re increases, shifting of streamlines towards left and right, wall is visible because viscosity decreases and movement is easy for flow particles near the walls as compared to lower Re.

At Re =100, effect of viscosity was larger and velocity was close to zero near the walls due to tackiness, i. e. tackiness was not helping streamlines to move fast. So velocities were slow towards left and the top. Also as Re number increases the effect of viscosity on the bottom and left right walls decreases, especially u velocity is affected the most. U velocity in the top left and bottom right corners become visible. The bend in curve shows the transition is slow, and now the phenomena will occur near the top wall.

When viscosity decreases, wave like behavior increases and FTCS scheme become unstable. For very high Re number, viscosity will become close to zero, and this scheme will not work. The same effect is depicted from the error graph as the oscillations increases by increasing the Re. When Re increases, vorticities shift from center to the right corners. As Primary vortices occur along the streamlines of fluid flow, they move toward top right corner. Secondary vortices are generated at a place from which the fluid departs and vacuum causes the flow to circulate back and hence generate secondary vortices.

More secondary vortices are formed at high Re and they are located at the top left and bottom left and right corners of the graphs. Conclusion: A Matlab code for 2-d incompressible navier stokes equation with artificial compressibility is developed using FTCS scheme and the results are compared with a paper by Ghia, Ghia ; Shin. The values were not comparable because of difference in the solution schemes. The behaviors of Graphs follow the same pattern as given by various papers on internet, so the results are validated. Velocity and vortices counters, Grid independence and variation in tolerance studies are also carried out.

It is found that behavior of error is non-monotonic. Convergence with respect to time is studied and plotted on Microsoft excel chart. The effect of variation in Re number is also discussed in detail, it is found that by increasing Re vector plot shift towards top right, it effect convergence, results get detoriated and further increase make FTCS scheme inapplicable. References www. wekipedia. com Ionut Danaila, Pascal Jolly, “ An introduction to scientific computing, Springer, pp 251-2 Klaus A. Hoffmann, Steve T. Chiang “Computational Fluid Dynamics,” Engineering eduaction system,Wichita, Kansas, 67208-1078, USA, vol.

I, pp. 316–317,August 2001. Rodolfo salvi, the navier stokes equations, theory and numerical methods, pp 7 Benjamin Seibold, A compact and fast Matlab code,MIT, pp9. Klaus A. Hoffmann, Steve T. Chiang “Computational Fluid Dynamics,” Engineering eduaction system,Wichita, Kansas, 67208-1078, USA, vol. I, pp. 316–317, August 2001. U. Ghia, K. N. Ghia, And C. T. Shin, Journal of Computational Physis(1982),pp 396-397 U. Ghia, K. N. Ghia, And C. T. Shin, Journal of Computational Physis(1982),pp 396-397 U. Ghia, K. N. Ghia, And C. T. Shin, Journal of Computational Physis(1982),pp 396-397 Dale A. Anderson, John C.

Tannehill, Richard H. Pletcher,” Computational Fluid Mechanics and Heat Transfer”, Hemisphere publishing corporation, Newyork, H. P. 1984, pp 239 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Appendix % Matlab code for the Solution of 2-D incompressible Navier Stokes % equations using FTCS scheme with artificial compressibility. clc clear all format short Lx = 1. 0; % Length Ly = 1. 0;% Width t = 1; % time Re = 100; %Reonald number nx = 81; % No. of Space (x) steps ny = 81; nt = 50000; % No. of Time steps delT = t/(nt-1. ); % DelatT in Time delX = Lx/(nx-1. ; % DelatT in Space(x) delY = Ly/(ny-1. ); % DelatT in Space(x) bita = sqrt(0. 6); % Artifical compressibility factor %t = 0:delT:t; x = 0:delX:Lx; y = 0:delY:Ly; r1 = delT/(2. *delX); r2 = delT/(2. *delY); r3 = delT/delX^2 r4 = delT/delY^2; U(1:nx,1:ny)= 0; % Initialization of Solution matrix V(1:nx,1:ny)= 0; % Initialization of Solution matrix P(1:nx,1:ny)= 0. 5; % Initialization for pressure P1(1:nx,1:ny)= 0. 5; % Initial Condition of solution matrix. %Boundary Conditions U(:,ny)= 1. 0 ; %Boundary Condition on U err=1; %Error initilization iter = 1; tic; error = zeros(1,2); % Matrix for storing error. hile( err ; 1e-5) % for n = 1:nt % Loop for Time for i = 2:nx-1; % Loop for x for j = 2:ny-1; % Loop for y P1(i,j) = P(i,j)-(r1/bita^2)*(U(i+1,j)-U(i-1,j))-(r2/bita^2)*(V(i,j+1 )-V(i,j-1 )); U(i,j ) = U(i,j )-r1*(U(i+1,j )^2-U(i-1,j )^2)… -r2*(U(i,j+1 )*V(i,j+1 )-U(i,j-1 )*V(i,j-1 ))… -r1*(P1(i+1,j )-P1(i-1,j ))+(r3/Re)*(U(i+1,j )-2*U(i,j )+U(i-1,j ))… +(r4/Re)*(U(i,j+1 )-2*U(i,j )+U(i,j-1 )); V(i,j )= V(i,j )-r2*(V(i,j+1 )^2-V(i,j-1 )^2)… -r1*(U(i+1,j )*V(i+1,j )-U(i-1,j )*V(i-1,j ))… -r2*(P1(i,j+1 )-P1(i,j-1 ))… +(r4/Re)*(V(i,j+1 )-2*V(i,j )+V(i,j-1 ))… (r3/Re)*(V(i+1,j )-2*V(i,j )+V(i-1,j )); end end %% Neuman Boundary Conditions for i = 1:nx; P1(i,1) = P1(i,2); P1(i,ny) = P1(i,ny-1); end for j = 1:ny; P1(1,j) = P1(2,j); P1(nx,j) = P1(nx-1,j); end err=max(max(abs(P1-P))); error(1,iter) = err; P = P1 ; iter = iter + 1; end Vr = zeros(nx,ny); % Vorticity initilization for i = 2:nx-1; % Loop for x for j = 2: ny-1; % Loop for y Vr(i,j) = (0. 5*(V(i+1,j)-V(i-1,j))/delX)-(0. 5*(U(i,j+1)-U(i,j-1))/delY); end end toc %analysis quiver(U’,V’) figure(2); contourf(U’,20) title(‘U velociy Contours’) figure(3); contourf(V’,20) title(‘V velociy Contours’) igure(4); contourf(Vr’,20) title(‘Vorticity Contour Graph’) uu=U((nx-1)/2,:); figure(5); plot(uu’,y); title(‘U velocity profile on center point along verticle’); xlabel(‘U component of velocity’); ylabel(‘Grid points along Y-axis’); vv=V(:,(ny-1)/2); figure(6); plot(x,vv’) title(‘V velocity profile on center point along horizontal’); xlabel(‘Grid points along X-axis’); ylabel(‘V component of velociy’); figure(7); plot(error) title(‘Error Graph’); [A I]=min(U); [A1 I1]=min(a) min(min(uu’)) figure(9); semilogx(error) title(‘Semilog plot x’); figure (10); semilogy(error) title(‘Semilog plot y’); % END