From b95df453aa695fd15460b111808733a9ac29b844 Mon Sep 17 00:00:00 2001 From: Stephen Smith Date: Thu, 12 Nov 2026 21:13:12 +0000 Subject: [PATCH 1/2] Add a 1D pipe flow example for conjugate heat transfer co-simulation examples/pipe_cht/pipe1d.py: co-simulation FMU of a pipe discretised into N axial segments, to be coupled to a (3D) solid solver. Per segment it takes the wall temperature and returns the bulk fluid temperature and heat transfer coefficient (Robin coupling), plus wall heat flux, outlet temperature, total heat, Re, Nu and pressure drop. - implicit upwind finite volume energy equation: one O(N) sweep per step, stable for any communication step; reverse flow supported - Nu = 3.66 laminar / Gnielinski turbulent, Petukhov friction - cht_demo.py couples it to an axisymmetric steel-wall conduction model (heated middle third) through the same set/doStep/get loop a CHT driver uses; also runs the built FMU through FMPy with --fmu - README with the coupling scheme, parameters, inputs/outputs and model Tests: analytic constant-wall-temperature solution (laminar and turbulent), steady-state energy conservation, reverse-flow symmetry, Nusselt values, demo energy balance, model description arrays, and the built FMU run through the FMI C API matching the Python model. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01X8Wmw797DJpgrMFmT6n8CG --- examples/pipe_cht/README.md | 96 ++++++++++ examples/pipe_cht/cht_demo.py | 185 ++++++++++++++++++ examples/pipe_cht/cht_profile.png | Bin 0 -> 51748 bytes examples/pipe_cht/pipe1d.py | 218 ++++++++++++++++++++++ pythonfmu3/tests/test_example_pipe_cht.py | 161 ++++++++++++++++ 5 files changed, 660 insertions(+) create mode 100644 examples/pipe_cht/README.md create mode 100644 examples/pipe_cht/cht_demo.py create mode 100644 examples/pipe_cht/cht_profile.png create mode 100644 examples/pipe_cht/pipe1d.py create mode 100644 pythonfmu3/tests/test_example_pipe_cht.py diff --git a/examples/pipe_cht/README.md b/examples/pipe_cht/README.md new file mode 100644 index 0000000..4d23986 --- /dev/null +++ b/examples/pipe_cht/README.md @@ -0,0 +1,96 @@ +# 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_SEGMENTS` 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 | +| --- | --- | --- | --- | +| `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 | +| `N_SEGMENTS` | 20 | – | number of axial segments: a constant at the top of `pipe1d.py`, so edit it and rebuild the FMU to change it | + +The fluid properties are constant; the defaults are water at about 40 °C. The parameters +are fixed, so set them before `fmi3ExitInitializationMode`. + +## 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 --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..9cb0c7d --- /dev/null +++ b/examples/pipe_cht/cht_demo.py @@ -0,0 +1,185 @@ +"""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 N axial segments matching 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 --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 N_SEGMENTS, 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, **parameters): + self.model = Pipe1D(instance_name="pipe") + 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, **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: (N_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() + 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("--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 + + params = dict(length=length, diameter=d_inner, T_initial=T_in) + pipe = FmuPipe(args.fmu, **params) if args.fmu else PythonPipe(**params) + pipe.set("m_dot", m_dot) + pipe.set("T_in", T_in) + + solid = SolidWall(length, d_inner / 2, d_outer / 2, N_SEGMENTS, T0=T_in) + third = N_SEGMENTS // 3 + solid.q_outer[third:2 * third] = heat_flux + heater_power = heat_flux * math.pi * d_outer * (length / N_SEGMENTS) * 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_SEGMENTS / 2, x[2 * third - 1] + length / N_SEGMENTS / 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 0000000000000000000000000000000000000000..261943e8cd76bcd7c3057186271c89d8205b8a25 GIT binary patch literal 51748 zcmeFYXH-;8w=UX%7(he?L_mTfO3pz*q9jo=O>T0No7_Z+zJQ2`N)nKqb52cEl9F?7 z5ReRwqy~wn+SmQD)m3^%ui+q-hQ@!_7iZN}?jwOvR-HN9Rv|I$H)TxbPN0Y0}g1J9>ho12MPI5;Cr zT&=AUCYEq6M|(>fYkL=06MHlGQyC9-4}_zLwY9kjkFW`!;3Eq-hk$^v2?xKq84t%J zZeBQtu!*Uuu!RX9p9zntoip4HbOZ)#W$kF?XyOb;p<)4n;LZOTx77cB+$LtfALMl7 z{nwyB8u;;8-VdU@BeLj6N$>y(I*P*$I?eJ=Vx z+Wfx`&4d2GHg#~ewzRhAFmtfC08N4LEs^%N_73j$w&rjbGiPf@oOnDep8W2S;Xit0 zaJ8~_VX&~agEN@8f*Ux4Gu#OYcX4HKuwXD@Kmz4qFtImhK$v(kSedxN8O)sFCa!RE z1_x&bgoC-Yg*Ese_t@_)z{_31lDTjqTpTQ1-NE`vS%Ujib9*4rKd*K+v9pur;V8XAxW?oFIsZ{M}BgWFrWT2bA44UB{MpLsm|Jr7

nT3z|2$l47i&v<6IY}&+(YT#H0`%)as59s5Pxgqzqfmf-ZExX_ z{#$n(r2b~gGCYyurlTa`736-zEyVMPpHDzg58OTiw}L#vzi+u^ z{`aq9;ePn{s~G<6RexKm|9GPY*cZSWbL(=GUA**r`1p8#=0SjWi9i~1;iDIkZI|k8 zLh=#*2I)*;eVNs9^jO(1#SO~->a1IAczo{W?x)k8k7Ws5nHM=le0(Hc-{OZRE%`Jl zQZ@)xN2cp=YjIx# zhR+S*rpMheUbxMD%b&bgZT>=lvbjFv;X-R)y!fZjV+LRTUE&0m5N>mPZXOp}@+3U3 z!7k|AUb#83q}K^{7nLpiCh;qGRPZF>JNXiiWhZx6?vL&Nl1zLFZO7+!G%`T%o8*}) zW^qLY<&&s>V_qNQHOhz-JJv&!UWj0Ci4fr}3iVjgIMP&`YmmQ^inDyWBsin^&(A*+ z_(uZ&NZ=m{{3C&XB=C;}{*k~x68J{~|485;3H*O90j!Q>sIln%%RdK88$2KAFrSl) z)o>FmmHP#&B<%)kUwPE8j#vCw|458GR})%!*oSELa=hMC`$w`~DB;T?58ihoYSrX1 zqF^2(uCh9=8h>y%h?0q`qn(4RowX^$AJ~D9lZS_smz}}N)z#5Ogo_J6wtrf1IyhT$ z`G;e0!noq7pbMa4%A2@fJj>)kfO?l-d@3WQ>0!8mCj4kH@~dOp9mYr_gnvKa-GwUw zkk4o|xqBs73a)np9ctRAe zxGeWpQX0MBSu4GEo4Z*qoElDJR}uh+!)1MXl!Vacey8h3n;YeQRT5|05dk*kD3r+Y z;X=e_kI-zzuia&R3yUn#870AqRyQPtxl9Q*A2_=a80POd`rm`*&1KR8-1o>GxG{b9&&^5p`JkiEe$ zVP{v@=7eQ(|ZFpa&HJS{S+uzJRe zr$JUD>v^eVl4rX`U(#l8LaABDEKYcX%^8*5L<%JsVt$uL9j&B}Li0>_cb$@Hkz ztOV(ik&*CGA&;NN)pikdV$zkBl^ph5YW1hssnmgVc_Gi~h4ygbhF?nwYVPifozZiP zGu~Sbw@9tQ9OlOQk|iWg!0?+I+o~e4Gp%Jk*qppPqMjtNsnd#}YC9lI{aX?FJ3Bl5 z!pon?28x7cJj{j>;s?_%smM{?hil#?BU<`FS6WH?{=H&iswxU-kpG3+|*HymCgBkBBI`0jWrjwm8G0-D+ zAwptuE`;8`xB(mARe&!!eZvx8vg?E@|MlzFa>UHU$SFj1(BRQ>JFW_Y!CGy5AMZXtFwoM{TFFlK zwW&BE=BmGLe+~r0IRERzhY$US+Y9L3mJm|+30cUiYU61^LjMcdqD#xkevOGhZR}Of z&kmZZZ6l(YZ1S<}Y}5;B-SKBlK$oXL6HJ*jiSg>V4sHDG|EyY=DvB(-!B^EuipCFZPYKHY60m9VRzgF6Zc z+K!i)V86d46{HT^@AEqwLKtne&7NVa<+{!nx%`|WXkEH_u=z%(tn;4XIrVT?z7XtCR{aTEbw^Nek*sDp>W=IcH zLguPdvLkZ`k2EDym{^ZW&~8(;i_R`CXp!M?71*NQGg*O>qgGz$i7o{qFWgTU8nEnR<*w}(1JAc2CrarQ)}bjOIj7~`#INwHUn#HR34HZ zNavO4Tfh5ctMov}`{g`#(EGIi32CilcI^S|by!%-v8egicR2}$mWjyeg|J&3HhcAy z8+Y8Noz%h_bBI%t$SAUY4jm0cRnw4U6i7m-5Y+ma@jj+*3pG!3bZ6~F{}Fc9&+ix~ zuo;Vs?=mhIKvD?y*T?z_mE0#{IrMKCmUBu&OP`+~=$D$x9y-D6eSvAN+kQqpR^~V- zEu2COQ@Exx#Q40S9Im&nn3JT0U)%BWVOdaCK*Y`zw`nJr(Orl%$F3-(0T(|Rh zc=VLbdIPgEGFoX{s4PF3(7*T@(AB!)GrT~}l^teH2&1d4op~0%xXnB?^t04zSR4r8 zu~#yU_~^~b$~>O-o%G#@TGZLD^Pvyh_nZke3`fMmGhesoFXcEn``SDyp8ReSLOc2c z%P*%mEO%7k!QZHa$lctknzHBMEFP@*lvpgczTYlgaUec#qN5G7Pa zqmV9?=UuT}?v(zhchk1T4R*j%V~U!1tx9)Lh&DY{DawP+;gxK^^RrU>2POt?6Y!$( zS7kn4q{F^plw5A-2d^J6d+e8abY5K(+eUp3EmQG&E=wIB85y}?m`z^{c}*ng5m8wM zk<(T`&`vB)!zX9;n#-QQ*V)|JMrm|s`x8;4N7Sb6YgB^Q;tBLl`@(4>2DL6YWwB<# zyt~w6>r6x9v)x8nWSjn01?IIM7tdoBe0J2b78 zD3AS)hLQ~qigY|E8&J!ecjvEVCC*;&iq7HIJ<0*AKVa3DIHj@?BYOx%wtj-bnG?8VzhCBgT%*`4xlCUZ2x0h@q zhLlAM+a#5a-NBfuEcs|nNAjRpV14VHRMs%r-N{xb)c4T9w5C;DfL3F&72Jy(&*P+V2++bPMmnpSt0kUIfU71yQti;KZ*M7=q^NFIp5ER!Geb==J@C6iW$p6{y z-tDLJ*xL}p4KxEEzh^d9y~TWMJO1TUvuxLEc{F)t6(lUL)azI-t5>5s=>4XVG~|Bo zGTu07;5FfMAv=wdD7%vwMD z`t9}ybUgxsog3qgOj1pbsux(s0T~$z7eSZ>hLH;{oZW}-xu8rqw~;QwS`{u zSeTo2*;D?99|MPgv%{74I*qrU^?quRSLe+nueDsYu~aAh`5wkawTs$&wQ5UMxuqI2 zJEhYeRthj`XpgEyXWnY?>(?cZ_lGq&Z;v@0hb2NFU`yVI!0xbm>~P*b9S-VTk69)y z;OBfs-tBrH+M!V?L{X~Ryp-I+LJc#|Qo z=dvq-yhx}tpKQcR4Nq+KZjg`0f2jP)3ivWQkpal*TBi>GAvdLOsavl z6f}DX?HAGR#q5$*w^Q77<~0l-K|KAeY}-D*w>HA|*l+0cV0NqQkrMLZ&_`y6qI>Wc zf$jN+qnbiI;icBBj{BAzU0q%JJ}Waj2wlM&&p&d1-V6>*&qDpkRs24jtASppe)~*a zip=u>?Y{Z#VX~d8Li|(@+nxD_(RBZgF6Y&gFs*u7p-VBGu(iqB+Pnx7>*s{cp~`)&XGl=zT}mC7+g;<&iu+?*V%?R}C&;QQN845>(n z&F6M;HBE~FKOiK86Sb1|QxF7BUz)NrGlPoicRF*X*txh!6FryW=EiJb4PN_1Aolo+ z{&!@su9u&e*IB4Pb?~-f?!1z=A*-ycY+y~h$bicC6d2X9`K5;{Nq2ltPk!y{zW~9V z!TK%!^x$SO$R?>2O&N_K7Nm%-yM$;ATFH4j(Csv!((B$1|ASxy!y)1=z;R5 z@fh47y%BTZwU8o8+iuMvZ;_M#b+o@uSsIRuC!y;X z-cf%A2Jq_ero8*&{!q-$RHwa_E;QcuwPKd-r-+O5`8VU+PgmYH=6w)=>WK;ojlr{q z5sHi5`wnPbE_QZT-3RcPBpUl1~Bs7nfe*#3rUp(_2LB3!8dB*KznD`u68koQ#8{q>va+{ zUxwhYjwJk;zx}B!eLO3&sBZ-5@bppn5#gNf!WsI+ za*?s9#lhNQX|lpZ1>GQlVKFr-t{&a&f;isyB4-4okNYAUEmT*l;aU62TII>;DLr4P$DlEs51ZXg=Az^RPY}}UIqaPL^5S~bX3pDskBhHyge^355IB? zB4P-lC*=BQkJkYPK}AJH+Z1pc*Pnm0s;YoH#vnf@=W@3=!7DDb@jy?^;GjBfaDZ&5(y-0Af97gBnQvqSPRjf(uqu}?P~`q>*JRNkgwSn?g2!Rd{@`Xn zdc6>MZW@a5gPfeqUio`TudAAaIpG!%H@E-WtAa@-^6mkFYY+~b7&Y?q@!lWNZb>)? zMsO4DKD+rE{gjH{@1zv=V!O?4EjN1-e3enCM6=p4-<1jIaLotF6{B_d1GYNI1jI=g#cIxwUOb`i3%gFsqFHYx@6rhMVieu;@K0^g;5U0 z;U|*#xkW`)A3og7&&|EE7417_=EbRnJBm0A>&J^4uzFpcoz7t-06}ycH@O6?a2YU% z?I_lYF9Kr*3x>EzdBXR2(=skDl=9lj<3}4MEgPW^lS+UO)8{thxvX`^WbA2u%moN! zv#P$LQ~cP7LBO-unXRw@;g;CO1l^o@&`*uS9MD}6k9WTQG@Y23z(`YXmzuT_QFgxt zfRys%?bbuPjF&6*s>T9N^D_8gYxaBcE8DcE31C+XOk<;uOq`yU)6lC&S-$jMsvYah z+ZpF0QMd}QzM{HMWnwh;H72dfEZ+Q{({p$5B09)(nQa5q zOwVqtsQoq&G``~J{;656Qi+D023p@Szs>rGjb-xUaWn!sot+Bh+v^k)o^@BNLsV!T zNRv^~iDhW7IYuKjxpr17US3{102Wt;3(bs=0)(*in0hy(1HdL7av{hGcx>&mW^*&2 z!_4DqPGHPm4B_Ct^}^~TZ(MXwJHt=!v-h4@*>!vR9m*$HG-VcS6WhAYs)W;%-9C)8 z-*k2HO@NF~(;p z&T?9C7+!Xl(7+3$-E;&!WNA?G+6W}(ELl$S?k#jl_Q94)oUW>1dquZumh0Vc<3l&? z?}AInSSli~P!R2^QC;xX=X-O28vtcx#W*)Z-lpc{d<_qu@>(xMuXEM@`0%PIjgIiE zmr}C0&$MS2z|4ju4ExfJ?Qu9Hl_24JDM@WfupH9vAN=&xZ_=^tpMDt`HIQoZMLXNA zn>5(1?|V-yJqcsyVq=_F@cvOqJ9cl4RR^kJIT1JQLh^y&+M{Nuss>iR$s53ABnSCr zB#OWiU{2SwF%t#O_#Pw1dNz)? znj*g5L<1{isopK(>!h5WXS zL1fu8U|Ewk%epu@r!v1U0~FCcpb)qJBZ|MJ%r%sCNj#H0Dkv5PqOZ#_^O-eBLbuSNMLKOErLKYkX0c-Cv5x`B&G2Kq_{ zFAZNmQ3ZTqez@QHZuK!Z*h)a?KHG3Wv@(Lh`Qx^c!DzC|bFgo+M2$|C9nkZ^GAbn-PQ8K^kF~O(x%pDGT8%B z4lKJV&tsQs7T%UKJ&UPHQS4%=>sPB+bVYeYMB~3Rdn??228?J3u)TczVdVpqby zeAlk!+P^dussj*I2Vs;$!+s3gb`gU*`5gn~2$WMJ(=LXWJ8hqBJ$sY&wg*7y8Apqo z70vvw0WHKCk4k?n-X16Qa|ef5Mj|EIn&p>tri-1vpe9K2L;V4ckIy$_W7U9Y4n!c&%n~mS|v~ zPjE!wYruO!Byi9Dnm2-L=2z;0>atlH<8R3_Kl)f+-z|P>#IU~^2<_?R8_8#ws?sZt z?Rt(VwSP~s)OwaJ_K>Y*r<|{#zy@Iw77d36Tq_{gL|{WIWPUQjo%6aah|+MymwNDpx(!((0bS(?0M<89rh+e2 zqNp^FS$8veG}8{`0yfm-X(?`Kh6X++tnZxkIaz-^apNG>%D3t8wTdL(rHs~<(LTN^+nbl^h)hL2s0yZL@3(CsMZ@c(eBGZVCV47?S zgNJNzJ#r80={N0sxyCO82`KPX1@rGV$jK5%1*v7p2ZW${L(7b5pT?1G%C|6N39}B~ zC{zs!WK)xR6P$+Z3B!8PmNR=(7CuW2yg(?srj;Xzm0`0RxPz?Jc$Z~*RpN9nlGsrI zV;!j^Xzf-zYPBgYkWmna@7!=MLQ1}ogBZC;SKR!nyTJJ@VFa z$B{`%dpK1z-b)3_Wuk&8rboNoeCYG&Aaz+G*dZ%Tx+!6&`0295xzoh^;JCB|hsKHU zux;L_X%J}%5@I3O3vi1}f#{huSaWs7TdR@tG=#>;uhxWxTo%KhYt}tyVMJOv30dtz zEMNC|Uv(1YvVGzWD}-4-LBI?@hnIa?6JP6*&uT`C#zqzHGB=aM6HMnSivwSr@{cw8 z&XjsAKl`oDDUsbULkkDr>hYO#liVJo(yC{w8iFg@w|f;4{3&%b5)`CjJ8K;o7sSjm zA)4s5K)GWy9W_@b#+oQcGinZm&8`bZGK7S#PDj96>gnpnvEBK;b&J_RNLsu z&2z)1L9%s`$TCN)sPSm+)rdSXtL92-dM;5)9y4n~ZqZ+PmvgSB_Hq#Vo)x?wyH;xP zp?)m?ay!T=4E!isxA0Gn4cgVCi46nfau#IC#u3Pc72EJn(|% z+T}W@CHInA>xA>DscBF^{I6-j^eyR6E4qt(OyU!a*0vnZZwut8k2I!N-7&4cQyVAd zjj6&O6LgG=CtTgANKrgN6KM}naF;{{sb9|s#HX&#knBolcm{77AyI`vMwt5h*5^kY z)wKz!Nbo(Sb<8wrf^o% zW{s2)O+RBbLDjUl*zS5Aei1_n5j&yQ+F7j0Axo2_jgEiSDT|}+Z zVBw_$vL&7Qww&3)lzX(YaW@9sPT?JOuLEx@6&#N>osIM)W(;vPtrjT3t?I ztFA1Dhyf|i@^z>QYJoPY`I64#nP;ApqOh4E%@W_C*v1mIuf9#r#h_kiEzS;g)XYT^ zevBVd6e1wETdju%eGYDW>gsEy>rhlEBM+cENl|uPn2Cc=Q)_sldfcAXFxvh=?G({B zaEkC~I_lI5ZI*}~>D)1PCwvY^T^@UMo8rn(6rPu>S{|n^b<EtrN^ zJ!WMT_u8`(`erXbB))o2j}jA~QUjt`+J;cAIzT|M1DF7(>&2tdup|*8cql@Sx9z<&mkV_Eiv7(E=`THHz&B{K`hCGkF zV(U#;ES-t6$g8uxo)HwNZo$}G?@8}^a7ohi(9_zy`eKSeXTHs*goMrXu<5Z4MlnpQ z4>*uLRrCWL*xTq?=w6uib}07f1=hK4W3>|f;Ue~HWruR{mh`Az%Tc?K)iBJ<0dp6b zXybrom|!ML`~H-&e&$CLof36^Z6Uu?KB`_ey`bgCw@6rwBH1XbOF1r|c_tUcQ9roS z$!o={OW~<{gJ|FpLwK@UNd~8JX@6TKXFG@*3`I}(hFZcrA}@K9b-m^+gYJk;{`Fv_ zP|-(ZS~<%aQG~V6;|?Z*t{g^Bs2}G&OQahTXJvp|D49-V#~7vyYT=P4KI<&XjcS@v z0qCb}i+x5~$4I5x=URNmM4_3Fc~%u9g(}Cz-FfyzR&l{mBZEI$sNLNLpFXS{uP;3H zebuoq8~}azs)OcIiqGYOem#IDB=9<95@(4YI^{w&$71QQQ|6AOv2l~uEsa{^y*3Vs zbJ-Ok+<1Q;D=KGeMfaQh?NBMHFL==4W-hY*r1{%UXBe#N3xQ zmc>N^@jTfz!-N_RKGm%;F`9Bdiz!L9oFCQQo20~c@kf%*SiKfUDz$WBw!*sM4^J+L z#G}xOfCyn<@V4weD+%c<)Q!AhRQtvDo0?2SxQLos6F zF38L36_113IVkJfC3ov#P-g-e!d}&SO|hH8Uj7Kt`f9l4MZeJx1#es_@5g8rB|5rB zOhhx;;?b56^N# z=Opj5=5VPbvA^CaFRUn=rtN+0uX6J$1QrB z!^cB8MgB?Au7nr-1)nxF(Lww(glu(liL+)-d`l@}HD$$S?#g-=CCt*LA5Hnw<-!U+ z8#%*3su)s;GgS&iofY+bKA~(f+Y28))&k1%asFHefmG+wv zMM+v#FE}tbxY*Jeh43IxBXD9@t-rZSEwn$U(i5$)wHla9E3{$uYvM%MHUd-oeluH| z!BcE!vRo0OAbi4v54nFKTEWk3@}9Tt(FcoX8EhXXeO&!$`o${2>SS~SjWQ@K%5lnv z!kYsC7TC`U{`hm2&-*bAMkm?;buctLnvx@6eox26`g4vHpo4`t;jL(v#wg5MoElbj zd_1e_If_!)y80aG(I)Q8SE7EQd+ym>l^;4QVZGGU6HN*EEApXxz4G7of28#3fgfYi zteTqfEagNS9J)_JE5gs+u{|J9udAgMg8*sofC)?2`n|VC@G)BMr^{3`paS+}V2*u5YvXpGMD-{f%QtH55dIYM* z4812nEbU-BT2NZ+vOL*zw$DYwXYqWP6b|TRi7$4(;b6krIuZcDR%S+eVXk(DmDdN-LmKu5pK<2gJUs@`b@_Wnc~pwsK$#_q$zjGh=C%QsO0=LjA3%P&=*HpU^<#qgO)D!G* z)cPHgY!b8W6Jx))n?_yikjB| z>xlSp3fQ>xo2$tu1{;C7DwZ7Xuyg881(ou`u`Rf_`aQ9C{1MTXa8Ky>{H{!Vl7Rk4 z&Jp+zPQL;w>5RFtfai`m?xxo%)G7S}h@Ovd zc@(tZh>~i{!wQu`5_a*^J&o#X$Xwe_K;y(#NJ34%pLjrbU^sS@&}P8*s}|nnd)oI-V@KT9^(7{6xSx;tp3l2u zx4Y-%q(*{Ii8t$Ea^tP#Ri*2sGvjX4DaCNh2R+HLQFB{FIxw2Dvb-3!)KpuR6jc+z zSW2@CJu6I}Y48RqAkCjrQqk>hEVfyjtmzqff2}bMQz!eK)p~ojT#c7+YwZRh^uCc8&WYiKSEJ-IjBry zS361`hH^Kl+T_V=N^r3`72k&y+%P~pmp?D3_pDWbUS$7zWgITeY_3(x>fs=5Pg*Y6 zEtW~_d202dcI8}HSUA!71joXjU_n3Ds+M-39}ZHbz|$AKGy+Jy-?4SMCKLb-o3zfZ zT^Eh8Pg*}5$y6Kj!-d+@Z(#w6Uc4-hn%Ym%tr2oU*by_uj-irl@()JMjXVie=HG4} zWGH)Ney#@0wO6D`OId@Ll;IKwQw|abpW4-)YV4A&I7L(J8(xBZ%PYuh7e6$KBhG3} zs^(oYSu)S6v+?zQi*ZOG<#2|bzk9srYx5~u+XiN~dhwu}codH{=jhpj)9970GGPl( zF*0AFQ-BK)xQAz`uNbN>s@>ow6ayARySw8@IEn3=MfE2@PzWD&A2$inJwH2z&6Rl! zcog8M%;Aj3T}Ef7l!hKAkqJB4Sy8%-n2ak0M8w2s(QL9=dd*PvK~!X_C^ar#c~mlJ zK$;d+rf4qs13Iyd8@kLw2DwawIO4s};R0{Im1<8PpIahbRmTcm)lhJhf{Gh!ts%XR z8r5|rjnXSrVqUJpl(}WjwmG)RdmB_@&6*ahll??U!eN0=_IdfS>KZMglTQeFhoNe+ zyqt}P^A_4rP0b%kV>;{^Nnlv~Rf@ya0@PbCi-#*ao(YTqG5j& z1KSNu4fZ|zWqHsK@*9-Co%(W-SA^K}&)TR1&P4;@Pg)yoVW z%cT26^)k~1wnp*?L_a9hID%$AD#l@7crWwpGD!pNK#FP3U22gxI&ZX|l+|Od&n@?M za9&SNa*3ihzGT1VSf0$0vo~-+k26An-}<0(|EK{Lx_70;OQVOYZNJ+{|AG{LTI{hxxc7nj+>?4BT1{R;HS??E~*36jR_qsGC`YRW*m*tKvEe=JIV2$M4|ryusgnjU&X z6@7{;GST)$tBge&`TokJNh-sj@YvG=-JZq2iOx|BTT(qd-STq?q4TnmSpa?mCl<%r z7QOAehM#j*k0Z&`0;Nv;Oeq*?# z(U$lO;`<(k$71-tK=UstvNaplcpSy$F6r8>9DEX1%cAw17-sG|z$}7EaYsvOL(;c5 zh}ff~P{1GB1b(x>i$dGLX5c^TP7o&)Xs z%qjgY=sNd^N4E*3^UAV|UaU0=+KqElTCsQx?7wSrjOBf-C*M+=>o)BqPx%JsIv~qA z-oNFQxKyc*HHf^_Au_5kazf5+eT^NtHDw;D^yN6;gogS_=2iV-!-nZ_XLuMn1FItn z@Fl4!svZHr7ELZ1&;JWJ{O%JRg#Mjgv>-zqUZyqkkw>C3wB1agp2ZI^_3HeU+xn7T z>^v1+`e-(7dIY1s=&3+!eK+yfE4E@%OojE(Jc<@lsv&t(6PfzXN8GzsDgNXHho(I{ZM89oD`VA93t+>}fCSDThLKW9XN zGvPb=_L$I5s=WIecn@_%p{YlXl1p_6(&M*sZwbrm87%{#xTC1#bw0^H1+CmXF~ zDsM!5NYyXlXrVdDF;*wTFT5>}0P`bkqcYLLf-Lx;im~eUJnHm4R8oE=CSg_ar+0Qj znbrlZmFIwDSABein((jRq|)o*$wxnyO|{?6uX~1fs4QA!P0i4dmOSYosk$rs>q^D4 z!o^i<>rqCoF0H7-!pk|$Hi?Gc9L*uu{XjVW4gw$mv`izlgn*x{XS&#${ztPAJMpjB zDOU!n3$0K!UpF6P*>6b5D=IKj*VIf~eqWn_hP@tXeKPpT@Ygf+^kHtI;D%Y;n`G&S zE2}h?zZ=uy8V68?jn-^yY@L=21JritcI>9yUTI4QNon{ccgquUk!TdY7mjwf0pywz z$?SjRk=3r!Q|&-mU{f(tUvPEGG7`^iW$7hCRi!gC-8P5K_)*ExtPTKsx^OGei$jPi zb8S0RI1HczR<@_8#TrxcV|(6yB9}uwxv`yCgKRbiNbB3C?~v;c!?d->%Z)fnVw+%A znbDQuRlHs)?J9+4qIUfqyu;_D3J}>d9xd`>i2p|bu}EGz_MNjyc;wmW^K)VF=L<=? zHWg9({xOP-+hRH1mmH|9n&>=nKvLE33W}6%jq_hlaWoS9=r5qeP3K+W9O{Uc*S>^g zC20yz|H3@uJu|2lMRGVq@+L@D`@1pc7L;m(k@UT&Y`#@oWLJ?DAm)rd%=HnZGyXO4 zV)5Sgudc#k)?&9vE!%lxh~(}e8z}^la&_d%c_O0kMQw9u9(6?Wq`9NJkgv_uTmzLv zkmTy>I#Es>guf}i>4GSP0SEtCy*xju8Kyfmw{QKxDp|T441^Var@xTImJnqqr5p^; z%GT{UU((iam&CQkEpf%?jrp>l>Zt2t#zRtWUzc7MI#hoWEFe+1V{L6+^ux>P+vGMh z?`nd6V3ntV2y$KO<&oIg z=QBxO;}6CVD>2+t1m~@}XP$ok=GS$d8_pi!`yYaGwQ_){x4ARE!a;ggFJHZZ8fwvu zLhf_=l@g)wA51=i4$^262$@vn3(<4s}PX#9|vUBrAax&4{aC>0m=da zy>!0Rx|NCw44O&l<657N+H2}>?SAFWeq?yB+5qIf39@RqLe7qzov3n)9_lf1T-{*0 znf5JZyafP3L^^MNxQg2b(Tc1#)2+F5_-ywRUKRN_N^@Keo}AP>AY+3+)k{xH`vQlG zNZVGUC&z>`%gT6{&<@Wo=H}UKx)K&fg9!ae@KtL!o=p|o8KZ%5jxytI9B^6qqpAP7 z3Ljf{rnbO`XENH7ueC3{%Cn~;K$On2hKdrK1R>=#PrlfQ7OB{tYNXU6c*G0`OL&}$ zG~%k?v6#|tF!5xBY#Ix#Pu0}L8B z^KxAMZMb1}QLtVI$fwSg1>1;L;)*p9fSZL2i*V}5bbdz*Ri2y_GukNt3GFem_d$AA zjvQ%JC@7#l&H!HE^hY)%BaN_g1`cJY1ZU>}A+Il+c2NaDy=Kn#vL$M}I{+Qe3RJne zEJWkH_G1t}-ZMN0NL%ePz4r#slw%9Y<|mh4Cs{^N`YJXB41LRIeKocBG+X?xE2`D? z_%}|y0q1mXzhSpmw5siv{y3m2TGAjY-s8YKOxii1@0TM04Oq+2kggO#Xip@gnPQ4Y zjv+tCmvP~cJmMdwz~z*te3e5ygCoU_KOju)d%XL`HXS&2Or}uW8f*ewqDX&%vu?vQ zgxmp~)Q=&Mr60=!v3!|))CEr++CqK#Y*^5kS#D1s; zrEW1C`lHOC;4(87<Sba@fX3jFGlt-+D0wyilY46BC~A8^`_bePC>ZZF6>vvh#KM*qti8Q zPQHD>7h^IFv9FoCcnoM>o4C9j*j%HSrgPueHr4+0upK~nMUCr%$w2}FGzaTj`%X$WaZOCOn zma_Iu0`^-vM3L%_!}AZ=U!d;Qt+)xeP}60lULaR!2Z(+}9OUe$fa>HTHDYvrq)pjf zf3%iA2W7F%z8W9cjjOY6P?pNXry?{V>NNz2@6YvEfc^mnTTxe|_4F{1`HCSQa!r%+ zm38eb0CO?>=Ldf0mVO}3h~6S{9Tq?O)@g%66^v$|Y91>m2=*Y3K>;_C;w}X{;8R6n zvhOC6b!u5l(E04C9n7$N)H*f-QFf^ytzH1Y`1P!zo<6gl?b>&#K%PrSOU~y(48}iH zg%VAq{1_PdIR@F!k(O<*DRA}7b?ol~Z#VR&r_2G^)RpR*A7vAOs;1yCIB`6=VAW9> zT`EG6QdP$tTzOTB6{vo?!j)#%gQJ~wzB#U|%4JJJwSo=ci_6V{WGVW#(1dxUj^KJX zu5h~vaP)6hb-m3veh+i0@^f;+Zn`4;r2bNnO-iUsc&>4d-6Xb_Ps8y= z&S*cnlSAB+&M>-k=~-Q?_v3(%-eMl_T0W)TX>GgsQ#D+E)*E*(=K}A2P#BG~A7z-A zst&v^J%nnX8^+zi%27^$D|%KwW;Zz8sLgv13QQr7{GE+EoTymN*IX6Aq{qi2P52Xg{jxV0$5Brwl>X@X5#KS6zIUie)G* zLnQh`$-FkP-{~^&6$=}WE!R)p4bk(qIx5w4Z9Rmo1mH!85=7|EX|klYdT?&H4bXp^ zw6FHIw|;Q{%)32a0q+qp6zZGq=VLczNY4&vW`8fpxOSoR#xGa5ZFqkH+Z*anVlG?t zvd(~(O}%tsAe8m}PL4GQd|vL_n%YOUIf0AhIBmPsdh~!$t>8X_o@AC%fSv4)-RZ?=9ODBOPf5_??#->5t)4=a#q;xQskF!&qx~jIo)ckw0*TPX@3VQsE&ybbDv3g z5twkf#D&Nhn6~@0L;jUs<%Uecox~8cyeQ3=P9;%(2c@E_Z|tM80nIl%rS{Q}?i4R0 zpiG;a`C*>P3T?2B>N;^ef&Qnj(ZS!{7jf|gaq;n!3I-#_D4(Rv_LWXQjP9*X4Y0I? zTS%F5k)b9Xdbl&{k7!V8H6ozjfO+DlKY3WslYU7%)oqx!Wan(i-h;C7Q`~wo^Gy*s zM|ZBJ5H56Y$G0@G^^ZL7zg;XNv@Ju>Tn~tFVg@m9$}iT+jUJfUb2r2-EiN&~Xp6(m zP5`o1J4W7=TFk>({Ab3czS55cc9Fa zGTUfDp;=OSqSeIsOB@WiQj%>zeNvlbd%b%a-_^-_87W{E>>2 z^JA5FWIg?vi}+3c$lidUH<^CGJofN7a5t}i9clF(@kM<}-X-gzNVjnf+|irY^MF27 zBrTx5=_>h-9!10Yw!X)Oy_HJBk2tu@Rzh=x#+ER93lJhIzK+%{M4WpHp8cEK|1x$(AN-J&ebj{ znpB4;K27>jilxqR8k@^^N7lXwQ)q#ierbDoagSI1&VN&7Zp++$c8GmsXLPDhspPqs zFSQKD`PTcfFExLJ>?6z$%?fg#;5DZObpg6(*4_9-@ny|Poj};jmCGY#ekzPZ)J8y%qJyIQL{2aW^P}T_R#4{TA--M*B^S@ zxHI&QM^Emo2A<@?%lG?@6;$TG4VIozHSVsMZZZ z)8KxunwP%b4rv)_-EI7~_gZP%0grMne+=MCYQ@{N!*XeuA(;PQh_7WS-=cOW-KK2*%nyGh zXa9UHJu+__O~>dT_maGb@(of0ua|T5=)2Y#IL>6|&LvyhB2oW7`~JIqC6~|CX=El^ zZr2b(a9p|X3eoGIUN^uG49OEl=r{(1>)x1&?k;c41&9gck@D)6B%yuXRkm|81cemyuEDuFH75KAxeaEFXZJ?^3f)O|{8y)weRaG~rT!CN)1_cfE}! z*%JJ4&Go1VF(pK7aqjD;O-xXM0dFRpPSrTIt%_JH0=$y*f5N{Gq@`!%D&dv-4Vv)+ z>Eg*fs&L))&syiprwP_({6*gJ5ifvQoq%z*%4#;SQ$4H*ZK;JNb;s}aOc=S09JaC& z|8C6q9+49Z4k*7GRC^IbqE%Tt1a59AF4c7t#Sp8O`Ug`klfOJ9gsQRENMw`}^GM{F zO|=mY0aSTU)~q65<297FT{nUT3zjU7?#z;(n4;j}KNpeBfuZ@IOD&|i8b;HHl3!Le1`$QxDcY|#!!Vlqn(8VP;y#N^N;0a)wU9QyZ+f8|F ze1;w(TVx(jZ-5geP}CM$U~uRKHcD>hgFUzM$_!oKfIyar)xUPGDYm=Y3 z%d)WNT%Cr2H+M>cYmw2^P~bbYYs#P@g()q&R`uLpY!9Jr9Jh7x&L_H7=|P)H%6Iyu21A9*DM#NtNB~a)IM`C?L-p{-x<)m+_>%S8k z^M*J0aYB^Pf%b^VaZFoYv2*8Fkb`B9jg3ie9fuMukIAflbn`;uvEBr~C5JdbrjTfk z|HyYzH)xkM+bka)2rGVHKd9ZlDu<@Zm_2_BSA{tsEuE{CJRyR_t7)jCHE>qX|fz zJO^n8r*9dHNgAp5nxi$=DrTC9ykqG z#?9?XOY)Kd`a|*q1C3)GB)S&8!VyHdQ@?wxLialZ?iMmyiWoCy-+uy5^a*(`DZII} z<@$J2BXUWD9*o~}>SoexXuTn9TM+>1T|N2t*}t(?ExMq=^`ZJ}XO!JXnqoeG-Fj@H z0-V0Az((TCKWVc}Xc&0bA9(jE$F^<)#7vu&x4O=^Q16(zAZQ5hM2;xg6K2fwU$djgWeJkH}(ST%` z5nMqp8MV^EA{-+XvV5U1ash=|0cZ>jK}2TSQ3gVtypJ~+6=8k)ma|+%jov8gRa%`vG;0%W(&EUZc zD4VA)WlEpc0;QqQ+}>In98W^9+xQ|Np0{B|37F}EkSWb?X;cd7Bu%Sp6@eWB^>gop zy#HdT=(V@o5V%?NQa{1W&CM;j;Qmmt=EaL!%lS!TC}c#|#(u{wPd9U>%O>~u0?HM1 z!WCd>4Ue;Jcx|&a@brlePSQO1aT@W%l$4ZMC2afeUn-|M03206x>f&GGK=-$&SsIw zvEfY-ZWM#ojpXKk!|-1gt|+ty)we%B<43^=FJ*}yKWOv9#GG<^Oo>C_)H+z*VV_2c(t{8pk*I&rjqvz#fLI&~W=8 zKGg^tir)NQNzu{im%&(u+2SigL=q}aOy|`6BOU|G5154P{T5|C$d<$;-;ri=`qcgl zplpM(Z+ZfJAny6RzJ%5MWyMd_%;4OXwwVVD-%BP}#JRF-0La-61hQe%U&lP=OtLHk z5DEtwd_R4IqdBo0`2?~ezW`7-p0MkvZyutr$qm5p$Eq`_ zBuzDpu&QZ~Mqt)olA;Z}>(i84*genO8{vKb2a&*mwq>+oCI?HKSAMG?aji9MsLkN4 zFh;2lw+J!G;Z#|m++&lCC;oJI$p=Fx1nrUL98tF~H7mr0}M@|Mfbtii4`(dMTid|L}w*5Nf6)lue7= zQscaMF?L}2E=xP#So}8|_bOLx6H%V{CP%R)0}|ND=ha4#vT1mT3K`ZU9sk$)ov4@d z^*BdXI^c{24RxHPGm!Cj5N;`qkR)zN`VQH4Mk>*!?%{381!efAbC|C2-3=5#eRH9K zncT3&7dp-mRYu&lJ>t~cCJ?Oom07ahElJ{Ja10F(o4C5V7I*1CS{4!#DsU3Tn#@;2 z)q80hn_|{-vy+)+pZKs8IzO>+6;aio*U{Z}zZ~F3mxfrSuuP7oRB;995|Z?f-rCz? zKr{HLiRqf(<%U+DEJT;fBB-uNd9{Ip3>={WKy-{#1E_NzY&uNLdb>hSU=h+#1q@|z zBok-D)oh?&Gg5JY;+LA5>Iv{>7ceIghc&e8!*t1y`{+N{=Lpj%*IOQw2qW9h^~U;s zvz4m%!`*w^hrgOSt;mU`)88mkz;2)bQ0wkc`LfHXE06lsI^FXz*7vi)h3^eq>2SAQ zoQpFAg!CmbOyh4GCN*j?w+LY<|JX36fY^fs+LJI#=b^XI4Xnpq`SFFlo+j%wUQpW? zqTUQh;omId2!AqVGBdB9CX)!N{+5Gl^%#jh_C2f(`uC1NZuNw{^{qaTJI6B%xBYec z(_ascSXEV3r7s=BGz6Y$OAm-Y^kGsZIfxc2osd6M^y*?9u=I?RoM36l8+YjX3&gEK z+xy$r957=*Es|n!p+2A?o~y%1C6CE%Ae56X@l~LZ+_zj&jmYPS7*jo%6|3f7#uu^l zYr2o_DGGSsjo;AgFwWnNxgT6xpdVBmkE-LUBpviI5{Zg4@)yaDgxn?N75^5d@7&m6 ztRDA1W9zZj^r=Avlr>mkP5mp^4Lus`a>Cq$46x*TX`4rO#qnV2X>T|g2+aOD+R}Sb zR#Q{6=O*RJlP4^`KQZ@I+TdJj`{)LqgN+syfskbAn!)6BLmKcYDk|$~G@4R(*-i<+ zb|d@F5}zOI)Ys(XWPN`%%>96Wb;xTx3266FKk1s@-AIv&caMyFKGuDu+@q67UT*23 z`0=Kp=0)@~UuX%f`%l*GFHBi5!4bAC$n)R$cA$*+@_i_AAk|D^*LxqVfX9!3&4bkeTadgn89&v2s*+Bu}9_Q<;99Qu#AE7e^Q8qrjqw(1ILR$ za*$+TbZPG3j{&@2|Dt<=;A}<`wn^W-3|4$!dC6d?!QK>@=rbhpD`p7F3iU7YBIgZA zg!zivzHlK-SI1O=l7NJ#3I3}-!WzhMfx!Vh45dWnww)T zZzow?&Sd3>jIqBz(<$M&oKQ|ol9f+YAW&P|o&}I?ZBt7j7Q)MvL)@w69*uq6#z5)YG`|F zt059YLSVdCueMM|kb0A&p!$<2sG+%fjZy{mQA!rSsCDX?Kd^+-di~daE4zKuKv^qt zv&jKyz3Toy1Hn67?XcqlhYtOX#|0Q(e%dBaB=%^9YROm|abc0VXsSYpE6OI9FE(81 zDW_GjVUjqP|DM;X-M|%bWW@}3ItAkw@Ui(&|&&6fc@>jXv5pfWWq~yT5MNPiL z-~!g&may{y!t@M)8Kb%#0^1yjWM1>>-~~uB<4q6&1M4kZd(~XINGi_8Dtps*$-JVo!{Bu`WQ$8Wx)MO%+GpmCHE>bjez2FYDLs9<+oZtx2{OQ#bo|`Y z)8n}x=mM26YU7Kc252|(A!iKQ+UNP5!UV30?lc5J*G1eFe%P2_6nud4)Bd4kWs zIUzsB|Bu9BwSpiV;ojc>3e=0UalF<{>L(l96q8FAe%~^d+_WKr4y99tY_hnes~vmL z?jvneFs=MB0czb88Hm%`)pVw655GGnm$6iif+cxeXaiN^I3%}Q6ev)ZR8KNgnGsJx4?`x->pBZMv-Wjf(tMYSaO^3zQ<**!ukb1GXK~xcdEU zK&Pp#2#nhK*u*{z-=qrM4agEA=&Jj(e37UlQcEuT*00`<)Ax5f@CnIYby)@5yUxe^ z^mOxo<0N?e$t3?mNo+qXSQ?}yBu9GUNBIkyT3VB+#sYHO5i35(Yp42Ddfgwbw;aI@ zw-c(RZ7)kfgtJ%-Jp~SUI`31d;6(z4%uc@3A_E2tkhHkh5MZd$y2=CvBOG$umjQYeViYo)HumgV;S*EV6tOQfaMv_A( zfNh@N@I>+5v6b&{T5oK`+>6-S?`-h07OS1ynfPQ5fT@HQB zk6`){>)dDN*;#c|1&Sm6n_DsWY^6aZiZH$Jp4g(1S;-kvmQBh_qhCEsAGzydB>o`r zhB#RcdQ8^U)HoD4ly-n-ftk6aUg2ld1OD;eb?LMKb!4R%x6cwsz*i|{?No7n zn)>o(H#BIQox6{%t^Al;1iIdffg4y-sFJK{WzvnLxfO2N5ni74&&z;aR@FYq^RYG^ z;Og^AFzL(_IudIjFI#p96>CqUPn4z{!7r=8d&UZkup^HCY&yB`(+VaSpk_T-{tuTxaq!tdP?n*m ztJY+DO}kp$nxheO6KJ$`Z|evC5WZa#LXpG<$dWg9{d3ilS6Pj!F7PnB_W4!rCZTPQ zxS|!#{!MM&{tM`i-a#DFzaP_RofDn8W>K|u<~WjpF?ouZpKvkVBWXX>lag5+5xxVl z*-Ma=^u8?z>m;{Vwn`7$Mm|2@M8_mH!}6XKHo|A+{3dcysrQ^Lm7s-zRJ!=7)s5J_ z%VeY9y&jOax@>=#Z54*6?%2^$p*h&Dbl^KQQsV;8^AV&W1)l84>S0&waFS%;3i*(X ze1chLi&A&^)L~i|oXJ#Rwc1XQkL06)gN0ECY&AZ4u6 zI6zd>9r5SM*rzvt0Ef~J0Ofd1ueRT0+`dBoyz_-UvRW4v8Os*QntMBnhuxlYkj8@H zS8p!5sZ`g3#9Hlq(HNZ(dy^ji<0(lT`b4*-H? z{dr`YEC`vpKBS`eTNp{ltgE1Tk}e6ij_lAMNnO0(^#b+@Fd=lc=9ZcJRbTc!W*lw59(W4C^yg4qTv6_YAG}Aqw3vZ2cQImq5BowcHAr(W!j&lx8P8K(ayaiMQ zsWgtJ(Ul-#$_+6HNb|b3B>603bGC8?)N}$J6?C8xAz&i)#9=1rRP z75;Wx*beHSbO{)5AV4`!e9E2&R)})?dBY@+7?6{*}&t@<6kp3 zAc(@Z=(R0OD;@BrLM~-D94=kDBpZs6?ivop?T$GF)2Uy+$lMf|X;}JSa@NdsnITwU zktwWFJrA(eJ81Ma8B~lk+dT?4!F8@UlKqTuW|BCe>45Jdn}aSxz*0gWw*<%ITu=_) zwm1D1`0_HiUAaJUng>;S8?d)@Q+sg|7?g6V3FY40DBq5RudLf&rM<@2o1n06dKGuL?(?5UqS<%+*f6x!^t3jtJ>@B1a*n??_iPN#y7H0>y1<<4)4bIeNUyM`U ze=1g#WDv^a(SVTj*PLAarc720iX>=tqSFXK+H?cz`PXgOJ8+teMQYD2Dl-6BjrjrP z&OYSeVJb`L&|N&F;rDlR;v)_~y?4@-lX2|U(kF0N9jzvdx;GtXokHH6D&=i}M{T_U zu9v;u3@8PAb8>Tg4T;#uW=8F}{5x*DMT z>upXBQN{y#>(oZ$^DgnWc5gj?;vlbjJMM*2I?Lp2SKEFlXA#jLl1<{ulhI6c*!mNP(Q*%3)Kj>bOfgY61 zKb%)r7USDfVx`8BFYmS89SG_Lh-4mF=ZpD{{COD6eI~(XEf3~Q7cgKaRLU36znoO0 ze#p*Y)}AOk4+#B;Hww>{b%1OVxZf&e#j1MH!}q9tMu&#fvNECnKP5-eQ>zv8xDC_( zFvtpfj~F+EPwRG+BZaea-_7et<}0oYn9q&kqC$ib-GIe%*s2hNBz|!iQqm5WO{8}r zF%s*uTL9KMRHWf_oU)XPO5w@a6`%upRQ8}_6VOo?Zt*g#XUTYEn+1?QoRA=QeGVPF z60;b2z`9!9>;cEC{x^5M-e==3ms|DqUbkhq=WZmu15V;f;vr0f0(&~DC}0zCd0Eg4 z2WSlD=mXqZj5;7J>M6Yv);%DyhA%||wpxGm^Aoc_aBur~Ei6x^dk~Jnq=t7w;)ONa z)(Pi=Ool0b$eNgNMuus|Jte|_r*l{ylqvZbTe0qM-qE5j$w%t=Zt)Gk$C+&p*NdJV zO%4Gx`UI-cf*FJ5Gr#Fq?ht zFYME6)4!TL5vgE|ni+ZNMnGlneM$!OXdbwyr&f@WWJWIdi5@+6<#Wk*Wqu$p8wjof z7w!SO7ueGEM_D=JSmfPu6G7O&X8IcqYeE#*oE#di@Ds=qpu-mhrQD1Wk9@ErEf4M4 z2_ZfoS^@{h9vDpBf}Y#;#&RIV_mdG6g+h_=)%2J9s+01~3R zyT0xSAM(@He$GTC(BabF$LE)k)&1Xdr~(1Q&4KD0;k#5mKz(}>Awq|5C1P*EcU-^Y z&g)$xPtFv15X^mr+l8~BQA)KDHd!q3Y+Bgs(FVR<-~hlIZRG#+ze7&^&@2VT97{Dx^Vtm+drSuW8TIKxuK>FC^3emrc=vBi!yI@WH(aAqU#9Mz zuH-P2b&Iq_IiKsbtxWZ*J*^oRtTSF*Oq`OLOON~V{^4I{*Mbd)<2k7Fvl}S@5)IV# z;Ufm{SftrYuK!R+01piSgG%}V9cYrwXe&6d#rUnvoQXT^#oD!;SGi%=tD1jtNx3_s z+Z$6)0#qZvCQ-Q>w-bv9Q{9`T3oUD?zL6osX_r@Za>bv3e*4~32UD$b>%)9yI5=kT zM(cQ_MB%HA>&U-LL1#6vX`C0n_ai*AIOGN#&QXAjBj%wu)dlJ}C&Zm6`tV~jZ=o?W zJo;6>@B*2Hot2vkXaiR~Tfm7ufkxl;#6H1cbaYP7K7BL0Haih{NT5zD0NY?8E?*P3 zn`8iJO_L8C2c$$(`X6!`U&x6y-5XFno;8$?4(Nz0AA@;!pNB57e}4l~sTTw1(*+6b zq3UrTc*d!53QDjKU$-)E10*hy3%R6Bh8*N<*10aPVA-;V@waiLBc~plG>$Kvt--a{ zHPmLke!BGGz+=;m5%mK*4N(@$ngm;hPtLnOB&zitRfDN*;|f`4(kGR(`^dO}g&zhR zOewDcYKSy|L1GlkOgM~CDqts zNlGgLwCY4( zIbSUFaG<1^cxYjfzL~W0b0$1nhC$Dtr5o4PeF*N`fJ+L69AgP4lBss8?ii}k=D?;o zlz%%Wo3ox7#wVv$o(~&n*4!1@5U!8aNU_Xt{?K0Nu!=1YN#ccs`Mn_o+jwzcL9vJ1 z``+0?`Bx(YAb&Oa?CA~r{4uje0QR9Z^x^9?zXr(S&>!2NcXkD7Fc9zrZp6_{7ajrz zQ3Am&-WTB&#=oL7G;*5uYeU2P?xVIx6@~1;EmXS~lBSV69SLd2XQYzgckbNA6>gd(?2jlM>KV&_X zWlM9$ry3l;;H|g1_zfhG@Jftx?ZO#6uR!*JxE|{ajr#f*Du;k8VdH3mZF;N5degUm z6KQ2`A2abYqU1!9VFdu)`W*^BD#yPAfj?pWsq_arTZ*Yrzkc|(Ru2sYQcj4|Ob95H zI$Y72M>z$1s8oVAw|l&K+1WrqjtE9@+g%{M$N&c4Je(FG<+P0B34UKe_bb^*T^Rpb zC-^9v1JaM3hI130Er`^D%906iayi?s%q4)W;(o01+F}El2LB$QK&)Ba#ZA3{$~ok} zYz9R@6!iAx@;_%jVwjTaAp!8TfB+=LESh zRWsA_HnjagJvo&GRdWWL-35ocHGMli=7R;_1?ANuv>jFTbo?l9hYo%SsemwvKsGRu zGxqWE$;lU?Vjvy`11ERZe?CXjCO0_AiSM43{y6L<1>(d3RY0c&PpMbes(j`o4S(-z zN`j_~O>25U8}!qAdMIdL0b89{^Q{uEF@cNVPiA6d-q#j`CRazdgJ{36d9}Ef-R8)|N)qaan^;$w0m|>e@jHG%}*tu8P}0rKM`I$Ku5&i`4B{ z+)d?`oM2BGq?Z+7h;tk>}7kk;y5(^UUj8*0(VbCNg=5PC8{Vf zAOmoXgeeT31;|#`d$${RUnQx>$)t6mW-a0s-kST$tFG-rVT~?b9Tz+d>fKopGs4)r#(dn0Xn_ zF9b+IK)CK!YWyZ9m9-5(+U66hy75>97fjI_popMjBnnZd5CMI0$8NE2tN>Hu!EZQp zEgNcxPhefXUKW92xWr{$?4IO0lb$NqpZJc>RggB1>E7Z(KKlJ*OC#RXywn|o^5!~lPLlh*1RtZQ1Urhbaq;3nq|A6i zT-)Ad4r~U15lv*Fq;i*5n?bZrO%_5%z^8~my>)bSj@;jZ?e@0>PXEZx(NPB$aVlpm z+(%!ycW!xyd&38ewA)?Na&K(i7|8T4&kMBn#V*=n?s@Os#t(^>`T{3)BgGteGd1FX z2qE+I;)$HlIIdTRP<8-L^i&@v=KfD_mXQu>wQtjQ6D292X``!a1-X0z@xn()1cp8h zRE4^9=GfNsw8Mamh3r(({(RDalr;WRhMAaD9D`!|lUMHCz^2)l0khfE>Nje90U0Ol z{7biS-&?-E>~OK8sx`5~B!eGbob_OeBj&R?Y9m{7kHyR%{`nU$jSbf}F7n+_cFT$wIMT9l_MRrj>;W-d-CBuL%fXwRQL5n^dl;U=sfViYAqKGsh`zpWcgC z6^`gU0621LzepD!sY0E?#oW_}OZqW)?++xUAi3SsvK;&vlMbj#@mVS^iMa6Y+5#~W zH*#~E8BPnxyAsim&Y(Ujv@#_PLUBAhpne*FNN1{JZ0q6iXZ2AIjh)x_7xH6f4SNjn zsIJys9p!f%J5j(+W^4W0h^Y!f2uKlHhiOO#V+nD#Ps(bax^Lu)IESLOl84GHbh$>1 zW2Y(-54}J+G-~0dLlbJL_)px(&ecY6?YHfi;KJVG(3OS~i7%rwTNPme9E!#R|55@e zAkQ-Xm+ibegT`+0<T?z#2X8*$ZPR|KmYiF7Rbb%~jDV1%VO)31>owr+5_j_Q@`87XuHzPG z*Z=ZRgt2r@6k}7qNmL>vYxjS@@CR`%B?V&(LF-fEFZ7Dib6&K!>$l&25Pks0aO7*8 z#H$z-Pr1-*%}dTx9eBmVq`o~4^FjfbyAh8=w1{uskpbu_G7>rK^OPj8_|6cp@d#P; z1F5hn26dkG2_tDQzB{rcHM@aI0P=Ih<1`)~3D)N_8923JuG1FdQ9!?4(^ncCLM2%} z>wy7=4TG=`c>Iu$N(lN_JDxJiZ)jI2o8*{)3D}_LxDP$W%p)5J*S^`>8XVU~1lW5Lg#c`2&L2u~K? zk(RE6U?hRP9+m)3ejfXk0{^ttVH&%Ep~rRCAAycIF0*EC3WDeLE*z>t89Oj@5Vsdl z-QUhYK004!z*~U~w1_e?Ij+5IBM^_HP6*(JIZ5T^q}pe#ssg)j-Lz&(L-scUq}nBj zZq5+D3!s?*$8S!K2R=u&6>0+@AGTYt=ZAmJq!A6eP`L#{)#$;Yp{_sg?zTc~6xe5* z3;O(S#3|ZeXzH^jq{2v6A_p*CBWrwHer2Z*c9{TlP}E9;g2YVwC)6>WU}feChg?2G z$J3(TJNTkHKHvU@&S2#r4H@@;XV(IPDL3u`WacjD?wA7DcJ+T3mxirVzy>Uo(D@M{ zDR%~mc@kgN3gU6@Qu8YT0J~+VAx_xO5gjauj$Hv-V9#1h^Vr6sKG1Nlr}aSOv4bS_ zWW@j7^}|O`F7>)8yg7UhCdT`*@&*o9KvxQ^szM|moakLq7uny41j*SUiXZ-YlasSE z+U8jQ?d=omwPG|!Miktn$V79iY)L|(TOAmG2&8E@9_EW-@oT6cSpyU{K9c=P`|Ym4iTCbOH}jsJ2yQ8WP{`~ zR{A+1>ixnVl^SVSE(^Z+Z@0Z7Vbcqgb>L(G#!K8S07Qa1d=A)piHd(F;u#UpMSv@j z=||z8-Bg^0Gqfpzt}-X=j0x)Pv54-qEelh`xXw3wa+iQr&je5@@=-W%&{{uz7(vp2 z#rOZxxK23>k+b#J`5z{ubKsS;wQHF-B_)S$%wVRkg+XZ{nCLH?qV^T`rAK)p#D+72 zMUecufQyaSDQR?=_F=#!d-eX^xYcvb0#N$6HNZ>dVx9jW4LE%YPm?acrn3{fR|S#Z z;^qb{P|`t>;~}S*sI$`<<~O!PT)O6&L5QhJdQ~yL)fbP zI-O&34mf5om;tn*mrDpJC|~ZKMtz5mCT_@w=N;tq0)DPQuXA5pMF%6nu}2|>LkQ`- zw_S=o6^Q(psC`xhj|DBIU$yo?a^mrY?j0#oBlZAA=0Y}9pT&f;ZENf2halbIpjvna z;?Fq4smz|R%GFCe{1@KraYo0G07tSJECCd{q4zKeaAX0j-cEcb0-X0NDmuh&(r_f! zjwv?pU}1rk79}URFM~mgtXf=oF%~H0P*wzG)+C&|kcFTCpnrUmpF_kDLm&e^w_E$9GbtNejgl4Zp*om|J*Y@ z9M3VyjT-``)@yycF33s?f6%9jE@ZBJ5#^yij@A#L*#qr)>ZKCop%b61AGkD#V}fDj zJ?!$^iqyC5FaMU>ZHQL{S_oA$7SKF_kN7Mq)SKE1#Zyb6O&BQx1)ZG_MSv|ykOn0hVKCG21_1TO0lZ3tZpN)0pYV)%M@MAZVUmg{#P|FXr`lda$If z0q9&;7dVLk+58taEBz0E(8>XbfM>6nl4M4Z7xh^XnUHj$ODJG^>t>Z*2x1W=p*Re6v_@{mn?pCG8rbp2N+Z( z`EjDHgTvf%dCBb}CdP&8t+SPC>Cz%Zf}w(9bB4V+#rtJ6jj~}J!$4w4waM=B`tvmO z7ra5dBnjb{3Tzb(UdSH;8~jhEpiT{jlAj&G`4X~|sycv6PkAtXy*TXu|3Q&*=6?zeg!vvzz#+>A}qfybc%2!%gd>{<#7r;Pyp^!bUTV0nR!Jj;had;(mQcZrBc@29DB_~IM)^CAy3I*S#v$a z`M03!ny}T&_qMg#MGfH zyPLb#eoaP=--!1YZOD+Xy>6BQ;)0M^Vs9dd=`{hjQXeMFWyV zvMLvnCf!2O{lR|W$e0^bj5C_%m>v+&IB!(dlF=RZ5?-vmY_7Xb>*=(qG8QfDNPWAr zcUs#|=hj75!_xeI=N~!uNTvKg)_%1}ZS&+IC85DN=&wTBqA;ECu4AX}#nYX-#D$qd zAYbxKdcgf%E7Q<5NUiJ(2fe4t;2gW3wZCqh$**4R74IXObsiHbIg(dtvLo>%EBomY&E z%{-oS@}{PmEUF$huDNs9Xf&&I|Kn6uPmm6-e-W$%@BB(U#1znY7BMKh>gHllz9+FT ztd^#;^?A;(wEV9xl#hIJgN6{_Ka>`K zGj0C;H^b)nn+$qmOp=%kmcoVJl#o$LU`@coPt7t0CjIr%8QOcB{TC(1KD66@Gw#eP zb!mA2CgWzjXnpgd-WKe5MWGUv^jwA@Xl@oh*s$E#bxlvL+-@`rG>P`-I-hMl${0Je;{Wvpe}3%LNx{q0 zU`3C@L6#=N5tJ!%C2ODzHmH*}8ROXJF;{1#ynk$}b?B!%x3qp+X3pA>or8l}X*Zms zLPB4nJZ_IRjs_4Z;HVInf${mE!A9RV+{y8-714$HPmh z&e9RXG?dW_uCv&dW$zbdYxnmd8P$85De|!s>?Qol&tzNujQi~V5I<{paepyw-Ud*j zI~b2V`GI7CR*BGc!)Lj+tc3BQxG|W|wiyYc9zT?VI9(d<+%3x1Yb#YAnK0wR+>t~! zfNZMUhtKlDnK1gJG~rJ8`$Gy*Z=UBG;u`FuS7 zLyS$3+7YKGwpz)oEc#{6=qEBFK51=2{ElmBjf?fb7xEk_B_r+LTL zr;hK9l}z;|IXe;kP#98vG!W{RG;w|LTm8fBF4Ypmh0#-WX~v>+6~v`0QyWRmd~-rEQc0 z-hn#eq99j7Pp=C_^@-|_Tv5b_s{tNy$3D*5mmWvowp_nw;Cr(PxVp!j__quJV2TDmF8`k{kI0!&5qOHJlgo5wPeco7A?v@C-*k1g|h$hc(O=QwNf)#e-DCy z{0RKBIzt9ZP*h_4wpZ3X2Wc9K%Vl*JU?KN5I}vw@UU=Wky!n2flJjvJMb&u^pY7nr zr3e`6O(^d`ctM0ge+YnE6{05v(95H~K(A-FaZQlYKtGpP3-0RvqQ_fKjMf>96B)&>?&8GtZ zyB|}$OWwFCNY<<0p2c$ttn%L#)$5N2b7sgDrv0nTkCDEwe-!9WPa@~vn<{m)8@#)? z^^F+Ar_IP7{gx28VwTh8c79A~z+n-2Xo&a&$jfi_#g8&7H1AH^!|Z%5v%$qWD(Z7t zI0qK-4SnPSq+`LynXr62E5qYtL<>~VNL=Pd6Y^xWjAu5SRsaW&UmUy37qk z-YAq@g-~+$SCQAB*36wD8S1|B!WtPO>9Ew%1fhDOy+``LyaNu;kWQ@7G`0}GC7c_^ zB$fWcFbRag1>~@ky8wmkzcpnRu|4hgw04HfL%-X!`~Oe{M=|Uz8V3@LLn8%QFa9cN zUEUd51oQB2CLAGh-8<>CHHHOB&#^QxuGsRcr*w`NfbNca$=fQwKz%6`EuuQJ1)tU` z>4)(0I?%$91_V%q@ww*Mp)8pC34AYL;$Irgr-grJXdAO>TXpXJe3Rk+>6lP$Enn;M zCXAE=sGWdxy0+ud8>9Nfrw3|y{k_{2>Y`8G@~5v~c0Gto*J+LCDwjq|ThnydpSX2- z)7(PJA1o!b4$T!rZsGU5et|sK8H1ilLI6GYfV}y(pDgP7B>G(PYyW=zj*(KHRC|-r zr5mqUezuP_He!TJXtW={XHjvo`C|0zRZ_5)^ZqaOn6x-x$*1BYkH(5;1G=W*HbeXd zd$%w3v|Q$YnbphvPIEnJF6O;1!NdnFD*eT97O7Mo!%~phWx&8zreG=2@|F71&wEkD zm)%-D%}vB_`K!3nt(w~H97Y3t)i7*bTFS8(F-^JZ4cjQiGrGL+&dNKT@40j9*^W=2 zhg}1U1d2Shf+y(2!z!X+D$!EfrnBya5o~!GU+xaa%02w-lmc1IZJM)p7Hp@R&lPiG z)F||HVa)(N(9<)YcC-z)gFiEd9YG7n+UG3U+oGkpo}EmzH1;Th1>3-Gh_UhcaV!b; z0qiQ8x2LnA96&vuhY9V6mRW+Me5>Pu>48t=O?*8ZS|@eUENn&L|8Mc(tPubH5ayto zaevs^_DPg-ev3|<{04i=lb>S=Z#Y|`Ff7S>{f-8jmsr%&uLloZ_4UcduP&7ErMX3) zyVs@cVH8AB7{G+jp`Bp85UdC`sFuli@fzO8U|{-6(X_F#+4<+nn^J8J7n|aHzLmoJ zad=i=Qjxa{!-Cb&OW(H5Fu$`4{XGdM zShK+khgxB1P7(dwA;ih(a?DNqtPnSQyEDw7!Ev0=!uv#&-QdZB>t?@q=m+KSakdD- zQI6RXww80C^C?OMjsGWT5UyK80{qqsVBF{DtX}{Ib|bDhB3+|{L5^<}wLy~JR0knl zY}2!(%q3p8hoL_v`R7oQCpEtg*4N5F9^mz}cX^ z*gtk{NFh_;w97DL>S~B4P`!# zE~Na4yCd(S``$S}defkMxXbyAp-H~%l0-R%VJ;8#x4nVJ=0qSJVAB?)w}msr(=v?n zYY`1f*bwoVk&sb6uLb0pV}W#laf`5ZRy^=4KY<0I3p$^L?*p9fQc#e$KO$1bX@>+k zbIYAo;T=WPnmrOl?-(5$v~9^0vR=nn!C17J3miva*LTE0d0fp-$g$*p*zVki%P}sV zVV%Lf^$M{cj5#qZ?~JiVVd$f`NO6IrT)m5+hQ%L~5gh;>{{#Y}1_lPcUkB~C^38>- zCT#a~IJmbL)T{B|6&MDY1?QDJksswydpgF*0QJIGp?GB%V-tgqa?jYRM0165vvHDC zXuqDw8!gkuS3|e03`38p4U&}#a+qC|C^CBue0bf~EX*jj#5#sEZPkKe%cy05`G@lE zZz)AO=bT2prq@1-{m9F!r{)cMQa%bRKby7Ajh!92tK$B;tG;#Q--g*_e@}mN_P@$(k)|r3J&3#*)1_p%_{@jqFSIh_Q?<5v7H(jRqr;vCd=}BV+kp z@9yvYxc|EEzwWO;j>mZ%&CzFhzpwXoy`#tlvEf7i48Yp^q!_&S-kKx9i zjh4F2ld|W+ruV^tKmhvUN(L`-L4VQ*eOb@&jomg}iCUqaY}HUfZ7tMaKIWE}?|k8f z%s!=Mxn3_F(+l-R1UY&4>n~!?Z}nsK-MU%&`Rm!o>$tcS<;MTYhrwtc00EIUXaUHn z8tF{^mC5e>%=D|VFVIZc8utrHT)v!-0?2snGhzVDqBf~s2yNOT}zHT zE`IBEIWO4fIC| z&j4Wo`WQEuZNzy(-V3y`DOHyYpUy3XA}ua91{^~wGdKPHtE6wwfWXqu)3Xp6Xk^DM zy5B;hxf}YKs`t@HDM}s7ocP;NTc-<*@Pj`d|H&dAn}sfO2%X0v>X1-TQIT(bzVr|X zWg%65S~_{`i%WOcDH|@e!8Ri!wa^UhC&E1y@*m?sErOnB&Vh&x7PhuyZ%mK9kiUF7 zBcUIrA&ePBE-k;*uzs?y;x&@I`$BN9a^G`)^gA4c2%!fuV%v}3T>xa2H)w~s!VUmQ z%}@yhH1I`Q|L5XzmF@g5eZOD_eLf;i)K#YfI{@8vDj-^D9Y*WCpG;6mCe$5|)ZQpm zI&a;#2gx7Mm>v`<+xmt+Z?oS?urz}h=BzAF?Ec!B)0cK@m3b{o+6f;X=Dw6p^ioL2 z563CTYEe~JU56JjH-)8hI+Owo4!8(#xirUk1)lVzmT^sa%9ERC41`<6Jw-y{kN=;4 zdC=QDG6nAlS85V$oBz-6IUr_$uW96ng)iy@@XP-H{4Q0R^eJPn*4lC@I1GIjfJsY1 z>M1nT8(bb$DF+2+dENlhid7_ot_gbgtsWV8X6$PV(CgbOEOJV`|SXJf1jzzb146VZzEC}+2v)6~~2BxFOc{K$vLzJOFgv&9{A zvPd)R3p#q0-%H@IM%;+bvESt{IXX(KI@GDM#b2Wm5PywGVlDc?1UrY=T6jSWMwduU}hXG6Kz!ObCT9Tv3Zb`0}L2 zB~6<|x0AssAWZ*iEL-nAZoyv3Y>TaXs?vW;R~-y3(=f`3b4%1=!Lrc{I+$@^XL{n- zYTGVQHUm}q^6Q$}hb!zA`QDMV zV6cSv0h(0I+JO6@qt##s0tj~ureQJuiD*oE=NtSLBH_j} zdY*!$Z>eu^b+?a(QM}9+&yx4M89tDhmxiL!Y9X@(5D%iTZ$|rb9(WAByPt;K7NRyz zPEM|Hiy3|^DA)YB_zh!YV^=ik`}C7ek8F+=6%Jv-($^pOtsO;OIi(~gGiaR=TI$nf z7Y?X6euk7BdVAD%iT)A*uSA4W>fvc@{+EIS;DAz}6Pl%E_z-{g^W}jvNV*F+h7|oJ z+8#w$v&EnwVB^@IDggbt{Dl8lQi7T(kpW~O`I@4ipAsVLjnb1E0}5c;{@!!bdtgwV zMjEURFXvV&VpJ$$OS1!&~t25^z}3VD>7T zo&*xvSJ(MaN+EV3u(ukT*VD_O5+|&PF~}JP3W_)Wf$K>V#Q`V771GAH@a|POxXfmN zpf2k79e4=Z_flFmSbLcPll%`~kslck8SflX^q5~-Q40dl?5=RIw->H^)V-+$mf&N6 zzwY{Vp90Sf6-0ScC6VfeZ71n3&A>ep7=nh;8+wL8)fgDvw=GQ6w|9R~Nse*?D9tLh zEn56RdDiT`2VSS}z!mha+x~j%XcnKH8Gi4HW>Gr45G`il zBTVmq#GDY&>qtAm%Wze);Bp2kDRjTQ*N>)5ci*HI8m~sM%G!iYy}Q*VA9y~iDJYS7 zh>2U7w&c%i$rvb@4O1QX_bfM04~upxPb!8&ui%$1z}YrhsO zx{RCO@2>mnG*r{r5-DoH|GuHOP=e)7#Eo_>TL2T<{?X{YOd=Ri9Wq#BaisDR*b59! zsd!Qq*^5xjz7yn12DxPiq~W<1&ZimWCajaGZlKl-NcjRXgVU$An}jRy8D}`iP;WCfJrA=c zenH~QF@6MiC(&i=muHp5X<*PCmgU!nDCn7_@5Bl&#$p2zk=13cyQRqOHU)1i&E!qz3U@15@psQH)Gm9vh~ zLloNCh1l#G@HsGU)!fa|xYrHy^u4T&_bAO1-H&;^b+zf=LGd9Z%7By1&t%vqGmlX= zC;0@Qza)SXNr~zw!57e{kFj}Hh~?$EmG1a~<1fe}1@6GhKJ#wKthN>7#!-DQd?WnD zcdA`u4aBZ3Lm$-R*y}vp|RXC zR>WF8dLnVk`{~_;a9~g?AW3|6HhrI@hz9~31mP!G(Y!jw4QjEC;WGb`=_du zzr7o;(rWyNbjPFKy2Bcc4>&>htLG#II_ThzwiDxvs?y24e-ZD8eTgsoDU?=P6qn1h zXA^Q=LI7j0K)&k{4UGZV8~IcgFaC8pF@rTir}u&PWDn~a3fK`IDWdq-46+;)X~S?< za1@EHne!g;qQQn+q^he(OSKa4_=P0YfpQ$yBz@yMjP^3z@Ebs}o2t=|IzjZhemxu6 zjK3rstt0g4Yd4*>hPf50h_2+>a92SY+e`KyC1<)0NC_3!!+LAyJe1jwm!m^n4sp5& z*Ms@zD@haNTMac^C9f$EI>1o>OMEbUEH;r5R_?1o4S~EbLd~AoYkfMYsVD6AdYcBu z;h+^#q*$x4aVdQGmDZTM3%3a=FsIsxH2LaHE9p_k{Rvc9{BzyoiC=gsk?yRhmhfDF zEbcSV5Q-%0ptAGHvSRhZraY`5nOvUC(djZD6eOccXaO2!)moP1*54D$-Qpfe!{Fwr zfO_4Yv37Jz5&LX|RUYOOCxaOY#gX5K8;85cuS_C)l*%qFD9@rs+u8rTkt2C!C=Aiq zt#H;JsNbX_$S5s~?Q>~S4C4SD4RAZV18kChV77&r0ykB}}j&IK?Z9(Utl1V@=c}`Y| z>&=|9W-?sMJf7=Sk6GFjp%N;o5t%hVQJd8^nJc)%CE6;R|r@=R5}Kq7o+!8=j8_oY&`>V1K(Io>0{(;#=~$ueS$R3d`h2TOav>?K?sV zwfk=t=;4Cz_y|@N7Z-O=k{5i^-B15;Mb~3Y)rWF=g{eNP>SJRn$XN<6EINovtgo8f1CF1eLDyP6YStWZ0%I(6^_>ufvzsE>Uh08t4dH~M$suaPU&YB!a{tTXl@ z6x;1A%lnii7zAax!NF@A;Dn0vlf4Nx|3Y^$AIXJg5#af73hk7ScqgH77AgRsrYh=J z3b#(VkpJ`66i-~Dap%f@R!sS z3f1(d)#^;%7VGy1E`A+au&^=2553R`^nI%0m0qi*7HY+>t>XUU_cufR3H{uy<^v7Z z2@~bMhCfBHNfFcCuoSSCYd2C4ylA`ZD4uVW7p9MwLu?fLZavQ`{jq{cm-YFo?UlkS zI(XlT3JBj}x%BK2<*oFL?a5#qB+d0dT6%<8*27|ps->rAj&5}0zoa-}bIh|)=Gl=e zCrup9@*Y~qbbXyrNN?33VQ#}=PqFG>Mm9<=UKolL*AK?~wSYFTC;D3F79-CBtZLp- z$|H6q<10L?qhG|r8DV*~gx2q8Ft5?!c_cW*Dhmg`D*QGbm2F*EA@(lJ{w;VTe>kf< z?C16|#_KavmfmOq(Z(=%o8N|{=+>N7FcZ6#sSS9WB^X}wTj9nMI3X5&lbnl=8?AVE zu%3-m?nniTS#HRwO0+*Jl*tspOU4|!SQn|@3eWo__0#0ucC!nEjfEQQlK7%F@TG61 z{H)Ez1_?2<*?O;yiMU0;o%yvCI_bEF_t3aJJ}I4+q2N0hlf&`u?=x66FUr_djGT&j z`e4fH{9&6{Yw(D8HVp?1Kzg!-7)liyz$k4pYc$+hYxomoSw~$%jfK$<9YF$!8*|45 zXXAguX;aeZzwP~668bT&isJu{zU=UFm(7oO7#%IFk^TE#tj2 z!IqKr)1Xrl#d2FAQ1m(36z&(Jmpi@C8i*f#;EzJcaqblnLMECgp(DYjf@e3IoY@^9 z4!B!gkY#VjIs7?}qVTGsU=y|RO#R63#(D#{X3r-1NlHqZ^IaTNw;KdA^q|L6CXA7-OB?1NsM{j$kJvgSsmzJVZg=;j7SK3}qwSF3hi9J*g z{pS8c)2G)Dmu8Qabg(gHCq;Zr0_XVLL$c!mQX*I+in4?CnbD(yT|vf6ShwEQkD0mt z1>VWZHOY6N!{Qi-ntphKHcGP_T3z87EU#dof{P&#TF%j}%eEBy#Kr-Wk5~~(D{=L% zeUy+OXjggba#LpjU+r&(BA>Iu%vbxBLiano0>B?K;WsDV9_uB@2na4Sj0$Dt zJ6Cze0gWQ~#YEg^D3bXnAU+z3!dmjx3LPJmraau$6iCRmZ^?b|n6*hsLL!5F@B`tO zbqh+tG!U7E-U+kk1PmPncVV7!-pN2U&`QL|*;<)UrKn=#i(y@PW8X}MweDgqbZRDGrD@LJNz9&Z>>GmH? z;K>1JgHw}}bL5JGd7Mha0*rV@5Baab70^W!D^CB|+#K}2B#m(T z-IkZKw9gK1C;Zbr*i6f8b;ks!7;warl`zOf&-zr~EE+NroevXU@C(amCU4)M4o za{E;6E@efxL+w=k;5{p8T1rliZ!K62)(KsNYvUyThSG>M=QQ1WYsUkXqQl1~j+0sS_^5MPTNSoB zAnXo0NF63E&d45OrOkwRKhBIkR1W~zN^Qf}RnX8pklZP{{bq#QS^$Et%BH826Nca< zCNL9lEMIGEDdVZk{Jl=`&7sv}e%vqRz}%S8Xt~?84dUs)4bh)zAqz?(Ld^y?2GP;ERC5>o)Ap0afOK5B0K1Y=XZg zl{*M-aB7EmCv;IN&D9ay+HrZ&|AihggK>o8PalDk>YfW19H8ks#1vcOpkAv%|FWSn zd3`=i2Q~XQg2&jQX}|Q=d|_2ytxZs+Zvk$RlKb4430l@^u6wzgpap62M#)Wryk464 zAl7+0=WSPpa?H=bMd)J0^~MHp=NK3?UK6Nv#C=gZwxGG02R*)f~`wH*CxlO1@0POjP0Zd)Am|IdWGLQl{VYDC_RCN)9jo z&7eDQSAOoTHM!Pa%htAXQ}i=6Qj9ge)*u`J1CfDB8GLnz#}IfLjhgAw%18wG_f~g$y41z~Cs(v%eyZ=3W z^K(@5*>MMJ?#u{9Ym4DP#W6}^l-+w9G)%S8ZqXXLZDe`%3^TeebUrElc7>2&J5`gZ zlg^9*vr2=N?i^24RP)ZfJ7I@1x?dg*St@F0S2Y)pbAq7L?wb;xeMm}BpqISoBEF)V zlmFl%F=o{HoUf4ZBEm$Os88=oR=v0iCm53Q{OFcl0w4XMv0?AA18?P~inUIVhe{*Z zcu+dP+0xxH^zhT9=9O{I>E|9Sc1uRRkx*%_a?h^aAxoMR)<5x`#CtpkU4KV(J}hJR zn;s(gp6M)JJe=F(|LV7(Pwl8&ar$vL#F7)lTFXeIezH4s8dSRY|{RgEUhDC*CubSzUUERP-{{Hi&(Ua!s z=9St9tBXP3MQJvwbdCCk1UdT2wr>@3&;XZSjhUF7%&4!|iVK_{BR`Ded)vQ+0T`1z z_SJ^Fzo#1``AT5*a>Zbx!L__LBbRuNciMKHj0eF6Zo;n2UHo@DDsgm^Q=zM;5`KBY zD5aSM{doZxgte-hgRlGdAjX9EaFJ83kgy6iCQ03|nTVmx*5uCU85kHqOXyJC3Z$Fv z5f#O^rD(PE6xzEmX|?87KW_-87A!|@N>#Kjv=gZIW5&XT9K8!>@oZFiVtR>^u()S6w@;YKIFc{34omfe|$+)nSuXTZ@H)ash~x9h|zF-T($1sP{r3=sg@R|+3@9%><7PFPJaHU*VjJ9 zzXw*uEVd`~x7FN3%O45ncx{Rt_R2Zcd?@U#H=2u$i!(4aH5Jhc30j(aWt!3qp0JFp zEJ)H=8IdPsFi zmbRH6vzu>gD-<^2YBefepZW~NgVz>0R{K7`jpzC-9wxw*H3x_6-ulqc0>^6QgdtHi z?<1Na%eEw+_41*ct+}3keJo+{1Ac@|p(NO+K;u9&;ji^?Qx4d!Gdt?y;&NP3ahJO9 z)CWHsC1qv(z`#J@p96_Kq7w(4_TU*O5nGqR@@w%?KLFi!$iB|0+)Q`v-Pkbp*H~L jtNR*YT>4-6XKMXATBAbalHe61wCow4G0_7(ec1m1Qjio& literal 0 HcmV?d00001 diff --git a/examples/pipe_cht/pipe1d.py b/examples/pipe_cht/pipe1d.py new file mode 100644 index 0000000..696a783 --- /dev/null +++ b/examples/pipe_cht/pipe1d.py @@ -0,0 +1,218 @@ +"""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). + +Build: pythonfmu3 build -f pipe1d.py +""" +import math + +import numpy as np + +from pythonfmu3 import ( + Dimension, + Fmi3Causality, + Fmi3Initial, + Fmi3Slave, + Fmi3Variability, + Float64, + Unit, +) + +# Number of axial segments. Array sizes are fixed when the FMU is built: change this +# value and rebuild for a different resolution. +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" + + n = N_SEGMENTS + self.time = 0.0 + + # 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 + 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 + self.T_wall = np.full(n, 293.15) # K, wall temperature per segment + + # Outputs + self.T_fluid = np.full(n, self.T_initial) # K, bulk temperature per segment + self.h = np.zeros(n) # W/(m2.K) + self.q_wall = np.zeros(n) # W/m2, wall -> fluid + self.Q_wall = np.zeros(n) # W per segment, wall -> fluid + self.x = (np.arange(n) + 0.5) * self.length / n # m, segment centres + 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), + ]) + + segments = [Dimension(start=str(n))] + 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)) + + 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)) + + # --- FMI callbacks ----------------------------------------------------------------- + + def exit_initialization_mode(self): + n = N_SEGMENTS + 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(N_SEGMENTS, 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 = N_SEGMENTS + 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(N_SEGMENTS) + self.q_wall = self.h * (T_wall - self.T_fluid) + self.Q_wall = self.q_wall * self._perimeter * self.length / N_SEGMENTS + 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..19d5d12 --- /dev/null +++ b/pythonfmu3/tests/test_example_pipe_cht.py @@ -0,0 +1,161 @@ +"""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(**inputs): + pipe = Pipe1D(instance_name="test") + 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 * N_SEGMENTS))) + # 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_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")} + for name in ("T_wall", "T_fluid", "h"): + dims = variables[name].findall("Dimension") + assert [d.get("start") for d in dims] == [str(N_SEGMENTS)] + 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")} + + T_wall = np.linspace(300.0, 380.0, N_SEGMENTS) + 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"): + getattr(lib, f).argtypes = [C.c_void_p] + 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.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_SEGMENTS}), "h": get("h", {N_SEGMENTS}), "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(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 + 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) From a8bcc863370a8c858d4230ef09f6bbc41ecf08f3 Mon Sep 17 00:00:00 2001 From: Stephen Smith Date: Thu, 12 Nov 2026 21:17:31 +0000 Subject: [PATCH 2/2] pipe_cht: make the number of segments N a structural parameter N (UInt64, causality structuralParameter) sizes every per-segment array via . Setting it (in Configuration Mode, or on the Python class) reallocates the arrays through a property setter, so the resolution no longer requires editing the script and rebuilding. - cht_demo.py: --segments option; sets N through Configuration Mode for --fmu - README: N parameter and the configuration-mode call sequence - tests: N through the slave API, convergence with N against the analytic solution, model description dimensions, and the built FMU run through the FMI C API with N = 35 set in Configuration Mode Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01X8Wmw797DJpgrMFmT6n8CG --- examples/pipe_cht/README.md | 23 +++++++- examples/pipe_cht/cht_demo.py | 28 ++++++---- examples/pipe_cht/pipe1d.py | 67 +++++++++++++++++------ pythonfmu3/tests/test_example_pipe_cht.py | 53 +++++++++++++++--- 4 files changed, 132 insertions(+), 39 deletions(-) diff --git a/examples/pipe_cht/README.md b/examples/pipe_cht/README.md index 4d23986..46748ac 100644 --- a/examples/pipe_cht/README.md +++ b/examples/pipe_cht/README.md @@ -8,8 +8,8 @@ around it, is resolved by that solver, typically in 3D. Only the fluid is in the ## Coupling -The pipe is split into `N_SEGMENTS` equal axial segments. The FMU and the solid solver -exchange one value per segment, every communication step: +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 | | --- | --- | --- | @@ -33,6 +33,7 @@ One coupling step (what `cht_demo.py` does): | 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 | @@ -40,11 +41,26 @@ One coupling step (what `cht_demo.py` does): | `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 | -| `N_SEGMENTS` | 20 | – | number of axial segments: a constant at the top of `pipe1d.py`, so edit it and rebuild the FMU to change it | 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 `` 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 | @@ -80,6 +96,7 @@ 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) ``` diff --git a/examples/pipe_cht/cht_demo.py b/examples/pipe_cht/cht_demo.py index 9cb0c7d..c46cf90 100644 --- a/examples/pipe_cht/cht_demo.py +++ b/examples/pipe_cht/cht_demo.py @@ -3,7 +3,7 @@ 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 N axial segments matching the pipe, and M radial cells. +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) @@ -12,6 +12,7 @@ 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 @@ -26,7 +27,7 @@ import numpy as np sys.path.insert(0, str(Path(__file__).parent)) -from pipe1d import N_SEGMENTS, Pipe1D # noqa: E402 +from pipe1d import Pipe1D # noqa: E402 class SolidWall: @@ -80,8 +81,9 @@ def energy(self): class PythonPipe: """Drives the Pipe1D class directly, through the same calls an importer would make.""" - def __init__(self, **parameters): + 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() @@ -100,16 +102,20 @@ def do_step(self, t, dt): class FmuPipe: """Drives the built FMU through FMPy (pip install fmpy).""" - def __init__(self, fmu_path, **parameters): + 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: (N_SEGMENTS if v.dimensions else 1) 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) @@ -131,21 +137,23 @@ def main(): 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, **params) if args.fmu else PythonPipe(**params) + 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_SEGMENTS, T0=T_in) - third = N_SEGMENTS // 3 + 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_SEGMENTS) * third + heater_power = heat_flux * math.pi * d_outer * (length / n) * third h, T_fluid = pipe.get("h"), pipe.get("T_fluid") t = 0.0 @@ -171,7 +179,7 @@ def main(): 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_SEGMENTS / 2, x[2 * third - 1] + length / N_SEGMENTS / 2, + 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]") diff --git a/examples/pipe_cht/pipe1d.py b/examples/pipe_cht/pipe1d.py index 696a783..5a94142 100644 --- a/examples/pipe_cht/pipe1d.py +++ b/examples/pipe_cht/pipe1d.py @@ -22,6 +22,10 @@ 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 @@ -35,11 +39,11 @@ Fmi3Slave, Fmi3Variability, Float64, + UInt64, Unit, ) -# Number of axial segments. Array sizes are fixed when the FMU is built: change this -# value and rebuild for a different resolution. +# Default number of axial segments (the structural parameter N can override it). N_SEGMENTS = 20 RE_LAMINAR = 2300.0 @@ -82,8 +86,8 @@ def __init__(self, **kwargs): self.author = "PythonFMU3 examples" self.description = "1D pipe flow with wall heat transfer for CHT coupling with a solid solver" - n = N_SEGMENTS self.time = 0.0 + self._n = N_SEGMENTS # Parameters self.length = 2.0 # m @@ -94,17 +98,15 @@ def __init__(self, **kwargs): self.mu = 6.53e-4 # Pa.s self.T_initial = 293.15 # K, initial fluid temperature - # Inputs + # 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 - self.T_wall = np.full(n, 293.15) # K, wall temperature per segment - - # Outputs - self.T_fluid = np.full(n, self.T_initial) # K, bulk temperature per segment - self.h = np.zeros(n) # W/(m2.K) - self.q_wall = np.zeros(n) # W/m2, wall -> fluid - self.Q_wall = np.zeros(n) # W per segment, wall -> fluid - self.x = (np.arange(n) + 0.5) * self.length / n # m, segment centres + + # 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 @@ -125,13 +127,18 @@ def __init__(self, **kwargs): Unit(name="Pa.s", kg=1, m=-1, s=-1), ]) - segments = [Dimension(start=str(n))] 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)) @@ -159,10 +166,34 @@ def __init__(self, **kwargs): 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 = N_SEGMENTS + 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() @@ -189,13 +220,13 @@ 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(N_SEGMENTS, self.Nu * self.k_fluid / self.diameter) + 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 = N_SEGMENTS + 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 @@ -211,8 +242,8 @@ def _advance(self, T_old, dt): return T_new def _update_outputs(self): - T_wall = np.asarray(self.T_wall, dtype=float).reshape(N_SEGMENTS) + 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 / N_SEGMENTS + 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 index 19d5d12..e650f4a 100644 --- a/pythonfmu3/tests/test_example_pipe_cht.py +++ b/pythonfmu3/tests/test_example_pipe_cht.py @@ -20,8 +20,9 @@ from pipe1d import N_SEGMENTS, Pipe1D # noqa: E402 -def make_pipe(**inputs): +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(): @@ -44,12 +45,40 @@ def test_constant_wall_temperature_matches_analytic(m_dot): # 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 * N_SEGMENTS))) + 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)) @@ -86,9 +115,11 @@ def test_build_fmu(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")} - for name in ("T_wall", "T_fluid", "h"): + 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("start") for d in dims] == [str(N_SEGMENTS)] + assert [d.get("valueReference") for d in dims] == [n.get("valueReference")] assert variables["T_wall"].get("causality") == "input" assert variables["h"].get("causality") == "output" @@ -107,7 +138,8 @@ def test_fmu_matches_python_model(tmp_path): md = fromstring(z.read("modelDescription.xml")) vr = {v.get("name"): int(v.get("valueReference")) for v in md.find("ModelVariables")} - T_wall = np.linspace(300.0, 380.0, N_SEGMENTS) + 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 @@ -119,8 +151,9 @@ def test_fmu_matches_python_model(tmp_path): 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"): + 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] @@ -128,6 +161,9 @@ 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): @@ -143,7 +179,7 @@ def get(name, n): 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_SEGMENTS}), "h": get("h", {N_SEGMENTS}), "T_out": get("T_out", 1)}})) + 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) @@ -151,11 +187,12 @@ def get(name, n): import json result = json.loads(proc.stdout.strip().splitlines()[-1]) - reference = make_pipe(m_dot=0.03, T_in=290.0, T_wall=T_wall) + 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)