shaia · shaia · Jul 9, 2026 · Jul 5, 2026 · Jul 7, 2026 · Jul 8, 2026
diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -173,6 +173,9 @@ if(BUILD_EXAMPLES)
     add_executable(lid_driven_cavity_direct examples/lid_driven_cavity_direct.c)
     target_link_libraries(lid_driven_cavity_direct PRIVATE CFD::Library)
 
+    add_executable(turbulent_channel examples/turbulent_channel.c)
+    target_link_libraries(turbulent_channel PRIVATE CFD::Library)
+
     # Poiseuille flow with stretched grids
     add_executable(poiseuille_stretched_grid examples/poiseuille_stretched_grid.c)
     target_link_libraries(poiseuille_stretched_grid PRIVATE CFD::Library)
@@ -309,6 +312,14 @@ if(BUILD_TESTS)
     # Energy equation tests
     add_executable(test_energy_solver tests/solvers/energy/test_energy_solver.c)
     add_executable(test_natural_convection tests/validation/test_natural_convection.c)
+
+    # Turbulence model tests
+    add_executable(test_turbulence_kepsilon tests/solvers/turbulence/test_turbulence_kepsilon.c)
+    add_executable(test_turbulence_sa tests/solvers/turbulence/test_turbulence_sa.c)
+    add_executable(test_turbulence_wall_functions tests/solvers/turbulence/test_turbulence_wall_functions.c)
+    add_executable(test_turbulence_disabled_regression tests/solvers/turbulence/test_turbulence_disabled_regression.c)
+    add_executable(test_turbulent_channel tests/validation/test_turbulent_channel.c)
+    add_executable(test_turbulence_cross_backend tests/solvers/turbulence/test_turbulence_cross_backend.c)
     # Release-only: also run the Ra=1e4 de Vahl Davis tier (finer grid, slower)
     if(CAVITY_FULL_VALIDATION)
         target_compile_definitions(test_natural_convection PRIVATE CAVITY_FULL_VALIDATION=1)
@@ -429,6 +440,12 @@ if(BUILD_TESTS)
     target_link_libraries(test_poiseuille_3d PRIVATE CFD::Library unity)
     target_link_libraries(test_energy_solver PRIVATE CFD::Library unity)
     target_link_libraries(test_natural_convection PRIVATE CFD::Library unity)
+    target_link_libraries(test_turbulence_kepsilon PRIVATE CFD::Library unity)
+    target_link_libraries(test_turbulence_sa PRIVATE CFD::Library unity)
+    target_link_libraries(test_turbulence_wall_functions PRIVATE CFD::Library unity)
+    target_link_libraries(test_turbulence_disabled_regression PRIVATE CFD::Library unity)
+    target_link_libraries(test_turbulent_channel PRIVATE CFD::Library unity)
+    target_link_libraries(test_turbulence_cross_backend PRIVATE CFD::Library unity)
     target_link_libraries(test_finite_differences PRIVATE unity $<$<NOT:$<PLATFORM_ID:Windows>>:m>)
     target_include_directories(test_finite_differences PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/lib/include)
     target_link_libraries(test_finite_differences_3d PRIVATE unity $<$<NOT:$<PLATFORM_ID:Windows>>:m>)
@@ -580,6 +597,16 @@ if(BUILD_TESTS)
     add_test(NAME EnergySolverTest COMMAND test_energy_solver)
     add_test(NAME NaturalConvectionTest COMMAND test_natural_convection)
 
+    # Turbulence model tests
+    add_test(NAME TurbulenceKEpsilonTest COMMAND test_turbulence_kepsilon)
+    add_test(NAME TurbulenceSpalartAllmarasTest COMMAND test_turbulence_sa)
+    add_test(NAME TurbulenceWallFunctionsTest COMMAND test_turbulence_wall_functions)
+    add_test(NAME TurbulenceDisabledRegressionTest COMMAND test_turbulence_disabled_regression)
+    add_test(NAME TurbulentChannelTest COMMAND test_turbulent_channel)
+    set_tests_properties(TurbulentChannelTest PROPERTIES LABELS "validation" TIMEOUT 600)
+    add_test(NAME TurbulenceCrossBackendTest COMMAND test_turbulence_cross_backend)
+    set_tests_properties(TurbulenceCrossBackendTest PROPERTIES LABELS "cross-arch")
+
     # Mathematical accuracy tests
     add_test(NAME FiniteDifferencesTest COMMAND test_finite_differences)
     add_test(NAME FiniteDifferences3DTest COMMAND test_finite_differences_3d)

diff --git a/README.md b/README.md
@@ -14,7 +14,8 @@ A production-grade computational fluid dynamics (CFD) library in C for solving 2
 - 🔧 **Pluggable Solvers**: Explicit Euler, RK2 (Heun), Projection Method (Chorin's algorithm)
 - 📊 **Linear Solvers**: Jacobi, SOR, Red-Black SOR, CG/PCG, BiCGSTAB
 - 🌡️ **Heat Transfer**: Energy equation (advection–diffusion) + Boussinesq buoyancy + thermal BCs
-- 🎯 **Validated**: Ghia lid-driven cavity, Taylor-Green vortex, Poiseuille flow, natural convection benchmarks
+- 🌀 **Turbulence (RANS)**: Standard k-ε and Spalart-Allmaras with log-law wall functions; validated vs turbulent channel flow at Re_τ = 395
+- 🎯 **Validated**: Ghia lid-driven cavity, Taylor-Green vortex, Poiseuille flow, natural convection, turbulent channel flow
 - 📈 **VTK/CSV Output**: Ready for ParaView, VisIt visualization
 - 💾 **Restart/Checkpoint**: Portable, versioned binary save/restore of complete simulation state
 - ⚡ **Performance**: SIMD-optimized with runtime CPU detection
@@ -153,6 +154,15 @@ int main(void) {
 | `rk4_omp` | OpenMP | Multi-threaded RK4 |
 | `rk4_gpu` | GPU | CUDA-accelerated RK4 |
 
+### Turbulence backend coverage
+
+| Feature | Scalar | AVX2 | NEON | OMP | GPU |
+| ------- | ------ | ---- | ---- | --- | --- |
+| k-ε + SA + wall functions | done | done | — | done | — |
+
+GPU NS solvers return `CFD_ERROR_UNSUPPORTED` when a turbulence model is enabled.
+Turbulence is limited to 2D uniform grids in this release.
+
 ## Project Structure
 
 ```text

diff --git a/ROADMAP.md b/ROADMAP.md
@@ -47,6 +47,7 @@ The single source of truth for backend gaps. Each algorithm targets scalar (CPU)
 |                     | RK2 (Heun)     | done | done     | —        | done     | done |
 |                     | RK4 (classical)| done | done     | —        | done     | done |
 | **Energy Eq.**      | Advec-diff + Boussinesq + thermal BCs | done | done | — | done | done |
+| **Turbulence**      | k-ε / SA + wall functions             | done | done | — | done | —    |
 | **Linear Solvers**  | Jacobi         | done | done     | done     | —        | done |
 |                     | SOR            | done | done     | done     | —        | done |
 |                     | Red-Black SOR  | done | done     | done     | done     | done |
@@ -84,7 +85,7 @@ Genuine constraints to be aware of (not backlog items):
 | Phase | Theme | Priority | Status — what remains |
 | ----- | ----- | -------- | --------------------- |
 | 1 | Core Solver Improvements | P0–P3 | GMRES, multigrid, implicit integrators, SIMPLE/PISO, nonlinear & eigenvalue solvers |
-| 2 | Physics Extensions | P1–P3 | Turbulence (RANS), compressible, species, multiphase; energy-eq. extensions |
+| 2 | Physics Extensions | P1–P3 | RANS k-ε/SA ✅; realizable k-ε, k-ω SST remain; compressible, species, multiphase; energy-eq. extensions |
 | 3 | Geometry & Mesh | P1–P2 | Unstructured meshes, mesh I/O, adaptive refinement (3D ✅) |
 | 4 | Scalability & Performance | P1–P2 | MPI, GPU improvements, profiling tools (modular libs ✅) |
 | 5 | I/O & Post-processing | P1–P3 | HDF5, modern VTK XML, in-situ viz (CSV ✅) |
@@ -226,13 +227,17 @@ source via host callback). See CHANGELOG.
 
 ### 2.2 Turbulence Models (P1)
 
-- [ ] Spalart-Allmaras (1-equation)
-- [ ] k-epsilon standard
+- [x] Spalart-Allmaras (1-equation)
+- [x] k-epsilon standard
 - [ ] k-epsilon realizable
 - [ ] k-omega SST
-- [ ] Wall functions
+- [x] Wall functions
 - [ ] Low-Reynolds number treatment
 
+**Done (2D, uniform grids):** standard k-ε and Spalart-Allmaras with log-law wall functions
+on scalar/OMP/AVX2 backends; validated against turbulent channel flow at Re_τ = 395
+(k-ε: u_τ error 2.9%; SA: 3.1%). GPU turbulence not yet implemented.
+
 ### 2.3 Compressible Flow (P2)
 
 - [ ] Density-based solver
@@ -507,7 +512,7 @@ easier to support new architectures.
 
 | Milestone | Target |
 | --------- | ------ |
-| **v0.4.0 — Turbulence** | At least one RANS model (k-ε or SA), wall functions, turbulent channel-flow validation |
+| ✅ **v0.4.0 — Turbulence** | *(complete: k-ε + SA, log-law wall functions, turbulent channel-flow validation at Re_τ = 395)* |
 | **v0.5.0 — Parallel Computing** | MPI parallelization, scalability benchmarks, HDF5 parallel I/O |
 | **v0.6.0 — Unstructured Meshes** | Unstructured mesh support, Gmsh import, complex-geometry examples |
 | **v1.0.0 — Production Ready** | All Phase 1–6 features, comprehensive validation, complete docs, stable API, performance optimized |

diff --git a/docs/guides/examples.md b/docs/guides/examples.md
@@ -828,6 +828,56 @@ Part 3: Grid Refinement (Explicit Euler, dt=5e-04, T=0.5)
 
 ---
 
+### 17. turbulent_channel.c
+
+**Purpose:** RANS turbulence model demonstration — turbulent channel flow at Re_τ = 395
+
+**What it demonstrates:**
+
+- Enabling k-ε or Spalart-Allmaras turbulence via `params.turb_model`
+- Configuring wall-function walls with `BC_TYPE_NOSLIP` faces in `params.turb_bc`
+- Calling `turbulence_init_uniform()` before time-stepping
+- Recovering the friction velocity u_τ via `turbulence_wall_u_tau()`
+- Comparing the computed u+ profile against the log law
+- VTK output with the four turbulence scalar fields (`turbulent_kinetic_energy`,
+  `dissipation_rate`, `nu_tilde`, `turbulent_viscosity`)
+
+**Problem setup:**
+
+- Domain: 4 × 2 (channel half-height δ = 1), flow is streamwise-uniform
+- Re_τ = 395, ν = 1/395, ρ = 1; constant body force f_x = u_τ²/δ = 1, so the
+  exact steady friction velocity is u_τ = 1
+- Grid: 16×21 uniform, first-node y+ ≈ 40 (wall-function window 30–100)
+- Bottom/top faces: `BC_TYPE_NOSLIP` (wall functions); left/right periodic
+- Direct solver interface (registry + `solver_step`), marched to a
+  kinetic-energy residual below 1e-6
+
+**Run:**
+```bash
+./turbulent_channel        # k-epsilon (default)
+./turbulent_channel ke     # k-epsilon explicitly
+./turbulent_channel sa     # Spalart-Allmaras
+```
+
+**Expected output (k-ε):**
+```
+  Converged at step 24183 (KE residual 9.98e-07)
+
+Recovered u_tau = 0.9712 (exact force balance: 1.0000)
+
+        y+        u+   log-law       err
+      39.5     14.10     14.17      0.5%
+      79.0     16.24     15.86      2.4%
+     118.5     17.26     16.85      2.5%
+     ...
+     395.0     20.10     19.78      1.6%
+
+Wrote turbulent_channel.vtk (open in ParaView to inspect nu_t/k).
+```
+The SA variant converges similarly (u_τ ≈ 0.969, u+ within ~4% of the log law).
+
+---
+
 ## Visualization
 
 ### VTK Files (ParaView/VisIt)

diff --git a/docs/reference/api-reference.md b/docs/reference/api-reference.md
@@ -260,11 +260,33 @@ typedef struct {
     double rho;             // Density
     int max_iterations;     // Maximum iterations
     double tolerance;       // Convergence tolerance
+    // Turbulence
+    turbulence_model_t        turb_model;  // TURB_MODEL_NONE (default), K_EPSILON, SPALART_ALLMARAS
+    ns_turbulence_bc_config_t turb_bc;    // Per-face turbulence BC types
 } ns_solver_params_t;
 
 ns_solver_params_t ns_solver_params_default(void);
 ```
 
+`turbulence_model_t` is defined in `cfd/solvers/navier_stokes_solver.h`:
+
+```c
+typedef enum {
+    TURB_MODEL_NONE             = 0,  // Laminar (default, zero overhead)
+    TURB_MODEL_K_EPSILON        = 1,  // Standard k-epsilon (Launder-Spalding)
+    TURB_MODEL_SPALART_ALLMARAS = 2,  // Spalart-Allmaras 1-equation (no-ft2)
+} turbulence_model_t;
+```
+
+`ns_turbulence_bc_config_t` holds one BC type per domain face.  Available types:
+
+| Constant | Meaning |
+|----------|---------|
+| `BC_TYPE_PERIODIC` | Periodic (default for all faces) |
+| `BC_TYPE_NEUMANN` | Zero-gradient (outlet) |
+| `BC_TYPE_DIRICHLET` | Fixed inlet values via per-face `k_values` / `eps_values` / `nu_tilde_values` |
+| `BC_TYPE_NOSLIP` | Log-law wall-function treatment |
+
 ### Solver Statistics
 
 ```c
@@ -273,6 +295,7 @@ typedef struct {
     double residual;        // Final residual
     double elapsed_time;    // Execution time (seconds)
     cfd_status_t status;    // Solver status
+    double max_nu_t;        // Maximum turbulent viscosity (0.0 when laminar)
 } ns_solver_stats_t;
 
 ns_solver_stats_t ns_solver_stats_default(void);
@@ -358,13 +381,18 @@ void flow_field_destroy(flow_field* field);
 
 ```c
 typedef struct {
-    size_t nx, ny, nz;  // Grid dimensions (nz=1 for 2D)
-    double* u;          // x-velocity [nx * ny * nz]
-    double* v;          // y-velocity [nx * ny * nz]
-    double* w;          // z-velocity [nx * ny * nz] (zero for 2D)
-    double* p;          // Pressure [nx * ny * nz]
-    double* rho;        // Density [nx * ny * nz]
-    double* T;          // Temperature [nx * ny * nz]
+    size_t nx, ny, nz;    // Grid dimensions (nz=1 for 2D)
+    double* u;            // x-velocity [nx * ny * nz]
+    double* v;            // y-velocity [nx * ny * nz]
+    double* w;            // z-velocity [nx * ny * nz] (zero for 2D)
+    double* p;            // Pressure [nx * ny * nz]
+    double* rho;          // Density [nx * ny * nz]
+    double* T;            // Temperature [nx * ny * nz]
+    // Turbulence fields (always allocated; zero when TURB_MODEL_NONE)
+    double* turb_k;       // Turbulent kinetic energy k [nx * ny * nz]
+    double* turb_eps;     // Dissipation rate epsilon   [nx * ny * nz]
+    double* turb_nu_tilde;// SA transported variable ν̃  [nx * ny * nz]
+    double* nu_t;         // Turbulent viscosity ν_t    [nx * ny * nz]
 } flow_field;
 ```
 
@@ -465,6 +493,97 @@ typedef struct {
 poisson_solver_stats_t poisson_solver_stats_default(void);
 ```
 
+## Turbulence API
+
+```c
+#include "cfd/solvers/turbulence_solver.h"
+```
+
+### turbulence_init_uniform
+
+```c
+cfd_status_t turbulence_init_uniform(flow_field* field,
+                                     const ns_solver_params_t* params,
+                                     double k0, double eps0, double nu_tilde0);
+```
+
+Initialize turbulence fields to spatially uniform values.  **Call this once before
+time-stepping whenever `params->turb_model != TURB_MODEL_NONE`.**
+
+- `k0` — initial turbulent kinetic energy (k-ε only; ignored for SA)
+- `eps0` — initial dissipation rate (k-ε only; ignored for SA)
+- `nu_tilde0` — initial SA transported variable (SA only; ignored for k-ε)
+
+Returns `CFD_SUCCESS`, or `CFD_ERROR_INVALID` if `field` or `params` is NULL.
+
+### turbulence_step_explicit
+
+```c
+cfd_status_t turbulence_step_explicit(flow_field* field,
+                                      const grid* g,
+                                      const ns_solver_params_t* params,
+                                      double dt, double time);
+```
+
+Advance the turbulence transport equations by one explicit time step.  This function is
+**called automatically** by all CPU/OMP/AVX2 NS solvers after the velocity and energy steps
+— users normally do not call it directly.
+
+Returns `CFD_SUCCESS` (no-op when `TURB_MODEL_NONE`), `CFD_ERROR_UNSUPPORTED` on a 3D grid
+or non-uniform grid spacing, or `CFD_ERROR_INVALID` for NULL arguments.
+
+### turbulence_apply_bcs
+
+```c
+cfd_status_t turbulence_apply_bcs(flow_field* field,
+                                  const grid* g,
+                                  const ns_solver_params_t* params);
+```
+
+Apply turbulence boundary conditions (periodic, Neumann, Dirichlet, or wall-function) to all
+faces as configured in `params->turb_bc`.  Also called automatically by the NS solvers.
+
+### turbulence_wall_u_tau
+
+```c
+double turbulence_wall_u_tau(double u_p, double y_p, double nu);
+```
+
+Compute the friction velocity u_τ from the first-node parallel velocity `u_p`, wall-normal
+distance `y_p`, and kinematic viscosity `nu` using the log-law (`u+ = ln(y+)/κ + B`) via
+Newton iteration.  Below y+ = 11.63 the linear viscous-sublayer law is used instead.  Used
+internally by the wall-function BC; exposed for testing and post-processing.
+
+### VTK and CSV turbulence output
+
+When a turbulence model is active, `write_vtk_flow_field()` automatically includes four
+additional point-data scalars:
+
+- `turbulent_kinetic_energy` (k)
+- `dissipation_rate` (ε)
+- `nu_tilde` (ν̃)
+- `turbulent_viscosity` (ν_t)
+
+CSV centerline files written by `write_centerline_to_csv()` gain three additional columns:
+`turb_k`, `turb_eps`, `nu_t`.
+
+### Minimal usage example
+
+```c
+sim->params.mu = 1.0 / 395.0;
+sim->params.turb_model = TURB_MODEL_K_EPSILON;
+sim->params.turb_bc.bottom = BC_TYPE_NOSLIP;  /* wall-function walls */
+sim->params.turb_bc.top    = BC_TYPE_NOSLIP;  /* other faces stay PERIODIC */
+
+turbulence_init_uniform(sim->field, &sim->params, k0, eps0, 0.0);
+
+/* run_simulation_step advances k/eps/nu_t automatically */
+for (int step = 0; step < n_steps; step++) {
+    cfd_status_t st = run_simulation_step(sim);
+    if (st != CFD_SUCCESS) { /* handle */ break; }
+}
+```
+
 ## I/O API
 
 ### VTK Output