Table of Contents Introduction3 Development of numerical schemes5 Partial Differential Equations5 Initial and Boundary condition5 Modelling Approaches6 Numerical Methods6 Explicit method8 Implicit method8 Numerical Coding10 Explicit method10 Final code11 Implicit Method15 Final Code16 Numerical results18 Analysis of the Numerical results23 Conclusion24 References25 Introduction Over the years the importance of fluid dynamics has grown exponentially.
It represents the theoretical and physical aspects of the fluid in motion, as it flows naturally or when effected by a force.This application can be applied to liquids and gases providing a deeper understanding of pure sciences such as atmospheric, geophysics and oceanography. However the main application of this study is in industries such as turbo machinery, aerospace, civil engineering, automotive, water industry and more as these are realising the importance of this field.
Over the last century there were two different approaches used in this field. The theoretical-analytical approach requires the uses of partial differential equations which consist of continuity, Euler and Navier-Stokes equations.These helped to study and understand the behaviour of the fluid as it flows and predict and changes which may occur. However there are some issues which occur in mathematical analysis, calculating the data can be extremely difficult complicated therefore assumptions are made in order to simplify these equations. Due to this the results produced are not completely accurate and contain errors. On the other hand experimental approaches were used to physically model the flow in a test or a lab to understand the behaviour.
It can show different aspects of the flow during separation or the transition of laminar to turbulent, but as the data is only qualitative it doesn’t provide any information on velocity or pressure of the fluid (J. M. McDonough, 2009).
(Anderson. Jr, 1995) As the technology advanced it allowed the theoretical and experimental approach to be combined together though the use of computer simulation techniques, this is known as computational fluid dynamics. This uses mathematical equations, numerical methods and software to produce quantitative and qualitative data for the flow.The advantages of CFD compared to the previous approach are shown in the table blow (Kuzmin, D ). Experiment| CFD| * Calculations made for a signal quantity * Model size limited to test or laboratory size * Range of problems are limited * Expensive to carry out experiments * Slow and extremely time consuming | * For all the quantities at the same time * Can represent the actual flow area * Can be applied to any problem and conditions * Relatively cheap to model flows * Fast and can complete many iterations |The report shows a detailed study of the different methods used in CFD, the process or creating a program and analysis of the results gained. This will help to gain deeper understanding of how fluids work and how CFD can be applied to different scenarios to predict the flow outcome. Aim The purpose of this report is to create a model of a couette flow and analyse the results.
This consists of a vicious fluid between two parallel plates. As the upper plate moves and the lower plate remain stationary and the fluid in the middle flows due to the shear stress created by the upper plate.This creates a linear velocity across the flow, which can be calculated by using the Navier-Stokes equation and Finite Difference Method (FDM).
Objectives * A detail research on physical and theoretical aspects of the CFD * To develop a numerical schemes * Numerical coding using a programming language such as MatLab * Graphic output of the numerical results * Analysis of numerical results Development of numerical schemes Partial Differential Equations Fluid dynamics is bases on the calculation of hydrodynamic equations as these represent the behaviour of the fluid.These consist of partial differential equations which represent the variation of the same point over a range of time. Continuity equation is one of the main equations which is used in the CFD, as show below it calculates the velocity U,V and W of the point in a three dimensional space X,Y Z. This can be altered and applied to the incompressible and compressible flows. (Zuo, W. 2005) Euler’s equation is another main equation used for inviscid fluids, which is the flow of an ideal fluid with no viscosity. This allows simplification of problems which have fluids with low viscosity making the problem easier to solve, this is shown below.
However Navier-Stokes equation is the main one used in CFD during calculation. This is derived in the same process as the Euler’s equation but has additional shear stresses from Newton’s second law. Using this Momentum equation which shows the acceleration and continuity equations which shows the mass can be derived, these are shown bellow. Combing these two produces four equations, using numerical methods these can be used to determine the four unknowns which are the pressure (P) and the velocity in three direction( U, V and W). (CHEM 520)Initial and Boundary condition During computation the limits of the model need to be defined. These include the initial conditions which determine the starting point which include the time as 0, velocity as 0 and the pressure as 0. The boundary conditions on the other hand remain the same throughout the computation.
For the Navier-Stokes equation these can be simplified into three types which are Hyperbolic, Parabolic and Elliptic. All three of these have different solving methods as two of them have a marching direction and the Elliptic type doesn’t.However for this report only parabolic type will be used which is determined by the equation shown bellow. Where (T) is used represent the variable such as temperature or velocity, (x) is the direction or displacement to the next point and (t) is the time step taken. What is alpha Modelling Approaches For modelling techniques it is important to first understand the characteristics of the fluid, to see how it would behave during motion, this is determined by the Reynolds number which is the ratio of inertial to viscous forces in the flow.
Laminar flow represents the flow in a smooth motion where the Reynolds number is blow 2000, this usually occurs in slow moving fluids. On the other hand if the Reynolds number is above 4000 it represent a turbulent flow. At this point the motion of the fluid is fast, highly irregular and constantly changing directions. It produces eddy currents in the fluid which causes it to mixing of the fluid particles. In most practical cases the flow is mainly turbulent so the modelling techniques are usually based on this.Direct numerical Method (DNS) is one of the main ones which provides an accurate result for the Navier-Stokes equation by not making any assumptions. Due to this the computation of this can become extremely complicated and expensive therefore it is usually limited to problems which a small simulation or the Reynolds number.
LARGE Eddy Simulation (LES) is another main one which is used when the problem is too complicated for DNS approach, it usually deals with large simulations with a higher Reynolds number. Reynolds Averaged Navier?Stokes (RANS) is another method which uses the time as a function to average out the steady and unsteady flows. As LES and DNS provide instantaneous results with excessive detail which may not be necessary, in these cases the RANS approach is used which provides an average flow properties. Numerical Methods Solving the partial differential equations results in a answer which describes the whole area of the problem, however using numerical methods it is possible to transform the equation to gain the answer at a certain point in the area.These include methods such as Finite volume (FVM) which provides a fast solution for large problems with turbulent flow, this is known to be the main one used due to its broad applicability. Finite element method (FEM) which is mainly used for structural analysis but can be modified for fluid analysis, it produces a conservative solution but can have a higher processing time. Finite difference method (FDM) provides a simple accurate solution and can be used for the complex problems and geometries.
(Bakker. A, 2002).Finite difference method is based on the Taylor series expansion.
It uses a grid to define the different points, it transforms the partial different equation into an algebraic equation to solve for each of these points. Every problem has a size or domain which is limited by the boundary conditions, this domain can be broken down into grid points as shown in the figure below. The Y direction represents the time interval to determine the amount of iterations and the X direction shows the difference in space as the fluid moves form one point to the next.The algebraic equation derived from the Taylor series expansion can be used to calculate a solution for each of these points. In this method a centre point is used to calculate the pervious point and the one after, this is known as the first order backwards difference and first order forward difference. The information from both of the side points is used to determine the centre point known as the second order central difference. (Urroz.
E. G, 2004) The model equation is a parabolic type with regards to time, therefore it require the first derivative as shown above and the second derivative.This is shown below as the second order central second difference. This allows the partial differential equation to be transformed into the finite difference equations. As shown in the figure below the first derivative is represented by the first order forward difference, calculating the time step in the Y direction. The second derivative is used to calculate the spatial difference in the X direction using the second order central difference as this will require information from each of the side points. The figure above also shows a Truncation error also known as the discretization error.
As the domain is separated into a grid of space and time, the calculations only occur at these points so the error represents the difference between the numerical solution and the exact solution. Explicit method The finite difference equation can be used in two different methods to calculate the unknown points in the grid. The explicit one is the simple and straightforward one, it uses the information from the known points in the present time to calculate the value of the point in the future time. It uses the known values of three points such as j1, j2 and j3 at time N and calculates the value of j2 at N+1 in the future.Using initial values as 0 and final value as 1 the methods repeats over and over again to calculate the points at each time step until the final value is reached. (Anderson, j. Etal 2009) Tjn+1= Tjn+? dytdxt Tj+1n-2 Tjn + Tj-1n Implicit method The implicit method on the other hand uses the finite difference equation in a different way.
The explicit equation only has one unknown where as the implicit equation has three unknowns, a matrix is required and the values are solved simultaneously. This allows a larger time step to be used reducing the computation time however this method is more complicated and difficult to program.As shown in the figure below the method uses three known values in the present time step to calculate the three unknowns in the future time step, this takes a longer computation time on each step but due to less steps the overall time is lower than compared to the explicit method. In these methods the error analysis is extremely important. As mentioned earlier the discretization error and the round off error needs to controlled, if this increases from one step to the next then the error will keep on growing and accumulate during the iterations.
To stop this a stability criteria is used which is shown in the equation below, for a defined step size in the special direction requires a step size in time direction to be small enough to satisfy this equation. (LeVeque. R. J,2005). ?dt(dx)2? 12 Numerical Coding The software required to calculate the finite difference equation needs to be powerful such as Fortran, MatLab or C which have the ability to handle large amount of iterations and calculate simulations equations. In this report MatLab is used to create a numerical program, in which the finite difference equation is transformed into a code to create the model.The following explains each part of the code in detail and how the final code is constructed.
Explicit method clc clear Used to clear all information displayed from previous programs Clear is a input required at the end which is used to delete the variables saved in the directory, if the code is changed it ensures that the previous variables do not intervene in the new calculation. n=10 iter=10 n represents the number of grid points used to determine the domain. Iter represents the amount of iterations carried out in this program.
re=5000 y=1/n dt=0. 5*re*dy*dy re is the Reynolds number used to define the flow properties as turbulent. dy is the gird spacing between each point based on the grid size. dt is the time step in the grid based on the stability criteria. for i=1:n u(i)=0 end This is a loop function which determines the initial conditions from 1 to n all values of i will be equal to 0. u(1)=0 u(n+1)=1 This represents the boundary conditions as the bottom plate has a velocity of 0 and the top plate has a velocity of 1.
These will stay constant throughout the program. or i=1:iter for j=2:n unew(j)= u(j)+1/re*dt/dy/dy*(u(j+1)-2*u(j)+u(j-1)) end for k=2:n u(k)=unew(k) end end The outer loop determines the i values between 1 it the number of iterations required, The first inner loop uses the grid points from 2 to n and applies the finite difference equation. The second loop saves the new calculated values as the old ones and uses them as the base point in the next iteration. y=0:dy:1; This determines the y values from 0 to1 in the step size of dy.
plot(u,y) title(‘Explicit 10 grid points’) xlabel(‘Velocity’) label(‘Distance’) legend Iteration-10 clear Finally the graph is plotted of u as the velocity and y as the time step. The title and the labels are defined as names with the legends for each of the curves. clear is a input required at the end which is used to delete the variables saved in the directory, if the code is changed it ensures that the previous variables do not intervene in the new calculation. Final code The code above only produces a graph of one iteration, to see the effects of different iterations the code needs to be modified.The code below calculate 5 different iterations for the same grid points, the same loop as above is carried out 5 times and the values are saved at the end of each loop as a different variable. These variables can then be plotted in the same graphs.
clc clear n=100 aiter=10 biter=20 citer=50 diter=200 eiter=500 fiter=1000 giter=2000 hiter=3000 iiter=5000 re=5000 dy=1/n dt=0. 5*re*dy*dy for i=1:n u(i)=0 end u(1)=0 u(n+1)=1 for i=1:aiter for j=2:n unew(j)= u(j)+1/re*dt/dy/dy*(u(j+1)-2*u(j)+u(j-1)) end for k=2:n u(k)=unew(k) end a=u end for i=1:n u(i)=0 end u(1)=0 u(n+1)=1 or i=1:biter for j=2:n unew(j)= u(j)+1/re*dt/dy/dy*(u(j+1)-2*u(j)+u(j-1)) end for k=2:n u(k)=unew(k) end b=u end for i=1:n u(i)=0 end u(1)=0 u(n+1)=1 for i=1:citer for j=2:n unew(j)= u(j)+1/re*dt/dy/dy*(u(j+1)-2*u(j)+u(j-1)) end for k=2:n u(k)=unew(k) end c=u end for i=1:n u(i)=0 end u(1)=0 u(n+1)=1 for i=1:diter for j=2:n unew(j)= u(j)+1/re*dt/dy/dy*(u(j+1)-2*u(j)+u(j-1)) end for k=2:n u(k)=unew(k) end d=u end for i=1:n u(i)=0 end u(1)=0 u(n+1)=1 for i=1:eiter for j=2:n unew(j)= u(j)+1/re*dt/dy/dy*(u(j+1)-2*u(j)+u(j-1)) end for k=2:n u(k)=unew(k) end e=u end for i=1:n (i)=0 end u(1)=0 u(n+1)=1 for i=1:fiter for j=2:n unew(j)= u(j)+1/re*dt/dy/dy*(u(j+1)-2*u(j)+u(j-1)) end for k=2:n u(k)=unew(k) end f=u end for i=1:n u(i)=0 end u(1)=0 u(n+1)=1 for i=1:giter for j=2:n unew(j)= u(j)+1/re*dt/dy/dy*(u(j+1)-2*u(j)+u(j-1)) end for k=2:n u(k)=unew(k) end g=u end for i=1:n u(i)=0 end u(1)=0 u(n+1)=1 for i=1:hiter for j=2:n unew(j)= u(j)+1/re*dt/dy/dy*(u(j+1)-2*u(j)+u(j-1)) end for k=2:n u(k)=unew(k) end h=u end for i=1:n u(i)=0 end u(1)=0 u(n+1)=1 for i=1:iiter for j=2:n unew(j)= u(j)+1/re*dt/dy/dy*(u(j+1)-2*u(j)+u(j-1)) end for k=2:n u(k)=unew(k) nd i=u end y=0:dy:1; plot(a,y,b,y,c,y,d,y,e,y,f,y,g,y,h,y,i,y) title(‘Explicit 100 grid points’) xlabel(‘Velocity (u)’) ylabel(‘Distance (dy)’) legend Iteration-10 Iteration-20 Iteration-50 Iteration-200 Iteration-500 Iteration-1000 Iteration-2000 Iteration-3000 Iteration-5000 Implicit Method clc clear Used to clear all information displayed from previous programs Clear is a input required at the end which is used to delete the variables saved in the directory, if the code is changed it ensures that the previous variables do not intervene in the new calculation. =10; iter=20; re=5000; dy=1/n; dt=0.
5*re*dy*dy; n represents the number of grid points used to determine the domain. Iter represents the amount of iterations carried out in this program. re is the Reynolds number used to define the flow properties as turbulent. dy is the gird spacing between each point based on the grid size. dt is the time step in the grid based on the stability criteria. for i=1:n; u(i)=0; end u(1)=0; u(n+1)=1; This is a loop function which determines the initial conditions from 1 to n all values of i will be equal to 0.This represents the boundary conditions as the bottom plate has a velocity of 0 and the top plate has a velocity of 1.
These will stay constant throughout the program. A=((1/re)*dt)/(2*(dy*dy)); B=1+((1/re)*dt)/(dy*dy); This represents the two variables used in the matrix for i=1:iter; for j=2:n; unew(j)= -u(j)-(((1/re)*dt)/(2*(dy/dy)))*(u(j+1)-(2*u(j))+u(j-1)); end for k=2:n; u(k)=unew(k) end main_diagonal=-B*ones(n,1) off_diagonal=A*ones(n-1,1) t=diag(main_diagonal)+diag(off_diagonal,1)+diag(off_diagonal,-1) Transpose=t’ squareT=(Transpose*t) T=u(k)*inv(squareT) for i=1:n+1; T=u(i) nd end the outer loop is used to determine the values of the each iteration between 1 and the iter. The first inner loop uses the grid points from 2 to n and applies the finite difference equation in the implicit form, the know values are calculated the present time. The second loop saves the new calculated values as the old ones and uses them as the base point in the next iteration. The A and B variables are used to create a tri diagonal matrix which has the –B value in the middle and A values on each side, as these run down diagonally through the matrix the reset of the values are 0.
An inverse of the matrix is taken which is then multiplied by the u values, this results in the final T values which determine the 3 new points in the next time step. y=0:dy:1; plot(u,y) title(‘Implicit 10 grid points’) xlabel(‘Velocity’) ylabel(‘Distance’) legend Iteration-10 The y values are determined from 0 to1 in the step size of dy and the graph is plotted of u as the velocity and y as the time step. The title and the labels are defined as names with the legends for each of the curves. Final Code clc clear n=10; iter=20; re=5000; dy=1/n; dt=0. 5*re*dy*dy; for i=1:n; u(i)=0; end u(1)=0; (n+1)=1; A=((1/re)*dt)/(2*(dy*dy)); B=1+((1/re)*dt)/(dy*dy); for i=1:iter; for j=2:n; unew(j)= -u(j)-(((1/re)*dt)/(2*(dy/dy)))*(u(j+1)-(2*u(j))+u(j-1)); end for k=2:n; u(k)=unew(k) end main_diagonal=-B*ones(n,1) off_diagonal=A*ones(n-1,1) t=diag(main_diagonal)+diag(off_diagonal,1)+diag(off_diagonal,-1) Transpose=t’ squareT=(Transpose*t) T=u(k)*inv(squareT) for i=1:n+1; T=u(i) end end y=0:dy:1; plot(u,y) title(‘Implicit 10 grid points’) xlabel(‘Velocity’) ylabel(‘Distance’) legend Iteration-10 Numerical results Graph [ 1 ] (Explicit method, 10 grid points with 10 iterations)Graph [ 2 ] (Explicit method, 5 grid points with 10, 20, 30, 40, 50 iterations) Graph [ 3 ] (Explicit method, 10 grid points with 10, 20, 30, 40, 50 iterations) Graph [ 4 ] (Explicit method, 20 grid points with 10, 20, 30, 40, 50 iterations) Graph [ 5 ] (Explicit method, 100 grid points with 10, 20, 30, 40, 50 iteration) Graph [ 6 ] (Explicit method, 100 grid points with 20, 40, 60, 80, 100 iteration) Graph [ 7 ] (Explicit method, 100 grid points with 200, 400, 600, 800, 1000 iteration) Graph [ 8 ] (Explicit method, 100 grid points with 10, 20, 50, 200, 500, 1000, 2000, 3000 iteration)Graph [ 9 ] (Explicit method, 100 grid points with 10, 20, 50, 200, 500, 1000, 2000, 3000, 5000 iteration) Graph [ 10 ] (Implicit method, 10 grid points with 10 iteration) Analysis of the Numerical results The graph 1 which is shown earlier represent the explicit approach of the finite difference method, it uses a grid points in the domain at a spacing of 0. 1. The program carries out 10 iterations at these points to produce the graph.
This graph shows how the program works and the solution is calculated, the initial velocity is 0 at the bottom plate and the final velocity at the top plate is 1.The curve shows the increase in the velocity is slow during the start but becomes faster as it reaches the end, however this curve doesn’t converge to produce the right result therefore the amount of iterations need to be altered. Graph 2 shows the domain with 5 grid points at spacing of 0.
2, and the program is repeated for 5 different iteration starting from 10 up to 50. As shown in the graph the with only 20 iterations the graph converges which means that the solution of the finite difference method reaches the solution of the partial differential equation as the difference between the space and time step go to 0. he produces the right graph besides one issues which is the curvature of the line is not smooth. This means that there are not enough grid points in the domain to gain a accurate and reliable solutions as only certain few areas are being calculated. The smoothness of the curve is improved by increasing the grid points, graph 3 and 4 show the effects as the grid points are increased to 10 and 20.
In graph 3 with 10 grid points the curve finally converges with 50 iterations, and in graph 4 with 20 grid points even with 50 iterations the graph still doesn’t converge.What this shows is the fact that as the amount of grid points increase the curve becomes smoother, but also requires a larger amount of iterations to converge in order to gain the final solution. This results produced is more accurate but also increases the cost and the computation time. Graph 5 and 6 shows the grid points increased to 100, this produces an extremely smooth curve however the amount of iterations used do not produce the converged result.
Even doubling the iterations from 50 to 100 still does not make a major difference.The amount of iterations were then increased as shown in graph 6 and 7 while keeping the grid points the same. From graph 8 it is obvious that the curve nearly converges using 3000 iterations. Final results for the explicit method were represented in Graph 9, it uses 100 grid points which produce a smooth curve but requires 5000 iterations to converge. Using a large amount of grid points allows the grid spacing and time steps to become smaller and closer to 0 producing a more accurate answer with minimum error which is closer to the result of the partial differential equation.
The implicit method did not produce the correct results. As shown in graph 10 the same line was shown with different iterations, if the correct graph was produced it would have the curve would have converged to the correct result with less amount of iterations. However due to the complexity and problems in the implicit code the correct graph could not be produced. Conclusion The importance of computational fluid dynamics can be seen from the research and simulation above. The problem was correctly represented in a program using the FDM numerical technique, and by using an explicit method the result was produced.
The conditions and parameters were changed such as the iterations and the grid size to produce the correct and accurate result. The use of CFD allowed the problem to be easily solved without the need of an actual test in a test lab. The only problem occurred was between the two methods. The implicit method has a faster computing time but as shown in the report the programming of this method is extremely difficult, on the other hand the explicit method took a longer computing time to produce the results but due to its easier and less complicated programming this proved to be the preferred method.References * J. M. McDonough. (2009).
LECTURES IN ELEMENTARY FLUID DYNAMICS. * Jhon D. Anderson, Jr.
(1995). Computational Fluid Dynamics. The Basics with Application. University of Maryland.
* Kuzmin, D. Introduction to Computational Fluid Dynamics. Available: http://www. mathematik. uni-dortmund.
de/~kuzmin/cfdintro/lecture1. pdf. * Zuo, W. (2005). Introduction of Computational Fluid Dynamics. FAU Erlangen-Nurnberg * CHEM 520. Navier-Stokes Equations. Available: http://depts.
washington. du/chemcrs/bulkdisk/chem520A_aut05/notes_Week_05_Lecture_07. pdf.
* Andre Bakker. (2002). Applied Computational Fluid Dynamics. Lecture 5 – Solution Methods * Gilberto E. Urroz. (2004). Convergence, Stability, and Consistency of Finite Difference Schemes in the Solution of Partial Differential Equations. * Anderson, j.
Dick, E. Degrez, G. Grundmann, R. (2009).
In: John F. Wendt Computational Fluid Dynamics: An Introduction. 3rd ed: Technology & Engineering * Randall J.
LeVeque. (2005). Finite Difference Mehtods for Differential Equations.
