#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ ultra_navier.py — Système Foucaud–Claude in fabula, premier run. Modèle réduit d'un circuit hydraulique : un cœur périodique, un Y à fronce (hystérèse de Coandă), deux branches, une confluence qui dépose du grain à l'embouchure de la branche dominante. Horloge rapide = le battement ; horloge lente = le lit qui s'écrit entre deux battements. Dépendances : numpy, matplotlib. Rien d'autre. Sorties : ultra_navier_results.json (chiffres à lire, jamais à reconstruire) ultra_navier_run.png (partage au Y, rayons, bit s) ──────────────────────────────────────────────────────────────────────────── Les lignes maîtresses (mêmes lettres que dans le code) Segment (horloge rapide, Ohm + Navier + pile de gravité) : ΔP_i = 8 μ L_i / (π r_i^4 (1 + 4 b_i / r_i)) · Q_i − ρ g L_i sin α_i Cœur (Fourier, gain de Starling g) : P_in(t) = (1 + g (Q_ret − Q̄)/Q̄) · Σ_n a_n cos(n ω t + φ_n) Joint Y (fronce : partage + bit) : q = R_2/(R_1+R_2) + ε (2s − 1) s ← 1 − s si |R_2/(R_1+R_2) − 1/2| > h contre s Horloge lente (la Maine) : r_dom ← r_dom − k γ (v_dom − v_lent)_+ (le rapide freiné par le lent) Boucle : Q_in = Q_cœur + ρ_r Q_out(t − T_b) Tout bouton à zéro est présent pour être allumé plus tard, un à la fois. ──────────────────────────────────────────────────────────────────────────── """ import json import math import numpy as np # ═══════════════════════════════════ BOUTONS ══════════════════════════════════ # — Fluide (X) — MU = 1.0e-3 # viscosité [Pa·s] (eau ≈ 1e-3) RHO = 1000.0 # densité [kg/m³] GAMMA = 1.0 # grain : 0 = eau pure (aucune trace), 1 = chargée GAZ = 0.0 # fraction de gaz [ÉTEINT] T_INIT = 20.0 # température initiale [°C] [ÉTEINT] # — Cœur (Rodrigue) — P0 = 10.0 # pression moyenne [Pa] (bas exprès : on reste laminaire) A1, A2 = 0.6, 0.2 # amplitudes des deux harmoniques (fraction de P0) PHI2 = 0.8 # phase de la 2e harmonique [rad] — le cœur n'est pas une sinusoïde G_STARL = 0.0 # gain de Starling : 0 = cœur sourd [ÉTEINT] SUBSTEPS = 64 # échantillons par battement # — Espace (θ) : quatre segments — entrée, branche 1, branche 2, sortie — L_IN, R_IN = 0.30, 5.0e-3 # longueur [m], rayon [m] L_1, R_1 = 0.50, 5.2e-3 # branche 1 : un poil plus large → dominante au départ L_2, R_2 = 0.50, 5.0e-3 # branche 2 L_OUT, R_OUT = 0.30, 6.0e-3 # section aval (là où v_mix se calcule) B_SLIP = 0.0 # longueur de glissement de Navier [m] [ÉTEINT : Stokes] ALPHA = (0.0, 0.0, 0.0, 0.0) # angles [rad], 0 = horizontal [ÉTEINT] G_GRAV = 9.81 R_MIN = 1.0e-3 # un rayon ne descend pas sous ça (le banc ne bouche pas tout) # — Fronce du Y (Coandă) — EPS = 0.06 # décalage du partage vers la branche favorisée H = 0.05 # demi-largeur du pli : le bit tient tant que |r − ½| ≤ h S_INIT = 1 # bit initial : 1 = jet collé à la branche 1 # — Horloge lente (Hjulström réduit à une droite) — K_DEPOT = 5.0e-4 # [m par (m/s) par battement] : vitesse d'ensablement N_BEATS = 400 # nombre de battements simulés # — Boucle — RHO_RET = 0.0 # fraction qui revient à l'entrée (0 = robinet) [ÉTEINT] # ══════════════════════════════════════════════════════════════════════════════ def resistance(L, r, mu=MU, b=B_SLIP): """Poiseuille corrigé de Navier : 8 μ L / (π r⁴ (1 + 4b/r)).""" return 8.0 * mu * L / (math.pi * r**4 * (1.0 + 4.0 * b / r)) def emf(L, alpha, rho=RHO): """Pile de gravité le long de l'axe : ρ g L sin α (négatif si ça monte).""" return rho * G_GRAV * L * math.sin(alpha) def heart_pressure(t, omega=2.0 * math.pi): """Deux harmoniques, période 1 s. Toujours > 0 avec ces amplitudes.""" return P0 * (1.0 + A1 * math.cos(omega * t) + A2 * math.cos(2.0 * omega * t + PHI2)) def run(): r1, r2 = R_1, R_2 s = S_INIT hist = {"beat": [], "share_ohm": [], "share_real": [], "s": [], "r1": [], "r2": [], "v1": [], "v2": [], "v_mix": [], "depot": []} flips = [] A_in = math.pi * R_IN**2 A_out = math.pi * R_OUT**2 for n in range(1, N_BEATS + 1): # ── résistances avec la géométrie d'aujourd'hui ── R_i = resistance(L_IN, R_IN) Rb1 = resistance(L_1, r1) Rb2 = resistance(L_2, r2) R_o = resistance(L_OUT, R_OUT) R_par = Rb1 * Rb2 / (Rb1 + Rb2) R_tot = R_i + R_par + R_o E_tot = sum(emf(L, a) for L, a in zip((L_IN, L_1, L_2, L_OUT), ALPHA)) share_ohm = Rb2 / (Rb1 + Rb2) # Kirchhoff : part vers la branche 1 # ── fronce : le bit ne bascule que hors du pli, contre lui ── if s == 1 and share_ohm < 0.5 - H: s = 0; flips.append(n) elif s == 0 and share_ohm > 0.5 + H: s = 1; flips.append(n) share_real = min(1.0, max(0.0, share_ohm + EPS * (2 * s - 1))) # ── horloge rapide : un battement, échantillonné ── ts = np.linspace(0.0, 1.0, SUBSTEPS, endpoint=False) P = np.array([heart_pressure(t) for t in ts]) Q = (P + E_tot) / R_tot # débit total instantané [m³/s] Q1 = share_real * Q Q2 = (1.0 - share_real) * Q v1 = float(np.mean(Q1) / (math.pi * r1**2)) v2 = float(np.mean(Q2) / (math.pi * r2**2)) v_mix = float(np.mean(Q) / A_out) # ── horloge lente : la branche dominante s'ensable à sa propre embouchure ── # (v_mix est calculé et archivé, mais la chute qui fait lâcher le grain # est celle du rapide contre le lent — Kelvin-Helmholtz, pas la moyenne.) if v1 >= v2: depot = K_DEPOT * GAMMA * max(v1 - v2, 0.0) r1 = max(R_MIN, r1 - depot) else: depot = K_DEPOT * GAMMA * max(v2 - v1, 0.0) r2 = max(R_MIN, r2 - depot) for key, val in (("beat", n), ("share_ohm", share_ohm), ("share_real", share_real), ("s", s), ("r1", r1), ("r2", r2), ("v1", v1), ("v2", v2), ("v_mix", v_mix), ("depot", depot)): hist[key].append(val) return hist, flips def main(): hist, flips = run() # Reynolds moyen à l'entrée, pour vérifier qu'on a le droit d'être laminaire Q_mean = float(np.mean([(v * math.pi * r**2) for v, r in zip(hist["v1"], hist["r1"])])) \ + float(np.mean([(v * math.pi * r**2) for v, r in zip(hist["v2"], hist["r2"])])) v_in = Q_mean / (math.pi * R_IN**2) Re_in = RHO * v_in * (2 * R_IN) / MU results = { "parametres": {"EPS": EPS, "H": H, "K_DEPOT": K_DEPOT, "GAMMA": GAMMA, "R_1": R_1, "R_2": R_2, "P0": P0, "N_BEATS": N_BEATS}, "reynolds_entree_moyen": Re_in, "premiere_bascule_battement": flips[0] if flips else None, "toutes_les_bascules": flips, "nb_bascules": len(flips), "r1_final": hist["r1"][-1], "r2_final": hist["r2"][-1], "share_ohm_final": hist["share_ohm"][-1], } with open("ultra_navier_results.json", "w", encoding="utf-8") as f: json.dump(results, f, indent=2, ensure_ascii=False) print(f"Reynolds moyen à l'entrée : {Re_in:.0f} (laminaire si < ~2000)") if flips: print(f"Première bascule du Y au battement {flips[0]} ; {len(flips)} bascule(s) en {N_BEATS} battements.") print("Battements de bascule :", flips) else: print(f"Aucune bascule en {N_BEATS} battements (le pli h tient face à l'ensablement k).") print("→ ultra_navier_results.json écrit.") try: import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt fig, ax = plt.subplots(3, 1, figsize=(8, 8), sharex=True) b = hist["beat"] ax[0].plot(b, hist["share_ohm"], label="partage Kirchhoff (résistances)") ax[0].plot(b, hist["share_real"], label="partage réel (avec Coandă)") ax[0].axhline(0.5 + H, ls="--", lw=0.8); ax[0].axhline(0.5 - H, ls="--", lw=0.8) ax[0].set_ylabel("part vers branche 1"); ax[0].legend(loc="best") ax[1].plot(b, [r * 1e3 for r in hist["r1"]], label="r1 [mm]") ax[1].plot(b, [r * 1e3 for r in hist["r2"]], label="r2 [mm]") ax[1].set_ylabel("rayon d'embouchure"); ax[1].legend(loc="best") ax[2].step(b, hist["s"], where="post"); ax[2].set_ylabel("bit s"); ax[2].set_xlabel("battement") for fl in flips: for a in ax: a.axvline(fl, color="k", lw=0.5, alpha=0.4) fig.suptitle("Foucaud–Claude in fabula — premier run : cœur, fronce, dépôt") fig.tight_layout() fig.savefig("ultra_navier_run.png", dpi=120) print("→ ultra_navier_run.png écrit.") except ImportError: print("(matplotlib absent : pas de figure, le JSON suffit.)") if __name__ == "__main__": main()