Introduction

로터다이나믹 해석을 위하여 회전체를 모델링 하는 방법에는 전달행렬법, 유한요소법, 유한차분법 등이 있다.^[14] 본 연구에서는 1차원 빔 유한요소 모델을 사용하여 회전체를 모델링 하였으며, 각 절점의 자유도는 축 방향 자유도를 제외한 4자유도 (x,x_θ,y,y_θ)를 가진다. 아래 식 (1)을 통해 로터다이나믹 해석을 수행하였다.^[10]

여기서 [M] 행렬은 회전체의 질량 행렬, [C_G] 는 회전기계의 고유 특성인 자이로스코픽 행렬, [C_B, K_B] 행렬은 각각 저널 베어링의 감쇠, 강성 행렬, [K_S] 는 회전체의 강성행렬, F 는 불평형력, 외력과 같은 하중 벡터, q 는 변위 벡터를 의미한다. 요소 i 절점에 대한 행렬은 아래와 같다.

여기서 한 요소 i 절점에서의 변위벡터 {q_e,i}는 x 방향 변위, x 방향 각변위, y 방향 변위, y 방향 각변위 총 4 자유도(x,x_θ,y,y_θ)로 구성되어 있으며, 축 방향의 자유도는 고려되지 않았다. I_T,I_p 는 각각 한 요소에 대한 질량 관성 모멘트, 극 관성 모멘트를 의미하며, E, L은 각각 한 요소의 Young’s Modulus와 요소의 길이를 의미한다.

회전체를 모델링 하기 위하여 각 요소에 대한 행렬⁽²⁾을 만든 후 각 절점에 대하여 중첩법을 사용하였다. 베어링 동특성, 자이로 스코픽 효과, 추가질량 등을 포함한 회전체 통합 행렬로 로터다이나믹 해석을 진행하였다.

본 연구에 사용된 로터다이나믹 해석 모델은 Fig. 1과 같은 형태로 4개의 블레이드를 가지고 있으며, 두 개의 틸팅패드 저널베어링을 통해 회전체의 하중을 지지하며, 총 36개의 요소로 구성되어 있다.

Fig. 1 
Rotor-bearing system with 4 blades

4개의 블레이드 무게를 포함한 총 회전체의 질량은 약 1900kg, 길이는 3.36m이며 무게중심은 회전체 좌측 끝단에서 1.776m 떨어진 곳에 위치한다.

본 연구에 사용된 해석 모델에서 Thrust Collar는 존재하지만, 유한요소 모델링 시 축 방향 자유도를 고려하지 않았기 때문에 쓰러스트 베어링의 동특성과 영향은 무시하였다.

본 연구에서 사용된 베어링은 틸팅패드 저널베어링이며, 동적 안전성이 우수하여 고성능, 고출력, 고하중 회전기계에 자주 사용된다. 틸팅패드 저널베어링의 강성과 감쇠를 구하기 위하여 Reynolds 방정식(3)과 유한요소법을 통해 베어링 패드에 걸리는 유막의 압력 분포를 계산하였다. 윤활유 온도 변화에 따른 점도 변화를 무시하였으며, 윤활유는 비압축성으로 가정하였다.

여기서 h는 유막의 두께를 의미하며, 샤프트와 패드의 열변형, 탄성 변형을 무시하고 회전체와 패드의 움직임만 고려하게 된다면 유막 두께 h는 식(4)과 같이 표현할 수 있다.

여기서 C_P, C_B 는 패드와 베어링 간극을 의미하며, X, Y 는 각각 회전체의 중심 위치, δ, R 은 패드의 틸팅각, 패드의 틸팅 회전 반경을 의미한다.

식(3), 식(4)을 통해 유막의 압력을 계산하여 회전체의 평형위치를 찾은 후 Fig. 2와 같이 회전체와 패드의 미소 변위(Δx, Δy, Δδ), 미소 속도(ΔV_x, ΔV_y, ΔΩ_δ)에 대하여 섭동법을 진행하여 강성 및 감쇠 행렬(5)을 구하였다.

Fig. 2 
Prediction of TPJB dynamic characteristics using perturbation method

여기서 k, c 행렬은 평형위치에 있는 저널에 대하여 미소 변위(±Δx,± Δy,± Δδ), 미소 속도(±V_x,± V_y,± Ω_δ)로 섭동하여 힘의 변화량을 통해 아래 식(6)와 같이 구하였다.

식(6)의 하첨자 a,b는 각각 x, y, δ 방향을 의미하며, k_ab, c_ab의 의미는 b방향 미소 변위, 미소 속도에 의한 a방향 하중변화를 나타낸다.

행렬(5)를 통해 동기식 축약 동특성(Synchronously reduced dynamic coefficient)을 계산하였으며, 식(1)에 동기식 축약 동특성 [K_B, C_B]를 대입하여 로터다이나믹 해석을 수행하였다.

Fig. 3은 경사 환경에 놓인 로터-베어링 시스템을 나타내고 있다. 경사 환경(피칭, 요우잉)에 대한 운동방정식을 세우기 위하여 로터-베어링 시스템은 강체로 가정하였다.

Fig. 3 
Rotor-bearing system placed in the sloped environment

또한 전체 시스템의 회전중심에서 로터-베어링 시스템의 무게중심까지의 위치 벡터 r→ 은 고정되어있다고 가정하였다.

CR, CG는 각각 시스템의 회전중심, 회전체의 무게중심을 의미한다. 또한 여기서 θ는 피칭, 요우잉에 대한 회전 각도를 의미하며, eθ→ 는 위치 벡터 r→의 법선 벡터를 의미한다.

동적 경사 환경에 의한 강체 운동방정식은 식(7)과 같으며, 여기서 M 은 로터-베어링 시스템의 총 질량, ω 는 피칭, 요우잉의 회전 속도를 의미한다.

최대 경사각, 회전 속도를 통해 구한 동적 경사 하중(F_sloped)을 과도응답 해석 시 베어링에 작용하는 외력으로 대입하였다. 동적 경사 하중을 적용한 로터다이나믹 해석 운동 방정식은 식(8)과 같다.

Fig. 5는 본 연구에 사용된 로터-베어링 시스템에 대한 위험속도 선도로서, 가로축은 베어링 강성을 의미하며, 세로 축은 운전속도를 의미한다. 위험속도 선도는 회전체의 정격 운전속도에서 가로축으로 선을 긋고, 가로선과 그래프 상의 선이 만나는 지점이 위험속도가 된다.

Fig. 5 
Critical speed map

AP 684^[16]에서 권장하는 분리 여유 범위 내 그래프 상의 선이 만나는 지점이 없도록 베어링 강성 범위를 설정하였다.

본 연구에서는 API 648^[16]에 따라 16 이상의 분리 여유를 갖는 베어링 강성 범위를 선정하였으며, 그 결과 베어링 강성의 범위를 6×10⁷∼3×10⁸ N/m로 설계하였다.

Table 1은 설계된 베어링 결과를 나타내며 베어링 타입은 5 pad, LOP type이며, 윤활유는 ISO VG 32로 선정하였다.

육상 고정 상태에서의 좌측, 우측 베어링의 지지하중은 각각 8,300N, 10,500N 이다. 좌측 베어링의 x 방향 강성은 약 1.15×10⁸ N/m, y 방향 강성은 약 8.08×10⁷ N/m, x 방향 감쇠는 약 1.34×10⁵ N∙s/m이며, y 방향 감쇠는 약 1.11×10⁵ N∙s/m 이다. 우측 베어링의 x 방향 강성은 약 1.36×10⁸ N/m 이며, y 방향 강성은 약 8.41×10⁷ N/m 이며, x 방향 감쇠는 약 1.48×10⁵ N∙s/m 이며, y 방향 감쇠는 약 1.13×10⁵ N∙s/m 이다.

최대 경사각(피칭, θ_max=45°)에서 지지 하중이 감소함에 따라 베어링의 동특성이 감소하였다. 최대 경사 상태에서의 좌측 베어링 x 방향 강성은 약 8.26×10⁷ N/m, y 방향 강성은 약 7.44×10⁷ N/m, x 방향 감쇠는 약 1.12×10⁵ N∙s/m, y 방향 감쇠는 약 1.06×10⁵ N∙s/m 이다. 우측 베어링의 x 방향 강성은 약 8.86×10⁷ N/m, y 방향 강성은 약 7.58×10⁷ N/m, x 방향 감쇠는 약 1.16×10⁵ N∙s/m, y 방향 감쇠는 약 1.07×10⁵ N∙s/m 이다.

베어링 상세 설계 후, 고유값 해석을 통하여 캠벨 선도를 구하였다. 이때 육상에 고정되어 있는 상태에서의 베어링 동특성을 통해 예측된 캠벨 선도를 Fig. 6 (a)에 나타내었다. 또한 최대 경사 상태(피칭, θ_max=45°)에서의 저널 베어링 동특성을 통해 예측된 캠벨 선도를 Fig. 6 (b)에 나타내었다.

Fig. 6 
Campbell diagram (a) in the land-based state and (b) in the sloped state

해석 결과, 고정상태에서의 위험속도는 1,863 rpm, 2,012 rpm, 5,326 rpm, 6,765 rpm, 10,451 rpm이며, 최대 경사 상태에서의 위험속도는 1,825 rpm, 1,878 rpm, 4,981 rpm, 5,905 rpm, 9,964 rpm이다.

고정상태와 최대 경사 상태 모두 첫 번째 위험속도는 강체 모드 중 원통형(cylindrical) 모드이며, 두 번째 위험속도는 강체 모드 중 원추형(conical) 모드이며, 세 번째 위험속도의 모드 형상은 첫 번째 굽힘모드로 나타나였다. 정격 운전속도 9,000 rpm에 대한 육상에 고정상태에서의 분리 여유는 16.13 %이며, 최대 경사 상태에서의 분리 여유는 10.71 %이며, 경사 환경에 의하여 분리 여유가 5.42 %p 감소하였다.

불평형 질량과 경사 상태에 의한 회전체의 진동 응답 특성을 확인하기 위하여 조화응답 해석을 수행하였다. 조화응답 해석 시 회전체의 무게중심 위치에 불평형 질량을 부과하였다. 또한 불평형 질량 위치에 따른 진동응답을 비교하기 위하여 9,000rpm에서 가장 가까운 위험속도 모드형상의 최대변위 지점에 불평형 질량을 부과하였다. 이때 불평형 등급은 ISO G 2.5를 적용하여 조화응답 해석을 수행하였다. 정적 경사 상태에 대한 조화응답을 구하기 위하여 피칭 최대 경사각(θ=45°)에서의 베어링 동특성을 사용하였다.

Fig. 7은 무게중심에 불평형 질량을 부과하였을 때의 조화응답 해석 결과로서, (a)와 (b)는 각각 좌측 베어링에서의 진동 진폭, 위상 응답이며 (c)와(d)는 각각 우측 베어링에서의 진동 진폭, 위상 응답이다. 전 운전영역에 대하여 정적 경사 상태에서의 진동 응답이 육상에 고정되어 있는 상태에 비하여 큰 것을 확인하였다.

Fig. 7 
Unbalance response of (a) amplitude and (b) phase at the left journal bearing and (c) amplitude and (d) phase at the right journal bearing (Unbalance: 5101 gㆍmm, Location: gravity center)

정격 운전속도(9,000 rpm)에서의 해석 결과, 육상에 고정되어 있는 상태에서의 좌, 우측 베어링에서 측정되는 진동 실효값(rms)은 각각 4.325 μm, 4.354 μm이며, 최대 정적 경사 상태에서의 좌, 우측 베어링에서 측정되는 진동 실효값은 각각 5.207 μm, 4.889 μm이다.

육상에 고정되어 있는 상태의 불평형 응답이 최대 정적 경사 상태에서의 진동 응답보다 약 10∼20 % 작게 예측되었으며 육상 고정상태, 최대 정적 경사 상태 모두 ISO 10816-3^[15]에서 권장하는 진동 실효값 제한치인 37 μm 이하로 나오는 것을 확인하였다.

Fig. 8은 9,000rpm에서 가장 가까운 위험속도 모드형상의 최대변위 지점에 불평형 질량을 부과하였을 때의 조화응답 해석 결과로서, (a)와 (b)는 각각 좌측 베어링에서의 진동 진폭, 위상 응답이며 (c)와(d)는 각각 우측 베어링에서의 진동 진폭, 위상 응답이다. 불평형 질량을 무게중심에 부과하였을 때와 같이 전 운전영역에 대하여 정적 경사 상태에서의 진동 응답이 육상에 고정되어 있는 상태에 비하여 큰 것을 확인하였다.

Fig. 8 
Unbalance response of (a) amplitude and (b) phase at the left journal bearing and (c) amplitude and (d) phase at the right journal bearing (Unbalance: 5101 gㆍmm, Location: maximum displacement)

9,000rpm에서 가장 가까운 위험속도의 최대변위 지점에 불평형 질량 부과 시 정격 운전속도(9,000 rpm)에서의 해석 결과, 육상에 고정되어 있는 상태에서의 좌, 우측 베어링에서 측정되는 진동 실효값(rms)은 각각 28.58 μm, 23.62 μm이며, 최대 정적 경사 상태에서의 좌, 우측 베어링에서 측정되는 진동 실효값은 각각 33.20 μm, 27.42 μm이다.

육상에 고정되어 있는 상태의 불평형 응답이 최대 정적 경사 상태에서의 진동 응답보다 약 13∼16 % 작게 예측되었으며 육상 고정상태, 최대 정적 경사 상태 모두 ISO 10816-3^[15]에서 권장하는 진동 실효값 제한치인 37 μm 이하로 나오는 것을 확인하였다.

주기를 갖는 동적 경사 환경에 의한 회전체의 진동 응답을 확인하기 위하여 과도응답 해석을 수행하였다. 과도응답 해석을 수행하기 위하여 식(7)를 통해 동적 경사 환경에 의한 동적 경사 하중을 계산하였으며, 식(8)에 시간의 함수로 경사 하중을 고려하였다.

피칭 각속도에 의한 로터-베어링 시스템이 받는 진동 응답을 평가하기 위하여 동적 경사 하중을 베어링에 y축 방향으로 작용하는 하중으로 고려하여 해석을 수행하였다. 불평형 질량에 의한 불평형력는 무게중심 위치에 부과하였으며 크기는 ISO G 2.5를 적용하였다. Fig. 9는 동적 경사 하중과 불평형력을 비교한 그림으로 동적 경사 하중은 약 2,500 N의 크기, 10 초의 주기(=0.1 Hz)를 가지며, unbalance force는 약 2,000 N, 150 Hz로 적용하였다.

Fig. 9 
Dynamically sloping force and unbalance force

Fig. 10은 동적 경사 하중과 불평형력에 의한 과도응답 해석 결과로서, y축 방향으로 작용하는 동적 경사 하중에 의하여 x축 진동 응답보다 y축 진동 응답이 큰 것을 확인하였다. 해석 결과 동적 경사 환경에의 진동 실효값(rms)은 좌측과 우측 각각 35.89 μm, 34.02 μm이다.

Fig. 10 
Transient analysis response by dynamically slopping force with unbalance force: (a) 0∼10 sec. and (b) 0∼2 sec.

Table 2는 각 해석 조건에 따른 진동 실효값(rms)의 결과로서 (a)는 불평형력을 무게 중심에 부과하였을 때, (b)는 9,000rpm에서 가장 가까운 위험속도 모드형상의 최대변위 지점에 불평형 질량을 부과하였을 때의 결과를 나타낸다. 육상에 고정되어 있을 경우 분리 여유는 16.13 %, 피치각이 30°일 경우 분리 여유는 14.72 %, 피치각이 45°일 경우 10.71 %의 분리 여유를 가진다. 불평형력의 위치와 상관없이 고정상태 대비 정적 경사각(피치각)이 증가함에 따라 베어링의 지지하중이 감소하게 되고, 이에 따라 베어링 강성이 감소하여 분리 여유가 감소하는 것을 확인하였다.

불평형력이 무게 중심에 위치하고 피치각 45° 그리고 동적 경사 주기가 10초일 경우 베어링에서의 최대 진동 실효값은 35.89 μm이며, 동일한 조건에서 주기가 20 초로 바뀔 경우 베어링에서의 최대 진동 실효값은 14.79 μm로 감소하게 된다.

또한 불평형력이 최대변위 위치하고 피치각 45° 그리고 동적 경사 주기가 10초일 경우 베어링에서의 최대 진동 실효값은 48.51 μm이며, 동일한 조건에서 주기가 20 초로 바뀔 경우 베어링에서의 최대 진동 실효값은 36.86 μm로 감소하게 된다. 최대 경사각, 주기에 따라 ISO 10816-3^[15]에서 제한하는 진동 실효값인 37 μm보다 큰 값을 가질 수 있음을 확인하였다.