I. 서  론

원통형 실린더는 파이프라인, 해양플랜트, 풍력 타워 등 여러 목적의 구조물에 기본적으로 쓰이는 형상이다. 이러한 구조물에서의 진동 및 파동 전파 특성에 관해서는 많은 연구가 이루어졌다. Pochhammer^[1]와 Chree^[2]는 탄성 이론을 이용하여 꽉 찬 무한 원통형 실린더에서의 파동 전파 특성 연구를 수행하였고 McFadden,^[3] Ghosh^[4]는 쉘 형태의 원통형 실린더에서의 축대칭 모드에 대해 연구하였다. 축대칭 모드를 포함한 일반적인 경우에 관한 원통형 실린더의 3차원 운동은 Gazis^[5]가 탄성 이론을 이용하여 행렬 형태의 방정식으로 기술하였다. 이러한 이론들은 그 적용 목적에 따라 때때로 원통형 실린더의 안쪽 면이 유체와 접해 있는 경우, 바깥쪽 면이 유체에 접수되어 있는 경우, 혹은 안과 바깥쪽 양면이 동시에 같거나 혹은 서로 다른 유체에 접해 있는 경우까지 확장될 수 있다. 그러나 이런 탄성 이론들은 다루어야 할 행렬 요소의 개수가 많아 연산을 수행하는데 있어 때때로 많은 메모리를 필요로 하며 또한 일반해로 나타나게 되는 베셀함수를 모든 주파수와 파수영역에서 다루는데 있어 이를 정확한 값으로 수치화 하는데 어려움이 따른다.

한편 원통형 실린더의 자유진동에 관한 보다 간략화 된 이론들에 대해 많은 연구가 되어져 왔는데 이는 많은 경우에 있어 실린더의 두께가 반경에 비해 얇다는 사실로부터 도출된 이론들이다. 이 이론들은 기본적으로 Donnel-Mushtari 두께에 의한 영향을 해석함에 따라 그 방법 및 적용성에 따라 Love-Timoshenko, Goldenveizer-Novozhilov, Houghton-Johns의 모델이 존재하게 된다. 이러한 이론들에 대해서는 Leissa^[6]가 폭넓게 정리를 하였고 이에 따르면 얇은 원통형 실린더의 운동은 강성 행렬 형태의 방정식으로 기술된다. 이러한 얇은 원통형 쉘 기반 이론들은 앞선 탄성을 이용한 이론에 비해 간단한 행렬 요소를 이용해 빠른 연산 결과를 보여주는 장점이 있다. 본 논문에서는 얇은 원통형 쉘 이론 방식이 탄성이론 방식보다 10배 이상의 빠른 연산 속도를 보여주었다.

최근에는 컴퓨터의 연산속도의 획기적인 증대로 인해 유한요소법과 같은 수치해석 방법도 비교적 수월하게 수행할 수 있게 되었다. 이에 따라 원통형 실린더 뿐만 아니라 임의의 형상에 대한 실린더의 진동해석에 대한 수치해석 결과 정확도 및 신뢰성이 향상되었다. 본 연구에서는 탄성 이론과 얇은 원통형 쉘 이론을 적용한 결과를 검증하기 위해 상용 FEM 도구인 COMSOL 4.3b이 사용되었다. 물론 COMSOL을 이용한 결과가 탄성이론이나 얇은 원통형 쉘 이론 방식보다는 계산 시간이 많이 소요되지만 본 연구에서는 1분 이내의 비교적 빠른 연산 속도를 보여주었다.

본 논문에서는 무한 원통형 실린더에 대해 탄성 이론, 얇은 원통형 쉘 이론, 그리고 상용 소프트웨어인 COMSOL을 이용한 FEM을 통한 수치해석 방법을 적용하여 각각 자유 진동에서의 고유진동수를 구하였다. 고유진동수에서 형성되는 실린더의 고유 모드는 단면에 국한되어 진동하는 모드이다. 따라서 탄성 이론 및 얇은 원통형 쉘 이론에서는 축방향의 파수가 사라지는() 영역에서 고유주파수를 구하였고, 유한요소법에서는 2차원 공간에서 실린더의 단면을 구축하고 실린더 두께에 메쉬를 형성하여 고유주파수를 구하였다. 또한 주어진 반경의 실린더에 대해 그 두께를 각기 달리하여 각 방법을 통해 얻은 1 kHz 이하에 존재하는 원통형 실린더의 고유주파수 값이 어떻게 차이가 나는 지를 분석하였다.

II. 원통좌표계에서의 파동방정식

2.1 탄성 이론 접근법

외력이 작용하지 않는 탄성체에서의 진동운동은 변위 벡터 를 이용하여 다음과 같은 운동방정식으로 기술이 된다.^[7]

 (1)

여기서 는 탄성체의 밀도를 의미하며 는 각각 제 1, 제 2 Lam 상수이다. Eq.(1)은 어떠한 게이지 변환에서도 성립하며 Eq.(1)로부터 도출된 결과는 어느 변환 좌표계에서도 동일한 운동 특성을 기술하여야 한다. 이러한 특성을 게이지 불변성이라 부르고 이에 따라 임의의 좌표계에서의 변위 벡터 는 항상 스칼라 포텐셜 와 벡터 포텐셜 에 의해 Eq.(2)와 같은 관계로 표현 할 수 있다.

. (2)

Eq.(1)은 의 성질에 의해 다음과 같은 각각의 포텐셜에 대한 독립된 파동방정식을 얻게 된다.

, (3)

, (4)

여기서 와 는 각각 탄성체에서의 종파와 횡파의 속도를 의미하며 Lam 상수와 , 의 관계를 가진다. 위 파동 방정식이 시간에 대해 의 형태를 따른다고 가정하면 Eqs.(3)과 (4)는 각각 아래와 같은 헬름홀츠 방정식으로 표현된다.

, (5)

, (6)

여기서 , 이다. Fig. 1과 같은 원통형 좌표계에서 축방향( 방향)으로 전파하는 모드들이 공간에 대해 의 형태를 따른다고 가정하면 Eqs.(5)와 (6)에 대해 베셀함수의 무한합으로 표현이 되는 해를 구할 수 있게 된다. Eqs.(5)와 (6)에 대한 해의 표현식을 Eq.(2)에 대입함으로써 원통형 실린더에 대한 구체적인 변위 벡터를 구할 수 있고 구해진 변위 벡터를 아래 식에 대입함으로써 원통형 실린더의 응력성분을 구할 수 있게 된다.

Fig. 1. Configuration of a circular shell for elasticity theory.

, (7)

, (8)

. (9)

Fig. 1은 탄성 이론을 적용하기 위해 임의의 두께를 가진 원통형 실린더를 도식화한 것이다. 실린더의 내경 및 외경은 각각 , 로 표시를 한다. 이에 따라 원통형 실린더의 임의의 두께 는 로 주어진다.

Eqs.(2), (7) ~ (9)로부터 원통형 실린더의 변위 및 응력성분에 대해서는 6 개의 정해지지 않은 계수가 포함된 식을 얻게 되며 이러한 6 개의 정해지지 않은 계수들 사이의 관계는 원통형 실린더의 경계면인 와 에서의 경계조건으로부터 얻게 된다. 실린더의 응력 성분들에 각각 와 에서 과 같은 경계 조건을 적용하면 Eq.(10)과 같은 고유모드 특성 방정식을 얻는다.

, (10)

여기서 는 임의의 계수이며 행렬 요소 의 구체적인 표현식은 Reference [5]에 자세히 실려 있다. 각 행렬요소 는 각 모드 에 대한 함수이며 그 표현식은 제 1종, 2종 베셀 함수, 혹은 제 1종 2종 한켈 함수, 혹은 경우에 따라 변형(modified) 베셀 함수들의 선형결합으로 표시가 된다. 첫 번째 및 네 번째 열의 행렬 요소는(에서 ) 탄성 원통형 쉘의 종파의 운동과 관련된 항이며 나머지 네 열의 행렬 요소는(에서 ) 탄성 원통형 쉘의 횡파의 운동과 관련된 항이다. 또한 처음 세 행 및 나머지 세 행의 행렬 요소는(에서 과 ) 동일한 형태로 표현되나 각각 와 에서의 기술되는 차이만 있을 뿐이다. 따라서 이러한 차이가 운동방정식의 두께 방향에 대한 의존도를 생성하게 되는데 Eqs.(5)와 (6)의 해로 구해지는 베셀함수의 와 에서의 값들이 서로 선형적으로 연관되어 있지 않고 주기 함수 형태로 연관되어 있으므로 와 의 차이가 클수록 그 의존도가 커진다. 각 모드에 따른 탄성 원통형 쉘의 고유 진동수를 구하는 방법은 여러 가지가 있으나 본 논문에서는 각 모드 에 대한 고유주파수 은 이 되는 조건에서 을 만족시키는 주파수를 주파수에 따른 그래프를 통하여 구하였다. 여기서 는 Eq.(10)에서의 행렬의 determinant를 의미한다.^[5]

2.2 얇은 원통형 쉘 이론 접근법

Fig. 2는 얇은 원통형 쉘 이론을 기술할 때 적용되는 좌표계의 도식도이다. 앞선 탄성 이론 접근법과는 달리 이 이론에서는 원통형 쉘의 두께가 반경에 비해 작다고 가정하므로 실린더의 외경, 내경을 사용되지 않고 실린더의 평균반경 을 이용한다. 이러한 경우의 운동 방정식은 Eq.(11)과 같은 형태로 기술된다.

, (11)

여기서 은 강성 행렬, 는 변위 행렬, 는 외력 행렬이다. 여기서 강성 행렬 은 유도 방법에 따라 Love-Timoshenko, Goldenveizer-Novozhilov, Houghton-Johns 등의 여러 가지 모델이 존재하게 된다. 그러나 어떠한 모델에 대해서도 은 Eq.(12)와 같이 표현할 수 있다.^[6]

Fig. 2. Configuration of a circular thin-shell.

, (12)

여기서 은 Donnel-Mushtari 모델에 의한 강성 행렬, 는 쉘의 두께와 관계되는 무차원 패러미터이며 는 각 모델에 의해 확장되는 강성 행렬이다. 본 논문에서는 수치 해석적으로 우수하다고 알려진 Goldenveizer-Novozhilov 모델을 사용하여 기술하였다.^[6,8]

Fig. 2과 같은 좌표계에서 변위 행렬 의 성분은 축방향, 접선방향, 반경방향의 변위인 로 표시된다. 이때 각각의 변위 성분들이 의 공통된 항을 가진다고 가정한다면 Golden-veizer-Novozhilov 모델에 의한 강성 행렬 의 성분은 공간에서 Reference [8]의 식(9.1.6)에 자세하게 기술이 되어 있다.

 외력이 작용하지 않을 경우의 조건은 Eq.(11)에서 이 되며 얇은 원통형 쉘 이론을 이용한 각 모드 에 대한 고유주파수 은 이 되는 조건에서 을 만족시키는 주파수를 찾음으로써 얻을 수 있다. 여기서 는 Eq.(11)에서의 강성행렬의 determinant를 의미한다.

2.3 수치해석 접근법

본 논문에서는 유한요소법(FEM)을 통해 원통형 실린더의 고유주파수 을 구하였다. 이를 위해 상용 프로그램인 COMSOL Multiphysics 4.3b의 Structural Mechanics > Eigenfrequency analysis를 사용하였다. COMSOL에서 제공하는 외력이 작용하지 않는 경우 본 프로그램에서 제공되는 탄성체의 운동 방정식은 Eq.(15)와 같다.^[9]

, (13)

여기서 는 탄성체의 변위를 의미하는 텐서성분이며 는 탄성체에 작용하는 응력을 나타내는 텐서 성분이다. 고유진동수는 Eq.(13)을 만족시키는 주파수이며 이를 위해 2D 평면에서 원통형 실린더를 구성하여 계산하였다.

수치해석에 사용된 원통형 실린더의 물성치 및 크기는 Table 1과 같으며 1 kHz 이하에 존재하는 고유진동수를 수치해석을 통해 계산하였다. Fig. 3은 COMSOL을 이용한 수치해석 방법을 통해 구한 고유 모드의 형상의 예를 보여준다. 그림에서 사용된 원통형 실린더는 외경 1 m에 10 mm의 두께를 가지며 이와 같은 수치해석 방법으로 구해지는 고유주파수를 각각 탄성 이론적 방법 및 얇은 원통형 쉘 이론적 방법의 결과와 비교하였다.

Fig. 3. Mode shapes of the first four eigenmodes of the circular cylinder found in the frequency range below 1 kHz using COMSOL. Above are the calculations for the cylinder with the thickness of 10 T.

III. 무한 원통형 실린더의 고유주파수

무한 원통형 실린더의 고유주파수를 계산하기 위해 본 논문에서는 Steel AISI 4340 재질의 실린더를 사용하였다. 이에 대한 물성치는 Table 1과 같다. 각각의 계산 결과가 실린더의 두께와 어떤 관계가 있는 지를 파악하기 위해 4가지 두께(10 T, 20 T, 50 T, 100 T)를 선택하였고 외경은 1 m로 동일하다.

탄성 이론에 적용하기 위해서는 종파 속도 와 횡파속도 가 필요한데 이는 Table 1의 물성치로부터 Eqs.(14)와 (15)를 통해 구할 수 있다.

, (14)

. (15)

Table 2는 4가지 두께에 대해 무한 원통형 실린더의 고유주파수를 탄성 이론(A로 표시), 얇은 원통형 쉘 이론(B로 표시), COMSOL 수치해석(C로 표시)을 통해 얻은 결과를 정리한 결과이다. 고유주파수는 1 kHz 이하의 결과만 정리하였으며 각 모드수 n 별로 Hz의 단위로 소수점 첫째 자리까지 계산하였다.

계산 결과에 따르면 주어진 크기와 물성치의 원통형 실린더의 고유주파수가 1 kHz 이하인 경우에 축대칭모드(), 병진모드()에 대한 고유주파수는 나타나지 않으며 제일 낮은 고유주파수는 쌍극자모드()에 대한 고유주파수이다. 축대칭모드() 및 병진모드()에 대한 고유주파수는 1 kHz 이상에서 얻을 수 있다.

Fig. 4. Errors between eigenfrequencies obtained from three methods described in Section II. Open circles and triangles represent the deviation of the results by thin-shell theory and COMSOL numerical simulations, respectively, from the results by the elasticity theory.

Fig. 4는 Table 2에 정리된 각 고유 모드별 고유진동수에서 각각의 방법으로 계산한 결과 사이의 오차를 보여준다. 원형 마크는 수치해석 결과와 탄성 이론과의 오차(즉, C-A) 이며 삼각형 마크는 얇은 원통형 쉘 이론과 탄성 이론과의 오차(즉, B-A)를 보여준다.

얇은 원통형 쉘 이론에 의해 구한 원통형 실린더의 고유진동수는 탄성 이론 및 수치해석에 의해 구한 고유진동수와 다소 차이를 보이며 모드수 n 이 증가할 경우, 그리고 실린더의 반경 대비 두께의 비()가 증가할 경우 그 차이가 더 크게 일어남을 알 수 있다. 이는 두께가 증가할수록 얇은 원통형 쉘 이론이 탄성 이론 결과와 잘 일치하지 않음을 의미하는데 그 이유는 얇은 원통형 쉘 이론에서는 두께 사이에 존재하는 실린더의 변위(즉, 에서의 변위와 에서의 변위)는 늘 같다고 가정하기 때문이다.^[10] 반면 수치해석결과는 모드수 n 및 고유주파수가 커질 경우에도 탄성 이론과 거의 같은 결과를 얻을 수 있었다. 원통형 실린더의 두께 사이에 형성되는 주파수별 메쉬 사이즈를 최대 (는 파장) 이하로 설정하여 두께 사이에 존재하는 실린더의 에서의 변위와 에서의 변위차이를 충분히 보여줄 수 있기 때문이다.

Table 2에 정리된 각 모드별 고유진동수 분석을 통해 얇은 원통형 쉘 이론에서 구한 Steel AISI 4340쉘의 고유진동수와 탄성이론 결과 사이의 오차를 1 % 미만으로 얻기 위한 조건을 시뮬레이션을 통해 경험으로 아래와 같이 도출하였다.

, (16)

여기서 는 쉘의 평균 반경이며 은 모드수를 의미한다. Eq.(16)에 따르면 일 때 의 비가 각각 0.079, 0.051, 0.042, 0.038 보다 작을 경우 얇은 원통형 쉘 이론 결과와 탄성 이론 결과와의 오차가 1 % 이하인 것으로 예상된다. 실제 외경이 2 m 인 경우 50T의 두께를 가진 원통형 쉘()의 고유진동수를 에서 계산해 보면 얇은 원통형 쉘 이론 결과와 탄성이론은 각각 302.6 Hz와 299.5 Hz가 나오며 이때의 오차는 1.04 %가 된다. 따라서 1 % 미만의 오차는 의 조건에서 관찰된다. 이러한 방식으로 각각 0.5 m, 1 m, 2 m의 Steel AISI 4340 실린더의 고유 진동수를 계산해 보면 얇은 원통형 쉘 이론 결과와 탄성 이론 결과와의 오차가 Eq.(16)으로부터 계산된 Table 3에 정리된 조건일 경우에 1 % 이하인 것으로 확인이 되었다 .

IV. 결  론

무한 원통형 실린더에 대해 탄성이론, 얇은 원통형 쉘 이론, COMSOL을 이용한 수치해석 방법을 적용하여 1 kHz 이하에 존재하는 모든 고유주파수를 구하였다. 모든 두께 및 모드수에 대해 탄성 이론과 수치해석 결과가 0.4 Hz 미만의 오차를 보임으로써 이 두 해석 방법에 의한 결과가 거의 일치함을 알 수 있었다. 얇은 원통형 쉘 이론에 의해 구해진 결과는 실린더의 두께가 외경에 비해 작지 않거나 혹은 모드수 n이 증가할수록 탄성이론 결과와의 차이가 커짐을 확인 할 수 있었다. 또한 각 모드별 고유진동수 분석을 통해 얇은 원통형 쉘 이론에서 구한 Steel AISI 4340쉘의 고유진동수와 탄성이론 결과 사이의 오차를 1 % 미만으로 얻기 위한 경험적인 조건을 시뮬레이션을 통해 도출하였다. 향후 유체가 접해 있을 경우 이들 세 방법에 의한 고유 주파수의 변화를 관찰해 볼 예정이며 이를 통해 응용분야에 따른 효율적인 계산 방법을 도출해 낼 수 있을 것으로 기대된다.