材料力学Pythonコード


```
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')

```

```

```