To find the deformations on the quarter of a ring I found the following formulations:
ϵ=(θ′−1/R)∗t/2, θ(s)′′=−F4EI∗cosθ, θ(0)=π/2, θ(L/4)=0
Being L the medimum circumference, R the medium radius, t the thickness, E young modulus, I inertia and s the spatial coordinate following the arc.
I tried writing a code to solve this but the deformations I obtain are only negative while there should be a change in sign at some angle. Below is the code I tried:
``
import numpy as np
from scipy.integrate import solve_bvp
import matplotlib.pyplot as plt
# Costanti note
F = 10 # Forza massima applicata in Newton
E = 70000 # Modulo di Young in MPascal
R_e = 28 # Raggio esterno in millimetri
R_i = 27.5 # Raggio interno in millimetri
t = R_e - R_i # Spessore dell'anello in millimetri
h = 5 # Altezza della sezione in millimetri
I = 1/12 * h * t**3 # Momento d'inerzia della sezione in mm^4
R = (R_e + R_i) / 2 # Raggio medio in millimetri
L = 2 * np.pi * R # Circonferenza media in millimetri
# Equazioni differenziali
def equations(s, y):
theta = np.pi / 2 - s / R # Calcolo dell'angolo theta per ogni punto s
theta_prime = y[1]
dtheta_prime_ds = - (F / (4 * E * I)) * np.cos(theta)
return np.vstack((theta_prime, dtheta_prime_ds))
# Condizioni al contorno
def bc(ya, yb):
return np.array([ya[0] - np.pi / 2, yb[0]])
# Intervallo di integrazione
s = np.linspace(0, L / 4, 100)
# Condizioni iniziali stimate
y_init = np.zeros((2, s.size))
y_init[0] = np.pi / 2
# Risoluzione del problema ai valori al contorno
solution = solve_bvp(equations, bc, s, y_init)
# Verifica della soluzione
if not solution.success:
raise RuntimeError("La soluzione dell'equazione differenziale non ha avuto successo.")
# Calcolo della soluzione per questi valori di s
s_values = np.linspace(0, L / 4, 100)
theta_values = solution.sol(s_values)[0]
theta_prime_values = solution.sol(s_values)[1]
# Stampa dei valori di theta e theta'
print("theta(s):", theta_values)
print("theta'(s):", theta_prime_values)
# Calcolo della deformazione epsilon
epsilon_values = (theta_prime_values - (1 / R)) * t / 2
# Stampa dei risultati
print("s:", s_values)
print("theta(s):", theta_values)
print("theta'(s):", theta_prime_values)
print("epsilon(s):", epsilon_values)
# Visualizzazione dei risultati
plt.figure(figsize=(12, 8))
plt.subplot(3, 1, 1)
plt.plot(s_values, theta_values)
plt.xlabel('s (mm)')
plt.ylabel('theta(s)')
plt.title('Soluzione theta(s)')
plt.subplot(3, 1, 2)
plt.plot(s_values, theta_prime_values)
plt.xlabel('s (mm)')
plt.ylabel('theta\'(s)')
plt.title('Derivata theta\'(s)')
plt.subplot(3, 1, 3)
plt.plot(s_values, epsilon_values)
plt.xlabel('s (mm)')
plt.ylabel('epsilon(s)')
plt.title('Deformazione epsilon(s)')
plt.tight_layout()
plt.show()
# Grafico separato di epsilon vs theta (in gradi)
theta_degrees = np.degrees(theta_values)
plt.figure(figsize=(8, 6))
plt.plot(theta_degrees, epsilon_values)
plt.xlabel('theta (degrees)')
plt.ylabel('epsilon(s)')
plt.title('Deformazione epsilon vs theta (degrees)')
plt.grid(True)
plt.show()
``
