# Regenerative Cooling

This page explains the regenerative cooling model implemented in `pyskyfire.regen`. It traces the path through `coupled_solver.py`, shows how the solver calls the equations in `physics.py`, and explains how the hot-gas aerothermodynamic properties are prepared by `pyskyfire.skycea.aerothermodynamics`.

The model is a quasi-one-dimensional steady heat-exchanger calculation. At each axial station, it balances three heat-transfer paths:

1. heat transfer from the combustion gas to the hot wall,
2. conduction through one or more wall layers, and
3. heat transfer from the cold wall into the coolant flowing in the cooling channels.

The result is a streamwise solution for wall temperatures, coolant bulk temperature, coolant pressure, heat flux, heat-transfer coefficients, coolant velocity, and residuals.

---

## 1. Code path overview

The public entry point is:

```python
coupled_steady_heating_analysis(
    thrust_chamber,
    boundary_conditions,
    nodes=100,
    circuit_index=0,
    film="auto",
    solver="newton",
    output=True,
)
```

For the current implementation, the only accepted solver name is `"newton"`. This dispatches to:

```python
solve_coupled_heat_exchanger(
    thrust_chamber,
    boundary_conditions,
    nodes,
    circuit_index,
    output,
)
```

The coolant temperature and pressure are marched explicitly, but the two wall temperatures are obtained from a local nonlinear heat-balance solve using `scipy.optimize.least_squares`.

`nodes` may be an integer or three explicit, potentially ragged grids. An
integer constructs identical uniform grids. Explicit grids use
`[wall_nodes, heat_flux_nodes, coolant_nodes]`; the equivalent named form is:

```python
nodes = {
    "wall": [...],
    "heat_flux": [...],
    "coolant": [...],
}
```

Every coolant interval must contain at least one wall node. Multiple wall
nodes are attached to the same lumped coolant state, and their heat rates are
integrated before the next coolant state is evaluated. This is useful when a
detailed wall/heat-flux distribution is needed but coolant property flashes
dominate runtime.

The helper class:

```python
CoupledHeatExchangerPhysics(thrust_chamber, boundary_conditions, circuit_index)
```

collects the local physics calculations used by the marching solver:

- `hot_side_coefficients(x, T_hw)`
- `cold_side_coefficients(x, T_cw, T_cool)`
- `dQ_hot_dx(x, T_hw)`
- `dQ_cond_dx(x, T_hw, T_cw)`
- `dQ_cold_dx(x, T_cw, T_cool)`
- `coolant_temperature_rate(T_cool, p_cool, dQ_cold_dx)`
- `coolant_friction_rate(x, T_cool, p_cool)`
- `bulk_velocity(x, T_cool, p_cool)`
- `interface_temperatures(x, T_hw, T_cw)`

These methods call lower-level correlations from `physics.py`, mainly Bartz-type hot-gas heat transfer, Colburn coolant-side heat transfer, Reynolds number, Darcy friction factor, coolant velocity, curvature factor, and adiabatic-wall temperature.

The hot-gas property calls such as:

```python
combustion_transport.get_T(x)
combustion_transport.get_p(x)
combustion_transport.get_h(x)
combustion_transport.get_cp(x)
combustion_transport.get_mu(x)
combustion_transport.get_k(x)
combustion_transport.get_Pr(x)
combustion_transport.get_M(x)
combustion_transport.get_a(x)
combustion_transport.get_gamma(x)
```

are normally supplied by `skycea.aerothermodynamics.Aerothermodynamics`, which precomputes CEA-based equilibrium and temperature-pressure maps along the chamber/nozzle contour.

---

## 2. Symbols and sign conventions

The most important local variables are:

| Symbol | Meaning | Units |
|---|---:|---:|
| `x` | axial coordinate | m |
| `dx` | axial node spacing used by the marching solver | m |
| `T_hw` | hot-side wall temperature | K |
| `T_cw` | coolant-side wall temperature | K |
| `T_cool` | coolant bulk static temperature | K |
| `p_static` | coolant static pressure | Pa |
| `p_stagnation` | coolant stagnation pressure | Pa |
| `dA_hot/dx` | hot-side thermal area per unit axial length | m |
| `A_cool` | flow area of one coolant channel | m² |
| `D_h` | coolant hydraulic diameter | m |
| `D_hyd` | hot-gas-side hydraulic diameter approximation | m |
| `mdot_c` | coolant mass flow through one channel | kg/s |
| `mdot_g` | hot-gas mass flow | kg/s |
| `h_hot` | effective temperature-based hot-side coefficient | W/m²/K |
| `h_g` | enthalpy-based hot-side coefficient | kg/m²/s |
| `h_cold` | coolant-side heat-transfer coefficient | W/m²/K |
| `qpp_hot` | hot-side wall heat flux | W/m² |
| `dQ_hot_dx` | hot-side heat input per unit axial length | W/m |
| `dQ_cond_dx` | wall conduction per unit axial length | W/m |
| `dQ_cold_dx` | coolant heat pickup per unit axial length | W/m |

Heat flow is positive from the combustion gas into the wall and then into the coolant.

---

## 3. Boundary conditions

The `BoundaryConditions` class stores the inlet state of one cooling circuit:

```python
BoundaryConditions(T_coolant_in, p_coolant_in, mdot_coolant)
```

where:

$$
T_{c,in} = \texttt{T\_coolant\_in}
$$

$$
p_{c,in} = \texttt{p\_coolant\_in}
$$

$$
\dot{m}_{c,total} = \texttt{mdot\_coolant}
$$

The mass flow supplied here is the total mass flow through the selected cooling circuit. The solver divides it evenly across all geometric channel instances:

$$
N_{chan} = N_{positions} N_{channels/leaf}
$$

$$
\dot{m}_{c,chan} = \frac{\dot{m}_{c,total}}{N_{chan}}
$$

This single-channel mass flow is used for local velocity, Reynolds number, heat-transfer coefficient, coolant temperature rise, and pressure drop.

---

## 4. Aerothermodynamics from `skycea.aerothermodynamics`

The regenerative cooling solver needs the hot-gas state along the contour. In the CEA-based workflow, this is provided by the `Aerothermodynamics` class.

### 4.1 Initialization from thrust, area ratio, and characteristic length

One constructor is:

```python
Aerothermodynamics.from_F_eps_Lstar(
    fu, ox, MR, p_c, F, eps, L_star,
    T_fu_in=298.15,
    T_ox_in=298.15,
    p_amb=1.013e5,
)
```

The oxidizer-to-fuel mixture ratio is:

$$
MR = \frac{\dot{m}_{ox}}{\dot{m}_{fu}}
$$

The method builds CEA fuel and oxidizer objects at the specified inlet temperatures and runs a `CEA_Wrap.RocketProblem` at chamber pressure and exit area ratio. Internally, chamber pressure is converted from pascal to psi using:

$$
p[\mathrm{psi}] = 0.000145038\,p[\mathrm{Pa}]
$$

CEA returns design-point quantities such as characteristic velocity, vacuum specific impulse, chamber density, chamber temperature, throat temperature, and throat pressure. These are then used to derive engine-level quantities.

The total mass flow is calculated from the requested thrust and vacuum specific impulse:

$$
\dot{m} = \frac{F}{I_{sp,vac} g_0}
$$

using:

$$
g_0 = 9.81~\mathrm{m/s^2}
$$

The fuel and oxidizer mass flows are:

$$
\dot{m}_{fu} = \frac{\dot{m}}{1 + MR}
$$

$$
\dot{m}_{ox} = \dot{m} - \dot{m}_{fu}
$$

The throat area follows from the definition of characteristic velocity:

$$
c^* = \frac{p_c A_t}{\dot{m}}
$$

so that:

$$
A_t = \frac{c^* \dot{m}}{p_c}
$$

The throat radius is:

$$
r_t = \sqrt{\frac{A_t}{\pi}}
$$

The exit area and exit radius are:

$$
A_e = \varepsilon A_t
$$

$$
r_e = \sqrt{\frac{A_e}{\pi}}
$$

The chamber residence time is computed from the characteristic length, throat area, chamber density, and total mass flow:

$$
t_{stay} = \frac{L^* A_t \rho_c}{\dot{m}}
$$

The chamber volume estimate is then:

$$
V_c = \frac{\dot{m} t_{stay}}{\rho_c}
$$

which reduces to:

$$
V_c = L^* A_t
$$

when substituting the previous expression.

The vacuum thrust coefficient is derived from vacuum specific impulse:

$$
C_{F,vac} = \frac{I_{sp,vac} g_0}{c^*}
$$

Ambient thrust coefficient is modeled by subtracting the pressure-thrust penalty:

$$
C_{F,amb} = C_{F,vac} - \frac{p_{amb}}{p_c}\frac{A_e}{A_t}
$$

and the corresponding ambient specific impulse is:

$$
I_{sp,amb} = \frac{C_{F,amb} c^*}{g_0}
$$

The sea-level thrust coefficient and sea-level specific impulse use:

$$
p_{SL} = 1.01325\times10^5~\mathrm{Pa}
$$

$$
C_{F,SL} = C_{F,vac} - \frac{p_{SL}}{p_c}\frac{A_e}{A_t}
$$

$$
I_{sp,SL} = \frac{C_{F,SL} c^*}{g_0}
$$

### 4.2 Initialization from thrust and exit pressure

The alternative constructor is:

```python
Aerothermodynamics.from_F_pe_Lstar(
    fu, ox, MR, p_c, F, p_e, L_star,
    T_fu_in=298.15,
    T_ox_in=298.15,
    p_amb=1.013e5,
)
```

Instead of specifying exit area ratio directly, it specifies exit pressure. CEA is called with the pressure ratio:

$$
\frac{p_c}{p_e}
$$

The area ratio is then obtained from the CEA result:

$$
\varepsilon = \frac{A_e}{A_t}
$$

After that, the same mass-flow, area, radius, residence-time, and thrust-coefficient equations are used.

### 4.3 Attaching the contour

After construction, the method:

```python
compute_aerothermodynamics(contour)
```

stores the contour on the object and clears the station cache. No property
table is built here: every gas state is solved live by CEA when a getter asks
for it, so there is no resolution to choose at this point. `attach_contour` is
an alias for the same method.

### 4.4 Station states along the contour

With no temperature or enthalpy argument, a getter such as `get_T(x)` solves
the rocket problem at the local area ratio:

$$
\varepsilon(x) = \frac{A(x)}{A_t}
$$

The station is passed to CEA as a subsonic or supersonic area ratio depending
on the sign of `x`:

$$
x < 0 \Rightarrow \text{subsonic solve at } \varepsilon(x)
$$

$$
x \ge 0 \Rightarrow \text{supersonic solve at } \varepsilon(x)
$$

Stations inside the tolerance band $\varepsilon(x) \le 1 + 10^{-10}$ are solved
as the chamber state. CEA returns the state in its own units, which the code
converts to base SI: bar to pascal, kJ to J, millipoise to Pa·s, and
mW/(cm·K) to W/(m·K).

Because a station state depends only on `x`, it is cached on the object and
reused by every later getter at the same coordinate. Temperature- and
enthalpy-conditioned states are not cached: they are Newton trial points that
almost never recur, and caching them made the cache grow without bound during
a coupled solve. The coupled solver clears the cache when a simulation
finishes.

The cache is what makes post-processing cheap. A report tab that plots eleven
properties over the same axial grid pays for one CEA solve per station in the
first plot; the remaining ten read cached states. Visualization therefore
chooses its own resolution, and `PlotTransportProperty` accepts exactly one of:

```python
PlotTransportProperty(transport, prop="T", results=cooling_data)  # run grid
PlotTransportProperty(transport, prop="T", nodes=200)             # uniform
PlotTransportProperty(transport, prop="T", x=my_x)                # explicit
```

Passing `results` reuses the axial stations of the solved run, which keeps the
plots on the same grid as the simulation and hits the cache the run already
populated.

### 4.5 Imposed temperature and imposed enthalpy

When a temperature is supplied, for example:

```python
get_h(x, T=T_wall)
```

the local static pressure is taken from the station state at `x`, and CEA is
solved as a TP problem at $(T,\ p(x))$. An imposed enthalpy is solved the same
way as an HP problem at $(h,\ p(x))$.

Pressure is special: `get_p(x, T=..., h=...)` ignores `T` and `h` and always
returns the station pressure $p(x)$.

Composition is returned directly from the solve as a species mole-fraction
dictionary; no interpolation between stations is involved.

### 4.6 Continuation below the CEA temperature boundary

CEA becomes unreliable at low temperature, so states below
`minimum_cea_temperature` (200 K by default, configurable on the constructor)
are not solved directly. Instead, the properties are continued along their
right-hand tangent at the boundary $T_b$. Three TP solves are taken at
$T_b,\ T_b + \Delta T,\ T_b + 2\Delta T$, and a second-order one-sided
derivative is formed:

$$
\left.\frac{dZ}{dT}\right|_{T_b} = \frac{-3Z_0 + 4Z_1 - Z_2}{2\Delta T}
$$

Strictly positive properties ($\rho$, $c_p$, $\gamma$, $a$, $\mu$, $k$, $Pr$,
$MW$) are continued in log space, which keeps them positive:

$$
Z(T) = Z(T_b)\exp\left[\left.\frac{d\ln Z}{dT}\right|_{T_b}(T - T_b)\right]
$$

Enthalpy is continued linearly:

$$
h(T) = h(T_b) + \left.\frac{dh}{dT}\right|_{T_b}(T - T_b)
$$

This is $C^1$ at the boundary, so a Newton solve that crosses $T_b$ sees no
kink in the properties or their first derivatives. An HP query that fails, or
that lands below $T_b$, inverts the enthalpy tangent to recover the
temperature and then evaluates the same continuation. The tangents are cached
per (station, pressure) pair and cleared together with the station cache.

---

## 5. Marching grid and unknowns in `solver.py`

The cooling solver reads the selected cooling circuit:

```python
circuit = thrust_chamber.cooling_circuits[circuit_index]
```

and builds a marching grid from `circuit.x_domain`. If the circuit direction is positive:

$$
x_i = \mathrm{linspace}(x_0, x_N, n_{nodes})
$$

If the circuit direction is negative, the order is reversed:

$$
x_i = \mathrm{linspace}(x_N, x_0, n_{nodes})
$$

The axial spacing used in the solver is:

$$
\Delta x = \left|\frac{x_{N-1}-x_0}{n_{nodes}-1}\right|
$$

The stored solution arrays are:

- `T_hw_arr`: hot-side wall temperature,
- `T_cw_arr`: coolant-side wall temperature,
- `T_cool_arr`: coolant bulk static temperature,
- `p_static_arr`: coolant static pressure,
- `p_stagnation_arr`: coolant stagnation pressure,
- `dQ_dA_arr`: local heat flux,
- `velocity_arr`: coolant velocity,
- `T_stagnation_arr`: coolant stagnation temperature.

The inlet values are:

$$
T_{cool,0} = T_{c,in}
$$

$$
p_{static,0} = p_{c,in}
$$

$$
p_{0,0} = p_{c,in}
$$

The initial guesses for both wall temperatures at the first station are:

$$
T_{hw,guess} = \frac{T_g(x_0)+T_{c,in}}{2}
$$

$$
T_{cw,guess} = \frac{T_g(x_0)+T_{c,in}}{2}
$$

At each station, the unknowns solved by least squares are:

$$
T_{cw}
$$

and:

$$
\Delta T_w = T_{hw} - T_{cw}
$$

so that:

$$
T_{hw} = T_{cw} + \Delta T_w
$$

The lower bound enforces:

$$
T_{cw} \ge T_{cool}
$$

and:

$$
\Delta T_w \ge 0
$$

so the wall ordering remains:

$$
T_{hw} \ge T_{cw} \ge T_{cool}
$$

The hot-wall upper bound is based on the local gas temperature:

$$
T_{hw,max} = \max(T_{cool,i}+1,\ T_{g,i}-10^{-3})
$$

---

## 6. Hot-gas-side heat transfer

The hot-side calculation is performed by:

```python
CoupledHeatExchangerPhysics.hot_side_coefficients(x, T_hw)
```

### 6.1 Local gas-side geometry and state

The hot-gas-side hydraulic diameter is approximated as twice the local chamber radius:

$$
D_{hyd,g}(x) = 2r(x)
$$

The local hot-gas flow area is:

$$
A_g(x) = A(x)
$$

The hot-gas mass flow is:

$$
\dot{m}_g = \texttt{combustion\_transport.mdot}
$$

The local gas static temperature is:

$$
T_g = T_g(x)
$$

The temperature used in the current Bartz-property correction is the arithmetic mean of wall and gas temperature:

$$
T_{gr} = \frac{T_{hw}+T_g}{2}
$$

The gas enthalpy at the equilibrium state is:

$$
H_g = H_g(x)
$$

The hot-wall enthalpy is requested from the combustion-transport model at the wall temperature:

$$
H_{hw} = H(x,T_{hw})
$$

If this lookup fails, the implementation falls back to:

$$
H_{hw} = H_g
$$

The local Mach number and speed of sound are:

$$
M_g = M(x)
$$

$$
a_g = a(x)
$$

A reference enthalpy is computed:

$$
H_{gr} = \frac{1}{2}(H_{hw}+H_g) + 0.18\left(\frac{1}{2}M_g^2 a_g^2\right)
$$

Implementation note: 

Currently `H_gr`is calculated, but not used to query reference condition gas properties. The actual properties used in the Bartz correlation are currently retrieved at `x` as shown below. The reference-condition properties should ideally be evaluated at `H_gr`, but the current aerothermodynamics module does not yet support that workflow robustly. This is not trivial to solve, as it generally requires us to move away from NASA CEA as a reference lookup. There are currently no other programs that support as wide a library of propellants, so ditching it would hurt propellant compatibility strongly.

$$
c_{p,gr} \leftarrow c_p(x)
$$

$$
\mu_{gr} \leftarrow \mu(x)
$$

$$
k_{gr} \leftarrow k(x)
$$

$$
Pr_{gr} \leftarrow Pr(x)
$$

### 6.2 Bartz-type enthalpy-driven coefficient

The hot-side coefficient is calculated by `physics.h_gas_bartz_enthalpy_driven`:

$$
h_{gr} = 0.026\frac{k_{gr}}{D_{hyd,g}}
\left(\frac{c_{p,gr}}{k_{gr}\mu_{gr}}\right)^{0.4}
\left(\frac{\dot{m}_g D_{hyd,g}}{A_g}\right)^{0.8}
\left(\frac{T_g}{T_{gr}}\right)^{0.8}
$$

The solver multiplies this by a user-supplied hot-side correction factor:

$$
h_{gr,eff} = h_{gr}\,C_{hot}
$$

where:

$$
C_{hot} = \texttt{thrust\_chamber.h\_gas\_corr}
$$

The enthalpy-based heat-transfer coefficient used in the heat-flux equation is:

$$
h_g = \frac{h_{gr,eff}}{c_{p,gr}}
$$

Because `h_gr` has units of W/m²/K and `c_p` has units of J/kg/K, `h_g` has units of kg/m²/s. This is why the hot-side heat flux is written in terms of enthalpy difference rather than temperature difference.

### 6.3 Adiabatic-wall temperature and adiabatic-wall enthalpy

Adiabatic-wall temperature is computed using `physics.T_aw`:

$$
r = Pr^{1/3}
$$

$$
T_{aw} = T_\infty\left(1 + r\frac{\gamma-1}{2}M_\infty^2\right)
$$

In the hot-side call, this becomes:

$$
T_{aw} = T_g\left(1 + Pr_{gr}^{1/3}\frac{\gamma-1}{2}M_g^2\right)
$$

The heat flux itself is driven by adiabatic-wall enthalpy, not directly by `T_aw`. We compute:

$$
H_{aw} = H_g + \frac{1}{2}Pr_{gr}^{1/3}M_g^2 a_g^2
$$

### 6.4 Hot-side heat flux and heat per unit length

The hot-side wall heat flux is:

$$
q''_{hot} = h_g\left(H_{aw}-H_{hw}\right)
$$

The effective temperature-based hot-side coefficient reported for plotting and diagnostics is:

$$
h_{hot} = \frac{q''_{hot}}{T_{aw}-T_{hw}}
$$

The heat transfer per unit axial length is then:

$$
\frac{d\dot{Q}_{hot}}{dx} = q''_{hot}\frac{dA_{hot}}{dx}
$$

where `dA_hot/dx` comes from:

```python
cooling_circuit.dA_dx_thermal_exhaust(x)
```

---

## 7. Wall conduction

Wall conduction is calculated by:

```python
dQ_cond_dx(x, T_hw, T_cw)
```

The solver supports a stack of wall layers. Each wall layer has local thickness:

$$
\delta_j(x)
$$

and thermal conductivity evaluated at the mean wall temperature:

$$
k_j = k_j\left(\frac{T_{hw}+T_{cw}}{2}\right)
$$

The hot-side thermal area per unit axial length is:

$$
\frac{dA_{hot}}{dx}
$$

The thermal resistance per unit axial length of wall layer `j` is:

$$
R_j(x) = \frac{\delta_j(x)}{k_j\left(\frac{T_{hw}+T_{cw}}{2}\right)\frac{dA_{hot}}{dx}}
$$

The total wall resistance per unit length is the series sum:

$$
R_{tot}(x) = \sum_j R_j(x)
$$

The conduction heat flow per unit axial length is:

$$
\frac{d\dot{Q}_{cond}}{dx} = \frac{T_{hw}-T_{cw}}{R_{tot}}
$$

Interface temperatures are reconstructed after solving. Starting with:

$$
T_0 = T_{hw}
$$

and using:

$$
q'_x = \frac{d\dot{Q}_{cond}}{dx}
$$

for each wall layer:

$$
T_{j+1} = T_j - q'_x R_j
$$

The resulting list is:

$$
[T_{hw},\ T_1,\ T_2,\ \ldots,\ T_{cw}]
$$

The output array reverses this order so that the saved temperature vector at each station is:

$$
[T_{cool},\ T_{cw},\ \ldots,\ T_{hw}]
$$

---

## 8. Coolant-side heat transfer

The coolant-side coefficient is calculated by:

```python
cold_side_coefficients(x, T_cw, T_cool)
```

### 8.1 Coolant film and bulk properties

The film temperature is:

$$
T_{cf} = \frac{T_{cool}+T_{cw}}{2}
$$

The current implementation obtains a pressure for these property calls from the combustion-transport model:

$$
p \leftarrow p_g(x)
$$

and then evaluates:

$$
k_{cf} = k_c(T_{cf},p)
$$

$$
c_{p,cr} = c_{p,c}(T_{cf},p)
$$

$$
\mu_{cf} = \mu_c(T_{cf},p)
$$

The coolant bulk density and viscosity are evaluated at the coolant bulk temperature:

$$
\rho_c = \rho_c(T_{cool},p)
$$

$$
\mu_c = \mu_c(T_{cool},p)
$$

Implementation note: this pressure choice is specific to `cold_side_coefficients`. The coolant temperature and pressure marching functions use the coolant pressure array. If coolant properties are strongly pressure-dependent in a particular case, this distinction should be checked.

### 8.2 Coolant velocity and Reynolds number

The local single-channel flow area is:

$$
A_c(x) = \texttt{cooling\_circuit.A\_coolant(x)}
$$

The local hydraulic diameter is:

$$
D_c(x) = \texttt{cooling\_circuit.Dh\_coolant(x)}
$$

Coolant velocity is computed by `physics.u_coolant`:

$$
u_c = \frac{\dot{m}_{c,chan}}{\rho_c A_c}
$$

The Reynolds number is computed by `physics.reynolds`:

$$
Re_c = \frac{\rho_c u_c D_c}{\mu_c}
$$

### 8.3 Curvature factor

The local channel radius of curvature is:

$$
R_{curv}(x) = \texttt{cooling\_circuit.radius\_of\_curvature(x)}
$$

The curvature factor computed by `physics.phi_curv` is:

$$
\phi_{curv} = \left[Re_c\left(\frac{0.5D_c}{R_{curv}}\right)^2\right]^{0.05}
$$

For a straight section:

$$
R_{curv}=\infty \Rightarrow \phi_{curv}=1
$$

Implementation note: the current `cold_side_coefficients` method computes `phi_curv`, but then calls the Colburn correlation with `phi_curv=1`. Therefore, curvature is reported in the returned dictionary but is not currently applied to the coolant-side heat-transfer coefficient.

### 8.4 Colburn coolant-side heat-transfer coefficient

The coolant-side heat-transfer coefficient is calculated by `physics.h_coolant_colburn`:

$$
h_c = 0.023\frac{k_{cf}}{D_c}
\left(\frac{c_{p,cr}}{k_{cf}\mu_{cf}}\right)^{0.4}
\left(\frac{\dot{m}_{c,chan}D_c}{A_c}\right)^{0.8}
\phi_{curv}
$$

In the current solver call:

$$
\phi_{curv}=1
$$

The result is multiplied by a user-supplied coolant-side correction factor:

$$
h_{c,eff} = h_c C_{cold}
$$

where:

$$
C_{cold} = \texttt{thrust\_chamber.h\_cold\_corr}
$$

### 8.5 Coolant heat pickup per unit length

The coolant-side thermal resistance per unit length is delegated to the cooling-circuit geometry object:

```python
R_coolant_per_len(x, h_c=h_c, T_wall_rep=T_rep)
```

with representative wall/coolant temperature:

$$
T_{rep} = \frac{T_{cw}+T_{cool}}{2}
$$

The heat pickup per unit axial length is:

$$
\frac{d\dot{Q}_{cold}}{dx} = \frac{T_{cw}-T_{cool}}{R_{cool/len}}
$$

This lets the geometry object decide the effective coolant-side area and fin efficiency details, instead of hard-coding them directly in `solver.py`.

---

## 9. Local wall-temperature solve

At each axial station `i`, the solver forms three per-cell heat flows:

$$
Q_{hot,i} = \left(\frac{d\dot{Q}_{hot}}{dx}\right)_i \Delta x
$$

$$
Q_{cond,i} = \left(\frac{d\dot{Q}_{cond}}{dx}\right)_i \Delta x
$$

$$
Q_{cold,i} = \left(\frac{d\dot{Q}_{cold}}{dx}\right)_i \Delta x
$$

The local steady heat balance is:

$$
Q_{hot,i} = Q_{cond,i} = Q_{cold,i}
$$

The nonlinear residual vector is:

$$
R_1 = Q_{hot,i} - Q_{cond,i}
$$

$$
R_2 = Q_{cond,i} - Q_{cold,i}
$$

The residuals are scaled by:

$$
Q_{ref} = \max\left(|Q_{hot,i}|, |Q_{cond,i}|, |Q_{cold,i}|, 1\right)
$$

so the least-squares residual vector is:

$$
\mathbf{r} =
\begin{bmatrix}
R_1/Q_{ref}\\
R_2/Q_{ref}
\end{bmatrix}
$$

The solver uses `least_squares` with:

- trust-region reflective method, `method="trf"`,
- robust `soft_l1` loss,
- bounds enforcing physical wall-temperature ordering,
- `xtol = ftol = gtol = 1e-10`,
- `max_nfev = 200`.

After convergence, the solved values are:

$$
T_{cw,i} = T_{cw,sol}
$$

$$
T_{hw,i} = T_{cw,sol}+\Delta T_{w,sol}
$$

The converged wall temperatures become the initial guesses for the next axial station.

---

## 10. Coolant temperature marching

After the wall temperatures are solved at station `i`, the coolant temperature at station `i+1` is updated using the heat absorbed by the coolant in the current cell.

The continuous energy equation is:

$$
\frac{dT_{cool}}{dx} = \frac{1}{\dot{m}_{c,chan}c_p}\frac{d\dot{Q}_{cold}}{dx}
$$

In the implemented marching step, the solver first computes the per-cell coolant heat pickup:

$$
Q_{cold,i}=\left(\frac{d\dot{Q}_{cold}}{dx}\right)_i \Delta x
$$

and then applies:

$$
\Delta T_{cool,i} = \frac{Q_{cold,i}}{\dot{m}_{c,chan}c_p(T_{cool,i},p_i)}
$$

so:

$$
T_{cool,i+1} = T_{cool,i} + \Delta T_{cool,i}
$$

In the source, the function name and argument name still refer to `dQ_cold_dx`, but in the marching call the quantity passed has already been multiplied by `dx`. The implemented update is therefore a finite-volume energy update using heat per cell.

---

## 11. Coolant pressure marching

The solver marches **static** pressure, because static pressure is what the
momentum equation governs and what test data and the reference RL10 models
report. Stagnation pressure is carried alongside as
$p_0 = p + \tfrac{1}{2}\rho_c u_c^2$.

The irreversible part of the update is calculated by:

```python
coolant_friction_rate(x, T_cool, p_cool)
```

At the current station, the coolant density is:

$$
\rho_c = \rho_c(T_{cool},p_{cool})
$$

The single-channel velocity is:

$$
u_c = \frac{\dot{m}_{c,chan}}{\rho_c A_c}
$$

The Reynolds number is:

$$
Re_{D_h} = \frac{\rho_c u_c D_h}{\mu_c}
$$

The Darcy friction factor is then computed by `physics.f_darcy`.

### 11.1 Darcy friction factor

The following laminar threshold is used:

$$
Re_{lam}=2300
$$

and turbulent threshold:

$$
Re_{turb}=3500
$$

For laminar flow:

$$
f = \frac{64}{Re_{D_h}}
$$

For turbulent smooth-wall flow, a Petukhov-type expression is used:

$$
f = \left(0.79\ln Re_{D_h} - 1.64\right)^{-2}
$$

If roughness is supplied, the Colebrook-White equation is solved iteratively:

$$
\frac{1}{\sqrt{f}} + 2\log_{10}\left(
\frac{\epsilon(x)}{3.71D_h}
+ \frac{2.51}{Re_{D_h}\sqrt{f}}
\right)=0
$$

where:

$$
\epsilon(x) = \texttt{roughness(x)}
$$

In the transitional regime, the solver linearly blends the laminar and turbulent friction factors:

$$
f = \mathrm{interp}\left(Re_{D_h};\ [2300,3500],\ [f_{lam},f_{turb}]\right)
$$

### 11.2 Friction gradient

A geometric path-length factor is computed:

```python
circuit.ds_dx(x)
```

The friction gradient, evaluated at the segment midpoint on the upstream
state, is:

$$
\left(\frac{dp}{dx}\right)_{fric} = -\frac{f}{D_h}\frac{\rho_c u_c^2}{2}\frac{ds}{dx}
$$

where $f$ already carries the Ito curvature multiplier when
`pressure_curvature_correction` is enabled.

### 11.3 Acceleration term

The full one-dimensional momentum equation is:

$$
\frac{dp}{dx} = \left(\frac{dp}{dx}\right)_{fric} - \rho_c u_c\frac{du_c}{dx}
$$

The second term is not a purely local quantity: the coolant accelerates both
because the channel tapers and because heating lowers its density. It is
therefore applied segment to segment, trapezoidally in the mass flux
$G = \rho_c u_c$:

$$
\Delta p_{acc,i} = \tfrac{1}{2}\left(G_i + G_{i+1}\right)\left(u_{c,i+1} - u_{c,i}\right)
$$

This form is exact in both limits it has to span. At constant density it
reduces to Bernoulli, $\tfrac{1}{2}\rho_c(u_{c,i+1}^2 - u_{c,i}^2)$; at
constant area it reduces to $G^2\left(1/\rho_{i+1} - 1/\rho_i\right)$, which
for a heated duct is **twice** the change in dynamic head. Reconstructing
static pressure as $p_0 - \tfrac{1}{2}\rho_c u_c^2$ from a friction-only
stagnation march therefore counts the heating-driven acceleration at half
strength. Equivalently, a friction-only $dp_0/dx$ omits the heat-addition
(Rayleigh) stagnation loss $+\tfrac{1}{2}u_c^2\,d\rho_c/dx$.

### 11.4 Segment update

The static update is:

$$
p_{i+1} = p_i + \left(\frac{dp}{dx}\right)_{fric,i}\Delta x - \Delta p_{acc,i}
$$

$$
p_{0,i+1} = p_{i+1} + \tfrac{1}{2}\rho_{i+1}u_{c,i+1}^2
$$

The downstream density, velocity and pressure are mutually dependent through
the equation of state, so each segment is swept to a fixed point starting from
the friction-only guess. `PRESSURE_SWEEPS` caps the sweeps and
`PRESSURE_SWEEP_TOL` sets the convergence tolerance in pascals; a segment that
fails to settle raises a `RuntimeWarning`. The downstream enthalpy is fixed by
the segment's heat load and so does not participate in the sweep.

The inlet boundary condition `p_coolant_in` is interpreted as a **stagnation**
pressure, so node 0 starts at $p_0 = p_{coolant,in}$ and
$p = p_{coolant,in} - \tfrac{1}{2}\rho_c u_c^2$.

---

## 12. Derived outputs

After all nodes have been marched, the solver calculates several reporting quantities.

### 12.1 Heat flux

At each station:

$$
Q_{hot,i} = \left(\frac{d\dot{Q}_{hot}}{dx}\right)_i\Delta x
$$

$$
A_{hot,i} = \left(\frac{dA_{hot}}{dx}\right)_i\Delta x
$$

The reported heat flux is:

$$
q''_i = \frac{Q_{hot,i}}{A_{hot,i}}
$$

If `A_hot` is zero, the code stores zero to avoid division by zero.

### 12.2 Coolant velocity

The velocity reported at each node is:

$$
u_i = \frac{\dot{m}_{c,chan}}{\rho_i A_{c,i}}
$$

where density is evaluated using the solved coolant temperature and static pressure.

### 12.3 Stagnation temperature

The stagnation temperature is computed from the static coolant temperature and kinetic-energy term:

$$
T_{0,i} = T_i + \frac{u_i^2}{2c_{p,i}}
$$

where:

$$
c_{p,i}=c_p(T_i,p_i)
$$

### 12.4 Wall-interface temperature array

The returned `T` array has shape:

$$
(n_{nodes},\ 1+n_{walls}+1)
$$

The first column is coolant bulk temperature. The remaining columns are wall-interface temperatures ordered from cold side to hot side:

$$
T[i,:] = [T_{cool},\ T_{cw},\ T_{interface,1},\ldots,\ T_{hw}]
$$

### 12.5 Heat-transfer coefficients and adiabatic-wall temperature

For each node, the solver recomputes and stores:

$$
h_{hot}
$$

$$
h_g
$$

$$
h_{cold}
$$

$$
T_{aw}
$$

These are returned as:

```python
"h_hot"
"h_hot_enthalpy"
"h_cold"
"T_aw_hot"
```

---

## 13. Residual logging

When residual logging is enabled, every local wall-temperature solve records:

$$
(cell,\ iteration,\ R_1,\ R_2)
$$

The residual magnitude is:

$$
|R| = \sqrt{R_1^2+R_2^2}
$$

The global residual history is aggregated iteration-by-iteration. For finite `p`, the implemented norm is:

$$
R_{global,k} = \left(\mathrm{mean}\left(|R|^p\right)\right)^{1/p}
$$

For `p = \infty`, it is:

$$
R_{global,k} = \max(|R|)
$$

The final per-cell residual is the last recorded residual magnitude for each cell.

---


## 14. Assumptions and implementation details

The current regenerative cooling solver is an engineering heat-exchanger model. Its main assumptions are:

- The solution is steady in time.
- The solver marches along one spatial coordinate.
- Axial conduction in the wall is neglected.
- Wall conduction is treated as one-dimensional through the wall stack.
- Each coolant channel in the selected circuit receives an equal share of total circuit mass flow.
- Hot-gas heat transfer is modeled by an enthalpy-driven Bartz-style correlation.
- Coolant-side heat transfer is modeled by a Colburn-style turbulent internal-flow correlation.
- Coolant pressure drop is modeled with Darcy friction and a separate area-change correction.
- Hot-gas properties are solved live by CEA at the requested station and cached per axial coordinate.

There are also a few important implementation details to keep in mind:

- The Bartz gas properties are evaluated at the computed reference enthalpy `H_gr`.
- `phi_curv` is computed in the coolant-side model, but the current Colburn call passes `phi_curv=1`, so curvature does not currently modify `h_cold`.
- Gas states below `minimum_cea_temperature` are continued along a $C^1$ tangent instead of being solved by CEA.
- `get_p(x, T=..., h=...)` always returns the station pressure and ignores the supplied temperature or enthalpy.
- In `cold_side_coefficients`, coolant properties are evaluated using the marched coolant pressure.

These are limitations in the current implementation that is being worked on for future versions. Nevertheless, the exact behavior is reported here so results using the program are interpreted correctly. 

---

## 15. Returned data structure

`coupled_steady_heating_analysis` returns a `RegenResult` produced by `solve_coupled_heat_exchanger`. Its core fields are:

```python
{
    "x": x_wall,
    "x_wall": x_wall,
    "x_heat_flux": x_heat_flux,
    "x_coolant": x_coolant,
    "T": T_full,
    "T_static": T_cool_arr,
    "T_stagnation": T_stagnation_arr,
    "p_static": p_static_arr,
    "p_stagnation": p_stagnation_arr,
    "dQ_dA": dQ_dA_arr,
    "velocity": velocity_arr,
    "h_hot": h_hot_arr,
    "h_hot_enthalpy": h_hot_enthalpy_arr,
    "h_cold": h_cold_arr,
    "T_aw_hot": T_aw_hot_arr,
    "residuals": (global_R, final_R),
    "wall_residual_scaled": wall_residual_scaled,
    "wall_converged": wall_residual_scaled <= RESIDUAL_TOL,
    "film_regime": film_regime,
    "liquid_film": liquid_film_or_none,
    "gaseous_film": gaseous_film_or_none,
}
```

`x` remains an alias for `x_wall`. `T` is defined on `x_wall`; `dQ_dA`,
`qpp_hot`, `h_hot`, `h_hot_enthalpy`, `T_aw_hot`, and `T_drive` are defined on
`x_heat_flux`; coolant temperature, pressure, velocity, phase, quality, and
`h_cold` are defined on `x_coolant`.

`PlotWallTemperature` displays a red `x` at every wall node where
`wall_converged` is false. Pass `mark_nonconverged=False` to suppress these
markers. Their hover labels report the final scaled residual.

---

## 16. Summary

The regenerative cooling solver couples a CEA-based hot-gas property model to a one-dimensional cooling-channel heat-exchanger calculation. The solver walks along the selected cooling circuit, solves a local nonlinear balance for the hot- and cold-side wall temperatures, and then marches coolant temperature and pressure downstream.

At each station, the central balance is:

$$
Q_{hot}=Q_{cond}=Q_{cold}
$$

with:

$$
Q_{hot}=q''_{hot}A_{hot}
$$

$$
Q_{cond}=\frac{T_{hw}-T_{cw}}{R_{wall}}
$$

$$
Q_{cold}=\frac{T_{cw}-T_{cool}}{R_{coolant}}
$$

The hot side is enthalpy-driven, the wall is modeled as a stack of thermal resistances, and the coolant side is temperature-driven. This makes the solver fast enough for design iteration while retaining the main physical couplings needed for regenerative cooling analysis.