회전체의 하중을 지지하는 베어링의 성능 예측은 터보기계의 작동 신뢰성을 검증하는 단계에서 중요한 요구 사항 중 하나이다. 최근 터보기계의 고효율화와 경량화로 인해 베어링의 발열이 중요한 문제가 되고 있다. 더불어 개인 컴퓨터의 비약적인 성능향상과 수치해석 소프트웨어의 발전을 기반으로 회전체-베어링 시스템의 온도를 예측하기 위한 다양한 연구가 시도되고 있다.

Ha 와 Kim[1]은 대형 틸팅패드 베어링의 정적특성을 예측하고자, 와점성계수를 이용한 난류윤활이론을 사용하였다. Kim과 Palazzolo[2,3]는 2차원 유한요소모델을 이용하여 틸팅패드 저널베어링 수치해석 모델을 개발하였다. 특히 베어링 패드의 탄성변형을 고려하는 데에 따른 계산의 효율성 향상을 위하여 모드좌표 변환을 이용하였다. 본 연구는 베어링 패드의 탄성변형 예측을 위하여 유한요소법과 모드합성법을 적용한 최초의 연구였다. Gadanngi와 Palazzolo[4]는 시간과도응답 해석을 통하여 시간에 따른 베어링 시스템의 응답 변화를 해석하였다. Monmousseau 와Fillon[5]은 동적인 하중 하에서의 틸팅패드 저널베어링의 시간과도응답을 해석하였다. Nicholas[6]는 시스템 온도를 낮추기 위한 틸팅패드 저널베어링 설계에 대하여 연구하였다. Tschoepe 와 Childs[7]는 험 장비에서 고온 및 저온 상태에서의 유막 두께를 측정하였다. Suh와 Palazzolo[8,9]는3차원 유한요소모델을 이용하여 베어링 패드의 탄성변형과 열변형을 동시에 고려한 TEHD윤활모델을 개발하였으며, 시간에 따른 패드의 탄성변형을 고려할 수 있는 고정밀 수치해석 모델을 개발하였다. 특히 Kim과 Palazzolo[2,3]에 의해 최초 제안 되었던 유한요소법과 모드합성법을 이용하여 3차원 베어링 패드를 모델링 하였으며, 베어링 동특성계수 행렬 유도에 모드좌표와 물리좌표를 결합한 방법을 최초로 제안하였다. El-Butch와 Ashour[10]는 TEHD 윤활 모델을 사용하여 비정상 부하 조건에서 각정렬 불량 상태의틸팅패드저 저널베어링의 성능을 연구하였다. Suh와 Choi[11]는 틸팅패드 저널베어링의 피봇 형상과 각정렬불량에 따른 정특성 및 동특성을 연구하였다. Suh와 Hwang[12]은 베어링 패드의 열변형을 예압량 변화의 관점에서 정량화 시키는 연구를 수행하였다. Lee와 Sun[13]은 터빈 시뮬레이터용 틸팅패드 저널베어링에 대하여 열윤활 해석을 수행하고, 예측된 패드온도와 실험결과를 비교하였다. Lee[14]는 초임계 CO₂ 파워 터빈 회전축을 지지하는 틸팅패드 저널베어링의 열윤활 수치해석 결과를 실험적으로 측정한 패드 온도와 비교하였다. Yang과 Palazzolo[15-17]는 상용 소프트웨어인 CFX를 이용하여 틸팅패드 저널베어링의 성능변화를 해석하였으며, 패드의 열전달, 열변형 및 피봇의 탄성변형을 동시 고려하였다. 특히 그루브 에서의 오일혼합을 정밀하게 예측하였다.

위에서 언급한 연구 외에도 현재까지 수많은 저널베어링 수치해석 모델이 개발되었지만 베어링-유막-회전체 시스템의 온도 경계 조건 변화에 따른 성능 변화와 관련된 깊이 있는 연구는 많지 않다. 이는 베어링 패드 및 샤프트 주변의 열 전달을 모델링하고 그에 따른 열 변형을 예측하기 위한 수치해석 모델링의 어려움과 복잡성에 기인한다. 샤프트 축 길이 방향의 온도 경계 조건은 2차원 베어링 모델에서는 고려할 수 없다. 이는 2차원 모델에서 샤프트가 유막에 의해 둘러 쌓여 있기 때문이다. 다음은 현재까지 다양한 연구자들이 개발한 수치적 접근 방식이 정확한 온도 경계 조건을 고려하는 데 어려움을 겪는 이유를 요약한 것이다.

(a) 회전하는 저널과 유막 사이의 열 전달 모델링 (boundary A, Fig. 1a) (b) 유막과 베어링 패드 사이의 열 전달 모델링 (boundary B, Fig. 1a) (c) 샤프트 축 방향 열 경계 조건 모델링 (boundary C, Fig.1b) (d) 저널 및 베어링 패드의 열 변형으로 인한 유막두께 변화 해석 (Fig. 1b) (e) 베어링 동역학, 윤활, 열전달, 열변형 문제를 동시 고려한 다중물리 기반 모델링

유막과 베어링 패드의 온도는 베어링 수명과 밀접한 관련이 있으며, 구조물의 온도구배는 열 변형을 일으킨다. 열변형은 운전 상태에서의 베어링 유막 두께를 변화시키고, 변화된 두께는 유막의 점성 전단을 변화시켜 베어링 시스템의 온도 분포를 다시 변화시킨다. 이러한 복잡한 상호 작용은 베어링 시스템의 정확한 온도 예측을 어렵게 만든다.

본 연구에서는 3차원 유한요소모델에 기반한 베어링 수치 모델링 방법 및 알고리즘을 제시하고 온도 경계 조건의 변화에 따른 베어링의 성능변화를 분석하였다. 본 연구에서 사용하는 수치해석 모델은 Suh와 Palazzolo[8,9]가 2014년 수행한 연구를 기반으로 하였다. 본 연구에서 사용된 수치해석 모델은 상용수치해석 소프트웨어인 MATLAB 2019b 버전을 사용하여 인하우스 코드로 개발되었다. 본 수치모델은 Suh 와 Palazzolo [8,9]의 연구에서 타 연구결과[18,19] 와의 비교를 통하여 이미 검증되었다.

본 연구의 목적은 샤프트와 베어링 패드 주변의 온도 조건 변화에 따른 베어링의 정특성 및 동특성의 변화를 예측하는 것이다. 본 연구에서는 대류경계조건 (convection boundary condition)과 등온 조건 (constant temperature condition) 두가지를 베어링 패드와 샤프트에 적용하였다.

이러한 온도경계 조건의 변화는 베어링 성능변화에 큰 영향을 주는 것을 확인하였으며 대류경계조건에 비하여 등온조건이 더 큰 영향을 미치는 것을 확인하였다.

본 연구에서는 유막두께 방향의 점성변화를 고려할 수 있는 가변점성 레이놀즈방정식 (varying viscosity Reynolds equation)을 사용하여 유막에서 발생하는 압력을 계산하였다. 방정식 (1)은 가변점성레이놀즈 방정식을 나타내며, 유한요소법을 이용하여 수치모델링 하였다.

여기서 D₁과 D₂는 다음과 같이 정의된다.

유한요소모델에서의 노드 별 윤활유의 국부적인 점도는 온도만의 함수로 가정하였으며 방정식 (4)는 점성과 온도와의 관계를 나타낸다. 여기서 μ₀는 기준온도 T₀에서의 기준점도를 나타내며 β는 점도상수이다. 방정식 (5)에서 U는 윤활유 속도벡터를 나타낸다.

이러한 속도벡터는 압력의 구배(∇p)와 저널의 선속도(U)를 통해 계산할 수 있다. 윤활유의 속도는 유막두께방향으로 다른 값을 가지며 방정식(6)은 이러한 속도 구배로 인해 유막 내에서 발생하는 점성전단열, 열대류 및 열전도를 계산하는 에너지방정식을 나타낸다. 본 연구 에서는 베어링 축방향의 열구배를 고려할 수 있는 3차원 에너지방정식을 유한요소법을 이용하여 모델링 하였으며, 대류항으로 인해 발생하는 수치적 에러를 제거하고자 업윈드방법 (up-winding scheme)을 사용하여 가중함수를 수정하였다.

샤프트와 베어링 패드의 국부적인 열구배로 인한 구조물의 열변형은 유막의 두께를 변화시킨다. 변화된 유막두께는 베어링에서 발생하는 유압을 변화시키고, 변화된 유압은 베어링의 동적 거동을 변화시킴과 동시에 이로 인한 유속의 변화는 에너지방정식의 경계조건 및 유막발생 열을 변화시킨다. 본 연구에서는 방정식 (7)의 열전도 모델을 3차원 유한요소로 이산화 하여 모델링 하였다. 수렴과정에서 발생할 수 있는 발산 현상을 해결하고자 온도의 시간미분항을 고려하였다. 방정식 (8)은 베어링 구조물 주변의 온도경계조건을 나타내며, 등온조건 (S_T), 열속경계조건 (S_q) 및 대류경계조건 (S_h)을 사용할 수 있다. 아무 조건이 주어지지 않을 경우 단열조건이 성립된다.

방정식 (9)는 열전달 방정식 (7)의 유한요소 행렬모델을 나타내며, 앞에서 언급한 바와 같이 베어링 수렴특성에 따라 온도의 시간미분항을 고려하기도 한다. 방정식 (10)은 베어링 패드와 샤프트 구조물의 열변형 모델을 나타낸다.

본 연구에서는 베어링 시스템의 정적 수렴조건 해석을 위하여 일반적으로 사용하는 뉴튼랍슨 (Newton-Raphson) 방법을 사용하지 않고, 시간에 대하여 수치적분 하여 수렴조건을 계산하였다. 이는 본 연구에서 제안한 베어링 해석 알고리즘의 특징 중 하나이며 특히, 3차원베어링 모델에서 고려하는 패드의 틸트, 피치, 요우 운동, 그리고 피봇의 탄성변형 운동을 모두 고려할 경우 기존의 뉴튼랍슨 방법보다 계산시간 면에서 효율적인 수렴특성을 보인다. 본 연구에서 사용된 방법은 아니지만 패드의 탄성변형을 고려하는 경우 시간에 대한 수치적분 방법이 베어링 시스템의 수렴을 위한 유일한 방법으로 알려져 있다[8,9].

Fig. 2는 저널과 베어링 패드의 자유물체도를 나타내며, 식 (11)과 (12)는 각각 베어링패드와 샤프트의 운동방정식을 나타낸다. 패드는 레이놀즈 방정식에서 계산된 유압(F_Fluid)과 피봇에서의 지지력을 외력(F_Pivot)으로 받는다. 샤프트는 자중 및 기타 외력에 의한 하중(F_Load)과 유압(F_Fluid)을 외력으로 받는다.

Fig. 3은 본 연구에서 사용된 수치해석 알고리즘을 나타낸다. 총 3개의 폐쇄 루프(closed loop)가 존재하며, 이러한 루프는 수치해석 코드의 수치적분 알고리즘 내에서 매 단계마다 영향을 주고받으며 베어링의 정적 및 동적 특성을 수렴시킨다.

본 연구에서 성능예측에 사용된 베어링모델의 특성이 아래의 표와 그림에 나타나 있다. 윤활유 내에서의 점성 전단 열과 구조물 내에서의 열전달 특성을 고려하기 위하여 Table 1과 같이 온도에 따른 점도변화를 고려하였다. 또한 샤프트와 베어링 패드 구조물 내에서의 열전달 및 이로 인한 열변형 특성을 고려하기 위하여 Table 2와 같이 재료의 여러가지 특성을 고려하였고 수치모델에서 사용하였다. 소재의 국부적인 열변형은 온도구배에 의하여 발생하며, 온도변화에 따른 소재 물성치(밀도, 탄성계수, 열팽창계수 등) 변화는 없는 것으로 가정하였다.

Table 3은 베어링 모델의 기하학적 특성 및 운전조건을 나타낸다. Fig. 4(a)에서 보는 바와 같이 5개의 패드를 가지며 16,000N의 하중이 -y축 방향으로 패드 사이에 가해진다. 회전속도는 12,000 rpm으로 일정하며 오프셋은 0.5, 각 패드의 원주방향 길이는 60°를 가진다. Fig. 4(b)는 베어링의 3차원 유한요소 모델의 단면도를 나타낸다.

피봇의 강성을 고려하기 위하여 피봇을 경화스프링(hardening spring)으로 가정하였고 피봇에서 발생하는 힘(F_Pivot)을 식(13)의 형태로 모델링 하였다. x_p는 피봇에서 반경방향으로의 탄성변형량을 나타내며, k₁은 선형스프링 강성계수, k₃는 비선형 스프링강성을 나타낸다. Table 3에서 k₁과 k₃는 각각 선형, 비선형 강성값을 나타낸다. 본 값들은 대상 베어링의 피봇특성과 Hertzian 접촉이론 결과를 회귀분석하여 계산하였다.

본 연구에서는 Table 4에서 보는 바와 같이 총 4개의 경우에 대하여 매개변수 해석을 수행하였다. 먼저 1번 모델에서는 샤프트 주변의 대류온도조건(T_Shaft at Fig. 5(a))을 열대류계수 200 W/(m²K)에 주변온도를 50°C로 가정하였다. 윤활유 공급온도(T_Lube)는 50°C로 일정하게 유지시켰으며 패드 주변의 대류온도조건을 열대류계수 200 W/(m²K)에 주변온도를 10°C에서 90°C로 20°C간격으로 증가시키며 베어링 특성 변화를 해석하였다. 2번모델의 경우 샤프트 주변의 온도조건과 윤활유 공급온도는 1번과 같으며 패드 주변의 등온조건(constant temperature)을 10°C에서 90°C로 20°C간격으로 증가시키며 베어링 특성 변화를 해석하였다.

3번모델과 4번모델의 경우는 Fig. 5(a)의 T_Shaft와 같이 샤프트 주변의 온도를 열대류조건과 등온조건으로 각각 해석하였다. 이 때 패드 주변(T_Pad)은 열대류계수 200 W/(m²K)에 주변온도를 50°C로 가정하였다. 뒤에서 살펴볼 결과 그래프에서 사용될 선 모양(Line style)을 Table 4에 기입하여 독자의 이해를 돕고자 하였다.

먼저 베어링 정특성 해석을 수행하였다. Fig. 6는 저널과 패드의 정적 평형 상태를 나타낸다. Fig. 6(a)와 같이 저널 편심량의 경우 샤프트 주변의 등온조건(constant temperature condition)변화에 따른 변화가 가장 컸으며 패드 주변의 대류경계조건(convection temperature condition)에 따른 변화가 가장 작았다. 자세각의 경우는 이와 다르게 패드 주변의 등온조건 변화에 따른 특성이 가장 많이 바뀌었다. 패드의 정적 특성을 알아보고자 Fig. 6(c)와 같이 하중방향의 패드 중 하나인 5번 패드(see Fig. 6(b))의 틸팅각과 피봇탄성변형량의 변화를 살펴보았다. 패드와 샤프트의 등온조건 변화에 따른 특성변화는 비슷함을 볼 수 있다. 피봇의 탄성변형량은 패드주변의 온도 변화 보다 샤프트 주변의 온도변화에 민감함을 확인할 수 있다.

Fig. 7은 패드 내의 최대온도 변화를 나타낸다. 패드 주변의 등온조건 변화에 따른 패드의 최대온도 변화가 가장 뚜렷하였다. 이는 샤프트 주변의 온도변화가 패드의 온도에 미치는 영향이 작음을 의미하며, 패드 주변의 대류경계조건이 10°C에서 90°C까지 크게 변하더라도 그 영향이 등온조건에 비하여 작음을 알 수 있다.

Fig. 8는 무차원화 된 최소유막두께(h/c_b)를 나타낸다. 최소유막두께는 베어링의 운전건전성을 판단하기 위해 사용하는 중요 상태 중 하나이며 본 연구에서 사용된 수치해석 모델의 경우 샤프트의 반경방향 열팽창량 및 베어링 패드의 열변형에 큰 영향을 받는다. 패드 주변의 등온조건 변화 따른 최소유막두께의 변화가 가장 두드러졌으며 이는 패드 바깥 부분의 온도를 능동적으로 제어함으로써 베어링 건전성을 확보할 수 있음을 나타낸다. 이와 같이 패드 바깥부분의 등온조건 온도를 낮춤으로써 패드의 최대온도를 낮추고 최소유막두께를 증가시킬 수 있음을 알 수 있다.

Fig. 9는 베어링 축방향 중심부에서의 무차원화 된 저널의 반경방향 열팽창량을 나타낸다. 만약 베어링 패드의 열변형이 없다면 1.0의 열팽창은 저널과 베어링 패드의 표면이 닿는 것을 의미한다. 샤프트 주변의 온도변화가 샤프트의 열팽창에 가장 큰 영향을 미친다는 것을 볼 수 있으며 특히 샤프트 축방향 바깥쪽의 대류 경계조건의 변화에 비하여 등온조건 변화에 따른 샤프트의 열팽창량 변화가 크다는 것을 알 수 있다. 이는 패드와 마찬가지로 대류경계조건에 비하여 등온조건이 온도변화에 더 큰영향을 미치며 이로 인한 열변형특성 또한 큰 것임을 알 수 있다.

Fig. 10은 베어링 주변의 온도경계조건 변화에 따른 동력손실량 변화를 나타낸다. 주목할 만한 점은 다른 온도 조건변화에 비하여 패드 바깥면의 등온조건 변화에 따른 동력손실량 변화가 매우 크다는 것이다. 본 논문의 앞에서 살펴 보았던 Fig. 8과 Fig. 9에서와 같이 패드 바깥면 등온조건 온도가 낮을 경우 이로 인한 최소유막두께 증가에 따른 베어링 건전성 확보의 관점에서는 유리하지만, 동력손실량은 약 50%가량 증가함을 볼 수 있다. 이는 유막의 최대온도 감소로 인한 유막 점도의 증가를 그 원인으로 들 수 있다.

Fig. 11과 12는 각각 Table 4의 3,4번 모델에서 T_shaft: 90°C인 경우 베어링-샤프트 단면의 온도분포를 나타낸다. 등온조건의 경우 비하여 대류온도조건에서 패드의 원주방향 온도변화가 뚜렷하게 관찰되며, 최대온도는 5번 패드에서 발생한다. Fig. 13과 14는 각 베어링 패드에서의 유막두께 분포를 나타낸다. 주목할 만한 점은 유막두께의 경우 축방향으로의 변화가 발생한다는 것이다. 원주방향의 유막두께 변화량에 비해서는 작지만 베어링 축방향의 유막두께 변화량도 뚜렷하게 관찰된다.

Fig. 15는 온도경계조건 변화에 따른 베어링 동특성 변화를 나타낸다. 연성항 (cross coupled terms)의 경우 직결항 (direct terms)에 비하여 그 값이 매우 작기 때문에 본 논문에 나타내지 않았다. 먼저 Fig. 15(a), (b)는 각각 수평방향 및 수직방향 강성계수의 변화를 나타낸다. 두 경우 모두 변화하는 경향이 매우 비슷함을 볼 수 있다. 특히 샤프트 주변의 열경계조건의 변화가 베어링 패드에 비하여 강성계수의 변화에 큰 변화를 미친다는 것을 볼 수 있다. 모든 경우에 대하여 주변 온도가 증가할수록 강성계수가 증가함을 볼 수 있다. 이는 주변 온도증가에 따른 유막의 점도강하 효과보다, 구조물의 열팽창에 의하여 유막의 강성이 증가하는 효과가 더 크다는 것을 의미한다. Fig. 15(c), (d)는 감쇠계수의 변화를 나타낸다. 앞서 살펴본 강성계수의 변화와 다른 경향을 볼 수 있다. 패드 주변의 등온조건 온도 증가에 따른 감쇠계수는 감소하는 경향을 보인다. 샤프트 주변의 온도 증가는 감쇠계수를 증가시킴을 볼 수 있다. 샤프트 주변의 온도 변화는 유막의 온도증가 보다는 유막두께의 감소를 통해 감쇠계수를 증가시키는 것을 볼 수 있다.

본 연구에서는 저널 지름 105 mm, 회전속도 12,000 rpm, 16 kN의 하중을 지지하는 LBP타입의 5패드 틸팅패드 저널베어링 주변의 온도경계조건 변화에 따른 정특성 및 동특성 변화를 알아보았다. 이는 본 연구에서 제시된 특정 베어링에서 나타나는 현상일 뿐, 다른 형상과 운전조건을 지니는 베어링에 대해서는 다른 경향의 결과를 보일 수 있다.

본 연구에서는 베어링패드와 샤프트 주변의 온도경계 조건의 변화에 따른 베어링의 정특성 및 동특성의 변화에 대하여 살펴보았다. 매개변수 해석을 통하여 제시된 결과를 요약하면 다음과 같다.

(a) 베어링 패드와 샤프트 주변의 열경계조건을 열대류조건 (convection boundary condition)과 등온조건 (constant temperature)로 나누어 그 영향을 각각 살펴보았다.

(b) 베어링 패드와 저널의 정적 변화량은 등온조건이 열대류조건에 비하여 큰 영향을 미쳤다.

(c) 패드의 등온조건 온도의 감소는 최소유막두께 증가 및 패드 온도 감소로 인한 베어링 건전성 확보에 좋은 영향을 미쳤다. 하지만 이는 동력손실량을 증가시키는 원인으로 작용하였다.

(d) 강성계수는 패드와 샤프트 주변 온도의 증가로 인하여 함께 증가하는 경향을 보였다.

(e) 감쇠계수는 패드 주변온도 변화보다 샤프트의 주변온도 변화에 큰 영향을 받았다.

본 연구에서는 베어링 주변의 열경계조건을 예측하기 보다는 열경계조건의 변화에 따른 성능변화의 경향을 제시하고자 한다. 본 연구의 결과는 베어링 동작의 신뢰성을 예측하기를 원하는 베어링 설계자에게 설계지침을 제공할 수 있다.

h: Film thickness μ: Lubricant viscosity U: Shaft surface linear velocity t: Time T: Temperature x: Circumferential position y: Radial position z: Axial position r: Radial position ρ: Density c: Heat capacity k: Heat conductivity Cb: Bearing clearance δtilt: Pad tilt angle Ppvt: Pivot displacement θ: Circumferential nodal position Cp: Pad clearance θp: Pivot angular position R: Journal radius