WillTools 2026-09-09 -- clean snapshot of the development repo for review
ci / test (push) Successful in 6m18s

History-free snapshot of the private development repository at cbe45c5
(2026-09-09). Code, docs and public catalogue only; no plant or job
knowledge layers, no client or job names (tools/scrub.py clean). The
intro deck is left out of this snapshot because its slide images could
not be checked for job drawings.
This commit is contained in:
2026-09-09 08:32:15 -05:00
commit 8e6ee31596
406 changed files with 223435 additions and 0 deletions
+13
View File
@@ -0,0 +1,13 @@
"""beamcalc -- the CPython engine behind WillTools' Beam Analysis.
A precast beam (rectangular, L / spandrel, inverted-T; mild and/or
prestressed) on two bearings with optional overhangs, checked the way the
reference tool (Eriksson Beam v5.9) reports it: stripping stresses, in-place
strength (flexure, shear, torsion), in-place service (stresses, deflection)
and bearing. Input is one JSON job file (`beamcalc/1`, see `model.py`);
output is an `Analysis` and an HTML report. Nothing here imports Revit.
Units inside the engine: inches, kips, ksi. Feet appear only at the
boundary (the job file, the report).
"""
SCHEMA = 'beamcalc/1'
+3
View File
@@ -0,0 +1,3 @@
from .cli import main
import sys
sys.exit(main())
+70
View File
@@ -0,0 +1,70 @@
# -*- coding: utf-8 -*-
"""`python -m beamcalc job.json [-o report.html] [--json summary.json]`
Reads one `beamcalc/1` job, runs every check, writes the HTML report and
prints a one-line verdict. Exit 0 on PASS, 1 on FAIL, 2 when the job could
not be read (the message names the key)."""
from __future__ import annotations
import argparse
import json
import sys
from . import model, report, solve
def summary(a: solve.Analysis) -> dict:
worst_f = max(a.flex, key=lambda f: f.util) if a.flex else None
worst_s = max(a.shr, key=lambda s: s.r.util) if a.shr else None
out = {
'mark': a.beam.mark,
'ok': a.ok,
'stages': [{'name': s.name, 'ok': s.ok, 'util': s.util} for s in a.stages],
'flexure': {'util': worst_f.util, 'mu': worst_f.mu, 'phi_mn': worst_f.phi_mn,
'x_ft': worst_f.x_ft, 'combo': worst_f.combo} if worst_f else None,
'shear': {'util': worst_s.r.util, 'vu': worst_s.r.vu, 'phi_vn': worst_s.r.phi_vn,
'vc': worst_s.r.vc, 'x_ft': worst_s.x_ft,
'av_req_max': max(s.r.av_req for s in a.shr)} if worst_s else None,
'deflection': {k: a.deflection[k] for k in ('total_down', 'total_up', 'total_limit',
'live_down', 'live_up', 'live_limit')},
'reactions': {k: list(v) for k, v in a.reactions.items()},
'bearing': {k: {kk: vv for kk, vv in v.items() if kk != 'note'}
for k, v in a.bearing.items()},
'stripping': ({k: v for k, v in a.stripping.items() if k != 'result'}
if a.stripping else None),
'stations': [{'x_ft': f.x_ft, 'mu': f.mu, 'phi_mn': f.phi_mn, 'mcr12': f.mcr12,
'util': f.util, 'c': f.c, 'as': f.as_t, 'as_c': f.as_c, 'd': f.d,
'vu': s.r.vu, 'phi_vn': s.r.phi_vn, 'vc': s.r.vc, 'vs': s.r.vs_prov,
'bwd': s.r.bwd, 'av_req': s.r.av_req, 'av_prov': s.r.av_prov,
'f_act': s.r.f_act, 'f_alw': s.r.f_alw}
for f, s in zip(a.flex, a.shr)],
}
return out
def main(argv=None) -> int:
ap = argparse.ArgumentParser(prog='beamcalc')
ap.add_argument('job')
ap.add_argument('-o', '--out', help='HTML report path')
ap.add_argument('--json', help='write a JSON summary here')
args = ap.parse_args(argv)
try:
beam = model.load_job(args.job)
except (model.JobError, OSError, ValueError) as e:
print('COULD NOT READ JOB: {0}'.format(e), file=sys.stderr)
return 2
a = solve.analyse(beam)
if args.out:
with open(args.out, 'w', encoding='utf-8') as fh:
fh.write(report.render(a))
if args.json:
with open(args.json, 'w', encoding='utf-8') as fh:
json.dump(summary(a), fh, indent=1)
print('{0}: {1}'.format(beam.mark or 'beam', 'PASS' if a.ok else 'FAIL'))
for s in a.stages:
print(' {0:<20} {1} {2:.0%}'.format(s.name, 'PASS' if s.ok else 'FAIL', s.util))
return 0 if a.ok else 1
if __name__ == '__main__':
sys.exit(main())
+273
View File
@@ -0,0 +1,273 @@
# -*- coding: utf-8 -*-
"""Camber and deflections the way the reference tool tabulates them.
Reverse-engineered from RBA5 p.14 ("Camber and Deflections"), which prints
one row per load case per stage with a multiplier, then the total:
* **Multipliers are ACI 318-19 24.2.4.1's long-term factor, not PCI's.**
Each row's multiplier is 1 + ξ / (1 + 50·ρ') with ρ' = A's / (b·d) at
midspan. RBA5 (ρ' = 1.76 / (16 × 20) = 0.0055): ξ = 2.0 gives 2.57 and
ξ = 1.0 gives 1.78 -- exactly the printed 2.57 (SW, Dc, sustained live,
final) and 1.78 (SW at start of construction; topping, final). The
settings page's "Self Weight 1 / 1.3 / 2" and "Add'l Dead - / 0 / 2" are
therefore ξ per stage (erection / construction / end), and live load
always carries 1.00.
* **Stiffness is the effective moment of inertia per station**, ACI 318-19
24.2.3.5 (Bischoff): Ie = Icr / (1 - (2/3·Mcr/Ma)²·(1 - Icr/Ig)) where
Ma > 2/3·Mcr, else Ig -- with Ig and Icr both TRANSFORMED (the tool's
"Transform Reinforcement: All"), Mcr = fr·Sb on the transformed section,
and Ma the service moment of the STAGE's total load at that station.
At the start of construction only self weight is on, nothing cracks, and
SW deflects on the uncracked section with Eci (printed 0.029 / 1.78 =
0.0163" -- ours 0.0160"). At the final stage the same self weight
deflects on the cracked stiffness of the FULL service load (0.113 / 2.57
= 0.044").
* Self weight uses Eci when the setting says so; every other case Ec.
* Live load is patterned; the reported live deflection is the worst
pattern's, and the total takes the worst pattern too.
"""
from __future__ import annotations
import math
from typing import Dict, List, Optional, Sequence, Tuple
from . import loads as ld, section as sec_mod, statics
from .model import Beam
def cracked_inertia(beam: Beam, x_in: float, positive: bool = True,
ec_ksi: Optional[float] = None) -> Tuple[float, float, float]:
"""(Icr, kd, A_cr) of the cracked transformed section at x, elastic
(n·As tension, (n-1)·A's compression), solving the NA depth kd from
first moments. Positive bending: compression at the top."""
mats = beam.materials
ec = ec_ksi or mats.Ec
n_bar = mats.rebar.es_ksi / ec
n_str = mats.strand.ep_ksi / ec
sec = beam.section
h = sec.height_in
layers = [] # (depth from compression face, area, n)
for bx, by, area, db, e0, e1, row in beam.bars_at(x_in):
depth = (h - by) if positive else by
layers.append((depth, area, n_bar))
for g in beam.strands:
for sx, sy, area in g.strands(sec):
depth = (h - sy) if positive else sy
layers.append((depth, area, n_str))
if not any(d > 0 for d, a, nn in layers):
return sec.ixx, h, sec.area
block = sec.compression_area if positive else sec.compression_area_from_bottom
def first_moment(kd):
a, ybar = block(kd)
q = a * (kd - ybar) # concrete about the NA (above it)
for depth, area, nn in layers:
eq = (nn - 1.0) * area if depth < kd else nn * area
q += eq * (kd - depth) # positive above the NA
return q
lo, hi = 1e-4, h
for _ in range(100):
mid = (lo + hi) / 2.0
if first_moment(mid) > 0:
hi = mid
else:
lo = mid
kd = (lo + hi) / 2.0
a_c, ybar = block(kd)
# concrete: I about its own centroid + shift to the NA
icr = 0.0
top = h if positive else 0.0
for x0, y0, w, hh in sec.rects:
# part of this rect within the compression block
if positive:
lo_y = max(y0, h - kd)
hi_y = min(y0 + hh, h)
else:
lo_y = max(y0, 0.0)
hi_y = min(y0 + hh, kd)
if hi_y > lo_y:
d_h = hi_y - lo_y
centre = (lo_y + hi_y) / 2.0
depth_c = (h - centre) if positive else centre
icr += w * d_h ** 3 / 12.0 + w * d_h * (kd - depth_c) ** 2
a_cr = a_c
for depth, area, nn in layers:
eq = (nn - 1.0) * area if depth < kd else nn * area
icr += eq * (depth - kd) ** 2
a_cr += eq
return icr, kd, a_cr
def cracking_moment_kin(beam: Beam, x_in: float, positive: bool = True) -> float:
mats = beam.materials
fr = 7.5 * mats.lam * math.sqrt(mats.fc_ksi * 1000.0) / 1000.0
p = sec_mod.props_at(beam, x_in)
return fr * (p.sb if positive else p.st)
def effective_inertia(beam: Beam, x_in: float, ma_kin: float) -> float:
"""ACI 318-19 24.2.3.5 at one station for the service moment `ma`."""
positive = ma_kin >= 0
ig = sec_mod.props_at(beam, x_in).ixx
mcr = cracking_moment_kin(beam, x_in, positive)
ma = abs(ma_kin)
if ma <= (2.0 / 3.0) * mcr or mcr <= 0:
return ig
icr, kd, a_cr = cracked_inertia(beam, x_in, positive)
if icr >= ig:
return ig
denom = 1.0 - ((2.0 / 3.0) * mcr / ma) ** 2 * (1.0 - icr / ig)
return min(ig, icr / denom) if denom > 0 else ig
def long_term_multiplier(beam: Beam, xi: float) -> float:
"""1 + ξ / (1 + 50·ρ'), ρ' = compression steel ratio at midspan
(ACI 318-19 24.2.4.1.1 / 24.2.4.1.3)."""
if xi <= 0:
return 1.0
x_mid = beam.length_in / 2.0
h = beam.section.height_in
bars = beam.bars_at(x_mid)
# compression steel for positive bending = bars in the top half
as_c = sum(a for bx, by, a, db, e0, e1, row in bars if by > h / 2.0)
tens = [(a, h - by) for bx, by, a, db, e0, e1, row in bars if by <= h / 2.0]
d = (sum(a * dd for a, dd in tens) / sum(a for a, _ in tens)) if tens else 0.9 * h
b = beam.section.width_at(h - 1e-6) or beam.section.width_in
rho_c = as_c / (b * d) if b * d > 0 else 0.0
return 1.0 + xi / (1.0 + 50.0 * rho_c)
def _deflect(beam: Beam, loads, grid, ei) -> List[float]:
if not loads:
return [0.0] * len(grid)
xa, xb = beam.supports[0].x_in, beam.supports[1].x_in
res = statics.solve(loads, xa, xb, grid)
return statics.deflect(res, ei, xa, xb)
def analyse(beam: Beam, grid: Sequence[float]) -> dict:
s = beam.settings
mats = beam.materials
xa, xb = beam.supports[0].x_in, beam.supports[1].x_in
span = xb - xa
grid = list(grid)
n = len(grid)
e_sw = mats.Eci if s.ec_for_sw_defl.lower() == 'eci' else mats.Ec
ec = mats.Ec
sw = ld.self_weight(beam)
dnc = ld.case_loads(beam, 'deadnc')
dt = ld.case_loads(beam, 'topping')
dc = ld.case_loads(beam, 'deadc')
sus = ld.case_loads(beam, 'sus_live')
roof = ld.case_loads(beam, 'roof')
live_full = ld.case_loads(beam, 'live')
# multipliers (ξ per stage from the settings)
m_sw_start = long_term_multiplier(beam, s.sw_mult.get('erection', 1.0))
m_sw_final = long_term_multiplier(beam, s.sw_mult.get('end', 2.0))
m_dead_final = long_term_multiplier(beam, s.addl_dead_mult.get('end', 2.0))
m_top_final = long_term_multiplier(beam, s.sw_mult.get('erection', 1.0))
m_sus_final = m_dead_final
# --- start of construction: self weight (+ Dnc + Dt) on the uncracked
# section, judged for cracking under that load alone
start_loads = sw + dnc + dt
ma_start = statics.solve(start_loads, xa, xb, grid).moment
ie_start = [effective_inertia(beam, x, ma_start[i]) for i, x in enumerate(grid)]
d_sw_start = _deflect(beam, sw, grid, [e_sw * ie for ie in ie_start])
d_dnc_start = _deflect(beam, dnc, grid, [ec * ie for ie in ie_start])
d_dt_start = _deflect(beam, dt, grid, [ec * ie for ie in ie_start])
start_total = [m_sw_start * d_sw_start[i] + d_dnc_start[i] + d_dt_start[i]
for i in range(n)]
# --- final: worst live pattern for the total and for the live row
best = None
for pat in ld.patterns(beam):
live = ld.live_loads_for(beam, pat, 'live')
roof_p = ld.live_loads_for(beam, pat, 'roof')
# Two stiffnesses: the SUSTAINED rows (self weight, dead, topping,
# sustained live) deflect on Ie judged under the sustained load
# alone; the transient rows (live, roof) on Ie under the full
# service load. Measured on RBA5: with one Ie for everything the
# dead rows came out 5 % too soft (0.519 vs 0.495) while the live
# row was exact; Ie from the dead-only moment closes that gap.
sust_loads = sw + dnc + dt + dc + sus
final_loads = sust_loads + live + roof_p
ma_s = statics.solve(sust_loads, xa, xb, grid).moment
ma = statics.solve(final_loads, xa, xb, grid).moment
ie_s = [effective_inertia(beam, x, ma_s[i]) for i, x in enumerate(grid)]
ie = [effective_inertia(beam, x, ma[i]) for i, x in enumerate(grid)]
ei_sw = [e_sw * v for v in ie_s]
ei_s = [ec * v for v in ie_s]
ei_c = [ec * v for v in ie]
rows = {
'SW': (m_sw_final, _deflect(beam, sw, grid, ei_sw)),
'Dnc': (m_dead_final, _deflect(beam, dnc, grid, ei_s)),
'Dt': (m_top_final, _deflect(beam, dt, grid, ei_s)),
'Dc': (m_dead_final, _deflect(beam, dc, grid, ei_s)),
'Live': (1.0, _deflect(beam, live, grid, ei_c)),
'Sust. Live': (m_sus_final, _deflect(beam, sus, grid, ei_s)),
'Roof': (1.0, _deflect(beam, roof_p, grid, ei_c)),
}
total = [sum(m * d[i] for m, d in rows.values()) for i in range(n)]
sustained = [sum(m * d[i] for k, (m, d) in rows.items()
if k not in ('Live', 'Roof')) for i in range(n)]
live_row = [rows['Live'][1][i] + rows['Roof'][1][i] for i in range(n)]
key = max(abs(t) for t in total)
if best is None or key > best[0]:
best = (key, rows, total, sustained, live_row, ie, pat)
_, rows, total, sustained, live_row, ie, pat = best
idx = list(range(n)) # the tool reports over the whole length
# --- "prior / after installation of non-structural elements" split.
# RBA1/2/3/7 print the final stage in two blocks -- SW + Dnc (with their
# long-term multipliers) as "prior", Dt + Dc + live + sustained live +
# roof as "after" -- and judge the AFTER block against the total-load
# limit (RBA1: table total 0.423", reported "Defl. Total" 0.023" = the
# after block, limit L/480 -> PASS). RBA5 (composite, no Dnc load) and
# RBA6 ("Precast Only", Dnc present) print one block and judge the
# whole total. Inferred rule from those six packages: the split is
# made when the section behaves compositely AND a non-composite dead
# load exists. ACI 318-19 Table 24.2.2's "part of the total deflection
# occurring after attachment of nonstructural elements" is the idea.
composite = (s.composite_shear or '').lower() not in ('precast_only', 'precast only', 'none')
split = composite and any(abs(v) > 1e-12 for v in rows['Dnc'][1])
prior = [rows['SW'][0] * rows['SW'][1][i] + rows['Dnc'][0] * rows['Dnc'][1][i]
for i in range(n)]
after = [total[i] - prior[i] for i in range(n)]
judged = after if split else total
t_down = max(judged[i] for i in idx)
t_up = min(judged[i] for i in idx)
l_down = max(live_row[i] for i in idx)
l_up = min(live_row[i] for i in idx)
lim_t = span / s.total_defl_ratio
lim_l = span / s.live_defl_ratio
cracked = sum(1 for i in idx if ie[i] < sec_mod.props_at(beam, grid[i]).ixx - 1e-6)
return {
'span_in': span,
'total_down': t_down, 'total_up': -t_up, 'total_limit': lim_t,
'live_down': l_down, 'live_up': -l_up, 'live_limit': lim_l,
'total_util': max(t_down, -t_up) / lim_t if lim_t else 0.0,
'live_util': max(l_down, -l_up) / lim_l if lim_l else 0.0,
'total_ok': max(t_down, -t_up) <= lim_t + 1e-9,
'live_ok': max(l_down, -l_up) <= lim_l + 1e-9,
'rows': {k: {'multiplier': m, 'elastic': d} for k, (m, d) in rows.items()},
'total': total, 'sustained': sustained, 'live': live_row,
'prior': prior, 'after': after, 'split': split, 'judged': judged,
'table_total_down': max(total), 'table_total_up': -min(total),
'start_total': start_total, 'start_sw_multiplier': m_sw_start,
'ie': ie, 'pct_cracked': 100.0 * cracked / n if n else 0.0,
'multipliers': {'self_weight': m_sw_final, 'addl_dead': m_dead_final,
'topping': m_top_final, 'live': 1.0},
'e_sw': e_sw, 'e': ec,
'i': sec_mod.props_at(beam, (xa + xb) / 2.0).ixx,
'icr_mid': cracked_inertia(beam, (xa + xb) / 2.0)[0],
'pattern': pat,
'max_start': min(start_total) if start_total else 0.0,
'max_sustained': min(sustained) if sustained else 0.0,
}
+127
View File
@@ -0,0 +1,127 @@
# -*- coding: utf-8 -*-
"""Development of bars and strands along the beam.
A bar contributes only the stress it can develop at a station: fy times
(embedment / ld), capped at fy -- ACI 318-19 25.4.2.1 read as a linear
build-up, which is how the reference tool prints a partial `As` near the
ends (RBA5 p.12: As = 0.28 in² at 4" from the end for 4.4 in² of #6).
`ld` is Table 25.4.2.3 with ψt = ψe = ψg = 1. The tool's own setting "Bar
Confinement Provided" picks the column: unchecked (its default in every
package read) is the longer "other cases" length -- for a #6 at f'c = 5 ksi
that is 31.8", which is what RBA5's As build-up implies (2" of embedment ->
0.063 of the bar). Top-bar ψt = 1.3 is NOT applied: the printed d at the
end of RBA5 (12.80") is the area-weighted centroid of top and bottom bars
at EQUAL development, which ψt = 1.3 on the top bars would have broken.
"""
from __future__ import annotations
import math
from typing import Optional
from .model import Beam, BarRow, Materials
PSI_T_TOP = 1.0 # see module docstring; the reference tool's choice
def ld_bar(row: BarRow, mats: Materials, confined: bool = False,
mult: float = 1.0, top: bool = False) -> float:
"""Straight-bar TENSION development length (in) as the reference tool
prints it in its `As` column: fy·ψt·ψe·db / (20·λ·√f'c) for EVERY bar
size -- measured on RBA5 (#6, f'c 5 ksi: 31.8") and RBA7 (#7, f'c
8 ksi: 29.4"), both of which the size-split Table 25.4.2.3 would have
put in a different column. `confined` is kept for the day the tool's
"Bar Confinement Provided" box is understood; today it changes nothing.
"""
fy = mats.rebar.fy_ksi * 1000.0
fc = mats.fc_ksi * 1000.0
lam = mats.lam if mats.normalweight else 0.75
psi_t = PSI_T_TOP if top else 1.0
psi_e = 1.0
psi_g = psi_g_for(mats.rebar.fy_ksi)
ld = fy * psi_t * psi_e * psi_g / (20.0 * lam * math.sqrt(fc)) * row.db
return max(12.0, ld * mult)
def psi_g_for(fy_ksi: float) -> float:
"""ACI 318-19 Table 25.4.2.5 reinforcement grade factor: 1.0 for
Grade 40/60, 1.15 for Grade 80, 1.3 for Grade 100. Owner, 2026-09-06,
on how to treat Grade 80 development: "参考ACI吧". On RBA1 (Grade 80)
this closes part of the +72 % gap in the reported As against the
reference tool; the remaining ~x1.5 has no ACI text behind it and is
left as a recorded difference, not fitted (docs/beam-analysis-state.md)."""
if fy_ksi >= 100.0 - 1e-9:
return 1.3
if fy_ksi >= 80.0 - 1e-9:
return 1.15
return 1.0
def ld_capacity(row: BarRow, mats: Materials, mult: float = 1.0) -> float:
"""The length over which the reference tool lets a bar's FORCE build up
in the strength solve -- distinct from `ld_bar`, which only drives its
`As` column. Measured on RBA5: φMn = 36.7 k-ft at 2" of embedment and
148.3 at 8.8" fit a bar force rising linearly to fy over about 13.5",
which is ACI 318-19 25.4.9.2's COMPRESSION development length
max(0.02·fy·db/(λ√f'c), 0.0003·fy·db) = max(12.7, 13.5) for a #6 at
5 ksi. Whether that is the tool's intent or its accident is not known;
it is what its numbers do."""
fy = mats.rebar.fy_ksi * 1000.0
fc = mats.fc_ksi * 1000.0
lam = mats.lam if mats.normalweight else 0.75
db = row.db
ldc = max(0.02 * fy * db / (lam * math.sqrt(fc)), 0.0003 * fy * db)
return max(8.0, ldc * mult)
def ldh_bar(row: BarRow, mats: Materials, mult: float = 1.0) -> float:
"""Hooked-bar development length (in), ACI 318-19 25.4.3.1 with every
ψ = 1 (ψr, ψo, ψc, ψe) -- the conservative reading; a hooked bar end is
'developed' once it has ldh of embedment from the hook's outside."""
fy = mats.rebar.fy_ksi * 1000.0
fc = mats.fc_ksi * 1000.0
lam = mats.lam if mats.normalweight else 0.75
db = row.db
ldh = fy / (55.0 * lam * math.sqrt(fc)) * db ** 1.5
return max(8.0 * db, 6.0, ldh * mult)
def bar_fraction(embed_start: float, embed_end: float, row: BarRow,
mats: Materials, confined: bool = False, mult: float = 1.0,
top: bool = False, for_capacity: bool = False) -> float:
"""Developed fraction (0..1) of a bar at a station with the given
embedment towards each end. `for_capacity` picks the strength-solve
length (`ld_capacity`); otherwise the reported-`As` length."""
ld = ld_capacity(row, mats, mult) if for_capacity else \
ld_bar(row, mats, confined, mult, top)
f_start = embed_start / (ldh_bar(row, mats, mult)
if row.start_cond == 'hooked' else ld)
f_end = embed_end / (ldh_bar(row, mats, mult)
if row.end_cond == 'hooked' else ld)
return max(0.0, min(1.0, f_start, f_end))
def transfer_length_in(db: float) -> float:
"""ACI 318-19 25.4.8.1 / 22.5.9.1: lt = 50 db."""
return 50.0 * db
def strand_dev_length_in(db: float, fse_ksi: float, fps_ksi: float) -> float:
"""ACI 318-19 25.4.8.1: ld = (fse/3000) db + ((fps - fse)/1000) db (psi)."""
return (fse_ksi * 1000.0 / 3000.0) * db + ((fps_ksi - fse_ksi) * 1000.0
/ 1000.0) * db
def strand_available_ksi(embed: float, db: float, fse_ksi: float,
fps_ksi: float) -> float:
"""The stress a strand can carry at `embed` from its (bonded) end:
linear to fse over the transfer length, then linear to fps at ld."""
lt = transfer_length_in(db)
if embed <= 0:
return 0.0
if embed <= lt:
return fse_ksi * embed / lt
ld = strand_dev_length_in(db, fse_ksi, fps_ksi)
if embed >= ld:
return fps_ksi
return fse_ksi + (fps_ksi - fse_ksi) * (embed - lt) / (ld - lt)
+246
View File
@@ -0,0 +1,246 @@
# -*- coding: utf-8 -*-
"""Flexural strength by strain compatibility with a Whitney stress block.
At a station every bar and strand present is a layer with a developed
stress limit (`develop.py`); the neutral-axis depth c is found so that the
concrete block plus every compression layer equals every tension layer,
then Mn is the moment of all of it about the compression face. Positive
bending puts the block at the top; negative bending at the bottom -- the
same routine on the flipped section.
What is reported per station mirrors the reference tool's Flexure Summary
Table: c, εt, φ, φMn, the developed tension area `As`, the compression
area `A's`, `Aps`, and `d` (developed-area-weighted centroid of the
tension layers from the compression face).
Strand stress: PCI Design Handbook 8th Ed. power formula for 270 ksi
low-relaxation strand, εps = εse + εce(≈0) + 0.003·(dp - c)/c.
"""
from __future__ import annotations
import math
from dataclasses import dataclass, field
from typing import List, Optional, Tuple
from . import develop
from .model import Beam, Materials
@dataclass
class Layer:
depth: float # from the compression face (in)
area: float # in²
kind: str # 'bar' | 'strand'
f_limit: float # developable stress, ksi (bars: f·fy; strands: available fps)
es: float # modulus, ksi
fy: float = 60.0 # bars only
fse: float = 0.0 # strands only: effective prestress after losses
fpu: float = 270.0
dev: float = 1.0 # developed fraction (bars)
@dataclass
class FlexureResult:
phi_mn: float # kip-ft
mn: float # kip-ft
phi: float
c: float # in
eps_t: float
as_tension: float # developed tension bar area, in²
as_comp: float # bar area in compression, in²
aps: float # strand area, in²
d: float # in, from the compression face
controls: str # 'T' tension-controlled, 'C' compression, 'TR' transition
strain_ok: bool = True
note: str = ''
def strand_stress_ksi(eps: float, fpu: float = 270.0,
low_relaxation: bool = True) -> float:
"""PCI 8th Ed. Design Aid 15.3.3 power formula (270 ksi strand)."""
if eps <= 0:
return 0.0
if low_relaxation:
f = eps * (887.0 + 27613.0 / (1.0 + (112.4 * eps) ** 7.36) ** (1.0 / 7.36))
else:
f = eps * (885.0 + 27645.0 / (1.0 + (118.0 * eps) ** 6.0) ** (1.0 / 6.0))
return min(f, fpu)
def _forces(layers: List[Layer], c: float, section_block, beta1: float,
fc: float, ecu: float = 0.003) -> Tuple[float, float, list]:
"""(net force = C - T, moment about compression face of all forces,
per-layer (force, stress)) for a trial c. Compression positive."""
a = beta1 * c
ac, ybar = section_block(a)
cc = 0.85 * fc * ac
net = cc
mom = -cc * ybar # compression above → negative lever from face
per = []
for L in layers:
eps = ecu * (c - L.depth) / c
if L.kind == 'bar':
fs = max(-L.f_limit, min(L.f_limit, L.es * eps))
if fs > 0 and L.depth < c:
# displaced concrete: bar within the neutral-axis depth c
# (the reference tool's convention; the Whitney block is a
# resultant device, not the bar-side boundary)
fs -= 0.85 * fc
force = fs * L.area # + compression, - tension
else:
# eps is the concrete strain at the strand (+ compression);
# the strand's own strain is its effective prestrain minus that
eps_ps = L.fse / L.es - eps
fps = min(strand_stress_ksi(max(eps_ps, 0.0), L.fpu), L.f_limit)
force = -fps * L.area
net += force
mom += -force * L.depth # tension (negative force) below → positive moment
per.append((force, force / L.area if L.area else 0.0))
return net, mom, per
def solve_section(layers: List[Layer], section_block, height: float,
beta1: float, fc: float, phi_t: float, phi_c: float,
prestressed: bool, fy: float = 60.0, es: float = 29000.0,
strain_limit: float = 0.004) -> FlexureResult:
"""Find c by bisection so the section is in equilibrium; report."""
tension = [L for L in layers if L.f_limit > 0]
if not tension:
return FlexureResult(0.0, 0.0, phi_t, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
'NC', True, 'no developed reinforcement')
lo, hi = 1e-4, height
net_lo, _, _ = _forces(layers, lo, section_block, beta1, fc)
net_hi, _, _ = _forces(layers, hi, section_block, beta1, fc)
if net_lo > 0:
# even a hairline block out-pulls the steel: no tension capacity
c = lo
elif net_hi < 0:
c = hi
else:
for _ in range(80):
mid = (lo + hi) / 2.0
net_mid, _, _ = _forces(layers, mid, section_block, beta1, fc)
if net_mid > 0:
hi = mid
else:
lo = mid
if hi - lo < 1e-7:
break
c = (lo + hi) / 2.0
net, mom, per = _forces(layers, c, section_block, beta1, fc)
mn_kin = mom
# extreme tension steel strain
d_t = max(L.depth for L in layers)
eps_t = 0.003 * (d_t - c) / c
eps_ty = 0.002 if prestressed else fy / es
if eps_t >= 0.005:
phi = phi_t
ctl = 'T'
elif eps_t <= eps_ty:
phi = phi_c
ctl = 'C'
else:
phi = phi_c + (phi_t - phi_c) * (eps_t - eps_ty) / (0.005 - eps_ty)
ctl = 'TR'
as_t = sum(L.area * L.dev for (L, (f, s)) in zip(layers, per)
if L.kind == 'bar' and f < 0)
as_c = sum(L.area for (L, (f, s)) in zip(layers, per)
if L.kind == 'bar' and f >= 0)
aps = sum(L.area for L in layers if L.kind == 'strand')
# Reported d: "computed to the center of tension force in bars" (the
# reference tool's own footnote). Each tension layer weighs in with its
# developed AREA times the STRAIN-COMPATIBILITY stress at the solved c,
# capped at yield only -- not at the development stress cap the
# equilibrium solve uses. With the cap, two layers developing together
# both sit at the cap and d collapses to their plain area centroid
# (12.80 flat on RBA5); Eriksson's d ramps 12.80 -> 14.53 -> 19.49 ->
# 20.00 as the top layer's strain dies off approaching c. Probed on both
# fixtures 2026-09-05: As, c and phiMn unchanged bit for bit.
wt = []
for L, (f, s) in zip(layers, per):
if f < 0 and c > 0:
stress = L.es * 0.003 * (L.depth - c) / c
cap = L.fy if L.kind == 'bar' else L.f_limit
stress = max(-cap, min(cap, stress))
wt.append((L.area * L.dev * stress, L.depth))
d = sum(a * y for a, y in wt) / sum(a for a, _ in wt) if wt and sum(a for a, _ in wt) > 0 else 0.0
strain_ok = (eps_t >= strain_limit) if (not prestressed) else True
return FlexureResult(phi * mn_kin / 12.0, mn_kin / 12.0, phi, c, eps_t,
as_t, as_c, aps, d, ctl, strain_ok)
def _bar_layers(beam: Beam, x_in: float, positive: bool,
confined: bool = False) -> List[Layer]:
h = beam.section.height_in
mats = beam.materials
fy = mats.rebar.fy_ksi
out = []
for bx, by, area, db, e0, e1, row in beam.bars_at(x_in):
depth = (h - by) if positive else by
top = (by > h / 2.0)
f_cap = develop.bar_fraction(e0, e1, row, mats, confined,
beam.settings.mild_dev_mult, top=top,
for_capacity=True)
f_rep = develop.bar_fraction(e0, e1, row, mats, confined,
beam.settings.mild_dev_mult, top=top)
out.append(Layer(depth, area, 'bar', f_cap * fy, mats.rebar.es_ksi,
fy, dev=f_rep))
return out
def _strand_layers(beam: Beam, x_in: float, positive: bool,
fse_ksi: float, fps_nominal: float) -> List[Layer]:
h = beam.section.height_in
mats = beam.materials
out = []
for g in beam.strands:
db = _strand_db(g.size)
for sx, sy, area in g.strands(beam.section):
depth = (h - sy) if positive else sy
embed = min(x_in, beam.length_in - x_in)
avail = develop.strand_available_ksi(embed, db, fse_ksi, fps_nominal)
out.append(Layer(depth, area, 'strand', avail, mats.strand.ep_ksi,
fse=min(fse_ksi, avail), fpu=mats.strand.fpu_ksi))
return out
def _strand_db(size: str) -> float:
s = size.replace(' Sp', '').replace('"', '').strip()
return {'3/8': 0.375, '7/16': 0.4375, '1/2': 0.5, '9/16': 0.5625,
'0.6': 0.6}.get(s, 0.5)
def capacity_at(beam: Beam, x_in: float, positive: bool = True,
fse_ksi: float = 0.0, fc_ksi: Optional[float] = None
) -> FlexureResult:
"""φMn at x for positive (sagging) or negative bending."""
mats = beam.materials
fc = fc_ksi if fc_ksi is not None else mats.fc_ksi
beta1 = mats.beta1_at(fc)
sec = beam.section
layers = _bar_layers(beam, x_in, positive)
if beam.strands:
fps_nom = 0.98 * mats.strand.fpu_ksi
layers += _strand_layers(beam, x_in, positive, fse_ksi, fps_nom)
block = sec.compression_area if positive else sec.compression_area_from_bottom
return solve_section(layers, block, sec.height_in, beta1, fc,
beam.settings.phi_tension, beam.settings.phi_compression,
bool(beam.strands), mats.rebar.fy_ksi,
mats.rebar.es_ksi, beam.settings.strain_limit)
def cracking_moment(beam: Beam, x_in: float, positive: bool = True,
p_eff_k: float = 0.0, e_in: float = 0.0,
fc_ksi: Optional[float] = None) -> float:
"""Mcr (kip-ft) on the transformed section: (fr + P/A + P·e/S)·S with
fr = 7.5·λ·√f'c (ACI 19.2.3.1). `e_in` is the strand eccentricity
below the centroid (positive = below, which helps positive bending)."""
from . import section as sec_mod
mats = beam.materials
fc = fc_ksi if fc_ksi is not None else mats.fc_ksi
fr = 7.5 * mats.lam * math.sqrt(fc * 1000.0) / 1000.0 # ksi
p = sec_mod.props_at(beam, x_in)
s = p.sb if positive else p.st
sign = 1.0 if positive else -1.0
f_pre = p_eff_k / p.area + sign * p_eff_k * e_in / s if p_eff_k else 0.0
return (fr + f_pre) * s / 12.0
+138
View File
@@ -0,0 +1,138 @@
# -*- coding: utf-8 -*-
"""Load cases, patterns and combinations.
The reference tool carries nine load-case columns on every load row
(`model.CASES`) and combines them as D / L / R / W / E:
D = self weight + deadnc + topping + deadc
L = live (and sus_live rides inside live for strength; it is its own
case only for the sustained-stress and long-term checks)
R = roof
W = pressure or suction, whichever governs, each sign tried
E = seismic
`case_loads(beam, case)` returns the `statics.Load` list for one case.
`patterns(beam)` returns the live-load patterns: with "Pattern Live Loads"
on, live load is applied to each of the beam's three regions (left
overhang, span, right overhang) in every on/off combination, and every
station reports the worst of them. With it off, live is all-on.
"""
from __future__ import annotations
import itertools
from typing import Dict, List, Sequence, Tuple
from .model import Beam, CASES
from .statics import Load
SELF = 'self'
DEAD_CASES = ('deadnc', 'topping', 'deadc')
def self_weight(beam: Beam) -> List[Load]:
w = beam.self_weight_klf / 12.0 # k/in
return [Load('line', 0.0, beam.length_in, w0=w, w1=w, tag='self weight')]
def case_loads(beam: Beam, case: str) -> List[Load]:
"""Every load in one case, in kips / kip-in / k-in."""
if case == SELF:
return self_weight(beam)
out = []
for p in beam.point_loads:
v = p.cases.get(case, 0.0)
if v:
out.append(Load('point', p.x_in, p=v, tag=p.id))
for l in beam.line_loads:
w0, w1 = l.cases.get(case, (0.0, 0.0))
if w0 or w1:
out.append(Load('line', l.x0_in, l.x1_in, w0=w0 / 12.0,
w1=w1 / 12.0, tag=l.id))
for m in beam.point_moments:
v = m.cases.get(case, 0.0)
if v:
out.append(Load('moment', m.x_in, m=v * 12.0, tag=m.id))
return out
def clip(loads: Sequence[Load], x0: float, x1: float) -> List[Load]:
"""The part of each load that lies within [x0, x1] -- how a pattern
keeps only the live load on one region of the beam."""
out = []
for l in loads:
if l.kind in ('point', 'moment'):
if x0 - 1e-9 <= l.x0 <= x1 + 1e-9:
out.append(l)
else:
a = max(l.x0, x0)
b = min(l.x1, x1)
if b > a + 1e-9:
L = l.x1 - l.x0
wa = l.w0 + (l.w1 - l.w0) * (a - l.x0) / L if L else l.w0
wb = l.w0 + (l.w1 - l.w0) * (b - l.x0) / L if L else l.w1
out.append(Load('line', a, b, w0=wa, w1=wb, tag=l.tag))
return out
def regions(beam: Beam) -> List[Tuple[float, float]]:
"""The three pattern regions, dropping any of zero length."""
xa, xb = beam.supports[0].x_in, beam.supports[1].x_in
regs = [(0.0, xa), (xa, xb), (xb, beam.length_in)]
return [(a, b) for a, b in regs if b - a > 1e-6]
def patterns(beam: Beam) -> List[Tuple[Tuple[float, float], ...]]:
"""Every non-empty on/off pattern of the regions (all-on first), or just
all-on when patterning is off. The all-on pattern is always included so
a result can never be less than the un-patterned answer."""
regs = regions(beam)
if not beam.pattern_live or len(regs) == 1:
return [tuple(regs)]
out = [tuple(regs)]
for n in range(1, len(regs)):
for combo in itertools.combinations(regs, n):
out.append(tuple(combo))
return out
def live_loads_for(beam: Beam, pattern: Sequence[Tuple[float, float]],
case: str = 'live') -> List[Load]:
full = case_loads(beam, case)
out = []
for a, b in pattern:
out.extend(clip(full, a, b))
return out
def dead_loads(beam: Beam, include_self: bool = True) -> List[Load]:
out = []
if include_self:
out.extend(self_weight(beam))
for c in DEAD_CASES:
out.extend(case_loads(beam, c))
return out
def scaled(loads: Sequence[Load], f: float) -> List[Load]:
if f == 1.0:
return list(loads)
if f == 0.0:
return []
return [Load(l.kind, l.x0, l.x1, l.p * f, l.w0 * f, l.w1 * f, l.m * f, l.tag)
for l in loads]
def has_case(beam: Beam, case: str) -> bool:
return bool(case_loads(beam, case))
def wind_variants(beam: Beam) -> List[Tuple[str, List[Load]]]:
"""W as the reference tool takes it: pressure and suction are separate
cases and either may govern; each is tried with both signs."""
out = []
for c in ('pressure', 'suction'):
loads = case_loads(beam, c)
if loads:
out.append((c, loads))
out.append((c + ' (reversed)', scaled(loads, -1.0)))
return out
+771
View File
@@ -0,0 +1,771 @@
# -*- coding: utf-8 -*-
"""The job file (`beamcalc/1`) and the dataclasses it becomes.
The job file is the whole input: a beam re-run from it on a machine with no
Revit and no licence gives the same numbers. Its blocks mirror the reference
tool's own input pages one for one (Design Criteria, Concrete Extents,
Structural Model, Stripping & Handling, Reinforcement) so a package printed
by that tool can be transcribed into a job by hand -- that is exactly how the
parity fixtures in `beamcalc/fixtures/` were made.
Lengths in the job file are INCHES unless the key says otherwise
(`length_ft` is the one exception, because everyone says a beam length in
feet). Forces are kips, moments kip-ft, stresses ksi, line loads klf, area
loads psf. `x` positions are measured from the START of the beam unless a
`from: "end"` accompanies them.
"""
from __future__ import annotations
import json
import math
from dataclasses import dataclass, field
from typing import Dict, List, Optional, Sequence
SCHEMA = 'beamcalc/1'
# Load-case columns, in the reference tool's own order and naming. Every load
# row carries one value per case (a line load carries a start and end value
# per case, so it can taper).
CASES = ('deadnc', 'topping', 'deadc', 'live', 'sus_live', 'roof',
'pressure', 'suction', 'seismic')
# Bar areas (in²) and diameters (in) by size, ACI 318 Appendix A.
BAR_AREA = {'#3': 0.11, '#4': 0.20, '#5': 0.31, '#6': 0.44, '#7': 0.60,
'#8': 0.79, '#9': 1.00, '#10': 1.27, '#11': 1.56, '#14': 2.25,
'#18': 4.00}
BAR_DIA = {'#3': 0.375, '#4': 0.500, '#5': 0.625, '#6': 0.750, '#7': 0.875,
'#8': 1.000, '#9': 1.128, '#10': 1.270, '#11': 1.410,
'#14': 1.693, '#18': 2.257}
# Strand areas (in²) by nominal size (7-wire, ASTM A416).
STRAND_AREA = {'3/8': 0.085, '7/16': 0.115, '1/2': 0.153, '1/2 Sp': 0.167,
'9/16': 0.192, '0.6': 0.217}
class JobError(ValueError):
"""The job file cannot be read into a beam. The message names the key."""
def _num(d, key, default=None, required=False):
v = d.get(key, default)
if v is None or v == '':
if required:
raise JobError('missing required value: {0}'.format(key))
return default
try:
return float(v)
except (TypeError, ValueError):
raise JobError('{0} is not a number: {1!r}'.format(key, v))
@dataclass
class Rebar:
fy_ksi: float = 60.0
es_ksi: float = 29000.0
@dataclass
class StrandMaterial:
fpu_ksi: float = 270.0
ep_ksi: float = 28500.0
low_relaxation: bool = True
@dataclass
class Materials:
fc_ksi: float
fci_ksi: float
wc_pcf: float = 150.0
ec_ksi: Optional[float] = None # None -> ACI 19.2.2.1(a)
eci_ksi: Optional[float] = None
normalweight: bool = True
lam: float = 1.0
rebar: Rebar = field(default_factory=Rebar)
trans: Rebar = field(default_factory=Rebar)
strand: StrandMaterial = field(default_factory=StrandMaterial)
@staticmethod
def modulus_ksi(wc_pcf: float, fc_ksi: float) -> float:
"""ACI 318-19 19.2.2.1(a): Ec = wc^1.5 * 33 * sqrt(f'c) (psi)."""
return wc_pcf ** 1.5 * 33.0 * math.sqrt(fc_ksi * 1000.0) / 1000.0
@property
def Ec(self) -> float:
return self.ec_ksi or self.modulus_ksi(self.wc_pcf, self.fc_ksi)
@property
def Eci(self) -> float:
return self.eci_ksi or self.modulus_ksi(self.wc_pcf, self.fci_ksi)
@property
def beta1(self) -> float:
"""ACI 318-19 Table 22.2.2.4.3."""
fc = self.fc_ksi * 1000.0
if fc <= 4000.0:
return 0.85
if fc >= 8000.0:
return 0.65
return 0.85 - 0.05 * (fc - 4000.0) / 1000.0
def beta1_at(self, fc_ksi: float) -> float:
fc = fc_ksi * 1000.0
if fc <= 4000.0:
return 0.85
if fc >= 8000.0:
return 0.65
return 0.85 - 0.05 * (fc - 4000.0) / 1000.0
@dataclass
class Combo:
D: float
L: float
R: float
W: float
E: float
include: bool = True
label: str = ''
def name(self) -> str:
if self.label:
return self.label
parts = []
for f, s in ((self.D, 'D'), (self.L, 'L'), (self.R, 'R'),
(self.W, 'W'), (self.E, 'E')):
if f:
parts.append('{0:g}{1}'.format(f, s))
return ' + '.join(parts) or '0'
@dataclass
class Settings:
total_defl_ratio: float = 240.0 # L/240
live_defl_ratio: float = 360.0
ec_for_sw_defl: str = 'Eci' # 'Ec' | 'Eci'
sw_mult: Dict[str, float] = field(default_factory=lambda: {
'erection': 1.0, 'construction': 1.3, 'end': 2.0})
addl_dead_mult: Dict[str, float] = field(default_factory=lambda: {
'construction': 0.0, 'end': 2.0})
composite_shear: str = 'varying_strength'
dv_method: str = 'cg_tension'
min_shear_method: str = 'beam_9.6.3'
flexure_method: str = 'whitney'
ignored_distance_in: float = 0.0
strain_limit: float = 0.004
transfer_comp_ends: float = 0.7 # x f'ci
transfer_comp_other: float = 0.7
transfer_tens_ends_sqrt: float = 6.0 # x sqrt(f'ci) psi
transfer_tens_other_sqrt: float = 6.0
service_comp_sustained: float = 0.45 # x f'c
service_comp_total: float = 0.60
phi_tension: float = 0.9
phi_shear: float = 0.75
phi_compression: float = 0.65
ktr: float = 0.0
mild_dev_mult: float = 1.0
limit_shear_near_supports: bool = True
transform_reinforcement: bool = True
@dataclass
class Section:
"""A section built of axis-aligned rectangles. `kind` is what the user
chose; `rects` is what the engine uses -- (x0, y0, w, h) in inches with
y measured UP from the bottom of the section and x from its left edge."""
kind: str # 'rectangular' | 'L' | 'inverted_t'
width_in: float # overall bounding width
height_in: float # overall height
rects: List[tuple]
stem_width_in: Optional[float] = None
stem_height_in: Optional[float] = None
flange_width_in: Optional[float] = None
flange_thickness_in: Optional[float] = None
side: str = 'left'
@classmethod
def rectangular(cls, b: float, h: float) -> 'Section':
return cls('rectangular', b, h, [(0.0, 0.0, b, h)])
@classmethod
def l_beam(cls, stem_w: float, stem_h: float, flange_w: float,
flange_t: float, side: str = 'left') -> 'Section':
"""Spandrel / L-beam the way the reference tool defines it: a stem of
`stem_w` x (full height) with a ledge `flange_w` wide and `flange_t`
tall at the bottom on one side. `stem_h` is the height ABOVE the
ledge, so overall height is stem_h + flange_t."""
h = stem_h + flange_t
b = stem_w + flange_w
if side == 'left':
rects = [(flange_w, 0.0, stem_w, h), (0.0, 0.0, flange_w, flange_t)]
else:
rects = [(0.0, 0.0, stem_w, h), (stem_w, 0.0, flange_w, flange_t)]
s = cls('L', b, h, rects, stem_w, stem_h, flange_w, flange_t, side)
return s
@classmethod
def inverted_t(cls, stem_w: float, stem_h: float, flange_w: float,
flange_t: float) -> 'Section':
"""Inverted tee: a ledge `flange_w` wide each side of the stem."""
h = stem_h + flange_t
b = stem_w + 2 * flange_w
rects = [(flange_w, 0.0, stem_w, h), (0.0, 0.0, b, flange_t)]
# the bottom rect overlaps the stem; split it to avoid double count
rects = [(flange_w, 0.0, stem_w, h),
(0.0, 0.0, flange_w, flange_t),
(flange_w + stem_w, 0.0, flange_w, flange_t)]
return cls('inverted_t', b, h, rects, stem_w, stem_h, flange_w,
flange_t, 'both')
# -- gross properties -----------------------------------------------
@property
def area(self) -> float:
return sum(w * h for _, _, w, h in self.rects)
@property
def cy(self) -> float:
"""Centroid height above the bottom."""
return sum(w * h * (y + h / 2.0) for _, y, w, h in self.rects) / self.area
@property
def cx(self) -> float:
return sum(w * h * (x + w / 2.0) for x, _, w, h in self.rects) / self.area
@property
def ixx(self) -> float:
c = self.cy
return sum(w * h ** 3 / 12.0 + w * h * (y + h / 2.0 - c) ** 2
for _, y, w, h in self.rects)
@property
def iyy(self) -> float:
c = self.cx
return sum(h * w ** 3 / 12.0 + w * h * (x + w / 2.0 - c) ** 2
for x, _, w, h in self.rects)
@property
def st(self) -> float:
return self.ixx / (self.height_in - self.cy)
@property
def sb(self) -> float:
return self.ixx / self.cy
def width_at(self, y: float) -> float:
"""Total concrete width at height y above the bottom."""
return sum(w for _, y0, w, h in self.rects if y0 <= y < y0 + h)
def web_width(self) -> float:
"""bw for shear: the stem width (the narrowest full-height part)."""
if self.stem_width_in:
return self.stem_width_in
return self.width_in
def compression_area(self, depth_from_top: float) -> tuple:
"""(area, centroid depth from top) of the concrete within
`depth_from_top` of the top face -- the Whitney block for positive
bending. Handles a flange at the bottom by just not being there."""
top = self.height_in
a_tot = 0.0
m_tot = 0.0
for _, y0, w, h in self.rects:
y1 = y0 + h
lo = max(y0, top - depth_from_top)
hi = min(y1, top)
if hi > lo:
a = w * (hi - lo)
a_tot += a
m_tot += a * (top - (lo + hi) / 2.0)
return a_tot, (m_tot / a_tot if a_tot else 0.0)
def compression_area_from_bottom(self, depth_from_bottom: float) -> tuple:
"""Same, for negative bending (block at the bottom)."""
a_tot = 0.0
m_tot = 0.0
for _, y0, w, h in self.rects:
y1 = y0 + h
lo = max(y0, 0.0)
hi = min(y1, depth_from_bottom)
if hi > lo:
a = w * (hi - lo)
a_tot += a
m_tot += a * ((lo + hi) / 2.0)
return a_tot, (m_tot / a_tot if a_tot else 0.0)
@dataclass
class Support:
x_in: float # centreline of bearing, from the start
length_in: float # along the beam
width_in: float # across the beam
kind: str = 'pin'
@property
def face_start_in(self) -> float:
return self.x_in - self.length_in / 2.0
@property
def face_end_in(self) -> float:
return self.x_in + self.length_in / 2.0
@dataclass
class PointLoad:
id: str
x_in: float
ecc_in: float
cases: Dict[str, float]
description: str = ''
@dataclass
class LineLoad:
id: str
x0_in: float
x1_in: float
ecc0_in: float
ecc1_in: float
cases: Dict[str, tuple] # case -> (w0, w1) klf
description: str = ''
@dataclass
class PointMoment:
id: str
x_in: float
cases: Dict[str, float] # kip-ft
@dataclass
class BarRow:
"""One row of the reference tool's Rebar table: `cols` bars of `size`
in `rows` rows, the group's edges `left`/`right` from the section's
left/right faces, `top`/`bottom` the distance from the top/bottom faces
to the group's outermost bar CENTRES, running from `start_in` after the
beam start to `end_in` before the beam end."""
size: str
rows: int
cols: int
left_in: float
right_in: float
top_in: float
bottom_in: float
start_in: float = 0.0
end_in: float = 0.0
start_cond: str = 'straight'
end_cond: str = 'straight'
area_in2: Optional[float] = None # override; else BAR_AREA[size]
@property
def bar_area(self) -> float:
return self.area_in2 if self.area_in2 else BAR_AREA[self.size]
@property
def db(self) -> float:
return BAR_DIA.get(self.size, math.sqrt(4 * self.bar_area / math.pi))
def bars(self, section: Section) -> List[tuple]:
"""[(x, y, area, db), ...] every bar centre, y up from the bottom."""
out = []
h = section.height_in
b = section.width_in
ys = ([h - self.top_in] if self.rows == 1 else
[h - self.top_in - i * (h - self.top_in - self.bottom_in)
/ (self.rows - 1) for i in range(self.rows)])
xs = ([self.left_in] if self.cols == 1 else
[self.left_in + j * (b - self.left_in - self.right_in)
/ (self.cols - 1) for j in range(self.cols)])
for y in ys:
for x in xs:
out.append((x, y, self.bar_area, self.db))
return out
@dataclass
class StrandGroup:
"""One row of the reference tool's Prestress table."""
ids: str
size: str # '1/2', '1/2 Sp', '0.6', ...
pull_pct: float # % of fpu jacked
rows: int
cols: int
left_in: float
right_in: float
first_row_in: float # bottom row height above the bottom face
spacing_in: float # row spacing
force_k: Optional[float] = None # printed total jacking force, if known
@property
def strand_area(self) -> float:
return STRAND_AREA[self.size]
@property
def count(self) -> int:
return self.rows * self.cols
def strands(self, section: Section) -> List[tuple]:
"""[(x, y, area), ...] every strand centre, y up from the bottom."""
out = []
b = section.width_in
xs = ([self.left_in] if self.cols == 1 else
[self.left_in + j * (b - self.left_in - self.right_in)
/ (self.cols - 1) for j in range(self.cols)])
for i in range(self.rows):
y = self.first_row_in + i * self.spacing_in
for x in xs:
out.append((x, y, self.strand_area))
return out
@dataclass
class TransverseZone:
size: str
legs: int
first_in: float # first bar from the start (or end)
from_end: bool
spacing_in: float
length_in: float
qty: int
symmetric: bool = True
use: str = 'Av'
area_in2: Optional[float] = None # override total area per set
@property
def area(self) -> float:
if self.area_in2:
return self.area_in2
return self.legs * BAR_AREA[self.size]
@dataclass
class Stripping:
multiplier: float = 1.4
face_down: str = 'bottom'
strength: str = 'initial' # 'initial' | 'final'
limits: str = 'transfer'
lift_from_start_in: float = 0.0
lift_from_end_in: float = 0.0
pick_points: int = 2
@dataclass
class Beam:
mark: str
length_in: float
section: Section
materials: Materials
supports: List[Support]
settings: Settings
strength_combos: List[Combo]
service_combos: List[Combo]
point_loads: List[PointLoad] = field(default_factory=list)
line_loads: List[LineLoad] = field(default_factory=list)
point_moments: List[PointMoment] = field(default_factory=list)
bars: List[BarRow] = field(default_factory=list)
strands: List[StrandGroup] = field(default_factory=list)
transverse: List[TransverseZone] = field(default_factory=list)
stripping: Optional[Stripping] = None
self_weight_mult: float = 1.0
pattern_live: bool = True
bearing_analysis: bool = True
bearing_mu: float = 1.4
bearing_nu_ratio: float = 0.2
description: str = ''
# -- derived ---------------------------------------------------------
@property
def length_ft(self) -> float:
return self.length_in / 12.0
@property
def self_weight_klf(self) -> float:
return (self.section.area / 144.0) * self.materials.wc_pcf / 1000.0 \
* self.self_weight_mult
@property
def weight_k(self) -> float:
return self.self_weight_klf * self.length_ft
@property
def prestressed(self) -> bool:
return bool(self.strands)
@property
def span_in(self) -> float:
a, b = self.supports[0].x_in, self.supports[1].x_in
return b - a
def bars_at(self, x_in: float) -> List[tuple]:
"""Every bar present at x: [(x, y, area, db, embed_start, embed_end,
row), ...] with the embedment lengths from x to each bar end."""
out = []
for row in self.bars:
x0 = row.start_in
x1 = self.length_in - row.end_in
if x0 - 1e-9 <= x_in <= x1 + 1e-9:
for bx, by, a, db in row.bars(self.section):
out.append((bx, by, a, db, x_in - x0, x1 - x_in, row))
return out
# ---------------------------------------------------------------- parsing
def _x_from(d, key_x, key_from, length_in, default_from='start'):
x = _num(d, key_x, 0.0)
frm = (d.get(key_from) or default_from).lower()
if frm == 'end':
return length_in - x
return x
def _combos(rows, default_rows) -> List[Combo]:
rows = rows if rows is not None else default_rows
out = []
for r in rows:
out.append(Combo(D=float(r.get('D', 0)), L=float(r.get('L', 0)),
R=float(r.get('R', 0)), W=float(r.get('W', 0)),
E=float(r.get('E', 0)),
include=bool(r.get('include', True)),
label=r.get('label', '')))
return out
# ACI 318-19 5.3.1 as the reference tool prints them (Roof as its own case).
DEFAULT_STRENGTH = [
dict(D=1.4, L=0, R=0, W=0, E=0, include=True),
dict(D=1.2, L=1.6, R=0.5, W=0, E=0, include=True),
dict(D=1.2, L=1.0, R=1.6, W=0, E=0, include=True),
dict(D=1.2, L=0, R=1.6, W=0.5, E=0, include=False),
dict(D=1.2, L=1.0, R=0.5, W=1.0, E=0, include=False),
dict(D=1.2, L=1.0, R=0.2, W=0, E=1.0, include=False),
dict(D=0.9, L=0, R=0, W=1.0, E=0, include=False),
dict(D=0.9, L=0, R=0, W=0, E=1.0, include=False),
]
DEFAULT_SERVICE = [
dict(D=1, L=0, R=0, W=0, E=0, include=True),
dict(D=1, L=1, R=0, W=0, E=0, include=True),
dict(D=1, L=0, R=1, W=0, E=0, include=True),
dict(D=1, L=0.75, R=0.75, W=0, E=0, include=True),
dict(D=1, L=0, R=0, W=0.6, E=0, include=False),
dict(D=1, L=0, R=0, W=0, E=0.7, include=False),
dict(D=1, L=0.75, R=0.75, W=0.45, E=0, include=False),
dict(D=1, L=0.75, R=0.75, W=0, E=0.525, include=False),
dict(D=0.6, L=0, R=0, W=0.6, E=0, include=False),
dict(D=0.6, L=0, R=0, W=0, E=0.7, include=False),
]
def _section(d: dict) -> Section:
kind = (d.get('type') or 'rectangular').lower().replace(' ', '_')
if kind in ('rectangular', 'rect'):
return Section.rectangular(_num(d, 'width_in', required=True),
_num(d, 'height_in', required=True))
if kind in ('l', 'l_beam', 'spandrel', 'spandrel_/_l-beam'):
return Section.l_beam(_num(d, 'stem_width_in', required=True),
_num(d, 'stem_height_in', required=True),
_num(d, 'flange_width_in', required=True),
_num(d, 'flange_thickness_in', required=True),
(d.get('side') or 'left').lower())
if kind in ('inverted_t', 'it', 'inverted_tee'):
return Section.inverted_t(_num(d, 'stem_width_in', required=True),
_num(d, 'stem_height_in', required=True),
_num(d, 'flange_width_in', required=True),
_num(d, 'flange_thickness_in', required=True))
raise JobError('unknown section type: {0!r}'.format(kind))
def _materials(d: dict) -> Materials:
rb = d.get('rebar') or {}
tr = d.get('trans') or rb
st = d.get('strand') or {}
return Materials(
fc_ksi=_num(d, 'fc_ksi', required=True),
fci_ksi=_num(d, 'fci_ksi', _num(d, 'fc_ksi')),
wc_pcf=_num(d, 'wc_pcf', 150.0),
ec_ksi=_num(d, 'ec_ksi'),
eci_ksi=_num(d, 'eci_ksi'),
normalweight=(d.get('type') or 'normalweight').lower().startswith('normal'),
lam=_num(d, 'lambda', 1.0),
rebar=Rebar(_num(rb, 'fy_ksi', 60.0), _num(rb, 'es_ksi', 29000.0)),
trans=Rebar(_num(tr, 'fy_ksi', 60.0), _num(tr, 'es_ksi', 29000.0)),
strand=StrandMaterial(_num(st, 'fpu_ksi', 270.0),
_num(st, 'ep_ksi', 28500.0),
bool(st.get('low_relaxation', True))))
def _settings(d: dict) -> Settings:
s = Settings()
if not d:
return s
defl = d.get('deflection') or {}
for key, attr in (('total', 'total_defl_ratio'), ('live', 'live_defl_ratio')):
v = defl.get(key)
if v is not None:
if isinstance(v, str) and '/' in v:
v = float(v.split('/')[1])
setattr(s, attr, float(v))
if defl.get('ec_for_sw'):
s.ec_for_sw_defl = defl['ec_for_sw']
mult = defl.get('multipliers') or {}
if mult.get('self_weight'):
s.sw_mult.update({k: float(v) for k, v in mult['self_weight'].items()
if v is not None})
if mult.get('addl_dead'):
s.addl_dead_mult.update({k: float(v) for k, v in mult['addl_dead'].items()
if v is not None})
phi = d.get('phi') or {}
s.phi_tension = _num(phi, 'tension', s.phi_tension)
s.phi_shear = _num(phi, 'shear', s.phi_shear)
s.phi_compression = _num(phi, 'compression', s.phi_compression)
tl = d.get('transfer_limits') or {}
s.transfer_comp_ends = _num(tl, 'comp_ends', s.transfer_comp_ends)
s.transfer_comp_other = _num(tl, 'comp_other', s.transfer_comp_other)
s.transfer_tens_ends_sqrt = _num(tl, 'tens_ends_sqrt', s.transfer_tens_ends_sqrt)
s.transfer_tens_other_sqrt = _num(tl, 'tens_other_sqrt', s.transfer_tens_other_sqrt)
sl = d.get('service_limits') or {}
s.service_comp_sustained = _num(sl, 'comp_sustained', s.service_comp_sustained)
s.service_comp_total = _num(sl, 'comp_total', s.service_comp_total)
fx = d.get('flexure') or {}
s.ignored_distance_in = _num(fx, 'ignored_distance_in', 0.0)
s.strain_limit = _num(fx, 'strain_limit', 0.004)
sh = d.get('shear') or {}
s.dv_method = sh.get('dv_method', s.dv_method)
s.min_shear_method = sh.get('min_method', s.min_shear_method)
dv = d.get('development') or {}
s.ktr = _num(dv, 'ktr', 0.0)
s.mild_dev_mult = _num(dv, 'mild_mult', 1.0)
s.transform_reinforcement = bool(
(d.get('section_properties') or {}).get('transform', True))
return s
def _cases(d: dict, pair=False) -> dict:
out = {}
for c in CASES:
v = d.get(c, 0)
if pair:
if isinstance(v, (list, tuple)):
out[c] = (float(v[0]), float(v[-1]))
else:
out[c] = (float(v or 0), float(v or 0))
else:
out[c] = float(v or 0)
return out
def from_job(job: dict) -> Beam:
"""Build a `Beam` from a parsed `beamcalc/1` job dict. Raises JobError
naming the first thing wrong; never returns a half-read beam."""
if job.get('schema') != SCHEMA:
raise JobError('schema is {0!r}, this engine reads {1!r}'
.format(job.get('schema'), SCHEMA))
ext = job.get('extents') or {}
length_in = _num(ext, 'length_in')
if length_in is None:
length_ft = _num(ext, 'length_ft', required=True)
length_in = length_ft * 12.0
section = _section(ext.get('section') or {})
mats = _materials(job.get('materials') or {})
settings = _settings(job.get('settings') or {})
sup = job.get('supports') or {}
s0 = sup.get('start') or {}
s1 = sup.get('end') or {}
supports = [
Support(_num(s0, 'from_start_in', 0.0), _num(s0, 'length_in', 0.0),
_num(s0, 'width_in', section.width_in), s0.get('type', 'pin')),
Support(length_in - _num(s1, 'from_end_in', 0.0),
_num(s1, 'length_in', 0.0),
_num(s1, 'width_in', section.width_in), s1.get('type', 'pin')),
]
if supports[1].x_in <= supports[0].x_in:
raise JobError('supports cross: start at {0}", end at {1}"'
.format(supports[0].x_in, supports[1].x_in))
combos = job.get('combos') or {}
strength = _combos(combos.get('strength'), DEFAULT_STRENGTH)
service = _combos(combos.get('service'), DEFAULT_SERVICE)
loads = job.get('loads') or {}
pls = []
for i, p in enumerate(loads.get('point') or []):
pls.append(PointLoad(p.get('id') or 'P{0}'.format(i + 1),
_x_from(p, 'x_in', 'from', length_in),
_num(p, 'ecc_in', 0.0), _cases(p),
p.get('description', '')))
lls = []
for i, l in enumerate(loads.get('line') or []):
x0 = _x_from(l, 'x0_in', 'x0_from', length_in, 'start')
x1 = _x_from(l, 'x1_in', 'x1_from', length_in, 'end')
if x1 < x0:
x0, x1 = x1, x0
lls.append(LineLoad(l.get('id') or 'L{0}'.format(i + 1), x0, x1,
_num(l, 'ecc0_in', _num(l, 'ecc_in', 0.0)),
_num(l, 'ecc1_in', _num(l, 'ecc_in', 0.0)),
_cases(l, pair=True), l.get('description', '')))
pms = []
for i, m in enumerate(loads.get('moments') or []):
pms.append(PointMoment(m.get('id') or 'M{0}'.format(i + 1),
_x_from(m, 'x_in', 'from', length_in),
_cases(m)))
reinf = job.get('reinforcement') or {}
bars = []
for r in reinf.get('bars') or []:
bars.append(BarRow(r['size'], int(r.get('rows', 1)), int(r.get('cols', 1)),
_num(r, 'left_in', 2.0), _num(r, 'right_in', 2.0),
_num(r, 'top_in', 2.0), _num(r, 'bottom_in', 2.0),
_num(r, 'start_in', 0.0), _num(r, 'end_in', 0.0),
r.get('start_cond', 'straight'),
r.get('end_cond', 'straight'),
_num(r, 'area_in2')))
strands = []
for g in reinf.get('strands') or []:
strands.append(StrandGroup(g.get('ids', ''), g['size'],
_num(g, 'pull_pct', 70.0),
int(g.get('rows', 1)), int(g.get('cols', 1)),
_num(g, 'left_in', 2.0), _num(g, 'right_in', 2.0),
_num(g, 'first_row_in', 2.0),
_num(g, 'spacing_in', 2.0),
_num(g, 'force_k')))
trans = []
for z in reinf.get('transverse') or []:
trans.append(TransverseZone(z['size'], int(z.get('legs', 2)),
_num(z, 'first_in', 0.0),
(z.get('from') or 'start').lower() == 'end',
_num(z, 'spacing_in', required=True),
_num(z, 'length_in', 0.0),
int(z.get('qty', 0)),
bool(z.get('symmetric', True)),
z.get('use', 'Av'), _num(z, 'area_in2')))
hd = job.get('handling') or {}
st = hd.get('stripping')
stripping = None
if st:
stripping = Stripping(_num(st, 'multiplier', 1.4),
st.get('face_down', 'bottom'),
st.get('strength', 'initial'),
st.get('limits', 'transfer'),
_num(st, 'lift_from_start_in', 0.0),
_num(st, 'lift_from_end_in', 0.0),
int(st.get('pick_points', 2)))
br = job.get('bearing') or {}
return Beam(
mark=str(job.get('mark') or ''),
length_in=length_in, section=section, materials=mats,
supports=supports, settings=settings,
strength_combos=strength, service_combos=service,
point_loads=pls, line_loads=lls, point_moments=pms,
bars=bars, strands=strands, transverse=trans, stripping=stripping,
self_weight_mult=_num(loads, 'self_weight_mult', 1.0),
pattern_live=bool(loads.get('pattern_live', True)),
bearing_analysis=bool(br.get('analyse', True)),
bearing_mu=_num(br, 'mu', 1.4),
bearing_nu_ratio=_num(br, 'nu_ratio', 0.2),
description=str(job.get('description') or ''))
def load_job(path: str) -> Beam:
with open(path, 'r', encoding='utf-8') as fh:
return from_job(json.load(fh))
+383
View File
@@ -0,0 +1,383 @@
# -*- coding: utf-8 -*-
"""The results HTML (the Beam Analysis 后端 page).
Same shape as the Panel Analysis results page: verdict banner, one card per
stage, then collapsible sections that echo the inputs and print every
per-station table the reference tool prints, with SVG diagrams drawn by
`panellift.report.Plot` so the two reports read alike. Lengths print as
feet-inches through `panellift.units`.
"""
from __future__ import annotations
import datetime
import html
import os
import sys
from typing import List, Optional
_HERE = os.path.dirname(os.path.abspath(__file__))
_PANELLIFT = os.path.abspath(os.path.join(_HERE, '..', '..', '..', 'panellift', 'src'))
if _PANELLIFT not in sys.path:
sys.path.insert(0, _PANELLIFT)
from panellift.report import Plot, CSS as _PL_CSS, hover_script # noqa: E402
from panellift.units import fmt_ft_in, fmt_in # noqa: E402
from .solve import Analysis
from .model import Beam, CASES
def esc(s) -> str:
return html.escape(str(s), quote=True)
def n(v, d=1) -> str:
if v is None:
return '-'
if isinstance(v, float) and (v != v or v in (float('inf'), float('-inf'))):
return 'NC'
return ('{0:,.%df}' % d).format(v)
def ratio(v) -> str:
if v is None or v != v:
return '-'
if v == float('inf'):
return 'inf'
cls = 'fail' if v > 1.0 + 1e-9 else 'ok'
return '<span class="{0}">{1}</span>'.format(cls, n(v, 2))
CSS = _PL_CSS + """
.cards{display:flex;flex-wrap:wrap;gap:10px;margin:12px 0}
.card{border:1px solid #d8dde3;border-left:5px solid #2a7;border-radius:6px;padding:8px 12px;min-width:180px;background:#fff}
.card.bad{border-left-color:#c33}
.card .t{font-weight:600;font-size:13px}.card .u{font-size:22px;font-weight:700;margin:2px 0}
.card .s{font-size:12px;color:#555}
.verdict{padding:8px 12px;border-radius:6px;margin:10px 0;font-weight:600}
.verdict.ok{background:#e6f5ea;color:#1a6b34}.verdict.bad{background:#fbe6e6;color:#8a1c1c}
span.fail{color:#b00;font-weight:600}span.ok{color:#1a6b34}
table.st{font-size:12px}table.st td,table.st th{padding:2px 6px;white-space:nowrap}
.note{font-size:12px;color:#555;margin:4px 0 8px}
details{margin:8px 0}summary{cursor:pointer;font-weight:600;padding:6px 8px;background:#f2f4f6;border-radius:4px}
"""
# ---------------------------------------------------------------- pieces
def _sec(title: str, body: str, open_: bool = False, tag: str = '') -> str:
return ('<details{o}><summary>{t}{g}</summary><div class="body">{b}</div>'
'</details>').format(o=' open' if open_ else '', t=esc(title),
g=(' <span class="tag">' + tag + '</span>') if tag else '',
b=body)
def _table(headers: List[str], rows: List[List[str]], cls='st') -> str:
h = ''.join('<th>{0}</th>'.format(x) for x in headers)
b = ''.join('<tr>' + ''.join('<td>{0}</td>'.format(c) for c in r) + '</tr>'
for r in rows)
return '<div class="scroll"><table class="{0}"><thead><tr>{1}</tr></thead>' \
'<tbody>{2}</tbody></table></div>'.format(cls, h, b)
def _kv(pairs) -> str:
return '<table class="kv">' + ''.join(
'<tr><th>{0}</th><td>{1}</td></tr>'.format(esc(k), v) for k, v in pairs
) + '</table>'
def _cards(a: Analysis) -> str:
out = []
for st in a.stages:
bad = not st.ok
sub = ''
if st.name == 'In-Place Strength':
f = st.data.get('flexure')
s = st.data.get('shear')
bits = []
if f:
bits.append('flexure {0} at {1}'.format(n(f.util, 2), fmt_ft_in(f.x_ft)))
if s:
bits.append('shear {0} at {1}'.format(n(s.r.util, 2), fmt_ft_in(s.x_ft)))
sub = '; '.join(bits)
elif st.name == 'In-Place Service':
d = a.deflection
sub = 'total {0}" / {1}", live {2}" / {3}"'.format(
n(d['total_down'], 3), n(d['total_limit'], 3),
n(d['live_down'], 3), n(d['live_limit'], 3))
elif st.name == 'Stripping':
sub = '; '.join(st.lines[1:])
out.append('<div class="card{b}"><div class="t">{t} <span class="{c}">{v}</span></div>'
'<div class="u">{u}</div><div class="s">{s}</div></div>'.format(
b=' bad' if bad else '', t=esc(st.name),
c='fail' if bad else 'ok', v='FAIL' if bad else 'PASS',
u='{0:.0%} utilised'.format(st.util) if st.util < 10 else 'NC',
s=esc(sub)))
return '<div class="cards">' + ''.join(out) + '</div>'
def _inputs(b: Beam) -> str:
m = b.materials
sec = b.section
mat = _kv([
("f'c / f'ci", '{0} / {1} ksi'.format(n(m.fc_ksi, 2), n(m.fci_ksi, 2))),
('Ec / Eci', '{0} / {1} ksi'.format(n(m.Ec, 0), n(m.Eci, 0))),
('unit weight', '{0} pcf, {1}'.format(n(m.wc_pcf, 0),
'normalweight' if m.normalweight else 'lightweight')),
('rebar fy / Es', '{0} / {1} ksi'.format(n(m.rebar.fy_ksi, 0), n(m.rebar.es_ksi, 0))),
('strand fpu / Ep', '{0} / {1} ksi'.format(n(m.strand.fpu_ksi, 0), n(m.strand.ep_ksi, 0))
if b.strands else '-'),
])
geo = _kv([
('length', fmt_ft_in(b.length_ft)),
('section', '{0}, {1} wide x {2} high'.format(sec.kind, fmt_in(sec.width_in),
fmt_in(sec.height_in))),
('area / Ixx', '{0} in² / {1} in⁴'.format(n(sec.area, 1), n(sec.ixx, 0))),
('St / Sb (gross)', '{0} / {1} in³'.format(n(sec.st, 0), n(sec.sb, 0))),
('self weight', '{0} klf, {1} kips total'.format(n(b.self_weight_klf, 3),
n(b.weight_k, 2))),
('bearings', 'start CL at {0} ({1} long), end CL at {2} from the end ({3} long)'
.format(fmt_in(b.supports[0].x_in), fmt_in(b.supports[0].length_in),
fmt_in(b.length_in - b.supports[1].x_in), fmt_in(b.supports[1].length_in))),
('span (CL to CL)', fmt_ft_in(b.span_in / 12.0)),
])
load_rows = []
for p in b.point_loads:
load_rows.append([esc(p.id), 'point', fmt_ft_in(p.x_in / 12.0), '-'] +
[n(p.cases.get(c, 0.0), 2) for c in CASES])
for l in b.line_loads:
load_rows.append([esc(l.id), 'line', fmt_ft_in(l.x0_in / 12.0),
fmt_ft_in(l.x1_in / 12.0)] +
['{0}'.format(n(l.cases.get(c, (0, 0))[0], 2)) for c in CASES])
loads = _table(['id', 'kind', 'from', 'to'] + list(CASES), load_rows) if load_rows \
else '<p class="note">no applied loads</p>'
combos = _table(['#', 'D', 'L', 'R', 'W', 'E', 'include'],
[[str(i + 1), n(c.D, 2), n(c.L, 2), n(c.R, 2), n(c.W, 2), n(c.E, 2),
'yes' if c.include else 'no'] for i, c in enumerate(b.strength_combos)])
bars = _table(['size', 'rows x cols', 'from top', 'from bottom', 'left', 'right',
'starts', 'ends', 'area'],
[[esc(r.size), '{0} x {1}'.format(r.rows, r.cols), fmt_in(r.top_in),
fmt_in(r.bottom_in), fmt_in(r.left_in), fmt_in(r.right_in),
fmt_in(r.start_in) + ' from start', fmt_in(r.end_in) + ' from end',
n(r.bar_area * r.rows * r.cols, 2) + ' in²'] for r in b.bars]) \
if b.bars else '<p class="note">no mild bars</p>'
strands = _table(['ids', 'size', 'pull', 'rows x cols', 'first row', 'spacing'],
[[esc(g.ids), esc(g.size), n(g.pull_pct, 0) + '%',
'{0} x {1}'.format(g.rows, g.cols), fmt_in(g.first_row_in),
fmt_in(g.spacing_in)] for g in b.strands]) \
if b.strands else '<p class="note">no strands</p>'
trans = _table(['size', 'legs', 'area', 'first', 'spacing', 'zone length', 'qty', 'symmetric'],
[[esc(z.size), str(z.legs), n(z.area, 2), fmt_in(z.first_in) +
(' from end' if z.from_end else ' from start'), fmt_in(z.spacing_in),
fmt_in(z.length_in), str(z.qty), 'yes' if z.symmetric else 'no']
for z in b.transverse]) if b.transverse else '<p class="note">no transverse steel</p>'
return (_sec('Materials', mat) + _sec('Concrete extents', geo) +
_sec('Loads (kips, klf)', loads) + _sec('Strength combinations', combos) +
_sec('Mild reinforcement', bars) + _sec('Prestress', strands) +
_sec('Transverse reinforcement', trans))
def _diagram(a: Analysis) -> str:
"""Shear and moment envelopes over the grid, plus deflection."""
b = a.beam
xs = [x / 12.0 for x in a.grid]
vmax = [0.0] * len(xs)
vmin = [0.0] * len(xs)
mmax = [0.0] * len(xs)
mmin = [0.0] * len(xs)
for tag, res in a.diagrams.items():
for i in range(len(xs)):
vmax[i] = max(vmax[i], res.shear[i])
vmin[i] = min(vmin[i], res.shear[i])
mmax[i] = max(mmax[i], res.moment[i] / 12.0)
mmin[i] = min(mmin[i], res.moment[i] / 12.0)
out = []
for title, hi, lo, unit in (('Shear envelope Vu (kips)', vmax, vmin, 'k'),
('Moment envelope Mu (kip-ft)', mmax, mmin, 'k-ft')):
p = Plot(width=760, height=200)
y0 = min(lo) if lo else 0.0
y1 = max(hi) if hi else 1.0
pad = (y1 - y0) * 0.08 or 1.0
p.domain(0.0, b.length_ft, y0 - pad, y1 + pad)
p.area(list(zip(xs, hi)), cls='fill')
p.area(list(zip(xs, lo)), cls='fill')
p.path(list(zip(xs, hi)))
p.path(list(zip(xs, lo)))
p.hline(0.0)
for s in b.supports:
p.vline(s.x_in / 12.0, cls='support')
p.axis_x()
p.hline(y1, label='{0} {1}'.format(n(y1, 1), unit))
p.hline(y0, label='{0} {1}'.format(n(y0, 1), unit))
out.append(p.render(title))
d = a.deflection
p = Plot(width=760, height=200)
tot = d['total']
live = d['live']
y1 = max(max(tot), max(live), 0.0)
y0 = min(min(tot), min(live), 0.0)
pad = (y1 - y0) * 0.1 or 0.1
p.domain(0.0, b.length_ft, -(y1 + pad), -(y0 - pad))
p.area(list(zip(xs, [-t for t in tot])), cls='fill')
p.path(list(zip(xs, [-t for t in tot])))
p.path(list(zip(xs, [-l for l in live])), cls='series2')
p.hline(0.0)
p.hline(-d['total_limit'], cls='limit', label='L/{0:.0f}'.format(b.settings.total_defl_ratio))
for s in b.supports:
p.vline(s.x_in / 12.0, cls='support')
p.axis_x()
out.append(p.render('Deflection, total (filled) and live (line), inches, down is down'))
return ''.join(out)
def _flexure_table(a: Analysis) -> str:
rows = []
for f in a.flex:
rows.append([fmt_ft_in(f.x_ft), ratio(f.util), n(f.mu, 2),
'{0} {1}'.format(n(f.phi, 2), f.controls), n(f.phi_mn, 2),
n(f.mcr12, 2), n(f.c, 2), n(f.as_t, 2), n(f.as_c, 2), n(f.aps, 3),
n(f.d, 2), esc(f.combo)])
return _table(['station', 'util', 'Mu (k-ft)', 'φ', 'φMn (k-ft)', '1.2 Mcr', 'c (in)',
'As (in²)', "A's (in²)", 'Aps (in²)', 'd (in)', 'combination'], rows)
def _shear_table(a: Analysis) -> str:
rows = []
for s in a.shr:
r = s.r
rows.append([fmt_ft_in(s.x_ft), ratio(r.util), n(r.vu, 1), n(r.phi_vn, 1), n(r.vc, 1),
n(r.vs_prov, 1), n(r.bwd, 1), n(r.av_min, 3), n(r.av_prov, 3),
n(r.av_req, 3), n(r.s, 2) if r.s else '-', n(r.s_max, 2),
n(r.f_act, 1), n(r.f_alw, 1), esc(r.note)])
return _table(['station', 'util', 'Vu (k)', 'φVn (k)', 'Vc (k)', 'Vs prov (k)', 'bw·d (in²)',
'Av min', 'Av prov', 'Av req (in²/ft)', 's (in)', 's max', 'f act (psi)',
'f alw (psi)', ''], rows)
def _service_table(a: Analysis) -> str:
b = a.beam
fc = b.materials.fc_ksi * 1000.0
lim_tot = b.settings.service_comp_total * fc
lim_sus = b.settings.service_comp_sustained * fc
rows = []
for r in a.service:
rows.append([fmt_ft_in(r.x_ft), n(r.m, 2), n(r.f_top, 0), n(r.f_bot, 0),
n(r.f_top_sus, 0), n(r.f_bot_sus, 0),
ratio(max(r.f_top, r.f_bot) / lim_tot),
ratio(max(r.f_top_sus, r.f_bot_sus) / lim_sus), esc(r.combo)])
head = ('<p class="note">Stresses on the transformed section, psi, compression positive. '
'Limits: total {0} psi (0.{1:.0f} f\'c), sustained {2} psi.</p>'
.format(n(lim_tot, 0), b.settings.service_comp_total * 100, n(lim_sus, 0)))
return head + _table(['station', 'M (k-ft)', 'f top', 'f bot', 'f top sus', 'f bot sus',
'comp / total', 'comp / sustained', 'combination'], rows)
def _deflection(a: Analysis) -> str:
d = a.deflection
return _kv([
('span', fmt_ft_in(d['span_in'] / 12.0)),
('total deflection down / up', '{0}" / {1}" (limit {2}")'.format(
n(d['total_down'], 3), n(d['total_up'], 3), n(d['total_limit'], 3))),
('live deflection down / up', '{0}" / {1}" (limit {2}")'.format(
n(d['live_down'], 3), n(d['live_up'], 3), n(d['live_limit'], 3))),
('multipliers', 'self weight x {0}, additional dead x {1}, live x 1'.format(
n(d['multipliers']['self_weight'], 2), n(d['multipliers']['addl_dead'], 2))),
('stiffness', 'E {0} ksi (self weight on {1} ksi), I {2} in⁴ (uncracked, transformed)'
.format(n(d['e'], 0), n(d['e_sw'], 0), n(d['i'], 0))),
])
def _reactions(a: Analysis) -> str:
rows = []
for case, (ra, rb) in a.reactions.items():
if ra or rb:
rows.append([esc(case), n(ra, 2), n(rb, 2)])
br = a.bearing
bear = _table(['bearing', 'Vu (k)', 'A1 (in²)', 'plain φVn (k)', 'note'],
[[k, n(v['vu'], 1), n(v['a1'], 1), n(v['phi_vn_plain'], 1), esc(v['note'])]
for k, v in br.items()])
return (_table(['case', 'start (k)', 'end (k)'], rows) +
'<h4>Bearing</h4>' + bear)
def _stripping(a: Analysis) -> str:
s = a.stripping
if not s:
return '<p class="note">stripping not analysed</p>'
b = a.beam
body = _kv([
('lift points', '{0} and {1} from the start'.format(
fmt_in(b.stripping.lift_from_start_in),
fmt_in(b.length_in - b.stripping.lift_from_end_in))),
('factored weight', '{0} kips ({1} x {2})'.format(n(s['factored_weight_k'], 2),
n(s['weight_k'], 2), b.stripping.multiplier)),
('per pick point', '{0} kips vertical; sling tension {1} k, shear {2} k at 60°'.format(
n(s['per_point_k'], 2), n(s['sling_tension_k'], 2), n(s['sling_shear_k'], 2))),
("strength used", "f'{0} = {1} ksi".format('ci' if b.stripping.strength == 'initial' else 'c',
n(s['fc_ksi'], 2))),
('tension', '{0} / {1} psi at {2}'.format(n(s['f_tension'], 0), n(s['lim_tension'], 0),
fmt_ft_in(s['f_tension_x'] / 12.0))),
('compression', '{0} / {1} psi at {2}'.format(n(s['f_comp'], 0), n(s['lim_comp'], 0),
fmt_ft_in(s['f_comp_x'] / 12.0))),
('flexure', 'Mu {0} vs φMn {1} k-ft, utilisation {2}'.format(
n(s['mu'], 1), n(s['phi_mn'], 1), ratio(s['util_flex']))),
])
res = s['result']
p = Plot(width=760, height=180)
xs = [x / 12.0 for x in res.xs]
ms = [m / 12.0 for m in res.moment]
y0, y1 = min(ms), max(ms)
pad = (y1 - y0) * 0.1 or 1.0
p.domain(0.0, b.length_ft, y0 - pad, y1 + pad)
p.area(list(zip(xs, ms)))
p.path(list(zip(xs, ms)))
p.hline(0.0)
p.axis_x()
return body + p.render('Stripping moment (kip-ft), factored self weight')
def render(a: Analysis, title: Optional[str] = None) -> str:
b = a.beam
ok = a.ok
title = title or '{0} — beam analysis'.format(b.mark or 'beam')
verdict = ('<div class="verdict {0}">{1}</div>'.format(
'ok' if ok else 'bad',
'PASS — every solved check satisfied' if ok else 'FAIL — one or more solved checks not satisfied'))
worst_f = max(a.flex, key=lambda f: f.util) if a.flex else None
worst_s = max(a.shr, key=lambda s: s.r.util) if a.shr else None
summary = _kv([
('flexure', '{0}: Mu {1} vs φMn {2} k-ft at {3} ({4})'.format(
ratio(worst_f.util), n(worst_f.mu, 1), n(worst_f.phi_mn, 1),
fmt_ft_in(worst_f.x_ft), esc(worst_f.combo)) if worst_f else '-'),
('shear', '{0}: Vu {1} vs φVn {2} k at {3}; Av/s required {4} in²/ft'.format(
ratio(worst_s.r.util), n(worst_s.r.vu, 1), n(worst_s.r.phi_vn, 1),
fmt_ft_in(worst_s.x_ft),
n(max(s.r.av_req for s in a.shr), 3)) if worst_s else '-'),
('service compression', '{0} psi total, {1} psi sustained'.format(
n(max(max(r.f_top, r.f_bot) for r in a.service), 0),
n(max(max(r.f_top_sus, r.f_bot_sus) for r in a.service), 0)) if a.service else '-'),
('service tension', '{0} psi (Class U limit {1}, Class T limit {2})'.format(
n(max(max(-r.f_top, -r.f_bot) for r in a.service), 0),
n(7.5 * (b.materials.fc_ksi * 1000.0) ** 0.5, 0),
n(12.0 * (b.materials.fc_ksi * 1000.0) ** 0.5, 0)) if a.service else '-'),
])
body = [
'<h1>{0}</h1>'.format(esc(title)),
'<p class="sub">{0} · generated {1} · beamcalc/1</p>'.format(
esc(b.description), datetime.datetime.now().strftime('%Y-%m-%d %H:%M')),
verdict, _cards(a),
_sec('Summary', summary, open_=True),
_sec('Inputs', _inputs(b)),
_sec('Stripping', _stripping(a), open_=bool(a.stripping)),
_sec('In-place strength — diagrams', _diagram(a), open_=True),
_sec('In-place strength — flexure, per station', _flexure_table(a), open_=True),
_sec('In-place strength — shear, per station', _shear_table(a), open_=True),
_sec('In-place service — stresses', _service_table(a)),
_sec('In-place service — deflection', _deflection(a), open_=True),
_sec('Reactions and bearing', _reactions(a)),
]
if a.notes:
body.append(_sec('Notes', '<ul>' + ''.join('<li>{0}</li>'.format(esc(x)) for x in a.notes) + '</ul>'))
return ('<!doctype html>\n<html><head><meta charset="utf-8"><title>{t}</title>'
'<style>{css}</style></head><body><main>{b}</main>{js}</body></html>'
.format(t=esc(title), css=CSS, b=''.join(body), js=hover_script()))
+80
View File
@@ -0,0 +1,80 @@
# -*- coding: utf-8 -*-
"""Section properties: gross, and transformed for the reinforcement.
The reference tool's "Transform Reinforcement: All" setting puts (n - 1)·A
of every bar and strand into the section before it computes the stresses
and the cracking moment -- measured on RBA5 (16" x 22", 6-#6 bottom,
4-#6 top): gross Sb = 1291 in³ gives 1.2·Mcr = 68.5 k-ft, transformed
Sb = 1493 in³ gives the printed 79.2 k-ft. So `props_at(x)` is the one to
use for anything the reference tool computes from a section modulus.
"""
from __future__ import annotations
from dataclasses import dataclass
from typing import List, Optional, Tuple
from .model import Beam, Section
@dataclass
class Props:
area: float
cy: float # centroid above the bottom
ixx: float
height: float
@property
def st(self) -> float:
return self.ixx / (self.height - self.cy)
@property
def sb(self) -> float:
return self.ixx / self.cy
@property
def ct(self) -> float:
return self.height - self.cy
def gross(section: Section) -> Props:
return Props(section.area, section.cy, section.ixx, section.height_in)
def transformed(section: Section, steel: List[Tuple[float, float, float]],
) -> Props:
"""`steel` is [(y_above_bottom, area, n), ...]; each contributes
(n - 1)·A at its position. With no steel this is the gross section."""
a = section.area
q = a * section.cy
for y, area, n in steel:
a_eq = (n - 1.0) * area
a += a_eq
q += a_eq * y
cy = q / a
ixx = 0.0
for x0, y0, w, h in section.rects:
ixx += w * h ** 3 / 12.0 + w * h * (y0 + h / 2.0 - cy) ** 2
for y, area, n in steel:
ixx += (n - 1.0) * area * (y - cy) ** 2
return Props(a, cy, ixx, section.height_in)
def props_at(beam: Beam, x_in: float, ec_ksi: Optional[float] = None,
transform: Optional[bool] = None) -> Props:
"""Section properties at x: transformed for every bar and strand
present there (a bar counts fully once it exists at x -- the reference
tool does not scale the transformed area by development), else gross."""
if transform is None:
transform = beam.settings.transform_reinforcement
if not transform:
return gross(beam.section)
ec = ec_ksi or beam.materials.Ec
steel = []
n_bar = beam.materials.rebar.es_ksi / ec
n_strand = beam.materials.strand.ep_ksi / ec
for bx, by, area, db, e0, e1, row in beam.bars_at(x_in):
steel.append((by, area, n_bar))
for g in beam.strands:
for sx, sy, area in g.strands(beam.section):
steel.append((sy, area, n_strand))
return transformed(beam.section, steel)
+178
View File
@@ -0,0 +1,178 @@
# -*- coding: utf-8 -*-
"""One-way shear, ACI 318-19 Chapter 22.5, reported per station the way
the reference tool's Shear Summary Table does: Vu, φVn, Vc, Vs provided,
bw·d, Av,min, Av provided, Av required (in²/ft), s, s max, and the shear
stresses `fact = Vu / (bw·d)` and `falw = φ·(Vc + 8·√f'c·bw·d) / (bw·d)`
(the last two reverse-engineered from RBA5 p.12: 44.1 k / 147.2 in² =
299.6 psi; 0.75·(25.8 + 8·√5000·147.2/1000) / 147.2 = 555.9 psi).
Vc:
non-prestressed, Av ≥ Av,min: Table 22.5.5.1(a): 2·λ·√f'c·bw·d
prestressed: the lesser of Vci and Vcw (22.5.6.3), with the
approximate 22.5.6.2 as the reported fallback when the detailed data
are not available.
Vs = Av·fyt·d / s, capped at 8·√f'c·bw·d (22.5.1.2).
Av,min (9.6.3.4): max(0.75·√f'c, 50)·bw / fyt.
s max (9.7.6.2.2): d/2 and 24", halved when Vs > 4·√f'c·bw·d.
"""
from __future__ import annotations
import math
from dataclasses import dataclass
from typing import List, Optional, Tuple
from .model import Beam, TransverseZone
@dataclass
class ShearResult:
vu: float # kips (factored demand, absolute)
phi_vn: float # kips
vc: float # kips
vs_prov: float # kips
bwd: float # in²
d: float # in
av_min: float # in²/ft
av_prov: float # in²/ft
av_req: float # in²/ft
s: Optional[float] # in, provided spacing at x
s_max: float # in
f_act: float # psi
f_alw: float # psi
util: float
note: str = ''
vs_over_max: bool = False # raw Vs > 8·√f'c·bw·d (22.5.1.2)
def zone_at(beam: Beam, x_in: float) -> Tuple[Optional[TransverseZone], float]:
"""(zone, spacing) of the transverse steel covering x, or (None, 0)."""
L = beam.length_in
best = None
for z in beam.transverse:
if z.use.lower() != 'av':
continue
spans = []
x0 = (L - z.first_in) if z.from_end else z.first_in
length = z.length_in if z.length_in else (z.qty - 1) * z.spacing_in
if z.from_end:
spans.append((x0 - length - z.spacing_in / 2.0, x0 + z.spacing_in / 2.0))
else:
spans.append((x0 - z.spacing_in / 2.0, x0 + length + z.spacing_in / 2.0))
if z.symmetric:
a, b = spans[0]
spans.append((L - b, L - a))
# a zone counts once at x if x lies in ANY of its spans (its own
# direct span or the mirrored span -- a symmetric zone at midspan
# where the two overlap must not be counted twice)
if any(a - 1e-9 <= x_in <= b + 1e-9 for a, b in spans):
if best is None or z.spacing_in < best.spacing_in:
best = z
return (best, best.spacing_in) if best else (None, 0.0)
def av_provided_per_ft(beam: Beam, x_in: float) -> Tuple[float, Optional[float]]:
"""Total Av at x (in²/ft) summing every zone that covers x -- two
overlapping zones (the 'overlap stirrups' detail) add up -- and the
tightest spacing among them."""
L = beam.length_in
total = 0.0
s_min = None
for z in beam.transverse:
if z.use.lower() != 'av':
continue
x0 = (L - z.first_in) if z.from_end else z.first_in
length = z.length_in if z.length_in else (z.qty - 1) * z.spacing_in
spans = []
if z.from_end:
spans.append((x0 - length - z.spacing_in / 2.0, x0 + z.spacing_in / 2.0))
else:
spans.append((x0 - z.spacing_in / 2.0, x0 + length + z.spacing_in / 2.0))
if z.symmetric:
a, b = spans[0]
spans.append((L - b, L - a))
# a zone counts once at x if x lies in ANY of its spans (its own
# direct span or the mirrored span -- a symmetric zone at midspan
# where the two overlap must not be counted twice). Distinct zones
# stay additive: an inner/outer doubled stirrup is two zone rows
# and each counts once.
if any(a - 1e-9 <= x_in <= b + 1e-9 for a, b in spans):
total += z.area * 12.0 / z.spacing_in
s_min = z.spacing_in if s_min is None else min(s_min, z.spacing_in)
return total, s_min
def vc_nonprestressed(beam: Beam, bw: float, d: float, fc_ksi: float,
nu_k: float = 0.0) -> float:
lam = beam.materials.lam
fc_psi = fc_ksi * 1000.0
vc = 2.0 * lam * math.sqrt(fc_psi) / 1000.0 * bw * d
if nu_k:
vc += nu_k / (6.0 * beam.section.area) * bw * d
return max(0.0, min(vc, 5.0 * lam * math.sqrt(fc_psi) / 1000.0 * bw * d))
def vc_prestressed(beam: Beam, bw: float, d: float, dp: float, fc_ksi: float,
vu: float, mu_kin: float, mmax_kin: float, vd: float,
fpe_psi_bot: float, fd_psi_bot: float, fpc_psi: float,
yt: float, i_in4: float, vp: float = 0.0) -> Tuple[float, str]:
"""Detailed Vci / Vcw (22.5.6.3.1 / 22.5.6.3.2); returns (Vc, which)."""
lam = beam.materials.lam
sq = math.sqrt(fc_ksi * 1000.0)
dv = max(dp, 0.8 * beam.section.height_in)
mcre = (i_in4 / yt) * (6.0 * lam * sq + fpe_psi_bot - fd_psi_bot) / 1000.0
if abs(mmax_kin) > 1e-9:
vci = 0.6 * lam * sq / 1000.0 * bw * dv + vd + vu * mcre / abs(mmax_kin)
else:
vci = float('inf')
vci = max(vci, 1.7 * lam * sq / 1000.0 * bw * dv)
vcw = (3.5 * lam * sq + 0.3 * fpc_psi) / 1000.0 * bw * dv + vp
if vci <= vcw:
return vci, 'Vci'
return vcw, 'Vcw'
def check(beam: Beam, x_in: float, vu: float, d: float, fc_ksi: float,
prestress: Optional[dict] = None, mu_kin: float = 0.0,
mmax_kin: float = 0.0, vd: float = 0.0) -> ShearResult:
mats = beam.materials
lam = mats.lam
fc_psi = fc_ksi * 1000.0
sq = math.sqrt(fc_psi)
bw = beam.section.web_width()
fyt = mats.trans.fy_ksi
phi = beam.settings.phi_shear
d = max(d, 1e-6)
bwd = bw * d
av_prov, s = av_provided_per_ft(beam, x_in)
av_min = max(0.75 * sq, 50.0) * bw / (fyt * 1000.0) * 12.0
if prestress:
vc, which = vc_prestressed(beam, bw, d, prestress['dp'], fc_ksi,
abs(vu), mu_kin, mmax_kin, vd,
prestress['fpe_bot'], prestress['fd_bot'],
prestress['fpc'], prestress['yt'],
prestress['i'], prestress.get('vp', 0.0))
note = which
else:
vc = vc_nonprestressed(beam, bw, d, fc_ksi)
note = ''
vs_cap = 8.0 * lam * sq / 1000.0 * bwd
vs_prov = av_prov / 12.0 * fyt * d if av_prov else 0.0
vs_over_max = vs_prov > vs_cap
if vs_over_max:
note = (note + '; ' if note else '') + "Vs exceeds 8*sqrt(f'c)*bw*d (22.5.1.2)"
phi_vn = phi * (vc + vs_prov)
vu_a = abs(vu)
# required Av/s (in²/ft): (Vu/φ - Vc)/(fyt d) * 12, not below min when
# Vu > φVc/2 (9.6.3.1), zero when Vu <= φVc/2
vs_req = max(0.0, vu_a / phi - vc)
av_req = vs_req / (fyt * d) * 12.0 if vs_req > 0 else 0.0
if vu_a > 0.5 * phi * vc:
av_req = max(av_req, av_min)
s_max = min(d / 2.0, 24.0)
if vs_prov > 4.0 * lam * sq / 1000.0 * bwd:
s_max = min(d / 4.0, 12.0)
f_act = vu_a * 1000.0 / bwd if bwd else 0.0
f_alw = phi * (vc + vs_cap) * 1000.0 / bwd if bwd else 0.0
util = vu_a / phi_vn if phi_vn > 0 else (0.0 if vu_a == 0 else float('inf'))
return ShearResult(vu_a, phi_vn, vc, vs_prov, bwd, d, av_min, av_prov,
av_req, s, s_max, f_act, f_alw, util, note, vs_over_max)
+573
View File
@@ -0,0 +1,573 @@
# -*- coding: utf-8 -*-
"""Run every check on a beam and return one `Analysis`.
Order, and what feeds what:
1. stations -- L/20 grid + supports, support faces, critical sections,
lift points, load and bar ends, the maximum-moment location.
2. envelope -- every included strength combo × every live pattern × every
wind variant, solved with `statics.solve`; per station the max sagging
Mu, max hogging Mu and max |Vu| with the combo that made each.
3. flexure per station (`flexure.capacity_at`) for both signs.
4. shear per station (`shear.check`) with d from the flexure solve.
5. service -- every service combo, stresses on the transformed section.
6. deflection -- elastic per case, then the reference tool's multipliers.
7. stripping -- a two-point pick on the initial strength.
8. reactions, bearing.
Every number in the result carries the station it belongs to so the report
and the parity tests can address it by `(station_ft, column)`.
"""
from __future__ import annotations
import math
from dataclasses import dataclass, field, asdict
from typing import Dict, List, Optional, Tuple
from . import deflection as defl_mod, flexure, loads as ld, shear, section as sec_mod, statics, develop
from .model import Beam, Combo
@dataclass
class StationStrength:
x_in: float
mu_pos: float = 0.0 # kip-ft, max sagging
mu_neg: float = 0.0 # kip-ft, max hogging (negative number)
vu: float = 0.0 # kips, max |V|
vu_signed: float = 0.0
combo_pos: str = ''
combo_neg: str = ''
combo_v: str = ''
mu_at_vmax: float = 0.0 # moment coincident with the governing shear
vd_at_vmax: float = 0.0 # unfactored dead shear (prestressed Vci)
@dataclass
class StationFlexure:
x_in: float
x_ft: float
mu: float
phi_mn: float
mcr12: float
util: float
phi: float
c: float
eps_t: float
as_t: float
as_c: float
aps: float
d: float
controls: str
positive: bool
combo: str
strain_ok: bool
@dataclass
class StationShear:
x_in: float
x_ft: float
r: shear.ShearResult
combo: str
@dataclass
class StationService:
x_in: float
x_ft: float
f_top: float # psi, + compression
f_bot: float
f_top_sus: float
f_bot_sus: float
combo: str
m: float # kip-ft governing
@dataclass
class Stage:
name: str
ok: bool
util: float
lines: List[str] = field(default_factory=list)
data: dict = field(default_factory=dict)
@dataclass
class Analysis:
beam: Beam
stations_in: List[float]
grid: List[float]
strength: List[StationStrength]
flex: List[StationFlexure]
shr: List[StationShear]
service: List[StationService]
stages: List[Stage]
reactions: dict
deflection: dict
stripping: Optional[dict]
bearing: dict
diagrams: dict # combo name -> statics.Result (for plots)
notes: List[str] = field(default_factory=list)
@property
def ok(self) -> bool:
return all(s.ok for s in self.stages)
def flex_at(self, x_ft: float) -> Optional[StationFlexure]:
return _nearest(self.flex, x_ft)
def shear_at(self, x_ft: float) -> Optional[StationShear]:
return _nearest(self.shr, x_ft)
def _nearest(rows, x_ft, tol=0.02):
best = None
for r in rows:
if abs(r.x_ft - x_ft) <= tol and (best is None or
abs(r.x_ft - x_ft) < abs(best.x_ft - x_ft)):
best = r
return best
# ---------------------------------------------------------------- stations
def stations(beam: Beam, d_at_face: Optional[Dict[str, float]] = None
) -> List[float]:
L = beam.length_in
pts = set()
n = 20
for i in range(n + 1):
pts.add(L * i / n)
for s in beam.supports:
pts.add(s.x_in)
pts.add(s.face_start_in)
pts.add(s.face_end_in)
if d_at_face:
pts.add(beam.supports[0].face_end_in + d_at_face.get('start', 0.0))
pts.add(beam.supports[1].face_start_in - d_at_face.get('end', 0.0))
for p in beam.point_loads:
pts.add(p.x_in)
for l in beam.line_loads:
pts.add(l.x0_in)
pts.add(l.x1_in)
for row in beam.bars:
pts.add(row.start_in)
pts.add(L - row.end_in)
if beam.stripping:
pts.add(beam.stripping.lift_from_start_in)
pts.add(L - beam.stripping.lift_from_end_in)
return sorted(round(p, 4) for p in pts if 0.0 <= p <= L)
def breakpoints(beam: Beam) -> List[float]:
pts = [s.x_in for s in beam.supports]
pts += [p.x_in for p in beam.point_loads]
for l in beam.line_loads:
pts += [l.x0_in, l.x1_in]
pts += [m.x_in for m in beam.point_moments]
return pts
# ---------------------------------------------------------------- envelope
def _combo_loadsets(beam: Beam, combo: Combo):
"""Every (label, loads) the combo expands to over patterns and wind."""
dead = ld.dead_loads(beam)
out = []
pats = ld.patterns(beam) if combo.L or combo.R else [tuple(ld.regions(beam))]
winds = ld.wind_variants(beam) if combo.W else [('', [])]
for pat in pats:
live = ld.live_loads_for(beam, pat, 'live') if combo.L else []
roof = ld.live_loads_for(beam, pat, 'roof') if combo.R else []
for wname, wl in winds:
loads = (ld.scaled(dead, combo.D) + ld.scaled(live, combo.L)
+ ld.scaled(roof, combo.R) + ld.scaled(wl, combo.W)
+ ld.scaled(ld.case_loads(beam, 'seismic'), combo.E))
tag = combo.name()
if len(pats) > 1:
tag += ' [pattern {0}]'.format(len(pat))
if wname:
tag += ' ' + wname
out.append((tag, loads))
return out
def envelope(beam: Beam, combos: List[Combo], xs: List[float]
) -> Tuple[List[StationStrength], Dict[str, statics.Result]]:
xa, xb = beam.supports[0].x_in, beam.supports[1].x_in
rows = [StationStrength(x) for x in xs]
diagrams = {}
dead_only = statics.solve(ld.dead_loads(beam), xa, xb, xs)
for combo in combos:
if not combo.include:
continue
for tag, loads in _combo_loadsets(beam, combo):
res = statics.solve(loads, xa, xb, xs)
diagrams[tag] = res
for i, x in enumerate(xs):
m = res.moment[i] / 12.0
vl, vr = res.shear_jump(x)
v = _station_shear(beam, x, vl, vr)
r = rows[i]
if m > r.mu_pos:
r.mu_pos, r.combo_pos = m, tag
if m < r.mu_neg:
r.mu_neg, r.combo_neg = m, tag
if abs(v) > r.vu:
r.vu, r.vu_signed, r.combo_v = abs(v), v, tag
r.mu_at_vmax = m
dl, dr = dead_only.shear_jump(x)
r.vd_at_vmax = dl if abs(vl) >= abs(vr) else dr
return rows, diagrams
# ---------------------------------------------------------------- solve
def analyse(beam: Beam, extra_stations_in: Optional[List[float]] = None) -> Analysis:
"""`extra_stations_in`: more x (in) to report at -- the comparison tool passes
the reference package's printed stations so every row lines up."""
notes: List[str] = []
mats = beam.materials
xa, xb = beam.supports[0].x_in, beam.supports[1].x_in
L = beam.length_in
# --- effective depth at the support faces, for the critical sections
d_face = {}
for key, x in (('start', beam.supports[0].face_end_in),
('end', beam.supports[1].face_start_in)):
fr = flexure.capacity_at(beam, x, positive=True)
d_face[key] = fr.d if fr.d > 0 else 0.8 * beam.section.height_in
xs = stations(beam, d_face)
if extra_stations_in:
xs = sorted(set(xs) | {round(x, 4) for x in extra_stations_in if 0.0 <= x <= L})
grid = statics.grid(L, breakpoints(beam) + xs, step=0.5)
# --- strength envelope on the grid (so plots are smooth), sampled at xs
strength_grid, diagrams = envelope(beam, beam.strength_combos, grid)
by_x = {round(r.x_in, 4): r for r in strength_grid}
strength = [by_x[round(x, 4)] for x in xs]
# --- flexure per station
flex_rows: List[StationFlexure] = []
fse = _effective_prestress_ksi(beam) if beam.strands else 0.0
for r in strength:
x = r.x_in
best = None
for positive, mu, combo in ((True, r.mu_pos, r.combo_pos),
(False, r.mu_neg, r.combo_neg)):
cap = flexure.capacity_at(beam, x, positive=positive, fse_ksi=fse)
mcr = flexure.cracking_moment(beam, x, positive=positive)
if cap.phi_mn > 1e-9:
util = abs(mu) / cap.phi_mn
elif abs(mu) <= max(0.1, 0.01 * abs(1.2 * mcr)):
# no developed steel and next to no moment: the reference
# tool prints NC here (the beam ends, the bearing zone)
util = 0.0
else:
util = float('inf')
row = StationFlexure(x, x / 12.0, mu, cap.phi_mn if positive else -cap.phi_mn,
1.2 * mcr if positive else -1.2 * mcr, util,
cap.phi, cap.c, cap.eps_t, cap.as_tension,
cap.as_comp, cap.aps, cap.d, cap.controls,
positive, combo, cap.strain_ok)
if best is None or util > best.util or (
util == best.util and positive):
best = row
flex_rows.append(best)
# --- shear per station
shr_rows: List[StationShear] = []
for r, f in zip(strength, flex_rows):
x = r.x_in
d = _shear_depth(beam, x, f)
# The reference tool prints the actual Vu at every station and
# lets the critical section (face + d) be one more station; it
# does not clip the end-region demand to it. Same here.
vu = r.vu
res = shear.check(beam, x, vu, d, mats.fc_ksi)
shr_rows.append(StationShear(x, x / 12.0, res, r.combo_v))
# --- service stresses
service_rows = _service(beam, xs, grid)
# --- deflections
defl = defl_mod.analyse(beam, grid)
# --- stripping
strip = _stripping(beam) if beam.stripping else None
# --- reactions per case (unfactored)
reactions = _reactions(beam)
# --- bearing
bearing = _bearing(beam, diagrams)
# --- stage summaries
stages = []
if strip:
stages.append(Stage('Stripping', strip['ok'], strip['util'],
strip['lines'], strip))
worst_f = max(flex_rows, key=lambda f: f.util) if flex_rows else None
worst_s = max(shr_rows, key=lambda s: s.r.util) if shr_rows else None
flex_ok = all(f.util <= 1.0 + 1e-9 and f.strain_ok for f in flex_rows)
mcr_ok = all(abs(f.phi_mn) >= abs(f.mcr12) - 1e-9 or abs(f.mu) < 1e-9
for f in flex_rows)
stages.append(Stage('In-Place Strength',
flex_ok and (worst_s is None or worst_s.r.util <= 1.0 + 1e-9),
max(worst_f.util if worst_f else 0.0,
worst_s.r.util if worst_s else 0.0),
[], {'flexure': worst_f, 'shear': worst_s, 'mcr_ok': mcr_ok}))
svc_ok = defl['total_ok'] and defl['live_ok'] and _service_ok(beam, service_rows)
stages.append(Stage('In-Place Service', svc_ok,
max(defl['total_util'], defl['live_util']), [], {}))
return Analysis(beam, xs, grid, strength, flex_rows, shr_rows, service_rows,
stages, reactions, defl, strip, bearing, diagrams, notes)
def _station_shear(beam: Beam, x: float, vl: float, vr: float) -> float:
"""Which side of a shear jump a station reports. The reference tool
prints the OVERHANG side at a bearing centreline (RBA5 p.12: -0.9 k at
the 2" support, not the 43.9 k just inside it) -- the span-side value
belongs to the next station in. Everywhere else the larger magnitude."""
xa, xb = beam.supports[0].x_in, beam.supports[1].x_in
if abs(x - xa) < 1e-6:
return vl
if abs(x - xb) < 1e-6:
return vr
return vl if abs(vl) >= abs(vr) else vr
def _vu_at(rows: List[StationStrength], x: float) -> float:
best = min(rows, key=lambda r: abs(r.x_in - x))
return best.vu
def _shear_depth(beam: Beam, x: float, f: StationFlexure) -> float:
"""d for shear: the reference tool's 'CG of Tension Rein.' -- the
flexure solve's d (from the compression face) but not less than
0.8h for a prestressed member (22.5.2.1). Where no bar has any
development yet (the bearing centreline, As = 0) the reference tool
still prints bw·d with d at the plain area centroid of every bar
present, measured from the compression face of the governing moment's
sign at that station -- RBA5 p.12: bwd = 204.8 = 16 × 12.80 at x = 0
(from the top), 147.2 = 16 × 9.20 at x = 2.04 in just inside the
bearing where the hogging overhang moment governs (from the bottom),
never 16 × 0.8h. The bars are taken where they exist (at x, else at
midspan: at x = 0 the bars have not started yet and the tool still
prints 12.80); 0.8h only when the beam has no bars at all."""
h = beam.section.height_in
d = f.d
if d <= 0:
bars = beam.bars_at(x) or beam.bars_at(beam.length_in / 2.0)
tot = sum(b[2] for b in bars)
if tot > 0:
depth = (lambda by: (h - by) if f.positive else by)
d = sum(b[2] * depth(b[1]) for b in bars) / tot
else:
d = 0.8 * h
if beam.strands:
d = max(d, 0.8 * h)
return d
def _effective_prestress_ksi(beam: Beam) -> float:
"""Effective strand stress after losses. Until the PCI loss routine
lands, a flat 18 % of the jacking stress (PCI 8th Ed. typical for
pretensioned members), applied to the highest pull."""
if not beam.strands:
return 0.0
fpu = beam.materials.strand.fpu_ksi
pull = max(g.pull_pct for g in beam.strands) / 100.0
return fpu * pull * 0.82
def _service(beam: Beam, xs: List[float], grid: List[float]
) -> List[StationService]:
xa, xb = beam.supports[0].x_in, beam.supports[1].x_in
rows = []
dead = ld.dead_loads(beam)
sus = ld.case_loads(beam, 'sus_live')
results = []
for combo in beam.service_combos:
if not combo.include:
continue
for tag, loads in _combo_loadsets(beam, combo):
results.append((tag, statics.solve(loads, xa, xb, xs)))
sus_res = statics.solve(dead + sus, xa, xb, xs)
for i, x in enumerate(xs):
p = sec_mod.props_at(beam, x)
best = None
for tag, res in results:
m = res.moment[i] # kip-in
f_top = m / p.st * 1000.0 # psi, + compression at top for M>0
f_bot = -m / p.sb * 1000.0
key = max(-f_top, -f_bot, f_top, f_bot)
if best is None or key > best[0]:
best = (key, tag, m, f_top, f_bot)
m_s = sus_res.moment[i]
f_top_s = m_s / p.st * 1000.0
f_bot_s = -m_s / p.sb * 1000.0
if best is None:
best = (0.0, '', 0.0, 0.0, 0.0)
rows.append(StationService(x, x / 12.0, best[3], best[4], f_top_s,
f_bot_s, best[1], best[2] / 12.0))
return rows
def _service_ok(beam: Beam, rows: List[StationService]) -> bool:
fc = beam.materials.fc_ksi * 1000.0
lim_tot = beam.settings.service_comp_total * fc
lim_sus = beam.settings.service_comp_sustained * fc
for r in rows:
if max(r.f_top, r.f_bot) > lim_tot + 1e-6:
return False
if max(r.f_top_sus, r.f_bot_sus) > lim_sus + 1e-6:
return False
return True
def _deflections(beam: Beam, grid: List[float]) -> dict:
"""Elastic deflection per case on the gross (transformed) section,
combined with the reference tool's multipliers:
total = SW·m_sw + (deadnc + topping + deadc)·m_dead + live + roof
live = live (+ roof)
Self weight uses Eci when the setting says so; everything else Ec.
Limits are span / ratio with span = distance between bearings."""
xa, xb = beam.supports[0].x_in, beam.supports[1].x_in
s = beam.settings
mats = beam.materials
span = xb - xa
p = sec_mod.props_at(beam, (xa + xb) / 2.0)
ei_c = [mats.Ec * p.ixx] * len(grid)
e_sw = mats.Eci if s.ec_for_sw_defl.lower() == 'eci' else mats.Ec
ei_i = [e_sw * p.ixx] * len(grid)
def defl(loads, ei):
if not loads:
return [0.0] * len(grid)
res = statics.solve(loads, xa, xb, grid)
return statics.deflect(res, ei, xa, xb)
d_sw = defl(ld.self_weight(beam), ei_i)
d_dead = defl(sum((ld.case_loads(beam, c) for c in ld.DEAD_CASES), []), ei_c)
d_live = defl(ld.case_loads(beam, 'live'), ei_c)
d_roof = defl(ld.case_loads(beam, 'roof'), ei_c)
m_sw = s.sw_mult.get('end', 2.0)
m_dd = s.addl_dead_mult.get('end', 2.0)
total = [d_sw[i] * m_sw + d_dead[i] * m_dd + d_live[i] + d_roof[i]
for i in range(len(grid))]
live = [d_live[i] + d_roof[i] for i in range(len(grid))]
# report within the span (overhang tips separately)
idx = [i for i, x in enumerate(grid) if xa <= x <= xb]
t_down = max(total[i] for i in idx)
t_up = min(total[i] for i in idx)
l_down = max(live[i] for i in idx)
l_up = min(live[i] for i in idx)
lim_t = span / s.total_defl_ratio
lim_l = span / s.live_defl_ratio
return {
'span_in': span,
'total_down': t_down, 'total_up': -t_up, 'total_limit': lim_t,
'live_down': l_down, 'live_up': -l_up, 'live_limit': lim_l,
'total_util': max(t_down, -t_up) / lim_t if lim_t else 0.0,
'live_util': max(l_down, -l_up) / lim_l if lim_l else 0.0,
'total_ok': max(t_down, -t_up) <= lim_t + 1e-9,
'live_ok': max(l_down, -l_up) <= lim_l + 1e-9,
'elastic': {'self_weight': d_sw, 'dead': d_dead, 'live': d_live,
'roof': d_roof},
'total': total, 'live': live,
'multipliers': {'self_weight': m_sw, 'addl_dead': m_dd},
'e_sw': e_sw, 'e': mats.Ec, 'i': p.ixx,
}
def _stripping(beam: Beam) -> dict:
st = beam.stripping
mats = beam.materials
L = beam.length_in
xa = st.lift_from_start_in
xb = L - st.lift_from_end_in
grid = statics.grid(L, [xa, xb], step=0.5)
loads = ld.scaled(ld.self_weight(beam), st.multiplier)
res = statics.solve(loads, xa, xb, grid)
fc = mats.fci_ksi if st.strength == 'initial' else mats.fc_ksi
fc_psi = fc * 1000.0
s = beam.settings
lim_c = s.transfer_comp_other * fc_psi
lim_t = s.transfer_tens_other_sqrt * math.sqrt(fc_psi)
worst_t = (0.0, 0.0)
worst_c = (0.0, 0.0)
for i, x in enumerate(grid):
p = sec_mod.props_at(beam, x)
m = res.moment[i]
f_top = m / p.st * 1000.0
f_bot = -m / p.sb * 1000.0
for f in (f_top, f_bot):
if f > worst_c[0]:
worst_c = (f, x)
if -f > worst_t[0]:
worst_t = (-f, x)
# flexure at the max-moment station with the initial strength
xm, mm = res.max_abs_moment()
cap = flexure.capacity_at(beam, xm, positive=(mm >= 0), fc_ksi=fc)
mu = abs(mm) / 12.0
util_f = mu / cap.phi_mn if cap.phi_mn > 1e-9 else float('inf')
w = beam.weight_k * st.multiplier
per_point = w / st.pick_points
sling = math.radians(60.0)
tension = per_point
shear_k = per_point / math.tan(sling)
util = max(worst_t[0] / lim_t if lim_t else 0.0,
worst_c[0] / lim_c if lim_c else 0.0, util_f)
lines = [
'Lifted at {0:.1f}" and {1:.1f}" from the start, self weight x {2:g}'
.format(xa, xb, st.multiplier),
'tension {0:.0f} / {1:.0f} psi, compression {2:.0f} / {3:.0f} psi'
.format(worst_t[0], lim_t, worst_c[0], lim_c),
'flexure Mu {0:.1f} vs phiMn {1:.1f} k-ft'.format(mu, cap.phi_mn),
]
return {
'ok': util <= 1.0 + 1e-9, 'util': util, 'lines': lines,
'f_tension': worst_t[0], 'f_tension_x': worst_t[1], 'lim_tension': lim_t,
'f_comp': worst_c[0], 'f_comp_x': worst_c[1], 'lim_comp': lim_c,
'mu': mu, 'phi_mn': cap.phi_mn, 'util_flex': util_f, 'x_max': xm,
'weight_k': beam.weight_k, 'factored_weight_k': w,
'per_point_k': per_point, 'sling_tension_k': tension,
'sling_shear_k': shear_k, 'result': res, 'fc_ksi': fc,
}
def _reactions(beam: Beam) -> dict:
xa, xb = beam.supports[0].x_in, beam.supports[1].x_in
out = {}
for case in (ld.SELF,) + ld.CASES:
loads = ld.case_loads(beam, case)
if not loads:
out[case] = (0.0, 0.0)
continue
ra, rb = statics.reactions(loads, xa, xb)
out[case] = (ra, rb)
return out
def _bearing(beam: Beam, diagrams: dict) -> dict:
"""Factored reactions (max over the strength combos) and the PCI
plain-bearing capacity. Vu is the REACTION, read off each combo's own
diagram (`Result.r_start` / `r_end`), never the station shear at the
bearing centreline: that station reports the overhang side of the jump
(`_station_shear`), which is ~0 -- RBA5 prints 45.1 k here against a
-0.9 k shear row at the same x. The Cr reduction is not yet reproduced;
`phi_vn_plain` is the unreduced φ·0.85·f'c·A1 and is flagged as such."""
vu_l = max((abs(res.r_start) for res in diagrams.values()), default=0.0)
vu_r = max((abs(res.r_end) for res in diagrams.values()), default=0.0)
out = {'left': {}, 'right': {}}
for key, sup, vu in (('left', beam.supports[0], vu_l),
('right', beam.supports[1], vu_r)):
a1 = sup.length_in * sup.width_in
phi_vn = beam.settings.phi_compression * 0.85 * beam.materials.fc_ksi * a1
out[key] = {'vu': vu, 'a1': a1, 'phi_vn_plain': phi_vn,
'note': 'unreduced (Cr = 1); PCI 5.2.1 reduction pending'}
return out
+243
View File
@@ -0,0 +1,243 @@
# -*- coding: utf-8 -*-
"""Statics of a beam on two supports with overhangs.
Determinate, so reactions come from equilibrium and V(x), M(x) from a single
pass along a fine grid. Deflection is a double integration of M/EI with
w = 0 at both supports, so it is right for any load shape and for a
varying EI (the cracked-section case in `deflection.py`).
Sign convention: loads are positive DOWN, shear V(x) is the sum of the
upward forces to the left of x, moment M(x) is positive SAGGING (bottom in
tension). Deflection is positive DOWN. All lengths inches, forces kips,
moments kip-in.
"""
from __future__ import annotations
import bisect
import math
from dataclasses import dataclass, field
from typing import Dict, List, Optional, Sequence, Tuple
@dataclass
class Load:
"""One load on the beam in one case: a point force P at x, a linearly
varying line load (w0 at x0, w1 at x1, kips/in), or a point moment M
at x (kip-in, positive clockwise = sagging to the right)."""
kind: str # 'point' | 'line' | 'moment'
x0: float
x1: float = 0.0
p: float = 0.0 # kips (point)
w0: float = 0.0 # k/in (line)
w1: float = 0.0
m: float = 0.0 # kip-in (moment)
tag: str = ''
def total(self) -> float:
if self.kind == 'point':
return self.p
if self.kind == 'line':
return (self.w0 + self.w1) / 2.0 * (self.x1 - self.x0)
return 0.0
def centroid(self) -> float:
if self.kind == 'point':
return self.x0
if self.kind == 'line':
L = self.x1 - self.x0
if L <= 0:
return self.x0
tot = self.total()
if abs(tot) < 1e-12:
return (self.x0 + self.x1) / 2.0
# first moment of a trapezoid about x0
m1 = self.w0 * L * L / 2.0 + (self.w1 - self.w0) * L * L / 3.0
return self.x0 + m1 / tot
return self.x0
def moment_about(self, x: float) -> float:
"""Moment of this load about x (kip-in, clockwise positive for a
downward load to the right of x)."""
if self.kind == 'moment':
return self.m
return self.total() * (self.centroid() - x)
@dataclass
class Result:
xs: List[float]
shear: List[float]
moment: List[float]
r_start: float
r_end: float
deflection: Optional[List[float]] = None
loads: Optional[List] = None
xa: float = 0.0
xb: float = 0.0
def at(self, x: float, series: Sequence[float]) -> float:
return interp(self.xs, series, x)
def moment_at(self, x: float) -> float:
return self.at(x, self.moment)
def shear_at(self, x: float) -> float:
return self.at(x, self.shear)
def shear_jump(self, x: float) -> Tuple[float, float]:
"""(V just left, V just right) -- a point load or support makes V
discontinuous and the printed 'V at x' is whichever side governs."""
if self.loads is not None:
return (_vm_at(self.loads, self.xa, self.xb, self.r_start,
self.r_end, x, inclusive=False)[0],
_vm_at(self.loads, self.xa, self.xb, self.r_start,
self.r_end, x, inclusive=True)[0])
# Fall back to the neighbouring grid nodes for a Result that was not
# built by `solve` and carries no loads to recompute from.
i = bisect.bisect_left(self.xs, x - 1e-9)
j = bisect.bisect_right(self.xs, x + 1e-9)
left = self.shear[max(i - 1, 0)]
right = self.shear[min(j, len(self.xs) - 1)]
return left, right
def deflection_at(self, x: float) -> float:
if self.deflection is None:
return 0.0
return self.at(x, self.deflection)
def max_abs_moment(self) -> Tuple[float, float]:
i = max(range(len(self.xs)), key=lambda k: abs(self.moment[k]))
return self.xs[i], self.moment[i]
def interp(xs: Sequence[float], ys: Sequence[float], x: float) -> float:
if x <= xs[0]:
return ys[0]
if x >= xs[-1]:
return ys[-1]
i = bisect.bisect_right(xs, x)
x0, x1 = xs[i - 1], xs[i]
if x1 == x0:
return ys[i]
t = (x - x0) / (x1 - x0)
return ys[i - 1] + t * (ys[i] - ys[i - 1])
def grid(length: float, breakpoints: Sequence[float], step: float = 0.5
) -> List[float]:
"""A fine grid from 0 to `length` that lands exactly on every
breakpoint (supports, load ends, bar ends), so V and M kinks are
represented and nothing is interpolated across a jump."""
pts = set([0.0, length])
for b in breakpoints:
if 0.0 <= b <= length:
pts.add(round(b, 6))
n = int(math.ceil(length / step))
for i in range(n + 1):
pts.add(round(min(length, i * step), 6))
return sorted(pts)
def reactions(loads: Sequence[Load], xa: float, xb: float) -> Tuple[float, float]:
"""Upward reactions at supports xa < xb for downward loads."""
tot = sum(l.total() for l in loads)
m_about_a = sum(l.moment_about(xa) for l in loads)
rb = m_about_a / (xb - xa)
ra = tot - rb
return ra, rb
def _line_w_at(l: Load, x: float) -> float:
if x < l.x0 or x > l.x1 or l.x1 <= l.x0:
return 0.0
t = (x - l.x0) / (l.x1 - l.x0)
return l.w0 + t * (l.w1 - l.w0)
def _vm_at(loads: Sequence[Load], xa: float, xb: float, ra: float, rb: float,
x: float, inclusive: bool) -> Tuple[float, float]:
"""Closed-form (V, M) at x for the loads and reactions ra, rb.
`inclusive=True` includes a reaction, point load or applied moment
located exactly at x (V taken just to the RIGHT of x); `inclusive=False`
excludes it, giving the LEFT limit. Line loads are continuous and are
treated the same either way."""
v = 0.0
m = 0.0
# reactions (upward)
if (xa <= x + 1e-9) if inclusive else (xa < x - 1e-9):
v += ra
m += ra * (x - xa)
if (xb <= x + 1e-9) if inclusive else (xb < x - 1e-9):
v += rb
m += rb * (x - xb)
for l in loads:
if l.kind == 'point':
if (l.x0 <= x + 1e-9) if inclusive else (l.x0 < x - 1e-9):
v -= l.p
m -= l.p * (x - l.x0)
elif l.kind == 'line':
if x > l.x0:
xe = min(x, l.x1)
L = xe - l.x0
if L > 0:
we = _line_w_at(l, xe)
tot = (l.w0 + we) / 2.0 * L
# centroid of the trapezoid from x0 to xe
if abs(tot) > 1e-12:
m1 = l.w0 * L * L / 2.0 + (we - l.w0) * L * L / 3.0
c = l.x0 + m1 / tot
else:
c = (l.x0 + xe) / 2.0
v -= tot
m -= tot * (x - c)
elif l.kind == 'moment':
if (l.x0 <= x + 1e-9) if inclusive else (l.x0 < x - 1e-9):
m -= l.m
return v, m
def solve(loads: Sequence[Load], xa: float, xb: float, xs: Sequence[float]
) -> Result:
"""V(x) and M(x) on the grid `xs` for the given loads and supports.
Both are evaluated exactly (closed form per load) at every grid point
-- the grid is for reporting, not for accuracy -- with V taken just to
the RIGHT of x (so at a support or point load the grid point carries
the post-jump value; `shear_jump` gives both sides)."""
ra, rb = reactions(loads, xa, xb)
shear = []
moment = []
for x in xs:
v, m = _vm_at(loads, xa, xb, ra, rb, x, inclusive=True)
shear.append(v)
moment.append(m)
return Result(list(xs), shear, moment, ra, rb, loads=list(loads),
xa=xa, xb=xb)
def deflect(result: Result, ei: Sequence[float], xa: float, xb: float
) -> List[float]:
"""Deflection (positive DOWN, inches) by double integration of M/EI
with w(xa) = w(xb) = 0. `ei` is EI (kip-in²) at every grid point --
constant, or varying for a cracked section."""
xs = result.xs
n = len(xs)
curv = [result.moment[i] / ei[i] if ei[i] else 0.0 for i in range(n)]
# theta = integral of curvature (trapezoid); w = integral of theta
theta = [0.0] * n
for i in range(1, n):
h = xs[i] - xs[i - 1]
theta[i] = theta[i - 1] + h * (curv[i] + curv[i - 1]) / 2.0
w = [0.0] * n
for i in range(1, n):
h = xs[i] - xs[i - 1]
w[i] = w[i - 1] + h * (theta[i] + theta[i - 1]) / 2.0
# w_true = w + c1*x + c2 ; impose zero at supports
wa = interp(xs, w, xa)
wb = interp(xs, w, xb)
c1 = -(wb - wa) / (xb - xa)
c2 = -wa - c1 * xa
# sagging curvature (M>0) -> downward deflection; the integration above
# gives upward-positive for positive M, so flip the sign.
return [-(w[i] + c1 * xs[i] + c2) for i in range(n)]