#068 - Structuring Engineering Calculations in Python | 01 - The Basics

Writing clear, maintainable, and easily reviewable Python engineering calculation scripts.

Greetings everybody, I hope you are well.

We are staring down the barrel of a federal election here in Canada. The political landscape is fluid, who knows what’s coming.

Politics is outside of my scope for this article, instead we will look at some of the fundamental concepts involved in structuring your calculations in Python.

Engineers typically rely on tools like Excel, SMath Studio or Mathcad for calculations. They often suffice for one-off tasks. Yet, revisiting a complex spreadsheet months later frequently involves deciphering cryptic cell references, recalling implicit assumptions, or tracing the origins of specific factors, an inefficient and potentially risky process.

Inheriting somebody elses complex spreadsheet is like inheriting a house where the plumbing runs through the toaster, which is also, somehow, the primary structural support for the second floor, and changing a cell risks spontaneously changing the primary language of all household appliances to Sumerian.

As project complexity and information volume increase, relying solely on memory or poorly documented spreadsheets becomes untenable. At least it is for my information addled brain.

Obscured logic invites error. Programmatic thinking, specifically using a tool like Python, offers a robust alternative by enforcing structure, clarity, and repeatability, even for routine calculations.

The goal is making calculations explicit and reviewable. I hear plenty of pushback from engineers, usually more senior ones, who do not like the idea of Python or ‘black box’ engineering. The intent here is to provide the opposite, full transparency.

I will admit there’s subjectivity and personal preference herein. Everyone has an opinion on how things should be done, which is totally fine. This is simply my version of it. There are many paths to the same destination.

One other caveat is the code blocks below, Substack refuses to implement code syntax highlighting. A travesty, yes I know. So please do not be deterred by the blocks of code below. You can also review the code here in this GitHub gist so that your eyes don’t explode.

We will discuss an example, an anchor shear capacity check based on ACI 318-19, Building Code Requirements for Structural Concrete – in Python, focusing on why this structured approach is recommended. This is based on the US design code for reinforced concrete but the same general principles apply to any engineering calculations.

There are many things inherent in this example which should be part of a broader discussion, topics like:

• How to handle units

• How much commentary to include

• How to connect to other data sources (spreadsheets, FE models etc.)

• How to create reports

• The limits of parametrization, how much is enough?

I will discuss these topics in the future.

The Blueprint: Structuring Calculations for Clarity

“If I had an hour to solve a problem, I’d spend 55 minutes thinking about the problem and 5 minutes thinking about solutions.” - Albert Einstein

A well-structured calculation script follows a deliberate pattern designed for transparency and maintainability:

1. Imports: Declare necessary libraries upfront. What you’ll be using for your calculation. For me it’s often numpy, scipy, pandas, polars, matplotlib, these cover the majority of standard calculations but even for simple calculations relying only on the built-in math module like this example, the principle holds.

2. Fundamental Inputs: Define all basic parameters, geometry, material properties, loads, and critical assumptions in one location. Resist the overwhelming urge to use pre-calculated values that hide their origins.

3. Geometry Calculations (Derived Inputs): Explicitly calculate intermediate geometric properties from the fundamental inputs. This is often a major source of confusion.

4. Core Engineering Calculations: Implement the design code logic step-by-step, using your fundamental or previously derived values.

5. Helper Functions (Optional): Encapsulate reusable code blocks for consistency. Ideally you won’t need to touch these once they are defined with a parametric philosophy.

6. Results: Format and present the outputs clearly and concisely.

This logical flow ensures a traceable path from basic inputs to final conclusions.

Programmatic thinking is a fascinating topic, I write about it in more detail here:

Section 1: Imports - Loading Your Tools

We begin by importing required libraries. For basic math, it's simple:

import math # Provides math constants (pi) and functions (sqrt)

• Purpose: Explicitly loading the math module makes its functions available. Placing imports first clearly declares the script's dependencies.

Section 2: Fundamental Input Parameters & Assumptions - The Control Panel

This section serves as the calculation's single source of truth for inputs. Define only basic, understandable parameters here.

# ---- FUNDAMENTAL INPUT PARAMETERS & ASSUMPTIONS ----
# Basic geometry, material properties, loads, factors per ACI 318-19

# Anchor Properties
da = 0.75           # Anchor diameter (in)
futa = 58           # Anchor steel ultimate tensile strength (ksi)
hef = 12            # Effective embedment depth (in)

# Concrete Properties
fc = 5              # Concrete compressive strength (ksi)
lambda_a = 1.0      # Modification for lightweight concrete (17.2.6) - 1.0 for normalweight

# Anchor Group Geometry (Rectangular Group Assumed)
Nx = 2              # Number of anchors in group along X-axis (direction of ca1/ca2)
Ny = 2              # Number of anchors in group along Y-axis (direction of ca1_p/ca2_p)
sx = 6.0            # Spacing between anchors along X-axis (in)
sy = 6.0            # Spacing between anchors along Y-axis (in)

# Edge Distances (Defined relative to the CL of outermost anchors)
ca1 = 17.5          # Max edge distance, X-dir (in) - Direction // to Shear Vu
ca2 = 14.5          # Min edge distance, X-dir (in) - Opposite ca1
ca1_p = 20.5        # Max edge distance, Y-dir (in) - Direction perp. to Shear Vu
ca2_p = 17.5        # Min edge distance, Y-dir (in) - Opposite ca1_p

# Applied Load
V_app = 7.96        # Applied factored shear load (kips) - Assumed acting along X-axis

# Design Assumptions & Factors (ACI 318-19 References)
phi_vs = 0.70       # Strength reduction factor, steel shear (Table 21.2.1, Sec 17.3.3)
phi_vc = 0.70       # Strength reduction factor, concrete breakout/pryout (Table 21.2.1, Sec 17.3.3)
cracked = True      # Assumed condition: True or False (Affects psi_c factors)

• Rationale: Isolating inputs here provides clarity and parametrization. Changes (e.g., different fc or anchor layout Nx, sx) are made only in this section, preventing errors deep within the calculation logic. Descriptive variable names and comments (including units and code references like ACI 17.2.6) are essential for future comprehension and review. This section defines the problem unambiguously.

Section 3: Geometry Calculations - Deriving Projected Areas

Crucially, we calculate derived geometric inputs explicitly, rather than treating them as opaque, pre-calculated values.

# ---- GEOMETRY CALCULATIONS (DERIVED INPUTS) ----
# Calculate projected areas based on fundamental geometry inputs
# Assumes rectangular group & edge distances to CL of outermost anchors
# NOTE: ACI projected area calculations can be complex. This is a simplified interpretation based on 1.5*hef boundaries. VERIFY FOR REAL DESIGNS.

# Width of failure prism in X-direction (influenced by ca1, ca2)
group_width_x = (Nx - 1) * sx if Nx > 1 else 0
width_boundary_x1 = min(1.5 * hef, ca1) # Boundary distance from outer anchor CL towards edge ca1
width_boundary_x2 = min(1.5 * hef, ca2) # Boundary distance from outer anchor CL towards edge ca2
projected_width_for_AVc = group_width_x + width_boundary_x1 + width_boundary_x2

# Height of failure prism in Y-direction (influenced by ca1_p, ca2_p)
group_width_y = (Ny - 1) * sy if Ny > 1 else 0
width_boundary_y1 = min(1.5 * hef, ca1_p)
width_boundary_y2 = min(1.5 * hef, ca2_p)
projected_height_for_AVc = group_width_y + width_boundary_y1 + width_boundary_y2

AVc = projected_width_for_AVc * projected_height_for_AVc # Projected area // to load direction (in^2)

# AVc_p: projected area resisting shear perpendicular to edge (hypothetical load in Y-dir)
# Uses similar logic but requires careful check of ACI intent for perp. load case area definition
projected_width_for_AVcp = group_width_x + width_boundary_x1 + width_boundary_x2
projected_height_for_AVcp = group_width_y + width_boundary_y1 + width_boundary_y2
AVc_p = projected_width_for_AVcp * projected_height_for_AVcp # Projected area perp to load direction (in^2)

# Projected Area for Tension Breakout (ANc) - Used for Pryout calc Vcp = kcp*Ncbg
c_min_overall = min(ca1, ca2, ca1_p, ca2_p)
# Failure width/height bounded by min edge distances or 1.5*hef from outer anchors
width_anc = group_width_x + 2 * min(1.5 * hef, ca1, ca2)
height_anc = group_width_y + 2 * min(1.5 * hef, ca1_p, ca2_p)
ANc = width_anc * height_anc # Projected area for tension breakout (in^2)

# Single anchor projected areas (calculated for use in ACI group formulas)
AVco = 4.5 * ca1**2     # Area single anchor // load direction (Fig R17.7.2.1(b))
AVco_p = 4.5 * ca1_p**2 # Area single anchor perp load direction
ANco = 9 * hef**2       # Area single anchor tension (Fig R17.6.2.1(b))

# Print calculated areas for verification during script development/review
print("--- Calculated Geometric Areas ---")
print(f"{'Projected Area // Load (AVc):':<35s} = {AVc:>7.2f} in^2")
print(f"{'Projected Area perp Load (AVc_p):':<35s} = {AVc_p:>7.2f} in^2")
print(f"{'Projected Area Tension (ANc):':<35s} = {ANc:>7.2f} in^2")
print("-" * 45)

• Rationale: This section bridges the gap between fundamental inputs and the parameters needed by ACI formulas. By calculating AVc, ANc, etc., here, the script becomes self-contained and the logic transparent. Stated assumptions (rectangular group, simplified boundary logic) are vital. While ACI details require careful study for production use, this illustrates the principle of deriving inputs explicitly and it’s the same for all design codes. You need to understand the logic and flow of your particular calculation.

Section 4: Core Engineering Calculations (ACI 318-19)

This section implements the design code formulas, referencing the inputs and derived geometric properties.

# ---- CORE ENGINEERING CALCULATIONS (ACI 318-19 Ch 17) ----
# Organized by failure mode, using previously defined or calculated variables

# --- Steel Strength (ACI 17.7.1) ---
N_anchors = Nx * Ny
Ase = math.pi * da**2 / 4
Vsa_single = 0.8 * Ase * futa # Nominal steel shear strength per anchor (kips) - Eq (17.7.1.2b) simplified
Vsa = N_anchors * Vsa_single # Group nominal steel strength (kips) - Assumes direct summation; verify ACI group effects
phiVsa = phi_vs * Vsa     # Design steel shear strength (kips)

# --- Concrete Breakout Strength (Direction of Load) (ACI 17.7.2) ---
le = min(hef, 8 * da)     # Effective length for breakout (in) - 17.7.2.3
Vb1 = 7 * (le / da)**0.2 * math.sqrt(da) * lambda_a * math.sqrt(fc * 1000) * ca1**1.5 / 1000 # Eq (17.7.2.2a)
Vb2 = 9 * lambda_a * math.sqrt(fc * 1000) * ca1**1.5 / 1000 # Eq (17.7.2.2b)
Vb = min(Vb1, Vb2)        # Basic concrete breakout strength (kips) - 17.7.2.2

psi_ed = 0.7 + 0.3 * ca2 / (1.5 * ca1) if ca2 < 1.5 * ca1 else 1.0 # Edge effect (psi_ed,V - 17.7.2.5a)
psi_ec = 1.0 # Eccentricity (psi_ec,V - 17.7.2.4): Assumed 1.0
psi_c = 1.0 if cracked and lambda_a == 1.0 else (1.25 if cracked and lambda_a < 1.0 else 1.4) # Cracked concrete (psi_c,V - 17.7.2.7) - Check Table carefully
psi_h = 1.0 # Member thickness (psi_h,V - 17.7.2.8): Assumed 1.0

Vcbg = (AVc / AVco) * psi_ec * psi_ed * psi_c * psi_h * Vb # Group breakout strength (kips) - Eq (17.7.2.1b)
phiVcbg = phi_vc * Vcbg   # Design group breakout strength (kips)

# --- Concrete Breakout Strength (Perpendicular to Load) (ACI 17.7.2) ---
Vb1_p = 7 * (le / da)**0.2 * math.sqrt(da) * lambda_a * math.sqrt(fc * 1000) * ca1_p**1.5 / 1000
Vb2_p = 9 * lambda_a * math.sqrt(fc * 1000) * ca1_p**1.5 / 1000
Vb_p = min(Vb1_p, Vb2_p)
psi_ed_p = 0.7 + 0.3 * ca2_p / (1.5 * ca1_p) if ca2_p < 1.5 * ca1_p else 1.0 # Edge effect factor based on ca2_p
# Factor of 2 applies under specific geometric conditions (ACI Fig R17.7.2.1(c)) - Verify applicability!
Vcbg_p = 2 * (AVc_p / AVco_p) * psi_ec * psi_ed_p * psi_c * psi_h * Vb_p # Factor of 2 applied
phiVcbg_p = phi_vc * Vcbg_p

# --- Concrete Pryout Strength (ACI 17.7.3) ---
Nb_single = 16 * lambda_a * math.sqrt(fc * 1000) * hef**(5/3) / 1000 # Basic tension breakout single anchor (kips) - Eq (17.6.2.2a)

psi_ec_N = 1.0 # Eccentricity factor tension (Assume concentric) - 17.6.2.5
psi_ed_N = 0.7 + 0.3 * c_min_overall / (1.5 * hef) if c_min_overall < 1.5 * hef else 1.0 # Edge factor tension (Eq 17.6.2.6a)
psi_c_N = 1.0 if cracked else 1.25 # Cracking factor tension (psi_c,N - 17.6.2.7)
psi_cp_N = 1.0 # Splitting factor tension (psi_cp,N - 17.6.2.8)

Ncbg = (ANc / ANco) * psi_ec_N * psi_ed_N * psi_c_N * psi_cp_N * Nb_single # Group tension breakout (kips) - Eq (17.6.2.1b)

kcp = 1.0 if hef < 2.5 else 2.0 # Pryout coefficient - 17.7.3.1
Vcpg = kcp * Ncbg           # Pryout strength (kips) - Eq (17.7.3.1b)
phiVcpg = phi_vc * Vcpg     # Design pryout strength (kips)

• Rationale: This section translates the engineering code into executable logic. Using intermediate variables (le, Vb, Nb_single) and relying on previously defined inputs or calculated geometries (AVc, ANc, etc.) makes the formulas readable and traceable. Conditional logic (if/else structure) directly implements code provisions based on inputs like cracked.

Section 5: Helper Functions - Promoting Consistency

Defining functions for repetitive tasks like output formatting improves code quality.

# ---- HELPER FUNCTIONS ----
def format_check(label, capacity_value, applied_load):
    """Compares capacity to load and returns a formatted status string."""
    status = "PASS" if capacity_value >= applied_load else "FAIL"
    # f-string provides clear control over alignment and precision
    return f"{label:<55s} = {capacity_value:>7.2f} kips | Status: {status}"

• Rationale: Encapsulating the PASS/FAIL check and string formatting logic within a function adheres to the DRY (Don't Repeat Yourself) principle. It ensures consistent output and makes future changes to the format easy (modify only the function).

Section 6: Results - Clear Communication

The final step is presenting the results effectively. This typically involves summarizing the key outputs in a readable format.

# ---- RESULTS ----
# Use the helper function for concise and uniform output

print(f"\n{'--- Design Check Summary ---':^75}") # Centered title
print(f"{'Applied Factored Shear Load (Vu):':<55s} = {V_app:>7.2f} kips")
print("-" * 75)

print(format_check("Design Steel Shear Strength (ΦVsa)", phiVsa, V_app))
print(format_check("Design Concrete Breakout (ΦVcbg - // load)", phiVcbg, V_app))
print(format_check("Design Concrete Breakout (ΦVcbg - perp load)", phiVcbg_p, V_app))
print(format_check("Design Concrete Pryout Strength (ΦVcpg)", phiVcpg, V_app))
print("-" * 75)

# Determine controlling failure mode and overall status
min_capacity = min(phiVsa, phiVcbg, phiVcbg_p, phiVcpg)
overall_status = "PASS" if min_capacity >= V_app else "FAIL"
print(f"\n{'Overall Shear Check:':<55s}   Minimum Capacity = {min_capacity:.2f} kips → {overall_status}")

# (Note: The full script includes an additional, more detailed printout section for verification)

• Rationale: This section focuses purely on presenting the key outcomes concisely. Using the format_check helper function keeps this block clean and ensures consistency. Formatted strings (f-strings) provide control over alignment and precision, making the summary easy to read and interpret, while the final check explicitly identifies the minimum capacity and overall status.

This final printout is a key part of the validation process. You can see below, it’s pretty easy to read and understand. You can add/remove detail as your objective dictates.

========================= FINAL VERIFICATION PRINTOUT =========================

--- Fundamental Inputs & Assumptions ---
Anchor Diameter (da):               = 0.75 in
Anchor Steel Strength (futa):       = 58 ksi
Embedment Depth (hef):              = 12 in
Concrete Strength (fc):             = 5 ksi
Lightweight Mod (lambda_a):         = 1.0
Number Anchors X (Nx):              = 2
Number Anchors Y (Ny):              = 2
Spacing X (sx):                     = 6.0 in
Spacing Y (sy):                     = 6.0 in
Edge Dist Max X (ca1):              = 17.5 in
Edge Dist Min X (ca2):              = 14.5 in
Edge Dist Max Y (ca1_p):            = 20.5 in
Edge Dist Min Y (ca2_p):            = 17.5 in
Applied Shear (V_app):              = 7.96 kips
Phi Factor Steel (phi_vs):          = 0.7
Phi Factor Concrete (phi_vc):       = 0.7
Cracked Concrete Assumed:           = True

--- Derived Geometric Inputs ---
Projected Area // Load (AVc):       = 1577.00 in^2
Projected Area perp Load (AVc_p):   = 1577.00 in^2
Projected Area Tension (ANc):       = 1435.00 in^2
Single Anchor Area // (AVco):       = 1378.12 in^2
Single Anchor Area perp (AVco_p):   = 1891.12 in^2
Single Anchor Area Tension (ANco):  = 1296.00 in^2
Min Overall Edge Dist (c_min):      = 14.5 in

--- Intermediate Calculation Results ---
Number of Anchors (N_anchors):      = 4
Anchor Steel Area (Ase):            = 0.442 in^2
Effective Length (le):              = 6.00 in
Basic Breakout Strength (Vb):       = 46.59 kips
Basic Breakout Strength Perp (Vb_p): = 59.07 kips
Basic Tension Breakout Single (Nb_s): = 71.16 kips
Group Tension Breakout (Ncbg):      = 74.20 kips

--- Final Nominal Capacities ---
Nominal Steel Shear (Vsa):          = 82.00 kips
Nominal Breakout // (Vcbg):         = 46.15 kips
Nominal Breakout perp (Vcbg_p):     = 85.78 kips
Nominal Pryout (Vcpg):              = 148.39 kips

--- Final Design Capacities & Status ---
Design Steel Shear (ΦVsa):          = 57.40 kips
Design Breakout // (ΦVcbg):         = 32.31 kips
Design Breakout perp (ΦVcbg_p):     = 60.05 kips
Design Pryout (ΦVcpg):              = 103.88 kips
Applied Load (Vu):                  = 7.96 kips
Minimum Design Capacity:            = 32.31 kips
Overall Status:                     = PASS

========================= END VERIFICATION PRINTOUT ===========================

By explicitly outputting all fundamental inputs, key derived values (like the calculated geometric areas), intermediate results, and final capacities together, the script provides a complete audit trail within its own output.

This allows anyone unfamiliar with the Python code itself to quickly verify the inputs used, check the key calculated values, and confirm the final results against the summary. It demystifies the process and offers a level of transparent self-documentation that is often difficult to achieve reliably in complex spreadsheets where intermediate values can be hidden or hard to trace.

Why Structure Matters

This deliberate, structured approach provides significant advantages over ad-hoc spreadsheets or monolithic scripts:

• Transparency: The logic flows clearly from fundamental inputs to final results. Assumptions are explicit.

• Reviewability: Easy for peers (or your future self) to understand, verify, and trust the calculation.

• Maintainability: Changes to inputs, assumptions, or specific code clauses are localized and easier to implement correctly.

• Reduced Risk: Clarity minimizes the potential for hidden errors or misinterpretations.

Conclusion: Programmatic Thinking as Engineering Practice

Adopting a structured, programmatic approach for calculations goes beyond immediate efficiency; it fundamentally improves engineering rigor. By thinking systematically about the process, prioritizing clarity, explicitness, and a logical flow from inputs to outputs, we create tools that are not only more reliable and trustworthy today but vastly more usable in the long run.

Well-structured, clearly documented calculations are easier to review, debug, and modify months or years later. They become valuable assets, facilitating knowledge transfer to colleagues and ensuring your future self can quickly understand and leverage past work, rather than starting from scratch. This mindset helps manage the inherent complexity of engineering design more effectively.

Just look at the structure of Pynite, the FE library, a great example of this iterative, logical development framework.

Start simple, focus on structure, demand clarity from your tools, and build confidence in your results. It's an investment that pays dividends in reliability and enduring usability.

On that note, applying these principles of clarity and structure is becoming even more critical as we integrate more advanced computational tools. I'm currently involved in several exciting AI-based projects and initiatives internally here at Knight Piesold, exploring how we can leverage these technologies in practical engineering contexts. I look forward to sharing some insights from that work soon.

See you in the next one.

James 🌊

---

(Note: The complete, consolidated Python script suitable for saving or sharing via a Gist is here, containing all the code blocks presented above in sequence.)