[ Article ]

Transactions of the Korean Society for Noise and Vibration Engineering - Vol. 29, No. 2, pp.175-184

ISSN: 1598-2785 (Print) 2287-5476 (Online)

Print publication date 20 Apr 2019

Received 30 Sep 2018 Revised 11 Dec 2018 Accepted 11 Feb 2019

DOI: https://doi.org/10.5050/KSNVE.2019.29.2.175

축-원판-블레이드 연성 계의 고유진동수 수렴특성 연구

김형희^* ; 유홍희^†

Study on the Natural Frequency Convergence Characteristics of a Shaft-disk-blade Coupled System

Hyung Hee Kim^* ; Hong Hee Yoo^†

*Dept. of Mech. Convergence Eng., Hanyang University

Correspondence to: ^† Fellow Member, School of Mechanical Engineering, Hanyang University E-mail: hhyoo@hanyang.ac.kr ‡ Recommended by Editor Gi-Woo Kim

Abstract

Shaft-disk-blade coupled systems are widely used in rotating machinery such as helicopters, wind turbines, and steam turbines. In this study, an analytical model of a shaft-disk-blade coupled system is developed by considering the coupling effects in shaft bending, shaft torsion, shaft shear, and blade bending. Next, the natural frequencies of the system are obtained using the assumed mode method that employs three types of mode functions (simple beam functions, beam functions considering the effect of lumped end-mass, and beam functions considering the effects of lumped end-mass and moment of inertia). The results obtained by the three mode function sets are compared to those obtained using a commercial FE program, and the natural frequency convergence characteristics are discussed. It was found that only one of the mode function sets guaranteed accurate natural frequency convergence.

Keywords:

Shaft-disk-blade Coupled System, Assumed Mode Method, Natural Frequency, Convergence Characteristics

키워드:

축-원판-블레이드 연성계, 가상모드 방법, 고유진동수, 수렴 특성

1. 서 론

축-원판-블레이드 연성계는 발전용 스팀터빈, 가스터빈 엔진, 헬리콥터 회전익 계, 그리고 풍력발전기 등 다양한 회전기기에 사용되고 있다. 회전기기에는 진동이 어느 정도 항상 발생하므로 그 진동 특성은 시스템을 평가하고 설계하는 데에 대단히 중요한 요소이다. 따라서 많은 연구자들이 축, 원판, 그리고 블레이드에 대하여 개별적으로 연구하고 있을 뿐만 아니라 축-원판-블레이드 간의 연성 효과를 고려한 고유진동수 변화에 대해서도 연구를 수행하여 왔다.

Lee⁽¹⁾는 substructure synthesis 방법과 가상모드법을 사용하여 블레이드의 장착 각과 비틀림 각에 따른 축-원판-블레이드 연성계의 고유진동수 특성변화에 대한 연구를 수행하였고 Al-bedoor⁽²^,³⁾는 축의 비틀림과 블레이드의 굽힘이 고려된 축-원판-블레이드 연성계 모델을 제안하고 축의 비틀림 강성에 따른 계의 고유진동수 변화에 대하여 연구하였다. 또한, Yang⁽⁴⁾은 원판의 강성의 영향에 따른 축-원판-블레이드 연성계의 고유진동수 변화에 대해 연구했으며 Turban⁽⁵⁾은 다단 축-원판-블레이드 연성 계의 모델을 제안하고 원판과 블레이드 개수 증가에 따른 고유진동 특성변화에 대하여 연구하였다.

축-원판-블레이드 연성 계에 대한 진동 해석에는 운동방정식을 유도하기 위해 가상모드법이 종종 사용되고 있다. 하지만 적절한 모드함수가 사용되지 않으면 구해지는 고유진동수가 수렴하지 못할 뿐만 아니라 계산에 많은 시간을 낭비하게 된다. 따라서 적절한 모드함수 선택은 진동해석에 매우 중요하다.

이 연구의 주요목적은 축의 비틀림, 굽힘, 전단과 블레이드 굽힘이 고려된 축-원판-블레이드 연성계를 모델링하고 운동방정식을 유도하여 축의 굽힘 변위 및 각도에 대한 세 종류의 모드 함수를 사용할 때 모드 수 증가에 따른 모드 함수 종류별 고유진동수 수렴 특성을 파악하는 것이다. 이와 같은 연성 구조 모델의 정확한 수렴특성은 모델의 안정적인 사용을 위해 반드시 연구되어야 하는 필수사항이다.

이 논문의 구성은 다음과 같다. 다음 2장에서는 복합변형변수를⁽⁶⁾ 이용 축-원판-블레이드 연성 계를 Kane 방법⁽⁷⁾을 사용하여서 선형 운동방정식을 유도하였다. 3장에서는 회전속도가 0일 때, 축의 굽힘 변형에 의한 변위 및 각도와 관련된 세 종류의 모드 함수를 사용할 경우 각 모드함수 별로 2장에서 유도된 운동방정식에 근거하여 모드 수 증가에 따른 해석결과를 구하고 ansys 결과와 비교하여 그 수렴 특성을 파악하였다. 마지막으로 4장에는 결론을 수록하였다.

2. 시스템 모델링

2.1 사용된 가정

여러 개의 블레이드를 갖는 경우 블레이드 간 연성 효과가 나타나게 되나 이 연구에서는 그러한 연성 효과는 고려하지 않았다. 첫째, 축은 직경이 일정한 원형 단면을 갖고 블레이드의 단면은 폭과 두께가 일정한 사각 단면이다. 원판은 원형 강체로 가정하였다. 둘째, 축과 블레이드 물성은 균일하다. 셋째, 축은 비틀림과 전단 및 단면관성효과를 고려하였다 (티모센코 보 이론을 사용)⁽⁸⁾. 그러나 블레이드는 블레이드의 길이에 비해 두께가 상대적으로 매우 작으므로 전단 및 단면관성효과를 고려하지 않았으며 블레이드의 비틀림 효과도 비틀림 고유진동수가 굽힘 고유진동수에 비해 훨씬 크므로 이 연구에서는 이를 고려하지 않았다(오일러-베르누이 보 이론을 사용)⁽⁸^,¹⁵⁾.

2.2 운동방정식 유도

Fig. 1에서 $a^1$ , $a^2$ , $a^3$ 는 기준틀 A에 부착된 서로 수직인 단위벡터들이다. 이 구조계는 $a^1$ 방향으로 ω 크기로 회전운동을 하는 외팔보 형태 축 끝단이 원판 중심 D점과 연결된다. 반경 r인 원판 원주 상에는 블레이드가 부착되는데 그 위치를 R이라 한다. 여기서 x는 고정점 O로부터 축 중심선의 변형 전 임의 점까지 거리를, s₁은 축 중심선 상 임의 점의 인장 변형량을, $u → 1$ 은 u₁, u₂, u₃들로 구성된 탄성 변위 벡터를, z는 점 R에서 블레이드 중립 축 상 변형 전 임의 점까지 거리를, s₂는 블레이드 중립 축 상 임의 점의 인장변형량을, 그리고 $u → 2$ 는 u₄, u₅, u₆로 이루어진 탄성 변위 벡터를 나타낸다.

Fig. 1

Configuration of a shaft-disk-blade system

s 1 = ∑ i = 1 μ 1 ϕ 1 i (x) q 1 i (t) i = 1,2, ⋯, μ 1 u 2 = ∑ i = 1 μ 2 ϕ 2 i (x) q 2 i (t) i = 1,2, ⋯, μ 2 u 3 = ∑ i = 1 μ 3 ϕ 3 i (x) q 3 i (t) i = 1,2, ⋯, μ 3 θ 1 = ∑ i = 1 μ 4 ϕ 4 i (x) q 4 i (t) i = 1,2, ⋯, μ 4 θ 2 = ∑ i = 1 μ 5 ϕ 5 i (x) q 5 i (t) i = 1,2, ⋯, μ 5 θ 3 = ∑ i = 1 μ 6 ϕ 6 i (x) q 6 i (t) i = 1,2, ⋯, μ 6 s 2 = ∑ i = 1 μ 7 ϕ 7 i (z) q 7 i (t) i = 1,2, ⋯, μ 7 u 4 = ∑ i = 1 μ 8 ϕ 8 i (z) q 8 i (t) i = 1,2, ⋯, μ 8 u 5 = ∑ i = 1 μ 9 ϕ 9 i (z) q 9 i (t) i = 1,2, ⋯, μ 9

(1)

θ₁, θ₂, θ₃는 축 중심선 상 임의 점에서의 $b^1$ 방향 비틀림각과 $b^2$ 및 $b^3$ 방향 두 굽힘각을 나타낸다.

연속체 방정식을 구하기 위해 가상모드법을 이용 s₁, u₂, u₃, θ₁, θ₂, θ₃, s₂, u₄, u₅를 다음 식과 같이 근사화를 한다. 여기서 ϕ₁_i(x), ϕ₇_i(z)는 축과 블레이드 인장변위, ϕ₂_i(x), ϕ₃_i(x)는 축의 $a^2$ , $a^3$ 방향 굽힘 변위, ϕ₄_i(x), ϕ₅_i(x), ϕ₆_i(x)은 각각 축의 $a^1$ 방향의 비틀림 각과 $a^2$ , $a^3$ 방향의 탄성 축 기울기 그리고 ϕ₈_i(z)과 ϕ₉_i(z)는 블레이드 $b^1$ , $b^2$ 방향 굽힘변위들을 근사화하기 위한 모드함수들이다. 또한, q₁_i(t) ~ q₉_i(t)는 모드 좌표들이며 μ₁~ μ₉는 사용된 모드함수별 모드 좌표의 수이며 모드 좌표의 총합은 시스템의 총 자유도와 같다.

Kane의 방법⁽⁷⁾을 이용하여 운동방정식은 다음의 식으로 나타낼 수 있다.

F i + F i * = 0 1,2, ⋯, μ 1 + μ 2 + ⋯ μ 8 + μ 9

(2)

여기서 F_i와 F_i^*는 일반작용력과 일반관성력으로 다음과 같이 구할 수 있다.

F i = - ∂ U ∂ q i F ∼ i * = - ∫ 0 L 1 ρ A 1 a ~ Q ⋅ v ~ i Q + ω ~ i B ⋅ α ~ B ⋅ I ⃡ B + ω ~ B × I ⃡ B ⋅ ω ~ B d x - ∫ 0 L 2 ρ * A 2 a ~ P ⋅ v ~ i P d z + ω ~ i D ⋅ T ~ * + v ~ i D ⋅ R *

(3)

식 (3)에서 ρ와 ρ^*는 축과 블레이드의 단위길이당 질량이며 $v ~ i Q$ , $a ~ i Q$ 는 축 중심 선상 임의 점 Q의 선형화된 편속도와 가속도이다. $w ~ i B$ 와 $α ~ B$ 는 축의 미소부분 B의 선형화된 편각속도와 각가속도이다. $I ⃡ B$ 는 축단면 미소부 B의 단위길이당 관성다이아딕, $v → i P$ , $a ~ P$ 는 블레이드 중립축 임의점 P의 선형화된 편속도와 가속도, $v → i D$ 와 $w ~ i D$ 는 D점의 선형화된 편속도와 편각속도들을 그리고 $T ~ *$ 와 $R ~ *$ 는 원판의 선형화된 관성 토크와 관성력을 각각 나타낸다.

축의 전단과 단면관성효과를 고려하되 블레이드 전단과 단면관성효과를 고려하지 않은 시스템의 탄성에너지는 다음과 같이 표현된다.

U = 12 ∫ 0 L 1 E A 1 ∂ s 1 ∂ x 2 d x + ∫ 0 L 1 E I 2 ∂ θ 2 ∂ x 2 d x + ∫ 0 L 1 E I 3 ∂ θ 3 ∂ x 2 d x + ∫ 0 L 1 K 2 G A 1 ∂ u 2 ∂ x - θ 3 2 d x + ∫ 0 L 1 K 3 G A 1 ∂ u 3 ∂ x + θ 2 2 d x + ∫ 0 L 1 G K ∂ θ 1 ∂ x 2 d x + ∫ 0 L 2 E A 2 ∂ s 2 ∂ z 2 d z + ∫ 0 L 2 E I 5 ∂ 2 u 4 ∂ z 2 2 d z + ∫ 0 L 2 E I 4 ∂ 2 u 5 ∂ z 2 2 d z

(4)

여기서, L₁, L₂는 축과 블레이드의 길이, E는 축과 블레이드의 탄성계수, I₂, I₃는 축의 면적관성모멘트, K₂, K₃는 축의 단면형상계수, G는 축의 전단계수이다. 또한, A₁, A₂는 축과 블레이드의 단면적, I₄, I₅는 블레이드의 면적 관성 모멘트, K는 축의 비틀림 계수이다.

$a^1$ 방향으로 ω의 각속도로 회전하는 축의 변형에 따른 임의 미소부 B의 각속도 $w → B$ 는 다음과 같이 나타낼 수 있다.

ω N→B = ω >N→A + ω >A→B = ω a ̯ 1 + θ ˙ 1 c 2 c 3 + θ ˙ 2 s 3 b ̯ 1 + - θ ˙ 1 c 2 s 3 + θ ˙ 2 c 3 b ̯ 2 + θ ˙ 1 s 2 + θ ˙ 3 b ̯ 3

(5)

또한, 점 D에서 축단면 각속도 $w → D$ 는 다음과 같이 나타낼 수 있다.

ω N → D = ω N → A + ω A → D = ω a ̯ 1 + θ 1 ¯ ˙ c 2 ¯ c 3 ¯ + θ 2 ¯ ˙ s 3 ¯ b ̯ 1 + θ 1 ¯ ˙ c 2 ¯ s 3 ¯ + θ 2 ¯ ˙ c 3 ¯ b ̯ 2 + - θ 1 ¯ ˙ s 2 ¯ + θ 3 ¯ ˙ b ̯ 3

(6)

여기서 ￣는 점 D에서의 각도, 각속도 등을 나타낸다. 축의 변형에 따른 중심선 상 임의점 Q의 속도 $v → Q$ , 점 D의 속도 $v → D$ , 점 R의 속도 $v → R$ , 그리고 블레이드 변형에 따른 그 중립축 상 임의점 P의 속도 $v → P$ 는 다음과 같이 구할 수 있다.

ν → Q = v → O + A v → Q + N ω → A × x → + u → 1 = u ˙ 1 a ̯ 1 + u ˙ 2 - ω u 3 a ̯ 2 + u ˙ 3 + ω u 2 a ̯ 3 ν → D = u ¯ ˙ 1 a ̯ 1 + u ¯ ˙ 2 - ω u ¯ 3 a ̯ 2 + u ¯ ˙ 3 + ω u ¯ 2 a ̯ 3 ν → R = ν → Q ¯ + N ω → D × r b ̯ 3 = u ¯ ˙ 1 a ̯ 1 + u ¯ ˙ 2 - ω u ¯ 3 a ̯ 2 + u ¯ ˙ 3 + ω u ¯ 2 a ̯ 3 - r ω c ¯ 2 c ¯ 3 + θ ¯ ˙ 1 C ¯ 2 c ¯ 3 + θ ¯ ˙ 2 s ¯ 3 b ̯ 2 + r - ω c ¯ 2 s ¯ 3 - θ ¯ ˙ 1 c ¯ 2 s 3 ¯ + θ ¯ ˙ 2 c ¯ 3 b ̯ 1 ν → P = ν → R + N ω → D × z → + u → 2 + B ν → P = u ¯ ˙ 1 a ̯ 1 + u ¯ ˙ 2 - ω u ¯ 3 a ̯ 2 + u ¯ ˙ 3 + ω u ¯ 2 a ̯ 3 + r + z - ω c ¯ 2 s ¯ 3 - θ ¯ ˙ 1 c ¯ 2 s ¯ 3 + θ ¯ ˙ 2 c ¯ 3 + u 6 - ω c ¯ 2 s ¯ 3 - θ ¯ ˙ 1 c ¯ 2 s ¯ 3 + θ ¯ ˙ 2 c ¯ 3 - u 5 ω s ¯ 2 + θ ¯ ˙ 1 s ¯ 2 + θ ¯ ˙ 3 + u ˙ 4 b ̯ 1 + {- r + z ω c ¯ 2 c ¯ 3 + θ ¯ ˙ 1 c ¯ 2 c ¯ 3 + θ ¯ ˙ 2 s ¯ 3 - u 6 ω c ¯ 2 c ¯ 3 + θ ¯ ˙ 1 c ¯ 2 c ¯ 3 + θ ¯ ˙ 2 s ¯ 3 + u 4 ω s ¯ 2 + θ ¯ ˙ 1 s ¯ 2 + θ ¯ ˙ 3 + u ˙ 5} b ̯ 2 + u 5 ω c ¯ 2 c ¯ 3 + θ ¯ ˙ 1 c ¯ 2 c ¯ 3 + θ ¯ ˙ c ¯ 2 c ¯ 3 - u 4 - ω c ¯ 2 s ¯ 3 - θ ¯ ˙ 1 c ¯ 2 s ¯ 3 + θ ¯ ˙ 2 c ¯ 3 + u ˙ 6 b ̯ 3

(7)

이상 식 (7)에 나타나는 u₁은 s₁, u₂, u₃간 근사화된 관계식이며(식 (6)), u₆는 s₂, u₄, u₅간 근사 관계식으로 다음과 같이 나타낼 수 있다.

u 1 ≅ s 1 - 12 ∫ 0 x ∂ u 2 ∂ σ 2 + ∂ u 3 ∂ σ 2 d σ u 6 ≅ s 2 - 12 ∫ 0 z ∂ u 4 ∂ σ 2 + ∂ u 5 ∂ σ 2 d σ

(8)

축 변형에 따른 임의 단면 미소부 B의 선형화된 각가속도 $α ~ B$ 와 점 D에서의 선형화된 각가속도 $α ~ D$ 는 다음과 같이 나타낼 수 있다.

α ∼ B = ω ˙ + θ ¨ 1 b ̯ 1 + - ω ˙ θ 3 - ω θ ˙ 3 + θ ¨ 2 b ̯ 2 + ω ˙ θ 2 + ω θ ˙ 2 + θ ¨ 3 b ̯ 3 α ∼ D = ω ˙ + θ ¯ ¨ 1 b ̯ 1 + - ω ˙ θ ¯ 3 - ω θ ¯ ˙ 3 + θ ¯ ¨ 2 b ̯ 2 + ω ˙ θ ¯ 2 + ω θ ¯ ˙ 2 + θ ¯ ¨ 3 b ̯ 3

(9)

축 중심선 상 임의점 Q에서의 선형화된 가속도 $α ~ Q$ , 블레이드의 중립축 상 임의점 P의 선형화된 가속도 $α ~ P$ , D점의 선형화된 가속도 $α → D$ 는 선형화된 속도 $v ~ Q$ , $v ~ D$ , $v ~ P$ 를 이용해 다음과 같이 나타낼 수 있다.

a ∼ Q = A d v ∼ Q d t + N ω ∼ A × v ∼ Q = s ¨ 1 a ̯ 1 + u ¨ 2 - ω ˙ u 3 - 2 ω u ˙ 3 - ω 2 u 2 a ̯ 2 + u ¨ 3 + ω ˙ u 2 + 2 ω u ˙ 2 - ω 2 u 3 a ̯ 3 a ∼ D = A d v ∼ D d t + N ω ∼ A × v ∼ D = s ¯ ¨ 1 a ̯ 1 + u ¯ ¨ 2 - ω u ¯ 3 ˙ - 2 ω u ¯ ˙ 3 - ω 2 u ¯ 2 a ̯ 2 + u ¯ ¨ 3 + ω ˙ u ¯ 2 + 2 ω u ¯ ¨ 2 - ω 2 u ¯ 3 a ̯ 3 a ∼ P = A d A → d t + N ω ∼ A × A → + B d B → d t + N ω ∼ B × B → = s ¯ ¨ 1 a ̯ 1 + u ¯ 2 ¨ - ω ˙ u ¯ 3 - 2 ω u ¯ 3 ˙ - ω 2 u ¯ 2 a ̯ 2 + u ¯ ¨ 3 + ω ˙ u ¯ 2 + 2 ω u ¯ ˙ 2 - ω 2 u ¯ 3 a ̯ 3 + [r + z - ω θ ¯ ˙ 3 - ω θ ¯ ˙ 3 + θ ¯ ¨ 2 + r + z ω 2 θ ¯ 2 + ω θ ¯ ˙ 3 + u ¨ 4 b ̯ 1 + - r + z ω ˙ + θ ¯ ¨ 1 - u 5 ω 2 - 2 s ˙ 2 ω - s 2 ω ˙ + u ¨ 5 b ̯ 2 + 2 u ˙ 5 ω + u 5 ω ˙ + s ¨ 2 - r + z ω 2 + 2 ω θ ¯ ˙ 1 - s 2 ω 2 b ̯ 3

(10)

여기서 $A →$ 와 $B →$ 는 벡터미분정리를 이용, $a ~ P$ 를 구하기 위해 $v ~ P$ 를 좌표계 $a^i$ 로 표기한 성분과 좌표계 $b^i$ 로 표기한 성분을 각각 나타낸 것으로 다음과 같다.

A → = s ¯ ˙ 1 a ̯ 1 + u ¯ ˙ 2 - ω u ¯ 3 a ̯ 2 + u ¯ ˙ 3 + ω u ¯ 2 a ̯ 3 B → = r + z - ω θ ¯ 3 + θ ¯ ˙ 2 + u ˙ 4 b ̯ 1 + - r + z ω + θ ¯ ˙ 1 - s 2 ω + u ˙ 5 b ̯ 2 + u 5 ω + s ˙ 2 b ̯ 3

(11)

축단면 미소부 B의 단위길이당 관성 다이아딕은 다음과 같이 나타낼 수 있다.

I ↔ B = I 11 b ̯ 1 b ̯ 1 + I 22 b ̯ 2 b ̯ 2 + I 33 b ̯ 3 b ̯ 3

(12)

$T ~ *$ , $R ~ *$ 는 원판에 작용하는 선형화된 관성 토크와 관성력이며 다음과 같이 나타낼 수 있다.

T ∼ * = - α ∼ D ⋅ I ↔ d - ω ∼ D × I ↔ d ⋅ ω ∼ D = - I d 1 ω ˙ + θ ¯ 1 ¨ b ̯ 1 + I d 3 - I d 1 ω 2 θ ¯ 2 + I d 3 - I d 1 ω θ ¯ ˙ 3 - I d 2 - ω ˙ θ ¯ 3 - ω θ ¯ 3 ˙ + θ ¯ 2 ¨ b ̯ 2 + - I d 1 - I d 2 ω 2 θ ¯ 3 + I d 1 - I d 2 ω θ ¯ 2 ˙ - I d 3 ω ˙ θ ¯ 2 + ω θ ¯ ˙ 2 + θ ¯ 3 ¨ b ̯ 3 R ∼ * = - M D a ∼ D = M D - s ¯ 1 ¨ a 1 ̯ u ¯ 2 ¨ - ω ˙ u ¯ 3 - 2 ω u ¯ 3 ˙ - ω 2 u ¯ 2 a ̯ 2 - u ¯ 3 ¨ + ω ˙ u ¯ 2 + 2 ω u ¯ 2 ˙ - ω 2 u ¯ 3 a ̯ 3

(13)

여기서 M_d는 원판의 질량이고 $I ⃡ d$ 는 원판의 관성다이아딕으로 다음과 같이 나타낼 수 있다.

I ↔ d = I 1 d b ̯ 1 b ̯ 1 + I 2 d b ̯ 2 b ̯ 2 + I 3 d b ̯ 3 b ̯ 3

(14)

식 (4) ~ (11)과 (13)을 식 (3)에 대입하면 운동방정식을 유도할 수 있으며 이를 식 (15)을 이용해 정리하면 다음과 같은 행렬 방정식 형태로 나타낼 수 있다.

M q ¨ i + G q ˙ i + K q i = F 1,2, ⋯, μ 1 + μ 2 + ⋯ μ 8 + μ 9

(15)

여기서

M α β (1) = ∫ 0 L 1 ϕ α i x ϕ β j x d x M α β (2) = ϕ α i L 1 ϕ β j L 1 M α β (3) = ∫ 0 L 2 ϕ α i z ϕ β j z d z M α β (4) = ∫ 0 L 2 ϕ α i L 1 ϕ β j z d z M α β (5) = ∫ 0 L 2 ϕ α i z ϕ β j L 1 d z M α β (6) = ∫ 0 L 2 ϕ α i L 1 ϕ β j L 1 d z M α β (7) = ∫ 0 L 2 r + z ϕ α i L 1 ϕ β j z d z M α β (8) = ∫ 0 L 2 r + z ϕ α i z ϕ β j L 1 d z M α β (9) = ∫ 0 L 2 r + z ϕ α i L 1 ϕ β j L 1 d z M α β (10) = ∫ 0 L 2 r + z 2 ϕ α i L 1 ϕ β j L 1 d z K α β (1) = ∫ 0 L 1 ϕ α i' x ϕ β i' x d x K α β (2) = ∫ 0 L 1 ϕ α i' x ϕ β i x d x K α β (3) = ∫ 0 L 1 ϕ α i x ϕ β i' x d x K α β (4) = ∫ 0 L 2 r L 2 - z ϕ α i' z ϕ β i' z d z K α β (5) = ∫ 0 L 2 12 L 2 2 - z 2 ϕ α i' z ϕ β i' z d z K α β (6) = ∫ 0 L 2 ϕ α i ″ z ϕ β i ″ z d z F α (1) = ∫ 0 L 2 r + z ϕ α i L 1 d z F α (2) = ∫ 0 L 1 ϕ α i x d x F α (3) = ∫ 0 L 2 r + z 2 ϕ α i L 1 d z F α (4) = ϕ α i L 1

(16)

여기서 M, G, K, F 행렬 등 각 행렬들의 부분행렬은 다음과 같다.

M = [m 11] 000 [m 51] 00 [m 18] 0 0 [m 22] 0 [m 42] 0000 [m 29] 00 [m 33] 000 [m 37] 00 0 [m 42] 0 [m 44] 0000 [m 49] [m 51] 000 [m 55] 00 [m 58] 0 00000 [m 66] 000 00 [m 73] 000 [m 77] 00 [m 81] 000 [m 85] 00 [m 88] 0 0 [m 92] 0 [m 94] 0000 [m 99] m 11 = ρ ∑ j = 1 μ 1 M 11 (1) + ∑ j = 1 μ 1 M 11 (4) + M D M 11 (2) m 15 = m 51 T = ρ ∑ j = 1 μ 5 M 15 (9) m 18 = m 81 T = ρ ∑ j = 1 μ 8 M 18 (5) m 22 = ρ ∑ j = 1 μ 2 M 22 (1) + ∑ j = 1 μ 2 M 22 (4) + M D M 22 (3) m 24 = m 42 T = - ρ ∑ j = 1 μ 4 M 24 (9) m 29 = m 92 T = ρ ∑ j = 1 μ 9 M 28 (5) m 33 = ρ ∑ j = 1 μ 3 M 33 (1) + ∑ j = 1 μ 3 M 33 (4) + M D M 33 (2) m 37 = m 73 T = ρ ∑ j = 1 μ 7 M 37 (5) m 44 = I 11 ∑ j = 1 μ 4 M 44 (1) + ρ ∑ j = 1 μ 4 M 44 (10) + I 1 d M 44 (2) m 49 = m 94 T = - ρ ∑ j = 1 μ 9 M 49 (7) m 55 = I 22 ∑ j = 1 μ 5 M 55 (1) + ρ ∑ j = 1 μ 5 M 55 (10) + I 2 d M 55 (2) m 58 = m 85 T = ρ ∑ j = 1 μ 8 M 58 (7) m 66 = I 33 ∑ j = 1 μ 6 M 66 (1) + I 3 d M 66 (2) m 77 = ρ ∑ j = 1 μ 7 M 77 (3) m 88 = ρ ∑ j = 1 μ 8 M 88 (3) m 99 = ρ ∑ j = 1 μ 9 M 99 (3)

(17)

G = 000000000 00 g 23 000 g 27 00 0 g 32 0 g 34 0000 g 39 00 g 43 000 g 47 00 000000000000000000 0 g 72 0 g 74 0000 g 79 000000000 00 g 93 000 g 97 00 g 23 = - g 32 T = - 2 ω ρ ∑ j = 1 μ 3 M 23 (1) + ∑ j = 1 μ 3 M 23 (4) + M D M 23 (2) g 27 = - g 72 T = - 2 ω ρ ∑ j = 1 μ 7 M 27 (5) g 34 = - g 43 T = - 2 ω ρ ∑ j = 1 μ 4 M 34 (9) g 39 = - g 93 T = 2 ω ρ ∑ j = 1 μ 9 M 39 (5) g 47 = - g 74 T = 2 ω ρ ∑ j = 1 μ 7 M 47 (7) g 79 = - g 97 T = 2 ω ρ ∑ j = 1 μ 9 M 79 (3)

(18)

K = k 11 00000000 0 k 22 k 23 k 24 0 k 26 k 27 0 k 29 0 k 32 k 33 k 34 k 35 0 k 37 0 k 39 0 k 42 k 43 k 44 00 k 47 00 00 k 53 0 k 55 k 56 000 0 k 62 00 k 65 k 66 000 0 k 72 k 73 000 k 77 0 k 79 0000 k 85 k 86 000 0 k 92 k 93 000 k 97 0 k 99 k 11 = E A 1 ∑ j = 1 μ 1 K 11 (1) k 22 = - ω 2 m 22 + K 2 G A 1 ∑ j = 1 μ 2 K 22 (1) k 23 = k 32 T = ω ˙ 2 ω g 23 k 24 = k 42 T = - ω 2 m 24 k 26 = k 62 T = - K 2 G A 1 ∑ j = 1 μ 6 K 26 (2) k 27 = k 72 T = ω ˙ 2 ω g 27 k 29 = k 92 T = - ω 2 ρ ∑ j = 1 μ 9 K 29 (5) k 33 = - ω 2 m 33 + K 3 G A 1 ∑ j = 1 μ 3 K 33 (1) k 34 = k 34 T = ω ˙ 2 ω g 34 k 35 = k 53 T = K 3 G A 1 ∑ j = 1 μ 5 K 35 (2) k 37 = k 73 T = - ω 2 ρ ∑ j = 1 μ 7 K 37 (5) k 39 = k 93 T = ω ˙ 2 ω g 39 k 44 = G K ∑ j = 1 μ 4 K 44 (1) k 47 = ω ˙ ω g 47 k 55 = E I 1 ∑ j = 1 μ 5 K 55 (1) + K 3 G A 1 ∑ j = 1 μ 5 M 55 (1) + ω 2 ρ ∑ j = 1 μ 5 M 55 (10) + I 1 d - I 3 d M 55 (2) k 56 = k 65 T = ω ˙ I 11 - I 22 ∑ j = 1 μ 6 M 56 (1) + ω ˙ I 1 d - I 2 d M 56 (2) k 58 = k 85 T = ω 2 m 58 k 66 = E I 3 ∑ j = 1 μ 6 K 66 (1) + K 2 G A 1 ∑ j = 1 μ 6 M 66 (1) + ω 2 I 1 d - I 2 d M 66 (2) k 68 = - k 86 T = ω ˙ ρ ∑ j = 1 μ 8 M 68 (2) k 77 = - ω 2 m 77 + E A 2 ∑ j = 1 μ 7 K 77 (1) k 79 = - k 97 T = ω ˙ 2 ω g 79 k 88 = ω 2 ρ ∑ j = 1 μ 8 K 88 (3) + ∑ j = 1 μ 8 K 88 (4) + E I 5 ∑ j = 1 μ 8 K 88 (5) k 99 = - ω 2 m 99 - ρ ∑ j = 1 μ 9 K 99 (4) + ∑ j = 1 μ 9 K 99 (5) + E I 4 ∑ j = 1 μ 9 K 99 (6)

(19)

F = 0 ω ˙ ρ F 2 (1) ω 2 ρ F 3 (1) ω ˙ I 11 F 4 (2) + ρ F 4 (3) + I 1 d F 4 (4) 00 ω 2 ρ F 7 (1) 0 ω ˙ ρ F 9 (1)

(20)

3. 수치해석

3.1 ANSYS 고유진동수 해석결과

Table 1은 이 연구에서 사용된 모델의 데이터이다. 축, 원판과 블레이드는 각 1개만 모델링하였다. 축, 원판과 블레이드 유한요소는 1차원 beam 188요소를 사용하였으며 축, 원판과 블레이드에 각각 10개, 2개 그리고 5개의 요소 사용하였다. 원판은 축의 자유단 끝에 두께 1 mm로 모델링하였다. 원판 중심노드와 블레이드 길이방향 첫 번째 노드 사이에 강체 경계조건을 적용하였다. Table 2에는 ANSYS⁽⁹⁾를 이용해 구한 축 굽힘모드 1차 ~ 4차 고유진동수, 인장모드 1차 고유진동수, 비틀림모드의 1차 고유진동수 그리고 블레이드 굽힘모드의 1, 2차 고유진동수의 결과가 나타나 있다.

Table 1

Numerical data of a shaft-disk-blade system

Table 2

Numerical frequencies obtained with ANSYS(unit: Hz)

3.2 축의 굽힘 변위 및 각도에 대한 모드 함수

이 절은 앞에서 유도된 운동방정식을 수치적으로 해석하기 위한 축의 굽힘 변위 및 각도와 관련된 세 종류 모드함수들을 소개하고자 한다. 첫 번째 경우(case 1)는 외팔보의 굽힘 변위 및 각도에 대한 일반화된 모드함수들을 축의 굽힘 변위 및 각도에 대한 모드함수로 사용하는 경우로 다음과 같다.

ϕ 2 i (x) = c o s h λ 2 i x L 1 - c o s λ 2 i x L 1 - c o s h λ 2 i + c o s λ 2 i s i n h λ 2 i + s i n λ 2 i s i n h λ 2 i x L 1 - s i n λ 2 i x L 1

(21)

c o s λ 2 i ⋅ c o s h λ 2 i + 1 = 0

(22)

여기서 식 (21)은 모드함수, 식 (22)는 특성방정식을 나타내는데 ϕ₂_i(x)는 ϕ₃_i(x)와 동일하다.

두 번째 경우(case 2)는 외팔보의 자유단의 끝에 집중질량을 가진 보의 굽힘 변위 및 각도에 대한 모드함수들을 축의 굽힘 변위 및 각도에 대한 모드함수로 사용한 경우로 다음 같이 나타낼 수 있다⁽¹⁰⁾.

ϕ 2 i (x) = c o s λ 2 i x - c o s h λ 2 i x - c o s λ 2 i L 1 + c o s h λ 2 i L 1 s i n λ 2 i L 1 + s i n h λ 2 i L 1 (s i n λ 2 i x - s i n h λ 2 i x)

(23)

1 + 1 c o s λ 2 i L 1 c o s h λ 2 i L 1 - M D ρ A L 1 λ 2 i L 1 t a n λ 2 i L 1 - t a n h λ 2 i L 1 = 0

(24)

여기서 식 (23)은 모드함수, 식 (24)는 특성방정식을 나타낸다.

세 번째 경우(case 3)는 외팔 보 자유단 끝에 집중질량과 집중질량의 질량 관성모멘트 효과를 동시에 고려한 굽힘 변위 및 각도 관련 모드함수들을 축의 굽힘 변위 및 각도를 위한 모드함수들로 사용하는 경우로 다음과 같다⁽¹¹⁾.

ϕ 2 i (x) = s i n λ 2 i x - s i n h λ 2 i x + c o s λ 2 i L 1 + c o s λ 2 i L 1 s i n λ 2 i L 1 + s i n h λ 2 i L 1 c o s λ 2 i x - c o s h λ 2 i x

(25)

1 + c o s λ 2 i L 1 c o s h λ 2 i L 1 - J D J S λ 2 i L 1 3 s i n h λ 2 i L 1 c o s λ 2 i L 1 + s i n λ 2 i L 1 c o s h λ 2 i L 1 = 0

(26)

여기서 J_D는 원판과 블레이드의 질량관성모멘트의 합이고 J_S는 축의 질량관성모멘트이다. 이외 블레이드의 굽힘 변위를 위한 모드함수는 식 (21), 축의 비틀림은 외팔보의 자유단 끝에 집중질량 및 그 집중질량의 질량관성모멘트를 고려한 모드함수로 사용했으며 축의 인장은 외팔보 자유단 끝에 집중질량이 있는 모드함수로 사용하였다⁽¹²^,¹³⁾. 수치해석 결과와 ANSYS 결과 간의 오차를 식 (27)과 같이 정의하였다.

P r e s e n t - A N S Y S A N S Y S × 100 = P % P > 0, E r r o r % = P p ≤ 0, E r r o r % = 0

(27)

여기서 P는 수치해석 결과와 ANSYS 결과 간의 상대오차를 나타낸다. 그런데 여기서 P < 0이면 수치해석 결과가 ANSYS의 결과보다 더 낮은 고유진동수 값을 갖는다는 것을 나타내는 것이다. 유한요소법에서 요소 수를 증가시키거나 가상모드법에서 사용모드 수를 증가시키면 해석결과가 점점 더 작은 값으로 수렴해 간다는 일관된 수렴특성은 잘 알려져 있다. 따라서 P < 0이면 사용된 모델의 수치해석 결과가 ANSYS 결과보다 더 정확한 값에 수렴했다는 것을 보여주는 것이므로 error를 0으로 정의한 것이다.

3.3 축에 대한 모드함수 별 수치해석 결과와 ANSYS 결과 간 오차

제안된 축-원판-블레이드 연성모델을 앞서 소개된 모드함수(case 1 ~ case 3)별로 모드 수(모드함수의 수)를 축의 인장, 굽힘, 비틀림 방향 그리고 블레이드 인장, 굽힘 방향별로 각각 동일한 수로 증가할 때에 고유진동수 수치해석 결과⁽¹⁴⁾를 구하고 수치해석 결과를 식 (27)로 ANSYS 고유진동수 해석결과와의 오차를 Figs. 2 ~ 4에 나타내었다.

Fig. 2

Error versus the number of mode (case 1)

Fig. 3

Error versus the number of mode (case 2)

Fig. 4

Error versus the number of mode (case 3)

Fig. 2는 모드 수 증가에도 불구하고 축의 3, 4차 굽힘 고유진동수와 블레이드 1차 굽힘 고유진동수 오차가 1 % 이하로 감소하지 않았음을 알 수 있다. Fig. 3은 Fig. 2와 달리 축 4차 굽힘 고유진동수와 블레이드 1차 굽힘 고유진동수가 모드 수 증가에 따라 1 % 이하로 감소하나 축 3차 굽힘 고유진동수는 1 % 이하로 감소하지 않음을 알 수 있다.

Fig. 4에는 모드 수가 3개일 때부터 축 1~4차 굽힘 고유진동수, 축 1차 인장 고유진동수, 축 1차 비틀림 고유진동수 그리고 블레이드 1, 2차 굽힘 고유진동수들이 모두 오차 0.11 % 이하임을 확인할 수 있다. Fig. 4에서 모드 수가 6개일 때를 충분히 수렴한 결과라 보고 이를 Table 3에 나타내었다.

Table 3

Natural frequencies obtained with the present method and ANSYS

3.4 모드함수 별 수렴 특성

Table 3의 present 결과를 기준으로 하여 3.2장에서 소개한 case 1과 case 2에 대한 오차들을 식 (27)를 이용하여 Figs. 5, 6과 같이 나타내었다. Fig. 5에는 Fig. 2와 달리 모드 수가 증가함에도 축 3, 4차 굽힘 고유진동수, 축 1차 인장 고유진동수 그리고 블레이드 1차 굽힘 고유진동수가 1 % 이하로 감소하지 않았음을 알 수 있다. Fig. 6에서는 Fig. 3과 달리 모드 수가 증가해도 축 3, 4차 굽힘 고유진동수와 블레이드의 1차 굽힘 고유진동수가 1 % 이하로 감소하지 않음을 알 수 있다.

Fig. 5

Error versus the number of mode (case 1) employing the reference values in Table 3

Fig. 6

Error versus the number of mode (case 2) employing the reference values in Table 3

4. 결 론

축의 굽힘 방향 변위 및 각도에 대해 세 종류의 외팔보의 굽힘 방향 모드함수로 사용한 결과, 일반화된 외팔보의 굽힘 방향 모드함수(case 1)와 외팔보 자유단 끝에 집중질량만 고려한 외팔 보 굽힘 방향 모드함수(case 2)를 축의 굽힘 방향 변위 및 각도에 대한 모드함수로 사용하면 모드 수가 증가하여도 특정 모드에 대한 고유진동수들의 오차를 감소시킬 수 없음을 확인할 수 있었다. 이에 반해 외팔보의 자유단 끝에 집중질량 및 그 집중질량의 질량관성모멘트를 동시 고려한 모드함수를 사용하면 모드 수의 증가에 따라 모든 모드의 빠른 수렴성을 달성할 수 있음을 확인할 수 있었으며 적은 모드 수만 사용하여도 ANSYS와 동등한 정확도의 결과를 얻을 수 있음을 확인하였다.

이 연구에서 유도된 축-원판-블레이드 연성계의 모델은 향후 이의 설계를 위한 진동 및 과도 특성 해석에 매우 유용하게 사용할 수 있을 것이다. 유한요소모델에 비해 소수의 일반좌표들만을 사용하는 이러한 모드좌표 기반 해석모델은 축-원판-블레이드 연성계의 기본설계 및 평가에 효율적으로 이용될 수 있을 것이다.

Acknowledgments

이 연구는 산업통상자원부(MOTIE)와 한국에너지기술평가원(KETEP)의 지원을 받아 수행한 연구 과제입니다(No. 2018201010636A).

References

Lee, C. W., and Chun, S. B., (1996), Vibration Analysis of Shaft-bladed Disk System by Using Substructure Synthesis and Assumed Modes Method, Journal of Sound and Vibration, 189(5), p857-608. [https://doi.org/10.1006/jsvi.1996.0038]
Al-Bedoor, B. O., (1999), Dynamic Model of Coupled Shaft Torsional and Blade Bending Deformation in Rotors, Computer Methods in Applied Mechanics and Engineering, 169(1-2), p177-190. [https://doi.org/10.1016/s0045-7825(98)00184-4]
Al-Bedoor, B. O., (2007), Natural Frequencies of Coupled Blade-bending and Shaft-torsional Vibration, Shock and Vibration, 14, p65-80. [https://doi.org/10.1155/2007/506165]
Yang, C. H., and Huang, S. C., (2007), The Influence of Disk’s Flexibility on Coupling Vibration of Shaft-disk-blades Systems, Journal of Sound and Vibration, 301(1-2), p1-17. [https://doi.org/10.1016/j.jsv.2006.01.053]
Turban, Ö., and Bulut, G., (2006), Linearly Coupled Shaft-torsional and Blade-bending Vibrations in Multi-stage Rotor-blade Systems, Journal of Sound and Vibration, 296(1-2), p292-318.
Yoo, H. H., Ryan, R. R., and Scott, R. A., (1995), Dynamics of Flexible Beam Undergoing Overall Motions, Journal of Sound and Vibration, 181(2), p261-278. [https://doi.org/10.1006/jsvi.1995.0139]
Kane, T., and Levinson, D., (1985), Dynamics: Theory and Applications, New York, McGraw-Hill Co.
Inman, D. J., (2014), Engineering Vibrations, 4th Ed., Boston, Pearson.
ANSYS, (1998), User’s Manual Structural Analysis Guide, Pennsylvania, ANSYS Inc.
Rao, S. S., (2007), Vibration of Continuous System, Wiley, New York.
Braun, S. G., Ewins, D. J., and Rao, S. S., (2001), Encyclopedia of Vibration, Academic Press, New York.
R. D. Blevins, (1979), Formulas for Natural Frequency and Mode Shape, Van Nostrand Rein-hold Ltd., New York.
Meirovitch, L., (1967), Analytical Methods in Vibrations, Macmillan Publishing Co., Inc., New York. [https://doi.org/10.1049/ep.1967.0353]
Yoo, H. H., and Shin, S. H., (1998), Vibration Analysis of Rotating Cantilever Beam, Journal of Sound and Vibration, 212(5), p807-828. [https://doi.org/10.1006/jsvi.1997.1469]
Shin, S. H., and Yoo, H. H., (1996), Modeling and Vibration Analysis of Rotating Cantilever Beams Considering Shear and Rotary Inertia Effects, Transactions of the Korean Society for Noise and Vibration Engineering, 6(2), p179-185.

Hyung Hee Kim received his B.S. degree in the Department of Mechanical Engineering in Youngnam University in 2012. He is working as a M.S. candidate in the Department of Mechanical Convergence Engineering in Hanyang University, Seoul, Korea. His research interests include structural vibration and dynamics.

Hong Hee Yoo received his B.S. and M.S. degrees in the Department of Mechanical Design in Seoul National University in 1980 and 1982. He received his Ph.D. degree in the Department of Mechanical Engineering and Applied Mechanics in the University of Michigan at Ann Arbor in 1989. He is a professor in the Department of Mechanical Engineering at Hanyang University, Seoul, Korea. His research interests include multi-body dynamics, structural vibration and statistical uncertainty analysis in Mechanics.

Description	Value
Material property	E = 200 GPa, ρ = 7850 kg/m³
Poisson ratio	ν = 0.25
Blade length	0.125 m
Blade cross section area	2.54 × 0.1 cm²
Blade mass per unit length	0.2 kg/m
Disk radius	0.05 m
Disk moment of inertia	0.012 331 25 kg·m²
Shaft length	0.5 m
Shaft diameter	1 cm
Shaft torsional stiffness	157.08 N·m/rad

Shaft 1st bending natural frequency (N.F)	2.4374
Shaft 2nd bending N.F	2.4378
Shaft 3rd bending N.F	53.776
Shaft 4th bending N.F	58.610
Shaft 1st torsional N.F	17.713
Shaft 1st longitudinal N.F	282.23
Blade 1st bending N.F	47.727
Blade 2nd bending N.F	328.43

	ANSYS [Hz]	Present [Hz]	Error (%)
Shaft 1st bending N.F	2.4374	2.4402	0.11
Shaft 2nd bending N.F	2.4378	2.4406	0.11
Shaft 3rd bending N.F	53.776	53.502	0
Shaft 4th bending N.F	58.610	58.443	0
Shaft 1st torsional N.F	17.713	17.714	0.01
Shaft 1st longitudinal N.F	282.23	279.81	0
Blade 1st bending N.F	47.727	47.614	0
Blade 2nd bending N.F	328.43	328.06	0