Analysis of Integration Factor Effect in Dynamic-Structure-Fluid-Heat Coupled Time Transient Staggered Integration Scheme for Morton Effect Analysis

모튼이펙트 해석을 위한 동역학-구조-유체-열전달 시간과도응답 연성해석 시차적분법에서 시상수 효과 분석

Abstract

The present study focuses on the effect of staggered integration factor (SIF) on Morton effect simulation results. The Morton effect is a synchronous rotordynamic instability problem caused by the temperature differential across the journal in fluid film bearings. Convection and conduction of heat in the thin film displaces the hot spot, which is the hottest circumferential position in the thin film, from -20 to 40 degrees ahead of the high spot, where the minimum film clearance is experienced. The temperature differential across the journal causes a bending moment and the corresponding thermal bow in the rotating frame acts like a distributed synchronous excitation in the fixed frame. This thermal bow may cause increased vibrations and continued growth of the synchronous orbit into a limit cycle. The SIF is developed assuming that the response of the rotor-lubricant-bearing dynamic system is much quicker than that of the bearing-journal thermal system, and it is defined as the ratio between the simulation time of the thermal system and the rotor-spinning period. The use of the SIF is unavoidable for efficient computing. The value of the SIF is chosen empirically by the software users as a value between 100 and 400. However, the effect of the SIF on Morton effect simulation results has not been investigated. This research produces simulation results with different values of SIF.

Nomenclature

M_R :회전체 질량 행렬

C_Gyro :회전체 자이로스코픽 댐핑 행렬

C_Brg :베어링 댐핑 행렬

K_Struc :회전체 강성 행렬

K_Brg :베어링 강성 행렬

X_R :회전체 변위 벡터

F_R :회전체 외력 벡터

m_G :베어링 패드 질량

p_pvt :베어링 패드 피봇 방향 변위

δ_tilt :베어링 패드 틸트각도

K_p :피봇 강성

J_G :베어링 패드 질량관성모멘트

F_f :유막이 패드에 가하는 외력

M_O :유막이 패드에 가하는 모멘트

p :유막에서 발생하는 압력

h :유막 두께

t :시간

U :유막 내 저널 원주방향 선속도

μ :유막 점성

μ0 :에서의 점성

β :윤활유 점성 상수

T :유막 온도

T₀ :유막 기준 온도

ρ :윤활유 밀도

c :윤활유 비열

u :윤활유 원주방향 속도

w :윤활유 축방향 속도

x :유막 원주방향 기준좌표

y :유막 두께방향 기준좌표

z :축방향 기준좌표

X :베어링 패드 원주 접선 방향 변위

Y :베어링 패드 반지름 방향 변위

k :유막 열전도계수

T :회전체의 회전 주기

 

1. 서 론

모튼이펙트(Morton effect)는 Fig. 1에서와 같이 유체 베어링에서 저널의 동기 휘돌림(synchronous whirling)발생 시 최소 유막두께(minimum film thickness)가 발생하는 고점(high spot) 부근에서 지속적으로 발생하는 최대 점성전단열(maximum viscous shearing heat)로 인해 발생한다. 이러한 점성전단 열은 원주방향(circum-ferential direction)으로 고점 근처에서 쌓이게 되고 회전하는 저널은 저널비대칭 가열(journal differential heat-ing) 된다. 원주방향으로 비대칭 가열된 저널은 열팽창계수에 따른 비대칭 열굽힘을 일으키고 이러한 열굽힘은 새로운 열편심(thermal bow induced eccentricity)을 일으켜 회전체의 진동을 악화시킨다.

 

Fig. 1. Journal synchronous whirling with high spot.

모튼이펙트의 정밀한 해석을 위해서는 베어링 유막에서 발생하는 힘, 윤활 및 열에 대한 정밀한 모델링이 필수적이며, 이에 따른 열전달, 열변형 및 동특성 변형해석이 필요하다. 또한 모튼이펙트는 시간에 따라 주기적으로 변화하는 진폭의 진동을 수반하므로 THD(thermo-hydro-dynamic) 윤활 모델을 기반으로 한 시간과도응답(time transient) 회전체-베어링 동역학-윤활 수치해석을 수행해야 한다.

본 연구에서는 다중물리 시스템의 시간과도응답 해석을 위한 여러 단계 중 동역학 모델과 열모델의 시차적분 과정에서 사용자가 임의로 또는 경험적으로 선택하게 되어 있는 시차적분 상수(staggered integration co-efficient)가 모튼이펙트 해석 결과에 미치는 영향에 대해 알아보고자 한다.

 

Table 1. Comparison of numerical approaches for Morton effect analysis

Table 1은 모튼이펙트와 관련되어 진행된 선행 연구들[1-7]의 특성을 비교한다. Keogh와 Morton[1]은 모튼이펙트 시스템을 해석적 방법(analytical method)으로 모델링한 최초의 연구이며 저자 중 한 사람인 Morton의 이름을 따라 현재까지도 모튼이펙트라고 부른다. Larsson은[2] Morton의 연구[1]에서 제시한 수치모델을 발전시켰으며 당시에는 다중물리 시스템을 수치적으로 모델링 및 해석하는 데에 한계가 있었으므로 해석적 방법으로 접근하였다.

Gomiaciaga와 Keogh[3]는 회전좌표계에 존재하는 저널과 정지좌표계에 있는 유막 사이의 열 경계조건을 구현하기 위하여 궤도 평균 열유속 경계 조건(orbit aver-aged heat flux boundary condition)을 처음으로 제안하였다. 하지만 해당 연구는 미리 주어진 저널 휘돌림이 존재할 때에 저널 주변에 발생하는 비대칭 가열현상을 수치적으로 보여주는 것에 그쳤다. 저널의 운동은 동역학 방정식이 아닌 미리 주어진 동기궤도를 그리는 것으로 가정하였으며, 베어링의 경우 고정패드 베어링(fixed pad journal bearing)을 사용하였다. 해당 연구는 본 연구에서 수행하는 저널 비대칭가열 해석 모델 개발에 있어비교 및 검증되었다.

Balbahadur와 Kirk[4]는 당시까지 발표된 모튼이펙트관련 실험 결과 예들을 수치적으로 모델링하였으며 시간에 따른 동적 거동을 해석하지 않고, 선형 모델을 도입하여 불안정 회전속도 영역을 예측하고자 하였다. Murphy[5]는 구체적인 수치해석 모델을 제시하지 않고해석적 관점에서 모튼이펙트 현상을 설명하려고 하였다.

하지만 해당 연구는 모튼이펙트의 원리에 대해 설명을 할 뿐 특정 시스템의 불안정 여부에 대해서는 예측이 불가능하다. Lee와 Palazzolo는[6] 모튼이펙트 현상을 연구하기 위하여 최초로 시차적분법(staggered integrationscheme)을 제안하였으나 시차적분상수의 구체적 의미와 영향에 대해 다루지 않았다. Suh와 Palzzolo[7]도 역시모튼이펙트 현상을 해석하기 위하여 시차적분법을 사용하였으며 모튼이펙트의 대표적인 현상 중 하나인 긴 주기를 가지는 진동의 진폭 변화와 나선형 진동을 처음으로 예측하였다. 하지만 해당 연구에서도 시차적분상수의 영향과 의미에 대해 구체적으로 다루지는 않았다.

본 연구는 Suh와 Palzzolo[7]의 연구에서 개발된 수치해석 소프트웨어를 기반으로 수행되었으며, 시차적분상수의 경우 그 값의 설정과 관련하여 선행 연구가 없었다. 본 연구에서는 지금까지 소프트웨어 사용자의 경험에 의해 사용되어 왔던 시차적분 상수의 구체적 의미와 그것이 모튼이펙트 수치해석결과에 미치는 영향에 대해 알아보고자 한다.

 

2. 수치해석 모델

 

2-1. 회전체-베어링 동역학 모델

회전체동역학 해석을 위해서 횡 방향 진동을 고려하였으며, 계산 시간 단축을 위해 모드 좌표 변환이 사용되었다. 모드 해석은 정확도의 손실 없이 자유도 감소를 위한 로터 베어링 시스템 분석에 유용한 도구이다. 시스템 응답은 직교하는 고유 벡터의 선형 조합이라고 가정할 수 있다. 물리좌표계에서 회전체-베어링 계의 운동방정식은 아래 식 (1)과 같이 표현될 수 있다.

\(M_{R} \ddot{X}_{R}+\left(C_{G y r o}+C_{B r g}\right) \dot{X}_{R}+\left(K_{S r u c}+K_{B r g}\right) X_{R}=F_{R}\)     (1)

 

Fig. 2. Dynamic model of titling pad bearing

본 연구에서는 틸팅패드 저널베어링에 지지되는 돌출질량(overhung mass)을 가지는 회전체의 열적 불안정성에 조사하기 위하여 오일러빔(Euler beam) 유한요소모델(finite element method, FEM)을 사용하였다. Fig. 2는 축방향에서 바라본 패드 동역학 모델을 나타낸다. 패드는 단순한 틸트운동(tilting motion) 뿐 아니라 피봇에서 접촉강성에 의한 변위를 고려하였다(Y방향). 피봇에서의 탄성변형을 헤르지안 접촉이론(Hertzian contact theory)을 이용하여 시간에 따라 변하는 비선형 스프링으로 모델링하였으며 패드는 틸트와 로터 지름 방향으로 병진운동을 가지는 2자유도 시스템으로 모델링하였다. 패드의 운동방정식은 식(2), (3)과 같다.

\(\begin{array}{rl} {m_{G} \ddot{p}_{p v t}-e_{G} m_{G} \ddot{\delta}_{t i l i}} & {\cos \left(\delta_{ti l t}+e_{G} m_{G} \dot{\delta}_{ti l t} {^2\sin }\left(\delta_{ti l t}\right)\right)} \\ {} & {+K_{p} p_{p v t}=F_{o}(t)} \end{array}\)       (2)

\(\begin{aligned} J_{G} \ddot{\delta}_{t i l t}+e_{G} m_{G} \cos \delta_{ {t i l t}} \ddot{\mathrm{p}}_{p v t}+e_{G}^{2} m_{G} \cos \delta_{{t i l t}} & {^2\ddot{\delta}_{{t i l t}}} \\ -e_{G}^{2} m_{G}^{2} \sin \delta_{{t i l t}} \cos \delta_{{t i l t}} &=M_{O}(t) \end{aligned}\)       (3)

 

2-2. 열전달-열변형 모델

저널의 열굽힘을 일으키는 온도 구배의 정밀한 해석을 위하여 회전체 축 방향 온도 구배를 고려하였으며 유한요소로 모델링하였다. 그러나, 회전체 열전달-열변형모델을 위하여 축 전체 길이를 고려하면 유한요소 행렬 사이즈가 커져 상당한 계산 시간이 소요된다. 축 열전달 모델의 궁극적인 목적은 베어링 위치 부근에서 축열 굽힘을 정밀하게 예측하는 것이다. 이러한 이유로 베어링 길이의 7배가 샤프트 열전달-열변형 유한요소 모델의 총 길이로 사용되었다[7].

 

2-3. 윤활모델

얇은 유막(thin film)에서의 압력 분포를 구하기 위한 지배 방정식은 Reynolds 방정식이며, 두 상대 이동하는 평면 사이의 두께 h를 가지는 얇은 유막을 모델링한다. 레이놀즈 방정식은 운동방정식과 연속방정식에서 유도되어 압력과 속도 분포를 계산할 수 있다. 운동방정식의 적분은 압력 구배의 관점에서 속도 분포를 제공하며 이 속도 프로파일을 유량 연속 방정식에 대입하면 레이놀즈 방정식을 얻을 수 있다. 유막 두께 방향으로 점성의 변화를 고려한 일반화된 레이놀즈 방정식은 다음 식 (4)와 같이 나타낼 수 있다.

\(\nabla \cdot\left(D_{1} \nabla p\right)+\left(\nabla D_{2}\right) \cdot\left(U_{2}-U_{1}\right)+(\nabla h) \cdot U_{1}+\frac{\partial h}{\partial t}=0\)       (4)

위의 일반화된 레이놀즈 방정식은 원주방향으로 회전하지 않는 베어링과 회전하는 회전축에 대해 다음과 같이 축약된 형태로 다시 쓸 수 있다.

\(\nabla \cdot\left(D_{1} \nabla p\right)+\left(\nabla D_{2}\right) \cdot U+\frac{\partial h}{\partial t}=0\)       (5)

THD 윤활 문제의 경우, 점도는 온도의 함수로 가정할 수 있으며, 상수 D1과 D2는 다음 식 (6), (7)과 같이 표현할 수 있다.

\(D_{1}=\int_{0}^{h} \int_{0}^{z} \frac{\xi}{\mu} d \xi d z-\frac{\int_{0}^{h} \frac{\xi}{\mu} d \xi}{\int_{0}^{h} \frac{1}{\mu} d \xi}\int^{h}_{0} \int_{0}^{z} \frac{1}{\mu} d \xi d z\)       (6)

\(D_{2}=\frac{\int_{0}^{h} \int_{0}^{z} \frac{1}{\mu} d \xi d z}{\int_{0}^{h} \frac{1}{\mu} d \xi}\)       (7)

일반화된 레이놀즈 방정식에서 유막은 층류, 필름 두께 방향의 일정한 압력, 무시할 수 있는 샤프트 곡률 효과, 무시할 수 있는 유체 관성, 일정한 유체 밀도 및 온도 영향 점도를 갖는 것으로 가정한다. 또한, 고체 및유체 경계면에서 비압축성 뉴튼유체 및 슬립이 발생하지 않는 것으로 가정한다. 점도와 유체의 온도 사이의 관계는 다음과 같이 표현할 수 있다[7].

\(\mu=\mu_{0} e^{-\beta\left(T-T_{0}\right)}\)       (8)

얇은 유체 필름의 온도 분포는 열 경계 조건, 압력, 속도 및 점도 분포가 필요한 에너지 방정식에 의해 결정된다. 레이놀즈 방정식과 에너지 방정식은 점도 프로파일을 통해 결합된다. 업데이트된 점도 분포는 레이놀즈 방정식의 해를 구하는 데에 사용된다. 층류, 비압축성 및 뉴톤 유체의 경우, 에너지 방정식은 다음 식 (9)와 같다.

\(\begin{array}{l} {\rho c\left(u \frac{\partial T}{\partial x}+w \frac{\partial T}{\partial z}\right)} \\ {=k\left(\frac{\partial^{2} T}{\partial x^{2}}+\frac{\partial^{2} T}{\partial y^{2}}+\frac{\partial^{2} T}{\partial z^{2}}+\mu\left[\left(\frac{\partial u}{\partial y}\right)^{2}+\left(\frac{\partial w}{\partial y}\right)^{2}\right]\right.} \end{array}\)       (9)

여기서 열전도는 x, y 및 z 방향으로 고려되며, 대류항은 레이놀즈 방정식에 의해 생성된 속도 프로파일(u,w)이 x 및 z 성분만을 갖기 때문에 x 및 z 방향으로 만고려된다. 점성전단항은 박막 두께가 단지 막 두께 방향이기 때문에 (∂u/∂y)²항 및 (∂w/∂y)²만으로 고려된다.

 

2-4. 모튼이펙트 해석 알고리즘

모튼이펙트 현상은 회전체-베어링 동역학, 유막에서의 열발생-전달, 베어링-회전체 구조물에서의 열전달 및 열변형, 그리고 이러한 현상들 간의 상호관계를 해석해야 하는 다중물리 문제이다. Fig. 3은 모튼이펙트 현상을 해석하기 위한 기본 알고리즘을 나타내며 각각의 블록에서 일어나는 연산 과정은 아래와 같다.

(1) 먼저 각 물리시스템에 대한 초기값이 필요하다. (베어링 패드 초기온도 Y_^OThermPad , 저널 초기온도 Y_^OThermShaf , 로터-베어링 초기변위 Y_^ODyn , 초기 유체점성 μ_O )

(2) 회전체-베어링 동역학 및 윤활 시스템을 Matlabode function을 이용하여 시간과도응답 해석한다. 본 시스템 안에서는 베어링 동역학-회전체동역학-레이놀즈 방정식이 같은 시간 영역에서 동시 해석된다. 레이놀즈 방정식에서 계산된 힘은 베어링 패드 및 회전체 저널부분에 외력을 가하게 되고 이들을 수치적분한다. 본 회전체-베어링 동역학 문제를 해석하는 자세한 과정은 아래와 같다.

(2-1) 로터-패드의 상태변수(Y_Dyn ), 베어링 패드 열변형(h_T.E.P ), 저널의 열변형(h_T.E.J ) 정보를 기반으로 유막두께및 유막두께의 시간변화량을 계산한다.

(2-2) 시스템에 저장되어 있는 유막에서의 점성 및(2-1)에서 계산된 유막두께를 기반으로 일반화된 레이놀즈 방정식을 이용하여 유막에서 발생하는 힘 및 윤활유의 속도를 계산한다.

(2-3) 유막에서 발생된 유압, 로터의 길이방향으로 주어진 편심력(eccentric force) 그리고 로터의 열굽힘에 기인한 편심력이 회전체에 작용하고 베어링 패드의 경우 유막만이 외력으로 작용한다.

(2-4) (2-5)에서 계산될 에너지 방정식의 전 단계로써베어링 패드 사이에 공급되는 윤활유 온도 및 바로 앞의패드에서 나오는 윤활유의 양 및 온도를 기반으로 뒤쪽패드에 들어가는 윤활유의 온도를 계산하며 본 온도는 에너지방정식에서 경계조건으로 사용된다.

(2-5) 에너지방정식에서는 윤활유의 속도, 유막두께, 점성, 패드온도, 샤프트 온도 및 입력단에서의 온도를 초기 및 경계조건으로 이용하여 점성전단 열발생 및 열전달을 계산한다.

(2-6) 본 회전체-베어링 동역학 시스템이 수렴될 때까지 본 (2)번에 포한된 과정을 반복한다.

(3-1) (2-5)에서 계산된 유막에서의 온도와 현재의 베어링 패드 온도를 기반으로 유막과 패드 사이의 열경계조건을 계산한다.

(3-2) (3-1)에서 계산된 열경계조건을 이용하여 패드의 온도를 시간과도응답법을 이용하여 계산한다.

(3-3) (3-2)에서 계산된 패드 온도구배(temperaturedistirbution)를 이용하여 패드의 열변형을 계산한다. 이때의 열변형량은 반지름방향으로의 변위벡터(h_T.E.P )로 계산되어 유막두께(h)에 사용된다.

(4-1) 베어링 패드에서와 마찬가지로 (2-5)에서 계산된 유막에서의 온도계산결과와 저널의 온도구배 정보를 바탕으로 유막과 회전하는 저널 사이의 열경계조건을 만들어낸다. 이 때 유막의 원주방향(circumferential dir-ection) 위치는 고정되어 있으나 저널이 회전하게 되어 (3-1)에서 사용된 온도경계조건을 사용할 수 없다. 이 때에는 Keoh와 Gormiciaga[3]가 제안한 궤도 평균 열유속경계 조건을 사용하여 회전하는 저널 주변의 온도조건을 계산한다. 단, 로터가 동기궤도(synchronous orbiting)를 그린다고 가정하며 수치해석 코드에서 로터가 동기궤도를 그리는지 여부에 대해 매 시간 단계마다 확인한다.

(4-2) (4-1)에서 계산된 열경계조건을 이용하여 저널에서의 저널의 온도구배를 계산한다.

(4-3) (4-2)에서 계산된 저널 온도구배를 기반으로 로터의 베어링 부분에서의 열변형 및 열굽힘량을 계산한다. 열변형량은 유막두께를 변화시키며, 열굽힘량의 경우 열편심량을 유발한다.

 

Fig. 3. Algorithm for Morton effect analysis.

 

2-5. 시차적분법 및 상수

 

Fig. 4. Schematic diagram for staggered integration scheme.

본 연구에서 제시된 모튼이펙트 해석 알고리즘에서 시간에 대한 적분과정은 몇가지 가정을 가지고 있다. 이러한 가정은 동역학-윤활-열전달-열변형 특성을 고려해야 하는 다중물리 시스템의 시간과도응답(변위, 온도, 열변형)을 구하고자 할 때 해석 시간을 줄이기 위해 고려되어야 한다. 그 중 시스템의 해석 결과에 가장 큰 영향을 미치는 가정은 동역학-윤활 시스템과 패드-저널 열전달-열변형 시스템 사이에서 고려된다.

Fig. 4는 본 연구에서 다루고자 하는 시차적분법에 대한 개념도이다. T는 동역학-윤활 시스템에서 저널이 한 바퀴 회전하는 시간 즉 회전주기를 의미한다. 실제 회전체-베어링-윤활 시스템에서는 동역학-윤활 시스템에서의 여러 상호작용은 물론이고 유막과 직접 맞닿아 있는 저널 및 베어링과의 열전달이 동시간대에 이루어져야 한다. 또한 열전달로 인한 열구배-열변형이 동시에 이루어진다. 하지만 이러한 현상을 동시간 대에서 수치적으로 적분하기 위해서는 상당한 시간이 소요되며 특히 긴 시간(최소 2분에서 길게는 10분 이상)을 해석해야 하는 모튼이펙트의 경우는 불가능하다. 현존하는 어떤 상용 다중물리 소프트웨어도 모튼이펙트를 다중물리 관점에서 해석한 예는 없으며 본 연구에서는 상용 수치해석 소프트웨어인 Matlab을 이용하여 모델링하였다. 동역학-윤활 시스템이 정상상태에 도달하기 위해서 폐쇄 루프(closed loop)는 m번 반복된다(Fig. 3 단계2) 이 때 정상상태에 도달한 회전체-윤활-베어링 동역학 시스템의 유막에 분포되어 있는 온도분포를 저장한다.

동역학 시스템에서 시스템이 정상상태에 도달한 m번째 루프에서 저장된 마지막 한 주기 동안의 유막 내 온도분포 데이터는 시스템 메모리에 저장되고, 저널-패드간의 열전달 경계조건으로 사용되며 본 경계조건을 이용하여 T·SIF시간 동안 열시스템의 시간과도응답을 해석한다. 본 연구에서는 SIF의 변화에 따른 모튼이펙트다중물리 해석 결과의 변화에 대해 고찰해보고자 한다.

Runge-Kutta 법과 같은 일반적인 수치적분에서도 시간단계(time step, Δt)가 작을수록 더 정확한 값에 수렴한다는 사실을 통해 작은 시차적분상수(n)를 사용하는 것이 실제와 가까운 현상 해석이 가능하다고 예상할 수 있다. 하지만 모든 수치해석에서 그렇듯이 효율적인 해석과 정밀한 해석 사이에서는 반드시 타협점이 존재해야 한다.

 

2-6. 회전체-베어링 모델

 

Table 2. Bearing configuration [7]

 

Table 3. Lubricant properties [7]

 

Fig. 5. Rotordynanic Model [7]

본 연구에서는 dejong [7]이 연구한 로터-베어링 시스템을 해석하였다. 입력변수는 Table 2~3에 제시되어 있으며, 로터의 형상은 Fig. 5와 같다. 주목할 점은 비작동측(non-driven end, N.D.E) 베어링 우측에 돌출질량(overhung mass)이 존재한다는 것이다. 일반적으로 모튼이펙트는 비작동측 베어링 바깥쪽에 돌출질량을 가지는 회전체 시스템에서 주로 발생한다고 알려져 있으며 있다. 이는 저널베어링의 열굽힘으로 인한 큰 열편심력을 일으키기 위해서는 편심위치에 무거운 질량이 위치해 있어야 하기 때문이다.

 

3. 결과 및 고찰

본 연구에서는 6,500 rpm 및 7,000 rpm 두 회전속도영역에서의 회전체 동적 거동을 해석하였다. 본 연구에서는 100, 200, 300, 400 네 가지의 시차적분 상수(n)를 사용하였을 때의 결과들을 비교하였다. Fig. 6은 비작동측 베어링에서의 로터 진폭의 시간에 따른 변화량을 나타내며, Fig. 7은 시간에 따른 저널 비대칭 온도(journaltemperature differential, ΔT)를 나타낸다. 6,500 rpm의 경우 모튼이펙트가 발생하지 않고 안정적인 경우이며, 약6분 후 진폭이 수렴하고 저널 비대칭 온도도 약 1.75° C에 수렴하는 것을 볼 수 있다. 7,000 rpm의 경우, 진폭이 수렴하지 않고 약 3분 주기로 진동하는 것을 볼 수 있으며, 저널 비대칭 온도 또한 같은 주기로 진동하는 것을 볼 수 있다. 이러한 저널 비대칭 가열 현상은 Fig. 8에서 보는 바와 같이 베어링-유막-저널의 단면 온도분포를 통해서도 확인 가능하다. 유막의 경우 온도 분포를 쉽게 확인하기 위하여 실제 두께보다 과장되게 표현하였다. Fig. 8(a)는 로터 회전속도가 6,500 rpm 및 해석시간 7분 42초일 때의 단면 온도분포를 나타내며, Fig. 6(a)에서와 같이 저널 비대칭온도가 약 1.75o C이다. 로터 회전속도가 7,000 rpm이며 해석 시간 4분 18초일경우 저널 비대칭온도가 약 12o C로 크며 이러한 온도분포가 회전체 진동을 악화시켰음을 알 수 있다.

 

Fig. 6. Vibration amplitude at N.D.E. bearing with staggered integration coefficient (pk-pk)-Morton case.

 

Fig. 7. Journal temperature differential

Fig. 6과 Fig. 7을 통해 시차적분상수 값의 변화에 따른 동적 거동의 변화가 눈에 띄게 달라지는 것을 확인할 수 있다. 예상한 것과 같이 작은 값의 시차적분상수를 사용할수록 더 빠른 응답을 보이는 것을 확인할 수 있다. 하지만 단순히 시스템의 응답만 지연되는 것이 아니라 시스템의 진폭 또한 변하는 것을 확인할 수 있다. 이는 1자유도 진동시스템의 시간응답을 수치해석할 때에 작은 시간단계를 사용할수록 응답이 빨라지고 진폭이 줄어드는 것과 비슷한 현상이다. 하지만 모튼이펙트와 같이 여러 물리시스템을 시간과도응답법으로 해석해야 하는 경우 무작정 작은 값을 선택하는 것은 때에 따라 불가능하거나 비효율적일 수 있다. 이에 대한 최적의 값을 찾는 것은 정확하고 효율적인 해석을 위해 필수적이다. 일반적인 모튼이펙트 해석에 있어 시차적분상수 값이 100이고, 로터 회전속도가 6,500 rpm 경우, 인텔 제온8코어 cpu, 64G 메모리, 그리고 Matlab 병렬처리 알고리즘 사용시 8분 해석을 위해 약 1주일이 소요되었다.

발생 원인과 해결책이 명확히 알려져 있지 않고, 고속회전기계에서 주로 발생하는 것으로 알려진 모튼이펙트의 해결책 제시를 위해서는 현재보다 더 정밀한 해석이 요구된다. 이에 더 작은 시차적분상수 사용이 필요하지만 해석의 경제성과 정밀도 사이의 균형이 필요하다. 본 연구에서는 지금까지 명확한 기준 없이 사용자의 경험에 따라 사용해 온 시차적분상수 값의 변화에 따라 다중물리 회전체-베어링 동역학 시스템의 해석결과가 영향을 미친다는 것을 확인하였다. 이에 후속 연구에서는 최적의 시차적분상수값을 찾기 위한 추가적인 수치모델및 이에 대한 고찰이 필요하다.

 

Fig. 8. Temperature distribution (n=100).

 

4. 결 론

큰 시상수를 사용하는 것은 빠른 해석의 장점이 있지만 모튼이펙트와 같이 다중물리 시스템의 시간에 따른 상태 변화량이 큰 경우 정밀한 예측이 어렵다는 것을 확인할 수 있었다. 한편으로 작은 시상수의 선택은 해석시간을 길게 만들어 연구의 효율성을 떨어뜨린다. 결국 해석결과에 큰 영향을 미치지 않는 범위 내에서 최적의 시상수 선택이 필요하다.

지금까지 연구자의 경험에 의해 선택되어 온 시차적분 상수가 모튼이펙트 해석결과에 미치는 영향에 대해 살펴본 것은 본 연구가 처음이며, 결코 연구자의 경험이 아닌 물리적 타당성을 바탕으로 선택되어야 한다.

본 연구를 바탕으로 다중물리 시스템의 시간과도응답해석에서 각 시스템 별 시상수(time constant)와의 관계를 바탕으로 최적의 시상수의 선택법에 대해 연구가 추가로 필요하다.