코딩일지(2024-05-12)

계측공학에서 배운걸로 5장 문제를 풀어보자!

---

저번에 5장 숙제인 9,14,19번 문제를 풀었고, 오늘은 나머지 문제들을 해결해보려고 한다.

5장 9번

5.9. 초기 부피 𝑉0를 가진 물이 담긴 탱크가 탱크 바닥의 구멍을 통해 물을 배출한다. 배출 속도가 탱크 안의 물의 부피에 비례한다면, 이 시스템의 시간 상수를 구하라.

1>배수율에 대한 관계 설정

2>미분방정식 풀이

3>시정수 구하기

시정수는 기존값의 36.8%가 되는 지점의 x값을 말하는 것이다. 보통 그림을 그려보면 다음과 같이 상승하거나 하강하는 양상의 그림이 그려진다.

이를 참고하여 다음과 같이 실제코드로 작성해보았다.

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# 주어지 초기 조건 및 상수
V0 = 1000 #초기 물의 양 (예를 들어, 1000리터)
k = 0.1   #비례 상수 (임의의 값, 단위: 1/초)

# 시정수 계산
tau = 1/k
print(f"시정수 (tau): {tau}초")

# 시간 범위 설정 (0부터 5*tau까지)
t_span = (0, 5 * tau)
t_eval = np.linspace(*t_span, 1000)

# 미분방정식 정의
def dVdt(t,V):
    return -k * V

# 미분 방정식 풀이
sol = solve_ivp(dVdt, t_span, [V0],  t_eval=t_eval)

# 결과 플로팅
plt.plot(sol.t, sol.y[0], label=f"V(t) = {V0} * exp(-{k}t)")
plt.axhline(0, color='grey', lw=0.5) # x축과 평핸한 점근선
plt.axvline(tau, color='red')
plt.xlabel('Time (s)')
plt.ylabel('Volume (L)')
plt.title('Volume of Water in Tank Over Time')
plt.legend()
plt.grid()
plt.show()

results>

5장 14번

5.14. 압력 변환기는 2차 시스템으로 작동한다. 감쇠되지 않은 고유 주파수가 4000 Hz이고 감쇠가 임계의 75%일 때, 측정 오차가 5%를 넘지 않는 주파수 범위를 구하라.

단순히 수식으로만 접근하면 다음과 같이 구현해볼 수 있다.

import sympy as sp

# 주어진 값들
omega_n = 4000
zeta = 0.75

# 주파수 변수를 정의
omega = sp.symbols('omega', real=True, positive=True)


# 정의
lhs = omega_n**4/0.9025 # lhs: left half side
rhs = (omega_n**2-omega**2)**2 + (2 * zeta * omega_n * omega)**2

# 해결
solution = sp.solve(lhs - rhs, omega)
solution

result>

[1904.31399610052]

근데 이 문제를 풀면서 한가지 의문이 들었던 부분이 있었는데, 바로 0.95 <=

H(jw)

<= 1.00에서 실질적으로 접근을 할때, 0.95 =

H(jw)

으로 두고 방정식을 푸는 것에 초점을 맞춘다는 것이다. 왜 1 =

H(jw)

에 대해서는 다루지 않을까 하는 의문이 들었다. 그리고 이 의문은 다음 그래프를 구현하여 해석해봄으로써 이해할 수 있었다.

실제 w와 H(jw)크기에 대한 관계 그래프를 나타내보면 다음과 같다.

import numpy as np
import matplotlib.pyplot as plt

# 주어진 값
omega_n = 4000 # 자연 진동수 (Hz)
zeta = 0.75    # 감쇠 비율

# 주파수 범위 설정
omega = np.linspace(0, 2 * omega_n, 10000)

# 주파수 응답 함수 정의
def H(omega, omega_n, zeta):
    numerator = omega_n**2
    denominator = np.sqrt((omega_n**2 - omega**2)**2 + (2 * zeta * omega_n * omega)**2)
    return numerator / denominator

# 주파수 응답 계산
H_omega = H(omega, omega_n, zeta)

# 5% 오차 범위
lower_bound = 0.95
upper_bound = 1.05

# 오차 범위 내의 주파수 찾기
valid_indices = np.where((H_omega >= lower_bound) & (H_omega <= upper_bound))
valid_frequencies = omega[valid_indices]

# 결과 출력
print(f"Measurement error not greater than 5% for frequencies in the range {valid_frequencies[0]:.2f} Hz to {valid_frequencies[-1]:.2f} Hz")

# 결과 플로팅
plt.plot(omega, H_omega, label='|H(jw)|')
plt.axhline(1, color='red', linestyle='--', label='|H(jw)| = 1')
plt.axhline(upper_bound, color='green', linestyle='--', label='Upper bound (1.05)')
plt.axhline(lower_bound, color='blue', linestyle='--', label='Lower bound (0.95)')
plt.fill_between(omega, lower_bound, upper_bound, where=((H_omega >= lower_bound) & (H_omega <= upper_bound)), color='grey', alpha=0.3)
plt.xlabel('Frequency (Hz)')
plt.ylabel('|H(jw)|')
plt.title('Frequency Response of the Second-Order System')
plt.legend()
plt.grid()
plt.show()

results>

5장 19번

a)시정수 구하기:

b)스위치가 닫힌 후 0.5초 뒤의 저항 양단 전압:

이를 코드를 풀어보면 다음과 같아지고,

import numpy as np

# 주어진 값
R = 10 * 10**3     # R: 저항(옴)
C = 200 * 10**(-6) # C: 커패시턴스(패럿)
E = 5              # E: 전원전압(테스트용 값으로 5V로 설정)

# 시정수 계산
tau = R * C
print(f"시정수 (tau): {tau:.2}초\n")

# 특정 시간 (0.5초) 후의 전압 계산
#1. 저항
voltage_across_resistor = E * (np.exp(-t/tau))
print(f"0.5초 후 저항 양단의 전압: {voltage_across_resistor:.2f} V")
ratio = (voltage_across_resistor/E)*100
print(f"5초 뒤의 저항 양단의 전압은 전원 전압E(5V)의 약{ratio:.2f}% 입니다.")

#2. 커패시터
t = 0.5
voltage_across_capacitor = E * (1 - np.exp(-t/tau))
print(f"0.5초 후 커패시터 양단의 전압: {voltage_across_capacitor:.2f} V")

ratio = (voltage_across_capacitor/E)*100
print(f"5초 뒤의 커패시터 양단의 전압은 전원 전압E(5V)의 약{ratio:.2f}% 입니다.\n")

result>

시정수 (tau): 2.0초

0.5초 후 저항 양단의 전압: 3.89 V
5초 뒤의 저항 양단의 전압은 전원 전압E(5V)의 약77.88% 입니다.
0.5초 후 커패시터 양단의 전압: 1.11 V
5초 뒤의 커패시터 양단의 전압은 전원 전압E(5V)의 약22.12% 입니다.

이를 그래프로 시각화해보면 다음과 같더라?(신기신기)

import numpy as np
import matplotlib.pyplot as plt

t = np.linspace(0, 10, 10000)
voltage_across_resistor = E * (np.exp(-t/tau))
plt.plot(t, voltage_across_resistor, label='V(t)')
plt.axvline(0.5, color='red',linestyle='--', lw=1.5, label='0.5s after the switch is closed')
plt.axhline(3.89, color='orange', linestyle='--', lw=1.5, label='voltage across the resistor 0.5s after') # x축과 평핸한 점근선
plt.xlabel('time(sec)')
plt.ylabel('V(t)')
plt.title('Time Response of the One-Order System(in Resistor)')
plt.legend()
plt.grid()
plt.show()

result>

추가적으로 커패시터에 걸리는 전압은 이거에 반비례하게 되니 이것도 확인해보았다.

import numpy as np
import matplotlib.pyplot as plt

t = np.linspace(0, 10, 10000)
voltage_across_resistor = E * (1 - np.exp(-t/tau))
plt.plot(t, voltage_across_resistor, label='V(t)')
plt.axvline(0.5, color='red',linestyle='--', lw=1.5, label='0.5s after the switch is closed')
plt.axhline(1.11, color='orange', linestyle='--', lw=1.5, label='voltage across the resistor 0.5s after') # x축과 평핸한 점근선
plt.xlabel('time(sec)')
plt.ylabel('V(t)')
plt.title('Time Response of the One-Order System(in Capacitor)')
plt.legend()
plt.grid()
plt.show()

result>

손풀이

다음에 나머지 문제들도 풀어보자!