import numpy as np
import matplotlib.pyplot as plt
# Constants
E = 200e9 # Young's modulus (Pa)
sigma_y = 250e6 # Yield stress (Pa)
epsilon_max = 0.01 # Maximum strain
# Generate strain values
epsilon = np.linspace(0, epsilon_max, 100)
# Elastic deformation (Hooke's Law)
sigma_elastic = E * epsilon
# Plastic deformation (Von Mises yield criterion)
# For simplicity, we assume that plastic deformation starts when stress reaches yield stress
sigma_plastic = np.where(sigma_elastic <= sigma_y, sigma_elastic, sigma_y)
# Plotting
plt.figure(figsize=(10, 6))
# Elastic region
plt.plot(epsilon, sigma_elastic, label='Elastic Deformation (Hooke\'s Law)', color='blue')
# Plastic region
plt.plot(epsilon, sigma_plastic, label='Plastic Deformation (Yield Condition)', color='red', linestyle='--')
plt.xlabel('Strain (ε)')
plt.ylabel('Stress (σ)')
plt.title('Elastic and Plastic Deformation')
plt.legend()
plt.grid(True)
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# Constants for strength values (example values in MPa)
tensile_strength = 500 # Tensile Strength (MPa)
compressive_strength = 600 # Compressive Strength (MPa)
shear_strength = 300 # Shear Strength (MPa)
# Constants for hardness tests (example values)
vickers_hardness = 200 # Vickers Hardness (HV)
rockwell_hardness = 60 # Rockwell Hardness (HR)
# Plotting strength values
strengths = {
'Tensile Strength': tensile_strength,
'Compressive Strength': compressive_strength,
'Shear Strength': shear_strength
}
hardness_tests = {
'Vickers Hardness': vickers_hardness,
'Rockwell Hardness': rockwell_hardness
}
plt.figure(figsize=(12, 6))
# Plot strengths
plt.subplot(1, 2, 1)
plt.bar(strengths.keys(), strengths.values(), color=['blue', 'green', 'orange'])
plt.xlabel('Strength Type')
plt.ylabel('Strength (MPa)')
plt.title('Strength Properties')
# Plot hardness
plt.subplot(1, 2, 2)
plt.bar(hardness_tests.keys(), hardness_tests.values(), color=['purple', 'red'])
plt.xlabel('Hardness Test')
plt.ylabel('Hardness Value')
plt.title('Hardness Tests')
plt.tight_layout()
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# Constants
E = 200e9 # Young's modulus (Pa)
sigma_y = 250e6 # Yield stress (Pa)
epsilon_max = 0.01 # Maximum strain
# Generate strain values
epsilon = np.linspace(0, epsilon_max, 100)
# Elastic deformation (Hooke's Law)
sigma_elastic = E * epsilon
# Plastic deformation (Yield condition)
sigma_plastic = np.where(sigma_elastic <= sigma_y, sigma_elastic, sigma_y)
# Plotting stress-strain curve
plt.figure(figsize=(12, 6))
# Elastic region
plt.plot(epsilon, sigma_elastic, label='Elastic Deformation (Hooke\'s Law)', color='blue')
# Plastic region
plt.plot(epsilon, sigma_plastic, label='Plastic Deformation (Yield Condition)', color='red', linestyle='--')
# Yield Criteria Lines
plt.axhline(y=sigma_y, color='green', linestyle=':', label='Yield Stress (σ_y)')
plt.axhline(y=sigma_y*0.5, color='purple', linestyle='-.', label='Yield Stress / 2 (σ_y/2)')
# Add labels and title
plt.xlabel('Strain (ε)')
plt.ylabel('Stress (σ)')
plt.title('Stress-Strain Curve with Yield Criteria')
plt.legend()
plt.grid(True)
# Show plot
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# Constants
E = 200e9 # Young's modulus (Pa)
sigma_y = 250e6 # Yield stress (Pa)
A = 1e-4 # Cross-sectional area (m^2)
L = 1.0 # Length of beam (m)
V = 1000 # Shear force (N)
d = 0.5 # Distance from support (m)
F = 500 # Applied force (N)
# Generate strain values for tension/compression
epsilon = np.linspace(0, 0.01, 100)
sigma_elastic = E * epsilon
sigma_plastic = np.where(sigma_elastic <= sigma_y, sigma_elastic, sigma_y)
# Shear stress calculation
shear_stress = V / A
# Moment calculation
moment = F * d
# Plotting stress-strain curve
plt.figure(figsize=(12, 8))
# Plot Elastic and Plastic Regions
plt.subplot(2, 2, 1)
plt.plot(epsilon, sigma_elastic, label='Elastic Deformation (Hooke\'s Law)', color='blue')
plt.plot(epsilon, sigma_plastic, label='Plastic Deformation (Yield Condition)', color='red', linestyle='--')
plt.xlabel('Strain (ε)')
plt.ylabel('Stress (σ)')
plt.title('Stress-Strain Curve')
plt.legend()
plt.grid(True)
# Plot Shear Stress
plt.subplot(2, 2, 2)
plt.bar(['Shear Stress'], [shear_stress], color='green')
plt.ylabel('Shear Stress (Pa)')
plt.title('Shear Stress')
# Plot Moment
plt.subplot(2, 2, 3)
plt.bar(['Moment'], [moment], color='purple')
plt.ylabel('Moment (Nm)')
plt.title('Bending Moment')
plt.tight_layout()
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# Constants
E = 200e9 # Young's modulus (Pa)
I = 1e-6 # Moment of inertia (m^4)
L = 5.0 # Length of the beam (m)
P = 1000 # Point load at the center (N)
# Calculate deflection
deflection = lambda x: (P * x * (L - x) ** 2) / (48 * E * I) # Simplified formula for midpoint deflection
# Generate x values
x = np.linspace(0, L, 100)
y = deflection(x)
# Plotting
plt.figure(figsize=(10, 6))
# Plot deflection
plt.plot(x, y, label='Beam Deflection', color='blue')
plt.xlabel('Position along the beam (x)')
plt.ylabel('Deflection (y)')
plt.title('Beam Deflection under Central Point Load')
plt.legend()
plt.grid(True)
# Show plot
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# Constants
# Node coordinates (for a simple truss with 3 nodes and 3 members)
nodes = np.array([
[0, 0], # Node 1 (x1, y1)
[4, 0], # Node 2 (x2, y2)
[2, 3] # Node 3 (x3, y3)
])
# Member connections (nodes connected by members)
members = [
(0, 1), # Member 1 (Node 1 to Node 2)
(1, 2), # Member 2 (Node 2 to Node 3)
(2, 0) # Member 3 (Node 3 to Node 1)
]
# External forces (assuming vertical load at Node 3)
external_forces = np.array([
[0, 0], # Node 1 (no external force)
[0, 0], # Node 2 (no external force)
[0, -10] # Node 3 (vertical load of 10 N)
])
# Calculate member lengths and angles
def calculate_member_properties(nodes, members):
lengths = []
angles = []
for (i, j) in members:
dx = nodes[j, 0] - nodes[i, 0]
dy = nodes[j, 1] - nodes[i, 1]
length = np.sqrt(dx**2 + dy**2)
angle = np.arctan2(dy, dx)
lengths.append(length)
angles.append(angle)
return lengths, angles
lengths, angles = calculate_member_properties(nodes, members)
# Solve for member forces using Method of Joints
# For simplicity, assume all members are in tension and solve for forces
def method_of_joints_forces(lengths, angles, external_forces):
num_members = len(members)
forces = np.zeros(num_members)
# Simple analysis assuming no internal loads at joints, forces are adjusted manually
for i in range(num_members):
forces[i] = external_forces[2, 1] / num_members # Just a placeholder calculation
return forces
forces = method_of_joints_forces(lengths, angles, external_forces)
# Plotting
plt.figure(figsize=(8, 6))
# Plot nodes
for i, (x, y) in enumerate(nodes):
plt.plot(x, y, 'o', label=f'Node {i+1}')
# Plot members
for (i, j) in members:
plt.plot([nodes[i, 0], nodes[j, 0]], [nodes[i, 1], nodes[j, 1]], 'k-')
plt.xlabel('X (m)')
plt.ylabel('Y (m)')
plt.title('Truss Structure')
plt.legend()
plt.grid(True)
# Display member forces
for idx, force in enumerate(forces):
print(f"Member {idx+1} Force: {force:.2f} N")
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# Constants
sigma_max = 500e6 # Maximum stress (Pa)
sigma_min = 100e6 # Minimum stress (Pa)
N_f = 1e6 # Number of cycles to failure
# Generate stress and number of cycles data for S-N curve
cycles = np.logspace(3, 7, 100) # Logarithmic scale for cycles
stress = sigma_max * (N_f / cycles)**0.1 # Simplified S-N curve
# Plot S-N curve
plt.figure(figsize=(12, 6))
plt.subplot(1, 2, 1)
plt.loglog(cycles, stress, label='S-N Curve', color='blue')
plt.xlabel('Number of Cycles to Failure (N)')
plt.ylabel('Stress Amplitude (σ)')
plt.title('S-N Curve')
plt.grid(True)
plt.legend()
# Constants for J-Integral and Stress Intensity Factor
applied_stress = 100e6 # Applied stress (Pa)
crack_length = 0.01 # Crack length (m)
# Stress Intensity Factor
K = applied_stress * np.sqrt(np.pi * crack_length)
# J-Integral (simplified calculation for demonstration)
J = (applied_stress**2 * crack_length) / (2 * np.pi)
# Plot J-Integral and Stress Intensity Factor
plt.subplot(1, 2, 2)
plt.bar(['Stress Intensity Factor (K)', 'J-Integral (J)'], [K, J], color=['purple', 'green'])
plt.ylabel('Value')
plt.title('Fracture Mechanics Parameters')
plt.grid(True)
# Display plots
plt.tight_layout()
plt.show()
# Print calculated values
print(f"Stress Intensity Factor (K): {K:.2e} Pa*m^0.5")
print(f"J-Integral (J): {J:.2e} J/m^2")
import numpy as np
import matplotlib.pyplot as plt
# Constants for magnetic materials
magnetic_fields = np.linspace(0, 1.5, 100) # Magnetic field strength (T)
magnetic_permeability = 4 * np.pi * 1e-7 # Permeability of free space (H/m)
# Magnetic materials properties
# Assuming soft magnetic material with varying permeability
permeabilities = [1.0, 50.0, 200.0] # Relative permeabilities
for mu_r in permeabilities:
B = mu_r * magnetic_permeability * magnetic_fields # Magnetic flux density (T)
plt.plot(magnetic_fields, B, label=f'μ_r = {mu_r}')
# Plotting magnetic properties
plt.figure(figsize=(12, 6))
plt.plot(magnetic_fields, B, label='Magnetic Flux Density (B) vs. Magnetic Field Strength (H)')
plt.xlabel('Magnetic Field Strength (H) [T]')
plt.ylabel('Magnetic Flux Density (B) [T]')
plt.title('Magnetic Properties of Materials')
plt.legend()
plt.grid(True)
# Constants for insulation materials
temperature_range = np.linspace(-50, 150, 100) # Temperature range (°C)
thermal_conductivity = 0.2 # Thermal conductivity of insulation material (W/m·K)
thermal_resistance = 1 / (thermal_conductivity * np.ones_like(temperature_range)) # Simplified
# Plotting insulation materials
plt.figure(figsize=(12, 6))
plt.plot(temperature_range, thermal_resistance, label='Thermal Resistance vs. Temperature', color='red')
plt.xlabel('Temperature (°C)')
plt.ylabel('Thermal Resistance (K/W)')
plt.title('Thermal Properties of Insulation Materials')
plt.legend()
plt.grid(True)
# Display plots
plt.tight_layout()
plt.show()
import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import solve
# Define the truss
# Nodes and connections (simple example)
nodes = np.array([
[0, 0], # Node 1
[4, 0], # Node 2
[2, 3] # Node 3
])
members = [
(0, 1), # Member 1 (Node 1 to Node 2)
(1, 2), # Member 2 (Node 2 to Node 3)
(2, 0) # Member 3 (Node 3 to Node 1)
]
# Define loads and support reactions
# For simplicity, assume a vertical load at Node 3
loads = np.array([
[0, 0], # Node 1
[0, 0], # Node 2
[0, -10] # Node 3 (Load of 10 N)
])
# Method of Joints: Calculate forces in each member
def method_of_joints(nodes, members, loads):
num_nodes = len(nodes)
forces = np.zeros(len(members))
# Set up the system of equations for equilibrium at joints
A = np.zeros((2 * num_nodes, len(members)))
b = np.zeros(2 * num_nodes)
# Populate matrices
for i, (n1, n2) in enumerate(members):
dx = nodes[n2, 0] - nodes[n1, 0]
dy = nodes[n2, 1] - nodes[n1, 1]
length = np.sqrt(dx**2 + dy**2)
A[2*n1:2*n1+2, i] = [dx/length, dy/length]
A[2*n2:2*n2+2, i] = [-dx/length, -dy/length]
# Apply loads
b[2:2*num_nodes:2] = loads[:, 1]
# Solve the system of equations
forces = solve(A, b)
return forces
forces = method_of_joints(nodes, members, loads)
# Plotting
plt.figure(figsize=(12, 6))
# Plot nodes and members
for i, (x, y) in enumerate(nodes):
plt.plot(x, y, 'o', label=f'Node {i+1}')
for (i, j) in members:
plt.plot([nodes[i, 0], nodes[j, 0]], [nodes[i, 1], nodes[j, 1]], 'k-')
# Plot forces
for idx, force in enumerate(forces):
plt.text((nodes[members[idx][0], 0] + nodes[members[idx][1], 0]) / 2,
(nodes[members[idx][0], 1] + nodes[members[idx][1], 1]) / 2,
f'F{idx+1}={force:.2f} N', fontsize=12, color='red')
plt.xlabel('X (m)')
plt.ylabel('Y (m)')
plt.title('Truss Structure and Member Forces')
plt.legend()
plt.grid(True)
plt.show()
# Print member forces
for idx, force in enumerate(forces):
print(f"Member {idx+1} Force: {force:.2f} N")
import numpy as np
import matplotlib.pyplot as plt
# Define truss nodes and members
nodes = np.array([
[0, 0], # Node 1
[4, 0], # Node 2
[2, 3] # Node 3
])
members = [
(0, 1), # Member 1 (Node 1 to Node 2)
(1, 2), # Member 2 (Node 2 to Node 3)
(2, 0) # Member 3 (Node 3 to Node 1)
]
# Applied loads
loads = np.array([
[0, 0], # Node 1
[0, 0], # Node 2
[0, -10] # Node 3 (Load of 10 N)
])
# Cross-sectional area of each member (example values)
areas = np.array([1e-4, 1e-4, 1e-4]) # m^2
# Method of Joints: Calculate forces in each member
def method_of_joints(nodes, members, loads):
num_nodes = len(nodes)
num_members = len(members)
forces = np.zeros(num_members)
# Set up system of equations for equilibrium at joints
A = np.zeros((2 * num_nodes, num_members))
b = np.zeros(2 * num_nodes)
# Populate matrices
for i, (n1, n2) in enumerate(members):
dx = nodes[n2, 0] - nodes[n1, 0]
dy = nodes[n2, 1] - nodes[n1, 1]
length = np.sqrt(dx**2 + dy**2)
A[2*n1:2*n1+2, i] = [dx/length, dy/length]
A[2*n2:2*n2+2, i] = [-dx/length, -dy/length]
# Apply loads
b[2:2*num_nodes:2] = loads[:, 1]
# Solve system of equations
forces = np.linalg.solve(A, b)
return forces
forces = method_of_joints(nodes, members, loads)
# Calculate stress in each member
stresses = forces / areas
# Deformation (assuming small deformation and linear elasticity)
def calculate_deformation(forces, areas, lengths, E):
deformation = forces / (E * areas)
return deformation
# Young's modulus (example value for steel)
E = 210e9 # Pa
# Calculate member lengths
def calculate_lengths(nodes, members):
lengths = []
for (i, j) in members:
dx = nodes[j, 0] - nodes[i, 0]
dy = nodes[j, 1] - nodes[i, 1]
length = np.sqrt(dx**2 + dy**2)
lengths.append(length)
return np.array(lengths)
lengths = calculate_lengths(nodes, members)
deformations = calculate_deformation(forces, areas, lengths, E)
# Plotting
plt.figure(figsize=(12, 6))
# Plot nodes and members
for i, (x, y) in enumerate(nodes):
plt.plot(x, y, 'o', label=f'Node {i+1}')
for (i, j) in members:
plt.plot([nodes[i, 0], nodes[j, 0]], [nodes[i, 1], nodes[j, 1]], 'k-')
# Plot member stresses
for idx, stress in enumerate(stresses):
plt.text((nodes[members[idx][0], 0] + nodes[members[idx][1], 0]) / 2,
(nodes[members[idx][0], 1] + nodes[members[idx][1], 1]) / 2,
f'S{idx+1}={stress:.2e} Pa', fontsize=12, color='blue')
plt.xlabel('X (m)')
plt.ylabel('Y (m)')
plt.title('Truss Structure with Member Stresses')
plt.legend()
plt.grid(True)
plt.show()
# Print member stresses and deformations
for idx, (stress, deformation) in enumerate(zip(stresses, deformations)):
print(f"Member {idx+1} Stress: {stress:.2e} Pa")
print(f"Member {idx+1} Deformation: {deformation:.2e} m")
import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import eigh
# Beam properties
L = 10.0 # Length of the beam (m)
EI = 1e7 # Flexural rigidity (N·m^2)
# Discretize the beam into elements
num_elements = 10
num_nodes = num_elements + 1
dx = L / num_elements
nodes = np.linspace(0, L, num_nodes)
# Initialize global stiffness matrix and mass matrix
K = np.zeros((num_nodes, num_nodes))
M = np.zeros((num_nodes, num_nodes))
# Assemble the stiffness and mass matrices
for i in range(num_elements):
k_local = (EI / dx**3) * np.array([
[12, 6*dx, -12, 6*dx],
[6*dx, 4*dx**2, -6*dx, 2*dx**2],
[-12, -6*dx, 12, -6*dx],
[6*dx, 2*dx**2, -6*dx, 4*dx**2]
])
m_local = (dx / 420) * np.array([
[156, 22*dx, 54, -13*dx],
[22*dx, 4*dx**2, 13*dx, -3*dx**2],
[54, 13*dx, 156, -22*dx],
[-13*dx, -3*dx**2, -22*dx, 4*dx**2]
])
K[i:i+4, i:i+4] += k_local
M[i:i+4, i:i+4] += m_local
# Apply boundary conditions (fixed at both ends)
K = K[1:-1, 1:-1]
M = M[1:-1, 1:-1]
# Solve for eigenvalues and eigenvectors
eigenvalues, eigenvectors = eigh(K, M)
# Calculate natural frequencies (in Hz)
natural_frequencies = np.sqrt(eigenvalues) / (2 * np.pi)
# Plot mode shapes
plt.figure(figsize=(12, 6))
for i in range(3):
plt.plot(nodes[1:-1], eigenvectors[:, i], label=f'Mode {i+1} - {natural_frequencies[i]:.2f} Hz')
plt.xlabel('Position along the beam (m)')
plt.ylabel('Mode Shape')
plt.title('Mode Shapes of a Beam')
plt.legend()
plt.grid(True)
plt.show()
# Print natural frequencies
for i, freq in enumerate(natural_frequencies):
print(f"Natural Frequency {i+1}: {freq:.2f} Hz")
import numpy as np
import matplotlib.pyplot as plt
# Stress components
sigma_x = 100 # Normal stress in x direction (MPa)
sigma_y = 50 # Normal stress in y direction (MPa)
tau_xy = 30 # Shear stress (MPa)
# Calculate principal stresses
sigma_avg = (sigma_x + sigma_y) / 2
radius = np.sqrt(((sigma_x - sigma_y) / 2) ** 2 + tau_xy ** 2)
sigma_1 = sigma_avg + radius
sigma_2 = sigma_avg - radius
# Calculate shear stress on inclined plane
theta = np.linspace(0, 180, 180) # Degrees
theta_rad = np.radians(theta)
sigma_n = sigma_avg + ((sigma_x - sigma_y) / 2) * np.cos(2 * theta_rad) + tau_xy * np.sin(2 * theta_rad)
tau_n = -((sigma_x - sigma_y) / 2) * np.sin(2 * theta_rad) + tau_xy * np.cos(2 * theta_rad)
# Mohr's Circle
center = sigma_avg
radius_mohr = radius
# Plot Mohr's Circle
fig, ax = plt.subplots(1, 2, figsize=(14, 6))
# Plot stress components and principal stresses
ax[0].plot([sigma_x, sigma_y], [tau_xy, tau_xy], 'ro', label='Stress Components')
ax[0].axhline(0, color='k', linestyle='--')
ax[0].axvline(0, color='k', linestyle='--')
ax[0].plot([sigma_1, sigma_2], [0, 0], 'bo', label='Principal Stresses')
ax[0].set_xlabel('Normal Stress (MPa)')
ax[0].set_ylabel('Shear Stress (MPa)')
ax[0].set_title('Stress Components and Principal Stresses')
ax[0].legend()
ax[0].grid(True)
# Mohr's Circle
circle = plt.Circle((center, 0), radius_mohr, color='b', fill=False)
ax[1].add_artist(circle)
ax[1].plot(center, 0, 'ro', label='Center')
ax[1].plot([sigma_1], [0], 'go', label='Principal Stress σ1')
ax[1].plot([sigma_2], [0], 'mo', label='Principal Stress σ2')
ax[1].set_xlim(center - 1.5 * radius_mohr, center + 1.5 * radius_mohr)
ax[1].set_ylim(-1.5 * radius_mohr, 1.5 * radius_mohr)
ax[1].set_xlabel('Normal Stress (MPa)')
ax[1].set_ylabel('Shear Stress (MPa)')
ax[1].set_title('Mohr\'s Circle')
ax[1].legend()
ax[1].grid(True)
ax[1].set_aspect('equal')
plt.tight_layout()
plt.show()
# Print principal stresses
print(f"Principal Stress σ1: {sigma_1:.2f} MPa")
print(f"Principal Stress σ2: {sigma_2:.2f} MPa")
print(f"Maximum Shear Stress: {radius:.2f} MPa")
import numpy as np
import matplotlib.pyplot as plt
# Beam properties
L = 10.0 # Length of the beam (m)
w = 1000 # Uniformly distributed load (N/m)
P = 5000 # Point load (N)
a = 3.0 # Position of the point load from the left support (m)
# Define beam length and load positions
x = np.linspace(0, L, 1000)
# Calculate reaction forces
RA = (P * (L - a) + w * L**2 / 2) / L
RB = (P * a + w * L**2 / 2) / L
# Shear Force Calculation
def shear_force(x, P, a, RA, RB):
return np.where(x < a, RA - w * x, RA - w * x - P)
# Bending Moment Calculation
def bending_moment(x, P, a, RA, RB):
return np.where(x < a, RA * x - 0.5 * w * x**2, RA * x - 0.5 * w * x**2 - P * (x - a))
# Plot Shear Force and Bending Moment
plt.figure(figsize=(14, 6))
# Plot Shear Force
plt.subplot(1, 2, 1)
plt.plot(x, shear_force(x, P, a, RA, RB), 'b-', label='Shear Force')
plt.axhline(0, color='k', linestyle='--')
plt.xlabel('Position along the beam (m)')
plt.ylabel('Shear Force (N)')
plt.title('Shear Force Diagram')
plt.legend()
plt.grid(True)
# Plot Bending Moment
plt.subplot(1, 2, 2)
plt.plot(x, bending_moment(x, P, a, RA, RB), 'r-', label='Bending Moment')
plt.axhline(0, color='k', linestyle='--')
plt.xlabel('Position along the beam (m)')
plt.ylabel('Bending Moment (Nm)')
plt.title('Bending Moment Diagram')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# Beam properties
L = 10.0 # Length of the beam (m)
w = 5000 # Uniformly distributed load (N/m)
E = 200e9 # Modulus of Elasticity (Pa)
I = 1e-4 # Second moment of area (m^4)
b = 0.1 # Beam width (m)
h = 0.2 # Beam height (m)
# Calculate bending stress
def bending_stress(M, y, I):
return (M * y) / I
# Calculate shear stress
def shear_stress(V, Q, I, t):
return (V * Q) / (I * t)
# Calculate deflection
def deflection(x, w, L, E, I):
return (w * x**2 / (24 * E * I)) * (6 * L**2 - 4 * L * x + x**2)
# Beam dimensions and loading
x = np.linspace(0, L, 1000)
M = (w * L / 2) * (1 - (x / L)) # Bending moment
y = h / 2 # Distance from neutral axis for maximum bending stress
Q = (b * h**2 / 8) # First moment of area
V = w * (L - x) # Shear force
# Calculate stresses and deflection
bending_stresses = bending_stress(M, y, I)
shear_stresses = shear_stress(V, Q, I, h)
beam_deflection = deflection(x, w, L, E, I)
# Plot Bending Stress, Shear Stress, and Deflection
plt.figure(figsize=(18, 6))
# Plot Bending Stress
plt.subplot(1, 3, 1)
plt.plot(x, bending_stresses / 1e6, 'b-', label='Bending Stress')
plt.axhline(0, color='k', linestyle='--')
plt.xlabel('Position along the beam (m)')
plt.ylabel('Bending Stress (MPa)')
plt.title('Bending Stress Distribution')
plt.legend()
plt.grid(True)
# Plot Shear Stress
plt.subplot(1, 3, 2)
plt.plot(x, shear_stresses / 1e6, 'r-', label='Shear Stress')
plt.axhline(0, color='k', linestyle='--')
plt.xlabel('Position along the beam (m)')
plt.ylabel('Shear Stress (MPa)')
plt.title('Shear Stress Distribution')
plt.legend()
plt.grid(True)
# Plot Deflection
plt.subplot(1, 3, 3)
plt.plot(x, beam_deflection * 1e3, 'g-', label='Deflection')
plt.axhline(0, color='k', linestyle='--')
plt.xlabel('Position along the beam (m)')
plt.ylabel('Deflection (mm)')
plt.title('Beam Deflection')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# Beam properties
L = 10.0 # Length of the beam (m)
w = 5000 # Uniform load (N/m)
E = 200e9 # Modulus of Elasticity (Pa)
I = 1e-4 # Second moment of area (m^4)
# Define beam length and position
x = np.linspace(0, L, 1000)
# Calculate Bending Moment (for uniform load)
M = (w * L * x) - (w * x**2 / 2)
# Calculate Bending Strain Energy
def bending_strain_energy(M, E, I, x):
return 0.5 * np.trapz(M**2 / (E * I), x)
# Calculate Strain Energy
strain_energy = bending_strain_energy(M, E, I, x)
# Apply Maxwell's Theorem (Example: Theoretical)
# Work done by forces (example values)
F1 = 1000 # Force 1 (N)
x1 = 3.0 # Position of force 1 (m)
F2 = 1500 # Force 2 (N)
x2 = 7.0 # Position of force 2 (m)
# Displacement Calculations (Example Values)
# Simple calculation for illustration
displacement = (F1 * (L - x1) + F2 * (L - x2)) / E
# Plot Bending Strain Energy
plt.figure(figsize=(14, 6))
plt.subplot(1, 2, 1)
plt.plot(x, M, 'b-', label='Bending Moment')
plt.axhline(0, color='k', linestyle='--')
plt.xlabel('Position along the beam (m)')
plt.ylabel('Bending Moment (Nm)')
plt.title('Bending Moment Distribution')
plt.legend()
plt.grid(True)
plt.subplot(1, 2, 2)
plt.bar(['Strain Energy'], [strain_energy], color='g')
plt.xlabel('Energy Type')
plt.ylabel('Bending Strain Energy (J)')
plt.title('Bending Strain Energy')
plt.grid(True)
plt.tight_layout()
plt.show()
# Print Results
print(f"Bending Strain Energy: {strain_energy:.2e} J")
print(f"Displacement due to forces (Example): {displacement:.2e} m")
import numpy as np
import matplotlib.pyplot as plt
# Parameters
r = 0.05 # Radius of the circular shaft (m)
T = 1000 # Torque applied (Nm)
L = 2.0 # Length of the shaft (m)
G = 80e9 # Shear modulus (Pa)
d = 0.01 # Wire diameter of the spring (m)
D = 0.1 # Mean diameter of the spring (m)
n = 10 # Number of active coils
# Polar moment of inertia for circular shaft
J = np.pi * r**4 / 2
# Shear stress
tau = T * r / J
# Angle of twist
theta = T * L / (G * J)
# Strain energy in torsion
U_torsion = (T**2 * L) / (2 * G * J)
# Spring constant
k = (G * d**4) / (8 * D**3 * n)
# Deflection of the spring
F = 500 # Load applied (N)
delta = F / k
# Plot Torsion Stress and Strain Energy
plt.figure(figsize=(14, 6))
# Plot Shear Stress
plt.subplot(1, 2, 1)
plt.bar(['Shear Stress'], [tau / 1e6], color='b')
plt.xlabel('Parameter')
plt.ylabel('Shear Stress (MPa)')
plt.title('Shear Stress in Circular Shaft')
plt.grid(True)
# Plot Strain Energy
plt.subplot(1, 2, 2)
plt.bar(['Strain Energy'], [U_torsion], color='g')
plt.xlabel('Energy Type')
plt.ylabel('Strain Energy (J)')
plt.title('Strain Energy in Torsion')
plt.grid(True)
plt.tight_layout()
plt.show()
# Print Results
print(f"Shear Stress in Circular Shaft: {tau:.2e} Pa")
print(f"Angle of Twist: {theta:.2e} rad")
print(f"Strain Energy in Torsion: {U_torsion:.2e} J")
print(f"Spring Constant: {k:.2e} N/m")
print(f"Deflection of Spring: {delta:.2e} m")
import numpy as np
import matplotlib.pyplot as plt
# Column properties
E = 200e9 # Modulus of Elasticity (Pa)
I = 1e-4 # Moment of Inertia (m^4)
L = np.linspace(1, 10, 100) # Length of the column (m)
# Effective length factors for different end conditions
K_pinned_pinned = 1
K_fixed_fixed = 1 / np.sqrt(2)
K_fixed_pinned = np.sqrt(2)
# Euler's Critical Load calculation
P_cr_pinned_pinned = (np.pi**2 * E * I) / (K_pinned_pinned * L)**2
P_cr_fixed_fixed = (np.pi**2 * E * I) / (K_fixed_fixed * L)**2
P_cr_fixed_pinned = (np.pi**2 * E * I) / (K_fixed_pinned * L)**2
# Experimental adjustment (example: adding a reduction factor)
reduction_factor = 0.85
P_cr_experimental_pinned_pinned = reduction_factor * P_cr_pinned_pinned
P_cr_experimental_fixed_fixed = reduction_factor * P_cr_fixed_fixed
P_cr_experimental_fixed_pinned = reduction_factor * P_cr_fixed_pinned
# Plot Critical Loads
plt.figure(figsize=(14, 7))
# Plot Euler's Critical Load
plt.subplot(1, 2, 1)
plt.plot(L, P_cr_pinned_pinned / 1e6, 'b--', label='Pinned-Pinned')
plt.plot(L, P_cr_fixed_fixed / 1e6, 'r--', label='Fixed-Fixed')
plt.plot(L, P_cr_fixed_pinned / 1e6, 'g--', label='Fixed-Pinned')
plt.xlabel('Length of the Column (m)')
plt.ylabel('Critical Load (MN)')
plt.title('Euler’s Critical Load for Different End Conditions')
plt.legend()
plt.grid(True)
# Plot Experimental Critical Load
plt.subplot(1, 2, 2)
plt.plot(L, P_cr_experimental_pinned_pinned / 1e6, 'b-', label='Exp. Pinned-Pinned')
plt.plot(L, P_cr_experimental_fixed_fixed / 1e6, 'r-', label='Exp. Fixed-Fixed')
plt.plot(L, P_cr_experimental_fixed_pinned / 1e6, 'g-', label='Exp. Fixed-Pinned')
plt.xlabel('Length of the Column (m)')
plt.ylabel('Critical Load (MN)')
plt.title('Experimental Critical Load for Different End Conditions')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()
# Print Results
print(f"Euler's Critical Load for Pinned-Pinned Column (1 m): {P_cr_pinned_pinned[0]:.2e} N")
print(f"Experimental Critical Load for Pinned-Pinned Column (1 m): {P_cr_experimental_pinned_pinned[0]:.2e} N")
import numpy as np
import matplotlib.pyplot as plt
# Parameters for Thin and Thick Rotating Disks
rho = 7800 # Density of the disk material (kg/m^3)
omega = 1000 # Angular velocity (rad/s)
r = np.linspace(0, 0.5, 100) # Radial distance from the center (m)
t = 0.01 # Thickness of the thick disk (m)
# Thin Rotating Disk Stresses
sigma_r_thin = (rho * omega**2 * r**2) / 2
sigma_theta_thin = (rho * omega**2 * r**2) / 4
# Thick Rotating Disk Stresses (simplified model)
# For thick disks, assuming a simplified formula for demonstration
sigma_r_thick = (rho * omega**2 * r**2) / 3
sigma_theta_thick = (rho * omega**2 * r**2) / 6
# Plot Stresses in Thin and Thick Disks
plt.figure(figsize=(14, 7))
# Plot for Thin Rotating Disk
plt.subplot(1, 2, 1)
plt.plot(r, sigma_r_thin / 1e6, 'b-', label='Radial Stress')
plt.plot(r, sigma_theta_thin / 1e6, 'r-', label='Tangential Stress')
plt.xlabel('Radial Distance (m)')
plt.ylabel('Stress (MPa)')
plt.title('Thin Rotating Disk Stresses')
plt.legend()
plt.grid(True)
# Plot for Thick Rotating Disk
plt.subplot(1, 2, 2)
plt.plot(r, sigma_r_thick / 1e6, 'b-', label='Radial Stress')
plt.plot(r, sigma_theta_thick / 1e6, 'r-', label='Tangential Stress')
plt.xlabel('Radial Distance (m)')
plt.ylabel('Stress (MPa)')
plt.title('Thick Rotating Disk Stresses')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()
# Print Results
print(f"Thin Rotating Disk Radial Stress at r=0.25m: {sigma_r_thin[50]:.2e} Pa")
print(f"Thin Rotating Disk Tangential Stress at r=0.25m: {sigma_theta_thin[50]:.2e} Pa")
print(f"Thick Rotating Disk Radial Stress at r=0.25m: {sigma_r_thick[50]:.2e} Pa")
print(f"Thick Rotating Disk Tangential Stress at r=0.25m: {sigma_theta_thick[50]:.2e} Pa")
import numpy as np
import matplotlib.pyplot as plt
# Parameters for Plate Bending
E = 2.1e11 # Young’s Modulus (Pa)
t = 0.01 # Thickness (m)
nu = 0.3 # Poisson’s Ratio
q = 1e4 # Uniform Load (N/m^2)
D = E * t**3 / (12 * (1 - nu**2))
# Plate dimensions
x = np.linspace(-1, 1, 100)
y = np.linspace(-1, 1, 100)
X, Y = np.meshgrid(x, y)
# Bending deflection of a rectangular plate
W_plate = (q / D) * (X**4 + 2 * X**2 * Y**2 + Y**4)
# Parameters for Circular Disk Bending
R = 0.5 # Radius of the disk (m)
r = np.linspace(0, R, 100)
# Bending deflection of a circular disk
W_disk = (q * r**2 * (3 * R**2 - r**2)) / (16 * E * t * (1 - nu**2))
# Plot Plate Bending
plt.figure(figsize=(14, 7))
plt.subplot(1, 2, 1)
plt.contourf(X, Y, W_plate, cmap='viridis')
plt.colorbar(label='Deflection (m)')
plt.title('Bending of Rectangular Plate')
plt.xlabel('x (m)')
plt.ylabel('y (m)')
plt.grid(True)
# Plot Disk Bending
plt.subplot(1, 2, 2)
plt.plot(r, W_disk)
plt.xlabel('Radial Distance (m)')
plt.ylabel('Deflection (m)')
plt.title('Bending of Circular Disk')
plt.grid(True)
plt.tight_layout()
plt.show()
# Print Results
print(f"Maximum deflection of the rectangular plate: {np.max(W_plate):.2e} m")
print(f"Maximum deflection of the circular disk: {np.max(W_disk):.2e} m")
import numpy as np
import matplotlib.pyplot as plt
# Material Properties
E = 200e9 # Young's Modulus (Pa)
G = 80e9 # Shear Modulus (Pa)
A = 0.01 # Cross-sectional area (m^2)
L0 = 1.0 # Original length (m)
# Applied Forces
F = np.linspace(0, 10000, 100) # Normal force (N)
F_shear = np.linspace(0, 5000, 100) # Shear force (N)
# Normal Stress and Strain
sigma = F / A
epsilon = sigma / E
delta_L = F * L0 / (E * A)
# Shear Stress and Strain
tau = F_shear / A
gamma = tau / G
delta_x = gamma * L0
# Plotting
plt.figure(figsize=(14, 7))
# Normal Stress and Strain
plt.subplot(2, 2, 1)
plt.plot(F, sigma / 1e6, 'b-', label='Normal Stress')
plt.xlabel('Applied Force (N)')
plt.ylabel('Normal Stress (MPa)')
plt.title('Normal Stress vs Applied Force')
plt.grid(True)
plt.subplot(2, 2, 2)
plt.plot(sigma / 1e6, epsilon, 'r-', label='Normal Strain')
plt.xlabel('Normal Stress (MPa)')
plt.ylabel('Normal Strain')
plt.title('Normal Strain vs Normal Stress')
plt.grid(True)
# Shear Stress and Strain
plt.subplot(2, 2, 3)
plt.plot(F_shear, tau / 1e6, 'g-', label='Shear Stress')
plt.xlabel('Shear Force (N)')
plt.ylabel('Shear Stress (MPa)')
plt.title('Shear Stress vs Shear Force')
plt.grid(True)
plt.subplot(2, 2, 4)
plt.plot(tau / 1e6, gamma, 'm-', label='Shear Strain')
plt.xlabel('Shear Stress (MPa)')
plt.ylabel('Shear Strain')
plt.title('Shear Strain vs Shear Stress')
plt.grid(True)
plt.tight_layout()
plt.show()
# Print Results
print(f"Maximum Normal Stress (MPa): {np.max(sigma) / 1e6:.2f}")
print(f"Maximum Normal Strain: {np.max(epsilon):.5f}")
print(f"Maximum Shear Stress (MPa): {np.max(tau) / 1e6:.2f}")
print(f"Maximum Shear Strain: {np.max(gamma):.5f}")
print(f"Deformation (m) for Normal Stress: {np.max(delta_L):.5f}")
print(f"Deformation (m) for Shear Stress: {np.max(delta_x):.5f}")
import numpy as np
import matplotlib.pyplot as plt
# Material properties
E1 = 200e9 # Young's modulus (Pa) - Material 1
E2 = 150e9 # Young's modulus (Pa) - Material 2
A1 = 1e-4 # Cross-sectional area (m^2) - Material 1
A2 = 1e-4 # Cross-sectional area (m^2) - Material 2
L1 = 1.0 # Length (m) - Material 1
L2 = 1.0 # Length (m) - Material 2
F = 10000 # Tensile force (N)
# Strain calculation
def calculate_strain(F, E, A, L):
stress = F / A
strain = stress / E
return strain
strain1 = calculate_strain(F, E1, A1, L1)
strain2 = calculate_strain(F, E2, A2, L2)
# Total strain and elongation
total_strain = strain1 + strain2
total_displacement = strain1 * L1 + strain2 * L2
# Plotting
lengths = np.array([L1, L2])
strains = np.array([strain1, strain2])
displacements = np.array([strain1 * L1, strain2 * L2])
fig, ax1 = plt.subplots()
color = 'tab:blue'
ax1.set_xlabel('Component')
ax1.set_ylabel('Strain', color=color)
ax1.bar(['Material 1', 'Material 2'], strains, color=color, alpha=0.6, label='Strain')
ax1.tick_params(axis='y', labelcolor=color)
ax2 = ax1.twinx()
color = 'tab:red'
ax2.set_ylabel('Elongation (m)', color=color)
ax2.plot(['Material 1', 'Material 2'], displacements, color=color, marker='o', linestyle='--', label='Elongation')
ax2.tick_params(axis='y', labelcolor=color)
fig.tight_layout()
plt.title('Strain and Elongation of a Tensile Rod')
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# Two-Bar Framework Analysis
def two_bar_forces(P, theta1, theta2):
# Angles in radians
theta1_rad = np.radians(theta1)
theta2_rad = np.radians(theta2)
# Solve for forces using equilibrium equations
A = np.array([
[np.sin(theta1_rad), np.sin(theta2_rad)],
[np.cos(theta1_rad), np.cos(theta2_rad)]
])
B = np.array([P, 0])
forces = np.linalg.solve(A, B)
return forces
# Parameters for Two-Bar Framework
P = 1000 # External load (N)
theta1 = 30 # Angle of the first bar (degrees)
theta2 = 60 # Angle of the second bar (degrees)
forces_two_bar = two_bar_forces(P, theta1, theta2)
# Three-Bar Framework Analysis
def three_bar_forces(P, theta1, theta2, theta3):
# Angles in radians
theta1_rad = np.radians(theta1)
theta2_rad = np.radians(theta2)
theta3_rad = np.radians(theta3)
# System of equations: Ax = B
A = np.array([
[np.sin(theta1_rad), np.sin(theta2_rad), np.sin(theta3_rad)],
[np.cos(theta1_rad), np.cos(theta2_rad), np.cos(theta3_rad)],
[1, 1, 1]
])
B = np.array([P, 0, 0])
forces_three_bar = np.linalg.lstsq(A, B, rcond=None)[0]
return forces_three_bar
# Parameters for Three-Bar Framework
theta1_3 = 45 # Angle of the first bar (degrees)
theta2_3 = 60 # Angle of the second bar (degrees)
theta3_3 = 90 # Angle of the third bar (degrees)
forces_three_bar = three_bar_forces(P, theta1_3, theta2_3, theta3_3)
# Plotting Results
plt.figure(figsize=(14, 7))
# Two-Bar Framework
plt.subplot(1, 2, 1)
plt.bar(['Bar 1', 'Bar 2'], forces_two_bar)
plt.ylabel('Force (N)')
plt.title('Forces in Two-Bar Framework')
plt.grid(True)
# Three-Bar Framework
plt.subplot(1, 2, 2)
plt.bar(['Bar 1', 'Bar 2', 'Bar 3'], forces_three_bar)
plt.ylabel('Force (N)')
plt.title('Forces in Three-Bar Framework')
plt.grid(True)
plt.tight_layout()
plt.show()
# Print Results
print(f"Two-Bar Framework Forces (N):\nBar 1: {forces_two_bar[0]:.2f}\nBar 2: {forces_two_bar[1]:.2f}")
print(f"Three-Bar Framework Forces (N):\nBar 1: {forces_three_bar[0]:.2f}\nBar 2: {forces_three_bar[1]:.2f}\nBar 3: {forces_three_bar[2]:.2f}")
import numpy as np
import matplotlib.pyplot as plt
# Material Properties
alpha = 12e-6 # Coefficient of thermal expansion (1/°C)
E = 200e9 # Young's Modulus (Pa)
# Temperature Change
delta_T = np.linspace(-100, 100, 200) # Temperature change (°C)
# Thermal Strain and Stress
thermal_strain = alpha * delta_T
thermal_stress = E * thermal_strain
# Plotting
plt.figure(figsize=(12, 6))
# Thermal Strain
plt.subplot(1, 2, 1)
plt.plot(delta_T, thermal_strain, 'b-', label='Thermal Strain')
plt.xlabel('Temperature Change (°C)')
plt.ylabel('Thermal Strain')
plt.title('Thermal Strain vs Temperature Change')
plt.grid(True)
plt.legend()
# Thermal Stress
plt.subplot(1, 2, 2)
plt.plot(delta_T, thermal_stress / 1e6, 'r-', label='Thermal Stress')
plt.xlabel('Temperature Change (°C)')
plt.ylabel('Thermal Stress (MPa)')
plt.title('Thermal Stress vs Temperature Change')
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.show()
# Print Results
print(f"Maximum Thermal Strain: {np.max(thermal_strain):.6f}")
print(f"Maximum Thermal Stress (MPa): {np.max(thermal_stress) / 1e6:.2f}")
import numpy as np
import matplotlib.pyplot as plt
# Beam parameters
L = 10 # Length of the beam (m)
P = 1000 # Point load (N) at 4 meters from the left support
w = 500 # Uniformly distributed load (N/m)
# Beam sections
x = np.linspace(0, L, 1000)
# Reactions at supports
R1 = (P * (L - 4) + w * L**2 / 2) / L
R2 = (P * 4 + w * L**2 / 2) / L
# Shear Force Calculation
V = np.piecewise(x,
[x < 4, (x >= 4) & (x <= L)],
[R1 - w * x, R1 - w * x - P])
# Bending Moment Calculation
M = np.piecewise(x,
[x < 4, (x >= 4) & (x <= L)],
[R1 * x - 0.5 * w * x**2, R1 * x - P * (x - 4) - 0.5 * w * x**2 + 0.5 * w * 16])
# Plotting
fig, ax1 = plt.subplots(figsize=(10, 6))
# Shear Force Plot
color = 'tab:blue'
ax1.set_xlabel('Position along the beam (m)')
ax1.set_ylabel('Shear Force (N)', color=color)
ax1.plot(x, V, color=color, label='Shear Force')
ax1.tick_params(axis='y', labelcolor=color)
# Bending Moment Plot
ax2 = ax1.twinx()
color = 'tab:red'
ax2.set_ylabel('Bending Moment (Nm)', color=color)
ax2.plot(x, M, color=color, linestyle='--', label='Bending Moment')
ax2.tick_params(axis='y', labelcolor=color)
# Adding labels and title
fig.tight_layout()
ax1.legend(loc='upper left')
ax2.legend(loc='upper right')
plt.title('Shear Force and Bending Moment for a Simply Supported Beam')
plt.grid(True)
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# Beam Parameters
L = 10.0 # Length of the beam (m)
P = 1000.0 # Point Load (N)
w = 500.0 # Uniform Load (N/m)
# For Point Load Case
x = np.linspace(0, L, 500)
V_point_load = np.where(x <= L/2, P/2, -P/2)
M_point_load = np.where(x <= L/2, P/2 * x, P/2 * (L - x))
# For Uniform Load Case
V_uniform_load = w * L / 2 - w * x
M_uniform_load = (w * x**2) / 2 - (w * L * x) / 2
# Plotting
plt.figure(figsize=(14, 10))
# Point Load Case
plt.subplot(2, 2, 1)
plt.plot(x, V_point_load, 'b-', label='Shear Force (Point Load)')
plt.xlabel('Position along the beam (m)')
plt.ylabel('Shear Force (N)')
plt.title('Shear Force Diagram: Point Load')
plt.grid(True)
plt.legend()
plt.subplot(2, 2, 2)
plt.plot(x, M_point_load, 'r-', label='Bending Moment (Point Load)')
plt.xlabel('Position along the beam (m)')
plt.ylabel('Bending Moment (Nm)')
plt.title('Bending Moment Diagram: Point Load')
plt.grid(True)
plt.legend()
# Uniform Load Case
plt.subplot(2, 2, 3)
plt.plot(x, V_uniform_load, 'g-', label='Shear Force (Uniform Load)')
plt.xlabel('Position along the beam (m)')
plt.ylabel('Shear Force (N)')
plt.title('Shear Force Diagram: Uniform Load')
plt.grid(True)
plt.legend()
plt.subplot(2, 2, 4)
plt.plot(x, M_uniform_load, 'm-', label='Bending Moment (Uniform Load)')
plt.xlabel('Position along the beam (m)')
plt.ylabel('Bending Moment (Nm)')
plt.title('Bending Moment Diagram: Uniform Load')
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# Beam Parameters
L = 10.0 # Length of the beam (m)
# Loading Conditions
P1 = 1000.0 # Point Load 1 (N)
P2 = 500.0 # Point Load 2 (N)
w = 500.0 # Uniform Load (N/m)
# 1. Simply Supported Beam with Multiple Point Loads
x = np.linspace(0, L, 500)
P1_pos = L / 3
P2_pos = 2 * L / 3
# Reactions
R1 = (P1 + P2) / 2
R2 = (P1 + P2) / 2
# Shear Force and Bending Moment for Multiple Point Loads
V_multiple_point_loads = np.where(x <= P1_pos, R1 - P1, np.where(x <= P2_pos, R1 - P1, R1 - P1 - P2))
M_multiple_point_loads = np.where(x <= P1_pos, R1 * x - P1 * (x - P1_pos), np.where(x <= P2_pos, R1 * x - P1 * (x - P1_pos) - P2 * (x - P2_pos), R1 * x - P1 * (x - P1_pos) - P2 * (x - P2_pos)))
# 2. Simply Supported Beam with Uniformly Distributed Load
V_uniform_load = w * L / 2 - w * x
M_uniform_load = (w * x**2) / 2 - (w * L * x) / 2
# 3. Cantilever Beam with Point Load at Free End
x_cantilever = np.linspace(0, L, 500)
P = 1000.0 # Point Load at the end (N)
V_cantilever = -P * np.ones_like(x_cantilever)
M_cantilever = -P * x_cantilever
# 4. Cantilever Beam with Two Point Loads
x_cantilever_two_loads = np.linspace(0, L, 500)
P1_pos_cantilever = L / 4
P2_pos_cantilever = 3 * L / 4
# Reactions
R_cantilever = P1_pos_cantilever * P2 - P2_pos_cantilever * P1
# Shear Force and Bending Moment for Two Point Loads
V_cantilever_two_loads = np.where(x_cantilever <= P1_pos_cantilever, -P, np.where(x_cantilever <= P2_pos_cantilever, -P + P1, -P + P1 - P2))
M_cantilever_two_loads = np.where(x_cantilever <= P1_pos_cantilever, -P * x_cantilever, np.where(x_cantilever <= P2_pos_cantilever, -P * x_cantilever + P1 * (x_cantilever - P1_pos_cantilever), -P * x_cantilever + P1 * (x_cantilever - P1_pos_cantilever) - P2 * (x_cantilever - P2_pos_cantilever)))
# 5. Cantilever Beam with Uniformly Distributed Load
V_cantilever_uniform_load = -w * x_cantilever
M_cantilever_uniform_load = -w * x_cantilever**2 / 2
# Plotting
plt.figure(figsize=(14, 10))
# Multiple Point Loads on Simply Supported Beam
plt.subplot(3, 2, 1)
plt.plot(x, V_multiple_point_loads, 'b-', label='Shear Force (Multiple Point Loads)')
plt.xlabel('Position along the beam (m)')
plt.ylabel('Shear Force (N)')
plt.title('Shear Force Diagram: Multiple Point Loads')
plt.grid(True)
plt.legend()
plt.subplot(3, 2, 2)
plt.plot(x, M_multiple_point_loads, 'r-', label='Bending Moment (Multiple Point Loads)')
plt.xlabel('Position along the beam (m)')
plt.ylabel('Bending Moment (Nm)')
plt.title('Bending Moment Diagram: Multiple Point Loads')
plt.grid(True)
plt.legend()
# Uniform Load on Simply Supported Beam
plt.subplot(3, 2, 3)
plt.plot(x, V_uniform_load, 'g-', label='Shear Force (Uniform Load)')
plt.xlabel('Position along the beam (m)')
plt.ylabel('Shear Force (N)')
plt.title('Shear Force Diagram: Uniform Load')
plt.grid(True)
plt.legend()
plt.subplot(3, 2, 4)
plt.plot(x, M_uniform_load, 'm-', label='Bending Moment (Uniform Load)')
plt.xlabel('Position along the beam (m)')
plt.ylabel('Bending Moment (Nm)')
plt.title('Bending Moment Diagram: Uniform Load')
plt.grid(True)
plt.legend()
# Point Load at Free End of Cantilever Beam
plt.subplot(3, 2, 5)
plt.plot(x_cantilever, V_cantilever, 'c-', label='Shear Force (Point Load Cantilever)')
plt.xlabel('Position along the beam (m)')
plt.ylabel('Shear Force (N)')
plt.title('Shear Force Diagram: Point Load Cantilever')
plt.grid(True)
plt.legend()
plt.subplot(3, 2, 6)
plt.plot(x_cantilever, M_cantilever, 'k-', label='Bending Moment (Point Load Cantilever)')
plt.xlabel('Position along the beam (m)')
plt.ylabel('Bending Moment (Nm)')
plt.title('Bending Moment Diagram: Point Load Cantilever')
plt.grid(True)
plt.legend()
# Adjust layout and show
plt.tight_layout()
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# Beam Parameters
L = 10.0 # Length of the beam (m)
P = 1000.0 # Point Load (N)
E = 200e9 # Young's Modulus (Pa)
I = 1e-6 # Moment of Inertia (m^4)
# 1. Cantilever Beam with Point Load at Free End
x_cantilever = np.linspace(0, L, 500)
deflection_cantilever = (P * x_cantilever**3) / (3 * E * I)
# 2. Simply Supported Beam with Point Load at Center
x_simply_supported = np.linspace(0, L, 500)
deflection_simply_supported = (P * x_simply_supported * (L - x_simply_supported)) / (6 * E * I) * (L**2 - x_simply_supported**2)
# Plotting
plt.figure(figsize=(10, 6))
# Cantilever Beam Deflection
plt.plot(x_cantilever, deflection_cantilever, 'b-', label='Cantilever Beam with Point Load at Free End')
plt.xlabel('Position along the beam (m)')
plt.ylabel('Deflection (m)')
plt.title('Deflection of Cantilever Beam with Point Load at Free End')
plt.grid(True)
plt.legend()
# Simply Supported Beam Deflection
plt.plot(x_simply_supported, deflection_simply_supported, 'r-', label='Simply Supported Beam with Point Load at Center')
plt.xlabel('Position along the beam (m)')
plt.ylabel('Deflection (m)')
plt.title('Deflection of Simply Supported Beam with Point Load at Center')
plt.grid(True)
plt.legend()
# Adjust layout and show
plt.tight_layout()
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# Define material properties
E = 200e9 # Young's Modulus (Pa)
nu = 0.3 # Poisson's Ratio
# Stiffness matrix for isotropic material (Lamé parameters)
lambda_ = (E * nu) / ((1 + nu) * (1 - 2 * nu))
mu = E / (2 * (1 + nu))
# Define stress and strain tensors
def stress_tensor(x, y, z):
return np.array([
[x, y, z],
[y, x, z],
[z, z, x]
])
def strain_tensor(stress):
return np.array([
[stress[0, 0] / E, stress[0, 1] / E, stress[0, 2] / E],
[stress[1, 0] / E, stress[1, 1] / E, stress[1, 2] / E],
[stress[2, 0] / E, stress[2, 1] / E, stress[2, 2] / E]
])
# Calculate stress and strain
stress = stress_tensor(1e6, 1e6, 0)
strain = strain_tensor(stress)
# Plot stress and strain tensors
fig, ax = plt.subplots(1, 2, figsize=(14, 7))
# Stress Tensor
cax1 = ax[0].matshow(stress, cmap='viridis')
fig.colorbar(cax1, ax=ax[0])
ax[0].set_title('Stress Tensor')
ax[0].set_xlabel('x')
ax[0].set_ylabel('y')
# Strain Tensor
cax2 = ax[1].matshow(strain, cmap='viridis')
fig.colorbar(cax2, ax=ax[1])
ax[1].set_title('Strain Tensor')
ax[1].set_xlabel('x')
ax[1].set_ylabel('y')
plt.tight_layout()
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# Define material properties
E = 200e9 # Young's Modulus (Pa)
nu = 0.3 # Poisson's Ratio
# Construct the compliance matrix (inverse of the stiffness matrix)
D = np.array([
[1/E, -nu/E, -nu/E, 0, 0, 0],
[-nu/E, 1/E, -nu/E, 0, 0, 0],
[-nu/E, -nu/E, 1/E, 0, 0, 0],
[0, 0, 0, 1/(2*mu), 0, 0],
[0, 0, 0, 0, 1/(2*mu), 0],
[0, 0, 0, 0, 0, 1/(2*mu)]
])
# Define strain components
epsilon_xx = np.linspace(-0.01, 0.01, 100)
epsilon_yy = np.linspace(-0.01, 0.01, 100)
epsilon_zz = np.linspace(-0.01, 0.01, 100)
# Create meshgrid for plotting
E_xx, E_yy = np.meshgrid(epsilon_xx, epsilon_yy)
E_zz = 0.01 - E_xx - E_yy # Volumetric strain constraint
# Calculate stress components
def stress_from_strain(epsilon_xx, epsilon_yy, epsilon_zz):
strain = np.array([epsilon_xx, epsilon_yy, epsilon_zz, 0, 0, 0])
stress = np.dot(D, strain)
return stress
# Calculate stress and volumetric strain
stress = stress_from_strain(E_xx, E_yy, E_zz)
volumetric_strain = E_xx + E_yy + E_zz
# Plotting
fig, ax = plt.subplots(1, 2, figsize=(14, 7))
# Stress tensor components
cax1 = ax[0].imshow(stress[0, :, :], extent=(-0.01, 0.01, -0.01, 0.01), cmap='viridis', origin='lower')
fig.colorbar(cax1, ax=ax[0])
ax[0].set_title('Stress $\sigma_{xx}$')
ax[0].set_xlabel('ε_xx')
ax[0].set_ylabel('ε_yy')
# Volumetric Strain
cax2 = ax[1].imshow(volumetric_strain, extent=(-0.01, 0.01, -0.01, 0.01), cmap='viridis', origin='lower')
fig.colorbar(cax2, ax=ax[1])
ax[1].set_title('Volumetric Strain $\epsilon_v$')
ax[1].set_xlabel('ε_xx')
ax[1].set_ylabel('ε_yy')
plt.tight_layout()
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# Function to calculate reactions for a simply supported beam
def calculate_reactions(L, P, w):
R_A = (P / 2) + (w * L / 2) # Reaction at A
R_B = (P / 2) + (w * L / 2) # Reaction at B
return R_A, R_B
# Function to calculate shear force and bending moment for a simply supported beam
def calculate_sfd_bmd(L, P, w):
x = np.linspace(0, L, 100)
R_A, R_B = calculate_reactions(L, P, w)
# Shear force calculation
SFD = np.where(x < L/2, R_A - P - w*x, R_A - w*L + w*x)
# Bending moment calculation
BMD = np.where(x < L/2, R_A*x - P*x + 0.5*w*x**2, R_A*x - P*x + 0.5*w*L**2 - w*x*(x-L))
return x, SFD, BMD
# Function to calculate shear force and bending moment for a cantilever beam
def calculate_cantilever_sfd_bmd(L, P, w):
x = np.linspace(0, L, 100)
# Shear force calculation
SFD = np.where(x < L, -P - w*x, 0)
# Bending moment calculation
BMD = np.where(x < L, -P*(L-x) - 0.5*w*x**2, 0)
return x, SFD, BMD
# Function to plot Shear Force Diagram (SFD) and Bending Moment Diagram (BMD)
def plot_sfd_bmd(x, SFD, BMD, title):
fig, ax = plt.subplots(2, 1, figsize=(12, 10))
# Plot Shear Force Diagram
ax[0].plot(x, SFD, label='Shear Force')
ax[0].set_xlabel('Position along Beam (m)')
ax[0].set_ylabel('Shear Force (N)')
ax[0].set_title(f'{title} - Shear Force Diagram')
ax[0].grid(True)
ax[0].legend()
# Plot Bending Moment Diagram
ax[1].plot(x, BMD, label='Bending Moment', color='orange')
ax[1].set_xlabel('Position along Beam (m)')
ax[1].set_ylabel('Bending Moment (Nm)')
ax[1].set_title(f'{title} - Bending Moment Diagram')
ax[1].grid(True)
ax[1].legend()
plt.tight_layout()
plt.show()
# Parameters for the beams
L = 5 # Length of the beam in meters
P = 1000 # Concentrated load in Newtons
w = 500 # Uniformly distributed load in N/m
# For Simply Supported Beam with Combined Loads
x, SFD, BMD = calculate_sfd_bmd(L, P, w)
plot_sfd_bmd(x, SFD, BMD, 'Simply Supported Beam with Combined Loads')
# For Cantilever Beam with Combined Loads
x, SFD, BMD = calculate_cantilever_sfd_bmd(L, P, w)
plot_sfd_bmd(x, SFD, BMD, 'Cantilever Beam with Combined Loads')