Thursday, October 24, 2019

Density-Based Approach vs Pressure-Based Approach

Let us consider unsteady Euler equations along with a compressible equation of state (eos) - for example ideal gas.


  • the solution at t=0 along with
  • appropriate boundary conditions at all physical boundaries (at all times), 
calculate the "time-accurate" solution for t>0.

Time Marching

Each transport equation (continuity + momentum + energy) has a time derivative. So the simplest approach is to march each transport equation in time - explicitly or implicitly. That gives us {rho, rho*u, rho*v, rho*w, rho*E} at next time level. {rho, rho*u, rho*v, rho*w, rho*E} can be translated to {rho, u, v, w, T}:
  • {rho, rho*u, rho*v, rho*w, rho*E} => {rho, u, v, w, E}
  • {rho, u, v, w, E} => {rho, u, v, w, e} using E = e + (u*u + v*v + w*w)/2
  • {rho, u, v, w, e} => {rho, u, v, w, T} using e = Cv*T
So we have solved for 5 primary variables by marching the transport equations in time. What about the 6th primary variable (p)? p is obtained using the eos. This numerical approach is often termed as density-based approach because density is "solved" for by time marching the continuity equation while pressure is "updated" using eos.

What if we have a Constant Density EOS?

With constant density eos, the time derivative term in continuity equation disappears. What can we do now? We still have time derivatives in momentum and energy equations. We can march these equations in time to get {u, v, w, E} at next time level. {u, v, w, E} can be translated to {u, v, w, T}:
  •  {u, v, w, E} => {u, v, w, e} using E = e + (u*u + v*v + w*w)/2
  •  {u, v, w, e} => {u, v, w, T} using e = Cp*T
Now we have solved for 4 primary variables by marching the transport equations in time. rho is already known but we still need to compute p. Unfortunately the eos (rho = constant) does not help here in updating p! But there is one transport equation that we haven't utilized yet - the continuity equation. Although there is no pressure term in continuity equation, it can be applied to momentum equation to obtain a Poisson equation for pressure. This pressure Poisson equation can be used to update p. This numerical approach is called pressure-based approach because pressure is the primary variable (as opposed to density) which is "solved" for using a transport equation. There, however, is an important thing to note here. When momentum equations are marched in time to compute velocity they have to use the old pressure. This velocity satisifies momentum equation but does not satisfy continuity equation. After this, the pressure Poisson equation (which is representative of continuity equation) is solved to update pressure. So the momentum equations and pressure Poisson equation need to be solved iteratively several times one after another at each time level.

Let us redefine Density-Based Approach

We said that density-based approach is called so because we solve for density (and not pressure) using the continuity transport equation. But that is not entirely true in modern perspective. Today, a method can be said to follow density-based approach if it has following features:
  • all the transport equations (continuity + momentum + energy) are marched in time (or pseudo-time if only interested in steady state solution) either explicitly or implicitly, and
  • convective fluxes are computed using a Riemann solver based on the eigenstructure of all the transport equations (continuity + momentum + energy) taken together.
Note that this new definition does not put any emphasis on solving for density using a transport equation and updating pressure using eos. Even if a method solves for pressure using a transport equation and updates density using eos it can be classified as density-based approach if it fulfills above criteria.

Low Mach Number Flows

Let us get back to compressible ideal gas eos. As mentioned above, a density-based approach needs to compute convective fluxes using a Riemann solver based on the eigenstructure of all the transport equations (continuity + momentum + energy) taken together. Even for multi-dimensional problems, usually a 1D Riemann solver is used. The eigenvalues of the system used in 1D Riemann solver are {un, un, un, un + c, un - c} where un is the face-normal velocity and c is the speed of sound. At low Mach number, un << max{|un + c|, |un - c|}. This huge disparity in eigenvalues forces extremely small time steps to be used for explicit time marching even when one is not interested in disturbances traveling at acoustic speeds. This makes the convergence to steady state extremely slow. Moreover, in the limit of vanishing Mach number, discretization based on density-based approach introduces error that grows inversely with Mach number. This is true regardless of whether one uses explicit or implicit time marching.

Although the pressure-based approach was originally designed for constant-density equation of state, over the years it has been modified and adapted to be used for compressible equations of state as well. Because pressure-based approach does not use convective fluxes based on eigenstructure of all the transport equations (continuity + momentum + energy) taken together, it does not face the same problems as density-based approach. Because of this, pressure-based approach is widely used to solve low Mach number flows for compressible equations of state.

Current State

As mentioned above, pressure-based approach has been adapted over years to be to be used for compressible equations of state as well. The state-of-art methods based on pressure-based approach can be applied to the whole range of flows involving compressible equations of state, not just limited to low Mach number flows.

The same can be said about state-of-art methods based density-based approach. Over the years, density-based approach has also evolved. Current methods based on state-of-art density-based approach can be applied to low Mach number flows with compressible eos as well as flows involving constant-density equation of state.

Thursday, October 17, 2019

Working Pressure in Momentum Equation

Momentum equation in (Navier-Stokes equations) can be written as

Pressure in the above equation is the thermodynamic/mechanical (assuming Stokes hypothesis). But CFD codes in general don't use thermodynamic pressure in momentum equation discretization. They instead define a working pressure and use that. Here are the steps that take us to the final equation that is used for discretization.

Step 1

Introduce a constant reference density.


Introduce a constant reference pressure such that




  • Step 1 helps in decreasing the magnitude of gravity source term that can otherwise have some destabilizing effects. Step 2 is done to avoid precision issues. In general, the variation in pressure is orders of magnitude smaller than the mean value of pressure. If we work with absolute magnitude of pressure in discretization, we will waste most of the significant digits and may not capture variation in pressure to desired accuracy.
  • In a stagnant (constant density) liquid column, if we chose reference density to be liquid's density then p_therm will increase linearly as we go deeper but p_working will be constant.

Wednesday, October 16, 2019

Incompressible Flow

Incompressible Flow Definition

Continuity equation can be written as

If the material derivative of density is 0 (or negligible) i.e. if the density of a material particle does not change (or does not change significantly) as it moves through the flow domain, then the divergence of velocity vector is equal to 0 (or close to 0). Divergence free velocity is generally taken as the definition of incompressible "flow".

Note that constant density fluid implies divergence free velocity but the reverse is not true. Flows involving constant density fluid are always incompressible. Flows involving non-constant density fluid can also be assumed to be incompressible or divergence free under appropriate operating conditions. There is a difference between incompressible "fluid" and incompressible "flow" and that must always be kept in mind. Incompressible "fluid" implies incompressible "flow" but incompressible "flow" does not always imply incompressible "fluid".

Compressible / Incompressible "Fluid" and Incompressible "Flow"

For any fluid, equation of state (eos) can be prescribed to express density in terms of other thermodynamic variables - pressure and temperature. Fluids with dependence of density on pressure are called compressible "fluids" and those with no dependence of density on pressure are called incompressible "fluids". Strictly speaking, no fluid is incompressible. There is always some amount of compressibility, however small that may be. But for many fluids, it is fair to neglect compressibility for a wide class of problems. For example, water's density does not change appreciably when subjected to a wide range of pressures. For these fluids, we can simply neglect the dependence of density on pressure. And in many cases, we can go a step further and assume density to be constant.
  • Constant Density "Fluid" As mentioned above, any flow involving constant density incompressible "fluid" can be termed incompressible "flow" (i.e. divergence free velocity).
  • Ideal Gas "Fluid" At low Mach numbers, density variation in flows involving ideal gas can be very small. So the assumption of incompressible "flow" (i.e. divergence free velocity) is fairly accurate for low Mach number flows involving ideal gas.
  • Other Compressible "Fluids" For any general eos, we have to look for flow regime where the density variation over the solution domain is negligible. Over that flow regime, the assumption of incompressible "flow" would hold in general.

OK! Now we can identify incompressible "flow" based on eos and operating condition but how does it help?
  • Simplified Theoretical Analysis With the assumption of incompressible flow, simple formula can be derived that are easy and simple to use for theoretical analysis e.g. Bernoulli's equation and simple definitions of total/stagnation temperature & pressure.
  • Faster/Robust Solution Approach For a compressible "fluid", all the flow + energy variables (rho, p, V, T) are tightly coupled via the set of Euler (or Navier-Stokes) equations and equation of state. In general, this set of tightly coupled equations are solved numerically using what is generally termed as density-based approach. The approach has inherent coupling of both convective and acoustics velocity scales. At moderate and high Mach numbers, the two velocity scales are comparable but low Mach numbers the convective velocity scale is much smaller than acoustic velocity scale. This wide separate between convective and acoustic scales presents several problems with density-based approach. So people started looking for alternatives and pressure-based approach was derived with the assumption of incompressible flow. Pressure-based approach decouples the solution of momentum equation (based on convective scales) from that of pressure (based on acoustic scales) and provides a much better (faster & robust) alternative to density-based approach for incompressible flows. So if we have identified that the problem under consideration involves incompressible "flow", the choice of which numerical approach to use is simple.

Thursday, September 19, 2019

Why do we need Prism Mesh in Boundary Layer?

Let us consider a cell-centered unstructured finite-volume CFD solver. In the core of simulation domain, away from wall boundaries, we can use hex, tet, or poly cells.

Why can't we use the same type of cells in boundary layer? In the core of simulation domain, gradients are mild or negligible. But in boundary layer close to a wall boundary, solution gradients are extremely high normal to the wall but almost negligible along directions parallel to the wall. To resolve these highly anisotropic gradients we have two options to choose from:
  1. either use isotropic cells (similar cell size in all directions) with cell length scale based on the requirements of resolving wall normal gradient,
  2. or use highly stretched anisotropic cells with small length scale (based on the requirements of resolving wall normal gradient) in wall normal direction and much larger length scale in other directions.
Option #1 leads to a much larger number of cells compared to option #2. These extra cells unnecessarily add to computational cost without improving accuracy in a noticeable way. That's why option #2 is preferred.

Can we use highly anisotropic hex, tet, or poly cells? Anisotropic tets and polys can have very bad skewness and non-orthogonality and solver has difficulty in handling them. (What do skewness & non-orthogonality mean and how is solver affected by them? We will leave that for a later post.) Anisotropic hexes, however, do not have the same problem. They can be highly squished along just one direction while maintaining zero skewness and full orthogonality. So
  • if the core mesh is hex, we use quad prims (stretched hexes - with quad bases) in boundary layer,
  • if the core mesh is tet, we use tri prisms (tri bases) in boundary layer, and
  • if the core mesh is poly, we use poly prisms (poly bases) in boundary layer.

Sunday, September 15, 2019

Brushy Creek Regional Trail

The full Brushy Creek Regional Trail (Cedar Park, TX) is around 10 miles long (one-way). Part of this trail extending from Brushy Creek Lake Park to Twin Lakes YMCA is around 6 miles round trip.

We cover parts of this trail often but today morning we did the whole 6 miles round trip. At a brisk pace it took us around 2 hours. This is a nice well-maintained concrete trail. Only the section near Brushy Creek Lake Park is paved. There are a lot of people walking up and down the trail over weekend mornings. Bikers share the same trail although there are several offroad shoots as well. The trail is mostly surrounded by trees on both sides. So even on a hot day you don't feel burnt down. The section close to Brushy Creek Lake Park, however, is directly exposed to the sun and can be difficult to walk after early mornings. Here are some photos from today's walk.

And here is a small video with soothing sound of flowing creek.

Saturday, March 30, 2019

Star Trek type Moon gif

Taken sometime in 2015 (or perhaps earlier) using Orion XT8i and Orion StarShoot Solar System Camera