diff --git a/examples/pipe_cht/README.md b/examples/pipe_cht/README.md
new file mode 100644
index 0000000..46748ac
--- /dev/null
+++ b/examples/pipe_cht/README.md
@@ -0,0 +1,113 @@
+# 1D pipe flow for conjugate heat transfer (CHT)
+
+`pipe1d.py` models flow and heat transfer in a pipe in 1D (along its axis), as a
+co-simulation FMU meant to be coupled to a solid solver. The pipe wall, or any solid
+around it, is resolved by that solver, typically in 3D. Only the fluid is in the FMU.
+
+![Wall and fluid temperatures from the demo](cht_profile.png)
+
+## Coupling
+
+The pipe is split into `N` equal axial segments. The FMU and the solid solver exchange one
+value per segment, every communication step:
+
+| Direction | Variable | Meaning |
+| --- | --- | --- |
+| solid → FMU | `T_wall[i]` | wetted-wall temperature of segment *i*, area-averaged by the solid solver [K] |
+| FMU → solid | `T_fluid[i]` | bulk fluid temperature of segment *i* [K] |
+| FMU → solid | `h[i]` | heat transfer coefficient of segment *i* [W/(m²·K)] |
+
+On the wetted surface of segment *i*, the solid solver applies a convective (Robin) boundary
+condition: `q'' = h[i] · (T_surface − T_fluid[i])`. This is the usual stable way to couple
+CHT. If your solver needs a heat flux instead, `q_wall[i]` (W/m², wall → fluid) and
+`Q_wall[i]` (W per segment) are also outputs.
+
+One coupling step (what `cht_demo.py` does):
+
+1. solid → FMU: set `T_wall`.
+2. `doStep(t, dt)` on the FMU.
+3. FMU → solid: get `T_fluid` and `h`.
+4. Advance the solid by `dt` with those boundary conditions.
+
+## Parameters
+
+| Name | Default | Unit | Meaning |
+| --- | --- | --- | --- |
+| `N` | 20 | – | number of axial segments (structural parameter, `UInt64`) |
+| `length` | 2.0 | m | pipe length |
+| `diameter` | 0.02 | m | inner (hydraulic) diameter |
+| `rho` | 992.2 | kg/m³ | fluid density |
+| `cp` | 4179 | J/(kg·K) | fluid specific heat |
+| `k_fluid` | 0.631 | W/(m·K) | fluid thermal conductivity |
+| `mu` | 6.53e-4 | Pa·s | fluid dynamic viscosity |
+| `T_initial` | 293.15 | K | initial fluid temperature |
+
+The fluid properties are constant; the defaults are water at about 40 °C. The parameters
+are fixed, so set them before `fmi3ExitInitializationMode`.
+
+`N` sizes every per-segment array (each has `<Dimension valueReference="…"/>` pointing at
+`N` in the model description). Set it in Configuration Mode, before initialisation:
+
+```
+fmi3InstantiateCoSimulation(...)
+fmi3EnterConfigurationMode(instance)
+fmi3SetUInt64(instance, {vr of N}, 1, {50}, 1)
+fmi3ExitConfigurationMode(instance)
+fmi3EnterInitializationMode(...)       # T_wall, T_fluid, h, ... now have 50 elements
+```
+
+In the Python class, `N` is a property whose setter reallocates the arrays. That's the
+pattern to use for any PythonFMU3 model whose array sizes depend on a structural
+parameter, since arrays are not resized automatically. Importers that don't support
+Configuration Mode use the default from the model description (20).
+
+## Inputs and outputs
+
+| Name | Causality | Unit | Meaning |
+| --- | --- | --- | --- |
+| `m_dot` | input | kg/s | mass flow; positive from x = 0 to x = L. Negative values (reverse flow) are supported |
+| `T_in` | input | K | temperature of the fluid entering at the upstream end |
+| `T_wall[N]` | input | K | wall temperature per segment |
+| `T_fluid[N]` | output | K | bulk fluid temperature per segment |
+| `h[N]` | output | W/(m²·K) | heat transfer coefficient per segment |
+| `q_wall[N]`, `Q_wall[N]` | output | W/m², W | heat from the wall into the fluid, per segment |
+| `T_out` | output | K | outlet temperature |
+| `Q_total` | output | W | total heat from the wall into the fluid |
+| `Re`, `Nu` | output | – | Reynolds and Nusselt numbers |
+| `dp` | output | Pa | frictional pressure drop |
+| `x[N]` | calculated parameter | m | segment centres, for mapping segments onto the solid mesh |
+
+## Model
+
+- **Energy:** finite volume, first-order upwind, implicit Euler:
+  `ρ·cp·A·Δx·dTᵢ/dt = ṁ·cp·(T_upstream − Tᵢ) + hᵢ·P·Δx·(T_wall,ᵢ − Tᵢ)`.
+  The system is solved in one O(N) sweep along the flow and is stable for any
+  communication step, so the coupling step can follow the solid solver's time step.
+- **Heat transfer:** fully developed flow. Nu = 3.66 when laminar (Re < 2300),
+  Gnielinski when turbulent (Re > 4000), blended linearly in between.
+- **Friction:** f = 64/Re when laminar, Petukhov when turbulent.
+- **Not modelled:** entrance effects, property variation with temperature, axial
+  conduction in the fluid, and phase change.
+
+## Running it
+
+```bash
+pip install pythonfmu3 numpy
+pythonfmu3 build -f pipe1d.py          # creates Pipe1D.fmu
+
+python cht_demo.py                      # Python class + simple solid model
+python cht_demo.py --segments 50        # finer axial resolution
+python cht_demo.py --fmu Pipe1D.fmu     # same coupling through the FMU (needs FMPy)
+python cht_demo.py --plot profile.png   # save the temperature profile (needs matplotlib)
+```
+
+The demo heats the middle third of a steel tube (outer diameter 30 mm) with 50 kW/m² on
+its outer surface, with water at 0.1 kg/s (Re ≈ 9700). The solid is a small axisymmetric
+conduction model standing in for your 3D solver. At steady state, the water's enthalpy
+rise matches the heater power.
+
+To use your own solver instead, replace the `SolidWall` class: area-average your solver's
+wetted-wall temperature over each axial segment (use `x` and `length` to map your mesh
+faces to segments) to get `T_wall`, and apply each segment's `h` and `T_fluid` as that
+segment's convective boundary condition. The model is unchanged when the pipe sits
+inside a larger solid, or is one of several channels.
diff --git a/examples/pipe_cht/cht_demo.py b/examples/pipe_cht/cht_demo.py
new file mode 100644
index 0000000..c46cf90
--- /dev/null
+++ b/examples/pipe_cht/cht_demo.py
@@ -0,0 +1,193 @@
+"""Conjugate heat transfer demo: the 1D pipe model coupled to a solid wall model.
+
+A heated, thick-walled steel tube: the heater applies a uniform heat flux on the outer
+surface of the middle third of the tube, the rest of the outer surface is insulated, and
+water flows inside. The solid stands in for your 3D solver: it is an axisymmetric (x, r)
+conduction model with the same N axial segments as the pipe, and M radial cells.
+
+Coupling (Robin / "h and T_fluid" exchange), every communication step dt:
+    1. solid -> pipe:  T_wall[i]   (inner-surface temperature of segment i)
+    2. pipe.do_step(t, dt)
+    3. pipe -> solid:  T_fluid[i], h[i]
+    4. solid advances dt with q'' = h[i] * (T_surface - T_fluid[i]) on the inner surface
+
+Run:  python cht_demo.py            (uses the Pipe1D class directly)
+      python cht_demo.py --segments 50
+      python cht_demo.py --fmu Pipe1D.fmu   (uses the built FMU through FMPy)
+
+With a 3D solver, replace SolidWall by your solver: area-average its wetted-wall
+temperature over each axial segment for T_wall, and apply each segment's h and T_fluid as
+its convective boundary condition.
+"""
+import argparse
+import math
+import sys
+from pathlib import Path
+
+import numpy as np
+
+sys.path.insert(0, str(Path(__file__).parent))
+from pipe1d import Pipe1D  # noqa: E402
+
+
+class SolidWall:
+    """Axisymmetric conduction in a tube wall (explicit finite volume, sub-stepped)."""
+
+    def __init__(self, length, r_inner, r_outer, n_axial, n_radial=6,
+                 k=16.0, rho=7900.0, cp=500.0, T0=293.15):
+        self.k, self.rho_cp = k, rho * cp
+        self.n, self.m = n_axial, n_radial
+        self.dx = length / n_axial
+        self.r_faces = np.linspace(r_inner, r_outer, n_radial + 1)
+        self.r = 0.5 * (self.r_faces[1:] + self.r_faces[:-1])
+        self.dr = self.r_faces[1] - self.r_faces[0]
+        self.T = np.full((n_axial, n_radial), T0)
+        self.q_outer = np.zeros(n_axial)   # W/m2 into the solid at r_outer
+        self.volume = math.pi * (self.r_faces[1:] ** 2 - self.r_faces[:-1] ** 2) * self.dx   # (m,)
+        alpha = k / self.rho_cp
+        self.dt_max = 0.2 * min(self.dx, self.dr) ** 2 / alpha
+
+    def wall_temperature(self, h, T_fluid):
+        """Inner-surface temperature, from a half-cell conduction / convection balance."""
+        g = self.k / (0.5 * self.dr)
+        return (g * self.T[:, 0] + h * T_fluid) / (g + h)
+
+    def advance(self, dt, h, T_fluid):
+        steps = max(1, math.ceil(dt / self.dt_max))
+        sub = dt / steps
+        A_ax = math.pi * (self.r_faces[1:] ** 2 - self.r_faces[:-1] ** 2)   # axial face areas (m,)
+        A_rad = 2.0 * math.pi * self.r_faces * self.dx                        # radial face areas (m+1,)
+        for _ in range(steps):
+            Q = np.zeros_like(self.T)
+            # radial conduction between cells
+            flux_r = self.k * A_rad[1:-1] * (self.T[:, 1:] - self.T[:, :-1]) / self.dr
+            Q[:, :-1] += flux_r
+            Q[:, 1:] -= flux_r
+            # axial conduction between segments (insulated ends)
+            flux_x = self.k * A_ax * (self.T[1:, :] - self.T[:-1, :]) / self.dx
+            Q[:-1, :] += flux_x
+            Q[1:, :] -= flux_x
+            # inner surface: convection to the fluid (series with the half-cell)
+            T_s = self.wall_temperature(h, T_fluid)
+            Q[:, 0] -= h * A_rad[0] * (T_s - T_fluid)
+            # outer surface: heater
+            Q[:, -1] += self.q_outer * A_rad[-1]
+            self.T += sub * Q / (self.rho_cp * self.volume)
+
+    def energy(self):
+        return float(np.sum(self.rho_cp * self.volume * self.T))
+
+
+class PythonPipe:
+    """Drives the Pipe1D class directly, through the same calls an importer would make."""
+
+    def __init__(self, segments, **parameters):
+        self.model = Pipe1D(instance_name="pipe")
+        self.model.N = segments          # structural parameter: resizes the arrays
+        for name, value in parameters.items():
+            setattr(self.model, name, value)
+        self.model.enter_initialization_mode()
+        self.model.exit_initialization_mode()
+
+    def set(self, name, value):
+        setattr(self.model, name, value)
+
+    def get(self, name):
+        return np.atleast_1d(np.asarray(getattr(self.model, name), dtype=float))
+
+    def do_step(self, t, dt):
+        self.model.do_step(t, dt)
+
+
+class FmuPipe:
+    """Drives the built FMU through FMPy (pip install fmpy)."""
+
+    def __init__(self, fmu_path, segments, **parameters):
+        from fmpy import extract, read_model_description
+        from fmpy.fmi3 import FMU3Slave
+
+        md = read_model_description(fmu_path)
+        self.vr = {v.name: v.valueReference for v in md.modelVariables}
+        self.size = {v.name: (segments if v.dimensions else 1) for v in md.modelVariables}
+        self.fmu = FMU3Slave(guid=md.guid, unzipDirectory=extract(fmu_path),
+                             modelIdentifier=md.coSimulation.modelIdentifier, instanceName="pipe")
+        self.fmu.instantiate()
+        # The structural parameter N is set in Configuration Mode
+        self.fmu.enterConfigurationMode()
+        self.fmu.setUInt64([self.vr["N"]], [segments])
+        self.fmu.exitConfigurationMode()
+        self.fmu.enterInitializationMode()
+        for name, value in parameters.items():
+            self.set(name, value)
+        self.fmu.exitInitializationMode()
+
+    def set(self, name, value):
+        values = np.atleast_1d(np.asarray(value, dtype=float))
+        self.fmu.setFloat64([self.vr[name]], list(values))
+
+    def get(self, name):
+        return np.asarray(self.fmu.getFloat64([self.vr[name]], nValues=self.size[name]), dtype=float)
+
+    def do_step(self, t, dt):
+        self.fmu.doStep(currentCommunicationPoint=t, communicationStepSize=dt)
+
+
+def main():
+    parser = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
+    parser.add_argument("--fmu", help="use this built Pipe1D FMU (via FMPy) instead of the Python class")
+    parser.add_argument("--end", type=float, default=1200.0, help="end time [s]")
+    parser.add_argument("--dt", type=float, default=1.0, help="communication step [s]")
+    parser.add_argument("--segments", type=int, default=20, help="number of axial segments N")
+    parser.add_argument("--plot", help="save a plot to this file (needs matplotlib)")
+    args = parser.parse_args()
+
+    length, d_inner, d_outer = 2.0, 0.02, 0.03
+    m_dot, T_in, heat_flux = 0.1, 293.15, 50e3      # kg/s, K, W/m2 on the heated outer surface
+
+    n = args.segments
+    params = dict(length=length, diameter=d_inner, T_initial=T_in)
+    pipe = FmuPipe(args.fmu, n, **params) if args.fmu else PythonPipe(n, **params)
+    pipe.set("m_dot", m_dot)
+    pipe.set("T_in", T_in)
+
+    solid = SolidWall(length, d_inner / 2, d_outer / 2, n, T0=T_in)
+    third = n // 3
+    solid.q_outer[third:2 * third] = heat_flux
+    heater_power = heat_flux * math.pi * d_outer * (length / n) * third
+
+    h, T_fluid = pipe.get("h"), pipe.get("T_fluid")
+    t = 0.0
+    print(f"{'t [s]':>7} {'T_out [K]':>10} {'max T_wall [K]':>15} {'Q to fluid [W]':>15}")
+    while t < args.end - 1e-9:
+        pipe.set("T_wall", solid.wall_temperature(h, T_fluid))
+        pipe.do_step(t, args.dt)
+        h, T_fluid = pipe.get("h"), pipe.get("T_fluid")
+        solid.advance(args.dt, h, T_fluid)
+        t += args.dt
+        if abs(t % 120.0) < 1e-6 or t >= args.end - 1e-9:
+            T_wall = solid.wall_temperature(h, T_fluid)
+            print(f"{t:7.0f} {pipe.get('T_out')[0]:10.2f} {T_wall.max():15.2f} {pipe.get('Q_total')[0]:15.1f}")
+
+    absorbed = m_dot * 4179.0 * (pipe.get("T_out")[0] - T_in)
+    print(f"\nheater power {heater_power:.1f} W, enthalpy rise of the water {absorbed:.1f} W "
+          f"({100 * absorbed / heater_power:.1f}% at t = {t:.0f} s)")
+    print(f"Re = {pipe.get('Re')[0]:.0f}, Nu = {pipe.get('Nu')[0]:.1f}, h = {h[0]:.0f} W/(m2.K), "
+          f"dp = {pipe.get('dp')[0]:.1f} Pa")
+
+    if args.plot:
+        import matplotlib.pyplot as plt
+        x = pipe.get("x")
+        plt.plot(x, solid.wall_temperature(h, T_fluid) - 273.15, label="wall (inner surface)")
+        plt.plot(x, T_fluid - 273.15, label="fluid (bulk)")
+        plt.axvspan(x[third] - length / n / 2, x[2 * third - 1] + length / n / 2,
+                    alpha=0.15, label="heated section")
+        plt.xlabel("x [m]")
+        plt.ylabel("temperature [°C]")
+        plt.legend()
+        plt.title(f"Pipe CHT at t = {t:.0f} s")
+        plt.savefig(args.plot, dpi=150, bbox_inches="tight")
+        print(f"plot saved to {args.plot}")
+
+
+if __name__ == "__main__":
+    main()
diff --git a/examples/pipe_cht/cht_profile.png b/examples/pipe_cht/cht_profile.png
new file mode 100644
index 0000000..261943e
Binary files /dev/null and b/examples/pipe_cht/cht_profile.png differ
diff --git a/examples/pipe_cht/pipe1d.py b/examples/pipe_cht/pipe1d.py
new file mode 100644
index 0000000..5a94142
--- /dev/null
+++ b/examples/pipe_cht/pipe1d.py
@@ -0,0 +1,249 @@
+"""1D pipe flow with wall heat transfer, for conjugate heat transfer (CHT) co-simulation.
+
+The fluid is resolved along the pipe axis only (N segments). The pipe wall belongs to a
+separate solid solver (e.g. a 3D conduction model) that exchanges, per axial segment:
+
+    solid -> FMU:  T_wall[i]   area-averaged wall temperature of segment i   [K]
+    FMU -> solid:  T_fluid[i]  bulk fluid temperature of segment i           [K]
+                   h[i]        heat transfer coefficient of segment i        [W/(m2.K)]
+
+so the solid applies a convective (Robin) condition  q'' = h[i] * (T_solid - T_fluid[i])
+on the wetted surface of segment i. This Robin coupling is the usual stable choice for CHT;
+q_wall[i] (heat flux into the fluid) is also provided for a Neumann-type coupling.
+
+Model, per segment (finite volume, first-order upwind, implicit Euler):
+
+    rho*cp*A*dx * dT_i/dt = m_dot*cp*(T_up - T_i) + h_i*P*dx*(T_wall_i - T_i)
+
+The implicit upwind system is (bi)diagonal, so each step is one O(N) sweep and is stable for
+any communication step size. Reverse flow (m_dot < 0) is supported: T_in then enters at x=L.
+
+Heat transfer: fully developed Nusselt number, Nu = 3.66 (laminar, Re < 2300) and
+Gnielinski (Re > 4000), blended linearly in between. Friction: f = 64/Re (laminar) and
+Petukhov (turbulent). Constant fluid properties (defaults: water at about 40 C).
+
+The number of segments N is a structural parameter: set it in Configuration Mode
+(fmi3EnterConfigurationMode, fmi3SetUInt64, fmi3ExitConfigurationMode) before
+initialisation, and every per-segment array is resized to N.
+
+Build:  pythonfmu3 build -f pipe1d.py
+"""
+import math
+
+import numpy as np
+
+from pythonfmu3 import (
+    Dimension,
+    Fmi3Causality,
+    Fmi3Initial,
+    Fmi3Slave,
+    Fmi3Variability,
+    Float64,
+    UInt64,
+    Unit,
+)
+
+# Default number of axial segments (the structural parameter N can override it).
+N_SEGMENTS = 20
+
+RE_LAMINAR = 2300.0
+RE_TURBULENT = 4000.0
+
+
+def friction_factor(re: float) -> float:
+    """Darcy friction factor of a smooth pipe."""
+    if re <= 0.0:
+        return 0.0
+    laminar = 64.0 / max(re, 1e-12)
+    if re <= RE_LAMINAR:
+        return laminar
+    turbulent = (0.790 * math.log(max(re, RE_TURBULENT)) - 1.64) ** -2  # Petukhov
+    if re >= RE_TURBULENT:
+        return turbulent
+    w = (re - RE_LAMINAR) / (RE_TURBULENT - RE_LAMINAR)
+    return (1.0 - w) * laminar + w * turbulent
+
+
+def nusselt(re: float, pr: float) -> float:
+    """Fully developed Nusselt number: 3.66 laminar, Gnielinski turbulent, blended in between."""
+    laminar = 3.66
+    if re <= RE_LAMINAR:
+        return laminar
+    re_t = max(re, RE_TURBULENT)
+    f = (0.790 * math.log(re_t) - 1.64) ** -2
+    turbulent = (f / 8.0) * (re_t - 1000.0) * pr / (1.0 + 12.7 * math.sqrt(f / 8.0) * (pr ** (2.0 / 3.0) - 1.0))
+    if re >= RE_TURBULENT:
+        return turbulent
+    w = (re - RE_LAMINAR) / (RE_TURBULENT - RE_LAMINAR)
+    return (1.0 - w) * laminar + w * turbulent
+
+
+class Pipe1D(Fmi3Slave):
+
+    def __init__(self, **kwargs):
+        super().__init__(**kwargs)
+
+        self.author = "PythonFMU3 examples"
+        self.description = "1D pipe flow with wall heat transfer for CHT coupling with a solid solver"
+
+        self.time = 0.0
+        self._n = N_SEGMENTS
+
+        # Parameters
+        self.length = 2.0          # m
+        self.diameter = 0.02       # m
+        self.rho = 992.2           # kg/m3
+        self.cp = 4179.0           # J/(kg.K)
+        self.k_fluid = 0.631       # W/(m.K)
+        self.mu = 6.53e-4          # Pa.s
+        self.T_initial = 293.15    # K, initial fluid temperature
+
+        # Inputs (T_wall is allocated in _allocate)
+        self.m_dot = 0.05                       # kg/s (positive: from x=0 to x=L)
+        self.T_in = 293.15                      # K, temperature of the entering fluid
+
+        # Per-segment arrays: T_wall (input, K), T_fluid (K), h (W/(m2.K)),
+        # q_wall (W/m2, wall -> fluid), Q_wall (W per segment), x (m, segment centres)
+        self._allocate(self._n)
+
+        # Scalar outputs
+        self.T_out = self.T_initial                 # K
+        self.Q_total = 0.0                          # W, wall -> fluid over the whole pipe
+        self.Re = 0.0
+        self.Nu = 0.0
+        self.dp = 0.0                               # Pa, frictional pressure drop
+
+        self.register_units([
+            Unit(name="K", K=1),
+            Unit(name="m", m=1),
+            Unit(name="kg/s", kg=1, s=-1),
+            Unit(name="kg/m3", kg=1, m=-3),
+            Unit(name="J/(kg.K)", m=2, s=-2, K=-1),
+            Unit(name="W/(m.K)", kg=1, m=1, s=-3, K=-1),
+            Unit(name="W/(m2.K)", kg=1, s=-3, K=-1),
+            Unit(name="W/m2", kg=1, s=-3),
+            Unit(name="W", kg=1, m=2, s=-3),
+            Unit(name="Pa", kg=1, m=-1, s=-2),
+            Unit(name="Pa.s", kg=1, m=-1, s=-1),
+        ])
+
+        parameter = dict(causality=Fmi3Causality.parameter, variability=Fmi3Variability.fixed)
+        output = dict(causality=Fmi3Causality.output, variability=Fmi3Variability.continuous,
+                      initial=Fmi3Initial.calculated)
+
+        self.register_variable(Float64("time", causality=Fmi3Causality.independent, variability=Fmi3Variability.continuous))
+
+        # Structural parameter: the arrays below have N elements (setting N resizes them)
+        n_var = UInt64("N", causality=Fmi3Causality.structuralParameter, variability=Fmi3Variability.fixed,
+                       start=self._n, description="number of axial segments")
+        self.register_variable(n_var)
+        segments = [Dimension(valueReference=str(n_var.value_reference))]
+
+        self.register_variable(Float64("length", start=self.length, unit="m", **parameter))
+        self.register_variable(Float64("diameter", start=self.diameter, unit="m", **parameter))
+        self.register_variable(Float64("rho", start=self.rho, unit="kg/m3", **parameter))
+        self.register_variable(Float64("cp", start=self.cp, unit="J/(kg.K)", **parameter))
+        self.register_variable(Float64("k_fluid", start=self.k_fluid, unit="W/(m.K)", **parameter))
+        self.register_variable(Float64("mu", start=self.mu, unit="Pa.s", **parameter))
+        self.register_variable(Float64("T_initial", start=self.T_initial, unit="K", **parameter))
+
+        self.register_variable(Float64("m_dot", causality=Fmi3Causality.input, variability=Fmi3Variability.continuous,
+                                       start=self.m_dot, unit="kg/s"))
+        self.register_variable(Float64("T_in", causality=Fmi3Causality.input, variability=Fmi3Variability.continuous,
+                                       start=self.T_in, unit="K"))
+        self.register_variable(Float64("T_wall", causality=Fmi3Causality.input, variability=Fmi3Variability.continuous,
+                                       unit="K", dimensions=segments))
+
+        self.register_variable(Float64("T_fluid", unit="K", dimensions=segments, **output))
+        self.register_variable(Float64("h", unit="W/(m2.K)", dimensions=segments, **output))
+        self.register_variable(Float64("q_wall", unit="W/m2", dimensions=segments, **output))
+        self.register_variable(Float64("Q_wall", unit="W", dimensions=segments, **output))
+        self.register_variable(Float64("x", causality=Fmi3Causality.calculatedParameter, variability=Fmi3Variability.fixed,
+                                       initial=Fmi3Initial.calculated, unit="m", dimensions=segments))
+        self.register_variable(Float64("T_out", unit="K", **output))
+        self.register_variable(Float64("Q_total", unit="W", **output))
+        self.register_variable(Float64("Re", **output))
+        self.register_variable(Float64("Nu", **output))
+        self.register_variable(Float64("dp", unit="Pa", **output))
+
+    # --- structural parameter ----------------------------------------------------------
+
+    @property
+    def N(self) -> int:
+        return self._n
+
+    @N.setter
+    def N(self, value):
+        value = int(value)
+        if value < 1:
+            raise ValueError(f"N must be at least 1, got {value}")
+        if value != self._n or len(self.T_fluid) != value:
+            self._allocate(value)
+
+    def _allocate(self, n: int):
+        """(Re)allocate every per-segment array for n segments."""
+        self._n = n
+        self.T_wall = np.full(n, self.T_initial)
+        self.T_fluid = np.full(n, self.T_initial)
+        self.h = np.zeros(n)
+        self.q_wall = np.zeros(n)
+        self.Q_wall = np.zeros(n)
+        self.x = (np.arange(n) + 0.5) * self.length / n
+
+    # --- FMI callbacks -----------------------------------------------------------------
+
+    def exit_initialization_mode(self):
+        n = self._n
+        self.T_fluid = np.full(n, self.T_initial)
+        self.x = (np.arange(n) + 0.5) * self.length / n
+        self._update_heat_transfer()
+        self._update_outputs()
+
+    def do_step(self, current_time, step_size):
+        self._update_heat_transfer()
+        self.T_fluid = self._advance(self.T_fluid, step_size)
+        self._update_outputs()
+        self.time = current_time + step_size
+        return True
+
+    # --- model -------------------------------------------------------------------------
+
+    @property
+    def _area(self):
+        return math.pi * self.diameter ** 2 / 4.0
+
+    @property
+    def _perimeter(self):
+        return math.pi * self.diameter
+
+    def _update_heat_transfer(self):
+        self.Re = abs(self.m_dot) * self.diameter / (self._area * self.mu)
+        pr = self.mu * self.cp / self.k_fluid
+        self.Nu = nusselt(self.Re, pr)
+        self.h = np.full(self._n, self.Nu * self.k_fluid / self.diameter)
+        velocity = abs(self.m_dot) / (self.rho * self._area)
+        self.dp = friction_factor(self.Re) * self.length / self.diameter * 0.5 * self.rho * velocity ** 2
+
+    def _advance(self, T_old, dt):
+        """One implicit-Euler step of the upwind finite-volume energy equation."""
+        n = self._n
+        dx = self.length / n
+        T_wall = np.asarray(self.T_wall, dtype=float).reshape(n)
+        storage = self.rho * self.cp * self._area * dx / dt     # W/K
+        convection = abs(self.m_dot) * self.cp                  # W/K
+        wall = self.h * self._perimeter * dx                    # W/K per segment
+
+        order = range(n) if self.m_dot >= 0.0 else range(n - 1, -1, -1)
+        T_new = np.empty(n)
+        T_up = self.T_in
+        for i in order:
+            T_new[i] = (storage * T_old[i] + convection * T_up + wall[i] * T_wall[i]) / (storage + convection + wall[i])
+            T_up = T_new[i]
+        return T_new
+
+    def _update_outputs(self):
+        T_wall = np.asarray(self.T_wall, dtype=float).reshape(self._n)
+        self.q_wall = self.h * (T_wall - self.T_fluid)
+        self.Q_wall = self.q_wall * self._perimeter * self.length / self._n
+        self.Q_total = float(np.sum(self.Q_wall))
+        self.T_out = float(self.T_fluid[-1] if self.m_dot >= 0.0 else self.T_fluid[0])
diff --git a/pythonfmu3/tests/test_example_pipe_cht.py b/pythonfmu3/tests/test_example_pipe_cht.py
new file mode 100644
index 0000000..e650f4a
--- /dev/null
+++ b/pythonfmu3/tests/test_example_pipe_cht.py
@@ -0,0 +1,198 @@
+"""Tests for examples/pipe_cht (only available in a source checkout)."""
+import math
+import subprocess
+import sys
+import textwrap
+import zipfile
+from pathlib import Path
+from xml.etree.ElementTree import fromstring
+
+import pytest
+
+np = pytest.importorskip("numpy")
+
+EXAMPLE = Path(__file__).parent.parent.parent / "examples" / "pipe_cht"
+if not EXAMPLE.exists():
+    pytest.skip("examples are only available in a source checkout", allow_module_level=True)
+
+sys.path.insert(0, str(EXAMPLE))
+import pipe1d  # noqa: E402
+from pipe1d import N_SEGMENTS, Pipe1D  # noqa: E402
+
+
+def make_pipe(n=N_SEGMENTS, **inputs):
+    pipe = Pipe1D(instance_name="test")
+    pipe.N = n
+    pipe.enter_initialization_mode()
+    pipe.exit_initialization_mode()
+    for name, value in inputs.items():
+        setattr(pipe, name, value)
+    return pipe
+
+
+def run_to_steady_state(pipe, steps=50, dt=1e4):
+    t = 0.0
+    for _ in range(steps):
+        pipe.do_step(t, dt)
+        t += dt
+    return pipe
+
+
+@pytest.mark.parametrize("m_dot", [0.005, 0.1])   # laminar and turbulent
+def test_constant_wall_temperature_matches_analytic(m_dot):
+    T_in, T_w = 293.15, 353.15
+    pipe = run_to_steady_state(make_pipe(m_dot=m_dot, T_in=T_in, T_wall=np.full(N_SEGMENTS, T_w)))
+
+    # Exact steady solution of the continuous problem: T_w - T = (T_w - T_in) exp(-h P x / (m_dot cp))
+    ntu_per_length = pipe.h[0] * math.pi * pipe.diameter / (m_dot * pipe.cp)
+    expected = T_w - (T_w - T_in) * np.exp(-ntu_per_length * (pipe.x + pipe.length / (2 * pipe.N)))
+    # first-order upwind: compare the outlet (end of the last segment) within 2% of the temperature rise
+    assert abs(pipe.T_out - expected[-1]) < 0.02 * (T_w - T_in)
+    assert np.all(np.diff(pipe.T_fluid) > 0)
+
+
+def test_resolution_converges_to_analytic():
+    """First-order upwind: the outlet error shrinks as N grows."""
+    T_in, T_w, m_dot = 293.15, 353.15, 0.005
+    errors = []
+    for n in (10, 40, 160):
+        pipe = run_to_steady_state(make_pipe(n=n, m_dot=m_dot, T_in=T_in, T_wall=np.full(n, T_w)))
+        ntu = pipe.h[0] * math.pi * pipe.diameter * pipe.length / (m_dot * pipe.cp)
+        errors.append(abs(pipe.T_out - (T_w - (T_w - T_in) * math.exp(-ntu))))
+    assert errors[0] > errors[1] > errors[2]
+    assert errors[2] < 0.25 * errors[0]
+
+
+def test_setting_N_resizes_arrays_through_the_fmi_api():
+    pipe = Pipe1D(instance_name="test")
+    n_vr = next(v.value_reference for v in pipe.vars.values() if v.name == "N")
+    wall_vr = next(v.value_reference for v in pipe.vars.values() if v.name == "T_wall")
+    fluid_vr = next(v.value_reference for v in pipe.vars.values() if v.name == "T_fluid")
+    pipe.set_uint64([n_vr], [35])
+    pipe.set_float64([wall_vr], list(np.linspace(300.0, 340.0, 35)))
+    pipe.enter_initialization_mode()
+    pipe.exit_initialization_mode()
+    pipe.do_step(0.0, 1.0)
+    assert len(pipe.get_float64([fluid_vr])) == 35
+    assert pipe.T_wall.shape == (35,)
+    with pytest.raises(ValueError):
+        pipe.N = 0
+
+
+def test_energy_is_conserved_at_steady_state():
+    T_wall = np.linspace(300.0, 380.0, N_SEGMENTS)
+    pipe = run_to_steady_state(make_pipe(m_dot=0.03, T_in=290.0, T_wall=T_wall))
+    enthalpy_rise = pipe.m_dot * pipe.cp * (pipe.T_out - pipe.T_in)
+    assert pipe.Q_total == pytest.approx(enthalpy_rise, rel=1e-9)
+
+
+def test_reverse_flow_mirrors_forward_flow():
+    T_wall = np.linspace(300.0, 380.0, N_SEGMENTS)
+    forward = run_to_steady_state(make_pipe(m_dot=0.03, T_in=290.0, T_wall=T_wall))
+    reverse = run_to_steady_state(make_pipe(m_dot=-0.03, T_in=290.0, T_wall=T_wall[::-1].copy()))
+    np.testing.assert_allclose(reverse.T_fluid[::-1], forward.T_fluid, rtol=1e-12)
+    assert reverse.T_out == pytest.approx(forward.T_out)
+
+
+@pytest.mark.parametrize("re, nu", [(1000, 3.66), (9749, 64.6)])
+def test_nusselt_correlation(re, nu):
+    assert pipe1d.nusselt(re, 4.33) == pytest.approx(nu, rel=0.01)
+
+
+def test_cht_demo_balances_heater_power():
+    demo = subprocess.run(
+        [sys.executable, str(EXAMPLE / "cht_demo.py"), "--end", "1800"],
+        capture_output=True, text=True, timeout=300,
+    )
+    assert demo.returncode == 0, demo.stderr
+    assert "(100.0% at t = 1800 s)" in demo.stdout or "(99.9% at t = 1800 s)" in demo.stdout, demo.stdout
+
+
+def test_build_fmu(tmp_path):
+    from pythonfmu3 import FmuBuilder
+
+    fmu = FmuBuilder.build_FMU(EXAMPLE / "pipe1d.py", dest=tmp_path)
+    with zipfile.ZipFile(fmu) as z:
+        md = fromstring(z.read("modelDescription.xml"))
+    variables = {v.get("name"): v for v in md.find("ModelVariables")}
+    n = variables["N"]
+    assert n.tag == "UInt64" and n.get("causality") == "structuralParameter" and n.get("start") == str(N_SEGMENTS)
+    for name in ("T_wall", "T_fluid", "h", "x"):
+        dims = variables[name].findall("Dimension")
+        assert [d.get("valueReference") for d in dims] == [n.get("valueReference")]
+    assert variables["T_wall"].get("causality") == "input"
+    assert variables["h"].get("causality") == "output"
+
+
+def test_fmu_matches_python_model(tmp_path):
+    """Run the built FMU through the FMI C API (ctypes) and compare with the Python class."""
+    from pythonfmu3 import FmuBuilder
+    from pythonfmu3.builder import FmuBuilder as _B
+
+    if not _B.has_binary():
+        pytest.skip("needs the compiled wrapper")
+    fmu = FmuBuilder.build_FMU(EXAMPLE / "pipe1d.py", dest=tmp_path)
+    target = tmp_path / "pipe"
+    with zipfile.ZipFile(fmu) as z:
+        z.extractall(target)
+        md = fromstring(z.read("modelDescription.xml"))
+    vr = {v.get("name"): int(v.get("valueReference")) for v in md.find("ModelVariables")}
+
+    n = 35   # not the default, set through Configuration Mode
+    T_wall = np.linspace(300.0, 380.0, n)
+    script = textwrap.dedent(f"""
+        import ctypes as C, json, sys
+        from pathlib import Path
+        root = Path({str(target)!r})
+        lib = C.CDLL(str(next(p for p in (root / "binaries").rglob("*.*")
+                              if p.suffix in (".so", ".dylib", ".dll") and "backend" not in p.name)))
+        lib.fmi3InstantiateCoSimulation.restype = C.c_void_p
+        lib.fmi3InstantiateCoSimulation.argtypes = [C.c_char_p, C.c_char_p, C.c_char_p, C.c_bool, C.c_bool, C.c_bool,
+            C.c_bool, C.c_void_p, C.c_size_t, C.c_void_p, C.c_void_p, C.c_void_p]
+        for f in ("fmi3EnterInitializationMode",):
+            getattr(lib, f).argtypes = [C.c_void_p, C.c_bool, C.c_double, C.c_double, C.c_bool, C.c_double]
+        for f in ("fmi3ExitInitializationMode", "fmi3FreeInstance", "fmi3EnterConfigurationMode", "fmi3ExitConfigurationMode"):
+            getattr(lib, f).argtypes = [C.c_void_p]
+        lib.fmi3SetUInt64.argtypes = [C.c_void_p, C.POINTER(C.c_uint32), C.c_size_t, C.POINTER(C.c_uint64), C.c_size_t]
+        lib.fmi3DoStep.argtypes = [C.c_void_p, C.c_double, C.c_double, C.c_bool, C.c_void_p, C.c_void_p, C.c_void_p, C.c_void_p]
+        for f in ("fmi3SetFloat64", "fmi3GetFloat64"):
+            getattr(lib, f).argtypes = [C.c_void_p, C.POINTER(C.c_uint32), C.c_size_t, C.POINTER(C.c_double), C.c_size_t]
+        def vrs(v): return (C.c_uint32 * 1)(v)
+        c = lib.fmi3InstantiateCoSimulation(b"pipe", b"", str(root / "resources").encode(), False, False, False, False,
+                                            None, 0, None, None, None)
+        assert c
+        assert lib.fmi3EnterConfigurationMode(c) == 0
+        assert lib.fmi3SetUInt64(c, vrs({vr}["N"]), 1, (C.c_uint64 * 1)({n}), 1) == 0
+        assert lib.fmi3ExitConfigurationMode(c) == 0
+        assert lib.fmi3EnterInitializationMode(c, False, 0, 0, False, 0) == 0
+        assert lib.fmi3ExitInitializationMode(c) == 0
+        def set_(name, values):
+            arr = (C.c_double * len(values))(*values)
+            assert lib.fmi3SetFloat64(c, vrs({vr}[name]), 1, arr, len(values)) == 0
+        def get(name, n):
+            arr = (C.c_double * n)()
+            assert lib.fmi3GetFloat64(c, vrs({vr}[name]), 1, arr, n) == 0
+            return list(arr)
+        set_("m_dot", [0.03]); set_("T_in", [290.0]); set_("T_wall", {T_wall.tolist()})
+        flags = (C.c_bool * 3)(); last = C.c_double()
+        t = 0.0
+        for _ in range(5):
+            assert lib.fmi3DoStep(c, t, 10.0, True, C.byref(flags, 0), C.byref(flags, 1), C.byref(flags, 2), C.byref(last)) == 0
+            t += 10.0
+        print(json.dumps({{"T_fluid": get("T_fluid", {n}), "h": get("h", {n}), "T_out": get("T_out", 1)}}))
+        lib.fmi3FreeInstance(c)
+    """)
+    proc = subprocess.run([sys.executable, "-c", script], capture_output=True, text=True, timeout=120)
+    assert proc.returncode == 0, proc.stdout + proc.stderr
+    import json
+    result = json.loads(proc.stdout.strip().splitlines()[-1])
+
+    reference = make_pipe(n=n, m_dot=0.03, T_in=290.0, T_wall=T_wall)
+    t = 0.0
+    for _ in range(5):
+        reference.do_step(t, 10.0)
+        t += 10.0
+    assert len(result["T_fluid"]) == n
+    np.testing.assert_allclose(result["T_fluid"], reference.T_fluid, rtol=1e-12)
+    np.testing.assert_allclose(result["h"], reference.h, rtol=1e-12)
+    assert result["T_out"][0] == pytest.approx(reference.T_out)