Incompressible Flow Analysis
Governing Equations, Pressure-Velocity Coupling & Solvers

It's definitely a practical guideline, not a sharp physics cutoff. The reasoning goes like this: for isentropic flow, the density ratio between a moving fluid and the stagnation point is \(\rho/\rho_0 = (1 + (\gamma-1)/2 \cdot Ma^2)^{-1/(\gamma-1)}\). Plug in Ma = 0.3 and you get about a 4.5% density variation — which most engineers are happy to ignore. The real test is whether your quantity of interest — force coefficients, pressure drop, flow rates — changes significantly when you switch from incompressible to compressible solver. For a ventilation fan spinning slowly, there's zero practical difference. For a turbofan blade tip moving at Ma = 0.7, you absolutely need compressibility. The 0.3 threshold is where errors become noticeable but not yet disastrous; it's been validated on so many benchmark cases that it became the standard rule of thumb.

Great point. The speed of sound in water is about 1,480 m/s — water flowing at 10 m/s has Ma ≈ 0.007, so for steady liquid flows the incompressible assumption is essentially exact. Water is so stiff elastically that it takes enormous pressures to measurably change its density. The exception is water hammer — rapid valve closure creates pressure waves that propagate at the acoustic speed through the pipe system, and there you genuinely need the compressible water hammer equations. Similarly, cavitation involves phase change and density jumps, so the simple incompressible model breaks down in those regions too. But for 99% of water hydraulics, pipe networks, and pump analysis, incompressible Navier-Stokes is perfectly correct.

Yes, it's just a convenient notation that groups the modified pressure. Since density \(\rho\) is constant in incompressible flow, you can define a "kinematic pressure" \(P = p/\rho\) [m²/s²] and rewrite the equation without \(\rho\) explicitly in the pressure gradient term. This is particularly useful in OpenFOAM — when you look at the output file p, you're actually reading kinematic pressure, not absolute pressure in Pascals. If you want the actual pressure in Pa, you multiply by density. It also clarifies that for incompressible flow, only pressure differences matter (the absolute level is arbitrary), which is why you need at least one pressure reference point in your boundary conditions — otherwise the pressure system is singular and your solver won't converge.

Pressure is absolutely a real physical force — it acts on every fluid surface. But in the incompressible limit, what determines its value is the mathematical requirement that the velocity field must stay divergence-free at every instant. Think of it this way: if you have a pipe that narrows, the flow must accelerate. The pressure gradient is what provides exactly the right amount of force to make that happen while also ensuring continuity. The pressure "adjusts itself" to whatever value is needed to prevent mass accumulation anywhere in the domain. In compressible flow, there's no such constraint — pressure is determined by the equation of state (\(p = \rho R T\)) and can evolve independently. In incompressible flow, density is fixed, so the equation of state drops out, and pressure is left doing double duty: both applying force and enforcing the constraint. That's the mathematical sense in which it acts as a Lagrange multiplier, and it's precisely why incompressible solvers need special pressure-velocity coupling algorithms.

The practical rule is: SIMPLE for steady state, PISO for time-accurate transient, PIMPLE when you need transient but want larger time steps. For a room air conditioning study where you just want the final steady airflow pattern — use simpleFoam (SIMPLE). No need for time accuracy; you iterate to convergence as fast as possible. For a vortex shedding problem behind a cylinder where you need to capture the shedding frequency accurately — use icoFoam or pimpleFoam with small time steps so the Courant number stays below 1 (PISO regime). For a complex transient HVAC startup with a big domain where small time steps would take forever — use pimpleFoam with nOuterCorrectors = 2 or 3 and push the Courant number to 5 or 10. You trade some temporal accuracy for feasibility. The most common mistake beginners make is using PIMPLE with nOuterCorrectors = 1 (which makes it pure PISO) but a Courant number of 20 — that's unstable because PISO isn't designed for that.

Under-relaxation damps the oscillations that come from the chicken-and-egg coupling between pressure and velocity. When you solve momentum with the current pressure guess, you get a velocity that doesn't satisfy continuity. When you correct pressure to enforce continuity, you change the driving force for momentum, so your velocity guess is now inconsistent with the new pressure. If you apply the full correction each time, these updates can overshoot and oscillate wildly — pressure swings high, drives velocity the wrong way, then pressure overcorrects, and so on. Under-relaxation says: "only move 30% of the way toward the corrected value this iteration." It's like a damped feedback loop rather than a bang-bang controller. Too much relaxation (low \(\alpha\)) means slow convergence; too little (high \(\alpha\)) means instability. The 0.7 / 0.3 split for velocity / pressure is a historical starting point that works well for most steady flows; for highly turbulent or buoyancy-driven flows you often need to reduce both further.

That's exactly what DNS is — but the required mesh resolution is brutal. The smallest turbulent eddies (Kolmogorov microscale) scale as \(\eta \sim Re^{-3/4} L\), so the number of grid points scales as \(Re^{9/4}\). For a pipe flow at Re = 10,000, you'd need roughly \(10,000^{9/4} \approx 10^9\) grid points, and you'd need to march in time with steps small enough to resolve the fastest turbulent oscillations. That's around \(10^{10}\) time steps for a flow-through event. On the world's fastest supercomputers today, DNS at Re = 10,000 in a pipe takes weeks. Your HVAC model at Re = 100,000? Completely infeasible as DNS with any foreseeable hardware. RANS reduces this to maybe a million cells solved in minutes on a laptop, by replacing the chaotic fluctuations with modeled turbulent stresses. That's the trade-off: extreme computational cost for accuracy vs. turbulence models with their uncertainties but tractable cost.

For wind engineering around buildings, the standard industry choice is k-ω SST (Shear Stress Transport). Here's why: building flows have large separated regions on the leeward faces and roof, recirculation zones in street canyons, and sharp-edge-triggered separation at corners. The k-ε model (standard or realizable) handles free-stream turbulence well but predicts separation point incorrectly and tends to underpredict separated regions. SST blends k-ω near walls (which it handles excellently with low-Re correction) and k-ε in the free stream, and activates a limiter that prevents over-prediction of turbulent shear stress in adverse pressure gradient regions — exactly where buildings make things hard. If you need really accurate pedestrian-level wind comfort data, consider SAS (Scale-Adaptive Simulation) or LES, but for design-level studies and wind load estimation, SST with a good mesh (y⁺ ≈ 1 or a wall function approach with y⁺ 30–300) is the practical standard.

Both are right — they're describing two different wall treatment approaches, and you need to be consistent with your turbulence model setup. If you're using a low-Reynolds-number turbulence model (like k-ω SST with the "low-Re correction" option, or the Spalart-Allmaras model in its standard form) that resolves the viscous sublayer directly, then you need \(y^+ \lesssim 1\) for your first wall cell. The model equations are valid all the way to the wall and require resolution of the viscous sublayer. If you're using wall functions (which model everything between the wall and the log-law region without resolving it), then you need \(y^+ \approx 30\mbox{–}300\) so that your first cell sits squarely in the log-law region where the wall function formula is valid. The dangerous zone is \(y^+ \approx 5\mbox{–}30\) — if you end up there, your first cell is in the buffer layer, neither the viscous sublayer formula nor the log-law is accurate, and your wall shear stress (and hence heat transfer) prediction will be poor. Pick one approach and stick to it throughout the mesh.

Use a flat-plate friction estimate — it's rough but good enough to set up the mesh before your first run. The steps are: (1) estimate the wall shear stress using a skin friction correlation, for example the flat-plate Schlichting formula \(C_f \approx 0.026\,Re_L^{-1/7}\) for turbulent flow; (2) compute the friction velocity \(u_\tau = \sqrt{C_f/2} \cdot U_\infty\); (3) compute the first-cell height from the target \(y^+\): \(\Delta y = y^+\,\nu/u_\tau\). For example, for air at 10 m/s over a 1 m plate: \(Re = 10/1.5\times10^{-5} \approx 667,000\), \(C_f \approx 0.026 \times 667000^{-1/7} \approx 0.0039\), \(u_\tau \approx \sqrt{0.0039/2} \times 10 \approx 0.14\) m/s, and for \(y^+ = 1\): \(\Delta y = 1 \times 1.5\times10^{-5}/0.14 \approx 0.11\) mm. After you run, check the y+ distribution on walls and iterate the mesh if needed. There are online y+ calculators that automate this — just make sure you understand the input parameters rather than blindly trusting the numbers.

Three things matter most. First, the mesh in the wake region: drag has large contributions from pressure recovery (or lack of it) in the separated wake behind the car. You need sufficient cell density in the wake to capture the size and structure of the recirculation region — not just around the car surface. A mesh that looks fine on the car body but has huge cells 1 metre downstream will give poor drag. Second, turbulence model choice: standard k-ε overestimates the turbulent kinetic energy in the stagnation region on the front hood and underestimates separation on the rear. k-ω SST with curvature correction, or Realizable k-ε, gives significantly better Cd predictions. For road car aerodynamics, industry uses DES or SAS increasingly for accurate time-averaged drag. Third, moving ground and rotating wheels: if you model the car with a stationary floor and stationary wheel surfaces, you're modeling something physically different from a car on a road. The ground boundary layer builds up under the car floor, pushing air outward, and stationary wheels dramatically underpredict wheel wake drag. Even for steady simulations, a moving floor BC and spinning wheel surfaces (using MRF or AMI patches) are essential for drag accuracy within ±5% of wind tunnel.

OpenFOAM is genuinely used in industry, perhaps more than most students realize. The automotive sector uses it extensively — BMW, Volkswagen, Toyota have large internal CFD teams running OpenFOAM for aerodynamics and thermal management, partly because the licensing cost scales from zero (for ESI OpenCFD version) and partly because it's fully customizable at the source level. Aerospace uses it for research-grade computations. Wind energy companies use it for large-domain wind farm simulations where the commercial license cost for thousands of CPU cores would be prohibitive. The main reason Fluent and Star-CCM+ dominate in smaller companies and consultancies is the GUI — setting up a case in Fluent's point-and-click interface takes 2 hours; doing the same in OpenFOAM with text-based dictionaries takes 2 days for a newcomer. Once you know OpenFOAM well, the speed difference shrinks significantly. The honest answer is: for specialized or large-scale work, OpenFOAM is excellent and widely used; for consultancy work where client deliverables need to be generated quickly by various skill levels, commercial tools with their GUIs win.

That depends entirely on what "stalls" means. If residuals reach 10⁻³ and level off flat — never going higher, never going lower — that often means the flow is physically unsteady (vortex shedding, oscillating separation) but you're trying to solve it with a steady-state solver. simpleFoam is looking for a steady solution that doesn't exist; the "converged" 10⁻³ result is a time-averaged nonsense. You should check whether the residuals oscillate slightly (a sign of unsteadiness) and whether key monitor quantities like drag or pressure drop are still changing iteration-to-iteration. If so, switch to pimpleFoam and run as transient, or use a time-averaging approach. If residuals genuinely stall because they drop smoothly to 10⁻³ and sit there stably, with monitor quantities flat, that's often acceptable — especially for turbulence quantities like k and ε which frequently don't converge below 10⁻³ in complex geometries. The key diagnostic is whether your engineering quantities of interest have converged, not whether the residuals have hit an arbitrary target.