Amaliyot — Qattiq Jism Mexanikasi

Python yordamida qattiq jism dinamikasini simulyatsiya qilish.

---

Loyiha: Quadcopter Dinamik Modeli

Maqsad

Quadcopterning roll, pitch, yaw dinamikasini modellashtirish va simulyatsiya qilish.

---

1-Qism: Muhitni Sozlash

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint, solve_ivp
from mpl_toolkits.mplot3d import Axes3D

plt.rcParams['figure.figsize'] = (12, 6)
plt.rcParams['font.size'] = 11
plt.rcParams['axes.grid'] = True

---

2-Qism: Inertsiya Momenti Hisoblash

Oddiy Shakllar

class InertiaCalculator:
    """Turli shakllar uchun inertsiya momenti."""
    
    @staticmethod
    def disk(mass, radius):
        """Disk (silindr) markazdan."""
        return 0.5 * mass * radius**2
    
    @staticmethod
    def sphere(mass, radius):
        """To'liq shar."""
        return 0.4 * mass * radius**2
    
    @staticmethod
    def rod_center(mass, length):
        """Tayoq markazdan."""
        return (1/12) * mass * length**2
    
    @staticmethod
    def rod_end(mass, length):
        """Tayoq uchidan."""
        return (1/3) * mass * length**2
    
    @staticmethod
    def hollow_cylinder(mass, r_inner, r_outer):
        """Ichi bo'sh silindr."""
        return 0.5 * mass * (r_inner**2 + r_outer**2)
    
    @staticmethod
    def point_masses(masses, distances):
        """Nuqtaviy massalar."""
        return np.sum(np.array(masses) * np.array(distances)**2)

# Misol: Quadcopter ramkasi
motor_mass = 0.05  # kg
arm_length = 0.15  # m
num_motors = 4

I_motors = InertiaCalculator.point_masses(
    [motor_mass] * num_motors,
    [arm_length] * num_motors
)
print(f"Motorlar inertsiyasi: {I_motors:.6f} kg·m²")

# Ramka (X shakl) — ikkita tayoq
frame_mass = 0.1  # kg (har bir tayoq)
frame_length = 0.3  # m (diagonal)
I_frame = 2 * InertiaCalculator.rod_center(frame_mass/2, frame_length)
print(f"Ramka inertsiyasi: {I_frame:.6f} kg·m²")

I_total = I_motors + I_frame
print(f"Jami Izz: {I_total:.6f} kg·m²")

3D Inertsiya Tensori

def quadcopter_inertia_tensor(motor_mass, arm_length, frame_mass, frame_width):
    """
    Quadcopter inertsiya tensorini hisoblash.
    X konfiguratsiya (motorlar diagonallarda).
    """
    # Motorlar pozitsiyalari (x, y, z)
    d = arm_length / np.sqrt(2)  # x va y komponentlar
    motor_positions = np.array([
        [ d,  d, 0],
        [-d,  d, 0],
        [-d, -d, 0],
        [ d, -d, 0]
    ])
    
    # Inertsiya tensori
    I = np.zeros((3, 3))
    
    for pos in motor_positions:
        x, y, z = pos
        I[0, 0] += motor_mass * (y**2 + z**2)  # Ixx
        I[1, 1] += motor_mass * (x**2 + z**2)  # Iyy
        I[2, 2] += motor_mass * (x**2 + y**2)  # Izz
        I[0, 1] -= motor_mass * x * y
        I[0, 2] -= motor_mass * x * z
        I[1, 2] -= motor_mass * y * z
    
    # Simmetriya
    I[1, 0] = I[0, 1]
    I[2, 0] = I[0, 2]
    I[2, 1] = I[1, 2]
    
    return I

# Hisoblash
I_tensor = quadcopter_inertia_tensor(0.05, 0.15, 0.1, 0.02)
print("Inertsiya tensori (kg·m²):")
print(I_tensor)
print(f"\nIxx = {I_tensor[0,0]:.6f}")
print(f"Iyy = {I_tensor[1,1]:.6f}")
print(f"Izz = {I_tensor[2,2]:.6f}")

---

3-Qism: Aylanma Dinamika Simulyatsiyasi

Oddiy Model (1 O'q)

def single_axis_rotation(y, t, I, tau_func):
    """
    Bir o'q atrofida aylanish.
    y = [theta, omega]
    """
    theta, omega = y
    tau = tau_func(t, theta, omega)
    alpha = tau / I
    return [omega, alpha]

# Misol: Qisqa moment impulsi
def torque_impulse(t, theta, omega):
    """t=0 da qisqa moment."""
    if 0 <= t < 0.1:
        return 0.5  # N·m
    else:
        return 0

I_test = 0.01  # kg·m²
y0 = [0, 0]  # theta=0, omega=0
t_span = np.linspace(0, 2, 500)

sol = odeint(single_axis_rotation, y0, t_span, args=(I_test, torque_impulse))

fig, axs = plt.subplots(1, 3, figsize=(15, 4))

axs[0].plot(t_span, [torque_impulse(t, 0, 0) for t in t_span], 'r-', linewidth=2)
axs[0].set_xlabel('Vaqt (s)')
axs[0].set_ylabel('Moment (N·m)')
axs[0].set_title('Kiritilgan Moment')

axs[1].plot(t_span, np.degrees(sol[:, 1]), 'b-', linewidth=2)
axs[1].set_xlabel('Vaqt (s)')
axs[1].set_ylabel('Burchak tezlik (°/s)')
axs[1].set_title('Burchak Tezlik')

axs[2].plot(t_span, np.degrees(sol[:, 0]), 'g-', linewidth=2)
axs[2].set_xlabel('Vaqt (s)')
axs[2].set_ylabel('Burchak (°)')
axs[2].set_title('Burchak')

plt.tight_layout()
plt.savefig('single_axis_rotation.png', dpi=150)
plt.show()

PD Kontrollerli Burchak Boshqaruvi

def angle_control(y, t, I, target_angle, Kp, Kd):
    """PD kontroller bilan burchak boshqaruvi."""
    theta, omega = y
    error = target_angle - theta
    tau = Kp * error - Kd * omega  # PD kontroller
    tau = np.clip(tau, -1, 1)  # Moment chegarasi
    alpha = tau / I
    return [omega, alpha]

# Maqsad burchak: 45 daraja
target = np.radians(45)
Kp, Kd = 2.0, 0.5
y0 = [0, 0]
t_span = np.linspace(0, 3, 500)

sol_pd = odeint(angle_control, y0, t_span, args=(0.01, target, Kp, Kd))

plt.figure(figsize=(12, 5))
plt.subplot(1, 2, 1)
plt.plot(t_span, np.degrees(sol_pd[:, 0]), 'b-', linewidth=2)
plt.axhline(y=45, color='r', linestyle='--', label='Maqsad')
plt.xlabel('Vaqt (s)')
plt.ylabel('Burchak (°)')
plt.title('PD Boshqaruv')
plt.legend()

plt.subplot(1, 2, 2)
plt.plot(t_span, np.degrees(sol_pd[:, 1]), 'g-', linewidth=2)
plt.xlabel('Vaqt (s)')
plt.ylabel('Burchak tezlik (°/s)')
plt.title('Burchak Tezlik')

plt.tight_layout()
plt.savefig('pd_control.png', dpi=150)
plt.show()

---

4-Qism: Quadcopter 3D Dinamikasi

Euler Burchaklari bilan Model

class QuadcopterDynamics:
    """Quadcopter 3D dinamik modeli."""
    
    def __init__(self, mass=1.0, arm_length=0.15):
        self.m = mass
        self.L = arm_length
        self.g = 9.81
        
        # Inertsiya (soddalashtirilgan)
        self.Ixx = 0.01
        self.Iyy = 0.01
        self.Izz = 0.02
        
        # Propeller koeffitsientlari
        self.k_thrust = 1e-5   # Thrust koeffitsienti
        self.k_torque = 1e-6   # Reaktiv moment koeffitsienti
    
    def dynamics(self, t, state, motor_speeds):
        """
        State: [x, y, z, vx, vy, vz, phi, theta, psi, p, q, r]
        motor_speeds: [omega1, omega2, omega3, omega4] rad/s
        """
        x, y, z, vx, vy, vz, phi, theta, psi, p, q, r = state
        w1, w2, w3, w4 = motor_speeds
        
        # Thrust kuchlari
        T1 = self.k_thrust * w1**2
        T2 = self.k_thrust * w2**2
        T3 = self.k_thrust * w3**2
        T4 = self.k_thrust * w4**2
        T_total = T1 + T2 + T3 + T4
        
        # Momentlar
        tau_phi = self.L * (T2 - T4)  # Roll
        tau_theta = self.L * (T3 - T1)  # Pitch
        tau_psi = self.k_torque * (w1**2 - w2**2 + w3**2 - w4**2)  # Yaw
        
        # Aylanish matritsasi (Z-Y-X konvensiya)
        c_phi, s_phi = np.cos(phi), np.sin(phi)
        c_theta, s_theta = np.cos(theta), np.sin(theta)
        c_psi, s_psi = np.cos(psi), np.sin(psi)
        
        # Chiziqli tezlanish (global frame)
        ax = (T_total/self.m) * (c_phi*s_theta*c_psi + s_phi*s_psi)
        ay = (T_total/self.m) * (c_phi*s_theta*s_psi - s_phi*c_psi)
        az = (T_total/self.m) * c_phi*c_theta - self.g
        
        # Burchak tezlanish
        p_dot = (tau_phi - (self.Izz - self.Iyy)*q*r) / self.Ixx
        q_dot = (tau_theta - (self.Ixx - self.Izz)*p*r) / self.Iyy
        r_dot = (tau_psi - (self.Iyy - self.Ixx)*p*q) / self.Izz
        
        # Euler burchaklari tezligi
        phi_dot = p + (q*s_phi + r*c_phi)*np.tan(theta)
        theta_dot = q*c_phi - r*s_phi
        psi_dot = (q*s_phi + r*c_phi) / c_theta if abs(c_theta) > 0.01 else 0
        
        return [vx, vy, vz, ax, ay, az, phi_dot, theta_dot, psi_dot, p_dot, q_dot, r_dot]

# Simulyatsiya
quad = QuadcopterDynamics(mass=1.0, arm_length=0.15)

# Boshlang'ich holat
state0 = [0, 0, 0,  # pozitsiya
          0, 0, 0,  # tezlik
          0, 0, 0,  # burchaklar
          0, 0, 0]  # burchak tezliklar

# Motor tezliklari (hovering + roll)
hover_speed = 1000  # rad/s
motor_speeds = [hover_speed, hover_speed*1.1, hover_speed, hover_speed*0.9]

t_span = (0, 2)
t_eval = np.linspace(0, 2, 200)

# Solve
sol = solve_ivp(quad.dynamics, t_span, state0, 
                args=(motor_speeds,), t_eval=t_eval, method='RK45')

# Vizualizatsiya
fig, axs = plt.subplots(2, 2, figsize=(14, 10))

# Pozitsiya
axs[0, 0].plot(sol.t, sol.y[0], label='x')
axs[0, 0].plot(sol.t, sol.y[1], label='y')
axs[0, 0].plot(sol.t, sol.y[2], label='z')
axs[0, 0].set_xlabel('Vaqt (s)')
axs[0, 0].set_ylabel('Pozitsiya (m)')
axs[0, 0].set_title('Pozitsiya')
axs[0, 0].legend()

# Burchaklar
axs[0, 1].plot(sol.t, np.degrees(sol.y[6]), label='Roll (φ)')
axs[0, 1].plot(sol.t, np.degrees(sol.y[7]), label='Pitch (θ)')
axs[0, 1].plot(sol.t, np.degrees(sol.y[8]), label='Yaw (ψ)')
axs[0, 1].set_xlabel('Vaqt (s)')
axs[0, 1].set_ylabel('Burchak (°)')
axs[0, 1].set_title('Euler Burchaklari')
axs[0, 1].legend()

# Burchak tezliklar
axs[1, 0].plot(sol.t, np.degrees(sol.y[9]), label='p (roll rate)')
axs[1, 0].plot(sol.t, np.degrees(sol.y[10]), label='q (pitch rate)')
axs[1, 0].plot(sol.t, np.degrees(sol.y[11]), label='r (yaw rate)')
axs[1, 0].set_xlabel('Vaqt (s)')
axs[1, 0].set_ylabel('Burchak tezlik (°/s)')
axs[1, 0].set_title('Burchak Tezliklar')
axs[1, 0].legend()

# 3D traektoriya
ax3d = fig.add_subplot(2, 2, 4, projection='3d')
ax3d.plot(sol.y[0], sol.y[1], sol.y[2], 'b-', linewidth=2)
ax3d.set_xlabel('X (m)')
ax3d.set_ylabel('Y (m)')
ax3d.set_zlabel('Z (m)')
ax3d.set_title('3D Traektoriya')

plt.tight_layout()
plt.savefig('quadcopter_dynamics.png', dpi=150)
plt.show()

---

5-Qism: Dumalanish Simulyatsiyasi

Nishablikdan Dumalanish

def rolling_simulation(shape='disk', h=1.0, angle_deg=30):
    """
    Turli shakllarning nishablikdan dumalanishini simulyatsiya qilish.
    """
    g = 9.81
    angle = np.radians(angle_deg)
    
    # Inertsiya koeffitsientlari
    shape_coeffs = {
        'disk': 0.5,
        'sphere': 0.4,
        'ring': 1.0,
        'hollow_sphere': 2/3
    }
    
    c = shape_coeffs.get(shape, 0.5)
    
    # Tezlanish (sirpanmasdan)
    a = g * np.sin(angle) / (1 + c)
    
    # Nishablik uzunligi
    L = h / np.sin(angle)
    
    # Vaqt va tezlik
    t_final = np.sqrt(2*L/a)
    v_final = np.sqrt(2*a*L)
    
    # Simulyatsiya
    t = np.linspace(0, t_final, 100)
    s = 0.5 * a * t**2
    v = a * t
    
    return t, s, v, v_final, t_final

# Turli shakllarni taqqoslash
shapes = ['disk', 'sphere', 'ring', 'hollow_sphere']
labels = ['Disk', 'Shar', 'Halqa', "Bo'sh shar"]
colors = ['blue', 'green', 'red', 'orange']

fig, axs = plt.subplots(1, 2, figsize=(14, 5))

for shape, label, color in zip(shapes, labels, colors):
    t, s, v, v_f, t_f = rolling_simulation(shape, h=1.0, angle_deg=30)
    axs[0].plot(t, s, color=color, label=f'{label} (v={v_f:.2f} m/s)', linewidth=2)
    axs[1].plot(t, v, color=color, label=label, linewidth=2)

axs[0].set_xlabel('Vaqt (s)')
axs[0].set_ylabel('Masofa (m)')
axs[0].set_title('Nishablikdan Dumalanish (h=1m, 30°)')
axs[0].legend()

axs[1].set_xlabel('Vaqt (s)')
axs[1].set_ylabel('Tezlik (m/s)')
axs[1].set_title('Tezlik')
axs[1].legend()

plt.tight_layout()
plt.savefig('rolling_comparison.png', dpi=150)
plt.show()

# Natijalar
print("Yakuniy tezliklar:")
for shape, label in zip(shapes, labels):
    _, _, _, v_f, t_f = rolling_simulation(shape)
    print(f"  {label}: v = {v_f:.3f} m/s, t = {t_f:.3f} s")

---

6-Qism: Giroskop Effekti

Pretsessiya Simulyatsiyasi

def gyroscope_precession(y, t, I_spin, I_perp, omega_spin, M, g, d):
    """
    Giroskop pretsessiyasi.
    y = [theta, phi, theta_dot, phi_dot]
    theta - nutatsiya burchagi
    phi - pretsessiya burchagi
    """
    theta, phi, theta_dot, phi_dot = y
    
    # Gravitatsion moment
    tau_g = M * g * d * np.sin(theta)
    
    # Harakatlar
    L = I_spin * omega_spin  # Spin impuls momenti
    
    # Soddalashtirilgan (kichik nutatsiya)
    phi_dot_new = tau_g / (L * np.sin(theta)) if np.sin(theta) > 0.01 else 0
    theta_dot_new = 0  # Barqaror pretsessiya
    
    return [theta_dot, phi_dot_new, 0, 0]

# Parametrlar
I_spin = 0.001  # kg·m²
omega_spin = 100  # rad/s (spin tezligi)
M = 0.1  # kg
g = 9.81
d = 0.05  # m

y0 = [np.radians(45), 0, 0, 0]  # 45 daraja burchak
t_span = np.linspace(0, 5, 500)

sol = odeint(gyroscope_precession, y0, t_span, 
             args=(I_spin, I_spin, omega_spin, M, g, d))

# Vizualizatsiya
fig, axs = plt.subplots(1, 2, figsize=(14, 5))

axs[0].plot(t_span, np.degrees(sol[:, 0]), 'b-', linewidth=2)
axs[0].set_xlabel('Vaqt (s)')
axs[0].set_ylabel('Nutatsiya (°)')
axs[0].set_title('Nutatsiya Burchagi')

axs[1].plot(t_span, np.degrees(sol[:, 1]), 'r-', linewidth=2)
axs[1].set_xlabel('Vaqt (s)')
axs[1].set_ylabel('Pretsessiya (°)')
axs[1].set_title('Pretsessiya Burchagi')

plt.tight_layout()
plt.savefig('gyroscope.png', dpi=150)
plt.show()

# Pretsessiya tezligi
omega_p = (M * g * d) / (I_spin * omega_spin)
print(f"Pretsessiya tezligi: {np.degrees(omega_p):.2f} °/s")

---

7-Qism: Topshiriqlar

Topshiriq 1: Robot Qo'li

1. 2 bo'g'inli robot qo'lini modellashtiring
2. Har bir bo'g'in uchun inertsiya momentini hisoblang
3. Yuklarni ko'tarish uchun kerakli momentlarni toping

Topshiriq 2: Quadcopter Tuning

1. Yuqoridagi modelda PID kontrollerlar qo'shing
2. Roll va Pitch uchun alohida tuning qiling
3. Step response va settling time ni o'lchang

Topshiriq 3: Raketa Barqarorligi

1. Spin stabilizatsiyali raketa modelini yarating
2. Tashqi buzilishlarga javobini simulyatsiya qiling
3. Optimal spin tezligini toping

---

Xulosa

Bu amaliyotda:

• Inertsiya momentlarini hisoblashni o'rgandik
• 3D aylanma dinamikani modellashtiridik
• Quadcopter dinamikasini simulyatsiya qildik
• Dumalanish va giroskopik effektlarni ko'rdik

Keyingi qadam: Fizika — Suyuqlik va gaz mexanikasi.