Category Archives: CFD

Rayleigh–Bénard Convection Using buoyantBoussinesqPimpleFoam

Here is an extremely simple simulation to set up that has a surprisingly beautiful output. In this post, we will simulation the classic Rayleigh–Bénard convection (see Wikipedia) in 3D using the buoyant solver, buoyantBoussinesqPimpleFoam.

buoyantBoussinesqPimpleFoam is a solver for buoyant, turbulent, incompressible flows. The incompressible part of the solver comes from the fact that it uses the Boussinesq approximation for buoyancy which is only suitable for flows with small density changes (aka incompressible). A good source for the equations in this solver is this blog post.

Simulation Set-up

The basic set-up for this case is simple: a hot bottom surface, a cold top surface, either cyclic or zero-gradient sides, and some initial conditions.

For this example, I used the properties of water, I set the bottom plate at 373 K (don’t ask me why… I know it’s close to boiling point of water), and the top plate at 273 K. For this case, we will not use any turbulent modeling and will simply use a laminar model (this simply means there is only molecular viscosity, there are no simplifications applied to the equations).

Geometry and Mesh

The geometry is shown below. As the geometry is so simple… I will not go over the blockMesh set up. The mesh discretization that I used was simplegrading (i.e. no inflation), with 200x200x50 cells.

GeometryLabelled

Constant

For this case, we will simulate water. The transportProperties file should look like:

transportModel Newtonian;

// Laminar viscosity
nu [0 2 -1 0 0 0 0] 1e-06;

// Thermal expansion coefficient
beta [0 0 0 -1 0 0 0] 0.000214;

// Reference temperature
TRef [0 0 0 1 0 0 0] 300;

// Laminar Prandtl number
Pr [0 0 0 0 0 0 0] 7.56;

// Turbulent Prandtl number (not used)
Prt [0 0 0 0 0 0 0] 0.7;

I don’t typically delete unused entries from dictionaries. This makes using previous simulations as templates much easier. Therefore note that the turbulent Prandtl number is in the dictionary but it is not used.

Selecting Reference Temperature TRef for buoyantBoussinesqPimpleFoam. To answer this recall that when the Boussinesq buoyancy approximation is used, the solver does not solve for the density. It solves the relative density using the linear function:

\frac{\rho}{\rho_0}=1-\beta \left(T-T_0\right)

Therefore, I think it makes sense that we should choose a temperature for T_{ref} that is somewhere in the range of the simulation. Thus I chose Tref=300. Somewhere in the middle!

And the turbulenceProperties file is:

simulationType laminar;

RAS
{
 RASModel laminar;

turbulence off;

printCoeffs off;
}

The g file tells the solver the acceleration due to gravity, as well as the direction:

dimensions [0 1 -2 0 0 0 0];
value (0 -9.81 0);

Boundary Conditions

In the “zero” folder, we need the following files: p, p_rgh, T, U, and alphat (this file needs to be present… however it is not used given the laminar simulationType.

T:

dimensions [0 0 0 1 0 0 0];

internalField uniform 273;

boundaryField
{
 floor
 {
 type fixedValue;
 value uniform 373;
 }
 ceiling
 {
 type fixedValue;
 value uniform 273;
 }
 fixedWalls
 {
 type zeroGradient;
 }
}

p_rgh:

dimensions [0 2 -2 0 0 0 0];

internalField uniform 0;

boundaryField
{
floor
{
type fixedFluxPressure;
rho rhok;
value uniform 0;
}

ceiling
{
type fixedFluxPressure;
rho rhok;
value uniform 0;
}

fixedWalls
{
type fixedFluxPressure;
rho rhok;
value uniform 0;
}
}

p:

dimensions [0 2 -2 0 0 0 0];

internalField uniform 0;

boundaryField
{
 floor
 {
 type calculated;
 value $internalField;
 }

ceiling
 {
 type calculated;
 value $internalField;
 }

fixedWalls
 {
 type calculated;
 value $internalField;
 }
}

U:

dimensions [0 1 -1 0 0 0 0];

internalField uniform (0 0 0);

boundaryField
{
 floor
 {
 type noSlip;
 }

ceiling
 {
 type noSlip;
 }

fixedWalls
 {
 type noSlip;
 }
}

Simulation Results

The results (as always) are the best part. But especially for this case since they are so nice to look at! I have made a couple animations of temperature fields and contours. Enjoy.

animation
3D Temperature Contours

 

cutanimation
Temperature Field – Slice Through xy Plane

 

Conclusion

This case demonstrated the simple set up of a case using buoyantBoussinesqPimpleFoam. The case (Rayleigh-Bénard convection) was simulated in 3D on a fine grid.

Comments and questions are welcome! Keep in mind I set this case up very quickly.

 

Time-Varying Cylinder Motion in Cross-flow: timeVaryingFixedUniformValue

This post is a simple demonstration of the timeVaryingFixedUniformValue boundary condition. This boundary condition allows a Dirichlet-type boundary condition to be varied in time. To demonstrate, we will modify the oscillating cylinder case.

Set-Up

Instead of using the oscillating boundary condition for point displacement. We will have the cylinder do two things:

  • Move in a circular motion
  • Move in a sinusoidal decay motion

The basics of this boundary condition are extremely simple. Keep in mind that although (here) we are modifying the pointDisplacement boundary condition for the cylinder, the basics of this BC would be the same if you were doing a time varying boundary condition for say pressure or velocity.

In the pointDisplacement file:

 cylinder
 {
 type timeVaryingUniformFixedValue;
 fileName "prescribedMotion";
 outOfBounds clamp;
 }

fileName points to the file where the time varying boundary condition is defined. Here we used a file called prescribedMotion however you can name it whatever you want. The outOfBounds variable dictates what the simulation should do if the simulation time progresses outside of the time domain defined in the file.

The additional file containing the desired motion prescribedMotion is formatted in the following way:

(
( 0 (0 0 0))
( 0.005 (-0.0000308418795848531 0.00392682932795517 0))
( 0.01 (-0.0001233599085671 0.00785268976953207 0))
( 0.015 (-0.000277531259507496 0.0117766126774107 0))
( 0.02 (-0.00049331789293211 0.0156976298823283 0))
...
( 9.99 (-0.0001233599085671 -0.00785268976953189 0))
( 9.995 (-0.0000308418795848531 -0.00392682932795479 0))
( 10 (0 -3.06161699786838E-016 0))
)

The first column is the time in seconds, and the vector defines the point displacement. In the present tutorial, these points were calculated in libreOffice and then exported into the text file.  I arbitrarily made up the motions purely for the sake of making this blog post.

The circular motion was defined as:

x=0.25\cos\left(\pi t\right)-0.25 and y=0.25\sin\left(\pi t\right)

Decaying sinusoidal motion was:

y=\sin(\pi t) \exp\left(-t/2\right)

The rest of the set-up is identical to the set-up in the oscillating cylinder example. The solver pimpleDyMFoam is then run.

Results

Circular Motion

animation

Sinusoidal Decay

animation

Conclusions

This post demonstrated how a more complicated motion can be prescribed by using a little math and the timeVaryingUniformFixedValue boundary condition. Always like to hear questions and comments! Has anybody else done something like this?

 

Oscillating Cylinder in Laminar Crossflow – pimpleDyMFoam

In this post I am going to simulate an oscillating cylinder in a cross-flow… just for fun… and to provide an additional tutorial case for those wishing to use some of the dynamic meshing features of OpenFOAM.

The case I am going to simulate is a cylinder in a Reynolds number 200 cross-flow (U=2 m/s, D=1 m, nu = 0.01 m^2/s), oscillating at a rate of 0.2 hz.

Tutorial Files

The tutorial files for this case can be downloaded from here:

Download Tutorial Files

Please let me know if the download does not work or if there is a problem running the tutorial files. Note: I ran this case in parallel on a relatively fast 6-core computer. It will take a long time if you try to run it as single core.

Mesh Generation

For this simulation, I built a simple two dimensional O-grid around the cylinder which joined to a series of blocks making a rectangular domain. I built this using the native OpenFOAM meshing utility blockMesh (which I like a lot).

mesh.0000
Fig: Grid

If you are wondering about the blockMesh set up… I intend to a blockMesh tutorial post… not that it’s all that necessary since the OpenFOAM manual covers it pretty well.

Case Set-up

dynamicMeshDict

When running a dynamic mesh case a solver runs as part of the solution and solves for the new grid at each timestep. Several are available in OpenFOAM and I am not really trying to do a full post on that right now. So I’ll just tell you that in this case I decided to use the displacementLaplacian solver.

Along with the solver one must define the coefficients that go with the solver. For most of the solvers this means setting the diffusivity for the Laplace solver. Since I am relatively new to these types of solutions I did what we all do… and turned to the great CFD online forum! A discussion in this post (http://www.cfd-online.com/Forums/openfoam-meshing/97299-diffusivity-dynamicmeshdict.htm) made me think that a good starting point would be the inverseDistance model for the mesh diffusivity.

In the end my dynamicMeshDict (which belongs in the constant folder) looked like:

dynamicFvMesh dynamicMotionSolverFvMesh;

motionSolverLibs ( "libfvMotionSolvers.so" );

solver displacementLaplacian;


displacementLaplacianCoeffs
{
      diffusivity inverseDistance 1(cylinder);
}

The solver produces the grid movements through time as the simulation is performed. Here is what the mesh looks like while it is moving! I added the stationary mesh and the moving mesh. Click on them to get a clearer picture 🙂

Boundary Conditions

The set up for this case is simple. An inlet boundary (on the far left) where I specified the incoming velocity and pressure (both are uniform fixedValue). The top and bottom boundaries are slip boundaries (but could also be freestream depending on your preferred set-up style).

The cylinder itself requires some special treatment because of the moving mesh. It must obey the no-slip condition right? So do we set it as uniform (0 0 0) ? No! The cylinder is moving! Luckily there is a handy boundary condition for this in the U file:

 cylinder
 {
      type movingWallVelocity;
      value uniform (0 0 0);
 }

For a typical simulation using pimpleFoam a p and U file would be all that are required. However, since we are doing a moving mesh simulation there is another parameter that must be solved for and requires boundary conditions… pointDisplacement.

For the pointDisplacement boundary conditions, we know that all of the outer edges should NOT move. Therefore they are all fixed with a type of fixedValue and  a value of uniform (0 0 0).

The cylinder however is moving and requires a definition. In this simulation we are simulating and oscillating cylinder. Since we are using the displacement solver the type is oscillatingDisplacement. We input and omega (rad/s) and an amplitude (m) in the following way:

 cylinder
 {
 type oscillatingDisplacement;
 omega 1.256; 
 amplitude (0 0.5 0); // Max piston stroke
 value uniform (0 0 0);
 }

Results

Yay! Now its time to look at the results. Well since I am not doing this for any particular scientific study… let’s look at some pretty pictures!

Here is an animation of vorticity:

vort
Wake of Oscillating Cylinder

Looks pretty nice! I personally have a big nerd love for vortex shedding…. I don’t know why.

Obviously if you intend to do any scientific or engineering work with this type of problem you would need to think very carefully about the grid resolution, diffusivity model, Reynolds number, oscillation frequency  etc. All of these were arbitrarily selected here to facilitate the blog post and to provide a nice tutorial example!

Conclusion

In this post I briefly covered the set-up of this type of dynamic meshing problem. The main difference for running a dynamic mesh case is that you require a dynamic mesh solver (you must specify in the dynamicMeshDict), and you also require boundary conditions for that solver.

Let me know if there are any problems with this blog post or with the tutorial files provided.

Shocktube – rhoCentralFoam TVD Schemes Test

Careful selection of TVD interpolation schemes is important for solving high speed compressible flows. In cases with discontinuities such as shockwaves and contact surfaces, these schemes help keep the simulations free of spurious oscillations. In this post I will use the shock tube tutorial case to test some of the available schemes in OpenFOAM (specifically rhoCentralFoam).

I am going to test the following schemes: vanLeer, vanAlbada, and Gamma (0, 0.5 and 1). To test them, I will first simulate the tutorial case without changing anything. I will then examine how each of the limiters respond to an increase in grid density and also a decrease in CFL (Courant-Freidrichs-Levy) number … also the amazing Canadian Football League :).

Here is a summary of what I’ve found… then you can read through all the plots after!

  1. In all cases the results at a cell count of 100 showed oscillatory and therefore unsatisfactory results. As you should already know, grid resolution is important.
  2. At a CFL of 0.2, vanLeer in rhoCentralFoam was inconsistent. What I mean by this is when using vanLeer, as you refine your grid the results actually get worse! This is not desirable in a numerical scheme for obvious reasons.
  3. At a CFL of 0.1 vanLeer showed good results. This indicates that one should probably avoid using vanLeer unless you are committed to using small timesteps (more computing time).
  4. vanAlbada is consistent at both CFL numbers tested. This indicates that by using vanAlbada, it is possible that a larger CFL could be tolerated (assuming sufficient grid density is used).
  5. Gamma 0 – Corresponds to central differencing and therefore we expect it to be the least stable (and not total variation diminishing). This is confirmed by the results which are poor in all cases! Decreasing the CFL number helps but it is clear that different schemes should be used.
  6. Gamma 0.5 – Corresponds to hybrid central and upwind (TVD) differencing. The results are much better than Gamma 0.  In fact, the results are very similar to the results of vanLeer.
  7. Gamma 1 – Corresponds to fully TVD linear upwind interpolation. This scheme should be the most stable.  However, because it is only first order, higher grid resolution is required to get the best results. This is the compromise that must be made between stability and accuracy when using TVD schemes.

In my opinion, the two best options are vanAlbada and Gamma 1. vanAlbada most likely provides the best mix between accuracy and stability (at least for the case examined here). Gamma 1 provides the best option if stability is the most important thing to you. However, because it is first order you will need more cells to get the most accurate solution.

Details of the test case and all of the results are listed below! Hopefully somebody finds this useful. If I have anything wrong please let me know! If you want me to try more schemes or anything different let me know as well!

Cheers,

curiosityFluids

Shock-tube Case

If you aren’t familiar with the shock tube case, it is the simple one dimensional problem where one side (the driver side) begins at a high pressure. The other side (the driven side) is at a lower pressure. When the simulation begins, the discontinuous initial conditions send a shock wave down the tube. Following the shock wave is a contact surface. Here is an animation of the simulation:

 

 

vanLeer

CFL=0.2

  • 100, 1000, 10000 cells

CFL=0.1

  • 100, 1000, 10000 cell

vanAlbada

CFL=0.2

  • 100, 1000, 10000 cells

CFL=0.1

  • 100, 1000, 10000 cells

Gamma

Gamma 0

  • CFL=0.2

    • 100, 1000, 10000 cells
  •  CFL=0.1

    • 100, 1000, 10000 cells

Gamma 0.5

  • CFL=0.2

    • 100, 1000, 10000 cells
  •  CFL=0.1

    • 100, 1000, 10000 cells

Gamma 1

  • CFL=0.2

    • 100, 1000, 10000 cells
  •  CFL=0.1

    • 100, 1000, 10000 cells

 

Converging-Diverging Nozzle v2- rhoCentralFoam

Recently I read a CFD Online forum post where the accuracy of rhoCentralFoam was called into question.

The test case was the simple example of a converging diverging nozzle flow. In fact, I have already done a blog post regarding this. However, I tackled the problem using a large reservoir. In fact, it would have been more prudent to just simulate the nozzle and attempt to match the example from the NASA cfd benchmarking website:

http://www.grc.nasa.gov/WWW/wind/valid/cdv/cdv.html

So this is what I’ve done!

The case file can be found here:

Case File Download

The geometry in the case is an axi-symmetric converging diverging nozzle defined by a cosine function. I built the geometry and mesh in PointWise (just for the sake of time… I mean, I already did a detailed nozzle post). The result is here:

geometry
Fig: Geometry (generated using pointwise)

Three cases are calculated for comparison on the nasa site given above:

  1. Subsonic flow at exit (Pe/Po=0.89)
  2. Supersonic flow with shock wave in expanding section (Pe/Po=0.75)
  3. Supersonic flow at the exit (Pe/Po=0.16)

Case 1- Subsonic

For the subsonic case, it is important to remember that the outlet pressure MUST be set. Why? Because for a given nozzle geometry and stagnation pressure there is more than one possible solution since the flow is not choked!

The two most important boundary conditions here are the inlet and outlet. The inlet is specified using the totalPressure boundary condition (in my set-up P0=10000 Pa). You can also specify the total temperature. However, this will not affect the mach number or pressure profiles, only the velocity, density and temperature profiles. I set the temperature to be fixed at 298 K. The nozzle itself is an adiabatic, slip wall. The front and back are wedge types. The outlet is set to zeroGradient for velocity and temperature and the pressure is fixed to 8900 Pa (P/Po=0.89)

The results for pressure distribution and Mach number are shown here. Remember that since we are comparing P/Po and Mach number, it didn’t matter that I used a different pressure and temperature than the test case. I also included gifs showing the unsteady part of the simulation at the start.

 

Clearly the results match! Woohoo! One down, two to go!

Case 2- Supersonic with shock in diverging section

Tip for supersonic flow solution: Initialize the problem as a shock tube! Then the throat is choked right away in the simulation and convergence to the solution is faster!

In this case we must also specify the outlet pressure. Why? Because there are multiple solutions possible! Changing the back pressure changes the strength and location of the shock in the diverging section! So the setup is identical to the subsonic case for the inlet and nozzle. The only difference is that the outlet is set to 7500 Pa (P/Po=0.75).  In order to speed convergence, we can use the setFields utility to set the LHS of the nozzle to the stagnation pressure, and the RHS of the nozzle to the outlet pressure.

The results for pressure distribution and Mach number are shown here:

 

Hurray! Two for two!

Case 3- Supersonic Flow

The fully supersonic case is in fact the most simple. Recall that a supersonic nozzle is actually an initial value problem! This is in contrast to the subsonic cases, or cases with shocks, where they are boundary value problems (the outlet pressure controls the solution). But in order to achieve the supersonic flow solution we just set the outlet pressure to be zeroGradient. Additionally, we again start the solution as a shock tube.

So you might ask, but we were given a specified outlet pressure! Well, we CAN set this outlet pressure. But this unnecissarily restricts our solution. Having a back pressure very close to, but not exactly the same as the actual supersonic outlet pressure could make things difficult. It is much easier to let the solution converge to where IT THINKS the right answer is. Then, how closely we match the outlet pressure is in fact an indicator of the quality of our solution!

Here are the results:

Three for three! rhoCentralFoam has passed the test.

As always comment and correct me where my mouth (or keyboard) has made error.

Cheers,

curiosityFluids

curiosityFluidsLogo_r2

Turbulent Zero Pressure Gradient Flat Plate – simpleFOAM, komegaSST

An excellent test case, and case to familiarize yourself with some of the turbulence models available in OpenFOAM is a 2D flat plate with zero pressure gradient. I will solve this problem using the solver simpleFoam and the komegaSST model.

Here are the sections of this post:

  1. Quick Overview: kω-SST (komegaSST) Boundary Conditions
  2. Case set-up and mesh
  3. Results
    • u+ vs. y+
    • Coefficient of friction
  4. Conclusions and useful resources

Download the case file here: kOmega Flat Plate Tutorial File

Quick Overview: kω-SST (komegaSST) Boundary Conditions

In this section I will describe the boundary set-up for komegaSST  where no wall functions are implemented. This requires that the y+ along the wall is less than or equal to one. For the komegaSST turbulence model the boundary conditions are as follows:

At the wall:

  • ω (omega) – specific dissipation rate
    • BC type: fixedValue
    • BC value: \omega_{wall}=10\frac{6\nu_\infty}{\beta_1\left(\Delta y_{wall}\right)^2}
  • k – turbulent kinetic energy
    • BC type: fixedValue
    • BC value: 0
  • nut – turbulent viscosity
    • BC type: fixedValue
    • BC value: 0

In the free-stream:

  • ω (omega) – specific dissipation rate
    • BC type: fixedValue
    • BC value: \frac{U_\infty}{L} < \omega_{\infty} < \frac{10 U_\infty}{L}
  • k – turbulent kinetic energy
    • BC type: fixedValue
    • BC value:\frac{10^{-5}U_\infty^2}{Re_L} < k_{\infty} < \frac{0.1U_\infty^2}{Re_L}
  • nut – turbulent viscosity
    • BC type: calculated
    • BC value: 0 (this is just an initial value)

where \beta_1=0.075, and \Delta y_{wall} is the normal distance from the wall to the first cell center.

Case set-up and mesh

Free-Stream Properties

For this case I have followed a similar set up to the 2D flat plate case used on the NASA turbulence modelling resource website. By doing this it gives me something to compare to! However, I am not going to use the exact set-up on that website. Since we are using simpleFoam, I am going to set up this case the way I prefer which is to use a velocity (U) of 1 m/s and scale all other properties accordingly. The simulation properties that I used are :

  • U∞=1 m/s
  • ν=4e-7 m2/s
  • L=2 m

These correspond to a Reynolds number at L=2m of 5 million .

Grid Generation

The grid used was generated in blockMesh. High inflation was used in the boundary layer region in order to achieve the desired y+ value of less than 1. For more details on grid generation using blockMesh see the OpenFOAM manual!

Grid_BoundaryTypes
Grid

Boundary Conditions

For the incompressible solver simpleFoam, the minimum boundary conditions required for a simulation are p and U. However, if the simulation is a RANS simulation additional boundary conditions are required. For the kω-SST model we need to have a boundary condition on k and ω as well. The boundary conditions I defined in the zero (0) folder can be found in the attached tutorial file.

The only boundary condition that really needs any comment is omega. We calculate omega using \omega_{wall}=10\frac{6\nu_\infty}{\beta_1\left(\Delta y_{wall}\right)^2}. In our case the wall distance to the first cell center is 5.48316E-06. Using our free stream viscosity of 4E-7 this gives a value of omega at the wall of 10643595.39.

Tip for fvSolution

If you find that the results you are getting are wrong, it could be that the residuals for the different properties are too high! Certain properties converge before others and therefore you need to ensure that they all converge to a sufficiently low value!

Results

First we compare the coefficient of friction to the .dat file available from the .DAT available from the NASA Turbulence modelling resource. NOTE: We simulated at Reynolds number of 5 million whereas the NASA setup is at 10 million. So the x coordinate in the following is plots is rescaled and in fact we are only using half of the data from the .DAT file!

kOmegaCf
Coefficient of Friction

We can see from the figure that the coefficient of friction from our simulation matches the expected data closely! Hurray!

Next let’s compare the u+ vs. y+ profile to the  universal profile for turbulent boundary layers. Recall that u+ is the velocity normalized by the friction velocity (u_*), and y+ is given by the following equation:

y^+ = \frac{u_* y}{\nu}

The u+ is vs y+ is plotted here:

kOmegauy
u+ vs y+

We can see from the figure that our solution is pretty good! The y+ value of the first node is located around a y+ of approximately 0.5, the viscous sublayer matches very closely and the log law layer is not significantly off!

Conclusions

In this post we simulated a zero pressure gradient flat plate at a Reynolds number of 5 million. We compared the results for shear stress to the NASA turbulence modelling resource expected results and showed close agreement. Then the u+ vs y+ profile was compared to the universal law of the wall and again the results were okay!

Some useful references

The NASA Turbulence Modelling Resource, http://turbmodels.larc.nasa.gov/