Wire BPHX channel pressure drop on both sides with selectable correlations.

Replace isobaric 4-port closures with SimplifiedChannel (default) and Martin1996 DP models so z_dp and UI dp_correlation actually affect the Newton solve.

Co-authored-by: Cursor <cursoragent@cursor.com>

crates/cli/examples/bphx_evaporator_condenser.json

@@ -33,6 +33,7 @@
           "target_subcooling_k": 5.0,
           "emergent_pressure": true,
           "correlation": "Longo2004",
+          "dp_correlation": "SimplifiedChannel",
           "ua": 2500.0
         },
         {
@@ -52,6 +53,7 @@
           "target_superheat_k": 5.0,
           "emergent_pressure": true,
           "correlation": "Longo2004",
+          "dp_correlation": "SimplifiedChannel",
           "ua": 2000.0
         },
         {

crates/cli/src/batch.rs

@@ -315,7 +315,7 @@ fn process_single_file(
                 initialization_diagnostics: None,
                 dof: None,
                 elapsed_ms: 0,
-                raw_state_vector: None,
+                raw_state_vector: None, solved_variables: Vec::new(),
             }
         }
     }
@@ -425,7 +425,7 @@ mod tests {
                 initialization_diagnostics: None,
                 dof: None,
                 elapsed_ms: 50,
-                raw_state_vector: None,
+                raw_state_vector: None, solved_variables: Vec::new(),
             },
             SimulationResult {
                 input: "test2.json".to_string(),
@@ -439,7 +439,7 @@ mod tests {
                 initialization_diagnostics: None,
                 dof: None,
                 elapsed_ms: 0,
-                raw_state_vector: None,
+                raw_state_vector: None, solved_variables: Vec::new(),
             },
         ];
 
@@ -483,7 +483,7 @@ mod tests {
                 initialization_diagnostics: None,
                 dof: None,
                 elapsed_ms: 50,
-                raw_state_vector: None,
+                raw_state_vector: None, solved_variables: Vec::new(),
             },
             SimulationResult {
                 input: "test2.json".to_string(),
@@ -497,7 +497,7 @@ mod tests {
                 initialization_diagnostics: None,
                 dof: None,
                 elapsed_ms: 0,
-                raw_state_vector: None,
+                raw_state_vector: None, solved_variables: Vec::new(),
             },
         ];
 
@@ -529,7 +529,7 @@ mod tests {
                 initialization_diagnostics: None,
                 dof: None,
                 elapsed_ms: 50,
-                raw_state_vector: None,
+                raw_state_vector: None, solved_variables: Vec::new(),
             },
             SimulationResult {
                 input: "test2.json".to_string(),
@@ -543,7 +543,7 @@ mod tests {
                 initialization_diagnostics: None,
                 dof: None,
                 elapsed_ms: 0,
-                raw_state_vector: None,
+                raw_state_vector: None, solved_variables: Vec::new(),
             },
         ];
 
@@ -571,7 +571,7 @@ mod tests {
                 initialization_diagnostics: None,
                 dof: None,
                 elapsed_ms: 50,
-                raw_state_vector: None,
+                raw_state_vector: None, solved_variables: Vec::new(),
             },
             SimulationResult {
                 input: "test2.json".to_string(),
@@ -585,7 +585,7 @@ mod tests {
                 initialization_diagnostics: None,
                 dof: None,
                 elapsed_ms: 1000,
-                raw_state_vector: None,
+                raw_state_vector: None, solved_variables: Vec::new(),
             },
         ];
 
@@ -643,7 +643,7 @@ mod tests {
                 initialization_diagnostics: None,
                 dof: None,
                 elapsed_ms: 100,
-                raw_state_vector: None,
+                raw_state_vector: None, solved_variables: Vec::new(),
             }],
         };
 
@@ -667,7 +667,7 @@ mod tests {
                 initialization_diagnostics: None,
                 dof: None,
                 elapsed_ms: 50,
-                raw_state_vector: None,
+                raw_state_vector: None, solved_variables: Vec::new(),
             },
             SimulationResult {
                 input: "ok2.json".to_string(),
@@ -681,7 +681,7 @@ mod tests {
                 initialization_diagnostics: None,
                 dof: None,
                 elapsed_ms: 60,
-                raw_state_vector: None,
+                raw_state_vector: None, solved_variables: Vec::new(),
             },
             SimulationResult {
                 input: "fail.json".to_string(),
@@ -695,7 +695,7 @@ mod tests {
                 initialization_diagnostics: None,
                 dof: None,
                 elapsed_ms: 0,
-                raw_state_vector: None,
+                raw_state_vector: None, solved_variables: Vec::new(),
             },
         ];
 

crates/cli/src/run.rs

@@ -76,6 +76,13 @@ pub struct SimulationResult {
     /// solve via the `ENTROPYK_INIT_STATE_JSON` environment variable).
     #[serde(default, skip_serializing_if = "Option::is_none")]
     pub raw_state_vector: Option<Vec<f64>>,
+    /// Named solver-computed unknowns (free actuators + calibration factors
+    /// like `z_ua`, `opening`, `z_dp`) read out of `raw_state_vector` and
+    /// labeled with their stable id, owning component, and `[min, max]`
+    /// bounds. Empty for systems with no bounded control variables; skipped
+    /// on serialization when empty so older consumers keep working.
+    #[serde(default, skip_serializing_if = "Vec::is_empty")]
+    pub solved_variables: Vec<entropyk::SolvedVariable>,
 }
 
 /// Compact DoF summary for CLI / web API consumers.
@@ -365,7 +372,7 @@ fn execute_simulation(
             initialization_diagnostics: None,
             dof: None,
             elapsed_ms,
-            raw_state_vector: None,
+            raw_state_vector: None, solved_variables: Vec::new(),
         };
     }
 
@@ -400,7 +407,7 @@ fn execute_simulation(
                 initialization_diagnostics: None,
                 dof: None,
                 elapsed_ms,
-                raw_state_vector: None,
+                raw_state_vector: None, solved_variables: Vec::new(),
             };
         }
     };
@@ -710,7 +717,7 @@ fn execute_simulation(
                             initialization_diagnostics: None,
                             dof: None,
                             elapsed_ms,
-                            raw_state_vector: None,
+                            raw_state_vector: None, solved_variables: Vec::new(),
                         };
                     }
                 },
@@ -730,7 +737,7 @@ fn execute_simulation(
                         initialization_diagnostics: None,
                         dof: None,
                         elapsed_ms,
-                        raw_state_vector: None,
+                        raw_state_vector: None, solved_variables: Vec::new(),
                     };
                 }
             }
@@ -775,7 +782,7 @@ fn execute_simulation(
                                     initialization_diagnostics: None,
                                     dof: None,
                                     elapsed_ms,
-                                    raw_state_vector: None,
+                                    raw_state_vector: None, solved_variables: Vec::new(),
                                 };
                             }
                         },
@@ -795,7 +802,7 @@ fn execute_simulation(
                                 initialization_diagnostics: None,
                                 dof: None,
                                 elapsed_ms,
-                raw_state_vector: None,
+                raw_state_vector: None, solved_variables: Vec::new(),
             };
                         }
                         None => 1,
@@ -816,7 +823,7 @@ fn execute_simulation(
                                     initialization_diagnostics: None,
                                     dof: None,
                                     elapsed_ms,
-                                    raw_state_vector: None,
+                                    raw_state_vector: None, solved_variables: Vec::new(),
                                 };
                             }
                         },
@@ -836,7 +843,7 @@ fn execute_simulation(
                                 initialization_diagnostics: None,
                                 dof: None,
                                 elapsed_ms,
-                raw_state_vector: None,
+                raw_state_vector: None, solved_variables: Vec::new(),
             };
                         }
                         None => 0,
@@ -860,7 +867,7 @@ fn execute_simulation(
                             initialization_diagnostics: None,
                             dof: None,
                             elapsed_ms,
-                            raw_state_vector: None,
+                            raw_state_vector: None, solved_variables: Vec::new(),
                         };
                     }
                 }
@@ -880,7 +887,7 @@ fn execute_simulation(
                         initialization_diagnostics: None,
                         dof: None,
                         elapsed_ms,
-                        raw_state_vector: None,
+                        raw_state_vector: None, solved_variables: Vec::new(),
                     };
                 }
             }
@@ -950,7 +957,7 @@ fn execute_simulation(
                 initialization_diagnostics: None,
                 dof: None,
                 elapsed_ms,
-                raw_state_vector: None,
+                raw_state_vector: None, solved_variables: Vec::new(),
             };
         }
     }
@@ -976,7 +983,7 @@ fn execute_simulation(
             initialization_diagnostics: None,
             dof: None,
             elapsed_ms,
-            raw_state_vector: None,
+            raw_state_vector: None, solved_variables: Vec::new(),
         };
         let find_component = |name: &str| {
             config
@@ -1054,7 +1061,7 @@ fn execute_simulation(
                         initialization_diagnostics: None,
                         dof: None,
                         elapsed_ms,
-                        raw_state_vector: None,
+                        raw_state_vector: None, solved_variables: Vec::new(),
                     };
                 }
                 system.add_saturated_controller(controller);
@@ -1072,7 +1079,7 @@ fn execute_simulation(
                     initialization_diagnostics: None,
                     dof: None,
                     elapsed_ms,
-                    raw_state_vector: None,
+                    raw_state_vector: None, solved_variables: Vec::new(),
                 };
             }
         }
@@ -1246,7 +1253,7 @@ fn execute_simulation(
                                 initialization_diagnostics: None,
                                 dof: None,
                                 elapsed_ms,
-                                raw_state_vector: None,
+                                raw_state_vector: None, solved_variables: Vec::new(),
                             };
                         }
                         system.add_free_actuator(var_id);
@@ -1268,7 +1275,7 @@ fn execute_simulation(
                             initialization_diagnostics: None,
                             dof: None,
                             elapsed_ms,
-                            raw_state_vector: None,
+                            raw_state_vector: None, solved_variables: Vec::new(),
                         };
                     }
                 }
@@ -1289,7 +1296,7 @@ fn execute_simulation(
             initialization_diagnostics: None,
             dof: None,
             elapsed_ms,
-            raw_state_vector: None,
+            raw_state_vector: None, solved_variables: Vec::new(),
         };
     }
 
@@ -1320,7 +1327,7 @@ fn execute_simulation(
                 initialization_diagnostics: None,
                 dof: Some(summary),
                 elapsed_ms,
-                raw_state_vector: None,
+                raw_state_vector: None, solved_variables: Vec::new(),
             };
         }
         summary
@@ -1517,7 +1524,7 @@ fn execute_simulation(
             initialization_diagnostics: None,
             dof: Some(dof_summary),
             elapsed_ms,
-            raw_state_vector: None,
+            raw_state_vector: None, solved_variables: Vec::new(),
         };
     }
 
@@ -1746,6 +1753,7 @@ fn execute_simulation(
                 dof: Some(dof_summary.clone()),
                 elapsed_ms,
                 raw_state_vector: Some(converged.state.clone()),
+                solved_variables: entropyk::extract_solved_variables(&system, &converged.state),
             }
         }
         Err(e) => {
@@ -1791,7 +1799,7 @@ fn execute_simulation(
                 initialization_diagnostics: Some(initialization_diagnostics),
                 dof: Some(dof_summary),
                 elapsed_ms,
-                raw_state_vector: None,
+                raw_state_vector: None, solved_variables: Vec::new(),
             }
         }
     }
@@ -2939,6 +2947,32 @@ fn bphx_correlation_from_params(
     Ok(Some(correlation))
 }
 
+/// Parse BPHX frictional `dp_correlation` / `DpCorrelation` (independent of HTC).
+///
+/// Default when absent: `SimplifiedChannel`. Unknown ids are rejected (no silent isobaric).
+fn bphx_dp_correlation_from_params(
+    params: &std::collections::HashMap<String, serde_json::Value>,
+    component_name: &str,
+) -> CliResult<entropyk_components::heat_exchanger::BphxDpCorrelation> {
+    use entropyk_components::heat_exchanger::BphxDpCorrelation;
+
+    let Some(value) = params
+        .get("dp_correlation")
+        .or_else(|| params.get("DpCorrelation"))
+        .or_else(|| params.get("dpCorrelation"))
+    else {
+        return Ok(BphxDpCorrelation::default());
+    };
+    let Some(raw) = value.as_str() else {
+        return Err(CliError::Config(format!(
+            "{component_name}: dp_correlation must be a string (got {value})"
+        )));
+    };
+    BphxDpCorrelation::parse(raw).map_err(|msg| {
+        CliError::Config(format!("{component_name}: {msg}"))
+    })
+}
+
 /// Extract calibration Z-factors for BphxEvaporator/BphxCondenser from JSON params.
 ///
 /// Errors if `ua_nominal == 0` and an explicit `ua` override is provided (geometry is
@@ -5095,7 +5129,11 @@ fn create_component(
                 .with_target_superheat(target_superheat)
                 .with_refrigerant(refrigerant)
                 .with_secondary_fluid(secondary_fluid)
-                .with_fluid_backend(Arc::clone(&backend));
+                .with_fluid_backend(Arc::clone(&backend))
+                .with_dp_correlation(bphx_dp_correlation_from_params(
+                    params,
+                    "BphxEvaporator",
+                )?);
             if let Some(correlation) =
                 bphx_correlation_from_params(params, "BphxEvaporator")?
             {
@@ -5151,7 +5189,11 @@ fn create_component(
                 .with_refrigerant(refrigerant)
                 .with_secondary_fluid(secondary_fluid)
                 .with_fluid_backend(Arc::clone(&backend))
-                .with_target_subcooling(target_subcooling);
+                .with_target_subcooling(target_subcooling)
+                .with_dp_correlation(bphx_dp_correlation_from_params(
+                    params,
+                    "BphxCondenser",
+                )?);
             if let Some(correlation) =
                 bphx_correlation_from_params(params, "BphxCondenser")?
             {
@@ -5630,7 +5672,7 @@ mod tests {
             initialization_diagnostics: None,
             dof: None,
             elapsed_ms: 50,
-            raw_state_vector: None,
+            raw_state_vector: None, solved_variables: Vec::new(),
         };
 
         let json = serde_json::to_string_pretty(&result).unwrap();
@@ -5787,6 +5829,56 @@ mod tests {
         );
     }
 
+    #[test]
+    fn test_bphx_dp_correlation_defaults_and_aliases() {
+        use entropyk_components::heat_exchanger::BphxDpCorrelation;
+        use serde_json::json;
+
+        assert_eq!(
+            bphx_dp_correlation_from_params(&std::collections::HashMap::new(), "BphxEvaporator")
+                .unwrap(),
+            BphxDpCorrelation::SimplifiedChannel
+        );
+        let params = std::collections::HashMap::from([(
+            "DpCorrelation".to_string(),
+            json!("Martin1996"),
+        )]);
+        assert_eq!(
+            bphx_dp_correlation_from_params(&params, "BphxCondenser").unwrap(),
+            BphxDpCorrelation::Martin1996
+        );
+    }
+
+    #[test]
+    fn test_bphx_dp_correlation_rejects_unknown() {
+        use serde_json::json;
+
+        let params = std::collections::HashMap::from([(
+            "dp_correlation".to_string(),
+            json!("Amalfi2016"),
+        )]);
+        let error = bphx_dp_correlation_from_params(&params, "BphxEvaporator").unwrap_err();
+        assert!(
+            error.to_string().contains("unsupported dp_correlation"),
+            "{error}"
+        );
+    }
+
+    #[test]
+    fn test_bphx_dp_correlation_rejects_non_string() {
+        use serde_json::json;
+
+        let params = std::collections::HashMap::from([(
+            "dp_correlation".to_string(),
+            json!(42),
+        )]);
+        let error = bphx_dp_correlation_from_params(&params, "BphxEvaporator").unwrap_err();
+        assert!(
+            error.to_string().contains("must be a string"),
+            "{error}"
+        );
+    }
+
     fn comp_from_json(v: serde_json::Value) -> crate::config::ComponentConfig {
         serde_json::from_value(v).expect("valid component JSON")
     }

crates/cli/tests/batch_execution.rs

@@ -92,7 +92,7 @@ fn test_simulation_result_statuses() {
             initialization_diagnostics: None,
             dof: None,
             elapsed_ms: 50,
-            raw_state_vector: None,
+            raw_state_vector: None, solved_variables: Vec::new(),
         },
         SimulationResult {
             input: "fail.json".to_string(),
@@ -106,7 +106,7 @@ fn test_simulation_result_statuses() {
             initialization_diagnostics: None,
             dof: None,
             elapsed_ms: 0,
-            raw_state_vector: None,
+            raw_state_vector: None, solved_variables: Vec::new(),
         },
         SimulationResult {
             input: "timeout.json".to_string(),
@@ -120,7 +120,7 @@ fn test_simulation_result_statuses() {
             initialization_diagnostics: None,
             dof: None,
             elapsed_ms: 1000,
-            raw_state_vector: None,
+            raw_state_vector: None, solved_variables: Vec::new(),
         },
     ];
 
@@ -163,7 +163,7 @@ fn test_batch_aggregator_csv_output() {
             initialization_diagnostics: None,
             dof: None,
             elapsed_ms: 150,
-            raw_state_vector: None,
+            raw_state_vector: None, solved_variables: Vec::new(),
         },
         SimulationResult {
             input: "scenario2.json".to_string(),
@@ -183,7 +183,7 @@ fn test_batch_aggregator_csv_output() {
             initialization_diagnostics: None,
             dof: None,
             elapsed_ms: 200,
-            raw_state_vector: None,
+            raw_state_vector: None, solved_variables: Vec::new(),
         },
         SimulationResult {
             input: "scenario3.json".to_string(),
@@ -197,7 +197,7 @@ fn test_batch_aggregator_csv_output() {
             initialization_diagnostics: None,
             dof: None,
             elapsed_ms: 0,
-            raw_state_vector: None,
+            raw_state_vector: None, solved_variables: Vec::new(),
         },
     ];
 
@@ -229,7 +229,7 @@ fn test_batch_aggregator_json_summary() {
             initialization_diagnostics: None,
             dof: None,
             elapsed_ms: 50,
-            raw_state_vector: None,
+            raw_state_vector: None, solved_variables: Vec::new(),
         },
         SimulationResult {
             input: "test2.json".to_string(),
@@ -243,7 +243,7 @@ fn test_batch_aggregator_json_summary() {
             initialization_diagnostics: None,
             dof: None,
             elapsed_ms: 75,
-            raw_state_vector: None,
+            raw_state_vector: None, solved_variables: Vec::new(),
         },
         SimulationResult {
             input: "test3.json".to_string(),
@@ -257,7 +257,7 @@ fn test_batch_aggregator_json_summary() {
             initialization_diagnostics: None,
             dof: None,
             elapsed_ms: 5000,
-            raw_state_vector: None,
+            raw_state_vector: None, solved_variables: Vec::new(),
         },
     ];
 
@@ -322,7 +322,7 @@ fn test_batch_summary_csv_with_convergence() {
         initialization_diagnostics: None,
         dof: None,
         elapsed_ms: 300,
-        raw_state_vector: None,
+        raw_state_vector: None, solved_variables: Vec::new(),
     }];
 
     let summary = BatchSummary {

crates/cli/tests/single_run.rs

@@ -29,7 +29,7 @@ fn test_simulation_result_serialization() {
         initialization_diagnostics: None,
         dof: None,
         elapsed_ms: 50,
-        raw_state_vector: None,
+        raw_state_vector: None, solved_variables: Vec::new(),
     };
 
     let json = serde_json::to_string_pretty(&result).unwrap();
@@ -70,7 +70,7 @@ fn test_error_result_serialization() {
         initialization_diagnostics: None,
         dof: None,
         elapsed_ms: 0,
-        raw_state_vector: None,
+        raw_state_vector: None, solved_variables: Vec::new(),
     };
 
     let json = serde_json::to_string(&result).unwrap();
@@ -96,7 +96,7 @@ fn test_error_result_serializes_failure_diagnostics() {
         initialization_diagnostics: None,
         dof: None,
         elapsed_ms: 0,
-        raw_state_vector: None,
+        raw_state_vector: None, solved_variables: Vec::new(),
     };
 
     let json = serde_json::to_string(&result).unwrap();
@@ -2303,7 +2303,7 @@ fn test_structural_failure_serializes_without_diagnostics() {
         initialization_diagnostics: None,
         dof: None,
         elapsed_ms: 0,
-        raw_state_vector: None,
+        raw_state_vector: None, solved_variables: Vec::new(),
     };
 
     let json = serde_json::to_string(&result).unwrap();
@@ -2348,7 +2348,7 @@ fn test_success_result_does_not_include_failure_diagnostics() {
         initialization_diagnostics: None,
         dof: None,
         elapsed_ms: 120,
-        raw_state_vector: None,
+        raw_state_vector: None, solved_variables: Vec::new(),
     };
 
     let json = serde_json::to_string(&result).unwrap();

crates/components/src/heat_exchanger/bphx_condenser.rs

@@ -24,6 +24,7 @@
 //! ```
 
 use super::bphx_correlation::{BphxCorrelation, CorrelationEvaluation};
+use super::bphx_dp::BphxDpCorrelation;
 use super::bphx_exchanger::BphxExchanger;
 use super::bphx_geometry::{BphxGeometry, BphxType};
 use super::correlation_registry::CorrelationSelectionError;
@@ -150,6 +151,17 @@ impl BphxCondenser {
         self
     }
 
+    /// Sets the channel frictional pressure-drop correlation (independent of HTC).
+    pub fn with_dp_correlation(mut self, correlation: BphxDpCorrelation) -> Self {
+        self.inner = self.inner.with_dp_correlation(correlation);
+        self
+    }
+
+    /// Returns the configured pressure-drop correlation.
+    pub fn dp_correlation(&self) -> BphxDpCorrelation {
+        self.inner.dp_correlation()
+    }
+
     /// Sets the target subcooling in Kelvin.
     ///
     /// # Panics

crates/components/src/heat_exchanger/bphx_correlation.rs

@@ -344,7 +344,11 @@ impl BphxCorrelation {
             CorrelationId::Cooper1984
             | CorrelationId::Mostinski1963
             | CorrelationId::Friedel1979
-            | CorrelationId::MullerSteinhagenHeck1986 => None,
+            | CorrelationId::MullerSteinhagenHeck1986
+            | CorrelationId::SimplifiedChannel
+            | CorrelationId::Martin1996
+            | CorrelationId::AmalfiThome2016
+            | CorrelationId::TaoInfanteFerreiraSurvey => None,
         }
     }
 

crates/components/src/heat_exchanger/bphx_dp.rs

@@ -0,0 +1,422 @@
+//! BPHX channel pressure-drop correlations (friction), separate from HTC.
+//!
+//! Pressure-drop form (Fanning-like `f`, matching existing SimplifiedChannel):
+//! `ΔP = z_dp · 2 · f · L · G² / (ρ · d_h)`.
+//!
+//! Martin (1996/1999) returns a Fanning friction factor from the reciprocal-
+//! square-root chevron model (same intermediate `f` as `fluids` /
+//! `ht.friction_plate_Martin_1999` before the ×4 Darcy conversion).
+
+use super::bphx_geometry::BphxGeometry;
+use crate::ComponentError;
+use std::f64::consts::PI;
+
+/// Reference dynamic viscosity [Pa·s] used for Re in channel DP models.
+/// Matches the historical SimplifiedChannel implementation (Newton-friendly constant).
+pub const BPHX_DP_MU_REF_PA_S: f64 = 0.0002;
+
+/// Selectable BPHX frictional pressure-drop correlation.
+#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
+pub enum BphxDpCorrelation {
+    /// Existing Blasius/laminar C¹ blend (default).
+    #[default]
+    SimplifiedChannel,
+    /// Holger Martin chevron-plate single-phase friction (1996 / 1999 form).
+    Martin1996,
+}
+
+impl BphxDpCorrelation {
+    /// Stable serialization name.
+    pub const fn as_str(self) -> &'static str {
+        match self {
+            Self::SimplifiedChannel => "SimplifiedChannel",
+            Self::Martin1996 => "Martin1996",
+        }
+    }
+
+    /// Parse CLI/UI identifiers. Rejects unknown strings (no silent fallback).
+    pub fn parse(name: &str) -> Result<Self, String> {
+        let key = name.trim().to_ascii_lowercase().replace(['-', '_', ' '], "");
+        match key.as_str() {
+            "simplifiedchannel" | "simplified" | "channel" | "default" => {
+                Ok(Self::SimplifiedChannel)
+            }
+            "martin1996" | "martin1999" | "martin" => Ok(Self::Martin1996),
+            _ => Err(format!(
+                "unsupported dp_correlation '{name}'. Use 'SimplifiedChannel' or 'Martin1996'."
+            )),
+        }
+    }
+}
+
+/// Evaluated channel pressure drop and ∂ΔP/∂G at fixed ρ and z_dp.
+#[derive(Debug, Clone, Copy)]
+pub struct BphxDpEval {
+    /// Pressure drop [Pa] (always ≥ 0 for the channel friction model).
+    pub delta_p_pa: f64,
+    /// ∂ΔP/∂G [Pa / (kg/(m²·s))].
+    pub d_delta_p_d_g: f64,
+}
+
+/// Fanning friction factor and ∂f/∂Re for the selected correlation.
+fn fanning_friction(
+    corr: BphxDpCorrelation,
+    re: f64,
+    chevron_angle_deg: f64,
+) -> Result<(f64, f64), ComponentError> {
+    match corr {
+        BphxDpCorrelation::SimplifiedChannel => Ok(simplified_channel_fanning(re)),
+        BphxDpCorrelation::Martin1996 => {
+            validate_martin_chevron(chevron_angle_deg)?;
+            Ok(martin1996_fanning(re, chevron_angle_deg))
+        }
+    }
+}
+
+fn validate_martin_chevron(angle_deg: f64) -> Result<(), ComponentError> {
+    if !angle_deg.is_finite() || !(0.0..90.0).contains(&angle_deg) {
+        return Err(ComponentError::InvalidState(format!(
+            "Martin1996 requires chevron_angle in (0, 90) degrees (got {angle_deg})"
+        )));
+    }
+    Ok(())
+}
+
+/// SimplifiedChannel: laminar 64/Re, Blasius 0.079 Re^−0.25, C¹ blend on [2300, 4000].
+///
+/// Note: laminar uses the Darcy 64/Re coefficient inside the Fanning ΔP form —
+/// historical Entropyk convention; kept identical for continuity.
+fn simplified_channel_fanning(re: f64) -> (f64, f64) {
+    const RE_LAMINAR: f64 = 2300.0;
+    const RE_TURBULENT: f64 = 4000.0;
+    let re_safe = re.max(1e-9);
+    let f_lam = 64.0 / re_safe;
+    let df_lam = -64.0 / (re_safe * re_safe);
+
+    if re <= RE_LAMINAR {
+        return (f_lam, df_lam);
+    }
+
+    let f_turb = 0.079 * re_safe.powf(-0.25);
+    let df_turb = 0.079 * (-0.25) * re_safe.powf(-1.25);
+
+    if re >= RE_TURBULENT {
+        return (f_turb, df_turb);
+    }
+
+    let f = entropyk_core::smoothing::cubic_blend(
+        f_lam,
+        f_turb,
+        re,
+        RE_LAMINAR,
+        RE_TURBULENT,
+    );
+    // d/dRe of blend(f_lam(Re), f_turb(Re), Re): product rule on cubic Hermite.
+    let df_blend_dx = entropyk_core::smoothing::cubic_blend_derivative(
+        f_lam,
+        f_turb,
+        re,
+        RE_LAMINAR,
+        RE_TURBULENT,
+    );
+    // cubic_blend treats endpoints as constants w.r.t x; add chain-rule terms for
+    // f_lam(Re) and f_turb(Re) via the complementary weights.
+    let t = ((re - RE_LAMINAR) / (RE_TURBULENT - RE_LAMINAR)).clamp(0.0, 1.0);
+    let w_b = t * t * (3.0 - 2.0 * t); // weight on b (= f_turb)
+    let w_a = 1.0 - w_b;
+    let df = df_blend_dx + w_a * df_lam + w_b * df_turb;
+    (f, df)
+}
+
+/// Martin 1996/1999 Fanning friction factor vs Re and chevron angle φ [deg].
+///
+/// ```text
+/// 1/√f = cosφ / √(0.045 tanφ + 0.09 sinφ + f0/cosφ) + (1−cosφ) / √(3.8 f1)
+/// ```
+/// with piecewise `f0`, `f1` (Shah/Martin 1999 appendix form).
+/// Laminar/turbulent `f0`/`f1` branches are C¹-blended on Re ∈ [1800, 2200]
+/// (same Newton rationale as SimplifiedChannel’s [2300, 4000] blend).
+fn martin1996_fanning(re: f64, chevron_angle_deg: f64) -> (f64, f64) {
+    const RE_LO: f64 = 1800.0;
+    const RE_HI: f64 = 2200.0;
+    let re_safe = re.max(1e-9);
+
+    let lam = |re_v: f64| -> (f64, f64, f64, f64) {
+        let f0 = 16.0 / re_v;
+        let df0 = -16.0 / (re_v * re_v);
+        let f1 = 149.0 / re_v + 0.9625;
+        let df1 = -149.0 / (re_v * re_v);
+        (f0, df0, f1, df1)
+    };
+    let turb = |re_v: f64| -> (f64, f64, f64, f64) {
+        let u = 1.56 * re_v.ln() - 3.0;
+        let u_safe = if u.abs() < 1e-12 {
+            if u >= 0.0 {
+                1e-12
+            } else {
+                -1e-12
+            }
+        } else {
+            u
+        };
+        let f0 = u_safe.powi(-2);
+        let df0 = -2.0 * u_safe.powi(-3) * (1.56 / re_v);
+        let f1 = 9.75 * re_v.powf(-0.289);
+        let df1 = 9.75 * (-0.289) * re_v.powf(-1.289);
+        (f0, df0, f1, df1)
+    };
+
+    let (f0, df0, f1, df1) = if re_safe <= RE_LO {
+        lam(re_safe)
+    } else if re_safe >= RE_HI {
+        turb(re_safe)
+    } else {
+        let (f0_l, df0_l, f1_l, df1_l) = lam(re_safe);
+        let (f0_t, df0_t, f1_t, df1_t) = turb(re_safe);
+        let f0 = entropyk_core::smoothing::cubic_blend(f0_l, f0_t, re_safe, RE_LO, RE_HI);
+        let f1 = entropyk_core::smoothing::cubic_blend(f1_l, f1_t, re_safe, RE_LO, RE_HI);
+        let t = ((re_safe - RE_LO) / (RE_HI - RE_LO)).clamp(0.0, 1.0);
+        let w_b = t * t * (3.0 - 2.0 * t);
+        let w_a = 1.0 - w_b;
+        let df0_blend = entropyk_core::smoothing::cubic_blend_derivative(
+            f0_l, f0_t, re_safe, RE_LO, RE_HI,
+        );
+        let df1_blend = entropyk_core::smoothing::cubic_blend_derivative(
+            f1_l, f1_t, re_safe, RE_LO, RE_HI,
+        );
+        let df0 = df0_blend + w_a * df0_l + w_b * df0_t;
+        let df1 = df1_blend + w_a * df1_l + w_b * df1_t;
+        (f0, df0, f1, df1)
+    };
+
+    let phi = chevron_angle_deg * PI / 180.0;
+    let cos_p = phi.cos();
+    let sin_p = phi.sin();
+    let tan_p = phi.tan();
+    // Guard grazing angles so /cosφ stays finite inside (0, 90).
+    let cos_safe = cos_p.signum() * cos_p.abs().max(1e-12);
+
+    let a = 0.045 * tan_p + 0.09 * sin_p + f0 / cos_safe;
+    let b = 3.8 * f1.max(1e-30);
+    let a_safe = a.max(1e-30);
+
+    let term0 = cos_p / a_safe.sqrt();
+    let term1 = (1.0 - cos_p) / b.sqrt();
+    let rhs = term0 + term1;
+    let rhs_safe = rhs.max(1e-30);
+    let f = rhs_safe.powi(-2);
+
+    let da_dre = df0 / cos_safe;
+    let db_dre = 3.8 * df1;
+    let dterm0 = cos_p * (-0.5) * a_safe.powf(-1.5) * da_dre;
+    let dterm1 = (1.0 - cos_p) * (-0.5) * b.powf(-1.5) * db_dre;
+    let drhs = dterm0 + dterm1;
+    let df = -2.0 * rhs_safe.powi(-3) * drhs;
+    (f, df)
+}
+
+/// Channel frictional ΔP and ∂ΔP/∂G for the given correlation.
+///
+/// `mass_flux` may be signed; Re uses `|G|`, ΔP uses `G²` (existing convention).
+pub fn evaluate_channel_pressure_drop(
+    corr: BphxDpCorrelation,
+    geometry: &BphxGeometry,
+    mass_flux: f64,
+    rho: f64,
+    z_dp: f64,
+) -> Result<BphxDpEval, ComponentError> {
+    if !mass_flux.is_finite() {
+        return Ok(BphxDpEval {
+            delta_p_pa: 0.0,
+            d_delta_p_d_g: 0.0,
+        });
+    }
+    if !rho.is_finite() || rho < 1e-10 || geometry.dh < 1e-10 {
+        return Ok(BphxDpEval {
+            delta_p_pa: 0.0,
+            d_delta_p_d_g: 0.0,
+        });
+    }
+    if !geometry.plate_length.is_finite() || geometry.plate_length <= 0.0 {
+        return Ok(BphxDpEval {
+            delta_p_pa: 0.0,
+            d_delta_p_d_g: 0.0,
+        });
+    }
+    if !z_dp.is_finite() || z_dp < 0.0 {
+        return Err(ComponentError::InvalidState(format!(
+            "BPHX z_dp must be finite and >= 0 (got {z_dp})"
+        )));
+    }
+
+    let g_abs = mass_flux.abs();
+    let re = g_abs * geometry.dh / BPHX_DP_MU_REF_PA_S;
+    let (f, df_dre) = fanning_friction(corr, re, geometry.chevron_angle)?;
+
+    // ΔP = z_dp · 2 f L G² / (ρ dh)
+    let k = z_dp * 2.0 * geometry.plate_length / (rho * geometry.dh);
+    let g2 = mass_flux * mass_flux;
+    let delta_p_pa = k * f * g2;
+
+    // dRe/dG = sign(G) · dh/μ  (0 at stagnation)
+    let d_re_d_g = if g_abs < 1e-30 {
+        0.0
+    } else {
+        mass_flux.signum() * geometry.dh / BPHX_DP_MU_REF_PA_S
+    };
+    // d(G²)/dG = 2 G
+    let d_delta_p_d_g = k * (df_dre * d_re_d_g * g2 + f * 2.0 * mass_flux);
+
+    Ok(BphxDpEval {
+        delta_p_pa,
+        d_delta_p_d_g,
+    })
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::heat_exchanger::BphxGeometry;
+
+    fn geo() -> BphxGeometry {
+        BphxGeometry::from_dh_area(0.003, 0.5, 20)
+    }
+
+    #[test]
+    fn parse_accepts_aliases_rejects_unknown() {
+        assert_eq!(
+            BphxDpCorrelation::parse("SimplifiedChannel").unwrap(),
+            BphxDpCorrelation::SimplifiedChannel
+        );
+        assert_eq!(
+            BphxDpCorrelation::parse("martin-1996").unwrap(),
+            BphxDpCorrelation::Martin1996
+        );
+        assert!(BphxDpCorrelation::parse("Amalfi2016").is_err());
+        assert!(BphxDpCorrelation::parse("friedel").is_err());
+    }
+
+    #[test]
+    fn simplified_matches_historical_scale() {
+        let g = 30.0;
+        let rho = 1100.0;
+        let eval = evaluate_channel_pressure_drop(
+            BphxDpCorrelation::SimplifiedChannel,
+            &geo(),
+            g,
+            rho,
+            1.0,
+        )
+        .unwrap();
+        assert!(eval.delta_p_pa > 0.0);
+        assert!(eval.delta_p_pa.is_finite());
+    }
+
+    #[test]
+    fn martin_differs_from_simplified_same_state() {
+        let g = 80.0;
+        let rho = 1000.0;
+        let simp = evaluate_channel_pressure_drop(
+            BphxDpCorrelation::SimplifiedChannel,
+            &geo(),
+            g,
+            rho,
+            1.0,
+        )
+        .unwrap();
+        let martin = evaluate_channel_pressure_drop(
+            BphxDpCorrelation::Martin1996,
+            &geo(),
+            g,
+            rho,
+            1.0,
+        )
+        .unwrap();
+        assert!(
+            (simp.delta_p_pa - martin.delta_p_pa).abs() > 1.0,
+            "Martin and Simplified should differ: {} vs {}",
+            simp.delta_p_pa,
+            martin.delta_p_pa
+        );
+    }
+
+    #[test]
+    fn z_dp_scales_linearly() {
+        let g = 40.0;
+        let rho = 1100.0;
+        let a = evaluate_channel_pressure_drop(
+            BphxDpCorrelation::SimplifiedChannel,
+            &geo(),
+            g,
+            rho,
+            1.0,
+        )
+        .unwrap();
+        let b = evaluate_channel_pressure_drop(
+            BphxDpCorrelation::SimplifiedChannel,
+            &geo(),
+            g,
+            rho,
+            0.5,
+        )
+        .unwrap();
+        assert!((b.delta_p_pa - 0.5 * a.delta_p_pa).abs() < 1e-9 * a.delta_p_pa.max(1.0));
+    }
+
+    #[test]
+    fn zero_flow_smooth() {
+        let eval = evaluate_channel_pressure_drop(
+            BphxDpCorrelation::SimplifiedChannel,
+            &geo(),
+            0.0,
+            1100.0,
+            1.0,
+        )
+        .unwrap();
+        assert_eq!(eval.delta_p_pa, 0.0);
+        assert!(eval.d_delta_p_d_g.is_finite());
+    }
+
+    #[test]
+    fn martin_rejects_invalid_chevron() {
+        let mut g = geo();
+        g.chevron_angle = f64::NAN;
+        let err = evaluate_channel_pressure_drop(
+            BphxDpCorrelation::Martin1996,
+            &g,
+            30.0,
+            1000.0,
+            1.0,
+        )
+        .unwrap_err();
+        assert!(format!("{err}").contains("chevron"));
+    }
+
+    #[test]
+    fn analytic_d_dg_matches_central_difference() {
+        let geo = geo();
+        let g0 = 55.0;
+        let rho = 1050.0;
+        for corr in [
+            BphxDpCorrelation::SimplifiedChannel,
+            BphxDpCorrelation::Martin1996,
+        ] {
+            let eval = evaluate_channel_pressure_drop(corr, &geo, g0, rho, 1.0).unwrap();
+            let h = 1e-4;
+            let plus =
+                evaluate_channel_pressure_drop(corr, &geo, g0 + h, rho, 1.0).unwrap();
+            let minus =
+                evaluate_channel_pressure_drop(corr, &geo, g0 - h, rho, 1.0).unwrap();
+            let fd = (plus.delta_p_pa - minus.delta_p_pa) / (2.0 * h);
+            let rel = (eval.d_delta_p_d_g - fd).abs() / fd.abs().max(1.0);
+            assert!(
+                rel < 1e-4,
+                "{:?}: analytic {} vs fd {} (rel {})",
+                corr,
+                eval.d_delta_p_d_g,
+                fd,
+                rel
+            );
+        }
+    }
+}

crates/components/src/heat_exchanger/bphx_evaporator.rs

@@ -24,6 +24,7 @@
 //! ```
 
 use super::bphx_correlation::{BphxCorrelation, CorrelationEvaluation};
+use super::bphx_dp::BphxDpCorrelation;
 use super::bphx_exchanger::BphxExchanger;
 use super::bphx_geometry::{BphxGeometry, BphxType};
 use super::correlation_registry::CorrelationSelectionError;
@@ -154,6 +155,17 @@ impl BphxEvaporator {
         self
     }
 
+    /// Sets the channel frictional pressure-drop correlation (independent of HTC).
+    pub fn with_dp_correlation(mut self, correlation: BphxDpCorrelation) -> Self {
+        self.inner = self.inner.with_dp_correlation(correlation);
+        self
+    }
+
+    /// Returns the configured pressure-drop correlation.
+    pub fn dp_correlation(&self) -> BphxDpCorrelation {
+        self.inner.dp_correlation()
+    }
+
     /// Returns the component name.
     pub fn name(&self) -> &str {
         self.inner.name()

crates/components/src/heat_exchanger/bphx_exchanger.rs

@@ -33,6 +33,7 @@ use super::bphx_correlation::{
     BphxCorrelation, CorrelationEvaluation, CorrelationParams, CorrelationResult,
     CorrelationSelector, ValidityStatus,
 };
+use super::bphx_dp::{evaluate_channel_pressure_drop, BphxDpCorrelation, BphxDpEval};
 use super::bphx_geometry::{BphxGeometry, BphxType};
 use super::correlation_registry::{
     CorrelationSelectionError, ExchangerGeometryType, FlowRegime, SelectionOutcome,
@@ -55,6 +56,7 @@ pub struct BphxExchanger {
     inner: HeatExchanger<EpsNtuModel>,
     geometry: BphxGeometry,
     correlation_selector: CorrelationSelector,
+    dp_correlation: BphxDpCorrelation,
     refrigerant_id: String,
     secondary_fluid_id: String,
     fluid_backend: Option<Arc<dyn entropyk_fluids::FluidBackend>>,
@@ -70,6 +72,7 @@ impl std::fmt::Debug for BphxExchanger {
             .field("ua", &self.ua())
             .field("geometry", &self.geometry)
             .field("correlation", &self.correlation_selector.correlation)
+            .field("dp_correlation", &self.dp_correlation)
             .field("refrigerant_id", &self.refrigerant_id)
             .field("secondary_fluid_id", &self.secondary_fluid_id)
             .field("has_fluid_backend", &self.fluid_backend.is_some())
@@ -108,6 +111,7 @@ impl BphxExchanger {
             inner: HeatExchanger::new(model, "BphxExchanger"),
             geometry,
             correlation_selector: CorrelationSelector::default(),
+            dp_correlation: BphxDpCorrelation::default(),
             refrigerant_id: String::new(),
             secondary_fluid_id: String::new(),
             fluid_backend: None,
@@ -125,6 +129,7 @@ impl BphxExchanger {
             inner: HeatExchanger::new(model, "BphxExchanger"),
             geometry,
             correlation_selector: CorrelationSelector::default(),
+            dp_correlation: BphxDpCorrelation::default(),
             refrigerant_id: String::new(),
             secondary_fluid_id: String::new(),
             fluid_backend: None,
@@ -151,6 +156,17 @@ impl BphxExchanger {
         self
     }
 
+    /// Sets the channel frictional pressure-drop correlation (independent of HTC).
+    pub fn with_dp_correlation(mut self, correlation: BphxDpCorrelation) -> Self {
+        self.dp_correlation = correlation;
+        self
+    }
+
+    /// Returns the configured pressure-drop correlation.
+    pub fn dp_correlation(&self) -> BphxDpCorrelation {
+        self.dp_correlation
+    }
+
     /// Sets the refrigerant fluid identifier.
     pub fn with_refrigerant(mut self, fluid: impl Into<String>) -> Self {
         self.set_refrigerant_id(fluid);
@@ -338,12 +354,10 @@ impl BphxExchanger {
         self.last_validity_warning.set(false);
     }
 
-    /// Computes the pressure drop using a simplified correlation.
+    /// Computes the channel frictional pressure drop for the configured correlation.
     ///
-    /// ΔP = f_dp × (2 × f × L × G²) / (ρ × d_h)
-    ///
-    /// where f is the friction factor, L is the plate length, G is mass flux.
-    /// The result is scaled by `calib().z_dp`.
+    /// ΔP = z_dp × (2 × f × L × G²) / (ρ × d_h) with Fanning-like `f`
+    /// (SimplifiedChannel default; Martin1996 when selected).
     ///
     /// # Arguments
     ///
@@ -352,46 +366,158 @@ impl BphxExchanger {
     ///
     /// # Returns
     ///
-    /// Pressure drop in Pa, scaled by f_dp calibration factor.
+    /// Pressure drop in Pa, scaled by `calib().z_dp`. Returns 0 on invalid inputs
+    /// or (for Martin) invalid chevron — callers that need errors should use
+    /// [`evaluate_channel_pressure_drop`].
     pub fn compute_pressure_drop(&self, mass_flux: f64, rho: f64) -> f64 {
-        if rho < 1e-10 || self.geometry.dh < 1e-10 {
-            return 0.0;
+        evaluate_channel_pressure_drop(
+            self.dp_correlation,
+            &self.geometry,
+            mass_flux,
+            rho,
+            self.calib().z_dp,
+        )
+        .map(|e| e.delta_p_pa)
+        .unwrap_or(0.0)
+    }
+
+    /// Live `z_dp`: free actuator from state when indexed, else calib value.
+    fn live_z_dp(&self, state: &StateSlice) -> f64 {
+        self.inner
+            .calib_indices_ref()
+            .z_dp
+            .and_then(|idx| state.get(idx).copied())
+            .unwrap_or_else(|| self.calib().z_dp)
+    }
+
+    fn side_channel_dp(
+        &self,
+        mass_flow: f64,
+        rho: f64,
+        z_dp: f64,
+    ) -> Result<BphxDpEval, ComponentError> {
+        let g = self.geometry.mass_flux(mass_flow);
+        evaluate_channel_pressure_drop(self.dp_correlation, &self.geometry, g, rho, z_dp)
+    }
+
+    /// Overwrites the two 4-port pressure rows with `P_out − P_in + ΔP = 0`.
+    fn overwrite_pressure_closures_with_dp(
+        &self,
+        state: &StateSlice,
+        residuals: &mut ResidualVector,
+    ) -> Result<(), ComponentError> {
+        let Some(edges) = self.inner.four_port_edges() else {
+            return Ok(());
+        };
+        let row0 = self.inner.model_n_equations();
+        if residuals.len() < row0 + 2 {
+            return Err(ComponentError::InvalidResidualDimensions {
+                expected: row0 + 2,
+                actual: residuals.len(),
+            });
         }
 
-        let re = mass_flux.abs() * self.geometry.dh / 0.0002;
-        // C¹ laminar→turbulent blend over [2300, 4000] (same rationale as the
-        // pipe friction factor and the Gnielinski correlation): a hard switch at
-        // Re = 2300 put a kink in the pressure-drop residual that the analytic
-        // Newton Jacobian could not see, hurting convergence near the transition.
-        // Reynolds uses |mass_flux| so transient reverse flow during iteration
-        // yields a physical friction factor instead of the flat `64` the old
-        // signed/`max(1.0)` form produced. The lower clamp only guards the exact
-        // div-by-zero at stagnation; at 1e-9 it is far below any reachable Re, so
-        // it introduces no visible kink (the previous `max(1.0)` kinked at Re = 1).
-        const RE_LAMINAR: f64 = 2300.0;
-        const RE_TURBULENT: f64 = 4000.0;
-        let f_laminar = 64.0 / re.max(1e-9);
-        let f = if re <= RE_LAMINAR {
-            f_laminar
-        } else {
-            let f_turbulent = 0.079 * re.powf(-0.25);
-            if re >= RE_TURBULENT {
-                f_turbulent
-            } else {
-                entropyk_core::smoothing::cubic_blend(
-                    f_laminar,
-                    f_turbulent,
-                    re,
-                    RE_LAMINAR,
-                    RE_TURBULENT,
-                )
-            }
+        let z_dp = self.live_z_dp(state);
+        let (m_h, p_h_in, h_h_in) = edges.hot_in;
+        let (_, p_h_out, _) = edges.hot_out;
+        let (m_c, p_c_in, h_c_in) = edges.cold_in;
+        let (_, p_c_out, _) = edges.cold_out;
+
+        let max_idx = [m_h, p_h_in, h_h_in, p_h_out, m_c, p_c_in, h_c_in, p_c_out]
+            .into_iter()
+            .max()
+            .unwrap_or(0);
+        if max_idx >= state.len() {
+            return Err(ComponentError::InvalidStateDimensions {
+                expected: max_idx + 1,
+                actual: state.len(),
+            });
+        }
+
+        let rho_hot = self.inner.side_density(
+            "hot",
+            self.inner.hot_fluid_id_str(),
+            state[p_h_in],
+            state[h_h_in],
+        )?;
+        let rho_cold = self.inner.side_density(
+            "cold",
+            self.inner.cold_fluid_id_str(),
+            state[p_c_in],
+            state[h_c_in],
+        )?;
+
+        let dp_hot = self.side_channel_dp(state[m_h], rho_hot, z_dp)?;
+        let dp_cold = self.side_channel_dp(state[m_c], rho_cold, z_dp)?;
+
+        residuals[row0] = state[p_h_out] - state[p_h_in] + dp_hot.delta_p_pa;
+        residuals[row0 + 1] = state[p_c_out] - state[p_c_in] + dp_cold.delta_p_pa;
+        Ok(())
+    }
+
+    /// Adds analytic ∂r/∂ṁ and ∂r/∂z_dp for BPHX pressure rows (P entries already set).
+    fn append_pressure_dp_jacobian(
+        &self,
+        state: &StateSlice,
+        jacobian: &mut JacobianBuilder,
+    ) -> Result<(), ComponentError> {
+        let Some(edges) = self.inner.four_port_edges() else {
+            return Ok(());
         };
+        let row0 = self.inner.model_n_equations();
+        let z_dp = self.live_z_dp(state);
+        let (m_h, p_h_in, h_h_in) = edges.hot_in;
+        let (_, p_h_out, _) = edges.hot_out;
+        let (m_c, p_c_in, h_c_in) = edges.cold_in;
+        let (_, p_c_out, _) = edges.cold_out;
 
-        let dp_base =
-            2.0 * f * self.geometry.plate_length * mass_flux.powi(2) / (rho * self.geometry.dh);
+        let max_idx = [m_h, p_h_in, h_h_in, p_h_out, m_c, p_c_in, h_c_in, p_c_out]
+            .into_iter()
+            .max()
+            .unwrap_or(0);
+        if max_idx >= state.len() {
+            return Err(ComponentError::InvalidStateDimensions {
+                expected: max_idx + 1,
+                actual: state.len(),
+            });
+        }
 
-        dp_base * self.calib().z_dp
+        let rho_hot = self.inner.side_density(
+            "hot",
+            self.inner.hot_fluid_id_str(),
+            state[p_h_in],
+            state[h_h_in],
+        )?;
+        let rho_cold = self.inner.side_density(
+            "cold",
+            self.inner.cold_fluid_id_str(),
+            state[p_c_in],
+            state[h_c_in],
+        )?;
+
+        let dp_hot = self.side_channel_dp(state[m_h], rho_hot, z_dp)?;
+        let dp_cold = self.side_channel_dp(state[m_c], rho_cold, z_dp)?;
+
+        let a_flow = self.geometry.channel_flow_area() * self.geometry.n_channels() as f64;
+        if a_flow <= 1e-30 {
+            return Err(ComponentError::InvalidState(
+                "BPHX channel flow area too small for pressure-drop Jacobian".into(),
+            ));
+        }
+        let d_g_d_m = 1.0 / a_flow;
+
+        // r = P_out − P_in + ΔP(G(ṁ), z_dp)  ⇒  ∂r/∂ṁ = ∂ΔP/∂G · ∂G/∂ṁ
+        jacobian.add_entry(row0, m_h, dp_hot.d_delta_p_d_g * d_g_d_m);
+        jacobian.add_entry(row0 + 1, m_c, dp_cold.d_delta_p_d_g * d_g_d_m);
+
+        if let Some(z_idx) = self.inner.calib_indices_ref().z_dp {
+            // ΔP = z_dp · ΔP_corr  ⇒  ∂r/∂z_dp = ΔP_corr = ΔP / z_dp
+            // Keep the column even near z_dp→0 so a free actuator stays visible to Newton.
+            let z_safe = z_dp.max(1e-30);
+            jacobian.add_entry(row0, z_idx, dp_hot.delta_p_pa / z_safe);
+            jacobian.add_entry(row0 + 1, z_idx, dp_cold.delta_p_pa / z_safe);
+        }
+        Ok(())
     }
 
     /// Updates UA based on computed HTC.
@@ -413,7 +539,9 @@ impl Component for BphxExchanger {
         state: &StateSlice,
         residuals: &mut ResidualVector,
     ) -> Result<(), ComponentError> {
-        self.inner.compute_residuals(state, residuals)
+        self.inner.compute_residuals(state, residuals)?;
+        // Replace isobaric P_out − P_in with channel ΔP on both sides.
+        self.overwrite_pressure_closures_with_dp(state, residuals)
     }
 
     fn jacobian_entries(
@@ -421,7 +549,9 @@ impl Component for BphxExchanger {
         state: &StateSlice,
         jacobian: &mut JacobianBuilder,
     ) -> Result<(), ComponentError> {
-        self.inner.jacobian_entries(state, jacobian)
+        self.inner.jacobian_entries(state, jacobian)?;
+        // Inner already wrote ∂r/∂P_out=+1, ∂r/∂P_in=−1; add ṁ and z_dp terms.
+        self.append_pressure_dp_jacobian(state, jacobian)
     }
 
     fn get_ports(&self) -> &[ConnectedPort] {
@@ -844,4 +974,131 @@ mod tests {
             "ΔP slope kink near Re = 1 ({s_below:.3e} vs {s_above:.3e})"
         );
     }
+
+    #[test]
+    fn test_bphx_both_sides_pressure_residuals_include_dp() {
+        use entropyk_fluids::TestBackend;
+        use std::sync::Arc;
+
+        // State layout per port triple (ṁ, P, h):
+        // hot_in 0..3, hot_out 3..6, cold_in 6..9, cold_out 9..12
+        let mut hx = BphxExchanger::new(test_geometry())
+            .with_fluid_backend(Arc::new(TestBackend::new()));
+        hx.set_hot_fluid("Water");
+        hx.set_cold_fluid("Water");
+        hx.set_port_context(&[
+            Some((0, 1, 2)),
+            Some((3, 4, 5)),
+            Some((6, 7, 8)),
+            Some((9, 10, 11)),
+        ]);
+
+        let mut state = vec![0.0; 12];
+        // Liquid water-ish P-h (TestBackend): ~2 bar, ~80 kJ/kg
+        state[0] = 0.25; // m_hot
+        state[1] = 200_000.0;
+        state[2] = 80_000.0;
+        state[3] = 0.25;
+        state[4] = 200_000.0; // isobaric guess — residual should show +ΔP
+        state[5] = 70_000.0;
+        state[6] = 0.20; // m_cold
+        state[7] = 200_000.0;
+        state[8] = 50_000.0;
+        state[9] = 0.20;
+        state[10] = 200_000.0;
+        state[11] = 60_000.0;
+
+        let n = hx.n_equations();
+        assert!(n >= 4, "4-port BPHX should include two pressure rows");
+        let mut residuals = vec![0.0; n];
+        hx.compute_residuals(&state, &mut residuals).unwrap();
+
+        let row_p = hx.inner.model_n_equations();
+        assert!(
+            residuals[row_p].abs() > 1.0,
+            "hot pressure residual should include ΔP, got {}",
+            residuals[row_p]
+        );
+        assert!(
+            residuals[row_p + 1].abs() > 1.0,
+            "cold pressure residual should include ΔP, got {}",
+            residuals[row_p + 1]
+        );
+
+        // z_dp scale: half calib → half residual contribution at same P_out=P_in
+        let mut hx_half = BphxExchanger::new(test_geometry())
+            .with_fluid_backend(Arc::new(TestBackend::new()));
+        hx_half.set_hot_fluid("Water");
+        hx_half.set_cold_fluid("Water");
+        hx_half.set_port_context(&[
+            Some((0, 1, 2)),
+            Some((3, 4, 5)),
+            Some((6, 7, 8)),
+            Some((9, 10, 11)),
+        ]);
+        let mut calib = Calib::default();
+        calib.z_dp = 0.5;
+        hx_half.set_calib(calib);
+        let mut r_half = vec![0.0; n];
+        hx_half.compute_residuals(&state, &mut r_half).unwrap();
+        assert!(
+            (r_half[row_p] - 0.5 * residuals[row_p]).abs() < 1e-6 * residuals[row_p].abs().max(1.0),
+            "z_dp should scale hot ΔP residual"
+        );
+        assert!(
+            (r_half[row_p + 1] - 0.5 * residuals[row_p + 1]).abs()
+                < 1e-6 * residuals[row_p + 1].abs().max(1.0),
+            "z_dp should scale cold ΔP residual"
+        );
+    }
+
+    #[test]
+    fn test_bphx_martin_dp_in_residuals_differs() {
+        use entropyk_fluids::TestBackend;
+        use std::sync::Arc;
+
+        let ports = [
+            Some((0, 1, 2)),
+            Some((3, 4, 5)),
+            Some((6, 7, 8)),
+            Some((9, 10, 11)),
+        ];
+        let mut state = vec![0.0; 12];
+        state[0] = 0.30;
+        state[1] = 200_000.0;
+        state[2] = 80_000.0;
+        state[3] = 0.30;
+        state[4] = 200_000.0;
+        state[5] = 70_000.0;
+        state[6] = 0.30;
+        state[7] = 200_000.0;
+        state[8] = 50_000.0;
+        state[9] = 0.30;
+        state[10] = 200_000.0;
+        state[11] = 60_000.0;
+
+        let mut simp = BphxExchanger::new(test_geometry())
+            .with_fluid_backend(Arc::new(TestBackend::new()));
+        simp.set_hot_fluid("Water");
+        simp.set_cold_fluid("Water");
+        simp.set_port_context(&ports);
+
+        let mut martin = BphxExchanger::new(test_geometry())
+            .with_dp_correlation(BphxDpCorrelation::Martin1996)
+            .with_fluid_backend(Arc::new(TestBackend::new()));
+        martin.set_hot_fluid("Water");
+        martin.set_cold_fluid("Water");
+        martin.set_port_context(&ports);
+
+        let n = simp.n_equations();
+        let mut r_s = vec![0.0; n];
+        let mut r_m = vec![0.0; n];
+        simp.compute_residuals(&state, &mut r_s).unwrap();
+        martin.compute_residuals(&state, &mut r_m).unwrap();
+        let row = simp.inner.model_n_equations();
+        assert!(
+            (r_s[row] - r_m[row]).abs() > 1.0,
+            "Martin residual should differ from Simplified"
+        );
+    }
 }

crates/components/src/heat_exchanger/correlation_registry.rs

@@ -38,6 +38,14 @@ pub enum CorrelationId {
     Friedel1979,
     /// Müller-Steinhagen-Heck (1986) two-phase pressure-drop correlation.
     MullerSteinhagenHeck1986,
+    /// BPHX SimplifiedChannel single-phase channel friction (default ΔP).
+    SimplifiedChannel,
+    /// Martin (1996/1999) chevron-plate single-phase friction.
+    Martin1996,
+    /// Amalfi–Thome (2016) boiling plate ΔP — metadata only (not implemented).
+    AmalfiThome2016,
+    /// Tao & Infante Ferreira condensation-in-BPHE DP survey — metadata only.
+    TaoInfanteFerreiraSurvey,
 }
 
 impl CorrelationId {
@@ -58,6 +66,10 @@ impl CorrelationId {
             Self::Mostinski1963 => "mostinski-1963",
             Self::Friedel1979 => "friedel-1979",
             Self::MullerSteinhagenHeck1986 => "muller-steinhagen-heck-1986",
+            Self::SimplifiedChannel => "simplified-channel",
+            Self::Martin1996 => "martin-1996",
+            Self::AmalfiThome2016 => "amalfi-thome-2016",
+            Self::TaoInfanteFerreiraSurvey => "tao-infante-ferreira-survey",
         }
     }
 }
@@ -337,6 +349,28 @@ const POOL_BOILING_GEOMETRIES: &[ExchangerGeometryType] = &[
     ExchangerGeometryType::FinnedTube,
 ];
 const DP_TUBE_GEOMETRIES: &[ExchangerGeometryType] = &[ExchangerGeometryType::SmoothTube];
+const BPHX_DP_BOUNDS: &[BoundedQuantity] = &[
+    bound(
+        BoundedQuantityKind::Reynolds,
+        1.0,
+        100_000.0,
+        SINGLE_PHASE_REGIMES,
+    ),
+    bound(
+        BoundedQuantityKind::MassFlux,
+        1.0,
+        2000.0,
+        SINGLE_PHASE_REGIMES,
+    ),
+];
+const AMALFI_BOUNDS: &[BoundedQuantity] = &[
+    bound(BoundedQuantityKind::MassFlux, 5.0, 600.0, EVAPORATION),
+    bound(BoundedQuantityKind::Quality, 0.0, 1.0, EVAPORATION),
+];
+const TAO_BOUNDS: &[BoundedQuantity] = &[
+    bound(BoundedQuantityKind::MassFlux, 5.0, 600.0, CONDENSATION),
+    bound(BoundedQuantityKind::Quality, 0.0, 1.0, CONDENSATION),
+];
 
 const fn bound(
     kind: BoundedQuantityKind,
@@ -481,6 +515,42 @@ pub fn correlation_metadata(id: CorrelationId) -> CorrelationMetadata {
             RefrigerantApplicability::EvidenceIncomplete,
             "Muller-Steinhagen, H., Heck, K. (1986). Chem. Eng. Process. 20, 297-308",
         ),
+        CorrelationId::SimplifiedChannel => (
+            "BPHX SimplifiedChannel (default)",
+            CorrelationPurpose::PressureDrop,
+            PLATES,
+            SINGLE_PHASE_REGIMES,
+            BPHX_DP_BOUNDS,
+            RefrigerantApplicability::FluidAgnostic,
+            "Entropyk channel Blasius/laminar C¹ blend; ΔP = z_dp·2 f L G²/(ρ dh)",
+        ),
+        CorrelationId::Martin1996 => (
+            "Martin (1996/1999) chevron plate friction",
+            CorrelationPurpose::PressureDrop,
+            PLATES,
+            SINGLE_PHASE_REGIMES,
+            BPHX_DP_BOUNDS,
+            RefrigerantApplicability::FluidAgnostic,
+            "Martin, H. (1996). Chem. Eng. Process. 35:301–310; Martin 1999 / Shah appendix form",
+        ),
+        CorrelationId::AmalfiThome2016 => (
+            "Amalfi–Thome (2016) boiling plate ΔP [not implemented]",
+            CorrelationPurpose::PressureDrop,
+            PLATES,
+            EVAPORATION,
+            AMALFI_BOUNDS,
+            RefrigerantApplicability::EvidenceIncomplete,
+            "Amalfi, Vakili-Farahani, Thome (2016). Int. J. Refrigeration — Parts 1–2; deferred",
+        ),
+        CorrelationId::TaoInfanteFerreiraSurvey => (
+            "Tao & Infante Ferreira condensation BPHE DP survey [not implemented]",
+            CorrelationPurpose::PressureDrop,
+            PLATES,
+            CONDENSATION,
+            TAO_BOUNDS,
+            RefrigerantApplicability::EvidenceIncomplete,
+            "Tao & Infante Ferreira — condensation-in-BPHE pressure-drop survey; deferred",
+        ),
     };
     CorrelationMetadata {
         id,
@@ -516,6 +586,10 @@ pub fn registered_correlations() -> Vec<CorrelationMetadata> {
         CorrelationId::Mostinski1963,
         CorrelationId::Friedel1979,
         CorrelationId::MullerSteinhagenHeck1986,
+        CorrelationId::SimplifiedChannel,
+        CorrelationId::Martin1996,
+        CorrelationId::AmalfiThome2016,
+        CorrelationId::TaoInfanteFerreiraSurvey,
     ]
     .into_iter()
     .map(correlation_metadata)
@@ -1369,16 +1443,23 @@ mod tests {
 
     #[test]
     fn registry_does_not_claim_unverified_fluid_agnostic_coverage() {
-        assert!(registered_correlations().iter().all(|metadata| matches!(
-            metadata.domain.refrigerants,
-            RefrigerantApplicability::EvidenceIncomplete
-        )));
+        // Geometry-based single-phase BPHX channel DP may declare FluidAgnostic;
+        // HTC / two-phase entries must stay EvidenceIncomplete until verified.
+        assert!(registered_correlations().iter().all(|metadata| {
+            matches!(
+                metadata.domain.refrigerants,
+                RefrigerantApplicability::EvidenceIncomplete
+            ) || matches!(
+                metadata.id,
+                CorrelationId::SimplifiedChannel | CorrelationId::Martin1996
+            )
+        }));
     }
 
     #[test]
     fn central_registry_contains_heat_transfer_and_pressure_drop_entries() {
         let registered = registered_correlations();
-        assert_eq!(registered.len(), 14);
+        assert_eq!(registered.len(), 18);
         assert!(registered
             .iter()
             .any(|entry| entry.id == CorrelationId::Friedel1979
@@ -1387,6 +1468,18 @@ mod tests {
             entry.id == CorrelationId::MullerSteinhagenHeck1986
                 && entry.domain.purpose == CorrelationPurpose::PressureDrop
         }));
+        assert!(registered.iter().any(|entry| {
+            entry.id == CorrelationId::SimplifiedChannel
+                && entry.domain.purpose == CorrelationPurpose::PressureDrop
+        }));
+        assert!(registered.iter().any(|entry| {
+            entry.id == CorrelationId::Martin1996
+                && entry.domain.purpose == CorrelationPurpose::PressureDrop
+        }));
+        assert!(registered.iter().any(|entry| {
+            entry.id == CorrelationId::AmalfiThome2016
+                && entry.domain.purpose == CorrelationPurpose::PressureDrop
+        }));
         assert!(registered
             .iter()
             .any(|entry| entry.id == CorrelationId::Cooper1984));

crates/components/src/heat_exchanger/exchanger.rs

@@ -163,6 +163,15 @@ impl HxSideConditions {
     }
 }
 
+/// Live four-port edge indices: each triple is `(ṁ_idx, P_idx, h_idx)`.
+#[derive(Debug, Clone, Copy)]
+pub(crate) struct FourPortEdgeIndices {
+    pub hot_in: (usize, usize, usize),
+    pub hot_out: (usize, usize, usize),
+    pub cold_in: (usize, usize, usize),
+    pub cold_out: (usize, usize, usize),
+}
+
 /// Generic heat exchanger component with 4 ports.
 ///
 /// Uses the Strategy Pattern for heat transfer calculations via the
@@ -516,7 +525,7 @@ impl<Model: HeatTransferModel + 'static> HeatExchanger<Model> {
     }
 
     /// `true` when all 4 edges are wired (Modelica-style 4-port mode).
-    fn edges_ready(&self) -> bool {
+    pub(crate) fn edges_ready(&self) -> bool {
         self.hot_in_idx.is_some()
             && self.hot_out_idx.is_some()
             && self.cold_in_idx.is_some()
@@ -525,6 +534,60 @@ impl<Model: HeatTransferModel + 'static> HeatExchanger<Model> {
             && !self.cold_fluid_id_str.is_empty()
     }
 
+    /// Live four-port edge index triples `(ṁ, P, h)` when wired.
+    pub(crate) fn four_port_edges(&self) -> Option<FourPortEdgeIndices> {
+        if !self.edges_ready() {
+            return None;
+        }
+        Some(FourPortEdgeIndices {
+            hot_in: self.hot_in_idx?,
+            hot_out: self.hot_out_idx?,
+            cold_in: self.cold_in_idx?,
+            cold_out: self.cold_out_idx?,
+        })
+    }
+
+    /// Thermal-model equation count (excludes the two pressure-closure rows).
+    pub(crate) fn model_n_equations(&self) -> usize {
+        self.model.n_equations()
+    }
+
+    /// Calibration state indices (e.g. free `z_dp` actuator).
+    pub(crate) fn calib_indices_ref(&self) -> &entropyk_core::CalibIndices {
+        &self.calib_indices
+    }
+
+    /// Hot-side fluid id used for live property queries.
+    pub(crate) fn hot_fluid_id_str(&self) -> &str {
+        &self.hot_fluid_id_str
+    }
+
+    /// Cold-side fluid id used for live property queries.
+    pub(crate) fn cold_fluid_id_str(&self) -> &str {
+        &self.cold_fluid_id_str
+    }
+
+    /// Inlet density [kg/m³] from the attached fluid backend at `(P, h)`.
+    pub(crate) fn side_density(
+        &self,
+        side: &str,
+        fluid_id: &str,
+        p_pa: f64,
+        h_jkg: f64,
+    ) -> Result<f64, ComponentError> {
+        self.query_live_property(side, fluid_id, Property::Density, p_pa, h_jkg)
+            .and_then(|rho| {
+                if rho.is_finite() && rho > 1e-10 {
+                    Ok(rho)
+                } else {
+                    Err(ComponentError::DomainViolation(DomainViolation {
+                        component: Some(self.name.clone()),
+                        detail: format!("{} {}-side density is invalid: {}", self.name, side, rho),
+                    }))
+                }
+            })
+    }
+
     fn live_state_required_error(&self) -> ComponentError {
         ComponentError::InvalidState(format!(
             "{} requires live four-port edge state (hot_inlet, hot_outlet, cold_inlet, cold_outlet); inlet-only boundary conditions cannot define outlet states",

crates/components/src/heat_exchanger/mod.rs

@@ -50,6 +50,7 @@
 pub mod air_cooled_condenser;
 pub mod bphx_condenser;
 pub mod bphx_correlation;
+pub mod bphx_dp;
 pub mod bphx_evaporator;
 pub mod bphx_exchanger;
 pub mod bphx_geometry;
@@ -83,6 +84,7 @@ pub use bphx_correlation::{
     BphxCorrelation, CorrelationEvaluation, CorrelationParams, CorrelationResult,
     CorrelationSelector, ValidityStatus,
 };
+pub use bphx_dp::{evaluate_channel_pressure_drop, BphxDpCorrelation, BphxDpEval};
 pub use bphx_evaporator::BphxEvaporator;
 pub use bphx_exchanger::BphxExchanger;
 pub use bphx_geometry::{BphxGeometry, BphxGeometryBuilder, BphxGeometryError, BphxType};

crates/components/src/lib.rs

@@ -132,11 +132,11 @@ pub use free_cooling_exchanger::{
 };
 pub use heat_exchanger::model::FluidState;
 pub use heat_exchanger::{
-    AirCooledCondenser, CoilGeometry, Condenser, CondenserCoil, CondenserRating, Economizer,
-    EpsNtuModel, Evaporator, EvaporatorCoil, EvaporatorRating, ExchangerType, FanCoilUnit,
-    FinCoilCondenser, FinType, FloodedCondenser, FloodedEvaporator, FloodedPoolBoilingConfig,
-    FlowConfiguration, GasCooler, HeatExchanger, HeatExchangerBuilder, HeatTransferModel,
-    HxSideConditions, LmtdModel, MchxCondenserCoil, ShellAndTubeHx, UaMode,
+    AirCooledCondenser, BphxDpCorrelation, CoilGeometry, Condenser, CondenserCoil, CondenserRating,
+    Economizer, EpsNtuModel, Evaporator, EvaporatorCoil, EvaporatorRating, ExchangerType,
+    FanCoilUnit, FinCoilCondenser, FinType, FloodedCondenser, FloodedEvaporator,
+    FloodedPoolBoilingConfig, FlowConfiguration, GasCooler, HeatExchanger, HeatExchangerBuilder,
+    HeatTransferModel, HxSideConditions, LmtdModel, MchxCondenserCoil, ShellAndTubeHx, UaMode,
 };
 pub use heat_source::HeatSource;
 pub use isenthalpic_expansion_valve::IsenthalpicExpansionValve;

crates/entropyk/src/lib.rs

@@ -174,8 +174,9 @@ pub use builder::{SystemBuilder, SystemBuilderError};
 
 mod result;
 pub use result::{
-    extract_simulation_result, ComponentResult, ConvergenceSummary, EdgeResult, EnergyResult,
-    PortState, SimulationOutcome, SimulationResult, SystemSummary,
+    extract_simulation_result, extract_solved_variables, ComponentResult, ConvergenceSummary,
+    EdgeResult, EnergyResult, PortState, SimulationOutcome, SimulationResult, SolvedVariable,
+    SystemSummary,
 };
 
 // =============================================================================

crates/entropyk/src/result.rs

@@ -173,6 +173,36 @@ pub struct SystemSummary {
     pub cop_heating: Option<f64>,
 }
 
+// ─────────────────────────────────────────────────────────────────────────────
+// Solved unknown (free actuator / calibration factor)
+// ─────────────────────────────────────────────────────────────────────────────
+
+/// A named solver-computed unknown with its converged value and physical bounds.
+///
+/// These correspond to bounded control variables that the Newton solve treats
+/// as unknowns — free actuators (e.g. expansion-valve `opening`, condenser
+/// `fan_speed`), hard inverse-control links, and saturated-controller
+/// actuators (calibration factors like `z_ua`, `z_dp`, `z_flow`). They ride
+/// inside the flat `raw_state_vector` returned by the solver; this struct
+/// attaches a stable id, owning component, short label, and `[min, max]`
+/// bounds so the UI can surface them without parsing the raw vector.
+#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
+#[serde(rename_all = "camelCase")]
+pub struct SolvedVariable {
+    /// Stable id from [`BoundedVariableId`] (e.g. `"exv__opening"`).
+    pub id: String,
+    /// Owning component name, when known (`None` for globals).
+    pub component: Option<String>,
+    /// Short human label, e.g. `"opening"`, `"z_ua"`, `"z_dp"`.
+    pub variable: String,
+    /// Converged value of the unknown.
+    pub value: f64,
+    /// Lower bound (inclusive).
+    pub min: f64,
+    /// Upper bound (inclusive).
+    pub max: f64,
+}
+
 // ─────────────────────────────────────────────────────────────────────────────
 // Top-level SimulationResult
 // ─────────────────────────────────────────────────────────────────────────────
@@ -195,6 +225,12 @@ pub struct SimulationResult {
     pub edges: Vec<EdgeResult>,
     /// Aggregated system performance summary.
     pub summary: SystemSummary,
+    /// Named solver-computed unknowns (free actuators + calibration factors).
+    ///
+    /// Empty for systems with no bounded control variables. Skipped on
+    /// serialization when empty so older consumers keep working.
+    #[serde(default, skip_serializing_if = "Vec::is_empty")]
+    pub solved_variables: Vec<SolvedVariable>,
 }
 
 impl SimulationResult {
@@ -422,6 +458,9 @@ pub fn extract_simulation_result(system: &System, converged: &ConvergedState) ->
         status,
     };
 
+    // --- Solved unknowns (free actuators + calibration factors) ---
+    let solved_variables = extract_solved_variables(system, state);
+
     SimulationResult {
         status,
         convergence,
@@ -429,6 +468,93 @@ pub fn extract_simulation_result(system: &System, converged: &ConvergedState) ->
         components,
         edges,
         summary,
+        solved_variables,
+    }
+}
+
+/// Extracts named solver unknowns (bounded variables / free actuators) from a
+/// solved state vector.
+///
+/// Every bounded variable registered on `system` whose value lives in the
+/// Newton state vector — as a physical free actuator, a hard inverse-control
+/// link, or a saturated-controller actuator — is materialized into a
+/// [`SolvedVariable`] carrying its converged value and `[min, max]` bounds.
+///
+/// Bounded variables with no resolvable state slot (purely registered, never
+/// wired into the solve) are silently skipped; a single `tracing::debug!`
+/// line summarizes how many were skipped so coverage gaps are diagnosable
+/// without panic risk. Out-of-range indices (state shorter than expected,
+/// e.g. a truncated warm-start vector) are likewise skipped.
+///
+/// # Arguments
+///
+/// * `system` - The solved system (must be finalized).
+/// * `state` - Converged Newton state vector slice (length =
+///   `system.full_state_vector_len()` when complete).
+pub fn extract_solved_variables(system: &System, state: &[f64]) -> Vec<SolvedVariable> {
+    let mut out = Vec::new();
+    let mut skipped = 0usize;
+
+    for bv in system.bounded_variables() {
+        let id = bv.id();
+        // Free actuators are laid out in their own block at the tail of the
+        // state vector; everything else (hard-control links and saturated-
+        // controller actuators) is resolvable via control_variable_state_index.
+        let idx = system
+            .free_actuators()
+            .position(|fid| fid == id)
+            .map(|i| system.free_actuator_index(i))
+            .or_else(|| system.control_variable_state_index(id));
+
+        match idx {
+            Some(i) if i < state.len() => {
+                let id_str = id.as_str();
+                let component = bv.component_id().map(|s| s.to_string());
+                let variable = derive_solved_variable_label(id_str, bv.component_id());
+                out.push(SolvedVariable {
+                    id: id_str.to_string(),
+                    component,
+                    variable,
+                    value: state[i],
+                    min: bv.min(),
+                    max: bv.max(),
+                });
+            }
+            _ => {
+                skipped += 1;
+            }
+        }
+    }
+
+    if skipped > 0 {
+        tracing::debug!(
+            skipped,
+            total = system.bounded_variable_count(),
+            "Bounded variables skipped when extracting solved variables (no resolvable state slot)"
+        );
+    }
+
+    out
+}
+
+/// Derives the short user-facing variable label from a bounded-variable id.
+///
+/// Bounded-variable ids follow the convention `"{component}__{factor}"`, e.g.
+/// `"exv__opening"`, `"cond__z_ua"`. When the component is known we strip the
+/// `"{component}__"` prefix; otherwise we fall back to the segment after the
+/// last `__`. The solver-internal `"actuator"` suffix (used for saturated
+/// controllers driving an `opening`/`injection` factor) is remapped to the
+/// user-facing `"opening"` label.
+fn derive_solved_variable_label(id: &str, component_id: Option<&str>) -> String {
+    let stripped = match component_id {
+        Some(c) => id
+            .strip_prefix(&format!("{}__", c))
+            .unwrap_or(id),
+        None => id.rfind("__").map(|i| &id[i + 2..]).unwrap_or(id),
+    };
+    match stripped {
+        "actuator" => "opening".to_string(),
+        other => other.to_string(),
     }
 }
 
@@ -572,6 +698,7 @@ mod tests {
             components: vec![],
             edges: vec![],
             summary: SystemSummary::default(),
+            solved_variables: vec![],
         };
         let json = result.to_json().unwrap();
         assert!(json.contains("\"status\": \"converged\""));
@@ -595,4 +722,65 @@ mod tests {
         assert_eq!(result.edges, de.edges);
         assert_eq!(result.summary, de.summary);
     }
+
+    #[test]
+    fn test_solved_variable_serialization() {
+        let sv = SolvedVariable {
+            id: "exv__opening".to_string(),
+            component: Some("exv".to_string()),
+            variable: "opening".to_string(),
+            value: 0.62,
+            min: 0.02,
+            max: 1.0,
+        };
+        let json = serde_json::to_string(&sv).unwrap();
+        let de: SolvedVariable = serde_json::from_str(&json).unwrap();
+        assert_eq!(sv, de);
+        assert!(json.contains("\"component\":\"exv\""));
+        assert!(json.contains("\"variable\":\"opening\""));
+    }
+
+    #[test]
+    fn test_solved_variables_skipped_when_empty() {
+        // Empty solved_variables should be omitted from serialized JSON
+        // (skip_serializing_if = "Vec::is_empty") for backward compatibility.
+        let result = SimulationResult {
+            status: SimulationOutcome::Converged,
+            convergence: ConvergenceSummary {
+                iterations: 1,
+                final_residual: 0.0,
+                converged: true,
+                status: SimulationOutcome::Converged,
+            },
+            metadata: entropyk_solver::SimulationMetadata::new("h".to_string()),
+            components: vec![],
+            edges: vec![],
+            summary: SystemSummary::default(),
+            solved_variables: vec![],
+        };
+        let json = result.to_json().unwrap();
+        assert!(!json.contains("solvedVariables"));
+    }
+
+    #[test]
+    fn test_derive_solved_variable_label() {
+        // Component known → strip "{component}__" prefix.
+        assert_eq!(
+            derive_solved_variable_label("exv__opening", Some("exv")),
+            "opening"
+        );
+        assert_eq!(
+            derive_solved_variable_label("cond__z_ua", Some("cond")),
+            "z_ua"
+        );
+        // Solver-internal "actuator" suffix remapped to user-facing "opening".
+        assert_eq!(
+            derive_solved_variable_label("exv__actuator", Some("exv")),
+            "opening"
+        );
+        // Component unknown → fall back to last "__"-separated segment.
+        assert_eq!(derive_solved_variable_label("glob__z_dp", None), "z_dp");
+        // No separator and no component → id itself.
+        assert_eq!(derive_solved_variable_label("global_var", None), "global_var");
+    }
 }