#!/usr/bin/env python3 """Fan Wall model, version 2 (29 September 2026). Every number on https://localsolarsystem.com/initiatives/green/wind-energy-research/ comes from this file. Output: web/src/data/wind/model.json (read by the page explainers) and a copy in web/public/wind/fan-wall-model.json; the method and sources are in docs/research/FAN_WALL.md. Run: python3 docs/research/fan_wall_model.py What changed from version 1 (28 September 2026): * Rotor efficiency is no longer a guess: a blade-element momentum (BEM) design of the optimum rotor (Glauert wake rotation, Prandtl tip and hub loss, drag from a lift-to-drag envelope that falls with Reynolds number) gives the ceiling for each rotor size; a realisation factor turns the design value into what a first prototype usually measures. * Wind direction: fixed, non-yawing rotors only see the wind component normal to them; version 1 assumed the wind always blows straight in. * Mounting: a rotor flat on a closed facade has almost no air passing through it (continuity bound below); the improved design sits where air can pass: a free-standing crown on the roof edge, or a barrier. * Ducts: power gain per rotor area from the wind-lens literature, divided by the extra frontal area the duct takes, times an array interaction factor. * Second rotor (contra-rotating) and stator vanes: modest gains, bounded by multiple actuator-disc theory (Newman 1986: 16/25 = 0.64 for two discs). * Electrical chain built from generator, rectifier/MPPT and collection efficiencies for each generator size. * Benchmarks computed on the same building: rooftop PV, facade PV and a small conventional turbine on a mast. """ import json, math, pathlib # --------------------------------------------------------------------------- # Constants RHO = 1.225 # kg/m3, ISA sea level NU = 1.5e-5 # m2/s, kinematic viscosity of air at 15 C HOURS = 8760 BETZ = 16 / 27 NEWMAN_2DISC = 16 / 25 # Newman (1986): two actuator discs in series, same area CASES = ("low", "mid", "high") # pessimistic / central / optimistic # Reference building: 40 m x 20 m footprint, 15 floors of 3 m (45 m). BUILDING = {"length_m": 40.0, "depth_m": 20.0, "floors": 15, "floor_h_m": 3.0} BUILDING["height_m"] = BUILDING["floors"] * BUILDING["floor_h_m"] FACADE_M2 = BUILDING["length_m"] * BUILDING["height_m"] # 1,800 m2, one long face ROOF_M2 = BUILDING["length_m"] * BUILDING["depth_m"] # 800 m2 PERIMETER_M = 2 * (BUILDING["length_m"] + BUILDING["depth_m"]) # 120 m # --------------------------------------------------------------------------- # 1. Aerodynamics of the blade section: lift-to-drag ratio versus Reynolds. # Envelope of the best measured sections (cambered thin plates below about # Re 5e4, low-Re airfoils such as SG6043 or SD7003 above): Lissaman 1983, # Laitone 1997, Mueller and DeLaurier 2003, Giguere and Selig 1998. # Values are deliberately at or below the best published data (moulded # blades, dust and rain roughness): SG6043 reaches 59 at Re 1e5 clean and # about 45 at Re 3e5 with leading-edge roughness (Giguere and Selig 1998). LD_ENVELOPE = [(5e3, 5.0), (1e4, 8.0), (3e4, 14.0), (6e4, 22.0), (1e5, 35.0), (2e5, 50.0), (5e5, 70.0), (1e6, 90.0)] def lift_to_drag(re): pts = LD_ENVELOPE if re <= pts[0][0]: return pts[0][1] if re >= pts[-1][0]: return pts[-1][1] for (r0, l0), (r1, l1) in zip(pts, pts[1:]): if r0 <= re <= r1: f = (math.log(re) - math.log(r0)) / (math.log(r1) - math.log(r0)) return l0 + f * (l1 - l0) CL_DESIGN = 0.8 # design lift coefficient (cambered low-Re section) ALPHA_DESIGN = 5.0 # degrees, design angle of attack def glauert_cp(lam, n=2000): """Ideal rotor with wake rotation (Glauert), no drag, infinite blades.""" s = 0.0 for i in range(n): lr = (i + 0.5) / n * lam phi = 2 / 3 * math.atan(1 / lr) a = 1 / (1 + 4 * math.sin(phi) ** 2 / (4 * (1 - math.cos(phi)) * math.cos(phi))) ap = (1 - 3 * a) / (4 * a - 1) s += lr ** 3 * ap * (1 - a) * lam / n return 8 / lam ** 2 * s def bem_design(radius, blades, lam, v_rotor, hub_ratio=0.2, stations=40, swirl=True): """Optimum rotor with wake rotation (Manwell, McGowan and Rogers, Wind Energy Explained, 2nd ed., sec. 3.11). Returns Cp referred to the rotor swept area (including hub), and the blade chord/twist table. v_rotor: wind speed the rotor sees (m/s), used for the Reynolds number. swirl=False gives the same rotor with its wake swirl fully recovered (the upper bound for a perfect stator or contra-rotating stage).""" rh = hub_ratio cp = 0.0 table = [] dl = (1 - rh) / stations for i in range(stations): x = rh + (i + 0.5) * dl # r/R lr = lam * x # local speed ratio phi = 2 / 3 * math.atan(1 / lr) chord = 8 * math.pi * x * radius / (blades * CL_DESIGN) * (1 - math.cos(phi)) sigma = blades * chord / (2 * math.pi * x * radius) a = 1 / (1 + 4 * math.sin(phi) ** 2 / (sigma * CL_DESIGN * math.cos(phi))) ap = (1 - 3 * a) / (4 * a - 1) if swirl else 0.0 w = v_rotor * math.hypot(1 - a, lr * (1 + ap)) re = w * chord / NU ld = lift_to_drag(re) f_tip = 2 / math.pi * math.acos(math.exp(-(blades / 2) * (1 - x) / (x * math.sin(phi)))) f_hub = 2 / math.pi * math.acos(math.exp(-(blades / 2) * (x - rh) / (rh * math.sin(phi)))) f = f_tip * f_hub if swirl: # dCp = 8/lam^2 * F lr^3 a'(1-a) [1 - cot(phi)/(L/D)] dlr dcp = 8 / lam ** 2 * f * lr ** 3 * ap * (1 - a) * (1 - 1 / (ld * math.tan(phi))) * lam * dl else: # swirl recovered: annulus reaches the local Betz value 4a(1-a)^2 at a = 1/3 # at a = 1/3 the inflow angle is atan((2/3)/lr): drag loss 1.5 lr/(L/D) dcp = 2 * x * dl * f * BETZ * (1 - 1.5 * lr / ld) cp += dcp table.append({"r_over_R": round(x, 3), "chord_mm": round(chord * 1000, 1), "twist_deg": round(math.degrees(phi) - ALPHA_DESIGN, 1), "re": round(re), "ld": round(ld, 1)}) return max(cp, 0.0), table def best_rotor(diameter, v_ref, speedup=1.0, hub_ratio=0.2, blade_options=range(2, 9), lam_options=None): """Search blade count and design tip-speed ratio for the highest Cp.""" lam_options = lam_options or [x / 4 for x in range(3, 33)] # 0.75 .. 8 best = None for b in blade_options: for lam in lam_options: cp, tab = bem_design(diameter / 2, b, lam, v_ref * speedup, hub_ratio) if best is None or cp > best["cp"]: best = {"cp": cp, "blades": b, "tsr": lam, "table": tab} cp_ns, _ = bem_design(diameter / 2, best["blades"], best["tsr"], v_ref * speedup, hub_ratio, swirl=False) best["cp_swirl_free"] = min(cp_ns, BETZ) best["glauert_ideal"] = glauert_cp(best["tsr"]) return best # Realisation factor: first-generation hardware measured in a tunnel versus # its BEM design value (laminar separation, tip gap, manufacturing, yaw in # gusts). An assumption to be replaced by Phase 1 tests. REALISE = {"low": 0.65, "mid": 0.78, "high": 0.90} # Calibration check: SWEPT (Kishore, Coudron and Priya 2013, 0.39 m rotor) # measured Cp about 0.32 where this BEM gives 0.39: ratio about 0.8. For # centimetre rotors (Re below 20,000) measured values fall much further below # design (Howey et al. 2011: about 0.09 for a 2 cm rotor, BEM here about 0.2). REALISE_TINY = {"low": 0.40, "mid": 0.55, "high": 0.70} # --------------------------------------------------------------------------- # 2. Ducts (diffuser-augmented rotors), stators and a second rotor. # Power gain relative to a bare rotor of the same diameter, for a compact # brimmed diffuser ("wind lens"): Ohya and Karasudani 2010 report 2 to 5 for # long diffusers; compact designs that fit a facade module give less. DUCT_GAIN = {"low": 1.4, "mid": 1.8, "high": 2.3} # Outer (brim) diameter over rotor diameter for a compact wind lens, and the # square footprint each module takes in an array. DUCT_OUTER = 1.35 # Array interaction: ducts packed side by side cannot all draw extra air from # around them (an isolated duct gains partly by pulling in air from outside # its own frontal area). Unmeasured; Phase 1 and 2 must measure it. ARRAY = {"low": 0.70, "mid": 0.80, "high": 0.90} # Stator vanes behind the rotor: recover this fraction of the swirl loss, # minus the drag of an extra blade row. # No wind-turbine measurement exists; swirl-recovery vanes behind propellers # gave +2.6 percent (Li et al. 2017, AIAA 2017-3571; literature 2 to 5 percent). STATOR_RECOVERY = {"low": 0.15, "mid": 0.30, "high": 0.50} STATOR_DRAG = 0.02 # fraction of rotor power lost to the stator blade row # Contra-rotating second rotor: gain over a single rotor. Measured range 0 to # 60 percent; well-controlled recent tests cluster at 2 to 10 percent (review: # Adema et al. 2026, Wind Energ. Sci. 11, 2037). Newman's 0.64 caps the pair. CR_GAIN = {"low": 1.00, "mid": 1.06, "high": 1.15} # Physical ceiling for power per unit of FRONTAL area of a ducted system # (van Bussel 2007: exit-area Cp has not exceeded about Betz). CP_FRONTAL_CEILING = BETZ # --------------------------------------------------------------------------- # 3. Electrical chain (generator x rectifier/MPPT x collection and inverter). ETA = { # 120 mm fan motor run as generator at 0.2 to 3 W: tiny copper, diode or # low-voltage active rectification, thousands of connectors. "fan": {"gen": {"low": 0.55, "mid": 0.65, "high": 0.75}, "conv": {"low": 0.80, "mid": 0.88, "high": 0.93}, "grid": {"low": 0.92, "mid": 0.94, "high": 0.96}}, # 50 to 300 W permanent-magnet generator per module, MPPT DC stage, string inverter. "module": {"gen": {"low": 0.78, "mid": 0.85, "high": 0.90}, "conv": {"low": 0.92, "mid": 0.95, "high": 0.97}, "grid": {"low": 0.92, "mid": 0.95, "high": 0.96}}, # 5 to 10 kW conventional small turbine (for the benchmark). "turbine": {"gen": {"low": 0.85, "mid": 0.90, "high": 0.93}, "conv": {"low": 0.94, "mid": 0.96, "high": 0.97}, "grid": {"low": 0.95, "mid": 0.97, "high": 0.98}}, } def eta(kind, case): e = ETA[kind] return e["gen"][case] * e["conv"][case] * e["grid"][case] # --------------------------------------------------------------------------- # 4. Wind at the device. SITES = { "city-rooftop": {"mean": 4.0, "label": "City rooftop", "note": "Carbon Trust (2008): few urban building sites reach 5 m/s"}, "highway": {"mean": 4.5, "label": "Highway barrier, open land", "note": "ambient 4.5 m/s at 5 to 10 m height"}, "coastal": {"mean": 6.5, "label": "Exposed coast or hilltop", "note": "good small-wind sites: 6 to 7 m/s annual mean"}, } # Roof-edge speed-up for the crown: Mertens (2006, TU Delft thesis, after # Mertens 2003): 1.06 to 1.25 at an upwind roof edge, but 0.38 when the edge is # downwind (the direction factor below keeps only upwind directions). CROWN_SPEEDUP = {"low": 1.0, "mid": 1.1, "high": 1.2} # Turbulence and gust loss in a built environment (urban TI 20 to 40 percent). TURB = {"low": 0.85, "mid": 0.90, "high": 0.95} AVAIL = {"low": 0.92, "mid": 0.95, "high": 0.97} # Wind rose: von Mises concentration around the prevailing direction the # array faces (0 = wind equally from all directions). ROSE_KAPPA = {"low": 0.5, "mid": 1.0, "high": 2.0} YAW_EXP = {"bare": 3.0, "ducted": 2.0} # power ~ cos(yaw)^n for a fixed rotor def direction_factor(kappa, n, offset=0.0): """Mean of max(cos t, 0)^n over a von Mises wind rose centred on the prevailing direction; `offset` is the angle between the prevailing wind and the direction the array faces.""" i0 = sum((kappa / 2) ** (2 * k) / math.factorial(k) ** 2 for k in range(30)) s, m = 0.0, 3600 for i in range(m): t = -math.pi + (i + 0.5) * 2 * math.pi / m f = math.exp(kappa * math.cos(t)) / (2 * math.pi * i0) s += f * max(math.cos(t - offset), 0) ** n * 2 * math.pi / m return s def rayleigh(mean): c = mean / (math.sqrt(math.pi) / 2) return lambda v: (2 * v / c ** 2) * math.exp(-(v / c) ** 2) def aep_per_m2(cp_frontal_fn, mean, v_in, v_rated, v_out, eff, factor): """kWh per m2 of frontal area per year. cp_frontal_fn(v) -> Cp on frontal area.""" f = rayleigh(mean) e, dv, v = 0.0, 0.1, 0.05 while v < 30: if v_in <= v <= v_out: ve = min(v, v_rated) e += 0.5 * RHO * ve ** 3 * cp_frontal_fn(ve) * eff * f(v) * dv v += dv return e * HOURS * factor / 1000 def mean_power_density(mean): f = rayleigh(mean); e, v = 0.0, 0.05 while v < 30: e += 0.5 * RHO * v ** 3 * f(v) * 0.1; v += 0.1 return e # --------------------------------------------------------------------------- # 5. Configurations. FAN_FRAME, FAN_RAIL = 0.120, 0.008 FAN_D = FAN_FRAME * 0.96 FAN_HUB = 0.35 FAN_SWEPT_FRACTION = (math.pi / 4 * FAN_D ** 2 * (1 - FAN_HUB ** 2)) / (FAN_FRAME + FAN_RAIL) ** 2 # Measured Cp of stock cooling fans run as turbines is low (Howey et al. 2011: # about 9 percent for a purpose-made 2 cm rotor); stock fans are designed to push air. STOCK_FAN_CP = {"low": 0.04, "mid": 0.08, "high": 0.12} MODULE_D = 0.50 # improved module rotor diameter (m) MODULE_PITCH = MODULE_D * DUCT_OUTER + 0.05 # square footprint incl. frame (m) MODULE_HUB = 0.18 MODULE_RATED = 13.0 # m/s, generator sized for this wind speed REF_V = 6.0 # m/s, design wind speed for blade Reynolds numbers # Closed-facade bound: rotors flat on a solid wall. The air entering them must # leave through the gap between modules and wall, open only at the edges. def closed_facade_bound(cavity_m=0.3): outlet = cavity_m * (2 * BUILDING["height_m"] + BUILDING["length_m"]) # two sides and the top ratio = outlet / FACADE_M2 # mean through-velocity / wind speed (upper bound) # free-standing optimum: the rotor plane sees 2/3 of the wind speed return {"cavity_m": cavity_m, "outlet_m2": round(outlet, 1), "velocity_ratio_max": round(ratio, 4), "power_ratio_max": round(min(1.0, (ratio / (2 / 3)) ** 3), 7)} # Costs, USD. Per m2 of frontal area unless stated. COST = { "fans": { # version 1 figures, per m2 of wall "fan_each": {"low": 3.0, "mid": 5.0, "high": 8.0}, "frame_mesh": {"low": 60, "mid": 100, "high": 160}, "electronics": {"low": 30, "mid": 60, "high": 100}, "install": {"low": 150, "mid": 300, "high": 600}, }, "module_each": { # one 0.5 m ducted module "rotor": {"low": 15, "mid": 25, "high": 40}, "generator": {"low": 40, "mid": 70, "high": 120}, "duct_diffuser": {"low": 40, "mid": 70, "high": 120}, "stator": {"low": 5, "mid": 10, "high": 20}, "guard_mesh": {"low": 10, "mid": 15, "high": 25}, "power_board": {"low": 25, "mid": 40, "high": 70}, "frame_connectors": {"low": 20, "mid": 35, "high": 60}, "second_rotor_and_generator": {"low": 50, "mid": 90, "high": 150}, }, "module_bos_m2": { # crown structure, installation, inverters, engineering share "structure_install": {"low": 250, "mid": 450, "high": 800}, "inverter_bos": {"low": 40, "mid": 60, "high": 100}, "engineering_cert": {"low": 30, "mid": 60, "high": 120}, }, } OM = {"low": 0.02, "mid": 0.03, "high": 0.05} # of capex per year RATE, YEARS = 0.07, 20 CRF = RATE * (1 + RATE) ** YEARS / ((1 + RATE) ** YEARS - 1) # Electricity the building owner avoids buying, USD/kWh: US commercial 2024 # average 12.75 c (EIA Table 4); EU non-household about EUR 0.18-0.19 (Eurostat, # H2 2025, band ID, excl. VAT); EU households about EUR 0.29. PRICE = {"low": 0.13, "mid": 0.19, "high": 0.29} def lcoe(capex, aep, om): return (capex * CRF + capex * om) / aep if aep > 0 else None def payback(capex, aep, om, price): """Simple payback in years; None when it never pays back within 40 years.""" net = aep * price - capex * om pb = capex / net if net > 0 else None return pb if pb is not None and pb <= 40 else None # --------------------------------------------------------------------------- def main(): out = {"version": 2, "date": "2026-09-29", "building": BUILDING, "facade_m2": FACADE_M2, "roof_m2": ROOF_M2, "perimeter_m": PERIMETER_M} kappa_dir = {c: {k: round(direction_factor(ROSE_KAPPA[c], YAW_EXP[k]), 3) for k in YAW_EXP} for c in CASES} out["direction_factor"] = kappa_dir out["closed_facade"] = closed_facade_bound() # --- Rotor design study: best Cp by rotor diameter (bare, BEM) ---------- sizes = [0.04, 0.08, 0.12, 0.25, 0.5, 1.0, 2.0, 7.0] study = [] for d in sizes: hub = 0.35 if d <= 0.14 else 0.18 if d <= 1 else 0.1 r = best_rotor(d, REF_V, hub_ratio=hub) rd = best_rotor(d, REF_V, speedup=DUCT_GAIN["mid"] ** (1 / 3), hub_ratio=hub) # ducted: faster flow at rotor re75 = [row["re"] for row in r["table"] if row["r_over_R"] >= 0.72][0] study.append({"d_m": d, "blades": r["blades"], "tsr": r["tsr"], "cp_bem": round(r["cp"], 3), "cp_bem_ducted_flow": round(rd["cp"], 3), "cp_swirl_free": round(r["cp_swirl_free"], 3), "glauert": round(r["glauert_ideal"], 3), "re_75": re75, "ld_75": round(lift_to_drag(re75), 1)}) out["rotor_study"] = study # --- Improved module rotor: blade design table -------------------------- v_rotor = REF_V * DUCT_GAIN["mid"] ** (1 / 3) mod = best_rotor(MODULE_D, REF_V, speedup=DUCT_GAIN["mid"] ** (1 / 3), hub_ratio=MODULE_HUB, blade_options=range(3, 8), lam_options=[x / 4 for x in range(8, 25)]) fan = best_rotor(FAN_D, REF_V, hub_ratio=FAN_HUB) mod_bare = best_rotor(MODULE_D, REF_V, hub_ratio=MODULE_HUB, blade_options=range(3, 8), lam_options=[x / 4 for x in range(8, 25)]) out["module_rotor"] = {"d_m": MODULE_D, "blades": mod["blades"], "tsr": mod["tsr"], "cp_bem": round(mod["cp"], 3), "cp_swirl_free": round(mod["cp_swirl_free"], 3), "v_rotor_at_6": round(v_rotor, 2), "rpm_at_6": round(mod["tsr"] * v_rotor / (MODULE_D / 2) * 60 / (2 * math.pi)), "tip_speed_at_10": round(mod["tsr"] * 10 * DUCT_GAIN["mid"] ** (1 / 3), 1), "table": mod["table"][::4] + [mod["table"][-1]]} out["fan_rotor_redesign"] = {"d_m": round(FAN_D, 4), "blades": fan["blades"], "tsr": fan["tsr"], "cp_bem": round(fan["cp"], 3), "table": fan["table"][::4] + [fan["table"][-1]]} # --- Configurations ----------------------------------------------------- configs = {} def cp_module(case, stator=True, second=False): base = mod["cp"] * REALISE[case] # bare-rotor Cp, realised swirl_loss = max(0.0, mod["cp_swirl_free"] - mod["cp"]) * REALISE[case] cp_rotor = base * DUCT_GAIN[case] # on rotor area, in the duct if stator: cp_rotor += swirl_loss * STATOR_RECOVERY[case] * DUCT_GAIN[case] - STATOR_DRAG * cp_rotor if second: cp_rotor *= CR_GAIN[case] frontal = cp_rotor * (math.pi / 4 * MODULE_D ** 2 * (1 - MODULE_HUB ** 2)) / MODULE_PITCH ** 2 * ARRAY[case] return min(frontal, CP_FRONTAL_CEILING), cp_rotor def add(key, label, cp_frontal, kind, yaw, mount_factor, v_in, v_rated, capex_m2, note, extra=None): row = {"label": label, "note": note, "cp_frontal": {}, "aep_m2": {}, "cf": {}, "capex_m2": capex_m2, "lcoe": {}, "payback": {}, "power_m2_at": {}} for c in CASES: eff = eta(kind, c) fac = kappa_dir[c][yaw] * TURB[c] * AVAIL[c] # direction, turbulence, availability row["cp_frontal"][c] = round(cp_frontal[c], 3) p_rated = 0.5 * RHO * v_rated ** 3 * cp_frontal[c] * eff # W per m2 at rated wind, straight on row["power_m2_at"][c] = {str(v): round(0.5 * RHO * v ** 3 * cp_frontal[c] * eff, 1) for v in (6, 10)} row["aep_m2"][c] = {} row["cf"][c] = {} for s, site in SITES.items(): mean = site["mean"] * mount_factor[c] # speed-up at the device e = aep_per_m2(lambda v, k=c: cp_frontal[k], mean, v_in, v_rated, 25, eff, fac) row["aep_m2"][c][s] = round(e, 1) row["cf"][c][s] = round(e / (p_rated * HOURS / 1000), 4) if p_rated > 0 else 0 # LCOE: best = high perf + low cost, mid = mid + mid, worst = low + high for s in SITES: row["lcoe"][s] = {"best": round(lcoe(capex_m2["low"], row["aep_m2"]["high"][s], OM["low"]), 3), "mid": round(lcoe(capex_m2["mid"], row["aep_m2"]["mid"][s], OM["mid"]), 3), "worst": round(lcoe(capex_m2["high"], row["aep_m2"]["low"][s], OM["high"]), 3)} pb = payback(capex_m2["mid"], row["aep_m2"]["mid"][s], OM["mid"], PRICE["mid"]) row["payback"][s] = round(pb, 1) if pb else None pbb = payback(capex_m2["low"], row["aep_m2"]["high"][s], OM["low"], PRICE["high"]) row.setdefault("payback_best", {})[s] = round(pbb, 1) if pbb else None if extra: row.update(extra) configs[key] = row fan_cost = {c: FAN_SWEPT_FRACTION / (math.pi / 4 * FAN_D ** 2 * (1 - FAN_HUB ** 2)) * COST["fans"]["fan_each"][c] + COST["fans"]["frame_mesh"][c] + COST["fans"]["electronics"][c] + COST["fans"]["install"][c] for c in CASES} fan_cost = {c: round(v) for c, v in fan_cost.items()} one = {c: 1.0 for c in CASES} add("fans-v1", "A. Version 1: stock 120 mm fans, free-standing wall", {c: STOCK_FAN_CP[c] * FAN_SWEPT_FRACTION for c in CASES}, "fan", "bare", one, 3.0, 15.0, fan_cost, "Stock cooling fans run backwards, 61 per m2, as if air could pass through the wall.") add("fans-redesign", "A2. Purpose-designed 120 mm rotors (best case for tiny rotors)", {c: fan["cp"] * REALISE_TINY[c] * FAN_SWEPT_FRACTION for c in CASES}, "fan", "bare", one, 3.0, 15.0, {c: fan_cost[c] + {"low": 20, "mid": 40, "high": 70}[c] for c in CASES}, "Same wall with custom blades optimised for Re about 10,000 (BEM design).") per_m2 = 1 / MODULE_PITCH ** 2 mcost = {c: sum(v[c] for k, v in COST["module_each"].items() if k != "second_rotor_and_generator") for c in CASES} bos = {c: sum(v[c] for v in COST["module_bos_m2"].values()) for c in CASES} crown_cost = {c: round(mcost[c] * per_m2 + bos[c]) for c in CASES} crown_cost2 = {c: round((mcost[c] + COST["module_each"]["second_rotor_and_generator"][c]) * per_m2 + bos[c]) for c in CASES} cpB = {c: cp_module(c)[0] for c in CASES} cpB2 = {c: cp_module(c, second=True)[0] for c in CASES} open_pitch = MODULE_D + 0.06 open_frac = math.pi / 4 * MODULE_D ** 2 * (1 - MODULE_HUB ** 2) / open_pitch ** 2 cpB0 = {c: min(CP_FRONTAL_CEILING, mod_bare["cp"] * REALISE[c] * open_frac * ARRAY[c]) for c in CASES} add("crown-open", "B0. 0.5 m open rotors behind a guard mesh, roof-edge crown", cpB0, "module", "bare", CROWN_SPEEDUP, 3.0, MODULE_RATED, {c: round((mcost[c] - COST["module_each"]["duct_diffuser"][c] - COST["module_each"]["stator"][c] + 10) / open_pitch ** 2 + bos[c]) for c in CASES}, "Same rotor without the duct: more rotors fit per m2, but no speed-up, less yaw tolerance, no acoustic or safety shroud.") add("crown", "B. Improved: 0.5 m ducted modules with stator, roof-edge crown", cpB, "module", "ducted", CROWN_SPEEDUP, 3.0, MODULE_RATED, crown_cost, "Wind-lens duct, 5-blade rotor designed for Re about 100,000, stator vanes, 1.8 modules per m2, free-standing above the parapet.", {"cp_rotor_area": {c: round(cp_module(c)[1], 3) for c in CASES}, "modules_per_m2": round(per_m2, 2), "module_cost": mcost}) add("crown-cr", "B2. As B plus a contra-rotating second rotor", cpB2, "module", "ducted", CROWN_SPEEDUP, 3.0, MODULE_RATED, crown_cost2, "Second rotor and generator behind the first, turning the other way.", {"cp_rotor_area": {c: round(cp_module(c, second=True)[1], 3) for c in CASES}}) add("barrier", "B3. Same modules as a free-standing barrier (no roof speed-up)", cpB, "module", "ducted", one, 3.0, MODULE_RATED, {c: round(mcost[c] * per_m2 + {"low": 180, "mid": 320, "high": 550}[c] + COST["module_bos_m2"]["inverter_bos"][c] + COST["module_bos_m2"]["engineering_cert"][c]) for c in CASES}, "Highway or coastal barrier: the wall structure exists or is cheaper than a roof crown.") out["configs"] = configs # --- Benchmarks on the same building ----------------------------------- pv = {} PV_KWP_M2 = 0.21 # module power density PV_YIELD = {"low": 950, "mid": 1150, "high": 1400} # kWh per kWp per year, tilted roof array PV_FACADE = {"low": 0.55, "mid": 0.62, "high": 0.70} # vertical south facade vs optimal tilt # NREL/DOE 2025 Q1 benchmark: commercial rooftop $1.95-1.98/W, residential $2.78-2.95/W. PV_COST_W = {"roof": {"low": 1.4, "mid": 1.95, "high": 2.9}, "facade": {"low": 2.2, "mid": 3.0, "high": 4.0}} PV_OM = 0.01 for kind in ("roof", "facade"): row = {"aep_m2": {}, "capex_m2": {}, "lcoe": {}, "payback": None} for c in CASES: y = PV_YIELD[c] * (PV_FACADE[c] if kind == "facade" else 1) * PV_KWP_M2 * 0.93 # 0.5%/yr degradation, 20-yr mean row["aep_m2"][c] = round(y, 1) row["capex_m2"][c] = round(PV_COST_W[kind][c] * 1000 * PV_KWP_M2) row["lcoe"] = {"best": round(lcoe(row["capex_m2"]["low"], row["aep_m2"]["high"], PV_OM), 3), "mid": round(lcoe(row["capex_m2"]["mid"], row["aep_m2"]["mid"], PV_OM), 3), "worst": round(lcoe(row["capex_m2"]["high"], row["aep_m2"]["low"], PV_OM), 3)} row["payback"] = round(payback(row["capex_m2"]["mid"], row["aep_m2"]["mid"], PV_OM, PRICE["mid"]), 1) row["cf"] = {c: round(PV_YIELD[c] * (PV_FACADE[c] if kind == "facade" else 1) * 0.93 / HOURS, 3) for c in CASES} pv[kind] = row out["pv"] = pv # Small conventional turbine: 7 m rotor, about 10 kW, on an 18 to 30 m mast # at the same site (not on the building). Cost: PNNL Distributed Wind reports # 2025 edition ($6,680/kW in 2024) and 2026 edition ($9,630/kW in 2025; ten-year # average $8,910/kW); fleet capacity factor 14 to 16 percent. tur = {"rotor_d_m": 7.0, "rated_kw": 10.0, "cp": {"low": 0.30, "mid": 0.35, "high": 0.40}, "cost_per_kw": {"low": 6680, "mid": 8910, "high": 9630}, "aep": {}, "cf": {}, "lcoe": {}, "payback": {}} area = math.pi / 4 * 7.0 ** 2 for s, site in SITES.items(): tur["aep"][s] = {}; tur["cf"][s] = {} for c in CASES: f = rayleigh(site["mean"] * 1.15) # mast height versus rooftop reference: +15 percent mean wind e, v = 0.0, 0.05 while v < 25: if v >= 3: e += min(0.5 * RHO * area * v ** 3 * tur["cp"][c] * eta("turbine", c), tur["rated_kw"] * 1000) * f(v) * 0.1 v += 0.1 e = e * HOURS * AVAIL[c] / 1000 tur["aep"][s][c] = round(e) tur["cf"][s][c] = round(e / (tur["rated_kw"] * HOURS), 3) cap = {c: tur["cost_per_kw"][c] * tur["rated_kw"] for c in CASES} tur["lcoe"][s] = {"best": round(lcoe(cap["low"], tur["aep"][s]["high"], OM["low"]), 3), "mid": round(lcoe(cap["mid"], tur["aep"][s]["mid"], OM["mid"]), 3), "worst": round(lcoe(cap["high"], tur["aep"][s]["low"], OM["high"]), 3)} pb = payback(cap["mid"], tur["aep"][s]["mid"], OM["mid"], PRICE["mid"]) tur["payback"][s] = round(pb, 1) if pb else None out["small_turbine"] = tur # --- Whole-building view ------------------------------------------------ crown_m2 = PERIMETER_M * 1.5 # 1.5 m tall crown around the whole roof edge # Per-m2 yields above are for modules facing the prevailing wind. Around a # whole roof, the long windward side faces it, the long leeward side faces # away and the two short sides face across it. k, n = ROSE_KAPPA["mid"], YAW_EXP["ducted"] sides = [(BUILDING["length_m"], 0.0), (BUILDING["length_m"], math.pi), (BUILDING["depth_m"], math.pi / 2), (BUILDING["depth_m"], -math.pi / 2)] ring = sum(L * direction_factor(k, n, off) for L, off in sides) / (PERIMETER_M * direction_factor(k, n)) out["whole_building"] = { "crown_area_m2": crown_m2, "crown_ring_factor": round(ring, 3), "crown_aep_kwh": {s: round(configs["crown"]["aep_m2"]["mid"][s] * crown_m2 * ring) for s in SITES}, "crown_windward_only_aep_kwh": {s: round(configs["crown"]["aep_m2"]["mid"][s] * BUILDING["length_m"] * 1.5) for s in SITES}, "crown_windward_only_capex": round(configs["crown"]["capex_m2"]["mid"] * BUILDING["length_m"] * 1.5), "crown_capex": round(configs["crown"]["capex_m2"]["mid"] * crown_m2), "fan_wall_aep_kwh": {s: round(configs["fans-v1"]["aep_m2"]["mid"][s] * FACADE_M2) for s in SITES}, "fan_wall_capex": round(configs["fans-v1"]["capex_m2"]["mid"] * FACADE_M2), "roof_pv_m2": round(ROOF_M2 * 0.6), "roof_pv_kwp": round(ROOF_M2 * 0.6 * PV_KWP_M2), "roof_pv_aep_kwh": round(ROOF_M2 * 0.6 * pv["roof"]["aep_m2"]["mid"]), "roof_pv_capex": round(ROOF_M2 * 0.6 * pv["roof"]["capex_m2"]["mid"]), } out["mean_power_density_w_m2"] = {s: round(mean_power_density(site["mean"]), 1) for s, site in SITES.items()} # --- Power curves for the chart (whole reference installation, kW) ------- speeds = list(range(3, 16)) curves = {} def curve(key, area_m2, v_rated, kind): cfg = configs[key] return {c: [round(0.5 * RHO * min(v, v_rated) ** 3 * cfg["cp_frontal"][c] * eta(kind, c) * area_m2 / 1000, 2) for v in speeds] for c in CASES} curves["fans-v1"] = {"area_m2": FACADE_M2, **curve("fans-v1", FACADE_M2, 15.0, "fan")} curves["fans-redesign"] = {"area_m2": FACADE_M2, **curve("fans-redesign", FACADE_M2, 15.0, "fan")} curves["crown"] = {"area_m2": crown_m2, **curve("crown", crown_m2, MODULE_RATED, "module")} curves["crown-cr"] = {"area_m2": crown_m2, **curve("crown-cr", crown_m2, MODULE_RATED, "module")} # Per square metre, straight-on wind (W/m2): the fair comparison. per_m2c = {} for key, kind, vr in (("fans-v1", "fan", 15.0), ("fans-redesign", "fan", 15.0), ("crown", "module", MODULE_RATED), ("crown-cr", "module", MODULE_RATED)): per_m2c[key] = {c: [round(0.5 * RHO * min(v, vr) ** 3 * configs[key]["cp_frontal"][c] * eta(kind, c), 1) for v in speeds] for c in CASES} out["speeds"] = speeds out["power_curve_w_m2"] = per_m2c out["power_curve_kw"] = curves # --- What would have to be true: LCOE for the crown by frontal Cp and cost targets = [] for site in ("coastal", "highway"): for cp in (0.15, 0.20, 0.25, 0.30): for cm2 in (300, 500, 800, 1100): e = aep_per_m2(lambda v: cp, SITES[site]["mean"] * CROWN_SPEEDUP["mid"], 3.0, MODULE_RATED, 25, eta("module", "mid"), kappa_dir["mid"]["ducted"] * TURB["mid"] * AVAIL["mid"]) targets.append({"site": site, "cp_frontal": cp, "cost_m2": cm2, "aep_m2": round(e, 1), "lcoe": round(lcoe(cm2, e, 0.02), 3)}) out["targets"] = targets # --- Benchmarks (USD per kWh) for the LCOE chart ------------------------- out["benchmarks"] = [ {"label": "Utility-scale onshore wind", "low": 0.037, "high": 0.099, "src": "Lazard LCOE+ June 2026"}, {"label": "Utility-scale solar PV", "low": 0.040, "high": 0.098, "src": "Lazard LCOE+ June 2026"}, {"label": "Commercial and industrial solar", "low": 0.088, "high": 0.197, "src": "Lazard LCOE+ June 2026"}, ] out["households"] = {"eu_kwh": 3600, "us_kwh": 10791} out["assumptions"] = { "rho": RHO, "nu": NU, "ld_envelope": LD_ENVELOPE, "cl_design": CL_DESIGN, "realisation": REALISE, "duct_gain": DUCT_GAIN, "duct_outer_ratio": DUCT_OUTER, "array_factor": ARRAY, "stator_recovery": STATOR_RECOVERY, "stator_drag": STATOR_DRAG, "contra_rotating_gain": CR_GAIN, "newman_two_disc": NEWMAN_2DISC, "betz": BETZ, "eta": ETA, "eta_total": {k: {c: round(eta(k, c), 3) for c in CASES} for k in ETA}, "crown_speedup": CROWN_SPEEDUP, "turbulence": TURB, "availability": AVAIL, "rose_kappa": ROSE_KAPPA, "yaw_exponent": YAW_EXP, "stock_fan_cp": STOCK_FAN_CP, "module_d_m": MODULE_D, "module_pitch_m": round(MODULE_PITCH, 3), "module_rated_ms": MODULE_RATED, "costs": COST, "om": OM, "discount_rate": RATE, "life_years": YEARS, "crf": round(CRF, 4), "price": PRICE, "pv": {"kwp_m2": PV_KWP_M2, "yield": PV_YIELD, "facade_ratio": PV_FACADE, "cost_w": PV_COST_W, "om": PV_OM}, "sites": SITES, } root = pathlib.Path(__file__).resolve().parents[2] txt = json.dumps(out, indent=1) (root / "web/src/data/wind/model.json").write_text(txt) (root / "web/public/wind/fan-wall-model.json").write_text(txt) (root / "web/public/wind/fan_wall_model.py").write_text(pathlib.Path(__file__).read_text()) # --- Console summary ------------------------------------------------------- print("direction factors", kappa_dir) print("closed facade", out["closed_facade"]) print("rotor study (bare, BEM ceiling at 6 m/s):") for r in study: print(" ", r) print("module rotor", {k: v for k, v in out["module_rotor"].items() if k != "table"}) for row in out["module_rotor"]["table"]: print(" ", row) print("fan redesign", {k: v for k, v in out["fan_rotor_redesign"].items() if k != "table"}) for k, r in configs.items(): print(f"\n{r['label']}\n Cp frontal {r['cp_frontal']} capex/m2 {r['capex_m2']}") if "cp_rotor_area" in r: print(" Cp rotor area", r["cp_rotor_area"]) for s in SITES: print(f" {s:13} kWh/m2 {r['aep_m2']['low'][s]:7.1f} {r['aep_m2']['mid'][s]:7.1f} {r['aep_m2']['high'][s]:7.1f} CF mid {r['cf']['mid'][s]:.3f} LCOE {r['lcoe'][s]} payback {r['payback'][s]}") print("\nPV", pv) print("turbine", {k: tur[k] for k in ("aep", "cf", "lcoe", "payback")}) print("building", out["whole_building"]) print("power density", out["mean_power_density_w_m2"]) for t in targets: print(t) if __name__ == "__main__": main()