changgong 2023 spring

analysis1.txt

## 공학적 분석 2023-05-06 18:30

모터의 중심을 xy평면의 원점에 놓고,
모터의 각속도가 ω라고 하자.
커넥팅 로드와 모터의 연결점으로부터
모터의 중심까지의 거리를 r이라고 하자.
그러면 모터와 커넥팅로드의 연결점의
좌표는 (r cos(ωt), r sin(ωt))이다.
커넥팅 로드의 길이가 ℓ일 때(단, ℓ > r),
커넥팅로드와 마사지건 헤드의 연결점이
x축 위에서만 왕복 운동을 하므로,
그것의 좌표를 (x_head(t), 0)이라고 하면,
        x_head(t) = r cos(ωt) + √( ℓ² - (r sin(ωt))² )
이다.
헤드가 타격하는 부위의 피부를
질량(mass)가 m인 판에
감쇠계수(damping coefficient)가 c인 댐퍼와
강성(stiffness)가 k인 용수철이 달려있는
시스템으로 모델링하자.
y_cr을 피부와 접하는 그 순간에서의
모터와 커넥팅로드의 연결점의 y좌표의
절댓값이라고 하자.
마사지건 헤드가 피부와 맞닿는 시간은
        { t ∈ ℝ : |r sin(ωt)| < y_cr and cos(ωt) > 0 }
이다.
마사지건 헤드의 질량이 m_head이고,
마사지건 헤드의 가속도가 a_head이면
마사지건 헤드가 판에 가하는 힘은
        F_head(t) = if |r sin(ωt)| < y_cr and cos(ωt) > 0 then - m_head * a_head(t) else 0
이다. 여기서, a_head = x_head''이다.
F_head가 판에 가해지는 외란이므로,
우리가 알아내고자 했던 판의 지배방정식은
        m x'' + c x' + k x = F_head
이다.

analysis2.py

# 2023-05-07 14:43

import numpy as np
import matplotlib.pyplot as plt
import sympy as sp

t = sp.Symbol('t')

r_motor = sp.Symbol('r_motor')
l_motor = sp.Symbol('l_motor')
omega_motor = sp.Symbol('omega_motor')

m_body = sp.Symbol('m_body')
r_body = sp.Symbol('r_body')
J_body = sp.Symbol('J_body')
c_body = sp.Symbol('c_body')
k_body = sp.Symbol('k_body')

theta_t = sp.Function('theta_body')
theta_body = theta_t(t)

x_head = r_motor * sp.cos(omega_motor * t) + sp.sqrt(l_motor**2 - (r_motor * sp.sin(omega_motor * t)))
x_body = (1 / 2) * x_head
a_body = sp.diff(x_body, t, 2)

T_disturbance = m_body * a_body / r_body

omega_body = sp.diff(theta_body, t)
alpha_body = sp.diff(omega_body, t)

# J_body alpha_body + c_body omega_body + k_body theta_body = m_body x_body'' / r_body
eq = sp.Eq(J_body * alpha_body + c_body * omega_body + k_body * theta_body, m_body * a_body / r_body)

sol = sp.dsolve(eq, theta_body, ics={theta_t(0): 0, sp.diff(theta_body, t).subs(t, 0): 0})

print(sol)

modeling.py

## 공학적 분석 2023.05.06 21:00

import numpy as np
import matplotlib.pyplot as plt
import sympy as sp

def RK4(phi, t0, y0, tf, dt): # The Runge–Kutta method (RK4)
    ret = []
    while t0 < tf:
        k1 = phi(t0, y0)
        k2 = phi(t0 + (dt/2.0), y0 + k1 * (dt/2.0))
        k3 = phi(t0 + (dt/2.0), y0 + k2 * (dt/2.0))
        k4 = phi(t0 + dt, y0 + k3 * dt)
        y1, t1 = y0 + (dt/6.0) * (k1 + 2*k2 + 2*k3 + k4), t0 + dt
        ret = ret + [(t1, y1)]
        t0, y0 = t1, y1
    return ret

# 기타 파라메터
rpm    = 1200
y_cr   = 0.01

# 마사지 건의 사양
l      = 57.25 * 1e-3             # 커넥팅 로드의 길이 [m]
r      = 6.5   * 1e-3             # 모터의 회전 반경   [m]
m_head = 50    * 1e-3             # 헤드의 질량        [kg]
omega  = rpm   * (2 * np.pi / 60) # 모터의 각속도      [rad/s]

# 살의 물성치
m      =   1.54  # 살의 질량     [kg]    (정당한지 모르겠음)
c      =   7.17  # 살의 감쇠계수 [N*s/m] (정당한지 모르겠음)
k      = 105.72  # 살의 강성     [N/m]   (정당한지 모르겠음)

t_sym = sp.Symbol('t')
# 헤드의 변위 [m]
x_head_sym = r * sp.cos(omega * t_sym) + sp.sqrt(l**2 - (r * sp.sin(omega * t_sym))**2)
# 헤드의 속도 [m/s]
v_head_sym = sp.diff(x_head_sym, t_sym)
# 헤드의 가속도 [m/s^2]
a_head_sym = sp.diff(v_head_sym, t_sym)
# 헤드가 가하는 힘 [N]
def F_head(t):
    if abs(r * np.sin(omega * t)) < y_cr and np.cos(omega * t) > 0:
        return - m_head * a_head_sym.subs(t_sym, t)
    else:
        return 0

# m x'' + c x' + k x = F_head
def phi(t, y):
    return np.array([y[1], - (c / m) * y[1] - (k / m) * y[0] + F_head(t) / m])

t0 = 0                     # 시작 시각               [s]
tf = 10                    # 종료 시각               [s]
n  = 10000                 # 점들의 갯수
dt = (tf - t0) / (n - 1.0) # 인접한 두 점 사이의 간격

x0 = 0                     # 처음 위치    [m]   (정당한지 모르겠음)
v0 = 0                     # 처음 속도    [m/s] (정당한지 모르겠음)

sol = RK4(phi, t0, np.array([x0, v0]), tf, dt)

ts   = list(map(lambda txv: txv[0], sol))           # [s]
xs   = list(map(lambda txv: 1000 * txv[1][0], sol)) # [mm]
imps = sum(map(lambda t: abs(F_head(t)) * dt, ts))  # [N*s] (0~10)[s]의 충격량의 총합
print("y_cr =", y_cr, "the sum of impulse =", imps)

F_head_max = max(map(lambda t: abs(F_head(t)), ts))
F_head_max_s = max([abs(F_head(t)) for t in ts if t > tf / 2])

print("F_head_max =", F_head_max)
print("F_head_max_s =", F_head_max_s)

# 플롯
plt.title("rpm = " + str(rpm))
plt.plot(ts, xs)
plt.xlabel('t [s]')
plt.ylabel('x [mm]')
plt.grid()
plt.show()