Transcript for:
Writing Code for CFD (Computational Fluid Dynamics)

alright so the next part in this video series is about how you can write a code for CFD so so far we have done some examples in ences so that was as part of the software so what happens when you want to write a code by yourself so you want to try a new scheme or you want to try something else so that's where coding comes into the picture and that is the need for the coding so you can write your own code starting from the scratch or starting from some other codes and then you can use them for the CFD so if you remember that the first example that we did was for the lid driven cavity and that's going to be the same problem that we are going to discuss here to write the code so this was the very problem statement that I have the square domain and my walls that is the left wall the right wall and the bottom wall they were completely stationary so on those walls in the u velocity and the V velocity that is the horizontal and the vertical velocity were zero while I said that the top wall is moving so because it was moving in the X direction so only the u component of the velocity is some value say u capital u while the V velocity that is the vertical velocity would be zero so we solve that problem in ANSYS we got some good results but we didn't compare it with the benchmark and we just say that we should always do that because that's the essence of the CFD it's not just some results it's the correct results that you should get so when we want to solve this kind of problem we discretize the problem domain or we say that we want to make the small parts of this domain or we want to create a mesh a mesh on which we will solve the equations so to create the mesh I'm going to be using the very simple Cartesian and very structured mesh over here as you can see on the right hand side so what I did is I just divided this whole square into a number of smaller squares and I have some points on the boundary which are in the red and I have some points inside the square which are in the green and I will call them as interior nodes so right now I have divided my one length into the six intervals so the number of points that I have in each direction is seven so I can say that I have a mesh of a size seven by seven this is only for the representation so as you might be aware of that we don't exactly solve the equations directly because of the so called checker board oscillations so what we do to avoid that kind of unphysical oscillation is that we do staggering of the grids struggling is that we displace the velocity components rather than being on the nodal points so if I say that I want to have the grid for the u velocity so this was my original mesh so if I want to stagger my u velocity that would mean that instead of defining my u velocity at the nodal point I would define it at the half distance from the node so if I say that I want to define my u velocity for this corner point at this place that is half the distance it could be in this case it's in the Y direction so you can see that it's only displaced in the Y direction so it could be displaced in the X direction also so you can try that but for now we are saying that we are displacing every X every X velocity or every u velocity in the Y direction so similarly for the other point it could be here for the third point it would be over its head and similarly for all these points on the boundary you can do it like this and for every other point you can do it like this but I'll also add one more additional row over here the pot which I will be explaining later so this would be my total grid for the u velocity and if you will count it it would be a seven by eight mesh because I have this additional u velocity layer so that would make it 7 by 8 and now the second part would be the grid for the V velocity so this was my original mesh and because I displace my u velocity in the Y direction and so I have to move my V velocity in the X direction so for this particular point I would be displacing it in the X so along this side so there comes my V velocity so similarly for this point I have my B velocity here and for the other points I could have my Green velocities like this so for every other point I can do it in the same way and I will have my V velocity distribution like this and similar to the u velocity I will be adding one layer over here and I will explain it to you later that why am i adding this layer and if you count it out it would be a eight by seven layer because of this extra you extra V velocity said that I added so this is right now the grids the different grids for the u velocity and the V velocity and it looks pretty complicated already so the third variable that we want to deal with is the pressure so for the pressure I will say that I have my pressure point something like this in the yellow and I will explain that why I had my pressure points over here and you can see that the pressure is 8 by 8 so my real mesh was 7 by 7 my u velocity turned out to be 7 by 8 and V velocity to be 8 by 7 and the pressure is actually 8 by 8 so if I put these all things together it would be very horrible and it looks something like this just don't get scared already our will be dealing it with in very nicely and very slowly so the reason why i define pressure over here is so that the control volume for the pressure would be represented at something like this that the u velocity comes in the u velocity goes out the V velocity comes in and the V velocity goes out so the only reason why I want these points at this particular position is such that the control volume can have their meanings so for the u velocity because the U velocity should only depend on the pressure gradient across the X direction so if I have these particular pressure points like this so my you velocity it is affected by this pressure gradient and similarly for the V velocity that V velocity is affected by this pressure gradient and that is what the checkerboard problem actually was that instead of these neighboring points it was affected by the alternate points so that was creating some non physical solutions that was never possible but mathematically they were true so with this grid set in mind let's deal with the control volumes one by one so let's talk about the X momentum equation so for the X momentum I want the X u velocity control volume so this is the u velocity in question and I have the control volume here and the pressure is over here and some of the U velocity velocities like this so I have numbered them 1 2 3 4 the pressure point as small a and small W that is east and west and for the u velocity I said that this is QP or just the P and this is east west north and south so you can cross-check it over here in the U velocity control volume that this is the U velocity P these four would be the surrounding east west north and south and these 1 2 3 4 and the pressure is 10 west so this is my momentum equation over here and I have added this particular term so that I can iterate in time and at the steady state this term would vanish so I will retain the original steady state jus momentum equation so there are some terms that I would be talking about so the first term is del u over del T so this is basically the time derivative of the U velocity so I cannot using a forward difference I can write it as the new velocity minus the old velocity or the current velocity divided by DT and your DT will depend if you're using an explicit scheme so you will have to take care of the CFM number and the second term is the U square over del X so if I walk if I use a central different scheme because this is del over del X of U square so if because I already know the velocities at different points in the X direction so I can use a central difference for this first derivative and I can simply write it as u e square minus UW square over 2 DX and this should be second-order accurate the next one is Del of del over del Y of UV that is the derivative first derivative in the Y direction of the multiplication of UV so just consider this control volume so if I want to have a derivative in the Y direction so this is the Y direction and this is the t-virus tense so if I consider this particular point that is center point over here then the multiplication of UV can be interpolated as first U is the yo u P plus u n divided by 2 and V is V 3 plus V 4 divided by 2 so this is the multiplication of UV the interpolated multiplication at this point and similarly you get two at this point which would be u P plus us over 2 and V 1 plus V 2 over 2 at this particular point and divided the whole thing divided by the distance that is dy similarly for Delta P over Delta X it shows the pressure gradient along the x direction which is PE minus PW over TX the last one is the viscous term which is the new and the laplacian of u velocity so this is basically the Delta square u over Delta X square plus Delta square u over Delta Y square and I would be using a second order central difference and with that I can write this particular terms in terms of the U at all these five points so that complete my description of how you can solve for the X momentum equation and the same thing you can do for the y momentum equation because these equations they are very similar the only thing is the convection and the viscous terms they they are in terms of V rather than Q so the next thing is the continuity equation so originally the continuity equation is only the divergence of velocity is zero so I added this particular term and when I write my equation in this form this is called as the artificial compressibility method I introduce an artificial term and originally it's in the form of Rho so that is it is like the compressible flow equations and that is why it's called as the artificial compressibility equation so because if you want to use this particular scheme it might not be available in your commercial software's and that's why we might need to write a code if you want to use this particular scheme and because this game is so damn simple so I love this scheme and moving on to the continuity equation so for the continuity I have my pressure control volume I have the U velocity coming in going out V velocity coming in and going out so for this particular term I have my the pressure and time derivative to be just like the last last page for the U velocity derivative that is P nu minus P o divided by DT and the delta is just a constant parameter that you can set and the Delta U over Delta X is simply u minus UW over DX v Delta V over Delta Y is V n minus V s over dy so all these terms we know it from either you can use the velocities that are at the current time step or you can also use the velocities that you just solve using the X&Y momentum equations so it would be faster if you use the updated velocities rather than the old velocities and using this equation you can get the pressure so one thing that this equation does is it gets the pressure directly so you don't need to use any pressure correction kind of schemes or a simple method or any other method so this equation directly gives you the pressure and that's one of the major use of this particular kind of equation that you get it so easily and that and the one thing about the incompressible flow equation is the decoupling between the velocities and the pressure and this particular scheme eliminates that pretty nicely so you get the new value you get the you get the new values for U and V and you keep doing that until so there is some criteria so usually the criteria is until the divergence of the velocity is zero so next thing and the very important thing is the boundary condition because we have the boundary so we have to have the boundary conditions so further you velocity I say that I can solve for the interior points but if I want to solve for the points at the boundary I may not have the all the values that are needed for from that are needed for the solution of the X momentum equation but because we know that on this well on this left and the right wall the you should be zero so on the left and the right wall I can simply say that these u velocities they should be zero the velocities that are in the red the velocity arrows that are in the red they should be zero so the first thing that we can do is to put them zero but we also have this layer in the bottom the yellow layer that we added extra so this is actually for the satisfaction of the boundary condition because we we will solve the X momentum equation for this particular u velocity but if if we don't have this velocity because in the physical situation the red node over here it should be zero so if if we want to make it zero we need to have something here so that these two things together they can make this zero so this happen in such a way that for this to be zero their average should be zero or in other words say if this velocity is one and if I say that if this velocity could be minus one then for this particular point it would be zero because we wanted zero the physical conditions say that it should be zero so that is why this yellow layer it's very important it it compensate for the change in the internal velocities and the same thing will happen at the top because at the top we want the U equal to U and because we can solve for this particular the blue u velocities but the physical conditions say that the red points at the top they should be one so that's why the average of the top and the next second top it should be u so that is why there is a layer at the top which should be averaged with the next layer and it should be equal to u so this Plus this divided by 2 should be equal to u so that the real lid or the real layer at the top should we have a velocity of U and the same goes for the boundary at the V velocities so again because the top and the bottom has the V equal to 0 so we can directly assign these V velocities because they are situated at the boundaries but the boundary belt but the velocities that are not at the boundaries we can do the same thing that we did over here that the average should be 0 or these yellow there would be negative of the brown ones so that their average is 0 and the physical condition is satisfied and the next part is the pressure boundary condition because you also have the pressure so you can solve for the interior pressure and for the exterior pressure it would be simply that exterior pressure across the interface or across the wall should be equal to the interior one because because you know that you cross this interface the V equal to CR the u equal to 0 so there is no transfer there is no transport across the wall so because the u equal to 0 there is no you and that suggests the pressure difference is should be zero because if there is a pressure difference there would be a flow so if the pressure difference is zero there would be no flow and same goes for here the V velocity is zero so there should be no flow across the interface these two should be equal these two should be equal and these two should be equal so that is all about how you can write the code so this was about to give you an idea about how you can write the code and in the next video I will show you the code that I have written myself using all the things that I taught you in this particular video so just hang on to it and click on the next one