Dynamic Analysis of Tie-rod-fastened Rotor Considering Elastoplastic Deformation

탄소성 변형을 고려한 타이로드 고정 회전체의 동역학 해석

Abstract

This study conducts numerical modeling and eigen-analysis of a rod-fastened rotor, which is mainly used in aircraft gas turbine engines in which multiple disks are in contact through curvic coupling. Nayak's theory is adopted to calculate surface parameters measured from the tooth profile of the curvic coupling gear. Surface parameters are important design parameters for predicting the stiffness between contact surfaces. Based on the calculated surface parameters, elastoplastic contact analysis is performed according to the interference between two surfaces based on the Greenwood-Williamson model. The equivalent bending stiffness is predicted based on the shape and elastoplastic contact stiffness of the curvic coupling. An equation of motion of the rod-fastened rotor, including the bending stiffness of the curvic coupling, is developed. Methods for applying the bending stiffness of a curvic coupling to the equation of motion and for modeling the equation of motion of a rotor that includes both inner and outer rotors are introduced. Rotordynamic analysis is performed through one-dimensional finite element analysis, and each element is modeled based on Timoshenko beam theory. Changes in bending stiffness and the resultant critical speed change in accordance with the rod fastening force are predicted, and the corresponding mode shapes are analyzed.

Nomenclature

An : Normal contact area of rough surface (m²) (거친 표면 수직 접촉 면적)

a_r : Reduction rate of area (면적 감소율)

a_p : Area of one side of tooth (치 한 면의 면적)

D_zero : The number of zero values per unit length (#/m) (단위 길이당 영점의 수)

D_extrema : The number of peak values per unit length (#/m) (단위 길이당 극 값의 수)

d : Separation of rough surface (간극)

[C] : Damping matrix(감쇠행렬)

E' : Equivalent Young’s modulus (Pa) (상당 탄성계수)

F_BW : Contact load of rough surface (N) (거친 표면 접촉력)

G : Gyroscopic moment (자이로스코픽 항)

H : Material hardness (Pa) (재료 강도)

I_d : Transverse mass of inertia (질량 관성 모멘트)

I_P : Polar moment of inertia (극 관성 모멘트)

[K] : Stiffness matrix (강성행렬)

K_z : Contact stiffness at z direction

k : Contact stiffness per unit area (N/m) (단위 면적에서의 접촉 강성)

[M] : Mass matrix (질량행렬)

m : Moment of power spectral density (파워 스펙트럼 밀도의 모멘트)

n : Number of tooth of coupling (커플링 치 개수)

R : Average asperity radius (m) (평균 돌기 반경)

Δ : Displacement (변위)

# : Number (개수)

α : Contact angle (°) (접촉각)

δ : Hertzian contact interference (헤르츠 접촉 간섭량)

ε : Interference of single asperity (m) (간섭량)

ε_c : Critical interference of single asperity (m) (임계 간섭량)

η : Spatial density of asperity (#/m²) (돌기의 공간 밀도)

λ : Eigenvalue (고유 값)

μ : Friction coefficient (마찰 계수)

σ : Stress (응력)

ν : Poisson’s ratio (푸아송 비)

Superscripts

B : Bearing (베어링)

d : Disk (디스크)

IB : Inter bearing (연결 베어링)

S : System

1. 서론

항공기 가스터빈 엔진은 무게, 가공, 조립, 제작 및 냉각 설계 용이성의 이유로 Rod-Fastened Rotor(RFR) 구조를 채택한다. RFR은 고정 방식에 따라 피치 원 직경의 여러 타이볼트로 고정되는 타이볼트 고정 회전체가 있고, 중앙 타이로드 고정 회전체가 있다. 그 중 중앙 타이로드 고정 회전체는 압축 너트를 통해 체결이 되며, 회전체의 디스크들은 물림 커플링(Claw Coupling)로 연결된다. 물림 커플링은 기어 치형에 따라 다양한 종류가 있다. 특히 커빅 커플링은 가공 및 교체가 간편하고 동력 전달 효율이 뛰어나 널리 채택된다.

하지만 이러한 조립 구조를 갖는 RFR은 비 연속적 구조로 인해 회전체 시스템의 굽힘 강성 및 동특성 예측에 어려움이 있다. 가공에 의해 형성된 금속의 표면이 육안으로 관측될 때에는 매끄러운 것처럼 보이나, 어떤 표면이든지 확대하면 산과 골로 형성된 표면 거칠기를 가진다. 연속체 회전체에서의 간단한 강성 계산식과 달리, RFR은 거칠기를 가지는 표면 간의 접촉으로 인해 복잡한 강성 예측 방법을 필요로 한다. 그러므로 커플링 접촉면에서 강성의 영향으로 전체 시스템의 강성과 동특성이 변할 수 있기 때문에 기존 전달 행렬법과 유한 요소법으로는 간단히 모델링할 수 없다.

접촉 표면의 프로파일은 접촉면 강성에 큰 영향을 미치므로 이를 고려한 해석 모델의 개발이 필요하다. RFR의 체결력과 접촉 강성을 고려한 선행 연구는 다음과 같다. Yuan[1] 등은 3차원 유한요소 접촉 모델을 이용하여 커빅커플링 간의 예압력을 정확하게 적용하는 방법을 연구하였다. 하지만 해당 연구에서는 접촉 인터페이스에서의 거친 표면에 대한 영향을 고려할 수 없었다. 또한 Yuan[2]은 GW 탄성 접촉 모델을 이용하여 커빅 커플링의 형상을 고려한 접촉 및 굽힘 강성을 예측하였고, 표면 거칠기 R_a값에 따른 회전체의 위험속도를 계산하였다. 하지만 굽힘 강성을 회전체 시스템에 어떻게 적용했는지 명확히 제시하지 않았다. Peng[3] 등은 디스크들이 접촉하는 타이 볼트 고정 회전체에서 디스크의 거칠기를 고려한 탄-소성 접촉해석 모델을 이용하여 접촉 및 굽힘 강성을 분석하고 이에 따른 위험속도를 예측하였다. Yu[4] 등은 GW 탄-소성 접촉 모델을 이용하여 커빅커플링 형상을 고려했을 뿐 아니라 고착(Stick)영역과 미끄러짐(Slip)영역을 구분하여 접촉 강성을 해석하였다. 해당 논문은 커빅커플링의 거칠기를 고려한 정확한 굽힘 강성 해석법을 제시하였지만, 커빅커플링 단품 자체에 대한 해석만을 하였고, 이를 이용한 회전체동역학 해석은 진행하지 않았다. Liu[5] 등은 커빅 커플링이 포함된 회전체동역학 모델을 제시하였는데, 체결력에 따른 축 방향 변형을 고려하여 커빅 커플링의 굽힘 강성을 계산하였다. 해당 연구에서는 계산된 굽힘 강성을 회전체의 시스템에 적용하는 방법을 제시하였다.

일반적인 회전체와는 달리 중앙 타이로드 체결 회전체는 중앙에 긴 로드(rod)와 바깥으로 여러 개의 디스크들이 압축 너트로 체결되는 구조를 가진다. 로드와 디스크 간의 간극이 존재하기 때문에 휘돌림(whirling) 궤적 크기가 서로 달라 마찰 충격 등이 발생할 수 있으므로 이를 단일 로터로 해석하면 여러 비선형 효과를 확인하지 못할 것이다. 따라서 로드와 디스크를 각각 안쪽 회전체(inner rotor)와 바깥쪽 회전체(outer rotor)로 모델링 할 수 있다. Sun[6]은 두 축의 마찰 충격과 관련된 비선형 동적 거동을 예측하기 위해 복 축 회전체의 운동방정식을 개발하였다. 본 연구에서는 Sun[6]이 도입한 복 축 회전체의 운동 방정식을 사용하였다. 대부분의 복 축 구조 회전체의 경우 간 베어링 단순 접촉 및 열 박음 구조를 가지는데, Hou[7]의 연구에서는 회전체 구조물 간에 연결되는 경우 회전체 구조 사이의 접촉강성을 계산하여 선형 베어링 강성으로 고려할 수 있다고 발표하였다.

본 연구에서는 Yu [4]의 연구에서 정립된 커빅커플링 접촉면에서의 탄-소성 접촉 및 이로 인한 굽힘강성해석을 수행하는 프로세스를 이용하였다. 커빅커플링 접촉면에서 고착 영역과 미끄러짐 영역을 구분하여 접촉 강성을 계산하였다. 실제 커빅커플링 표면의 거칠기를 측정하여 GW 모델에 필요한 표면 거칠기 변수를 계산하였고, 커빅커플링의 형상을 수치적으로 고려하는 방법을 제시하였다. 체결력에 따른 접촉 및 회전체 굽힘 강성 변화를 분석하였고, 고속회전에 따른 커빅커플링 원심 팽창으로 인한 겉보기 접촉 면적이 굽힘 강성에 미치는 영향을 분석하였다. RFR의 로드와 디스크 부분을 각각 안쪽 회전체와 바깥쪽 회전체의 복 축 회전체 시스템으로 고려하였고, Timoshenko beam모델을 이용하여 커빅커플링의 굽힘 강성이 회전체동역학 시스템에 미치는 영향을 분석하였다. 또한, 커빅커플링의 굽힘 강성을 회전체 시스템에 적용하는 방법과, 복 축 구조를 고려하는 방법에 대해 설명한 후 위험속도(Critical Speed)와 그에 따른 모드 형상(Mode shape)을 계산하였다.

2. 이론적 배경

2.1. 커플링 간의 접촉

커플링 접촉 해석을 위해 두 커플링 표면 사이의 거칠기 분포가 정규 분포를 따르는 것으로 가정하였다. 표면 거칠기가 정규분포를 따르면 GW모델을 이용해 접촉력과 접촉 강성을 구할 수 있으며, 이를 구하는데 필요한 표면 매개변수(parameter)는 Nayak[8,9]의 Random process 이론을 통해 구할 수 있다. GW 모델에서는 Fig. 1(a)와 같은 거친 표면을 Fig. 1(b)와 같이 총 N개의 동일한 곡률반경 R을 가지는 구형 돌기로 단순하게 모델링 된다.

Fig. 1.(a) Real rough surface, (b) Simplified rough surface.

2.1.1. 거친 표면의 Random Process

본 연구에서는 Random process이론을 기반으로 표면을 2차원 등방성 정규분포로 모델링하는 Nayak의 연구 결과[8,9]가 사용되었다. GW 모델을 이용한 접촉해석을 수행하기 위해 표면의 Power spectral density(PSD)의 모멘트로 세 가지 필수 표면 매개변수인 평균 돌기 반경(R), 돌기의 공간 밀도(η), 돌기 높이의 표준 편차(ξ_a)를 구할 수 있다. Longuet과 Higgins[10]의 연구결과를 바탕으로 PSD 모멘트와 D_zero, D_extrema는 다음과 같은 관계를 가진다[8,9].

\(\begin{align}D_{\text {zero }}=\frac{1}{\pi}\left(\frac{m_{2}}{m_{0}}\right)^{1 / 2}\end{align}\)       (1)

\(\begin{align}D_{\text {extrema }}=\frac{1}{\pi}\left(\frac{m_{4}}{m_{2}}\right)^{1 / 2}\end{align}\)       (2)

표면 매개변수는 PSD 모멘트로 다음과 같이 구할 수 있다.

\(\begin{align}\xi_{a}=\left(1-\frac{0.8968}{\beta}\right) \xi \text { where }\beta=\left(\frac{D_{\text {extrema }}}{D_{\text {zero }}}\right)^{2}\end{align}\)       (3)

\(\begin{align}R=\frac{3}{8} \sqrt{\frac{\pi}{m_{4}}}\end{align}\)       (4)

\(\begin{align}\eta=\frac{1}{6 \sqrt{3} \pi} \frac{m_{4}}{m_{2}}\end{align}\)       (5)

2.1.2. Greenwood Williamson 탄성 접촉 모델

GW 모델[11]에서 거친 표면은 N개의 돌기로 구성되며, 각각의 돌기들은 기준면으로부터 높이(z)를 갖는다. 모든 돌기들은 동일한 곡률 반경 R을 가진다고 가정하며, 돌기 높이 분포는 정규분포를 따른다. GW 모델에서의 접촉은 두 평면 간의 수많은 구와 강체 평면 간의 접촉이므로 각각의 돌기들 사이에 헤르츠(Hertz) 접촉 거동을 가정한다면 z와 dz 사이에 존재하는 돌기들에 의한 접촉력은 다음과 같이 표현된다.

\(\begin{align}\boldsymbol{d} F_{G W}=\frac{4}{3} E^{\prime} R^{\frac{1}{2}} \delta^{\frac{3}{2}} d n=\frac{4}{3} N E^{\prime} R^{\frac{1}{2}}(z-d)^{\frac{3}{2}} \varphi(z) d z\end{align}\)       (6)

GW모델에서는 무차원 정규분포를 이용한 접촉력 및 단위 면적당 수직 접촉 강성을 각각 아래와 같이 나타낸다.

\(\begin{align}F_{G W}=\frac{4}{3} \eta A_{G W} E^{\prime} R^{\frac{1}{2}} \sigma^{\frac{3}{2}} \int_{h}^{\infty}(s-h)^{\frac{3}{2}} \varphi(s) d s\end{align}\)       (7)

\(\begin{align}k_{G W}=2 \eta E^{\prime} R^{\frac{1}{2}} \sigma^{\frac{1}{2}} \int_{s}^{\infty}(s-h) \varphi(s) d s\end{align}\)       (8)

2.1.3. 탄-소성 접촉해석

하중에 의한 간섭량 ε이 임계 간섭량을 초과하면 돌기 내의 응력이 재료의 항복응력을 초과하여 소성 변형이 발생한다. 소성 영역의 영향까지 고려하기 위해 많은 연구가 수행되었는데, 유한 요소 해석을 통해 탄-소성 변형을 통합한 Kogut과 Etsion(KE 모델)[12]을 통해 소성 영역을 고려하였으며, Jackson과 Green (JG모델)[13]은 KE모델보다 큰 간섭량을 가질 때 재료 특성을 고려할 수 있다. 두 모델 모두 헤르츠 접촉 모델과 간섭량을 기반으로 한 회귀분석식을 제안하였다. 해당 탄-소성 접촉 모델에 관한 내용을 Table 1에 정리하였다.

Table 1. Elastic-plastic contact force equation

이와 같이 임의의 간섭량에서 접촉 하중은 헤르츠 접촉식과 간섭량에 대한 회귀분석식으로 표현이 가능하다. 각 간섭량 구간에 따른 접촉력 모델과 돌기 높이의 확률 분포를 고려하면, 미소 표면에서의 접촉력 및 접촉 강성을 다음과 같이 표현할 수 있다.

\(\begin{align}F_{E-P}=\frac{4}{3} \eta A_{n} R^{\frac{1}{2}} E^{\prime} \varepsilon_{c}^{\frac{3}{2}}\left[\begin{array}{c}\frac{\varepsilon}{\varepsilon_{c}}-1.5 \int_{h}^{h+\frac{\varepsilon}{\varepsilon_{c}}}(s-h)^{1.5} \phi^{*}(s) d s \\ +1.03 \frac{\varepsilon}{\varepsilon_{c}}-1.425 \int_{h+\frac{\varepsilon}{\varepsilon_{c}}}^{h+\frac{6 \varepsilon}{\varepsilon_{c}}}(s-h)^{1.425} \phi^{*}(s) d s \\ +1.40 \frac{\varepsilon}{\varepsilon_{c}}-1.263 \int_{h+\frac{6 \varepsilon}{\varepsilon_{c}}}^{h+\frac{110 \varepsilon}{\varepsilon_{c}}}(s-h)^{1.263} \phi^{*}(s) d s \\ +4.6 C(v) \int_{h+\frac{110 \varepsilon}{\varepsilon_{c}}}^{\infty}(s-h) \phi^{*}(s) d s\end{array}\right]\end{align}\)       (9)

\(\begin{align}=\frac{4}{3} \eta E^{\prime} R^{\frac{1}{2}} \varepsilon_{c}^{\frac{1}{2}}\left[\begin{array}{c}1.5 \frac{\varepsilon}{\varepsilon_{c}}-0.5 \int_{h}^{h+\frac{\varepsilon}{\varepsilon_{c}}}(s-h)^{0.5} \phi^{*}(s) d s \\ +1.03 \times 1.425 \frac{\varepsilon}{\varepsilon_{c}}-0.425 \int_{h+\frac{\varepsilon}{\varepsilon_{c}}}^{h+\frac{6 \varepsilon}{\varepsilon_{c}}}(s-h)^{0.425} \phi^{*}(s) d s \\ +1.40 \times 1.263 \frac{\varepsilon}{\varepsilon_{c}}-0.263 \int_{h+\frac{6 \varepsilon}{\varepsilon_{c}}}^{h+\frac{110 \varepsilon}{\varepsilon_{c}}}(s-h)^{0.263} \phi^{*}(s) d s \\ +4.6 C(v) \int_{h+\frac{110 \varepsilon}{\varepsilon_{c}}}^{\infty}(s-h) \phi^{*}(s) d s\end{array}\right]\end{align}\)       (10)

여기서 임계 간섭량(ε_c)은 다음과 같다.

\(\begin{align}\varepsilon_{c}=\left(\frac{\pi B H}{2 E^{\prime}}\right)^{2} R\end{align}\) , where B = 0.454 + 0.4ν       (11)

2.2. 커빅 커플링의 Geometry

Fig. 2는 커빅커플링 형상의 예를 나타낸다. 각 치는 α의 경사각을 가지며, 오목 치와 볼록 치가 접촉한다. 가공툴이 그리는 경로의 반경 r_tm을 중심으로 가공 폭 b(= r_o − r_I)으로부터 커플링 형상이 결정된다. 이런 형상을 고려하여 커빅커플링의 접촉 면적 및 관성 모멘트를 구해야 구할 수 있다.

Fig. 2. Geometry of curvic coupling.

2.2.1. 커빅 커플링의 치 형상과 접촉면적

Fig. 3은 커빅커플링 한 개의 기어치에 대한 정사영 단면을 나타낸다. 커빅커플링의 전체 접촉면적 A_n은 식(12)과 같이 나타낼 수 있다.

Fig. 3(a), (b). Tooth geometry of curvic coupling.

A_n = 2na_p / sin α       (12)

한 개의 기어치 접촉면을 x − y 평면에 정사영한 면적 a_p는 식(13-15)와 같다. 여기서 k는 커플링 치의 번호이다.

\(\begin{align}\begin{array}{l}a_{p}=\int_{a}^{\Box} d a=\int_{a}^{\Box} r \; d r d \theta \\=\int_{r_{i}}^{r_{o}} \int_{\theta_{1}^{k}}^{\theta_{2}^{k}} r \; d r d \theta+\int_{r_{i}}^{r_{o}} \int_{\theta_{3}^{k}}^{\theta_{4}^{k}} r \; d r d \theta \\ =2 \int_{r_{i}}^{r_{o}} \int_{\theta_{1}^{k}}^{\theta_{2}^{k}} r \; d r d \theta \end{array}\end{align}\)       (13)

\(\begin{align}\begin{array}{l}\theta_{m}^{k}=\frac{\pi}{n}\{1+2(k-1)\} \\ \theta_{m 1}^{k}=\theta_{m}^{k}-\frac{\pi}{2 n_{g}} \\ \theta_{m 2}^{k}=\theta_{m}^{k}+\frac{\pi}{2 n_{g}}\end{array}\end{align}\)       (14)

θ^k₁ = θ^k_m1 - θ

θ^k₂ = θ^k_m1 + θ

θ^k₃ = θ^k_m2 - θ

θ^k₄ = θ^k_m2 + θ       (15)

이 때, θ는 커플링의 기하학적 관계에 따라 아래와 같은 값을 가진다.

\(\begin{align}\theta=\frac{l \sin (\alpha)}{2 R_{c}}\end{align}\)       (16)

2.2.2. 커플링 형상을 고려한 접촉 및 굽힘 강성

커빅커플링 접촉면과 회전체 로드(rod)에 의한 체결력 F_z는 경사각을 가지므로, 접촉면에 대한 수직 및 접선방향 하중을 고려해야 한다. z방향 변위가 접촉면에 대해 가지는 수직 및 접선 방향 변위의 관계를 Fig. 4에 나타내었다.

Fig. 4. The direction of displacements.

단위 면적 당 접촉 강성 k_n으로부터 단위 면적 당 하중에 의한 변위를 계산할 수 있다. 하중의 방향과 접촉면이 경사를 이루고 있으므로 GW모델을 이용하여 접촉력을 계산할 때 접선방향 하중 역시 고려한다.

Δ_n = Δ_z sin(α)

Δ_t = Δ_z cos(α)       (17)

이때, GW 모델의 단위 면적당 접촉 압력을 σ_n이라 하면 식(18)이 성립한다.

\(\begin{align}\Delta_{n}=\frac{\sigma_{n}}{k_{n}}\end{align}\)       (18)

접선 방향 하중은 미끄러지는 영역과 미끄러지지 않는 영역에 따라 별도로 고려해야 한다. Mindlin[14]은 미끄러지지 않는 조건에서 헤르츠 구 접촉 모델에서 접촉 강성과 접선 방향 강성의 관계가 다음과 같음을 제시하였다.

\(\begin{align}\frac{k_{t}}{k_{n}}=\frac{2(1-v)}{2-v}\end{align}\)       (19)

접선방향 응력은 미끄러지지 않는 조건(no-slip)과 미끄러지는 조건(slip)으로 나눠지므로 아래와 같이 표현 가능하다.

σ_t = k_tΔ_t (σ_t ≤ μσ_n) : no - slip

σ_t = μσ_{n =} μk_nΔ_n (σ_t > μσ_n) : slip       (20)

따라서 z방향 응력은 식(21)와 같으며, 하중 F_s와 전체 수직 접촉강성은 식(22,23)와 같다.

σ_z = σ_nsin(α) + σ_tcos(α)       (21)

F_z = σ_zA_n

K_z = A_nk_z       (22)

\(\begin{align}=\left[\begin{array}{c}A_{n}\left(k_{n} \sin ^{2} \alpha+k_{t} \cos ^{2} \alpha\right) \\ A_{n}\left(k_{n} \sin \alpha(\sin \alpha+\mu \cos \alpha)\right)\end{array}\right.\end{align}\)       (23)

단위면적 당 접촉 강성 k_s에 관성 모멘트를 곱하면 미소 면적에서의 굽힘 강성 k_b를 구할 수 있고, 단위 면적에 대한 굽힘 강성을 접촉 면적 전체에 대해 적분하여 계산할 수 있다.

k_b = k_zy²dA_i       (24)

K_b = ∫k_bdA       (25)

2.3. 복 축 동역학 해석

2.3.1. Beam의 가정

Fig. 5는 본 연구에 사용된 회전체동역학 모델의 개략도를 나타낸다. 1차원 유한요소 티모섕코 빔(Timoshenko beam)모델을 사용하였다. 회전체는 4개의 디스크, 축, 베어링 및 5개의 커빅커플링으로 구성되며, 각 노드는 식(26)과 같이 4개의 자유도를 가진다.

Fig. 5. Rotor geometry.

q = {x θ_x y θ_y}^T       (26)

회전체동역학 해석을 위해 다음과 같이 가정하였다.

a) 디스크는 강체로 모델링 됨

b) 지지 베어링은 선형 강성과 감쇠를 가짐

c) 커빅커플링의 x방향 굽힘 강성은 굽힘 모멘트 M_x에만 영향을 받고, y방향 굽힘 강성은 M_y에만 영향을 받음

d) 복 축 회전체는 서로 선형 베어링으로 연결됨

2.3.2. 각 요소별 운동방정식

디스크의 운동방정식은 식(27, 28)과 같다.

\(\begin{align}\left[M^{d}\right]\left\{\ddot{q}^{d}\right\}-\Omega\left[G^{d}\right]\left\{\dot{q}^{d}\right\}=\left[F^{d}\right]\end{align}\)       (27)

\(\begin{align}\left[M^{d}\right]=\left[\begin{array}{cccc}M_{d} & 0 & 0 & 0 \\ 0 & I_{d} & 0 & 0 \\ 0 & 0 & M_{d} & 0 \\ 0 & 0 & 0 & I_{d}\end{array}\right],\left[G^{d}\right]=\left[\begin{array}{cccc}0 & 0 & 0 & 0 \\ 0 & 0 & 0 & -I_{p} \\ 0 & 0 & 0 & 0 \\ 0 & I_{p} & 0 & 0\end{array}\right]\end{align}\)       (28)

회전체 지지 베어링의 지지력 특성은 식 (29)과 같으며, 베어링 감쇠행렬과 강성행렬은 식 (30)과 같다.

\(\begin{align}-\left[C^{B}\right]\left\{\dot{q}^{B}\right\}-\left[K^{B}\right]\left\{q^{B}\right\}=\left\{F^{B}\right\}\end{align}\)       (29)

\(\begin{align}\left[C^{B}\right]=\left[\begin{array}{cccc}c_{x x} & 0 & c_{x y} & 0 \\ 0 & 0 & 0 & 0 \\ c_{y x} & 0 & c_{y y} & 0 \\ 0 & 0 & 0 & 0\end{array}\right],\left[K^{B}\right]=\left[\begin{array}{cccc}k_{x x} & 0 & k_{x y} & 0 \\ 0 & 0 & 0 & 0 \\ k_{y x} & 0 & k_{y y} & 0 \\ 0 & 0 & 0 & 0\end{array}\right]\end{align}\)       (30)

복 축 회전체를 연결하는 베어링의 특성은 식(31)과 같고, 감쇠 및 강성행렬은 식(32, 33)와 같다.

\(\begin{align}-\left[C^{I B}\right]\left\{\dot{q}^{I B}\right\}-\left[K^{I B}\right]\left\{q^{I B}\right\}=\left\{F^{I B}\right\}\end{align}\)       (31)

\(\begin{align}\left[C^{I B}\right]=\left[\begin{array}{cccc}c_{x x} & c_{x y} & -c_{x x} & -c_{x y} \\ c_{x y} & c_{y y} & -c_{y y} & -c_{y y} \\ -c_{x x} & -c_{x y} & c_{x x} & c_{x y} \\ -c_{y x} & -c_{y y} & c_{y x} & c_{y y}\end{array}\right]\end{align}\)       (32)

\(\begin{align}\left[K^{I B}\right]=\left[\begin{array}{cccc}k_{x x} & k_{x y} & -k_{x x} & -k_{x y} \\ k_{x y} & k_{y y} & -k_{y y} & -k_{y y} \\ -k_{x x} & -k_{x y} & k_{x x} & k_{x y} \\ -k_{y x} & -k_{y y} & k_{y x} & k_{y y}\end{array}\right]\end{align}\)       (33)

커빅커플링 굽힘 강성을 회전체 해석에 적용하기 위해 다음과 같이 가정한다.

a) 커플링 결합부를 길이가 0이고 질량이 없는 요소(element)로 가정함

b) 해당 요소는 강성행렬로만 표현됨

Fig. 6. Curvic coupling stiffness connection.

해당 요소에서 연결된 두 디스크가 반경방향으로 상대 변위가 없다고 가정하기 위해 반경방향 연결강성을 충분히 크게 설정함

\(\begin{align}\begin{array}{l} \left[K^{\text {Coupling }}\right]\\=\left[\begin{array}{cccccccc} k_{r} & 0 & 0 & 0 & -k_{r} & 0 & 0 & 0 \\ 0 & k_{b} & 0 & 0 & 0 & -k_{b} & 0 & 0 \\ 0 & 0 & k_{r} & 0 & 0 & 0 & -k_{r} & 0 \\ 0 & 0 & 0 & k_{b} & 0 & 0 & 0 & -k_{b} \\ -k_{r} & 0 & 0 & 0 & k_{r} & 0 & 0 & 0 \\ 0 & -k_{b} & 0 & 0 & 0 & k_{b} & 0 & 0 \\ 0 & 0 & -k_{r} & 0 & 0 & 0 & k_{r} & 0 \\ 0 & 0 & 0 & -k_{b} & 0 & 0 & 0 & k_{b} \end{array}\right] \end{array}\end{align}\)       (34)

복 축 회전체-베어링 시스템의 운동 방정식을 식(35)와 같이 나타낼 수 있다.

\(\begin{align}\begin{array}{l}{\left[\begin{array}{cc}M_{I} & 0 \\ 0 & M_{O}\end{array}\right]\left\{\begin{array}{l}\ddot{q}_{I} \\ \ddot{q}_{O}\end{array}\right\}+\left[\begin{array}{cc}C^{B}+C^{I B}+\Omega G_{I} & -C^{I B} \\ -C^{I B} & C^{I B}+\Omega G_{O}\end{array}\right]\left\{\begin{array}{l}\dot{q}_{I} \\ \dot{q}_{O}\end{array}\right\}} \\ +\left[\begin{array}{cc}K_{I}^{S}+K^{I B} & -K^{I B} \\ -K^{I B} & K_{O}^{S}+K^{I B}+K^{\text {Coupling }}\end{array}\right]\left\{\begin{array}{l}q_{I} \\ q_{O}\end{array}\right\}=\left\{\begin{array}{l}F_{I} \\ F_{O}\end{array}\right\} \\\end{array}\end{align}\)       (35)

3. 해석 결과 및 고찰

3.1. 해석 모델 및 운전 조건

조도측정기를 사용해 측정된 커빅커플링 치 면의 거칠기 프로파일은 Fig. 7과 같다. Table 2와 같이 왜도(Skewness)는 0에, 첨도(Kurtosis)는 3에 충분히 가까운 값을 가지기 때문에 Gaussian분포를 따른다 볼 수 있으며, GW 모델을 사용하기에 적합하다고 판단하였다.

Fig. 7. Surface profile of curvic coupling.

Table 2. Statistical surface properties

Table 3. Surface parameter

식(1-5)를 이용하여 해당 표면의 매개변수는 다음과 같이 계산되었다.

재질, 운전조건, 회전체 형상 및 커빅커플링 형상 정보는 Table 4에 제시되었다.

Table 4. Geometry of rotor system

3.2. 접촉 강성과 굽힘 강성

Fig. 8은 체결력에 따른 커빅커플링 접촉강성 및 회전체 굽힘강성을 나타낸다. 등가 굽힘 강성은 접촉 강성과 동일한 경향을 나타내며, 체결력에 따라 증가하는 경향을 보인다. 이는 굽힘 강성이 접촉 강성과 선형적인 관계를 가지기 때문이다.

Fig. 8. Contact and bending stiffness according to tightening force.

3.3. 커빅커플링이 포함된 회전체 동역학 해석

3.3.1. 체결력에 따른 위험속도 변화

2.3.2에서 언급된 커빅커플링의 등가 굽힘 강성을 회전체동역학 모델에 적용하여 RFR의 위험속도를 계산하였다. Fig. 9와 Fig. 10은 체결력의 크기에 따른 1-4차 위험속도를 나타낸다. 커빅 커플링의 영향을 평가하기 위해 해석에 사용된 회전체와 동일한 형상을 가지고 하나의 구조(uni-body)로 가정된 회전체의 위험속도를 계산하여 비교하였다.

Fig. 9. 1^st and 2^nd critical speed according to tightening force.

Fig. 10. 3^rd and 4^th critical speed according to tightening force.

체결력이 증가함에 따라 위험속도가 증가하며, 하나의 구조로 가정된 회전체의 위험속도로 수렴하는 경향을 보인다.

3.3.2. 위험속도에서의 모드형상

본 절에서는 다양한 위험 속도에서 회전체 시스템의 진동 모드를 분석하였다.

위험속도에서의 복 축 모델의 1-4차 모드 형상을 Fig. 11에 도시하였다. 안쪽 회전체는 파란색 점선, 바깥쪽 회전체는 붉은색 점선으로 표시되었다.

Fig. 11. Mode shapes at critical speed.

1차 위험속도에서는 원통형(cylindrical)모드가, 2차 위험속도에서는 원추형(conical)모드가 관찰되었다. 3차 위험속도에서는 2, 3차 굽힘 모드가 동시에 나타나는 경향을 보인다. 하지만 4차 위험속도에서는 여러 모드가 동시에 나타나면서 안쪽 회전체가 외부 회전체 보다 진폭이 큰 부분이 존재하는 모드 형상을 나타낸다.

4. 결론

본 연구에서는 커빅커플링의 표면을 측정하여 거친 표면 매개변수를 계산하여 이를 GW 모델에 적용시켜 거친 표면의 접촉 강성을 계산하였다. 간섭량에 따라 달라지는 접촉하중 모델을 사용하여 탄-소성 접촉해석을 수행하였고, 커빅커플링의 형상 및 미끄러짐을 고려하여 굽힘 강성을 계산하였다. 또한, 복 축 구조를 가지는 회전체의 티모섕코 빔 기반 유한요소 모델을 개발하였으며, 커빅커플링의 굽힘 강성을 회전체 강성 요소로 첨가하는 방법을 제안하였다. 시스템 모델링 후 고유 값 해석을 진행하여 위험속도를 예측하고 모드 형상을 예측하였다.

본 연구를 통해 아래의 결과를 도출하였다.

1. 커빅커플링의 접촉을 고려한 회전체의 위험속도는 체결력이 증가함에 따라 증가하는 양상을 보였다.

2. 커빅커플링의 접촉을 고려한 회전체의 위험속도는 단일 구조(uni-body)를 가지는 회전체 모델의 위험 속도 보다 항상 작았다.