Darcy–Forchheimer Porous Media Coefficients Calculator

Turn a measured pressure-drop curve into the viscous and inertial resistance coefficients your solver asks for: a closed-form least-squares fit of Δp/L = a·v + b·v² over up to five points, α, 1/α and C₂ in the units Fluent, OpenFOAM and Star-CCM+ want, the name of the box each goes in, and an honest statement of whether your data determined them at all. Plus the perforated-plate case from open-area ratio.

The model, the closed-form fit, and which box each number goes in

Which box each number goes in

What your data can and cannot determine

Perforated plates: the two limits, and how far apart they are

Fitting the two coefficients, and why R² does not tell you whether to trust them

Two numbers, and everybody gets one of them wrong. Every solver represents a porous region the same way: Δp/L = (μ/α)v + C₂·½ρv², a term linear in velocity for viscous drag on the pore walls and a term quadratic in velocity for the form drag and jet losses that dominate once the pore Reynolds number is past about 1. To use it you need the permeability α and the inertial resistance factor C₂, and what you actually have is a pressure-drop curve from a bench or a manufacturer’s datasheet. This page does the fit in closed form, prints both coefficients in the units every major solver asks for, names the box each one goes in, and — the part that matters — tells you whether your data actually determined them.

The three errors, in order of frequency. The first is forgetting to divide by the sample thickness. α and C₂ are per-unit-length quantities; if you fit Δp instead of Δp/L you scale both coefficients by the thickness, which for a 24 mm core is a factor of 42, and nothing in the solver will tell you. The second is using the interstitial velocity instead of the superficial one. Both coefficients are defined against the superficial velocity — volumetric flow divided by the whole frontal area, solid included — and swapping them scales α by the porosity and C₂ by its square. The third is fitting data taken in air and then using the coefficients as though they were properties of the medium measured in water: they are properties of the medium, but you extracted them by dividing out a particular μ and ρ, so the fluid the TEST ran in is the one to enter above.

Why the fit is fixed at five points, and why a blank row costs nothing. The expression language behind this page has no loops, so the number of data rows has to be fixed and the normal equations written out term by term. Five rows is enough for real bench data and the sums are short: a = (T₁S₄ − T₂S₃)/(S₂S₄ − S₃²) and b = (T₂S₂ − T₁S₃)/(S₂S₄ − S₃²), with S_k = Σv^k and the two cross-sums T₁ = Σ(Δp/L)v and T₂ = Σ(Δp/L)v². There is a pleasant accident in the algebra: every term of every one of those sums carries at least one factor of v, so a row entered with zero velocity contributes exactly nothing. Setting a velocity to zero does not approximately remove a point, it removes it algebraically. The point count, R² and the residuals are masked as well, so they count only the rows you used.

Two fits, because the standard one is not the better one. The fit every textbook and every solver manual describes is unweighted least squares on Δp/L against v, which is also what a spreadsheet trendline with the intercept forced to zero gives you. It has a known weakness: the residuals it minimises are absolute, so the highest-velocity point — where Δp is largest — dominates the fit, and the viscous coefficient, which is only visible at the low end, is fitted almost as an afterthought. The alternative is to divide through by v first and fit (Δp/L)/v = a + bv, an ordinary straight line that gives every point equal weight and is equally closed-form. Which is better is a measurable question, so: synthetic Darcy–Forchheimer data at five velocities were perturbed with proportional Gaussian noise and refitted 4,000 times at each of several noise levels. The equal-weight fit recovers the permeability 1.5 to 3 times more accurately at every noise level tested, recovers C₂ about 1.4 times more accurately, and roughly halves the chance of returning a negative coefficient. Both are offered; the unweighted one is the default because it is the one a reviewer will reproduce, and the gap between the two is a better uncertainty estimate than either R².

R² tells you about the curve, not about the coefficients, and the difference is enormous. This is the single most useful thing on the page. A five-point Darcy–Forchheimer fit can have R² = 0.9999 and a coefficient that is wrong by a factor of two, and the reason is structural rather than statistical: if one of the two terms accounts for only a few per cent of the pressure drop anywhere in your data, then the data contain almost no information about that term’s coefficient, and the fit will put it wherever the scatter pushes it while still passing beautifully through every point. The measured rule is simple enough to remember: the relative error in a coefficient is roughly your measurement scatter divided by that term’s share of the pressure drop. Two per cent scatter and a five per cent inertial share means C₂ is uncertain by about half its own value, and about one such fit in six returns a negative C₂ outright. This page prints the inertial share at your highest test velocity for exactly this reason, and the band above the result says plainly when a coefficient is not determined.

The Fluent manual’s own example fails this test, which is worth knowing before you copy it. The Fluent User’s Guide demonstrates the procedure on four points at 20, 50, 80 and 110 m/s and arrives at 1/α = −242,282 1/m² — a negative viscous resistance, which it flags and tells you to avoid. The cause is visible from the data: a pure quadratic through those four points, with no linear term at all, already reaches R² = 0.9938, and adding the linear term takes it only to 0.9998. The viscous coefficient is fitting the last six parts in a thousand of the variance, so its value is scatter and its sign is a coin toss. There is a second thing worth noticing. Refit those same four points here by unweighted least squares and you get a = −4.94333 and b = 0.28912, not the manual’s −4.33539 and 0.28296 — a 14 per cent difference in the viscous coefficient. The manual’s pair matches a fit of Δp/v against v (−4.358 and 0.28245) to within 0.6 per cent, so the published example used the equal-weight method without saying so. Three answers from four points, depending on how you fit them. That is not a criticism of the manual; it is the reason a page like this should show the fit quality and offer both fits rather than hand you one confident number.

Perforated plates and screens, where there is no curve to fit. A thin plate has essentially no viscous resistance: the loss is form drag on the jets issuing from the holes, so α is effectively infinite and all the resistance is inertial. The loss coefficient follows from a two-step argument that needs no empiricism at all — the flow accelerates into a vena contracta of area C_cβA, then suffers a Borda–Carnot sudden-expansion loss of ½ρ(v_vc − v)² on the way back out, which gives K = (1/(C_cβ) − 1)² on the approach velocity. Setting C_c = 1 is the thick-plate limit where the jet reattaches inside the hole. And here is the result that makes the section worth having: Idelchik’s published empirical form for a thin sharp-edged perforated plate, ζ = (0.707(1−β)^0.375 + 1 − β)²/β², is exactly algebraically identical to that same Borda–Carnot expression with C_c = 1/(1 + 0.707(1−β)^0.375). Not approximately: identically. So the handbook correlation is the elementary derivation in disguise, with a contraction coefficient that rises from 0.590 for an isolated small hole to 0.770 as the holes take over the plate — physically the right direction, and the isolated-hole limit 1/(1+1/√2) = 0.5858 sits close to the classical free-streamline value of 0.611. Convert to a solver input with C₂ = K/t, or better, use a porous-jump or thin-baffle boundary condition and give it K directly.

Ergun as a cross-check, never as a substitute. If your medium really is a packed bed, Ergun’s equation ties both coefficients to one particle diameter and one porosity: α = d_p²ε³/(150(1−ε)²) and C₂ = 3.5(1−ε)/(ε³d_p). This page inverts both and prints the two equivalent particle diameters. If they agree, your medium behaves like granular packing and you can extrapolate with some confidence. If they disagree by a large factor — 35 in the worked example above — it does not, and that is useful negative information rather than an error: a louvred-fin core, a wire-mesh screen or a perforated plate all produce far more inertial loss for their permeability than packing does. The wrong conclusion to draw is that one coefficient must be wrong; the right one is that you cannot use Ergun to fill in a coefficient you did not measure.

What to check in the model afterwards. Three things, in this order. That the porous zone thickness in your mesh equals the sample thickness you fitted with, because the coefficients are per unit length and this error survives every other check. That the direction is right: a real core is strongly anisotropic and the resistance across the fins is often a hundred times the resistance along them, so a single isotropic pair of coefficients models a core only if the flow is normal to the face. And that the velocity in the zone is near the range you fitted over, which you can read off the plot above. If your case runs at a velocity outside the data, you are extrapolating a two-parameter fit, which is usually safe upward and unreliable downward. For the approach flow itself, the pipe friction and pressure drop calculator gives the duct term this one does not, and the hydraulic diameter calculator handles the passage geometry if you want a Reynolds number for the pores rather than for the face.

Frequently asked questions

Related calculators

References

1. ANSYS Fluent User’s Guide, porous media conditions. The source of the model form used here: the momentum sink is S_i = −(ΣD_ij μv_j + ΣC_ij ½ρ|v|v_j), which for a homogeneous isotropic medium reduces to S_i = −(μ/α v_i + C₂ ½ρ|v|v_i), with the viscous resistance 1/α in m⁻² and the inertial resistance C₂ in m⁻¹ described as a loss coefficient per unit length. The guide also states that the porous jump boundary condition is a one-dimensional simplification of the same model and should be preferred where possible for robustness. Cited, not reproduced.
2. The worked example in that same guide — 28.4, 487.0, 1432.0 and 2964.0 Pa at 20, 50, 80 and 110 m/s, ρ = 1.225 kg/m³, μ = 1.7894×10⁻⁵ Pa·s, Δn = 1 m, quoted trendline Δp = 0.28296v² − 4.33539v giving C₂ = 0.462 and 1/α = −242,282 m⁻² — was used here as an independent test of this page’s closed-form fit, and the result is a finding rather than a confirmation. Unweighted least squares on those four points gives a = −4.94333 and b = 0.28912 (checked against an independent SVD least-squares solve, agreeing to eleven significant figures), hence 1/α = −276,256 and C₂ = 0.4720. The guide’s quoted pair is 14 per cent away from that in the viscous coefficient. It matches a least-squares fit of Δp/v against v (−4.35818 and 0.28245) to within 0.6 per cent, so the published trendline is the equal-weight fit rather than the unweighted one. Both fits are offered on this page and the difference is printed.
3. OpenFOAM, DarcyForchheimer.H and DarcyForchheimer.C in the porosity models. The header states the source as S = −(μd + (ρ|U|/2)f)U with d the Darcy coefficient in 1/m² and f the Forchheimer coefficient in 1/m, and the implementation carries the comment “leading 0.5 is from 1/2*rho” beside the line that halves f. Read directly from the published source rather than from secondary descriptions, because it establishes the one thing a user needs: OpenFOAM’s d and f are numerically identical to Fluent’s viscous and inertial resistances, so coefficients transfer between the two solvers unchanged.
4. S. Ergun, Fluid Flow through Packed Columns, Chemical Engineering Progress 48 (1952), 89–94, for Δp/L = 150μ(1−ε)²v/(ε³d_p²) + 1.75ρ(1−ε)v²/(ε³d_p). Inverted here into the two equivalent particle diameters printed as a cross-check: d_p = √(150α(1−ε)²/ε³) from the permeability, and d_p = 3.5(1−ε)/(ε³C₂) from the inertial factor, the second following from C₂ = 2b/ρ with b = 1.75ρ(1−ε)/(ε³d_p). Used only as a consistency check on a fitted pair, never to supply a coefficient. Cited, not reproduced.
5. I. E. Idel’chik, Handbook of Hydraulic Resistance, for the loss coefficient of a thin sharp-edged perforated plate, ζ = (0.707(1−β)^0.375 + 1 − β)²/β² on the approach velocity. Cited and not reproduced as a table, because this page does not need to: the expression was shown here to be algebraically identical to the elementary Borda–Carnot result K = (1/(C_cβ) − 1)² with C_c = 1/(1 + 0.707(1−β)^0.375), which is what the calculator evaluates. The identity was verified at β = 0.05 to 0.9 and reproduces Idel’chik’s values exactly (for instance β = 0.3 gives 19.3157 by both routes), and the implied contraction coefficient rises monotonically from 0.5905 to 0.7703, with the isolated-hole limit 1/(1 + 1/√2) = 0.58579 close to the classical free-streamline value π/(π+2) = 0.61101 for a slot.
6. L. Borda and L. Carnot’s sudden-expansion loss, ½ρ(v₁ − v₂)², which is classical and derivable from a control-volume momentum balance with no empirical content. Combined with a vena contracta of area C_cβA it gives the whole perforated-plate result on this page, including the thick-plate limit K = (1/β − 1)² at C_c = 1.
7. The uncertainty figures quoted on this page and in the second table were generated here by Monte Carlo, not taken from a reference: synthetic Darcy–Forchheimer data at five velocities, perturbed with proportional Gaussian noise at 2, 5 and 10 per cent, refitted 4,000 times at each level by both fit methods, and the median absolute error in each recovered coefficient reported. Three media were covered, spanning inertial shares from under 1 per cent to 38 per cent. The equal-weight fit recovered the permeability 1.5 to 3 times more accurately than the unweighted fit at every level tested, and roughly halved the incidence of a negative coefficient. The rule of thumb that falls out — relative coefficient error ≈ 1.3 × scatter / that term’s share of the pressure drop — fits the measured medians to within about 20 per cent across the cases tested.
8. National Institute of Standards and Technology, Special Publication 811, for the unit conversions used: 1 psi = 6 894.757 293 Pa, the conventional inch of water = 2.490 889×10² Pa (0.0254 m × 1000 kg/m³ × 9.80665 m/s² exactly), and 1 ft/min = 0.00508 m/s exactly. The darcy is taken as 9.869233×10⁻¹³ m², its definition from a viscosity of 1 cP, a gradient of 1 atm/cm and a flux of 1 cm/s. A US Government work.

Setup guidance, not validation. Correlations have ranges of validity and cell-count estimates are order-of-magnitude. A converged simulation is not a correct one. Full disclaimer at calcengines.com/disclaimer/