Finite Element Analysis of High-speed Rotating Disks Considering Impulsive Loading by the Clearance and Contact

간격 및 접촉에 의한 충격하중을 고려한 고속 회전 디스크의 유한요소 해석

Abstract

For the time integration solution of the impulsive dynamic contact problem of high-speed rotating disks formulated by the finite element technique, the velocity and acceleration contact constraints as well as the displacement contact constraint are imposed for the numerical stability without spurious oscillations. The solution of the present technique is checked by the numerical simulation using the concentric high-speed rotating disks with the clearance and impulsive loading. It is shown that the almost steady state solution agrees with the corresponding analytical solution of the elasticity and that the differentiated constraints are crucial for the numerical stability of such high-speed contact problems of the disks under impulsive loading.

1. 서 론

고속 회전 기계의 유한요소 모델에 동적인 접촉 조건을 부과여하여 안정된 해를 구하는 것은 동역학 및 진동 해석에서 중요한 문제의 하나이다. 특히 과도기 상태의 동역학 해석을 위하여는 운동방정식을 직접 시간적분하는 방법을 사용하여야 한다. 그런데 유한요소 모델의 동역학적 해석을 위하여 전통적인 Newmark 방법 등의 직접 시간 적분법을 활용하는 경우 수치해석 오차 누적으로 인하여 해의 안정성이 상실되므로, 일반적으로 에너지 또는 운동량 보존 법칙을 활용하여 시간적분법의 해를 보정하는 경우가 많다(1~3). 그러나 이러한 보존 법칙은 정확한 해의 필요조건일뿐 충분조건은 아니며, 특히 고속으로 회전 운동하는 접촉면에 적용되는 경우 해의 안정성을 보장할 수 없다. 최근 고속 회전 디스크에 접촉 조건을 부과한 동역학적 해석을 위하여 직접 시간적분법을 활용한 연구가 진행되고 있지만(4,5), 구속조건이 부과된 운동방정식의 안정된 시간적분에 대한 좀 더 체계적인 연구가 추가로 필요하다.

한편 고속회전 기계에 접촉 및 기구학적 구속조건을 부과하여 운동방정식을 적분하는 문제는 다물체 동역학에서 광범위하게 연구되었지만(6~8), 이 경우 접촉력을 국부적인 탄성변형 또는 윤활 현상을 이용하여 계산하므로 유한요소법을 이용한 해석에서는 활용할 수 없다 (즉 이 논문에서 접촉력은 구속조건의 Lagrange 승수 역할을 하므로, 침입거리 δ을 활용하여 접촉력을 kδn 등의 형식으로 표현하는 방법 등은 적용할 수 없음). 그런데 다물체 동역학에서 기구학적 구속조건의 경우에는 Lagrange 승수를 사용하여 구속조건을 표시하며 이를 미분한 속도 및 가속도 구속 조건을 운동방정식과 연립시켜 안정된 시간 적분을 가능하게 하고 있으므로, 이러한 기법은 유한요소법에서도 활용 가능할 것이다. 그러나 고속의 탄성체 사이에 접촉조건이 부과된 경우 접촉점이 변형된 접촉면 위를 고속으로 이동하므로 다물체 동역학 기법을 그대로 적용할 수는 없으며 접촉면의 탄성변형을 구속조건에 포함하는 해에 대한 추가 연구가 필요하다. 이러한 접촉문제와 관련하여 최근 저자는 접촉점에 변위뿐 아니라 속도 및 가속도 접촉조건을 부과하는 방법을 고안하여 안정된 시간 적분이 가능함을 입증하였다(9~11). 이 논문에서는 이러한 방법을 적용하여 초기 간격 및 충격 하중이 존재하는 고속회전 디스크 동접촉 문제를 풀기로 한다. 따라서 논문에서는 참고문헌(11)의 방정식 및 해법 중 필수적인 사항을 간단히 설명하고, 변위 접촉조건뿐 아니라 속도 및 가속도 접촉조건을 동시에 적용함으로써, 충격적인 하중이 가해지는 고속 디스크 접촉 문제의 경우에도 운동방정식을 시간적분하여 안정된 해를 구할 수 있음을 예제 계산을 통하여 설명한다.

 

2. 운동방정식 및 접촉조건

이 논문에서는 Fig. 1과 같이 고속 회전하는 디스크 A와 디스크 B 사이에 탄성 변형에 의하여 접촉이 발생하는 경우의 문제를 참고문헌(11)과 같이 2차원 유한요소 모델을 사용하여 해석하며, 특히 두 디스크 사이의 간격으로 인하여 충격적인 접촉력이 발생하는 경우를 고려한다. 디스크 A와 B는 평면 응력 요소를 사용하여 모델링하고, 디스크 A의 절점에서 디스크 B 표면에 수직선을 그어서 접촉점을 결정한다. 그런데 두 디스크는 강체 회전 운동은 크지만 내부의 탄성변형은 미소하므로 참고문헌(12) 등에서 설명된 바와 같은 물체에 부착되어 움직이는 좌표계를 활용하여 미소 탄성 변형을 계산하기로 한다. 이와 같은 목적으로 여기에서는 Fig. 1의 디스크 A의 경우 2점 I, J에 의하여 디스크 A에 부착된 x’ 축을 결정하고 여기에 수직 방향으로 y’ 축을 결정하며, 이러한 x’, y’축에서 미소 탄성변형을 계산한다 (물론 디스크 B에도 동일한 방법을 적용한다). 이 경우 참고문헌(11)에서 기술된 바와 같이 공간에 고정된 직각 좌표계 (x,y)에서 디스크 A의 유한요소법 운동방정식은 다음과 같이 표시된다.

Fig. 1Model of the two concentric disks having gap between them

위에서 M과 K는 일반적인 미소 변형 유한요소법에서 계산되는 간단한 질량 및 강성도 행렬이며, R은 물체에 부착된 좌표계 (x’,y’)와 공간에 고정된 (x,y) 좌표계 사이의 변환행렬, u는 대규모 강체운동과 미소 탄성변형을 합한 변위, ū는 물체에 부착된 좌표계의 대규모 강체운동으로 인한 변위, f는 외력, p는 (x,y) 좌표계에서의 접촉력을 나타낸다. 따라서 물체에 부착된 좌표계에서 측정한 미소 탄성 변위는 uA − 로 계산된다. 여기에서 물체에 부착된 좌표계 (x’,y’)의 강체회전 각도가 정확히 계산된다면 RA 와 u 는 정확히 계산할 수 있다. 이 논문에서는 다음 3장의 Steps 1.1~1.3에 의하여 물체에 부착된 좌표계의 정확한 회전각도를 2점 I, J를 잇는 선의 회전에 의하여 변위 u와 연립시켜 정확히 계산한다.

일반적인 유한요소법의 node-to-segment 접촉 기법에 의하여 디스크 A의 절점에 작용하는 접촉력이 디스크 B의 접촉점에 집중하중으로 작용하면, 디스크 B의 운동방정식은 다음과 같이 유도된다.

위에서 N은 접촉점에서의 형상 함수이다.

위의 방정식 (1)과 (2)의 해에 의하여 접촉점 i에서 두 물체 사이의 간격 (격리 거리)을 계산할 수 있는게 이 간격을 si라고 하며, 접촉점 i에서 접촉면에 수직인 접촉력을 pi라고 하면, 시간 t+Δt에서 접촉조건은 다음과 같이 된다.

 

3. 전체 방정식 해법

이 논문의 목적은 접촉조건 (3)을 만족시키는 방정식 (1)과 (2)의 해를 Newmark 수치적분 방법에 의하여 구하는 것이다. 그런데 1장에서 언급된 바와 같이 식 (1)과 (2)는 미분방정식이며 식 (3)은 대수부등식이므로 Newmark 방법 등 일반적인 시간적분법을 적용할 경우 그 해는 불안정하게 된다. 이러한 수치 불안정성을 해소하기 위하여, Lagrange 승수법을 사용하는 다물체 동역학에서와 같이, 이 논문에서는 접촉조건 (3)을 시간으로 미분하여 이를 운동방정식과 연립시켜 시간적분하는 방법을 사용한다. 특히 접촉조건이 등식 및 부등식으로 구성된 구속조건이며 또한 접촉점은 탄성변형이 발생하는 접촉면 위를 고속으로 이동하므로, 참고문헌(9~11) 기법을 활용하여 부록에 설명된 속도 및 가속도 접촉 오차 벡터를 0으로 만드는 방법을 사용하기로 한다.

이 논문에서는 다음에서 설명된 바와 같이 매 시간 단위마다 전체 운동방정식을 변위, 속도, 가속도 접촉조건을 각각 부과하여 수치적분 해를 구한다. 따라서 매 시간 단위에서 운동 방정식 해의 변위, 속도, 가속도는 각각 해당 접촉 조건을 부과시킨 경우의 해에 의하여 결정한다. 표준 Newmark 시간적분법 (𝛽1=0.25, 𝛽2=0.5)을 사용하며, 식 (1)과 (2)의 질량 및 강성도 행렬은 계산 초기에 한 번만 계산하면 된다. 이 경우 식 (1)과 (2)의 시간적분 해에서 대부분의 컴퓨터 계산 시간이 가우스 소거법에서 소요되지만, 시간 간격 Δt가 일정한 경우, 거대 행렬 K+1/(𝛽1Δt2)M의 forward reduction은 계산 초기에 1번만 시행하고, 그 후의 모든 계산에서는 그 행렬의 back substitution만 시행하면 되므로 대단히 경제적인 계산이 가능하다. 전체 계산 과정은 참고문헌(11)과 유사하며, 시간 t+Δt에서의 과정을 약술하면 다음과 같다.

Step 1.1 : 물체에 부착된 좌표계의 회전 각도 (θ A)t+Δt,k 에 의하여 (처음 k=1인 경우, (θ A)t+Δt,1은 (θ A)t 및 시간 t에서의 속도에 의하여 예측), 식 (1)의 좌표 변환 행렬 RA 및 강체 운동 변위 를 계산함. 디스크 B의 식 (2)에 대하여도 마찬가지로 계산함.

Step 1.2 : 부록의 반복 계산법 (A.3)에 의하여 접촉력 pt+Δt를 계산하고, 이것을 운동방정식(1)과 (2)에 대입하여 식 (1)과 (2)를 시간적분함. 계산 과정에서 식 (A.1)에 정의된 변위 접촉오차를 사용하며, 변위 접촉오차 크기 가 아주 작은 허용오차 이내로 될 때까지 반복 계산법에 의하여 정확한 접촉력 계산함 (이의 자세한 과정은 참고문헌(11) 참조). 변환 행렬 R은 직교 행렬이므로 (즉 R−1=RT), 가우스 소거법 계산과정에서 행렬 K+1/(𝛽1Δt2)M의 back substitution만이 필요함. 운동방정식 (1)과 (2)의 변위 ut+Δt 는 여기에서 결정됨.

Step 1.3: 위의 Step 1.2에서 계산된 운동방정식 (1)과 (2)의 해에 의하여 Fig. 1의 두 점 I와 J의 회전 각도를 계산하여 이 각도를 (θA)measured 라고 함. 이를 Step 1.1에서 사용된 회전 각도 (θA)t+Δt,k와 비교하여 만약 |(θA)measured − (θA)t+Δt,k|가 아주 작은 허용 오차 θtol 보다 작고, 또한 디스크 B에 대하여도 마찬가지 방법 적용하여 디스크 B의 회전각도가 결정된 경우, 다음의 Step 2 계산 진행함. 그렇지 않은 경우 디스크 A의 (k+1)번째 회전 각도를 다음 식에 의하여 보정한 후 다시 Step 1.1부터 반복 계산 시행함 (동일한 방법을 디스크 B에 대하여도 적용).

위에서 k는 강체 회전 각도를 계산하기 위한 반복계산법의 횟수를 나타내고, 𝛾는 상수인데 대부분의 계산에서 𝛾=1로 하면 θtol이 10−10 radin인 경우 1~2번의 반복계산에 의하여 θA와 θB가 결정됨.

Step 2: 위의 Step 1.2와 동일한 방법에 의하여, 속도 접촉 조건 (A.5)를 부과하기 위하여, 반복계산법 (A.6)에 의하여 접촉력을 계산하고 운동방정식 (1)과 (2)을 시간적분한다. 계산 과정에서 식 (A.4)에 정의된 속도 접촉오차를 사용하며, 속도 접촉오차 크기 가 아주 작은 허용오차 이내로 될 때까지 반복 계산법에 의하여 정확한 접촉력 계산함. 운동방정식 (1)과 (2)의 속도 는 여기에서 결정됨.

Step 3: 위의 Step 1.2와 동일한 방법에 의하여, 가속도 접촉 조건 (A.9)를 부과하기 위하여, 반복계산법 (A.10)에 의하여 접촉력을 계산하고 운동방정식 (1)과 (2)를 시간적분한다. 계산 과정에서 식(A.8)에 정의된 가속도 접촉오차를 사용하며, 가속도 접촉오차 크기 가 아주 작은 허용오차 이내로 될 때까지 반복 계산법에 의하여 정확한 접촉력 계산함. 운동방정식(1)과 (2_)의 가속도 는 여기에서 결정됨.

위의 반복계산법 (A.3), (A.6), (A.10)에 의하여 접촉력을 계산한 경우, 부록의 3개의 접촉오차 , , 는 모두 0을 향하여 단조 감소한다. 예로써 변위 접촉 오차 는 참고문헌 (11)에서 자세히 기술된 바와 같이 반복계산법 (A.3)에 의하여 단조 감소하며, 이를 다음에 간단히 약술한다. 디스크 A와 B 사이의 접촉면에서 각각의 유연성 행렬을 CA와 CB라고 하면, 수직 간격 s와 접촉력 p의 관계식은 다음과 같다 (여기에서 p는 접촉면에 수직방향인 성분으로 구성됨).

위의 식 및 변위 접촉 오차 (A.3)에 의하여, 참고문헌(11)에서 설명된 바와 같이, 반복 계산법 (A.3)의 α 가 작은 경우, 다음의 부등식이 성립한다.

위에서 m은 반복 횟수를 나타내고, 하첨자 non-zero 는벡터 의 0에 해당하는 성분을 모두 제거시킨 후 얻어지는 행과 열을 의미한다. 이제 λmax 와 λmin가 각각 행렬 (CA + NCBNT)non-zerot+Δt 의 최대 및 최소 고유치를 나타내기로 한다. 그러면 행렬 CA와 CB는 대칭이고 양정이므로, 반복계산법 (A.3)의 α가 0보다 크고 2 / λmax 보다 작은 경우, ∥et+Δt,m∥2 은 반복계산법 (A.3)에 의하여 0으로 수렴한다. 반복 계산법 (A.3)에서 최적의 α 는 2 / (λmax +λmin)이며 이를 위한 계산 과정은 참고문헌(11)에 설명되어 있지만, 비교적 간단한 모델의 경우 실용적으로 α 를 정도로 간단히 계산하여도 된다. 만약 위의 α 및 식 (A.3)에 의하여 가 얻어지지 않는 경우에는 (예로써 새로운 접촉점에서 접촉을 시작하는 경우 등) 참고문헌(11)에서 설명된 바와 같이 α 를 줄여가며 추가로 계산하여 를 얻게 된다.

 

4. 예제 계산

Fig. 1에서 보인 바와 같은 2개의 동심 디스크가 고속으로 회전하고 있으며, 디스크 A의 회전속도 ωA는 반시계 방향으로 2000 rpm, 디스크 B의 회전속도 ωB는 시계 방향으로 1000 rpm이다. 탄성계수 E는 70 GPa, 푸아송 비는 v는 0.3, 밀도는 2,600 kg/m3이며, Fig. 1에서 디스크 반경은 r0=15 cm, r1=15 cm, r2=r1+gap, r3=20 cm이다. 여기에서 gap은 변형되지 않은 상태에서 두 디스크 사이의 간격을 의미하며 이 예제에서는 gap=0.001 mm로 하여 계산을 수행한다. 또한 디스크 A의 내부압력이 초기 0에서 시간tloading에는 10 MPa까지 선형으로 증가하도록하며, 그 이후에는 10 MPa의 압력이 일정하게 유지되도록 한다. 유한요소법 계산을 위하여 각각의 디스크는 Fig. 2에서 보인 바와 같이 720개의 평면 응력 요소로 분할한다. 디스크 고속 회전으로 인한 초기 변형은 원심력을 외부의 체적력으로 환산하여 정역학적인 유한요소법 해에 의하여 구한다. 한편 부록에서 제시된 반복계산법 (A.3), (A.6), (A.10) 사용시 변위 접촉오차 허용치는 10−12 m, 속도 접촉오차 허용치는 10−12/Δt m/s, 가속도 접촉오차 허용치는 10−12/Δt2 m/s2으로 하여 계산한다(오차 한계를 위의 값보다 102 또는 10-2배로 하여도 실질적으로 동일한 결과가 얻어지며, 컴퓨터 수치 계산에서 속도 및 가속도 정밀도를 변위 정밀도에 비하여 1/Δt 및 1/Δt2배로 유지하여 계산 안정성을 확보함).

Fig. 2Finite element mesh of the two concentric rotating disks

만약 디스크 내부 압력이 증가하는 시간 tloading이 비교적 길다면 그 해는 정역학적인 해와 유사한 결과가 얻어질 것이다. 그런데 회전하는 디스크의 정역학적인 탄성 해는 축대칭 해석에 의하여 정확한 해석적인 해가 가능하다. 즉 정역학의 경우 회전하는 탄성체 디스크의 축대칭 해석에 의하여, 반경 방향 변위 ur 및 반경 방향 응력 𝜎r은 다음과 같다(13).

위에서 r은 반경, 𝞺는 밀도, C0, C1, D0, D1은 경계조건에 의하여 결정되는 상수이다. 이 예제와 같이 디스크 A와 B가 접촉하는 경우 경계조건은 다음과 같다.

위의 4개의 경계조건에 의하여 상수 C0, C1, D0, D1이 계산되며, 식 (7)에 의하여 반경 r1에서의 응력이 계산된다. 이 장에서는 식 (7)과 (8)에 의하여 해석적인 방법으로 계산된 반경 r1에서의 응력을 이 논문의 시간적분 수치해에 의한 접촉 응력과 비교하고자 한다.

먼저 디스크 내부 압력을 0.05 sec 동안 0에서 10 MPa로 증가 시키는 경우를 고려한다(즉 tloading = 0.05 sec). Newmark 시간 적분의 Δt는 0.01 msce로 하여, 이 논문에서 설명된 방법에 의하여 변위, 속도, 가속도 접촉조건이 부과된 운동방정식 시간적분 해를 구하였다. 이러한 수치 해에서 디스크 A 위의 절점 J에서 접촉 응력의 시간에 따른 변화가 Fig. 3에 그려져 있으며, 시간 0.07 sec에서 디스크 A 주변의 접촉 응력 분포가 Fig. 4에 그려져 있다. 이 경우 Fig. 3에서 보인 바와 같이 시간 0.05 sec 이후에 는 접촉력 변화는 거의 없으며 따라서 정상상태가 유지되고 있다. 정상상태에서 방정식 (7)의 해석적 해에 의한 접촉 응력이 2.572 MPa로 계산되므로 Fig. 3의 접촉 응력은 정확히 계산되었다고 보인다. 또 이 유한요소법에서는 (x,y)직각 좌표계에서 2차원 평면 응력 요소를 사용하여 계산을 수행하였지만 Fig. 4에서 보인 바와 같이 접촉 응력은 축대칭 형태로 분포하며, 역시 이 계산의 정확도를 뒷받침하고 있다. 특히 Fig. 4에서 보인 바와 같이, 변위, 속도, 가속도 접촉 조건에 의하여 계산되는 접촉 응력은 모두 동일하며, 이 계산 방법의 타당성을 입증한다(Figs. 4~9에서 변위, 속도, 가속도 해는 모두 색깔을 이용하여 표시했으며, 학회의 PDF 파일 참조 요망).

Fig. 3Contact pressure history on node J computed with displacement, velocity, and acceleration contact constraints(tloading : 0.05 sec)

Fig. 4Contact pressure distribution around the contact surface of disk A at time 0.07 sec computed with displacement, velocity, and acceleration contact constraints(tloading : 0.05 sec)

Fig. 5Contact pressure history on node J computed with displacement, velocity, and acceleration contact constraints(tloading : 0.005 sec)

Fig. 6Detail of contact pressure history on node J computed with displacement, velocity, and acceleration contact constraints(tloading : 0.005 sec)

Fig. 7Contact pressure history on node J computed with only displacement contact constraint without velocity and acceleration contact constraints( tloading : 0.005 sec)

Fig. 8Contact pressure history on node J computed with displacement and velocity contact constraints without acceleration contact constraint (tloading : 0.005 sec)

Fig. 9Contact pressure history on node J computed without Coriolis and centripetal accelerations in the acceleration contact constraint(tloading : 0.005 sec)

다음으로 디스크 내부 압력을 0.005 sec 동안 0에서 10MPa로 증가 시키는 경우를 고려한다(즉 tloading =0.005 sec). 이 경우에는 Fig. 5에서 보인 바와 같이 짧은 시간에 가해지는 충격적인 하중으로 인하여 진동이 발생한다. Fig. 5의 해는 앞 3장에서 설명된 바와 같이 변위, 속도, 가속도 접촉조건을 모두 가하여 얻어진 것이며, 이렇게 얻어진 해는 Fig. 6에서 자세히 보인 것처럼 접촉이 새롭게 시작되는 순간을 제외하고는 3개의 해가 일치한다. 그러나 변위 접촉조건만 부과하는 경우에는 (즉 속도 및 가속도 접촉조건을 부과하지 않으면) 그 시간적분 해는 Fig. 7처럼 발산한다. 또한 변위 및 속도 접촉조건을 부과하여 (즉 가속도 접촉조건을 부과하지 않으면) 얻어진 시간적분 해는 Fig. 8에서 보인 것처럼 상당히 안정된 해가 얻어지지만, Fig. 6과 비교하여 초기에 불안정 현상이 포함되어 있음을 알 수 있다. 또한 가속도 접촉조건을 부과하는 경우에도, 식 (A.7)의 Coriolis 가속도 및 원심 가속도를 포함하지 않으면, Fig. 9에서 보인 바와 같이 변위, 속도, 가속도 접촉조건을 부과한 해가 서로 일치하지 않는 부정확한 해가 얻어진다. 따라서 식 (A.7)의 Coriolis 가속도 및 원심 가속도는 가속도 접촉조건에서 필수적임을 알 수 있다.

 

5. 결 론

고속회전 디스크 사이에 접촉조건이 부과되는 경우에, 그 사이 간격으로 인하여 충격적인 접촉력이 발생하는 경우에도, 유한요소 운동방정식의 시간적분 수치 해에 변위 접촉조건뿐 아니라 속도 및 가속도 접촉조건을 부과하면 안정된 수치적분 해가 얻어질 수 있음을 설명하였다. 또한 비교적 정상해에 가까운 경우의 예제 계산에서 탄성학의 해석적인 해와 비교하여 이 계산 방법의 정밀도를 점검하였다. 이 논문에서 고속 회전 디스크의 접촉 문제와 관련된 수치적분 해의 정확도 및 안정성을 점검하였으므로, 추후 참고문헌(11) 등의 마찰 접촉 관련 식을 도입하여 고속 회전체 사이에 마찰이 존재하는 경우의 해도 마찬가지로 해석 가능하다고 예상된다. 예를 들어 자동차 및 기차 브레이크 등 각종 고속 회전체 사이에서 마찰에 의하여 유도되는 진동 해석에도 응용 가능 하리라고 예상되며 추후에 계속 연구가 필요하다.