Beam Design Structural Movement

Programming Language: Python

      
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from scipy.integrate import solve_ivp

# Convert constants to English units
ft_to_in = 12
slug_to_lb = 32.174

# Given parameters from calculations
m = (1663.63 + 1040.88) / slug_to_lb  # Convert force (lbf) to mass (slugs)
k = 1.85e7 / ft_to_in  # Convert stiffness from lbf/ft to lbf/in
c = float(input("Enter damping coefficient (lbf·s/in, enter 0 if negligible): "))
v = float(input("Enter wind speed (ft/s): "))
C_d = 1.2  # Drag coefficient (assumed)
A = 18  # Cross-sectional area in in²
rho = 0.002377  # Air density in slugs/ft³
F_wind = 0.5 * C_d * A * rho * (v**2)  # Wind force (lbf)
L = 109.5  # Beam length in inches
I = (1/12) * 6 * (8**3)  # Moment of inertia for I-Beam in in^4

# Calculating natural and damped frequencies
omega_n = np.sqrt(k / m)  # Natural frequency (rad/s)
zeta = c / (2 * np.sqrt(m * k))  # Damping ratio
omega_d = omega_n * np.sqrt(1 - zeta**2) if zeta < 1 else 0  # Damped frequency

# Time settings
time = np.linspace(0, 800, 1000)  # Time array (seconds)
A_disp = 2.625  # Amplitude in inches

# Define differential equation
def system(t, y):
    x, v = y
    M_wind = F_wind * L  # Wind-induced moment
    dxdt = v
    dvdt = (-c * v - k * x + M_wind / I) / m
    return [dxdt, dvdt]

# Solve using solve_ivp
y0 = [A_disp, 0]
sol = solve_ivp(system, [0, max(time)], y0, t_eval=time, method='RK45')

# Extract displacement
disp = sol.y[0]

# Plot results
plt.figure(figsize=(10, 6))
plt.plot(time / 60, disp, label="Displacement over time")
plt.xlabel("Time (minutes)")
plt.ylabel("Displacement (inches)")
plt.title("Oscillation of Refinery Stack")
plt.legend()
plt.grid()
plt.show(block=True)

# Streamline plot with user-input wind direction
angle_deg = float(input("Enter wind direction in degrees (0 to 360, step of 15): "))
angle_rad = np.radians(angle_deg)

x = np.linspace(-5, 5, 100)
y = np.linspace(-5, 5, 100)
X, Y = np.meshgrid(x, y)
U = v * np.cos(angle_rad) + 0.2 * np.sin(Y)  # Wind profile along the input direction
V = v * np.sin(angle_rad) + 0.1 * np.cos(X)

fig2, ax2 = plt.subplots(figsize=(8, 6))
ax2.streamplot(X, Y, U, V, density=1.5, color='b', linewidth=1)
ax2.add_patch(plt.Circle((0, 0), 0.5, color='gray', label="Stack 1"))
ax2.add_patch(plt.Circle((1, 0), 0.5, color='gray', label="Stack 2"))
ax2.add_patch(plt.Rectangle((-1, -0.25), 3, 0.5, color='black', alpha=0.5, label="Landing Platform"))
ax2.set_xlim([-5, 5])
ax2.set_ylim([-3, 3])
ax2.set_xlabel("")
ax2.set_ylabel("")
ax2.set_title(f"Streamlines Around Refinery Stacks (Wind {angle_deg}°)")
ax2.legend()
plt.grid()
plt.show(block=True