Comparison of Approximation and Rotordynamics Solutions for Design of a High Speed Air Spindle

고속 공기 스핀들 설계를 위한 근사해석과 회전체동역학의 비교

Abstract

This paper presents two methods for designing a high-speed air spindle operated over the rotational speed of 50,000 rpm. The first method is an approximate method, which assumes a symmetric spindle shape even though it is not symmetric in reality. The second is an analysis of rotordynamics using beam and solid models. The approximate method can be used to calculate the bearing load capacities, stiffness and damping coefficients, stability of the shaft system, and response of the forced excitation from the unbalanced mass. Designers can use this method to determine the dimensions of the desired spindle at the first stage of the design. The more detailed behavior of the spindle can be calculated using the rotordynamics theory using beam and solid models based on the Finite Element Method. In this paper, a spindle, with two air bearings, one motor at the end, and two air thrust bearings, is newly developed. The solutions from the two rotordynamics theories are compared with the solution obtained using the approximate method. The three calculations are in agreement, and the procedure for the design of a spindle system, supported on the externally pressurized air bearings, ispresented and discussed.

Nomenclature

A_s : Area of restrictor (m2 ) (자성형 저항 면적)

c : Clearance of bearing (m) (베어링 틈새)

C_ij : Damping coefficients (N/(m/s)) (감쇠계수)

C_d : Discharge Coefficients (확산계수)

H : Dimensionless Height (1) (무차원 공기막 두께)

I_T : Lateral mass moment of inertia (kg·m2 ) (횡방향 질량관성 모멘트)

I_P : Polar mass moment of inertia (kg·m2 ) (극관성모멘트)

K_ij : Stiffness coefficients (N/m) (강성계수)

P : Dimensionless pressure (1) (무차원 압력)

p_a : Ambient pressure (N/m2 ) (대기압)

Q : Dimensionless mass flow rate (무차원 유량)

r : Radius of bearing (m) (베어링 반경)

R : Gas constant (J/(mol·K)) (기체상수)

S_p : Length between mass center and bearing center (m) (베어링과 무게중심사이의 거리)

T : Absolute temperature (K) (절대온도)

X : Dimensionless coordinate along circumference(1) (원주방향 무차원 좌표)

Z : Dimensionless coordinate along shaft axis (1) (축방향 무차원 좌표)

α : Tilting along x axis (rad) (x 축 기울기)

β : Tilting along y axis (rad) (y 축 기울기)

κ : Adiabatic constant (단열지수)

Λ : Dimensionless bearing number (무차원 베어링 수)

μ : Absolute viscosity (Pa·s) (절대점성)

ω : Rotational Velocity (rad/s) (회전속도)

σ : Dimensionless squeeze number (무차원 스퀴즈 수)

 

1. 서 론

최근, OLED 나 휴대폰 액정 등에 사용되는 박막 유리나, 반도체 공정의 웨이퍼 절단 등에 사용되는 공작기계용 스핀들은 정밀가공을 위하여 고속 회전을 요구 받고 있다. 고속 회전을 하는 스핀들은 경제성을 고려하여 공기베어링을 주로 이용한 주축이 사용되고 있다.[1-4]

공기베어링에 관련한 기존의 연구는 베어링의 하중지지용량, 강성, 감쇠등의 베어링 특성수를 계산하는 연구가 주를 이루고 있다[5]. 하지만 스핀들에 적용하기 위해서는 베어링의 특성수 뿐만 아니라, 회전축과 연계된 회전축 진동의 연구가 함께 이루어져야 한다. 송[6] 등은 semi-implicit방법을 이용하여 회전축의 진동 특성을 정밀도와 연결하여 해석한 바 있다. 그 외에도 다양한 문헌들이 스핀들에 적용되어 해석 되어져 왔다[7-10].

하지만, 이전의 연구들은 주축의 형상이 모두 설계된 이후에, 해석을 진행하고 올바른 설계가 되었는지를 검증하는데 그치고 있기 때문에, 설계자들은 초기 설계를 진행하는데 어려움을 겪고 있다. 본 연구에서는 이러한 어려움을 극복하고자, 초기 설계를 위한 회전체 진동 근사해석 방법을 제시하고, 이를 통해 1차 고유진동수에 해당하는 위험속도, 베어링의 자려진동이 발생하는 불안정속도, 및 진폭 등을 계산하였다. 이러한 근사해석 방법과 회전체동역학을 적용한 수치해석과의 비교를 통해서 그 타당성을 확인해 보고자 한다.

 

2. 근사해석 방법

 

2-1. 해석 대상 스핀들

본 연구에서 해석 대상인 공기 베어링으로 지지된 고속 스핀들의 개략도를 Fig. 1에 나타내었다. 그림에서 보듯이, 가능한 두 저널 베어링이 전체 축의 대칭과 가깝게 설계 하기 위해서 무게가 가벼운 모터의 회전자를 두베어링의 왼쪽에, 스러스트 판을 두 베어링의 오른쪽에 배치하여 베어링과 베어링 사이에 무게중심이 위치할 수 있는 설계를 하였다. 이렇게 설계하는 이유는 대칭 구조가 됨으로써, 간이 계산을 통해서 초기 설계를 쉽게 할수 있는 장점을 활용하기 위해서이다.

 

Fig. 1. Schematics of air spindle.

 

2-2. 외부가압 공기베어링의 특성수

본 연구에서 사용한 외부가압 공기베어링은 공기 햄머링 현상이 나타나지 않도록 포켓이 없고 작은 구멍만으로 구성된 자성형 급기공을 갖는 베어링으로 선택하였다. 이 베어링의 형상에 관련한 좌표계는 Fig. 2에 나타내었다. 자성형 급기공을 통과하는 유량을 고려하여 무차원 레이놀즈 방정식을 나타내면 다음과 같다.

\(\frac{\partial}{\partial X}\left(P H^{3} \frac{\partial P}{\partial X}\right)+\frac{\partial}{\partial Z}\left(P H^{3} \frac{\partial P}{\partial Z}\right)+Q_{s i j}=\Lambda \frac{\partial(P H)}{\partial X}+\sigma \frac{\partial(P H)}{\partial \tau}\)       (1)

여기서 사용된 무차원 수들은 다음과 같다.

\(\Lambda=\frac{6 \mu \omega}{P_{a}}\left(\frac{r}{c}\right)^{2}, \sigma=\frac{12 \mu \omega}{P_{a}}\left(\frac{r}{c}\right)^{2}, Q_{s i j}=\frac{12 \mu r^{2} R T q_{s i j}}{A_{s} p_{a}^{2} c^{3}}\)       (2)

급기공이 없는 위치에서는 구멍으로 들어오는 유량이 0 으로 입력되며, 급기공이 있는 경우에만 유량이 더해지게 된다. 급기공을 통해 들어오는 유량은 등엔트로피 조건을 가정하여 다음과 같은 유량식을 적용하였다[11,12].

\(\begin{aligned} &q_{s i j}=c_{d} \pi d_{s} h \frac{p_{s}}{\sqrt{R T}}\left(\frac{2 \kappa}{\kappa+1}\right)^{\frac{1}{2}}\left(\frac{2}{\kappa+1}\right)^{\frac{1}{\kappa-1}}, \frac{p_{i j}}{p_{s}} \leq\left(\frac{2}{\kappa+1}\right)^{\frac{1}{\kappa-1}}\\ &=c_{d} \pi d_{s} h \frac{p_{s}}{\sqrt{R T}}\left(\frac{2 \kappa}{\kappa-1}\right)^{\frac{1}{2}} \sqrt{\left(\frac{p_{i}}{p_{s}}\right)^{\frac{2}{\kappa}}-\left(\frac{p_{i l}}{p_{s}}\right)^{\frac{\kappa+1}{\kappa}}}, \frac{p_{i j}}{p_{s}}>\left(\frac{2}{\kappa+1}\right)^{\frac{1}{\kappa-1}} \end{aligned}\)       (3)

 

Fig. 2. Configuration of an externally pressurized air bearing.

식 (3)에서 사용한 확산계수 (Discharge coefficient)는 0.6이다[13,14]. 본 연구에서는 상기 레이놀즈 방정식의 해를 구함에 있어서 유한요소해석법을 기반으로 한 상용 해석프로그램인 COMSOL 을 해를 구하기 위한 도구로 사용하였기 때문에 수치해석 기법에 대해서는 설명하지 않는다. 일반적으로 공기 스핀들을 해석하기 위해서는 베어링이 가지고 있는 강성 및 감쇠계수를 구하는 것이 매우 중요하다. 대부분의 선행 연구에서는 강성 및 감쇠계수를 섭동법에 의해 구해왔지만[1], 본 연구에서는 강성 계수와 감쇠계수를 정상상태의 하중과 그 위치에서 미소 변위와 미소 속도를 직접 베어링 틈새에 가함으로써 Fig. 3에 표시된 좌표계의 정의대로 강성 및 감쇠계수를 구하였다.

\(\begin{aligned} &\frac{\Delta F_{x}}{\Delta x}=K_{x x} \frac{\Delta F_{x}}{\Delta y}=K_{x y} \frac{\Delta F_{y}}{\Delta x}=K_{y x} \frac{\Delta F_{y}}{\Delta y}=K_{y y}\\ &\frac{\Delta F_{x}}{\Delta \dot{x}}=C_{x x} \frac{\Delta F_{x}}{\Delta \dot{y}}=C_{x y} \frac{\Delta F_{y}}{\Delta \dot{x}}=C_{y x} \frac{\Delta F_{y}}{\Delta \dot{y}}=C_{y y} \end{aligned}\)       (4)

 

Fig. 3. Stiffness and damping coefficients.

 

2-3. 회전축 진동의 고유진동수 근사해

공기 스핀들의 고유진동수를 구하기 위해, 우선 감쇠는 없으며 주강성만 고려하여 근사해를 구한다. 이렇게 가정할 수 있는 이유는 공기베어링은 불안정 속도에 도달하기 전 고유진동수 부근까지는 주강성이 연성강성에 비해 매우 크며, 감쇠는 매우 작기 때문이다. Fig. 1에 나타낸 바와 같이 x 축 방향으로의 기울어짐을 α, y 축 방향으로의 기울어짐을 β, 베어링의 주강성은 k라 할 때 계의 운동방정식은 다음과 같이 표현된다.

\(\left[\begin{array}{cccc} {m} & {0} & {0} & {0} \\ {0} & {m} & {0} & {0} \\ {0} & {0} & {I_{T}} & {0} \\ {0} & {0} & {0} & {I_{I}} \end{array}\right] \cdot\left[\begin{array}{c} {\ddot{x}} \\ {\ddot{y}} \\ {\ddot{\alpha}} \\ {\ddot{\beta}} \end{array}\right] + \left[\begin{array}{ccc} {000} & {0} & {0} \\ {00} & {0} & {0} \\ {00} & {0} & {-I_{p} \cdot \Omega} \\ {0} & {0 I_{p} \cdot \Omega} & {0} \end{array}\right] \cdot\left[\begin{array}{c} {\dot{x}} \\ {\dot{y}} \\ {\dot{\alpha}} \\ {\dot{\beta}} \end{array}\right] \\+ \left[\begin{array}{cccc} {2 k} & {0} & {0} & {0} \\ {0} & {2 k} & {0} & {0} \\ {0} & {0} & {\frac{1}{2} \mathrm{kL}^{2}} & {0} \\ {0} & {0} & {0} & {\frac{1}{2} \mathrm{kL}^{2}} \end{array}\right] \cdot\left[\begin{array}{c} {\sqrt{x}} \\ {y} \\ {\alpha} \\ {\beta} \end{array}\right]=0\)       (5)

이 계의 운동방정식을 풀면 고유진동수를 전방 휘돌림과 후방 휘돌림으로 나누어서 다음과 같이 구할 수 있다.

1) 전방 휘돌림 \(\frac{\Omega}{\omega}=1\)

\(\omega_{1}=\sqrt{\frac{2 k}{m}}, \omega_{2}=\sqrt{\frac{k L^{2}}{\left(I_{T}-I_{P}\right) \cdot 2}}\)       (6) 

2) 후방 휘돌림 \(\frac{\Omega}{\omega}=-1 \)

\(\omega_{1}=\sqrt{\frac{2 k}{m}}, \omega_{2}=\sqrt{\frac{k L^{2}}{\left(I_{T}+I_{P}\right) \cdot 2}}\)       (7)

축에 남아있는 불평형량에 의한 진동은 전방휘돌림에 대해서는 가진을 하며, 후방 휘돌림은 진폭을 줄이는 방향으로 작동하므로 고유진동수를 고려할 때는 전방 휘돌림만 고려하여 설계하면 된다.

 

2-4. 강제 진동에 의한 스핀들 진폭 예측

스핀들의 경우에는 고유진동수뿐만 아니라 공작기계의 회전 정밀도를 유지해야 하기 때문에, 스핀들의 진폭을 고려해서 설계해야 한다. 본 절에서는 스핀들의 강제진동 거동을 예측하기 위해서 연성강성 및 연성감쇠까지도 고려한 운동방정식에 대해서 설명한다. 밸런싱 작업을 하더라도, 축에는 불평형 질량이 남아 있게 되므로, 본 연구에서는 축의 불평형 질량을 고려한 다음과 같은 운동방정식을 사용하였다.

\(m \ddot{x}+2 c_{x x} \dot{x}+2 c_{x y} \dot{y}+2 k_{x x} x+2 k_{x y} y=m u \Omega^{2} \cos (\Omega t+\phi)\\m \ddot{y}+2 c_{y x} \dot{x}+2 c_{y y} \dot{y}+2 k_{y x} x+2 k_{y y} y=m u \Omega^{2} \sin (\Omega t+\phi)\\I_{T} \ddot{\alpha}+2\left(\frac{L}{2}\right)^{2} c_{y y} \dot{\alpha}+2\left(\frac{L}{2}\right)^{2} c_{y x} \dot{\beta}-I_{p} \cdot \Omega \cdot \dot{\beta}+2\left(\frac{L}{2}\right)^{2} k_{y y} \alpha\\+2\left(\frac{L}{2}\right)^{2} k_{y x}=-m u \Omega^{2} \rho \sin (\Omega t+\phi)\\I_{T} \ddot{\beta}+2\left(\frac{L}{2}\right)^{2} c_{x x} \dot{\beta}+2\left(\frac{L}{2}\right)^{2} c_{x y} \dot{\alpha}+I_{p} \cdot \Omega \cdot \dot{\alpha}+2\left(\frac{L}{2}\right)^{2} k_{x x} \beta\\+2\left(\frac{L}{2}\right)^{2} k_{x y} \alpha=m u \Omega^{2} \rho \cos (\Omega t+\phi)\)       (8)

식(8)의 해를 조화함수라 가정하면 아래와 같이 해를 유추할 수 있다.

\(\begin{aligned} &x=a \cdot \cos (\omega t)+b \cdot \sin (\omega t)\\ &y=c \cdot \cos (\omega t)+d \cdot \sin (\omega t)\\ &\alpha=e \cdot \cos (\omega t)+f \cdot \sin (\omega t)\\ &\beta=g \cdot \cos (\omega t)+h \cdot \sin (\omega t) \end{aligned}\)       (9)

무게 중심으로부터 베어링 중심까지의 거리를 Sp라 하고, 식(9)를 식(8)에 대입하여 정리하면 축의 진폭과 위상을 다음과 같이 구할 수 있다.

\(R_{\max }=\sqrt{\frac{A+B}{2}+\sqrt{D^{2}+C^{2}}} \gamma=\tan ^{-1}\left(\frac{D}{C}\right)\)       (10)

여기서 A,B,C,D는 다음과 같다.

\(\begin{aligned} &A=a^{2}+2 S_{p} a g+S_{p}^{2} g^{2}+c^{2}+2 S_{p} c e+S_{p}^{2} e^{2}\\ &B=b^{2}+2 S_{p} b h+S_{p}^{2} h^{2}+d^{2}+2 S_{p} d f+S_{p}^{2} f^{2}\\ &C=a b+S_{p}(a h+b g)+S_{p}^{2} h g+c d+S_{p}(c f+d e)+S_{p}^{2} e f \end{aligned}\\D=\frac{A-B}{2}\)       (11)

 

3. 근사해에 의한 결과

본 연구에서는 근사해와 회전체동역학을 이용한 수치해석 해의 비교를 목표로 하기 때문에, 다양한 해석 보다는 한 모델을 정하고 해석을 진행하였다. 본 연구에서 사용한 베어링 및 축의 제원은 Table 1에 나타내었다.

 

Table 1. Specification of the Spindle

우선 베어링의 압력분포가 계산이 되는지 확인을 위하여 베어링의 압력분포를 계산해 보았다. Table 1에 표기된 제원으로 압력분포를 계산한 결과 예를 Fig. 4에 나타내었다. 급기공 부근에는 압력이 높고 급기공에서 멀어질수록 압력이 줄어드는 전형적인 외부가압 공기베어링의 압력 분포가 잘 계산됨을 확인 하였다. 또한 계산된 베어링 강성 및 감쇠계수를 Fig. 5에 나타내었다. 다만 공기베어링은 가진 주파수에 따라 강성 및 감쇠계수가 달라질 수 있지만, 본 연구에서는 회전수에 동기된 성분이 가진 주파수라 가정하여 계산하였다.

 

Fig. 4. Example of bearing pressure.

 

Fig. 5. Dynamic coefficients of bearing.

식 (6)을 이용하여 축의 고유진동수를 예측한 결과를 Fig. 6(a)에 나타내었다. 일반적으로 베어링 지지점이 긴 경우에는 평행 모드(cylindrical mode)가 먼저 나타나고코니컬 모드(conical mode)에 의한 고유진동수가 나중에 나타난다. 하지만 베어링의 지지점이 짧아지게 되면 코니컬 모드가 먼저 나올 수도 있다. 본 연구에서는 베어링 뒤쪽에 모터가 들어가 있어 축은 길지만, 베어링 사이의 간극이 전체 길이에 비해 짧기 때문에 코니컬 모드에 의한 고유진동수가 82,000 rpm에서 먼저 나오고, 100,000 rpm 에서 평행 모드가 나오는 것으로 계산되었다.

 

Fig. 6. Results from Approximate method.

Fig. 6(b)에는 식(10) 에 의해서 구해진 축의 강제진동에 의한 응답을 나타내었다. 이 결과에서 살펴보면 코니컬 모드에 의한 변위는 매우 작아서 진폭을 키우지는 않는 것으로 판단되며, 100,000 rpm에서 가장 진폭이 크게 나타나는 것으로 계산된 것은 평행 모드에 해당하는 결과로 판단된다.

 

4. 회전체동역학 해석

본 장에서는 근사해가 아닌 수치해석적 기법에 의해 해를 구하고, 이를 근사해와 비교한다. 우선 회전체동역학은 로터를 보(beam)로 가정하여, 직선으로 모델링하고, 각 절점에 직경, 관성 등의 정보를 입력하여 해석하는 방법과 축 전체를 솔리드 요소로 모델링하여 3차원 형상을 모델링 하는 방법이 있다. 본 연구에서는 2가지 모두를 수행하고 근사해와 비교한다. Fig. 7에는 축의 절점을 9개로 나누어 보로 모델링한 축의 형상을 나타내었으며, 보 모델을 이용하여 계산된 Campbell 선도를 Fig. 8에 나타내었다. Fig. 8의 결과에서 보듯이 보 모델 해석을 통해 구한 위험속도가 코니컬 모드에서는 80,000 rpm, 원통 모드에서는 약 96,000 rpm 정도 나오는 것으로 계산되었으며, 이는 Fig. 6에 나타낸 근사해 보다 약 2.5% 정도만 낮은 결과로 매우 근사한 결과라 볼 수 있다.

 

Fig. 7. Stations for beam model

 

Fig. 8. Campbell diagram from Beam Model.

Fig. 9과 Fig. 10에는 솔리드 모델링을 하고, 베어링의 해석과 연동하여 다중 물리 해석을 한 예 및 Campbell 선도를 나타내었다. 솔리드 모델링의 경우에는 코니컬 모드가 80,000 rpm에서 발생하며 원통 모드는 95,000 rpm에서 발생하는 것으로 예측되었다. 이는 보 모델과 약 1000 rpm 만이 다른 결과이며, 거의 비슷한 결과를 나타내지만, Solid 모델링의 경우 해석시간이 매우 오래 걸리기 때문에 설계시에는 단순 계산을 통한 설계를 진행한 후에, 보 모델이나 솔리드 모델로 한번 더 확인 하는 절차만 거치면 충분히 설계 가능하다고 할 수 있다.

 

Fig. 9. Solid Model.

 

Fig. 10. Campbell diagram from solid modeling.

 

5. 결 론

본 연구에서는 공기베어링으로 지지된 고속 스핀들을 설계하기 위해, 초기 설계를 위한 근사식을 제시하고, 이 근사식에 대한 결과와 보 가정에 의한 회전체동역학 모델링 및 솔리드 가정에 의한 회전체동역학 모델링을 통한 결과를 비교하여 다음과 같은 결론을 얻었다.

1) 본 연구에서 제시한 근사식은 보 및 솔리드 모델과 일치하며 고유진동수는 거의 같은 값을 나타낸다.

2) 강제진동응답 및 고유진동수 해석은 근사해석과 수치해석이 거의 일치하는 결과를 나타낸다.

3) 솔리드 모델링의 경우에는 다른 모델에 비해 고유진동수가 조금 일찍 발생하는 것으로 해석되며 가장 정확한 해를 줄 것으로 예상된다.

4) 해석 시간의 효율을 생각해서 설계시에는 본 연구에서 제시한 근사식으로 설계를 진행하고, 보 또는 솔리드 모델링으로 검증하여 최종 설계를 마무리하는 과정으로 설계하는 것을 제안할 수 있다.