[기계공학/진동및소음] 05. 다자유도의 진동해석

교재: Julia 프로그래밍을 활용한 구조물의 진동해석 (곽문규, 김수민)

 

5.1 서론

 

 다자유진동계(Multi-Degree-of-Freedom System)는 독자적으로 움직일 수 있는 물체가 여러 개 있다는 것을 의미한다. 다자유도계의 짇농을 해석하기 위해서는 행렬의 사용이 필수적이다. 먼저 2 자유도계로 시작해보자.

 

5.2 비감쇠 다자유도계의 자유진동

 

 다음과 같은 2 자유도 진동계를 고려해보자.

 

 

 여기서 x₂ > x₁ 이라고 가정하고 자유물체도를 그려보자.

 

 

 뉴턴 제 2법칙을 적용하여 운동방정식을 유도해보자.

 

m₁x₁'' + (k₁ + k₂)*x₁ -k₂x₂ = 0

m₂x₂'' - k₂x₁+ (k₂ + k₃)*x₂ = 0

 

 이를 행렬식으로 간단히 표현하면 다음과 같이 쓸 수 있다.

 

Mx'' + Kx = 0

 

 여기서 x = [x₁ x₂]ᵀ 이고, M = [m₁ 0; 0 m₂],    K = [k₁+k₂ -k₂; -k₂ k₂+k₃] 이다. M과 K는 보통 질량 행렬과 강성행렬로 불린다. 이때 선형 상미분방정식 형태는 동일하기에 일반해를 다음과 같이 설정하자.

 

$x=X\times e^{st}=X\times e^{iwt}$

 

 일반해를 운동방정식에 대입하면 다음과 같은 식이 만족되어야 한다.

 

$KX=\lambda MX$

 

 여기서 λ = ω²이다. 해당 식은 만족시키기 위해서 결정자(Determinant)는 다음과 같이 유도된다.

 

$(m_1m_2)\times \lambda ^2-(m_1k_2+m_1k_3+m_2k_1+m_2k_2)\times \lambda +(k_1k_2+k_1k_3+k_2k_3)=0$

 

 근의 공식을 이용하여 위 2차방정식을 만족하는 2개의 λ, 즉 λ₁, λ₂를 구할 수 있다. 진동에 있어 가장 큰 영향을 주는 고유진동수는 제일 작은 고유진동수이기 때문에 λ₁ < λ₂ 순으로 정렬하는 것이 좋다. 그리고 고유진동수는 다음과 같이 계산한다.

 

$\omega _1=\sqrt{\lambda_1}, \quad \omega _2=\sqrt{\lambda_2}$

 

 해당 고유진동수에 맞는 고유벡터 또한 구할 수 있다.

 

Ex) $m_1 = m_2 = 1kg,\quad k_1 = k_2 = k_3 = 1N/m$ 라고 하자.

using LinearAlgebra


M = [1 0; 0 1]
K = [2 -1; -1 2]

val, vec = eigen(K, M)

display(val)
display(vec)

 

 val = [1 3]ᵀ 그리고 vec = [-0.707 -0.707; -0.707 0.707]로 나온다.

 

 자유진동 문제를 다루는데 초기조건, 즉 x₁(0) = x₁₀, x₁'(0) = v₁₀, x₂(0) = x₂₀, x₂'(0) = v₂₀일 경우 응답을 구할 수 있을까? 이는 많은 수학 연산을 요구한다. 이를 보다 간단하게 풀기 위해 수학자들이 모달 해석(Modal Analysis)이라는 방법을 개발했다.

모드해석의 이해

 

5.3 모달 해석에 의한 다자유도 비감쇠 진동계의 자유진동 응답

 위 예제 $m_1 = m_2 = 1kg,\quad k_1 = k_2 = k_3 = 1N/m$ 의 2 자유도 진동계에 대한 초기조건이 x₁(0) = 1, x₁'(0) = 0, x₂(0) = 0, x₂'(0) = 0일 경우의 자유진동 응답을 계산해보자. Mx'' + Kx = 0에서는 x₁과 x₂가 연성(couple)되어 있는 2차 상미분 연립방정식이다. 그래서 수학적으로 엄밀해를 구하기 어렵다. 이를 해결하기 위해 모달 좌표계(Modal Coordinate)를 도입하는 개념을 제시하였다.

 

Mx'' + Kx = 0,    x(0) = x₀,    x'(0) = v₀

 

 고유치 문제를 플어 얻은 고유치 행렬과 고유벡터 행렬을 각각 Λ와 U라고 하자. 고유치행렬 Λ는 대각행렬(Diagonal Matrix)인데 제일 작은 고유치부터 순차적으로 정렬되어 있는 행렬이다. 그리고 고유벡터 행렬의 각 열의 크기를 조정하면 다음과 같은 정상직교조건(orthonormality condition)을 만족하게 만들 수 있다.

 

$U^TMU=I, U^TKU=\Lambda$

 

 위 정상직교 조건은 진동계 해석에 있어 중요한 역할을 한다.

 

 다음과 같은 좌표 변환을 도입해보자.

 

x = Uq

 

 이 변환을 모달 변환(Modal Transformation)이라고 부른다. 여기서 x는 실제좌표(Physical coordinates), q는 모달좌표(Modal coordinates)이다. 위 식은 실제좌표와 모달좌표 사이의 관계를 나타내는 식인데 고유벡터 행렬 U가 좌표변환을 주도한다는데 그 특징이 있다.

 

 이 절의 첫 번째 식에 좌표변환을 도입하고 $U^T$ 를 앞에 곱하여 정상직교조건을 만족시키면 다음과 같다.

 

q'' + Λq = 0

 

 위 식을 2자유도계에 대해 풀어 쓰면 다음과 같이 쓸 수 있다.

 

$q_1''+\omega _1^2\times q_1=0, \quad q_2''+\omega _2^2\times q_2=0$

 

이 식들은 q₁ 과 q₂가 서로 연성되어 있지 않다. 이게 의미하는 바는 각각 풀어도 된다는 것이다. q₁ 과 q₂의 해는 다음과 같이 쓸 수 있다.

 

$q_1=A_1\times \sin (\omega _1 t)+B_1\times \cos (\omega _1 t), \quad q_2=A_2\times \sin (\omega _2 t)+B_2\times \cos (\omega _2t) $

 

A₁, B₁, A₂, B₂는 초기조건에 의 결정된다. 모달 좌표계에 대한 초기조건은 다음식을 이용해 계산하면 된다.

 

$q(0)=U^{-1}x(0), \quad q'(0)=U^{-1}x'(0)$

 

Ex) $m_1 = m_2 = 1kg,\quad k_1 = k_2 = k_3 = 1N/m$ 라고 하자. 초기조건 x₁(0) = 1, x₁'(0) = 0, x₂(0) = 0, x₂'(0) = 0에 대한 자유진동응답을 계산해보자.

using LinearAlgebra
using SymPy
using Plots


q1 = SymFunction("q1")
q2 = SymFunction("q2")
t = symbols("t")


M = [1 0; 0 1]
K = [2 -1; -1 2]
ω2, U = eigen(K, M)


#display(ω2);
ω = sqrt.(ω2);

#display(ω)
display(U)

# Orthonormality check

# x domain -> q domain

x0 = [ 1, 0 ];
xd0 = [ 0, 0 ]
q0 = U'*M*x0;
qd0 = U'*M*xd0;
q10 = q0[1]; q20 = q0[2];
qd10 = qd0[1]; qd20 = qd0[2];

diffeq=Eq(q1(t).diff(t,t)+q1(t),0)
ics=((q1,0.0,q10),(q1',0.0,qd10))
res1=dsolve(diffeq,ics=ics)
display(res1)

diffeq2=Eq(q2(t).diff(t,t)+3*q2(t),0)
ics=((q2,0.0,q20),(q2',0.0,qd20))
res2=dsolve(diffeq2,ics=ics)
display(res2)

#Back to original coordinates

x = U*[ res1.rhs, res2.rhs ]
display(x)
t=0:0.1:10

p2=plot(t,x[1],label="x1(t)")
p2=plot!(t,x[2],label="x2(t)")
display(p2)

 

5.4 다자유도 감쇠 진동게의 자유진동 응답

 

 다자유도 진동계에 댐퍼를 추가한 다음과 같은 진동계를 고려해보자.

 마찬가지로 x₂ > x₁ 이라고 가정하고 자유물체도를 그려보자.