1 · The problem
Two flat plates, one moving. What is the mean velocity between them?
Plane Couette flow is the simplest wall-bounded shear flow: a fluid of kinematic viscosity ν fills the gap H between two parallel plates, the lower one at rest and the upper one moving at Ue. There is no pressure gradient; the moving wall drags the fluid along by friction alone. Once the flow is fully developed, nothing varies along the plates, the wall-normal velocity is zero, and the streamwise velocity u depends on y only. Laminar, the answer is a straight line. Turbulent, it is not, and the question of the dissertation is how well a Reynolds-averaged turbulence model predicts that profile.
The steady, incompressible Reynolds-averaged momentum equation reduces to one line: the total shear stress, viscous plus turbulent, is the same at every height. With an eddy-viscosity closure νt for the turbulent part,
(1)d/dy [ (ν + νt) du/dy ] = 0, u(0) = 0, u(H) = Ue.
Everything interesting is in νt. The mixing-length model makes it algebraic, νt = ℓ² |du/dy| with the van Driest length ℓ = κ y (1 − e−y⁺/26) measured from the nearest wall. Transport models (k–ε, Spalart–Allmaras, k–ω SST) add one or two more differential equations for the turbulence itself. Results are compared in wall units, u⁺ = u/uτ and y⁺ = y uτ/ν, where uτ = √(τw/ρ) is the friction velocity, because in those units the near-wall profile is universal: u⁺ = y⁺ in the viscous sublayer and u⁺ = (1/κ) ln y⁺ + B in the log layer.
One flow, three Reynolds numbers
The 2000 report defines Re = Ue H/ν with the full gap and the moving-wall speed. The DNS papers it is compared against — Bech et al. (1995) and Lee & Moser (2018) — put the walls at ±Uw and use the half gap h: Re = Uw h/ν, which for the same flow is exactly one quarter of the report's number. Modellers prefer a third, the friction Reynolds number Reτ = uτ h/ν, which is what the 2026 code takes as input. Bech's “Re = 1300” is Re = 5200 in the report's convention and Reτ = 82.
2 · The conclusion
The 2000 conclusion holds — and it turns out to be stronger than the report claimed.
The report concluded that a mixing-length model “can provide a reasonable estimate of the mean velocity for a one-dimensional, incompressible turbulent Couette flow for a range of Reynolds numbers up to 20 000”, that the one-dimensional assumption is validated by the agreement with experiment, that the law of the wall is a good initial guess, and — with disarming honesty — that “unfortunately, the k–ε turbulence model code did not function properly before the end of the project and therefore its results could not be analysed.”
Re-run against DNS data that did not exist in 2000, with a solver that converges every model to round-off, the picture in 2026 is this:
The simplest model is the best one
Against Lee & Moser's DNS at Reτ = 93, 220 and 500, the van Driest mixing length predicts the wall speed Uw⁺ to within 1.2–2.5 % and the friction coefficient to within 2.4–4.9 %. Menter's k–ω SST is next (5–9 % on Cf), Launder–Sharma k–ε 12–13 %, Spalart–Allmaras 16–20 %.
The k–ε code was fixable, not wrong
The Launder–Sharma model the report could not stabilise converges in ~110 iterations and 0.2 s once the destruction terms are linearised into the matrix diagonal and the iteration is damped in pseudo-time. It then gives the answer it always would have given: a friction coefficient about 13 % too low.
The 2000 numbers were right, to a point
A line-by-line port of the Fortran program agrees with the new code to 0.007 % in Cf at equal κ, and its Thomas-algorithm test table is reproduced digit for digit. Its one numerical flaw — a one-sided velocity gradient that made the scheme first-order and slightly asymmetric — cost mesh points, not accuracy: 0.2 % on 250 points, 0.0005 % on the 10 000 it used.
All four closures over-predict the wall speed in wall units, i.e. under-predict friction, and the error grows with Reynolds number. The physical reason is visible in the DNS: Couette flow has no wake region, and its core is filled with very long streamwise rolls (Lee & Moser's “extreme-scale motions”) that carry more momentum than any local eddy-viscosity model allows. The mixing length wins here not because it is more sophisticated but because a near-wall-calibrated κ y happens to be a good description of a flow that is “all wall”.
3 · What was done in 2000
A Fortran program, a Poisson test, and graphs read off the page with a ruler
The report (48 pages, June 2000, Aeronautics 4th year, acknowledging Prof. Micha Wolfshtein and Dr. B. Aupoix for their help and Dr. R. E. Brown for his support) walks through the structure of a turbulent boundary layer, the two channel-flow archetypes, the mixing-length model and six low-Reynolds-number k–ε variants (Jones–Launder, Launder–Sharma, Chien, Nagano–Tagawa, Myong–Kasagi, So), the law of the wall derived from first principles, and a comparison against one DNS and two experimental data sets. Its numerical core is compact:
FormulationEquation (1) non-dimensionalised by Ue and H: d/dη[(1/Re + n) dz/dη] = 0 with n = κ²η²(1 − e−Re·γ·η/26)² |dz/dη|, γ = √(Cf/2).
DiscretisationCentral second-order on a uniform mesh of 10 000 points; coefficients averaged to the half-points; the velocity gradient inside νt by a backward difference.
IterationGuess Cf = 2/Re (laminar), compute νt, assemble the tridiagonal system, solve with the Thomas algorithm adapted from Numerical Recipes, recompute Cf from the wall gradient, repeat until |ΔCf| and max|Δz| < 10⁻⁸.
Output
SOLUTION.DAT with y⁺ and u⁺, plotted in Excel against the log law and the digitised data.What the report tested, and what those tests say now
Three checks were made in 2000. The Thomas solver was verified on z″ = (p+1)(p+2)ηp, whose exact solution is ηp+2, for meshes from 10 to 10 000 points (Table 1 of the report). The solution's symmetry about mid-gap was measured (Table 2). And the first mesh point was confirmed to lie inside the viscous sublayer, y⁺(1) < 1 (Table 3). Re-running the port:
| Check | 2000 report | 2026 port of the same program | Reading |
|---|---|---|---|
| Thomas test, max error, N = 10 000, p = 2 / 3 | 2.50E-9 / 6.41E-9 | 2.50E-9 / 6.42E-9 | Digit-for-digit. The solver was right; the residual error is the second-order truncation of the Poisson stencil, not the solver. |
| y⁺ of first mesh point, Re = 3 000 … 20 000 | 0.010, 0.03, 0.02, 0.028, 0.04, 0.026 | 0.010, 0.016, 0.021, 0.028, 0.041, 0.052 | Four of six agree; the entries for Re = 5 000 and 20 000 look like transcription slips in the 2000 table. All are far inside the sublayer. |
| Asymmetry, max|u⁺ − u⁺mirrored|, Re = 20 000 | 0.888 | 0.016 | Same trend (grows with Re), much smaller magnitude. The asymmetry is real and comes from the backward difference in νt; the 2000 figure may have been measured on the plotted, not the computed, data. |
| Iterations to converge, Re = 3 000 / 20 000 | “under a minute”, Pentium 133 | 219 / 1 316 iterations; at Re = 3 000, 3.4 s with a Fortran-style Thomas loop in Python and 0.16 s with LAPACK; 0.9 s at Re = 20 000 | Convergence count grows linearly with Re: Picard iteration on Cf is a fixed-point map whose contraction weakens as the wall layer thins. |
Reading it twenty-six years later
- The Reynolds number of the Bech comparison was mismatched by a factor of four. Graph 2 compares the program at “Re = 1300” with Bech's DNS at “Re = 1300”, but Bech's number uses the half gap and half the wall-speed difference. In the report's convention that DNS is Re = 5 200, Reτ = 82; the program run was at Reτ ≈ 26, barely above the point where plane Couette flow sustains turbulence at all (Re ≈ 360 in Bech's units, 1 440 in the report's). The comparison still looked good because in wall units the near-wall profile is universal — which is the law of the wall doing the work, not the model. Re-run at the right Reynolds number, the mixing length gives Uw⁺ = 16.0 against Bech's 15.8.
- The k–ε failure had a specific cause. The report notes that the ε destruction term was moved into the matrix diagonal for stability “but this didn't solve the problem”. The k equation's dissipation was left explicit on the right-hand side (equation 47 of the report), the ε destruction was linearised as εold·ε rather than Newton-style, and there was no under-relaxation or pseudo-time step. Any one of these can make the segregated iteration oscillate between two states, which is exactly the symptom described (“oscillating from iteration to iteration between two different kinds of solutions”). Section 4 shows the fix.
- Ten thousand uniform points was the wrong economy. A stretched mesh of 300 points reaches the same accuracy; with a first-order gradient, the uniform mesh needed 10 000 to bury the error. The “possible improvement” listed in the report — a variable step size — was the right instinct.
- The data were read off printed graphs by hand, and the report says so twice. Today all three DNS data sets used here are text files with a citation header; the experimental sets of 1987 and 2000 remain trapped in paper.
- What was right: the formulation, the discretisation of the diffusion operator, the Thomas solver, the convergence criteria, the choice of the mixing length as the workhorse, the insistence on testing the solver on a problem with a known solution, and the conclusion.
4 · The 2026 rewrite
Same equation, written so that the solver has less to do and the models have nowhere to hide
Formulation: the momentum equation is a quadrature
Integrating equation (1) once gives (ν + νt) du/dy = uτ², a constant. In wall units this is
(2)(1 + νt⁺) du⁺/dy⁺ = 1, so u⁺(y⁺) = ∫₀y⁺ ds / (1 + νt⁺(s)).
Given the eddy viscosity, the velocity is an integral — no tridiagonal system, no friction-coefficient iteration. The friction Reynolds number Reτ is the natural input; the wall speed Uw⁺ = u⁺(Reτ) and therefore Re and Cf = 2/Uw⁺² are outputs. If one wants a given Re instead, as the 2000 program did, a one-dimensional root find (Brent's method, about eight solves) maps Re to Reτ.
For the mixing length the closure is algebraic, νt⁺ = ℓ⁺² du⁺/dy⁺, and (2) becomes a quadratic in the gradient with the positive root
(3)du⁺/dy⁺ = 2 / (1 + √(1 + 4ℓ⁺²)).
The whole 2000 program — its 10 000-point matrices, its 200–1 300 Picard iterations, its convergence tests — collapses to two lines of NumPy and a cumulative trapezoid. The profile is symmetric to 10⁻¹³ by construction. That is not a criticism of 2000; it is the kind of thing one sees only after writing the general solver and noticing that momentum never needed one.
Numerics for the transport models
The three transport models — Launder–Sharma k–ε (the variant from the report's own table), Spalart–Allmaras with the 2012 S̃ modification, and Menter's 2003 k–ω SST — do need a solver, and the solver is where the 2026 code differs most from the 2000 one.
- Mesh. y⁺ = Reτ[1 + tanh(β(2ξ−1))/tanh β] with β = 3.5 and 300 points: first spacing 0.09 wall units at Reτ = 500, 20 wall units at the centre. The diffusion operator is the conservative second-order stencil on a non-uniform mesh, checked by a manufactured solution (observed order 2.00, 1.95 on the coarsest pair).
- Source linearisation. Every destruction term is split Patankar-style: the part that can only reduce the variable goes to the matrix diagonal, the rest to the right-hand side, so the update stays positive. Quadratic destruction (ε²/k, β ω², cw1 fw(ν̃/d)²) is linearised Newton-style, −2c φoldφ + c φold², which turns the oscillating map φ ← c/φ of the 2000 code into a contracting one.
- Stage 1, implicit pseudo-time. Each equation is stepped with backward Euler and a pseudo-time step that grows as the residual falls (switched evolution relaxation). The tridiagonal solves go to LAPACK's
gtsvthroughscipy.linalg.solve_banded; the eddy viscosity fed back to (2) is under-relaxed by 0.5. This is the 2000 algorithm made robust. - Stage 2, Newton. From a residual of 10⁻⁶, or as soon as stage 1 stops making progress, the coupled discrete system is solved by Newton's method with a sparse finite-difference Jacobian (seven residual evaluations per field, by column colouring) under pseudo-transient continuation. Spalart–Allmaras needs this: its Picard iteration settles into a limit cycle at a residual of 0.4 (the buffer-layer coupling between ν̃ and S̃), and Newton takes it from there to 10⁻¹³ in twelve steps.
- Residual. Convergence is judged on the steady residual of the discrete equations with all coefficients updated — not on the change between iterates, which goes to zero for any large pseudo-time step whether or not the equations are satisfied.
Verification
Verification asks whether the equations are solved right; validation, below, whether the right equations were solved. Three verification results, all in tests/: the report's Thomas-algorithm table is reproduced to the printed precision at every mesh size and polynomial degree; the manufactured-solution test on the stretched-mesh operator returns second order; and each model's wall speed converges as N⁻² to a Richardson-extrapolated value, so the default 300-point mesh carries an error of 0.06 % (mixing length, Spalart–Allmaras, k–ε) and 0.2 % (SST, whose ω wall condition depends on the first spacing).
Validation against DNS
Lee & Moser (2018) simulated plane Couette flow at Reτ = 93, 220 and 500 in boxes up to 100πh long, and publish the mean and fluctuation profiles. Those three cases span Re = 1 500 to 10 000 in their convention, 6 000 to 40 000 in the report's — the top of the 2000 range and twice beyond it.
| Model | Reτ = 93 Cf error | Reτ = 220 Cf error | Reτ = 500 Cf error | Iterations at 500 | Time at 500 | Character of the error |
|---|---|---|---|---|---|---|
| Mixing length (van Driest) | -2.5 % | -2.4 % | -4.9 % | 0 | 0 ms | Slightly too little mixing in the core; near-wall exact by calibration. |
| Launder–Sharma k–ε | -13.1 % | -11.9 % | -12.9 % | 111 | 156 ms | Turbulence switches on too late in the buffer layer; log-law intercept too high. |
| Spalart–Allmaras | -20.1 % | -15.6 % | -16.0 % | 166 | 300 ms | Same buffer-layer delay, plus a core where the wall-distance kink starves ν̃. |
| Menter k–ω SST | -5.5 % | -6.6 % | -9.4 % | 78 | 115 ms | Slightly too much mixing near the wall, too little in the core; best of the transport models. |
Then and now, in numbers
| Quantity | 2000 program (ported, LAPACK solve) | 2026 solver |
|---|---|---|
| Points in the wall-normal direction | 10 000, uniform | 300, tanh-stretched |
| Momentum solve | Tridiagonal system every iteration | One cumulative integral |
| Iterations, mixing length, Re = 20 000 | 1 316 | 0 (closed form); 8 solves to hit a given Re |
| Wall time, mixing length, Re = 20 000 | 0.9 s with LAPACK (a Fortran-style Thomas loop in Python is 20× slower) | 1 ms, including the root find for Re |
| Cf at Re = 20 000, κ = 0.40 | 1.36476 × 10⁻³ | 1.36467 × 10⁻³ (0.007 % apart) |
| Symmetry error, Re = 20 000 | 1.6 × 10⁻² | 1 × 10⁻¹³ |
| Order of accuracy | 1 | 2 |
| k–ε | Did not converge | 111 iterations, 0.16 s, residual 3 × 10⁻¹¹ |
| Everything on this page | — | 21 s, one command |
5 · The modernised dissertation
If the same project were set today
The problem is still a good one — it is the smallest wall-bounded flow in which a turbulence model can be wrong — and a modern version would keep its scope, one-dimensional and steady, but change what counts as a result. In 2000 the deliverable was a profile that agreed with a graph. In 2026 the deliverable is a verified solver, a validated set of models, and a quantified statement of how and why each model fails, all reproducible from one command. The proposed structure:
- Problem and reference dataPlane Couette flow; the three Reynolds-number conventions stated once; the DNS data sets of Bech (1995), Pirozzoli et al. (2014) and Lee & Moser (2018) as the ground truth, with their box sizes and sampling errors, since the core statistics of Couette flow converge slowly with box length.
- FormulationConstant total stress; velocity by quadrature; Reτ as the control parameter; the closures written in wall units so that ν = 1 and every model constant is dimensionless — mixing length, Launder–Sharma k–ε, Spalart–Allmaras, k–ω SST, and one Reynolds-stress model (Launder–Reece–Rodi or SSG) to see whether anisotropy helps in the core.
- Numerical methodStretched mesh, conservative operator, Patankar and Newton linearisation, pseudo-transient continuation, Newton with a sparse Jacobian; a section on why segregated iteration fails, with the 2000 k–ε as the worked example.
- VerificationManufactured solutions for every operator; order-of-accuracy tables; Richardson extrapolation and a grid-convergence index for each reported number; the pytest suite as an appendix rather than a paragraph.
- ValidationProfiles, Reynolds stress, friction law and the log-law constants (κ, B) fitted to each model and to the DNS, at three Reynolds numbers; error bars from the DNS sampling uncertainty; a clear statement of which model to use for what.
- Where the models fail and whyThe core of Couette flow: the DNS eddy viscosity at mid-gap exceeds κReτ, i.e. the large rolls transport momentum non-locally; a term-by-term budget of the k equation from the DNS against each model's; what a data-driven correction (field inversion on the DNS, a learned νt correction) can and cannot fix in one dimension.
- OutlookThe same four models in a three-dimensional finite-volume code (OpenFOAM or SU2) on the same case, to show that the 1-D result is the 3-D result; and the honest limit of RANS for this flow, set by the rolls.
What “appropriate today” means, beyond this problem
For a one-dimensional RANS problem the appropriate tool in 2026 is exactly what is used here: a few hundred lines of NumPy and SciPy, tests, a data directory with citation headers, and a script that regenerates every figure. Compiled languages are not needed until the problem is three-dimensional and unsteady, and then one does not write the solver either: OpenFOAM, SU2 and Nek5000 implement these same closures, and a dissertation's contribution moves to the case, the verification and the physics. Direct numerical simulation, which in 2000 was a single data set in one paper, is now a public database (the Oden Institute's and Johns Hopkins' collections) that any model can be tested against in an afternoon. And the frontier that did not exist in 2000 — closures augmented or replaced by machine learning, trained on that DNS — is best approached with precisely this kind of small, fully verified baseline, because the first question about a learned model is whether it beats a mixing length. Here, for this flow, nothing does.
The code is laid out for that next step. A closure is a class with three methods (initial fields, eddy viscosity, source terms); the solver does not know which one it is running. Adding a Reynolds-stress model or a learned correction is a new file, not a new program.
Appendix
Code, data and references
The package couette2026/ accompanies this page. python -m pytest -q runs the eight tests in 1.5 s; python run_all.py recomputes every number above in about 20 s; python figures.py redraws the figures. Requirements: Python 3.10+, NumPy, SciPy, matplotlib.
from couette_rans import solve, retau_from_re, MixingLength, MenterSST
r = solve(MenterSST(), retau=500) # friction Reynolds number on the half gap
r.ue_plus / 2, r.cf_half # wall speed U_w+ and Cf = 2 (u_tau / U_w)^2
r = retau_from_re(MixingLength(), 20000) # Re = U_e H / nu, the 2000 report's conventioncouette_rans/mesh.py · 31 lines
"""Wall-refined one-dimensional meshes.
In 2000 the mesh was uniform with 10 000 points so that the first point sat at
y+ ~ 0.03. A tanh stretching puts the same resolution at the wall with a few
hundred points and coarsens smoothly towards the centre where nothing happens.
"""
import numpy as np
def stretched_mesh(retau: float, n: int = 300, beta: float = 3.5) -> np.ndarray:
"""Return y+ coordinates of n+1 points on [0, 2*retau], clustered at both walls.
y+(xi) = Re_tau * (1 + tanh(beta*(2 xi - 1)) / tanh(beta)), xi in [0, 1].
beta ~ 3-4 gives a first spacing of order 1e-1 wall units for Re_tau ~ 500.
"""
xi = np.linspace(0.0, 1.0, n + 1)
return retau * (1.0 + np.tanh(beta * (2.0 * xi - 1.0)) / np.tanh(beta))
def uniform_mesh(retau: float, n: int) -> np.ndarray:
return np.linspace(0.0, 2.0 * retau, n + 1)
def wall_distance(y: np.ndarray) -> np.ndarray:
"""Distance to the nearest wall on the full gap."""
return np.minimum(y, y[-1] - y)
def gradient(f: np.ndarray, y: np.ndarray) -> np.ndarray:
"""Second-order derivative on a non-uniform mesh (central inside, one-sided at ends)."""
return np.gradient(f, y, edge_order=2)
couette_rans/models.py · 224 lines
"""Eddy-viscosity turbulence closures written in wall units (nu = 1).
Every transport model exposes the same small interface so the solver does not
know which model it is running:
fields names of the transported scalars
initial(y, dudy, nut0) initial fields from a mixing-length guess
eddy_viscosity(y, F, dudy) nu_t+ from the current fields
equations(y, F, u, dudy, nut)
-> {name: (gamma, S_C, S_P, (bc_left, bc_right))}
gamma : diffusivity at the nodes, so the equation reads
d/dy( gamma dphi/dy ) + S_C + S_P*phi = 0
S_C : explicit source (kept >= 0 wherever possible)
S_P : implicit source coefficient (kept <= 0) -- Patankar linearisation.
This is the single most important stabilisation trick and the one
the 2000 k-epsilon code was missing for the k equation.
"""
from dataclasses import dataclass
import numpy as np
from .mesh import wall_distance, gradient
EPS = 1e-30
# --------------------------------------------------------------------------
# Mixing length (van Driest 1956) — closed form for Couette flow
# --------------------------------------------------------------------------
@dataclass
class MixingLength:
kappa: float = 0.41
aplus: float = 26.0
name: str = "Mixing length (van Driest)"
fields: tuple = ()
def length_scale(self, y):
d = wall_distance(y)
return self.kappa * d * (1.0 - np.exp(-d / self.aplus))
def closed_form(self, y):
"""Because the total stress is uniform in Couette flow,
(1 + l+^2 |du/dy|) du/dy = 1 is a quadratic in du/dy with root
du+/dy+ = 2 / (1 + sqrt(1 + 4 l+^2)). No iteration is needed."""
l2 = self.length_scale(y) ** 2
dudy = 2.0 / (1.0 + np.sqrt(1.0 + 4.0 * l2))
nut = l2 * dudy
return dudy, nut
# --------------------------------------------------------------------------
# Launder & Sharma (1974) low-Reynolds-number k-epsilon
# --------------------------------------------------------------------------
@dataclass
class LaunderSharmaKE:
cmu: float = 0.09
ce1: float = 1.44
ce2: float = 1.92
sk: float = 1.0
se: float = 1.3
name: str = "Launder–Sharma k–ε"
fields: tuple = ("k", "eps")
def initial(self, y, dudy, nut0):
# local equilibrium P = eps with nu_t = Cmu k^2/eps => k = nu_t |u'| / sqrt(Cmu)
k = np.maximum(nut0 * np.abs(dudy) / np.sqrt(self.cmu), 1e-8)
eps = np.maximum(nut0 * dudy ** 2, 1e-8)
k[[0, -1]] = 0.0
eps[[0, -1]] = 0.0
return {"k": k, "eps": eps}
def eddy_viscosity(self, y, F, dudy):
k, eps = F["k"], F["eps"]
rt = k ** 2 / (eps + EPS)
fmu = np.exp(-3.4 / (1.0 + rt / 50.0) ** 2)
return self.cmu * fmu * k ** 2 / (eps + EPS)
def equations(self, y, F, u, dudy, nut):
k, eps = F["k"], F["eps"]
rt = k ** 2 / (eps + EPS)
f2 = 1.0 - 0.3 * np.exp(-rt ** 2)
P = nut * dudy ** 2
D = 2.0 * gradient(np.sqrt(np.maximum(k, 0.0)), y) ** 2 # extra dissipation at the wall
d2u = gradient(dudy, y)
E = 2.0 * nut * d2u ** 2 # Launder-Sharma E term
kk = np.maximum(k, EPS)
eq_k = (1.0 + nut / self.sk, P, -(eps + D) / kk, (0.0, 0.0))
# destruction -Ce2 f2 eps^2/k linearised Newton-style: -2 c eps_old*eps + c eps_old^2
c = self.ce2 * f2 / kk
eq_e = (1.0 + nut / self.se, self.ce1 * eps / kk * P + E + c * eps ** 2, -2.0 * c * eps, (0.0, 0.0))
return {"k": eq_k, "eps": eq_e}
# --------------------------------------------------------------------------
# Spalart & Allmaras (1992/1994), with the Allmaras-Johnson-Spalart 2012 S~ fix
# --------------------------------------------------------------------------
@dataclass
class SpalartAllmaras:
cb1: float = 0.1355
cb2: float = 0.622
sigma: float = 2.0 / 3.0
kappa: float = 0.41
cw2: float = 0.3
cw3: float = 2.0
cv1: float = 7.1
cv2: float = 0.7
cv3: float = 0.9
name: str = "Spalart–Allmaras"
fields: tuple = ("nuhat",)
@property
def cw1(self):
return self.cb1 / self.kappa ** 2 + (1.0 + self.cb2) / self.sigma
def initial(self, y, dudy, nut0):
nh = np.maximum(nut0, 1e-6)
nh[[0, -1]] = 0.0
return {"nuhat": nh}
def _fv1(self, chi):
return chi ** 3 / (chi ** 3 + self.cv1 ** 3)
def eddy_viscosity(self, y, F, dudy):
chi = np.maximum(F["nuhat"], 0.0) # nu = 1 in wall units
return chi * self._fv1(chi)
def equations(self, y, F, u, dudy, nut):
nh = np.maximum(F["nuhat"], 0.0)
d = np.maximum(wall_distance(y), 1e-12)
chi = nh
fv1 = self._fv1(chi)
fv2 = 1.0 - chi / (1.0 + chi * fv1)
omega = np.abs(dudy)
sbar = nh * fv2 / (self.kappa ** 2 * d ** 2)
# Allmaras et al. (2012) modification keeps S~ positive
stil = np.where(sbar >= -self.cv2 * omega,
omega + sbar,
omega + omega * (self.cv2 ** 2 * omega + self.cv3 * sbar)
/ np.maximum((self.cv3 - 2 * self.cv2) * omega - sbar, EPS))
stil = np.maximum(stil, 1e-12)
r = np.minimum(nh / (stil * self.kappa ** 2 * d ** 2), 10.0)
g = r + self.cw2 * (r ** 6 - r)
fw = g * ((1.0 + self.cw3 ** 6) / (g ** 6 + self.cw3 ** 6)) ** (1.0 / 6.0)
dnh = gradient(nh, y)
gamma = (1.0 + nh) / self.sigma
c = self.cw1 * fw / d ** 2 # destruction -c nh^2, Newton-linearised
# production cb1 S~ nh: the part cb1 (S~ - Omega) nh is negative where fv2 < 0
# (1 < chi < ~18, i.e. the buffer layer) and is treated implicitly there.
dS = self.cb1 * (stil - omega)
S_C = self.cb1 * omega * nh + np.maximum(dS, 0.0) * nh + self.cb2 / self.sigma * dnh ** 2 + c * nh ** 2
S_P = -2.0 * c * nh + np.minimum(dS, 0.0)
return {"nuhat": (gamma, S_C, S_P, (0.0, 0.0))}
# --------------------------------------------------------------------------
# Menter k-omega SST (2003 form)
# --------------------------------------------------------------------------
@dataclass
class MenterSST:
bstar: float = 0.09
a1: float = 0.31
kappa: float = 0.41
sk1: float = 0.85
sw1: float = 0.5
b1: float = 0.075
sk2: float = 1.0
sw2: float = 0.856
b2: float = 0.0828
name: str = "Menter k–ω SST"
fields: tuple = ("k", "omega")
@property
def g1(self):
return self.b1 / self.bstar - self.sw1 * self.kappa ** 2 / np.sqrt(self.bstar)
@property
def g2(self):
return self.b2 / self.bstar - self.sw2 * self.kappa ** 2 / np.sqrt(self.bstar)
def wall_omega(self, y):
dy1 = y[1] - y[0]
return 60.0 / (self.b1 * dy1 ** 2) # Menter's wall value (10 x Wilcox's 6/(b1 y1^2))
def initial(self, y, dudy, nut0):
k = np.maximum(nut0 * np.abs(dudy) / np.sqrt(self.bstar), 1e-8)
om = np.maximum(np.abs(dudy) / np.sqrt(self.bstar), 1e-6)
d = np.maximum(wall_distance(y), 1e-12)
om = np.maximum(om, 6.0 / (self.b1 * d ** 2) * (d < 5)) # near-wall asymptote
k[[0, -1]] = 0.0
om[[0, -1]] = self.wall_omega(y)
return {"k": k, "omega": om}
def _blend(self, y, k, om):
d = np.maximum(wall_distance(y), 1e-12)
dk, dw = gradient(k, y), gradient(om, y)
cd = np.maximum(2.0 * self.sw2 / (om + EPS) * dk * dw, 1e-10)
arg1 = np.minimum(np.maximum(np.sqrt(np.maximum(k, 0)) / (self.bstar * om * d + EPS),
500.0 / (d ** 2 * om + EPS)),
4.0 * self.sw2 * k / (cd * d ** 2 + EPS))
F1 = np.tanh(arg1 ** 4)
arg2 = np.maximum(2.0 * np.sqrt(np.maximum(k, 0)) / (self.bstar * om * d + EPS),
500.0 / (d ** 2 * om + EPS))
F2 = np.tanh(arg2 ** 2)
return F1, F2, dk, dw
def eddy_viscosity(self, y, F, dudy):
k, om = np.maximum(F["k"], 0.0), np.maximum(F["omega"], EPS)
_, F2, _, _ = self._blend(y, k, om)
return self.a1 * k / np.maximum(self.a1 * om, np.abs(dudy) * F2)
def equations(self, y, F, u, dudy, nut):
k, om = np.maximum(F["k"], 0.0), np.maximum(F["omega"], EPS)
F1, F2, dk, dw = self._blend(y, k, om)
sk = F1 * self.sk1 + (1 - F1) * self.sk2
sw = F1 * self.sw1 + (1 - F1) * self.sw2
beta = F1 * self.b1 + (1 - F1) * self.b2
gam = F1 * self.g1 + (1 - F1) * self.g2
P = nut * dudy ** 2
Pk = np.minimum(P, 10.0 * self.bstar * k * om)
eq_k = (1.0 + sk * nut, Pk, -self.bstar * om, (0.0, 0.0))
cross = 2.0 * (1 - F1) * self.sw2 / om * dk * dw
S_C = gam * dudy ** 2 + np.maximum(cross, 0.0) + beta * om ** 2 # -beta w^2 Newton-linearised
S_P = -2.0 * beta * om + np.minimum(cross, 0.0) / om
wo = self.wall_omega(y)
eq_w = (1.0 + sw * nut, S_C, S_P, (wo, wo))
return {"k": eq_k, "omega": eq_w}
couette_rans/solver.py · 304 lines
"""Steady 1-D RANS solver for plane Couette flow in wall units.
Strategy
--------
* Momentum never needs a matrix: (1 + nu_t+) du+/dy+ = 1 => u+ = ∫ dy+/(1+nu_t+).
* Stage 1 (globalisation): each turbulence transport equation is advanced in
pseudo-time with a fully implicit (backward-Euler) step and Patankar source
linearisation, which keeps k, eps, omega, nu~ positive. The tridiagonal
systems go through LAPACK's gtsv (scipy.linalg.solve_banded) — the modern
name for the Thomas algorithm that was hand-written in 2000. The pseudo-time
step grows as res0/res (switched evolution relaxation).
* Stage 2 (convergence): Newton's method on the coupled discrete system with a
sparse finite-difference Jacobian, again under pseudo-transient continuation.
Quadratic convergence to machine precision in ~10 iterations.
The two-stage design is the point: segregated Picard iteration — what the 2000
k-epsilon program did — is a fine globaliser but can limit-cycle (it does for
Spalart–Allmaras here); Newton needs a decent starting point but then finishes
the job.
"""
from dataclasses import dataclass, field
import numpy as np
from scipy.integrate import cumulative_trapezoid
from scipy.linalg import solve_banded
from scipy.optimize import brentq
from scipy.sparse import lil_matrix, csc_matrix, identity
from scipy.sparse.linalg import spsolve
from .mesh import stretched_mesh
from .models import MixingLength
@dataclass
class CouetteResult:
model: str
retau: float
y: np.ndarray # y+ on the full gap [0, 2 Re_tau]
u: np.ndarray # u+
dudy: np.ndarray
nut: np.ndarray # nu_t / nu
fields: dict
iterations: int
residual_history: list = field(default_factory=list)
converged: bool = True
newton_iterations: int = 0
# --- derived quantities -------------------------------------------------
@property
def ue_plus(self):
"""Velocity of the moving wall in wall units (the 2000 report's U_e / u_tau)."""
return float(self.u[-1])
@property
def re_gap(self):
"""Re = U_e H / nu, the 2000 report's definition (moving-wall speed, full gap)."""
return self.ue_plus * 2.0 * self.retau
@property
def re_half(self):
"""Re = U_w h / nu as in Bech et al. (1995) and Lee & Moser (2018) (= re_gap / 4)."""
return self.re_gap / 4.0
@property
def cf(self):
"""Skin friction on the moving-wall speed U_e: Cf = 2 (u_tau/U_e)^2 (2000 convention)."""
return 2.0 / self.ue_plus ** 2
@property
def cf_half(self):
"""Skin friction on U_w = U_e/2 (DNS convention): Cf = 2 (u_tau/U_w)^2."""
return 8.0 / self.ue_plus ** 2
@property
def reynolds_shear_stress(self):
"""-u'v'+ = nu_t+ du+/dy+ (modelled)."""
return self.nut * self.dudy
@property
def asymmetry(self):
"""max |u(y) - (U_e - u(2Re_tau - y))| — should be at round-off level."""
return float(np.max(np.abs(self.u - (self.u[-1] - self.u[::-1]))))
def velocity_from_nut(y, nut):
dudy = 1.0 / (1.0 + nut)
u = cumulative_trapezoid(dudy, y, initial=0.0)
return u, dudy
# ---------------------------------------------------------------------------
# discrete operator
# ---------------------------------------------------------------------------
def _operator(y, gamma):
"""Coefficients of the conservative second-order diffusion operator on a non-uniform mesh."""
n = len(y)
gf = 0.5 * (gamma[1:] + gamma[:-1]) # diffusivity at the faces
dy = np.diff(y)
h = np.empty(n)
h[1:-1] = 0.5 * (y[2:] - y[:-2])
aW = np.zeros(n)
aE = np.zeros(n)
aW[1:-1] = gf[:-1] / (dy[:-1] * h[1:-1])
aE[1:-1] = gf[1:] / (dy[1:] * h[1:-1])
return aW, aE
def apply_operator(y, gamma, phi):
"""d/dy(gamma dphi/dy) at interior points (used by the verification tests)."""
aW, aE = _operator(y, gamma)
return aW[1:-1] * phi[:-2] - (aW[1:-1] + aE[1:-1]) * phi[1:-1] + aE[1:-1] * phi[2:]
def _equation_residual(y, phi, gamma, S_C, S_P):
return apply_operator(y, gamma, phi) + S_C[1:-1] + S_P[1:-1] * phi[1:-1]
def steady_residual(y, phi, gamma, S_C, S_P):
"""Max-norm residual of d/dy(gamma dphi/dy) + S_C + S_P phi = 0, normalised by max|S_C|."""
r = _equation_residual(y, phi, gamma, S_C, S_P)
return float(np.max(np.abs(r)) / (np.max(np.abs(S_C[1:-1])) + 1e-30))
def _implicit_step(y, phi, gamma, S_C, S_P, bc, dtau):
"""Backward-Euler step of dphi/dtau = d/dy(gamma dphi/dy) + S_C + S_P phi.
S_P <= 0 goes to the diagonal (Patankar), which keeps the matrix strictly
diagonally dominant and the update positive for S_C >= 0.
"""
n = len(y)
aW, aE = _operator(y, gamma)
diag = 1.0 / dtau + aW + aE - np.minimum(S_P, 0.0)
rhs = phi / dtau + S_C + np.maximum(S_P, 0.0) * phi
diag[0] = diag[-1] = 1.0
aW[-1] = aE[0] = 0.0
rhs[0], rhs[-1] = bc
ab = np.zeros((3, n))
ab[0, 1:] = -aE[:-1]
ab[1, :] = diag
ab[2, :-1] = -aW[1:]
return solve_banded((1, 1), ab, rhs)
# ---------------------------------------------------------------------------
# stage 1: segregated implicit pseudo-time (Picard)
# ---------------------------------------------------------------------------
def picard(model, y, F, tol=1e-6, max_iter=3000, dtau0=1.0, dtau_max=1e4, relax=0.5,
stall_window=150, verbose=False):
dudy0, nut0 = MixingLength().closed_form(y)
nut = model.eddy_viscosity(y, F, dudy0)
u, dudy = velocity_from_nut(y, nut)
hist, res0, dtau, converged, it = [], None, dtau0, False, 0
for it in range(1, max_iter + 1):
eqs = model.equations(y, F, u, dudy, nut)
res = max(steady_residual(y, F[nm], g, sc, sp) for nm, (g, sc, sp, _) in eqs.items())
hist.append(res)
if verbose and (it % 25 == 0 or it == 1):
print(f"picard {it:5d} dtau={dtau:9.2e} res={res:9.2e} Ue+={u[-1]:.5f}")
if res < tol:
converged = True
break
if it > stall_window and res > 0.5 * min(hist[:-stall_window]):
break # limit cycle: hand over to Newton
res0 = res0 or res
dtau = float(np.clip(dtau0 * res0 / res, dtau0, dtau_max))
for nm, (gamma, S_C, S_P, bc) in eqs.items():
F[nm] = np.maximum(_implicit_step(y, F[nm], gamma, S_C, S_P, bc, dtau), 0.0)
nut = relax * model.eddy_viscosity(y, F, dudy) + (1.0 - relax) * nut
u, dudy = velocity_from_nut(y, nut)
return F, hist, converged, it
# ---------------------------------------------------------------------------
# stage 2: Newton with pseudo-transient continuation
# ---------------------------------------------------------------------------
def _pack(F, names):
return np.concatenate([F[n][1:-1] for n in names])
def _unpack(x, names, template):
F, m = {}, len(template[names[0]]) - 2
for i, n in enumerate(names):
f = template[n].copy()
f[1:-1] = x[i * m:(i + 1) * m]
F[n] = f
return F
def _state(model, y, F):
"""Velocity and eddy viscosity consistent with the fields (nu_t may depend on du/dy through
SST's limiter, so iterate that weak dependence twice)."""
dudy, _ = MixingLength().closed_form(y)
nut = model.eddy_viscosity(y, F, dudy)
for _ in range(50):
u, dudy = velocity_from_nut(y, nut)
nut_new = model.eddy_viscosity(y, F, dudy)
done = np.max(np.abs(nut_new - nut)) < 1e-12 * (np.max(nut_new) + 1e-30)
nut = nut_new
if done:
break
u, dudy = velocity_from_nut(y, nut)
return u, dudy, nut
def _residual_vector(model, y, x, names, template):
F = _unpack(x, names, template)
for n in names:
F[n] = np.maximum(F[n], 0.0)
u, dudy, nut = _state(model, y, F)
eqs = model.equations(y, F, u, dudy, nut)
out, scale = [], 0.0
for n in names:
gamma, S_C, S_P, bc = eqs[n]
phi = F[n].copy()
phi[0], phi[-1] = bc
out.append(_equation_residual(y, phi, gamma, S_C, S_P))
scale = max(scale, float(np.max(np.abs(S_C[1:-1]))))
return np.concatenate(out), scale, u, dudy, nut, F
def _fd_jacobian(fun, x, r0, stencil=3, nfields=1):
"""Sparse finite-difference Jacobian by column colouring: a residual only depends on
unknowns within +-stencil points (gradients use +-1, second derivatives +-2), so every
(2 stencil + 1)-th column of each field can be perturbed together."""
m = len(x) // nfields
ncol = 2 * stencil + 1
J = lil_matrix((len(x), len(x)))
for f in range(nfields):
for c in range(ncol):
cols = np.arange(f * m + c, (f + 1) * m, ncol)
dx = np.zeros_like(x)
h = 1e-7 * np.maximum(np.abs(x[cols]), 1e-3)
dx[cols] = h
r1 = fun(x + dx)
for col, hc in zip(cols, h):
i0 = col - f * m
lo, hi = max(0, i0 - stencil), min(m, i0 + stencil + 1)
rows = np.concatenate([np.arange(ff * m + lo, ff * m + hi) for ff in range(nfields)])
J[rows, col] = (r1[rows] - r0[rows]) / hc
return csc_matrix(J)
def newton(model, y, F, tol=1e-10, max_iter=100, dtau0=1e2, dtau_max=1e12, verbose=False):
names = list(model.fields)
template = {k: v.copy() for k, v in F.items()}
u, dudy, nut = _state(model, y, F)
eqs = model.equations(y, F, u, dudy, nut)
for nm in names: # boundary values are fixed by the model
template[nm][0], template[nm][-1] = eqs[nm][3]
x = _pack(F, names)
fun = lambda xx: _residual_vector(model, y, xx, names, template)[0]
hist, res0, converged, it = [], None, False, 0
for it in range(1, max_iter + 1):
r, scale, u, dudy, nut, Fc = _residual_vector(model, y, x, names, template)
res = float(np.max(np.abs(r))) / (scale + 1e-30)
hist.append(res)
if verbose:
print(f"newton {it:4d} res={res:9.2e} Ue+={u[-1]:.6f}")
if res < tol:
converged = True
break
res0 = res0 or res
dtau = min(dtau0 * res0 / res, dtau_max)
J = _fd_jacobian(fun, x, r, stencil=3, nfields=len(names))
dx = spsolve(identity(len(x), format="csc") / dtau - J, r)
alpha, xn = 1.0, x + dx
while np.any(xn < 0) and alpha > 1e-3: # damping to keep positivity
alpha *= 0.5
xn = x + alpha * dx
x = np.maximum(xn, 0.0)
r, scale, u, dudy, nut, Fc = _residual_vector(model, y, x, names, template)
return Fc, u, dudy, nut, hist, converged, it
# ---------------------------------------------------------------------------
# public entry points
# ---------------------------------------------------------------------------
def solve(model, retau, n=300, beta=3.5, y=None, tol=1e-10, picard_tol=1e-6, verbose=False):
"""Solve plane Couette flow at friction Reynolds number retau (based on the half gap)."""
if y is None:
y = stretched_mesh(retau, n, beta)
if isinstance(model, MixingLength):
dudy, nut = model.closed_form(y)
u, _ = velocity_from_nut(y, nut)
return CouetteResult(model.name, retau, y, u, dudy, nut, {}, 0, [], True)
dudy0, nut0 = MixingLength().closed_form(y)
F = model.initial(y, dudy0, nut0)
F, hist1, conv1, it1 = picard(model, y, F, tol=picard_tol, verbose=verbose)
F, u, dudy, nut, hist2, conv2, it2 = newton(model, y, F, tol=tol, verbose=verbose)
return CouetteResult(model.name, retau, y, u, dudy, nut, F, it1 + it2, hist1 + hist2,
conv2, newton_iterations=it2)
def retau_from_re(model, re_gap, n=300, beta=3.5, **kw):
"""Find Re_tau such that the solution has Re = U_e H / nu = re_gap (the 2000 convention).
The laminar limit (Re = 4 Re_tau^2) brackets from below; Brent's method needs ~8 solves.
"""
def f(rt):
return solve(model, rt, n=n, beta=beta, **kw).re_gap - re_gap
lo = 0.5 * np.sqrt(re_gap)
hi = max(2.0 * lo, re_gap / 20.0)
while f(hi) < 0:
hi *= 1.5
rt = brentq(f, lo * 0.9, hi, xtol=1e-6, rtol=1e-8)
return solve(model, rt, n=n, beta=beta, **kw)
couette_rans/legacy.py · 81 lines
"""A faithful NumPy port of the March-2000 Fortran mixing-length program (Appendix D).
Kept deliberately close to the original so that its behaviour — Picard iteration
on the friction coefficient, uniform mesh, one-sided (backward) difference for
the velocity gradient feeding the eddy viscosity, hand-written Thomas algorithm,
convergence test |ΔCf| < 1e-8 and max|Δz| < 1e-8 — can be measured against the
2026 solver. Variable names follow the Fortran where it helps (Z, DZDY, TV, G).
"""
import time
import numpy as np
def thomas(a, b, c, r):
"""Tridiagonal solve, as in Numerical Recipes 'tridag' (the SOLVE subroutine of 2000)."""
n = len(b)
u = np.empty(n)
gam = np.empty(n)
bet = b[0]
u[0] = r[0] / bet
for j in range(1, n):
gam[j] = c[j - 1] / bet
bet = b[j] - a[j] * gam[j]
if bet == 0.0:
raise ZeroDivisionError("TRIDAG FAILED")
u[j] = (r[j] - a[j] * u[j - 1]) / bet
for j in range(n - 2, -1, -1):
u[j] -= gam[j + 1] * u[j + 1]
return u
def legacy_mixing_length_2000(re, n=10000, kappa=0.4, aplus=26.0, epsilon=1e-8,
max_iter=10_000_000, use_lapack=False):
"""Solve Couette flow at Re = U_e H / nu exactly the way the 2000 program did.
Returns a dict with y+ , u+ , Cf, gamma, iterations, wall time and the
symmetry error the 2000 report tabulated in its Table 2.
"""
from scipy.linalg import solve_banded
t0 = time.perf_counter()
dy = 1.0 / n
y = np.arange(n + 1) * dy # eta = y / H
z = y.copy() # initial guess: laminar (linear) profile
kk = kappa * kappa
zpast = np.full(n + 1, 2.0)
cf, cfpast = 2.0 / re, 0.1
it = 0
for it in range(1, max_iter + 1):
dzdy = np.empty(n + 1)
dzdy[0] = (z[1] - z[0]) / dy
dzdy[1:] = (z[1:] - z[:-1]) / dy # backward difference, as in 2000
dudywall = (z[n] - z[n - 1]) / dy # gradient at the moving wall
cf = 2.0 * dudywall / re
g = np.sqrt(abs(cf) / 2.0)
# TURBULENTV subroutine
tv = np.zeros(n + 1)
ylo = y[1:-1]
d = np.where(1.0 - ylo >= ylo, ylo, 1.0 - ylo)
tv[1:-1] = kk * d ** 2 * (1.0 - np.exp(-re * g * d / aplus)) ** 2 * dzdy[1:-1]
a = tv + 1.0 / re
am = np.roll(a, 1) # A(I-1)
ap = np.roll(a, -1) # A(I+1)
aa = 0.5 * (a + am)
cc = 0.5 * (ap + a)
bb = -aa - cc
aa[0] = cc[0] = 0.0; bb[0] = 1.0
aa[n] = cc[n] = 0.0; bb[n] = 1.0
r = np.zeros(n + 1); r[n] = 1.0
if use_lapack:
ab = np.zeros((3, n + 1)); ab[0, 1:] = cc[:-1]; ab[1] = bb; ab[2, :-1] = aa[1:]
z = solve_banded((1, 1), ab, r)
else:
z = thomas(aa, bb, cc, r)
if abs(cf - cfpast) < epsilon and np.max(np.abs(zpast[1:-1] - z[1:-1])) < epsilon:
break
zpast = z.copy()
cfpast = cf
yplus = y * g * re
uplus = z * np.sqrt(2.0 / abs(cf))
sym_err = float(np.max(np.abs(uplus - (uplus[-1] - uplus[::-1]))))
return dict(re=re, n=n, yplus=yplus, uplus=uplus, cf=cf, gamma=g, retau_gap=g * re,
iterations=it, seconds=time.perf_counter() - t0, symmetry_error=sym_err)
couette_rans/verify.py · 70 lines
"""Verification: the 2000 tests, redone, plus the ones that were missing.
1. thomas_test_2000 — the report's Table 1: z'' = (p+1)(p+2) eta^p, exact z = eta^(p+2).
2. mms_operator — method of manufactured solutions for the non-uniform-mesh diffusion
operator: error must fall as h^2 (observed order ~ 2).
3. grid_convergence — U_w+ for each model on successively finer stretched meshes, with a
Richardson estimate of the mesh-converged value.
"""
import numpy as np
from scipy.linalg import solve_banded
from .legacy import thomas
from .solver import apply_operator, solve
from .mesh import stretched_mesh
def thomas_test_2000(points=(10, 100, 1000, 10000), ps=(0, 1, 2, 3)):
out = {}
for n in points:
row = {}
for p in ps:
eta = np.linspace(0, 1, n + 1)
d = 1.0 / n
a = np.ones(n + 1); b = np.full(n + 1, -2.0); c = np.ones(n + 1)
r = d ** 2 * (p + 1) * (p + 2) * eta ** p
a[0] = c[0] = a[-1] = c[-1] = 0; b[0] = b[-1] = 1.0
r[0], r[-1] = 0.0, 1.0
z = thomas(a, b, c, r)
row[p] = float(np.max(np.abs(z - eta ** (p + 2))))
out[n] = row
return out
def mms_operator(ns=(25, 50, 100, 200, 400, 800), L=1000.0, beta=3.5):
"""Solve d/dy(gamma dphi/dy) = f with phi = sin(pi y/L)(1+y/L), gamma = 1 + (y/L)^2 on a
tanh-stretched mesh; the exact source f is known, so the discrete error is measurable."""
errs = {}
for n in ns:
y = stretched_mesh(L / 2, n, beta)
s = y / L
phi = np.sin(np.pi * s) * (1 + s)
dphi = (np.pi * np.cos(np.pi * s) * (1 + s) + np.sin(np.pi * s)) / L
d2phi = (-np.pi ** 2 * np.sin(np.pi * s) * (1 + s) + 2 * np.pi * np.cos(np.pi * s)) / L ** 2
gamma = 1 + s ** 2
dgamma = 2 * s / L
f = dgamma * dphi + gamma * d2phi
# assemble and solve the linear system directly
from .solver import _operator
aW, aE = _operator(y, gamma)
m = len(y)
ab = np.zeros((3, m))
diag = aW + aE; diag[0] = diag[-1] = 1.0
rhs = -f.copy(); rhs[0], rhs[-1] = phi[0], phi[-1]
# sign: aW phi_{i-1} - (aW+aE) phi_i + aE phi_{i+1} = f -> multiply by -1
ab[1] = diag
ab[0, 1:] = -aE[:-1]; ab[0, 1] = 0.0
ab[2, :-1] = -aW[1:]; ab[2, -2] = 0.0
sol = solve_banded((1, 1), ab, rhs)
errs[n] = float(np.max(np.abs(sol - phi)))
ns_ = sorted(errs)
orders = {ns_[i + 1]: float(np.log(errs[ns_[i]] / errs[ns_[i + 1]]) / np.log(2)) for i in range(len(ns_) - 1)}
return errs, orders
def grid_convergence(model, retau=500, ns=(50, 100, 200, 400, 800), beta=3.5):
vals = {n: solve(model, retau, n=n, beta=beta).ue_plus / 2 for n in ns}
ns_ = sorted(vals)
# Richardson with the two finest meshes assuming second order
f1, f2 = vals[ns_[-1]], vals[ns_[-2]]
extrap = f1 + (f1 - f2) / 3.0
return vals, extrap
run_all.py · 143 lines
"""Reproduce every number and figure of the 2026 rewrite. `python run_all.py` writes results/results.json.
Reynolds-number conventions (they matter — see the page):
re_gap = U_e H / nu : 2000 report (moving-wall speed, full gap)
re_half = U_w h / nu : Bech (1995), Lee & Moser (2018) (= re_gap / 4)
Re_tau = u_tau h / nu : friction Reynolds number on the half gap
"""
import json, time, glob, os
import numpy as np
from couette_rans import (MixingLength, LaunderSharmaKE, SpalartAllmaras, MenterSST,
solve, retau_from_re, legacy_mixing_length_2000)
from couette_rans.verify import thomas_test_2000, mms_operator, grid_convergence
HERE = os.path.dirname(os.path.abspath(__file__))
DNS_DIR = os.path.join(HERE, "data", "lee_moser_2018")
OUT = os.path.join(HERE, "results")
os.makedirs(OUT, exist_ok=True)
MODELS = [MixingLength(), LaunderSharmaKE(), SpalartAllmaras(), MenterSST()]
def load_dns():
cases = {}
for f in sorted(glob.glob(os.path.join(DNS_DIR, "LM_Couette_R*_100PI_mean_prof.dat"))):
tag = os.path.basename(f).split("_")[2]
hdr = open(f).read()
nu = float(hdr.split("nu =")[1].split()[0])
utau = float(hdr.split("u_tau =")[1].split()[0])
retau = float(hdr.split("Re_tau =")[1].split()[0])
d = np.loadtxt(f, comments="%")
fl = np.loadtxt(f.replace("mean_prof", "vel_fluc_prof"), comments="%")
cases[tag] = dict(retau=retau, nu=nu, utau=utau, re_half=1.0 / nu,
uw_plus=1.0 / utau, cf_half=2.0 * utau ** 2,
yplus=d[:, 1].tolist(), uplus=d[:, 2].tolist(), dudy=d[:, 3].tolist(),
uv=(-fl[:, 5]).tolist() if fl.shape[1] > 5 else None,
fl_yplus=fl[:, 1].tolist())
return cases
def subsample(y, u, k=2):
return [float(v) for v in y[::k]], [float(v) for v in u[::k]]
def main():
R = {"conventions": __doc__}
t_all = time.time()
# ---- DNS reference data ------------------------------------------------
dns = load_dns()
R["dns"] = dns
R["bech1995"] = dict(re_half=1300, retau=82.2, cf_half=2 * (82.2 / 1300) ** 2, uw_plus=1300 / 82.2)
# ---- modern solver at the DNS friction Reynolds numbers -------------------
R["profiles"] = {}
R["summary"] = []
for tag, c in dns.items():
R["profiles"][tag] = {}
for m in MODELS:
t = time.time()
r = solve(m, c["retau"])
dt = time.time() - t
y, u = subsample(r.y, r.u)
_, uv = subsample(r.y, r.reynolds_shear_stress)
_, nut = subsample(r.y, r.nut)
R["profiles"][tag][m.name] = dict(yplus=y, uplus=u, uv=uv, nut=nut,
residual_history=[float(v) for v in r.residual_history])
R["summary"].append(dict(case=tag, model=m.name, retau=c["retau"], uw_plus=r.ue_plus / 2,
dns_uw_plus=c["uw_plus"], cf_half=r.cf_half, dns_cf_half=c["cf_half"],
cf_error_pct=100 * (r.cf_half / c["cf_half"] - 1),
re_half=r.re_half, iterations=r.iterations,
newton_iterations=r.newton_iterations, seconds=dt,
asymmetry=r.asymmetry, converged=r.converged,
final_residual=float(r.residual_history[-1]) if r.residual_history else 0.0))
print(f"{tag} {m.name:28s} Uw+={r.ue_plus/2:7.3f} (DNS {c['uw_plus']:.3f}) it={r.iterations:4d} {dt:5.2f}s")
# ---- friction law: Cf(Re) for each model -----------------------------------
R["friction"] = {}
for m in MODELS:
pts = []
for rt in [40, 50, 60, 82.2, 93, 120, 150, 220, 300, 400, 500, 700, 1000, 1500, 2000]:
r = solve(m, rt)
pts.append(dict(retau=rt, re_half=r.re_half, cf_half=r.cf_half, uw_plus=r.ue_plus / 2))
R["friction"][m.name] = pts
# laminar: Cf_half = 2 (u_tau/U_w)^2 with u_tau^2 = nu U_w / h -> Cf = 2/Re_half
R["friction_laminar"] = [dict(re_half=re, cf_half=2.0 / re) for re in [100, 200, 300, 400]]
# ---- the 2000 program, ported and re-run --------------------------------------
R["legacy"] = []
for re in [3000, 5000, 7000, 10000, 15000, 20000]:
L = legacy_mixing_length_2000(re, n=10000, use_lapack=True)
y, u = subsample(L["yplus"], L["uplus"], 25)
R["legacy"].append(dict(re_gap=re, n=L["n"], iterations=L["iterations"], seconds_lapack=L["seconds"],
cf=L["cf"], gamma=L["gamma"], retau_gap=L["retau_gap"], retau_half=L["retau_gap"] / 2,
yplus1=float(L["yplus"][1]), symmetry_error=L["symmetry_error"],
yplus=y, uplus=u))
print(f"legacy Re={re}: {L['iterations']} it, Cf={L['cf']:.6f}, Re_tau(half)={L['retau_gap']/2:.1f}, sym={L['symmetry_error']:.2e}")
# pure-python Thomas timing for one case (what 2000 actually did, loop by loop)
L = legacy_mixing_length_2000(3000, n=10000, use_lapack=False)
R["legacy_pure_python_seconds_re3000"] = L["seconds"]
R["legacy_pure_python_iterations_re3000"] = L["iterations"]
# mesh dependence of the legacy scheme (first-order backward difference in nu_t)
R["legacy_mesh"] = [dict(n=n, **{k: v for k, v in legacy_mixing_length_2000(10000, n=n, use_lapack=True).items()
if k in ("cf", "iterations", "symmetry_error", "seconds")}) for n in (250, 500, 1000, 2000, 4000, 10000)]
# ---- same Reynolds numbers with the 2026 mixing-length code (Re -> Re_tau by Brent) ------
R["modern_at_legacy_re"] = []
for re in [3000, 5000, 7000, 10000, 15000, 20000]:
t = time.time()
r = retau_from_re(MixingLength(), re)
R["modern_at_legacy_re"].append(dict(re_gap=re, retau_half=r.retau, cf=r.cf, seconds=time.time() - t,
asymmetry=r.asymmetry))
# apples to apples: the 2026 code with the 2000 value kappa = 0.4 (the legacy runs above use 0.4 too)
R["same_kappa"] = []
for L in R["legacy"]:
r = retau_from_re(MixingLength(kappa=0.4), L["re_gap"], n=800)
R["same_kappa"].append(dict(re_gap=L["re_gap"], legacy_cf=L["cf"], modern_cf_k040=r.cf,
diff_pct=100 * (L["cf"] / r.cf - 1)))
# the report's "Re = 1300" comparison with Bech: what it is in each convention
r_lo = retau_from_re(MixingLength(), 1300)
r_hi = retau_from_re(MixingLength(), 5200)
R["re1300"] = dict(retau_half_if_gap_convention=r_lo.retau, uw_plus_if_gap=r_lo.ue_plus / 2,
retau_half_if_bech_convention=r_hi.retau, uw_plus_if_bech=r_hi.ue_plus / 2,
bech_retau=82.2, bech_uw_plus=1300 / 82.2)
R["laminar_transition_note"] = "Plane Couette flow is linearly stable at all Re; experiments (Tillmark & Alfredsson 1992) find sustained turbulence above Re_half ≈ 360, i.e. Re_gap ≈ 1440."
# ---- verification --------------------------------------------------------------
R["verify"] = {}
R["verify"]["thomas_2000"] = {str(k): {str(p): v for p, v in row.items()} for k, row in thomas_test_2000().items()}
errs, orders = mms_operator()
R["verify"]["mms"] = dict(errors={str(k): v for k, v in errs.items()}, orders={str(k): v for k, v in orders.items()})
R["verify"]["grid"] = {}
for m in MODELS:
vals, extrap = grid_convergence(m, retau=500, ns=(50, 100, 200, 400, 800))
R["verify"]["grid"][m.name] = dict(uw_plus={str(k): v for k, v in vals.items()}, richardson=extrap)
R["total_seconds"] = time.time() - t_all
with open(os.path.join(OUT, "results.json"), "w") as f:
json.dump(R, f, indent=1)
print(f"done in {R['total_seconds']:.1f} s -> results/results.json")
if __name__ == "__main__":
main()
tests/test_solver.py · 58 lines
"""pytest -q (the tests that the 2000 report did by hand, plus the ones it could not do)"""
import numpy as np
from couette_rans import (MixingLength, LaunderSharmaKE, SpalartAllmaras, MenterSST,
solve, retau_from_re, legacy_mixing_length_2000)
from couette_rans.legacy import thomas
from couette_rans.verify import thomas_test_2000, mms_operator
def test_thomas_matches_lapack():
rng = np.random.default_rng(1)
n = 50
a, b, c = rng.random(n), 4 + rng.random(n), rng.random(n)
r = rng.random(n)
from scipy.linalg import solve_banded
ab = np.zeros((3, n)); ab[0, 1:] = c[:-1]; ab[1] = b; ab[2, :-1] = a[1:]
assert np.allclose(thomas(a, b, c, r), solve_banded((1, 1), ab, r))
def test_thomas_reproduces_2000_table_1():
t = thomas_test_2000(points=(10, 100), ps=(2, 3))
assert abs(t[10][2] - 2.49e-3) < 2e-5 and abs(t[10][3] - 6.39e-3) < 2e-5
assert abs(t[100][2] - 2.49e-5) < 2e-7
def test_operator_is_second_order():
_, orders = mms_operator(ns=(50, 100, 200, 400))
assert all(abs(o - 2) < 0.1 for o in orders.values())
def test_mixing_length_closed_form_satisfies_constant_stress():
r = solve(MixingLength(), 300)
assert np.allclose((1 + r.nut) * r.dudy, 1.0)
assert r.asymmetry < 1e-10
def test_laminar_limit():
# with nu_t = 0 the profile is linear: U_e+ = 2 Re_tau
r = solve(MixingLength(kappa=0.0), 20)
assert abs(r.ue_plus - 40) < 1e-9
def test_transport_models_converge_and_are_symmetric():
for m in (LaunderSharmaKE(), SpalartAllmaras(), MenterSST()):
r = solve(m, 220)
assert r.converged and r.residual_history[-1] < 1e-9
assert r.asymmetry < 1e-7
assert np.all(r.nut >= 0)
def test_re_to_retau_roundtrip():
r = retau_from_re(MixingLength(), 10000)
assert abs(r.re_gap - 10000) < 1e-2
def test_legacy_port_agrees_with_modern_code_at_same_kappa():
L = legacy_mixing_length_2000(5000, n=4000, use_lapack=True)
r = retau_from_re(MixingLength(kappa=0.4), 5000, n=600)
assert abs(L["cf"] / r.cf - 1) < 2e-3
Data
- M. Lee & R. D. Moser, “Extreme-scale motions in turbulent plane Couette flows”, J. Fluid Mech. 842 (2018) 128–145, doi:10.1017/jfm.2018.131. Profiles from
turbulence.oden.utexas.edu/couette2018(100π boxes), included indata/. - K. H. Bech, N. Tillmark, P. H. Alfredsson & H. I. Andersson, “An investigation of turbulent plane Couette flow at low Reynolds numbers”, J. Fluid Mech. 286 (1995) 291–325. Used for Re = 1300, Reτ = 82.2 in the friction law.
Models and methods
- E. R. van Driest, “On turbulent flow near a wall”, J. Aero. Sci. 23 (1956) 1007–1011.
- B. E. Launder & B. I. Sharma, “Application of the energy-dissipation model of turbulence to the calculation of flow near a spinning disc”, Lett. Heat Mass Transfer 1 (1974) 131–138.
- P. R. Spalart & S. R. Allmaras, “A one-equation turbulence model for aerodynamic flows”, Rech. Aérosp. 1 (1994) 5–21; S. R. Allmaras, F. T. Johnson & P. R. Spalart, “Modifications and clarifications for the implementation of the Spalart–Allmaras turbulence model”, ICCFD7-1902 (2012).
- F. R. Menter, M. Kuntz & R. Langtry, “Ten years of industrial experience with the SST turbulence model”, Turbulence, Heat and Mass Transfer 4 (2003) 625–632.
- S. V. Patankar, Numerical Heat Transfer and Fluid Flow, Hemisphere (1980) — source-term linearisation.
- C. T. Kelley & D. E. Keyes, “Convergence analysis of pseudo-transient continuation”, SIAM J. Numer. Anal. 35 (1998) 508–523.
- P. J. Roache, “Code verification by the method of manufactured solutions”, J. Fluids Eng. 124 (2002) 4–10.
- S. Pirozzoli, M. Bernardini & P. Orlandi, “Turbulence statistics in Couette flow at high Reynolds number”, J. Fluid Mech. 758 (2014) 327–343.
- N. Tillmark & P. H. Alfredsson, “Experiments on transition in plane Couette flow”, J. Fluid Mech. 235 (1992) 89–102 — the Re ≈ 360 threshold.
The 2000 report's own references, as cited there
- W. P. Jones & B. E. Launder, Int. J. Heat Mass Transfer 15 (1972) 301–314; K.-Y. Chien, AIAA J. 20 (1982) 33–38; Y. Nagano & M. Tagawa, J. Fluids Eng. 112 (1990) 33–39; H. K. Myong & N. Kasagi, JSME Int. J. 33 (1990) 63–72; R. M. C. So, NASA CR 189608 (1991); D. C. Wilcox, Turbulence Modeling for CFD (1993); Press et al., Numerical Recipes in Fortran 77; E. M. Aydin & H. J. Leutheusser, ASME FED 51 (1987); K. Nakabayashi, O. Kitoh & F. Nishimura, experimental study of turbulent Couette flow at low Reynolds number.