CFD Tutorial 13 – Turbulent Flow past NACA0012 Airfoil

Hi! Welcome to QuickerSim CFD Toolbox for MATLAB tutorial number 13. In this tutorial we will simulate a simple flow past a standard NACA0012 airfoil, using a single turbulence model. This tutorial consists of roughly five steps. First of all, we will have to read the mesh and generate the boundary layer, since it is necessary for running proper aerodynamic simulations. Next, we will properly initialize the solution and also set extra parameters for convergence such as under relaxation factors. After that, we will just solve the system as standard. Once the simulations are done, we will visualize pressure coefficient and also compute aerodynamic forces and coefficients. Okay, so let’s go! First of all, standrad thing that is reading the mesh file. The mesh file is provided together with the tutorial, that you can download from our website. The crucial thing as I’ve already said, was the generation of boundary layers. You use the structure that’s called ‘layers’. A few parameters need to be set up: mainly the height of the first element, the ratio of consecutive elements height and the number of elements in the boundary layer. Wall ID numbers also have to be specified. Here is only ID number 9 And finally we need to regenerate our mesh, by using function: extrudeLayers. Then just update match to the second order and visualize the mesh. So let’s execute the script. You can already see, that the mesh is visible. Let’s have a closer look at the boundary layer, how it was created. Yea, let’s look at another section. You can see. Maybe that’s not the finest mesh, but should be fair enough for performing a simple simulation. Okay, let’s see how it looks further on stream. Yea, you can see a few layers of the boundary layer present. Each one is thicker than the previous one. Okay, let’s continue. Let’s just set fluid viscosity. It will be air, so it’s like 1.5e to minus 5. Okay let’s initialize our solution. We do it in a standard way. Use the function initSolution, use zero velocities and pressures as the starter. We will also create the copy of the velocity pressure field, because it will be used in the under relaxation process Maximum revel of a level of residuals will be 10 to minus 6 and maximum number of iterations, let’s say, should be a hundred. We will also need the relaxation factor and let’s get it try 0.6. Another necessary thing, in the turbulent flow simulations, is the computation of wall distance. We will store the distances of cell centers to the airfoils walls in the variable ‘d’. We will get them using inbuilt function: computeWallDistance. WallIds are number nine. Now you’ll just need to define the value on direction of the inlet velocity. That’s pretty simple. 0.75 meters per second will give us Reynolds numer equal 50 thousands. We specified the angle of attack to be five degrees to recalculated into radians and also get proper components of inlet velocity. Okay we just check that everything is fine, just executing our script. So far so good. And now we set up the main interative loop. I will give it as is and then discuss the details. Okay, so iterate from 1 to max number of iterations. A few things should be already pretty well known to you. We assemble Navier-Stokes matrix, get NS and F. We impose boundary conditions at the inlet at the wall. Compute the residuals and plot the residuals. The crucial difference in here is that we include the turbulence model. The turbulence model actually only effects on turbulent viscosity that we call the parameter ‘nut’ and it’s an addition to standard viscosity in the standard Navier-Stockes formulation. We will be using CITM turbulence model Basically, it’s the hybrid of two models. If you want to learn more details about it just check out our manual turbulentViscosity. We will notice that CITM is constant intensity, just hybrid of Van Driest and Velocity Length Model. Another comment is in place. I mentioned, that will be using relaxation factor, that means that before we assemble Navier-Stokes matrix in the next iteration step, we combine the previous solution with the yet one step back that was safe at ‘us’ before we will solving in the system. That will help increase that convergence quality. Okay. Now, let’s solve the system We can observe how it starts converging. Great, the calculations are finished. The solution has converged. Let’s move to visualizing the solutions. Okay actually, before we have a look at the pressure coefficient. Let’s have a look at velocity. Just establish it another figure and execute the piece of code. That’s the solution displaySolution2D ‘x-velocity’ Let’s look how it looks close to the profile. There’s some sort of weight behind the profile. Let’s have a closer look just at the leading edge. Yea, so we can see stagnation point a little bit below the main line, because obviously it’s a positive angle of attack. And you can also watch how the boundary layer develops along the profile. Yea, so it starts growing as it should be. and there it goes. Okay, now we will move to calculating the pressure coefficient in the whole domain. First of all, since pressure is stored only on the first order nodes. We need to interpolate pressure data to neighboring nodes using pressure=generatePressureData, the domain. There we go to the new field. Furthermore we will need to find the maximum value of pressure which is obviously the stagnation pressure that appears somewhere on the airfoil. In order not to check for the maximum in the flow domain, we will just extract the indices of the nodes that lay on the airfoil and look for the maximum over there. We use function: extractNodeIdsOnEdges. We get the node edge wall number and then we look for the maximum over there, somewhere with this formula and we can see that there’s a value of pmax, which is larger than average zero pressure. Okay, furthermore we need to define the coefficient. what we do is simply calculate the ratio of local pressure minus the maximum pressure plus a half of the infinity squared divided by half of infinity squared. Please remember, that density is not included in our pressure therefore, it doesn’t appear here, nor here. Finally, we can set up the plot, so we need to establish yet another figure and display the whole solution. Let’s go! Here comes the plot. Let’s have a closer look on that. And exactly pressure coefficient equal to one is in the stagnation point, just below the leading edge of the airfoil. And the lowest pressure occurs just on the upper part of the airfoil in the dead zone. Lasted for today, is to calculate the forces and force coefficient that are on the profile. There is a built-in function in the Toolbox. That’s called computeForce, that would integrate the loads on specified wall. Here we specified by number nine. The resulting force obviously is giving XY reference frame horizontal and vertical. But we want to get the drag and lift. Therefore, we need to rotate this result in the direction that’s oripted with the flow. We use just rotation matrix of sines and cosines. We can execute the code to have a look at our results. Obviously, the sheer various of forces are little use when we want to compare different profiles We need to get the non dimensional coefficient simply by dividing lift and drag by half of velocity square. On the last step for today is to evaluate ‘cd’ and ‘cl’. Especially ‘cl is close the experimental value, while ‘cd’ is a little bit more complex. Therefore it may a little bit depart from the data that is given in the literature Anyway, that’s all that we wanted to present you in this tutorial. We hope, that you enjoy!

2 thoughts on “CFD Tutorial 13 – Turbulent Flow past NACA0012 Airfoil

Leave a Reply

Your email address will not be published. Required fields are marked *