실베스터-전달강성계수법에 의한 실습선 새동백호 추진축계의 비틀림 자유진동 해석

2.1 추진축계의 구성

Fig. 1은 새동백호 추진축계의 개략도이다. 추진축계는 주기관, 탄성 커플링, 감속기, 중간 축, OD(Oil Distributer) 축, 프로펠러 축 그리고 프로펠러로 구성된다. 주기관은 현대중공업의 7H32/40P, 탄성 커플링은 Centa사의 CX-181, 감속기는 Reintjes사의 LAF 5645이고, 축과 프로펠러의 제작사는 WÄRTSILÄ이다. 새동백호 추진축계의 구체적인 스펙을 Table 1에 제시한다.

Fig. 1

Propulsion shafting of training ship AEDONGBAEK

Table 1

Specification of propulsion shafting

2.2 추진축계의 모델링

새동백호 추진축계의 비틀림 자유진동을 해석하기 위하여 추진축계를 다수의 회전 질량과 비틀림 스프링 그리고 기어쌍과 기어축으로 모델링한 것이 Fig. 2이다. Fig. 2에서 절점 1은 캠 구동부이고 절점 2는 플랜지이다. 주기관은 7기통이므로 절점 3부터 절점 9는 순차적으로 크랭크 축과 커넥팅로드와의 연결점이다. 절점 10은 기관의 플라이 휠과 탄성 커플링의 결합점이고, 절점 11은 탄성 커플링과 감속기의 회전 질량과의 결합점이다. 절점 12는 클러치이고, 절점 13은 상부 기어와 하부 기어로 구성된 기어쌍이고, 절점 17은 프로펠러에 상당한다. 절점 14와 절점 15 사이에 중간 축이 있고, 절점 15와 절점 16 사이에 OD 축, 절점 16과 절점 17 사이에 프로펠러 축이 있다.

Fig. 2

Computational model

3.1 유한요소법³⁾

절점 1, 2, 4, 5에 각각 하나의 회전 질량이 존재하므로, 이들의 질량관성모멘트(J_i)를 가지고 크기 1×1의 요소 질량관성행렬을 구한다.

절점 3의 기어 쌍은 상부와 하부에 각각 질량관성모멘트 $$ {\widehat{J}}_{up}$$와 $$ {\widehat{J}}_{down}$$인 2개의 회전 질량으로 구성되므로, 상부 회전 질량을 중심으로 등가화하여 질량관성모멘트를 구하면

여기서 r은 아래 식과 같이 상부 기어와 하부 기어의 회전 속력비이다.

따라서 식 (3)의 등가 질량관성모멘트(J₃)를 식(1)에 적용하면 절점 3의 요소 질량관성행렬을 구할 수 있다.

인접한 두 절점을 연결하는 비틀림 스프링 1, 2, 4는 각각의 스프링상수(K_i)를 이용하여 식 (4)와 같은 요소 강성행렬을 유도할 수 있다.

유한요소해석³⁾으로부터 하부 기어($$ {\widehat{J}}_{down}$$)와 회전 질량(J₄)을 연결하는 기어 축의 요소 강성행렬은 기어 축의 비틀림 강성이 $$ {\widehat{K}}_i$$인 경우 다음과 같다.

요소행렬들을 모두 조립하면, Fig. 3의 계를 다음 식과 같은 전체 질량관성행렬(J)과 전체 강성행렬(K)로 나타낼 수 있다.

전체 계의 운동방정식을 행렬로 나타내면

여기서 ϴ는 각 절점의 각변위(ϴ_i)로 구성되는 전체 각변위벡터(ϴ = {ϴ₁, ϴ₂, ϴ₃, ϴ₄, ϴ₅}^T)이고, $$ \ddot{\theta }$$는 전체 각가속도벡터이다.

식 (8)을 일반 고유치 문제로 나타내면

가 되고, 여기서 λ는 고유치로서 각고유진동수 ω_n의 제곱이고, ψ는 고유벡터이다. 식 (9)의 우변 항을 좌변으로 넘기면

따라서 진동수 방정식은

가 되어, 식 (11)에서 고유치(λ_i)를 구하고, 식 (10)에서 고유벡터(ψ_i)를 구할 수 있다.

만약 계의 전체 자유도가 커지면 자코비 법(Jacobi method), 역방향 반복법(inverse iteration method), 축소공간 반복법(subspace iteration method) 등의 고유치 문제 해법⁶⁾으로 고유치와 고유벡터를 구하는 것이 효과적이다.

3.2 전달강성계수법⁵⁾

3.2.2 동강성계수의 전달

Fig. 4와 같이 절점 i에 질량관성모멘트가 J_i인 회전 질량이 존재하는 경우, 절점 i의 좌측에서 절점 i의 우측으로 동강성계수의 전달식은

$\[ S_i={\overline{S}}_i+\left(-{\omega }^2{J}_i\right) \left(i=\mathrm{1,2},\mathrm{4,5}\right) \]$

(14)

된다.

Fig. 4

Rotating masses and torsional spring

Fig. 4에서 절점 i와 절점 i+1 사이에 스프링상수가 K_i인 비틀림 스프링이 존재하는 경우, 절점 i와 절점 i+1의 각변위 사이의 관계는

$\[ {\theta }_i=V_i\times {\theta }_{i+1} \left(i=\mathrm{1,2},4\right) \]$

(15)

가 되고, 여기서

$\[ V_i=K_i/G_i, G_i=S_i+K_i \left(i=\mathrm{1,2},4\right) \]$

(16)

이고, 절점 i의 우측에서 절점 i+1의 좌측으로 동강성계수의 전달식은

$\[ {\overline{S}}_{i+1}=K_i\times \left(1-V_i\right) \left(i=\mathrm{1,2},4\right) \]$

(17)

이 된다.

Fig. 5에 나타낸 절점 i의 기어 쌍은 상부 회전 질량을 중심으로 한 등가 질량관성모멘트($$ J_i={\widehat{J}}_{up}+\left(r^2+{\widehat{J}}_{down}\right)$$)를 식 (14)에 적용하면 기어 쌍에 대해 동강성계수를 전달할 수 있다.

Fig. 5

Gear pair and gear shaft

Fig. 5에서 하부 기어($$ {\widehat{J}}_{down}$$)와 다음 회전원판(J_i+1) 사이에 연결된 기어 축의 비틀림 스프링상수가 $$ {\widehat{K}}_i$$라면, 이 스프링 요소의 좌측과 우측의 토크와 각변위 사이의 관계는 다음 식과 같다.

$\[ \left[\begin{array}{c}{T}_i\\ {\overline{T}}_{i+1}\end{array}\right]=\left[\begin{array}{cc}-r^2{\widehat{K}}_i& -r{\widehat{K}}_i\\ r{\widehat{K}}_i& {\widehat{K}}_i\end{array}\right]\left[\begin{array}{c}{\theta }_i\\ {\theta }_{i+1}\end{array}\right] \left(i=3\right) \]$

(18)

식 (18)의 T_i에 식 (13)을 대입하고 정리하면

$\[ {\theta }_i=V_i\times {\theta }_{i+1} \left(i=3\right) \]$

(19)

가 되고, 여기서

$\[ V_i=-r{K}_i/G_i,G_i=S_i+r^2{K}_i \left(i=3\right) \]$

(20)

식 (18)의 $$ {\overline{T}}_{i+1}$$에 식 (12)의 “i” 대신에 “i+1”을 대입한 식과 식 (19)을 대입하여 정리하면, 기어 축의 좌측에서 우측으로 동강성계수를 전달하는 아래 식을 유도할 수 있다.

$\[ {\overline{S}}_{i+1}={\widehat{K}}_i\times \left(1+r{V}_i\right) \left(i=3\right) \]$

(21)

Fig. 3의 계산 모델에 상기 식 (14), 식 (17), 식 (21)을 요소 별로 순차적으로 적용하면 최종적으로 마지막 절점인 절점 5의 우측의 동강성계수(S₅)를 구할 수 있다.

3.3 실베스터-전달강성계수법^1-2)

식 (9)의 일반 고유치 문제는 아래와 같은 직교성을 만족한다.

여기서 Φ = [ψ₁, ψ₁, ψ₂, ψ₃, ψ₄, ψ₅]는 1차부터 5차까지의 고유벡터로 구성된 모드행렬이다.

행렬 K - λJ의 앞과 뒤에 모드행렬을 식 (27)과 같이 곱하면

대각행렬이 되고 대각항의 요소들은 (λ₁ - λ), (λ₂ - λ), …, (λ₅ - λ)가 된다. 따라서 (K - λJ) 행렬의 고유치는 (λ_i - λ)가 된다. 그리고 이들 고유치 중에서 음수의 개수는 J 및 K 행렬의 고유치(λ_i) 중에서 λ보다 작은 고유치의 개수와 같다.

(K - λJ) 행렬의 앞뒤에 아래에 나타낸 행렬 L과 전치행렬인 L^T를 곱하면 다음과 같은 대각행렬이 된다.

여기서

실베스터의 관성이론으로부터 (K - λJ)행렬의 고유치 중에서 음의 고유치의 개수는 행렬 diag[G₁, G₂, G₃, G₄, G₅]의 대각요소의 음수의 개수와 동일하다. 그리고 식 (28)과 식 (29)의 G_i, V_i, S_i는 3.2절에 소개한 전달강성계수법의 G_i, V_i, S_i와 동일하다.

따라서 Fig. 3의 해석 모델에 대해 실베스터-전달강성계수법을 이용하여 고유진동수를 구하는 전산 알고리즘은 다음과 같다.

1) 임의의 각진동수 $$ \tilde{\omega }$$를 정한다.

2) 요소별 전달강성계수법의 전달식을 적용하여 모든 절점의 G_i와 S_i를 계산한다.

3) G₁부터 G₄ 그리고 S₅ 중에서 0보다 작은 값을 가지는 변수의 개수(m)를 센다.

4) 상기 3)의 개수(m)는 검색 진동수 $$ \tilde{\omega }$$ 보다 작은 고유각진동수의 개수이다. 진동수 검색 범위 안의 여러 $$ \tilde{\omega }$$에 대해 반복적으로 m을 구하면 빠뜨리지 않고 모든 고유각진동수를 구할 수 있다.

4.1 가변 피치 프로펠러 상태에 따른 추진축계의 비틀림 자유진동 해석

새동백호는 가변 피치 프로펠러가 장착되어 있으므로 프로펠러의 상태에 따라 프로펠러의 질량관성모멘트의 값이 달라진다.

탄성 커플링의 온도를 50℃로 설정하고, 프로펠러가 zero-pitch와 full-pitch일 때, 실베스터-전달강성계수법으로 추진축계의 고유진동수를 각각 계산하였다. 강체모드는 생략하고, 1차부터 5차까지의 고유진동수를 계산한 결과를 Table 2에 제시하였다.

Table 2

Natural frequencies of computational model (temperature of coupling : 50℃) [rpm]

Table 2로부터 프로펠러의 상태에 따른 고유진동수를 비교해 보면, 프로펠러가 zero-pitch에서 full-pitch로 변경되면 1차 고유진동수는 약 9% 정도 감소되고, 2차 고유진동수는 약 1% 정도 감소되었다. 그러나 3차 이상의 고유진동수에는 거의 변화가 없었다.

zero-pitch 및 full-pitch 2가지 경우에 대해 고유모드를 각각 계산하였다. Fig. 6에는 zero-pitch (o 표시)와 full-pitch(x 표시) 상태의 1차 고유모드를 제시하였다. 1차 고유모드를 살펴보면 탄성 커플링(절점 10과 절점 11사이)에서 진동을 하지 않는 절이 발생하고 있음을 알 수 있었다. 프로펠러의 상태에 따른 고유모드의 값을 비교해 보면 감속기 상부 기어의 좌측(절점 1에서 절점 13)에서 차이가 있음을 알 수 있었다.

Fig. 6

1st natural mode

Fig. 7은 2차 고유모드를 계산한 결과이다. 2차 고유모드는 탄성 커플링(절점 10과 절점 11 사이)과 프로펠러 축(절점 16과 절점 17 사이)에서 절이 발생하고 있음을 알 수 있었다. 그리고 프로펠러의 피치 값에 무관하게 2차 고유모드는 거의 동일함을 알 수 있었다.

Fig. 7

2nd natural mode

그리고 저자들은 프로펠러의 상태에 따른 3차부터 5차까지의 고유모드를 구한 후, 비교해 본 결과 거의 차이가 없음을 확인하였다.

4.2 탄성 커플링의 온도 변화에 따른 추진축계의 비틀림 자유진동 해석

새동백호 추진축계의 탄성 커플링은 온도에 따라 비틀림 강성이 변화한다. 커플링 제작사로부터 제공 받은 자료에 의하면 절점 10과 11 사이에 있는 커플링의 비틀림 스프링상수(K₁₀)의 값은 30℃에서 0.3162 MNm/rad이고, 70℃에서 0.2652 MNm/rad이었다.

프로펠러가 full-pitch이고, 커플링의 온도가 30℃와 70℃일 때, 실베스터-전달강성계수법으로 추진축계의 고유진동수를 각각 계산하였다. 강체모드는 생략하고, 1차부터 5차까지의 고유진동수를 계산한 결과를 Table 3에 제시하였다. 계산된 고유진동수를 비교해 보면, 1차 고유진동수는 약 4%의 변화가 있었고, 2차 고유진동수는 약 5%의 변화가 있었다. 그리고 탄성 커플링의 온도가 상승할 때 고유진동수가 감소함을 알 수 있었다. 그러나 3차, 4차, 5차 고유진동수는 거의 변화가 없었다.

Table 3

Natural frequencies of computational model (state of propeller : full-pitch) [rpm]

30℃와 70℃, 2가지 경우에 대해 고유모드를 각각 계산하였다. Fig. 8에는 탄성 커플링의 온도인 30℃(o 표시)와 70℃(x 표시)의 1차 고유모드를 제시하였다. 1차 고유모드를 살펴보면 탄성 커플링(절점 10과 절점 11사이)에서 절이 발생하고 있음을 알 수 있었다. 탄성 커플링의 온도에 따른 고유모드의 값을 비교해 보면 기어 박스에 해당하는 절점 11과 절점 14 사이에서 약간 차이가 있음을 알 수 있었다.

Fig. 8

1st natural mode

Fig. 9은 2차 고유모드를 계산한 결과이다. 2차 고유모드는 탄성 커플링(절점 10과 절점 11 사이)과 프로펠러 축(절점 16과 절점 17 사이)에서 절이 발생하고 있음을 알 수 있었다. 그리고 탄성 커플링의 온도와 무관하게 2차 고유모드는 거의 동일함을 알 수 있었다.

Fig. 9

2nd natural mode

그리고 저자들은 탄성 커플링의 온도에 따른 3차부터 5차까지의 고유모드를 구한 후, 비교해 본 결과 거의 차이가 없음을 확인하였다.