Lesson 4.3 - Engineering-Specific Charts

# --- [Electrical] Load-duration curve: every hour's load, sorted high -> low ---
np.random.seed(7)
# np.clip(a, lo, hi) caps values into a range; here just a lower bound of 800 kW
hourly_load = np.clip(np.random.normal(3800, 900, 8760)
                      + 700 * np.sin(np.arange(8760) / 8760 * 4 * np.pi), 800, None)
sorted_load = np.sort(hourly_load)[::-1]          # sort ascending, then [::-1] reverses it
pct_time = np.arange(1, 8761) / 8760 * 100        # x-axis: % of the year

fig, ax = plt.subplots(figsize=(7, 5))
ax.plot(pct_time, sorted_load, color="navy")
ax.fill_between(pct_time, sorted_load, alpha=0.15, color="navy")   # shade under the curve
ax.set_xlabel("Percent of year load is exceeded (%)")
ax.set_ylabel("Load (kW)")
ax.set_title("Load-Duration Curve")
plt.tight_layout(); plt.show()

# --- [Mechanical] Pump curve vs system curve: they cross at the operating point ---
Q = np.linspace(0, 0.03, 200)
pump_H = 45 - 6.0e5 * Q**2
system_H = 12 + 2.5e5 * Q**2
idx = np.argmin(np.abs(pump_H - system_H))        # index where the two curves are closest
fig, ax = plt.subplots(figsize=(7, 5))
ax.plot(Q * 1000, pump_H, label="Pump curve", color="darkred")     # *1000: m^3/s -> L/s
ax.plot(Q * 1000, system_H, label="System curve", color="steelblue")
ax.plot(Q[idx] * 1000, pump_H[idx], "ko")         # "ko" = black circle marker
ax.annotate("operating point", (Q[idx] * 1000, pump_H[idx]),
            textcoords="offset points", xytext=(10, 10))
ax.set_xlabel("Flow (L/s)"); ax.set_ylabel("Head (m)")
ax.set_title("Pump vs System Curve"); ax.legend()
plt.tight_layout(); plt.show()

# --- [Chemical/Mechanical] Moody diagram: friction factor vs Reynolds number ---
def colebrook_f(Re, rel_roughness, f_guess=0.02):
    f = f_guess
    for _ in range(30):                            # fixed-point iteration to solve for f
        f = (-2 * np.log10(rel_roughness / 3.7 + 2.51 / (Re * np.sqrt(f))))**-2
    return f

Re_range = np.logspace(3.7, 8, 100)               # log-spaced points from ~5e3 to 1e8
fig, ax = plt.subplots(figsize=(8, 6))
for rr in [0.05, 0.01, 0.001, 0.0001, 0.000001]:  # one curve per relative roughness
    ax.plot(Re_range, [colebrook_f(Re, rr) for Re in Re_range], label=f"eps/D = {rr}")
ax.set_xscale("log"); ax.set_yscale("log")        # both axes span orders of magnitude
ax.set_xlabel("Reynolds Number"); ax.set_ylabel("Friction Factor")
ax.set_title("Moody Diagram (Colebrook Equation)")
ax.legend(fontsize=8); ax.grid(True, which="both", alpha=0.3)   # which="both": major+minor gridlines
plt.tight_layout(); plt.show()

# --- Sankey energy balance: flow widths must sum to zero (in = out) ---
from matplotlib.sankey import Sankey
fig = plt.figure(figsize=(8, 5))
ax = fig.add_subplot(111, xticks=[], yticks=[])   # blank axes, no ticks
s = Sankey(ax=ax, unit="%", format="%.0f")
# positive flow = in, negative = out; orientations route each arm (0 = horizontal)
s.add(flows=[100, -55, -25, -12, -8],
      labels=["Grid input", "Process loads", "HVAC / cooling", "Lighting & aux", "Losses"],
      orientations=[0, 0, -1, 1, -1])
s.finish()
ax.set_title("Facility Energy Balance")
plt.show()

# Try the practice exercises here