This tutorial describes how to set up and run the first non-trivial flow problem. The lid-driven cavity flow is a standard test case for numerical schemes, and a number of results have been published in literature, see e. This tutorial assumes that you have completed the previous tutorial, know how to edit files and postprocess the solution with your favorite visualization tool, e.
Also, the later parts of the tutorial assume that you have access to a computer with an MPI-based parallelization with at least 4 computing cores — otherwise, it will just take a lot longer :. The tutorial is split into two parts: The basic part will teach you about setting up the code and running the simulations. The advanced part will build on this and give you a glimpse on how to make modification to code to accommodate more complex simulation and add features you might need.
If you are just interested in running the code as is, you may skip the advanced part, or just complete parts of it.
The standard settings are sufficient for this example. To check whether they are set, change to your build folder and open the cmake GUI. These settings are identical for both parts of the tutorial. Now continue with the basic or advanced tutorial! Advanced Tutorial. Skip to content. Lid-driven Cavity For the pdf version of this tutorial, see the userguide. The flow under consideration is essentially incompressible and two-dimensional, but we will use the three-dimensional code for the compressible Navier-Stokes equations to solve it here.
This is not the most efficient way to compute this flow, but it works well to show you how to set up and run a simulation in this tutorial. The wall of the cavity are modeled as isothermal walls, and a fixed flow is prescribed at the upper boundary, i. For the Reynolds numbers investigated here, this generates a steady, vortical flow field in the cavity.
Ghia U. Gao Z.Updated 10 Sep Cavity flow is simulated using the pressure correction method on a staggered grid using explicit differencing for the hyperbolic terms CD, MacCormack and Richtmyer method while both explicit and implicit methods are considered for the diffusive parabolic terms.Tutorial 2 Lid Driven Cavity IITP
For the implicit steps, preconditioned matrices are used using LU decomposition. The Pressure Poisson equation is also solved implicitly. Neumann boundary conditions are used for pressure and Dirichlet conditions for the velocity field. Suraj Shankar Retrieved April 15, The result seems incorrect, due to a strange recirculation in the upper right corner.
The code is meant to be pedagogical in nature and has been made in line with the steps to Navier-Stokes practical module, for which I would like to credit Lorena Barba and her online course on CFD. Learn About Live Editor. Choose a web site to get translated content where available and see local events and offers.
Based on your location, we recommend that you select:. Select the China site in Chinese or English for best site performance. Other MathWorks country sites are not optimized for visits from your location. Toggle Main Navigation. File Exchange.
Search MathWorks. Open Mobile Search. Trial software. You are now following this Submission You will see updates in your activity feed You may receive emails, depending on your notification preferences. Follow Download. Overview Functions. Cite As Suraj Shankar Comments and Ratings Luiz Miranda Luiz Miranda view profile. Radu Trimbitas Radu Trimbitas view profile. Sagar Sagar view profile. Shubham Maurya Shubham Maurya view profile. Seems to be caused by the projection step. Do you know how to fix that?
Hung Nguyen Hung Nguyen view profile.Engineers and Mathematicians -- Get Lost! This is a non-rigorous webpage intended to introduce principles to medical folks. You're welcome to images, videos etc. I'd also be happy to hear from you if I've made an obvious goof. The video above depicts an animation of the fluid flow in a Lid Driven Cavity that results when the "lid" or upper side of the square cavity is drawn to the right at a constant velocity, thereby inducing the fluid into a circular motion, i.
This doesn't have much to do with cardiology, but is here to impart some principles of fluid dynamics for a flow situation that's a little more complicated than what you may be used to thinking about.
The lid-driven cavity LDC is a common test or bench-mark problem in computational fluid dynamics CFD particularly as one that critically tests the accuracy of the advection convective acceleration scheme used for the computations. The figure that follows is an example calculation that illustrates the essential features of the computation which consists of a rectangular cavity, a square one in this case, where the Newtonian fluid inside is set in motion by dragging the upper edge of the cavity the "lid" to the right.
For a real viscous fluid, the fluid elements in contact with the lid move at the same velocity as the lid. The right, left, and lower boundaries of the cavity are stationary so the fluid at those borders has a velocity of zero. The test problem is invariably 2-dimensional as implied by the figure; it's as if we're taking a tomographic slice of a rectangular trough that extends infinitely into and out of the plane of the screen paper?
This particular solution was computed for a Reynolds number of ; this number constitutes the overall characteristic relationship between the inertial convective acceleration and viscous forces governing the flow.
In the figure as shown, velocity of the fluid is represented by arrows that indicate both magnitude and direction of velocity associate with each cell; color also represents velocity magnitude. The moving lid sets the fluid into a grossly circular motion, a vortex. If you look closely at the above figure you will see that the computational grid is also included in a light blue color.
The particular solution above was determined for a x cell, 2-dimensional grid. FYI, that means that the solution comes down to a set of equations with 10, variables!
Note that only a fraction of the computed velocities are shown so that the figure isn't saturated with vectors. The next figure shows more of the vectors still not all and you can see how it starts to become more difficult to interpret.
Note: you can see a full-size rendition of any of the figures in this article by right clicking the figure and selecting "Open image in new tab". The next 2 figures also show the velocity vectors, still colored by velocity magnitude, but in a format that might be easier to interpret.
Once again not all of the vectors are being shown, but this format has its advantages for some purposes. ALL of the vectors are shown that originate at a number of selected slices through the cavity. However this format tends to obscure the fact that there is a vortex whirling around a particular location inside the cavity. The next 2 figures are contour plots — essentially each is a "topographic map" showing you where the highs and lows are.
The figure on the left shows velocity magnitude contours. Note however that while the actual velocity in the solution ranges from 0. This was done because many of the velocities in the solution are lower than 0. Similarly pressure contours in the right-hand figure have also been adjusted.
Obviously there is a big pressure high point in the upper right-hand corner of the cavity and a low point in the upper left-hand corner.
There isn't a lot of variation in the pressure over the rest of the solution, but adjusting the color scale at least allows us to see that there is a relative pressure low at the center of the vortex. A fundamental aspect of fluid mechanics is that the characteristics of a flow are governed by the Reynolds number Re. The following sequence of figures is one of several locations of the website that illustrates and examines the fact.
The first group of 6 figures show solutions for increasing values of Re. The figures are essentially polluted with vectors, but this particular format allows us to readily see the location of the vortex center the "eye of the hurricane"which moves from a location somewhat towards the upper right quadrant towards the center of the cavity with increasing Re.Initially, the flow will be assumed laminar and will be solved on a uniform mesh using the icoFoam solver for laminar, isothermal, incompressible flow.
During the course of the tutorial, the effect of increased mesh resolution and mesh grading towards the walls will be investigated. Finally, the flow Reynolds number will be increased and the pisoFoam solver will be used for turbulentisothermal, incompressible flow. Here, the mesh must be 1 cell layer thick, and the empty patches planar. The cavity domain consists of a square of side length in the - plane.
Sample code for SIMPLE-algorithm and solving Lid-Driven Cavity flow test is uploaded
A uniform mesh of 20 by 20 cells will be used initially. The blockMesh mesh generator supplied with OpenFOAM generates meshes from a description specified in an input dictionary, blockMeshDict located in the system directory for a given case. For the remainder of the manual: For the sake of clarity and to save space, file headers, including the banner and FoamFile sub-dictionary, will be removed from verbatim quoting of case files The file first specifies the list of coordinates representing the block vertices ; These are in arbitrary units, and can be scaled to the real problem dimensions using the scale entry, e.
The final section defines the boundary patches. The mesh is generated by running blockMesh on this blockMeshDict file.
Any mistakes in the blockMeshDict file are picked up by blockMesh and the resulting error message directs the user to the line in the file where the problem occurred. There should be no error messages at this stage. The 0 sub-directory contains 2 files, p and Uone for each of the pressure and velocity fields whose initial values and boundary conditions must be set. For this case cavitythe boundary consists of walls only, split into 2 patches named: 1 fixedWalls for the fixed sides and base of the cavity; 2 movingWall for the moving top of the cavity.
The frontAndBack patch represents the front and back planes of the 2D case and therefore must be set as empty. In this case, as in most we encounter, the initial fields are set to be uniform. Here the pressure is kinematic, and as an incompressible case, its absolute value is not relevant, so is set to uniform 0 for convenience.
The dimensions are those expected for velocity, the internal field is initialised as uniform zero, which in the case of velocity must be expressed by 3 vector components, i. The boundary field for velocity requires the same boundary condition for the frontAndBack patch. The other patches are walls: a no-slip condition is assumed on the fixedWallshence a fixedValue condition with a value of uniform 0 0 0. For an icoFoam case, the only property that must be specified is the kinematic viscosity which is stored from the transportProperties dictionary.
The keyword for kinematic viscosity is nuthe phonetic label for the Greek symbol by which it is represented in equations.
Initially this case will be run with a Reynolds number of 10, where the Reynolds number is defined as: 2. The user should view this file; as a case control file, it is located in the system directory.
Therefore we set the startFrom keyword to startTime and then specify the startTime keyword to be 0. For the end time, we wish to reach the steady state solution where the flow is circulating around the cavity.
As a general rule, the fluid should pass through the domain 10 times to reach steady state in laminar flow. In this case the flow does not pass through this domain as there is no inlet or outlet, so instead the end time can be set to the time taken for the lid to travel ten times across the cavity, i. To specify this end time, we must specify the stopAt keyword as endTime and then set the endTime keyword to 0. Now we need to set the time step, represented by the keyword deltaT.
To achieve temporal accuracy and numerical stability when running icoFoama Courant number of less than 1 is required.The lid driven cavity is a classical benchmark for Navier Stokes solvers.
This is a demonstration of how the Python module shenfun can be used to solve the lid driven cavity problem with full spectral accuracy using a mixed coupled basis in a 2D tensor product domain. The demo also shows how to use mixed tensor product spaces for vector valued equations. Note that the regular lid driven cavity, where the top wall has constant velocity and the remaining three walls are stationary, has a singularity at the two upper corners, where the velocity is discontinuous.
We want to solve these steady nonlinear Navier Stokes equations with the Galerkin method, using the shenfun Python package.
Note that MPI for Python mpi4py is a requirement for shenfun, but the current solver cannot be used with more than one processor. With the Galerkin method we need basis functions for both velocity and pressure, as well as for the nonlinear right hand side. A Dirichlet basis will be used for velocity, whereas there is no boundary restriction on the pressure basis. And then we create two-dimensional basis functions like.
Special attention is required by the moving lid. In general, a nonzero boundary condition can be added on both sides of the domain using the following basis. To implement this boundary condition instead, we can make use of sympy and quite straight forward do.
Uncomment the last line to run the regularized boundary conditions. Otherwise, there is no difference at all between the regular and the regularized lid driven cavity implementations. The pressure basis that comes with no restrictions for the boundary is a little trickier. The reason for this has to do with inf-sup stability.
We have now created all relevant function spaces for the problem at hand. It remains to combine these spaces into tensor product spaces, and to combine tensor product spaces into mixed coupled tensor product spaces.
Tutorial 2 – Lid Driven Cavity
From the Dirichlet bases we create two different tensor product spaces, whereas one is enough for the pressure. These tensor product spaces are all scalar valued.
The mixed basis is created in shenfun as. The first obvious issue with Eq 11 is the nonlinearity. In other words we will need to linearize and iterate to be able to solve these equations with the Galerkin method.
All boundary integrals disappear since we are using test functions with homogeneous boundary conditions. Note that the bilinear form will assemble to a block matrix, whereas the right hand side linear form will assemble to a block vector. The bilinear form does not change with the solution and as such it does not need to be reassembled inside an iteration loop.
We will now implement the coupled variational problem described in previous sections. The test function v is using homogeneous Dirichlet boundary conditions even though it is derived from VQwhich contains W1. It is currently not and will probably never be possible to use test functions with inhomogeneous boundary conditions. With the basisfunctions in place we may assemble the different blocks of the final coefficient matrix. For this we also need to specify the kinematic viscosity, which is given here in terms of the Reynolds number:.
This is because the inner function returns a list of tensor product matrices of type TPMatrixand you cannot negate a list.
The assembled subsystems A, G and D are lists containg the different blocks of the complete, coupled, coefficient matrix. A actually contains 4 tensor product matrices of type TPMatrix. The first two matrices are for vector component zero of the test function v and trial function uthe matrices 2 and 3 are for components 1.
The first two matrices are as such for.Dear colleagues! It has a lot of shortcomings, but I think that it will be useful for all newcomers. I intend to develop it further, but now it very slowly converges. All notes will be useful, and hope it can be discussed. May 18,Simple-Algorithm.
Mijail Febres. It's not working properly, velocities converge in 2 or 3 iterations but pressures don't, I think I have a problem with BC's treatment constant known velocity on one side and constant known pressure on the other. Do you know any book that shows a clear example on boundary conditions treatment? The code you wrote is on a colocated grid, do you think it has advantages over a staggered grid for the 1-d problem I'm trying to solve?
Thanks in advance. May 18, Hi all. The code worked really well for Couette flow.
Sample code for solving Lid-Driven cavity test (Re=1000) - Fortran 90
However for lid driven cavity, it is not showing required results. I have tried using different Reynolds number like 10, Initially the code seems as if it will give the required result, however after a few iterations, oscillations begin and the solution fails.
Is there any stability criteria for this flow problem? The time step I had used was 0. Any kind of help in this regards will be highly appreciated. Thanks a lot. October 12, Join Date: Oct Hi Hamza We specify both velocity components on the wall u and v 0,0 for solid walls, 1,0 for moving wall Just fixed them they are not corrected during solution. December 27, Hi, thank you Michail for your answers they are very useful, i realized the code for Newtonian and Non-Newtionian fluids.
Now i am trayin the viscoelastic case, but i have a problem of numerical stability when increasing WE WE is a parameter in the equationsi want to know if you have an idea about this. December 28, First of all I would like to ask - how many equations do You solve? If 3 You are wrong. There are 3 additional equations for stresses.
They are solving the same way as momentum equations. But not too complicated. December 29, Hi michail, yes i see, we have three unknowns added so automatically we need 3 others equations which are given in the constitutive equations. Hallo the source term in equations of motion for U and V summarized from pressure gradient and stresses. Check for bugs.The contents of the folder should look like this picture.
In case you do not want to generate the mesh files yourselves, the meshes have already been provided. The next line maps the 6 faces of the cube to the boundary conditions following the CGNS notation, i. Again, the corresponding states can be set here, but overwritten later. For this basic tutorial, the simple meshes shown in this picture will be used. To get help on any of the parameters listed therein, you can run FLEXI from the command line by typing. The parameters in the file are grouped thematically, however, this is not mandatory.
NAnalyze determines the polynomial degree for the analysis routines, e. The boundary conditions are set via the BoundaryName identifier, which specifies a name that must be present in the mesh file see section on mesh generation above. Each line containing the boundary name must be followed by a line containing the BoundaryType. A list of types available for the Navier-Stokes equations can be found in the chapter of the userguide.
Note that the reference vectors are always in primitive variables, i. These boundaries are isothermal walls, so a wall temperature needs to be specified — this is computed from the primitive variables in the associated reference state. Note that the lines for the lower wall boundary are commented out. As has been mentioned above, the boundary conditions can be set in HOPR directly, and then be overwritten here later. This is just an example to show how both approaches work.
If you wish to run the code in parallel using MPI, the standard command is. One way to check if the solution has converged to a steady state is to check some characteristic quantities, e. Note that since the boundary conditions are applied weakly in a DG setting, a velocity slip at walls can occur, with its magnitude depending on the local wall resolution. For comparison with published data, 4 simulations where run on the provided meshes.
Details are listed in the table below. A contour plot of the velocity magnitude at the ent time is given in this figure. To generate this plot, convert the State files to a Paraview format by invoking. This figure shows the results for the 4 simulations run here, along with published data available in [ 1 ][ 2 ]. Skip to content. Also, we will only use one element in that direction to save computational costs. The identifier BoundaryType specifies what kind of boundary is to be applied to the face.
The second index indicates if the face is curved not for this test caseand the third indicates which reference state will be used at the boundary. The fourth entry is used to match periodic boundary condition sides, i. The parameter ProjectName is used to name all the files generated by the simulation, e. Turning it on will result in an on-the-fly output of the solution states, but this is not recommended as good practice — it will slow down the code considerably, in particular in parallel mode.
The recommended way is to use the flexi2vtk tool after the simulation to generate Paraview-readable files. The parameter N sets the degree of the solution polynomial, e. UseCurveds indicates whether the mesh is considered to be curved, i.