#!/usr/bin/env python3 """ ================================================================================ # BEGIN FILE: rigid_frame_analyzer.py # Date: June 26, 2025 # Author(s): Gary Argraves, Claude (Anthropic) # # Purpose: Structural analysis of single-bay, single-story rigid frame with # center peak-pitched roof, pinned supports, subject to wind and # gravity loads # # Usage Example: # python rigid_frame_analyzer.py # frame = RigidFrame(span=40, height=16, roof_pitch=30, frame_spacing=20) # results = frame.analyze(wind_load=25, dead_load=20, live_load=40) # # Features: # - Matrix structural analysis using stiffness method # - Wind load distribution on frame # - Gravity load analysis (dead + live) # - Member force calculations # - Deflection analysis # - Design check capabilities # - Load combination handling per building codes # ================================================================================ """ import numpy as np import math from dataclasses import dataclass from typing import Tuple, Dict, List import sys doc_main = lambda: f""" RIGID FRAME STRUCTURAL ANALYSIS PROGRAM ===================================== PURPOSE: Analyzes single-bay rigid frame structures with center-peaked roof under combined wind and gravity loading conditions. STRUCTURAL SYSTEM: - Single bay frame with pinned base supports - Center-peaked roof geometry - Fixed connections at eave and ridge - 20-foot frame spacing (tributary width) LOAD CASES ANALYZED: 1. Dead Load (gravity) 2. Live Load (gravity) 3. Wind Load (lateral) 4. Combined load cases per code requirements ANALYSIS METHOD: Matrix stiffness method with 2D frame elements STATUS: Ready for structural analysis """ @dataclass class LoadCase: """Container for different loading scenarios""" dead_load: float = 0.0 # psf on roof surface live_load: float = 0.0 # psf on roof surface wind_load: float = 0.0 # psf on vertical projection load_factor_dead: float = 1.0 load_factor_live: float = 1.0 load_factor_wind: float = 1.0 case_name: str = "Default" @dataclass class FrameGeometry: """Frame geometric properties""" span: float # ft - bay width eave_height: float # ft - column height roof_pitch: float # degrees frame_spacing: float = 20.0 # ft - spacing between frames @property def roof_rise(self) -> float: """Calculate roof rise at center""" return (self.span / 2.0) * math.tan(math.radians(self.roof_pitch)) @property def ridge_height(self) -> float: """Total height to ridge""" return self.eave_height + self.roof_rise @property def rafter_length(self) -> float: """Length of each rafter""" return (self.span / 2.0) / math.cos(math.radians(self.roof_pitch)) class RigidFrame: """ Rigid frame structural analysis using matrix stiffness method Node numbering: 1 ---- 3 ---- 2 | | | | | | 4 | 5 Element numbering: - Element 1: Node 4 to 1 (left column) - Element 2: Node 1 to 3 (left rafter) - Element 3: Node 3 to 2 (right rafter) - Element 4: Node 2 to 5 (right column) """ def __init__(self, geometry: FrameGeometry, section_props: Dict = None): """ Initialize rigid frame Args: geometry: FrameGeometry object with frame dimensions section_props: Dictionary with I, A, E for frame members """ self.geometry = geometry # Default section properties (can be overridden) default_props = { 'I_column': 1000.0, # in^4 - column moment of inertia 'A_column': 20.0, # in^2 - column cross-sectional area 'I_rafter': 800.0, # in^4 - rafter moment of inertia 'A_rafter': 18.0, # in^2 - rafter cross-sectional area 'E': 29000.0 # ksi - modulus of elasticity } self.section_props = section_props or default_props self._setup_nodes() self._setup_elements() def _setup_nodes(self): """Define node coordinates""" g = self.geometry # Convert to inches for calculations span_in = g.span * 12 eave_ht_in = g.eave_height * 12 ridge_ht_in = g.ridge_height * 12 self.nodes = { 1: (0.0, eave_ht_in), # Left eave 2: (span_in, eave_ht_in), # Right eave 3: (span_in/2, ridge_ht_in), # Ridge 4: (0.0, 0.0), # Left base 5: (span_in, 0.0) # Right base } def _setup_elements(self): """Define frame elements with properties""" self.elements = { 1: { # Left column 'nodes': (4, 1), 'I': self.section_props['I_column'], 'A': self.section_props['A_column'], 'E': self.section_props['E'] }, 2: { # Left rafter 'nodes': (1, 3), 'I': self.section_props['I_rafter'], 'A': self.section_props['A_rafter'], 'E': self.section_props['E'] }, 3: { # Right rafter 'nodes': (3, 2), 'I': self.section_props['I_rafter'], 'A': self.section_props['A_rafter'], 'E': self.section_props['E'] }, 4: { # Right column 'nodes': (2, 5), 'I': self.section_props['I_column'], 'A': self.section_props['A_column'], 'E': self.section_props['E'] } } def _element_stiffness_matrix(self, element_id: int) -> np.ndarray: """ Generate 6x6 element stiffness matrix for 2D frame element Returns: 6x6 stiffness matrix in global coordinates """ elem = self.elements[element_id] node_i, node_j = elem['nodes'] # Element properties E = elem['E'] * 1000 # Convert ksi to psi I = elem['I'] A = elem['A'] # Node coordinates xi, yi = self.nodes[node_i] xj, yj = self.nodes[node_j] # Element length and orientation L = math.sqrt((xj - xi)**2 + (yj - yi)**2) cos_theta = (xj - xi) / L sin_theta = (yj - yi) / L # Local stiffness matrix coefficients EA_L = E * A / L EI_L = E * I / L EI_L2 = E * I / (L * L) EI_L3 = E * I / (L * L * L) # Local stiffness matrix k_local = np.array([ [EA_L, 0, 0, -EA_L, 0, 0 ], [0, 12*EI_L3, 6*EI_L2, 0, -12*EI_L3, 6*EI_L2 ], [0, 6*EI_L2, 4*EI_L, 0, -6*EI_L2, 2*EI_L ], [-EA_L, 0, 0, EA_L, 0, 0 ], [0, -12*EI_L3, -6*EI_L2, 0, 12*EI_L3, -6*EI_L2 ], [0, 6*EI_L2, 2*EI_L, 0, -6*EI_L2, 4*EI_L ] ]) # Transformation matrix T = np.array([ [cos_theta, sin_theta, 0, 0, 0, 0], [-sin_theta, cos_theta, 0, 0, 0, 0], [0, 0, 1, 0, 0, 0], [0, 0, 0, cos_theta, sin_theta, 0], [0, 0, 0, -sin_theta, cos_theta, 0], [0, 0, 0, 0, 0, 1] ]) # Transform to global coordinates k_global = T.T @ k_local @ T return k_global def _assemble_global_stiffness(self) -> np.ndarray: """Assemble global stiffness matrix""" # 15 DOF total (5 nodes × 3 DOF each) K_global = np.zeros((15, 15)) # DOF mapping for each node (u, v, theta) dof_map = { 1: [0, 1, 2], # Node 1: DOF 0,1,2 2: [3, 4, 5], # Node 2: DOF 3,4,5 3: [6, 7, 8], # Node 3: DOF 6,7,8 4: [9, 10, 11], # Node 4: DOF 9,10,11 5: [12, 13, 14] # Node 5: DOF 12,13,14 } # Assemble each element for elem_id in self.elements: k_elem = self._element_stiffness_matrix(elem_id) node_i, node_j = self.elements[elem_id]['nodes'] # Get DOF indices for this element dof_i = dof_map[node_i] dof_j = dof_map[node_j] element_dofs = dof_i + dof_j # Add element stiffness to global matrix for i in range(6): for j in range(6): K_global[element_dofs[i], element_dofs[j]] += k_elem[i, j] return K_global def _apply_boundary_conditions(self, K: np.ndarray) -> Tuple[np.ndarray, List[int]]: """ Apply pinned support boundary conditions Returns: Reduced stiffness matrix and list of free DOFs """ # Pinned supports at nodes 4 and 5 (restrain u, v but allow rotation) constrained_dofs = [9, 10, 12, 13] # u4, v4, u5, v5 all_dofs = list(range(15)) free_dofs = [dof for dof in all_dofs if dof not in constrained_dofs] # Extract reduced stiffness matrix K_reduced = K[np.ix_(free_dofs, free_dofs)] return K_reduced, free_dofs def _calculate_loads(self, load_case: LoadCase) -> np.ndarray: """ Calculate applied loads for given load case Returns: Load vector for free DOFs """ # Initialize full load vector (15 DOF) F_full = np.zeros(15) # Frame spacing for tributary loading spacing = self.geometry.frame_spacing # Gravity loads on rafters (dead + live) if load_case.dead_load > 0 or load_case.live_load > 0: # Total gravity load per unit length of rafter w_dead = load_case.dead_load * load_case.load_factor_dead * spacing / 1000 # kips/ft w_live = load_case.live_load * load_case.load_factor_live * spacing / 1000 # kips/ft w_gravity = w_dead + w_live # Convert to component loads on inclined rafters rafter_len = self.geometry.rafter_length w_rafter = w_gravity * math.cos(math.radians(self.geometry.roof_pitch)) # Distributed load on each rafter converted to equivalent nodal loads # For uniformly loaded beam: F_end = wL/2, M_end = wL²/12 F_rafter = w_rafter * rafter_len / 2 # kips M_rafter = w_rafter * (rafter_len * 12)**2 / 12 / 12 # kip-in # Apply vertical components to nodes cos_pitch = math.cos(math.radians(self.geometry.roof_pitch)) sin_pitch = math.sin(math.radians(self.geometry.roof_pitch)) # Left rafter loads (element 2: nodes 1-3) F_full[1] += -F_rafter * cos_pitch # -Fy at node 1 F_full[6] += -F_rafter * cos_pitch # -Fy at node 3 (left side) F_full[2] += M_rafter # Mz at node 1 F_full[8] += -M_rafter # -Mz at node 3 # Right rafter loads (element 3: nodes 3-2) F_full[4] += -F_rafter * cos_pitch # -Fy at node 2 F_full[7] += -F_rafter * cos_pitch # -Fy at node 3 (right side) F_full[5] += -M_rafter # -Mz at node 2 F_full[8] += M_rafter # Mz at node 3 # Wind loads on frame if load_case.wind_load > 0: # Wind pressure on vertical projection w_wind = load_case.wind_load * load_case.load_factor_wind * spacing / 1000 # kips/ft # Distributed wind load on columns col_height = self.geometry.eave_height F_wind_col = w_wind * col_height / 2 # kips per column # Apply horizontal wind loads F_full[0] += F_wind_col # +Fx at node 1 (left eave) F_full[3] += F_wind_col # +Fx at node 2 (right eave) # Apply boundary conditions to get load vector for free DOFs _, free_dofs = self._apply_boundary_conditions(np.eye(15)) F_reduced = F_full[free_dofs] return F_reduced def analyze(self, load_case: LoadCase) -> Dict: """ Perform structural analysis for given load case Returns: Dictionary containing analysis results """ print(f"INFO: Analyzing load case: {load_case.case_name}") try: # Assemble global stiffness matrix K_global = self._assemble_global_stiffness() # Apply boundary conditions K_reduced, free_dofs = self._apply_boundary_conditions(K_global) # Calculate load vector F_reduced = self._calculate_loads(load_case) # Solve for displacements try: displacements_reduced = np.linalg.solve(K_reduced, F_reduced) except np.linalg.LinAlgError: print("ERROR: Singular stiffness matrix - check model stability") return {} # Reconstruct full displacement vector displacements_full = np.zeros(15) for i, dof in enumerate(free_dofs): displacements_full[dof] = displacements_reduced[i] # Calculate reactions reactions = K_global @ displacements_full # Extract nodal displacements nodal_disps = { 1: {'u': displacements_full[0], 'v': displacements_full[1], 'theta': displacements_full[2]}, 2: {'u': displacements_full[3], 'v': displacements_full[4], 'theta': displacements_full[5]}, 3: {'u': displacements_full[6], 'v': displacements_full[7], 'theta': displacements_full[8]}, 4: {'u': 0.0, 'v': 0.0, 'theta': displacements_full[11]}, 5: {'u': 0.0, 'v': 0.0, 'theta': displacements_full[14]} } # Calculate member forces member_forces = self._calculate_member_forces(displacements_full) # Calculate maximum deflections max_deflections = self._calculate_max_deflections(nodal_disps) results = { 'displacements': nodal_disps, 'reactions': { 'R4x': reactions[9]/1000, # kips 'R4y': reactions[10]/1000, # kips 'R5x': reactions[12]/1000, # kips 'R5y': reactions[13]/1000 # kips }, 'member_forces': member_forces, 'max_deflections': max_deflections, 'load_case': load_case.case_name } print("SUCCESS: Analysis completed successfully") return results except Exception as e: print(f"ERROR: Analysis failed - {str(e)}") return {} def _calculate_member_forces(self, displacements: np.ndarray) -> Dict: """Calculate internal member forces""" member_forces = {} # DOF mapping dof_map = {1: [0,1,2], 2: [3,4,5], 3: [6,7,8], 4: [9,10,11], 5: [12,13,14]} for elem_id in self.elements: elem = self.elements[elem_id] node_i, node_j = elem['nodes'] # Get element displacements dof_i = dof_map[node_i] dof_j = dof_map[node_j] elem_disps = np.concatenate([displacements[dof_i], displacements[dof_j]]) # Get element stiffness matrix k_elem = self._element_stiffness_matrix(elem_id) # Calculate element forces forces = k_elem @ elem_disps member_forces[f'Element_{elem_id}'] = { 'axial_i': forces[0]/1000, # kips 'shear_i': forces[1]/1000, # kips 'moment_i': forces[2]/12000, # kip-ft 'axial_j': -forces[3]/1000, # kips 'shear_j': -forces[4]/1000, # kips 'moment_j': forces[5]/12000 # kip-ft } return member_forces def _calculate_max_deflections(self, nodal_disps: Dict) -> Dict: """Calculate maximum deflections for design checks""" # Find maximum horizontal and vertical deflections max_horiz = max(abs(nodal_disps[node]['u']) for node in nodal_disps) max_vert = max(abs(nodal_disps[node]['v']) for node in nodal_disps) # Convert to more readable units max_horiz_in = max_horiz max_vert_in = max_vert # Calculate span-to-deflection ratios span_in = self.geometry.span * 12 height_in = self.geometry.eave_height * 12 horiz_ratio = span_in / max_horiz_in if max_horiz_in > 0 else float('inf') vert_ratio = span_in / max_vert_in if max_vert_in > 0 else float('inf') return { 'max_horizontal_in': max_horiz_in, 'max_vertical_in': max_vert_in, 'horizontal_span_ratio': f"L/{horiz_ratio:.0f}", 'vertical_span_ratio': f"L/{vert_ratio:.0f}" } def print_results(results: Dict, geometry: FrameGeometry): """Print formatted analysis results""" if not results: print("ERROR: No results to display") return print(f"\n{'='*60}") print(f"RIGID FRAME ANALYSIS RESULTS - {results['load_case']}") print(f"{'='*60}") print(f"\nFRAME GEOMETRY:") print(f" Span: {geometry.span:.1f} ft") print(f" Eave Height: {geometry.eave_height:.1f} ft") print(f" Roof Pitch: {geometry.roof_pitch:.1f}°") print(f" Ridge Height: {geometry.ridge_height:.1f} ft") print(f" Frame Spacing: {geometry.frame_spacing:.1f} ft") print(f"\nMAXIMUM DEFLECTIONS:") max_def = results['max_deflections'] print(f" Horizontal: {max_def['max_horizontal_in']:.3f} in ({max_def['horizontal_span_ratio']})") print(f" Vertical: {max_def['max_vertical_in']:.3f} in ({max_def['vertical_span_ratio']})") print(f"\nSUPPORT REACTIONS:") reactions = results['reactions'] print(f" Left Support (Node 4):") print(f" Horizontal: {reactions['R4x']:8.2f} kips") print(f" Vertical: {reactions['R4y']:8.2f} kips") print(f" Right Support (Node 5):") print(f" Horizontal: {reactions['R5x']:8.2f} kips") print(f" Vertical: {reactions['R5y']:8.2f} kips") print(f"\nMEMBER FORCES:") forces = results['member_forces'] for member, vals in forces.items(): print(f" {member}:") print(f" Axial (i-end): {vals['axial_i']:8.2f} kips") print(f" Shear (i-end): {vals['shear_i']:8.2f} kips") print(f" Moment (i-end): {vals['moment_i']:8.2f} kip-ft") def main(): """Main execution function""" # Display program information print(doc_main()) # Example frame geometry geometry = FrameGeometry( span=40.0, # 40 ft span eave_height=16.0, # 16 ft eave height roof_pitch=30.0, # 30 degree roof pitch frame_spacing=20.0 # 20 ft frame spacing ) # Create frame instance frame = RigidFrame(geometry) # Define load cases load_cases = [ LoadCase( dead_load=20.0, # 20 psf dead load live_load=40.0, # 40 psf live load wind_load=0.0, # No wind case_name="Dead + Live Load" ), LoadCase( dead_load=0.0, # No gravity loads live_load=0.0, wind_load=25.0, # 25 psf wind load case_name="Wind Load Only" ), LoadCase( dead_load=20.0, # Combined loading live_load=25.0, # Reduced live load per code wind_load=25.0, # Wind load load_factor_dead=1.2, load_factor_live=1.6, load_factor_wind=1.6, case_name="Ultimate Load Combination" ) ] # Analyze each load case all_results = [] for load_case in load_cases: print(f"\nNEXT: Analyzing {load_case.case_name}") results = frame.analyze(load_case) if results: all_results.append(results) print_results(results, geometry) print(f"\nSUCCESS: Analysis complete for {len(all_results)} load cases") return all_results if __name__ == "__main__": try: results = main() print(f"\nINFO: Program completed successfully") except KeyboardInterrupt: print(f"\nWARNING: Program interrupted by user") sys.exit(1) except Exception as e: print(f"ERROR: Program failed - {str(e)}") sys.exit(1) """ ================================================================================ # END FILE: rigid_frame_analyzer.py ================================================================================ """