The Korean Society for Power System Engineering

Journal Archive

Journal of Power System Engineering - Vol. 24 , No. 6

[ Research Paper ]
Journal of Power System Engineering - Vol. 24, No. 6, pp. 59-66
Abbreviation: The J. of the Korean Soc. for Power Syst. Eng.
ISSN: 2713-8429 (Print) 2713-8437 (Online)
Print publication date 31 Dec 2020
Received 31 Aug 2020 Revised 25 Nov 2020 Accepted 25 Nov 2020
DOI: https://doi.org/10.9726/kspse.2020.24.6.059

강성계수의 전달을 이용한 3차원 고체 구조물의 자유진동 해석
최명수*,
*교수, 전남대학교 기관시스템공학과

Free Vibration Analysis of Three Dimensional Solid Structure Using Transfer of Stiffness Coefficient
Myung-Soo Choi*,
*Professor, Department of Power System Engineering, Chonnam National University.
Correspondence to : Myung-Soo Choi : Professor, Department of Power System Engineering, Chonnam National University. E-mail : engine@jnu.ac.kr, Tel : 061-659-7137


Abstract

The author developed a computer program for the free vibration analyses of three dimensional solid structures by using the Sylvester-transfer stiffness coefficient method introducing the rectangular hexahedron element of the finite element analysis. In this paper, the computational algorithm for the free vibration analyses of three dimensional solid structures is formulated by the proposed method. The computer program is applied to the free vibration analyses of three dimensional solid structures which are a cube, an anvil, and a three dimensional cantilever beam structure. The trust of the proposed method is confirmed by comparing numerical results of the proposed method with those of the finite element method, experiment and analytical solution.


Keywords: Free vibration, Three dimensional solid structure, Sylvester-transfer stiffness coefficient method, Rectangular hexahedron element

1. 서 론

형상이 단순한 보나 판과 같은 구조물의 정적 및 동적 해석은 방정식을 세우고 이를 직접 풀어서 해를 구할 수 있지만 구조물의 형상이 복잡해지면 이와 같은 방법으로는 정도(accuracy)가 높은 해를 구하는 것은 현실적으로 곤란하다. 따라서 근래에는 컴퓨터 및 전산 프로그램을 이용하여 상기 문제점을 해결하기 위한 많은 연구가 진행되어 왔고, 그 결과 대표적인 전산 알고리즘으로 유한요소법과 전달행렬법 등이 있다.

전달행렬법은 축(shaft)과 같이 긴 구조물의 해석에 효과적인 알고리즘이라 알려져 있지만 구조물이 극단적으로 길어지면 해석 결과에 오류가 발생할 가능성이 높고, 자유진동 해석의 고유 진동수를 찾는 과정에서 일부 고유 진동수를 빠뜨릴 가능성이 존재하는 취약점이 있다.

반면에 유한요소법은 컴퓨터를 이용한 계산 과정에서 해석 모델의 전체 자유도에 상당하는 강성행렬 및 질량행렬을 일시에 이용해야 하므로 전달행렬법에 비해 계산 메모리의 활용 측면에서 부담이 많은 방법이다.

저자는 과거의 연구에서 강성계수의 전달을 이용하여 전달행렬법이 추구하는 목적을 유지하면서 전달행렬법의 취약점을 극복하기 위한 방안 연구를 수행한 바 있다.1),2) 그러나 이 방법은 적용 대상이 단면의 치수에 비해 길이가 긴 막대나 얇은 판과 같은 구조물에 제한되어 있어서 해석 방법으로서의 확장성이 부족하였다.

근래에는 복잡한 기계나 구조물이라 할지라도 3차원 캐드로 그린 후에 구조해석 소프트웨어의 자동 요소망(auto mesh)3) 생성 기능 및 3차원 고체 요소(solid element)4)를 활용하여 모델링할 수 있다. 따라서 이 연구에서는 강성계수의 전달을 이용하여 3차원 고체 구조물의 자유진동 해석을 수행할 수 있는 전산 프로그램을 개발한다. 그리고 정육면체 구조물, 앤빌(anvil), 3차원 외팔보 구조물을 대상으로 자유진동 해석을 수행하고 그 결과를 고찰한다.


2. 자유진동 해석 알고리즘

이 연구에서 제시하게 될 3차원 고체 구조물의 자유진동 해석을 위한 알고리즘을 명확히 설명하기 위하여 Fig. 1과 같이 탄성 지지된 직육면체를 해석 모델로 선정한다.


Fig. 1 
Analytical model

2.1 모델링

해석 모델을 Fig. 1과 같이 x축, y축, z축 방향과 수직하게 분할하면 해석 모델은 다수의 직육면체 요소(rectangular hexahedral element)5)로 분할된다. 해석 모델을 일차적으로 x축 방향과 수직하게 분할한 것을 블록(block)이라 정의한다(Fig. 2 참고). 해석 모델에서 블록을 나누는 분할 면과 해석 모델의 시작 면과 끝 면을 통칭하여 절면(nodal plane)이라 부른다. 따라서 해석 모델을 n개의 블록으로 분할하면 절면은 총 n+1개가 존재하며, 해석 모델의 시작 면에서 끝 면까지 각각의 절면을 순차적으로 절면 1, 절면 2, …, 절면 n+1이라고 정의한다.


Fig. 2 
Block I

Fig. 3에 나타낸 직육면체 요소는 8개의 꼭짓점에 절점을 하나씩 가지며, 각 절점은 x, y, z축 방향의 변위 u, v, w로 구성되는 3자유도를 갖는다. 따라서 직육면체 요소의 질량행렬(M(e))과 강성행렬(K(e))은 크기는 24×24가 된다.


Fig. 3 
Rectangular hexahedron element

블록의 질량행렬과 강성행렬은 블록을 구성하는 직육면체 요소의 질량행렬과 강성행렬을 조립하여 구할 수 있다. 블록의 질량행렬과 강성행렬을 각각 M(b)M(b)라고 하면, 블록의 동강성행렬 D(e)의 부분행렬(submatrix) A, B, C는 다음 식으로부터 유도할 수 있다.

Db=A BBT C=Kb-ω2Mb(1) 

여기서 ω는 고유 각진동수이다.

2.2 강성계수행렬의 정의 및 전달

임의의 절면 i의 변위벡터를 di라고 하고, 절면 i의 좌측과 우측의 힘벡터를 각각 f^ifˇi라 하면, 이들 변위벡터와 힘벡터의 관계를 다음 식과 같이 정의한다.

f^i=S^idi(2) 
fˇi=Sˇidi(3) 

여기서 S^iSˇi는 각각 절면 i 좌측 및 우측의 강성계수행렬이다.

절면 im개의 절점이 있고 이들 절점이 기초지지 스프링으로 지지될 경우, i번째 절면에서 힘의 평형식은

fˇi=f^i+Pidi(4) 

여기서

Pi=diagkx1,ky1,kz1,kx2,ky2,kz2,,kzmi(5) 

그리고 kxj, kyj, kzj는 절점 j를 기초로부터 지지하는 스프링의 x축, y축, z축 방향의 스프링상수이다.

식 (4)식 (2)(3)을 대입하여 정리하면 식 (6)과 같이 i번째 절면 좌우측의 강성계수의 전달식을 구할 수 있다.

Sˇi=S^i+Pi(6) 

i번째 블록 좌우측의 힘벡터와 변위벡터 사이의 관계식은 식 (1)로부터 다음 식과 같이 유도된다.

fˇif^i+1=-Ai -BiBiT Cididi+1(7) 

식 (3)(7)로부터 i번째 블록 좌우측의 변위벡터 사이의 관계식을 유도하면

di=Vidi+1(8) 

여기서

Vi=-Gi-1Bi, Gi=Sˇi+Ai(9) 

이다.

식 (2)i 대신에 i+1을 대입한 식과 식 (8)식 (7)에 대입하여 정리하면 식 (10)과 같이 i번째 블록 좌우측의 강성계수의 전달식을 구할 수 있다.

S^i+1=Ci+BiTVi(10) 
2.3 전달강성계수법

해석 모델의 시작 면(절면 1)의 경계조건을 식 (5)의 강성행렬(P1)로 모델링하면, 절면 1 좌측은 경계조건이 자유인 상태로 고려할 수 있다. 따라서 절면 1 좌측의 힘벡터는 영벡터(f^1=0)가 되고 식 (3)(4)로부터 절면 1 우측의 강성계수행렬을 다음 식과 같이 구할 수 있다.

Sˇ1=P1(11) 

절면 1을 제외한 나머지 절면의 강성계수행렬은 전달식인 식 (6)(10)을 통해 순차적으로 모두 구할 수 있다.

끝면(절면 n+1)의 경계조건도 식 (5)의 강성행렬(Pn+1)로 모델링하면 절면 n+1 우측은 해석적으로 경계조건이 자유(dn+10)인 상태로 고려할 수 있다. 따라서 절면 n+1 우측의 힘벡터는 영벡터(fˇn+1=0)가 되므로 식 (3)으로부터

Sˇn+1dn+1=0(12) 

가 되고, 따라서 다음 식과 같은 진동수 방정식을 유도할 수 있다.

Sˇn+1ω=0(13) 
2.4 실베스터-전달강성계수법

기존의 전달강성계수법은 진동수 방정식인 식 (13)을 통해 고유 진동수를 구할 수 있지만, 3차원 구조물의 형상이 대칭인 경우에는 고유 진동수가 중복근이 될 수 있어 계산 과정에서 이들을 빠뜨리는 문제점이 있다. 그리고 구조물의 일부 고유 진동수들이 아주 가까울 경우에도 이들을 찾아내지 못하고 빠뜨릴 가능성이 항상 존재한다. 이를 극복하기 위하여 실베스터-전달강성계수법(Sylvester-Transfer Stiffness Coefficient Method)6)이 개발되었고 상기 문제들은 완전히 해결되었다.

실베스터-전달강성계수법으로 고유 진동수를 구하는 진동수 방정식은 다음과 같다.

gω~=i=1nNGiω~+NSn+1ω~(14) 

여기서 함수 N은 입력된 행렬에 대한 음의 고유치를 세는 함수이다. 상기 식을 이용하면 임의의 시산 진동수ω~를 기준으로 그 직전의 고유 진동수의 개수를 파악할 수 있으므로 이를 이용하여 이분법(bisection method)에 적용하면 주어진 진동수 범위 내의 모든 고유 진동수를 찾아낼 수 있다.

2.5 고유 모드

진동수 방정식으로부터 고유 진동수를 구한 후, 식 (12)로부터 절면 n+1의 변위벡터 dn+1을 먼저 결정한다. 그리고 앞서 유도했던 식 (8)을 이용하면 직육면체의 끝 면에서 시작 면까지 순차적으로 모든 절면의 변위벡터를 계산할 수 있다.

그리고 변위벡터를 일차적으로 모두 구한 후에 전체 절점의 변위 중에서 가장 큰 변위 값을 찾아내어 그 값이 1이 되도록 전체 절점의 변위를 정규화한다.


3. 수치해석 및 고찰

이 연구에서 제안한 방법인 실베스터-전달강성계수법(S-TSCM)을 기반으로 매트랩(MATLAB)을 이용하여 전산 프로그램을 만든 후, 3차원 고체 구조물에 해당하는 정육면체 구조물, 직육면체 앤빌(anvil), 3차원 외팔보와 같은 3가지 계산 모델을 대상으로 자유진동 해석을 수행하고 그 결과를 고찰한다.

3.1 정육면체 구조물

Fig. 4에 나타낸 첫 번째 계산 모델은 한 변의 길이가 2 m인 정육면체 형상의 구조물이다. 이 계산 모델은 Y. Ma7) 등이 유한요소법(FEM)으로 3차원 구조물의 자유진동을 해석하면서 제시한 바 있다. 구조물은 탄성계수가 71 GPa, 푸아송 비 (Poisson's ratio)가 0.3, 밀도가 2700 kg/m3이다. 경계조건은 정육면체의 한 면만 구속이고 나머지 면은 모두 자유이다.


Fig. 4 
Cube with mesh 4×4×4

직육면체 유한요소를 이용하여 구조물을 4×4×4, 6×6×6, 8×8×8로 각각 분할하여 이 연구에서 제시한 방법(S-TSCM)으로 1차부터 6차까지 고유 진동수를 계산한 결과를 Table 1에 제시한다. 1차와 2차 고유 진동수 그리고 5차와 6차 고유 진동수는 진동수 방정식의 중복근(double root)에 해당하며 S-TSCM으로 구하는데 전혀 문제가 없었다.

Table 1 
Natural frequencies (Hz) of the three dimensional cube computed by using S-TSCM
Order 4×4×4 6×6×6 8×8×8
1 283.12 278.14 276.06
2 283.12 278.14 276.06
3 384.22 377.56 374.79
4 664.60 658.20 655.62
5 763.69 741.62 733.10
6 763.69 741.62 733.10

Table 2에 나타낸 Y. Ma7) 등이 동일한 모델에 대하여 유한요소법으로 고유 진동수를 계산한 결과와 비교해 보면, 8×8×8로 분할한 3차 고유 진동수를 제외하고 나머지 고유 진동수는 완전히 일치하였다. 따라서 이 연구에서 제시한 방법의 신뢰성을 확인할 수 있었다.

Table 2 
Natural frequencies (Hz) of the three dimensional cube computed by using FEM7)
Order 4×4×4 6×6×6 8×8×8
1 283.12 278.14 276.06
2 283.12 278.14 276.06
3 384.22 377.56 374.78
4 664.60 658.20 655.62
5 763.69 741.62 733.10
6 763.69 741.62 733.10

3.2 직육면체 앤빌

Fig. 5에 나타낸 두 번째 계산 모델은 고무(rubber) 층 위의 직육면체 앤빌(anvil)이다. 앤빌은 탄성계수가 207 GPa, 푸아송 비가 0.3, 밀도가 7,860 kg/m3이고 경계조건은 모두 자유로 모델링한다.


Fig. 5 
Anvil with mesh 6×6×2

직육면체 유한요소를 이용하여 앤빌을 3×3×1에서 18×18×6으로 각각 분할하여 이 연구에서 제시한 방법(S-TSCM)으로 강체모드를 제외하고 1차부터 4차까지 고유 진동수를 계산한 결과를 Table 3에 제시한다. Table 3에서 고유 진동수 아래에 밑줄을 그은 것은 실험치8)를 정해로 보고, 계산한 고유 진동수의 오차가 10% 이내의 범위에 들어온 것을 의미한다. 즉, 앤빌을 9×9×3 이상으로 분할하여 모델링한 경우 4차까지 고유 진동수의 백분율 오차는 10%이내에 있다는 것을 뜻한다.

Table 3 
Natural frequencies (kHz) of the anvil obtained by using S-TSCM and experiment
Order 3×3×1 6×6×2 9×9×3 12×12×4 15×15×5 18×18×6 Experiment8)
1 2.00 1.90 1.87 1.85 1.84 1.84 1.82
2 3.53 2.93 2.80 2.74 2.72 2.71 2.67
3 4.53 3.61 3.41 3.33 3.29 3.27 3.18
4 4.79 4.35 4.16 4.09 4.05 4.03 4.00

To8)의 연구에서 동일한 모델에 대해 실험을 통해 고유 진동수를 계측한 결과와 비교해 보면, S-TSCM으로 분할수를 증가하면서 계산한 고유 진동수가 실험 결과에 근접함을 확인할 수 있다.

구조물에 대한 요소 분할수를 증가하면 계산 결과가 좋아지는 것은 모델링을 이산화하는 유한요소 모델링의 일반적인 특징으로서 이 연구에서도 계산 결과를 통해 이를 확인할 수 있었다. 따라서 S-TSCM으로 3차원 구조물인 앤빌의 자유진동 해석도 충실히 수행됨을 알 수 있다.

S-TSCM으로 앤빌의 1차부터 3차까지 고유 모드를 계산하였고, 그 결과를 Fig. 6에 제시한다.


Fig. 6 
Natural modes of anvil with mesh 18×18×6

3.3 3차원 외팔보

Fig. 7에 나타낸 세 번째 계산 모델은 탄성계수가 206 GPa, 푸아송 비가 0.3, 밀도가 7,860 kg/m3사각단면을 갖는 3차원 외팔보이다.


Fig. 7 
Cantilever beam

외팔보를 y축과 z축 방향으로는 분할하지 않고 오직 길이방향(x축)으로 20, 50, 100, …, 400으로 직육면체 유한요소를 이용하여 각각 분할한 후 이 연구에서 제시한 방법(S-TSCM)으로 1차부터 10차까지 고유 진동수를 계산한 결과를 Table 4에 제시한다. 그리고 Table 4에는 외팔보를 연속계9)로 모델링하여 종·굽힘·비틀림 자유진동 해석을 통해 구한 10차 이내의 고유 진동수를 기준 값으로 제시한다.

Table 4 
Computed natural frequencies (Hz) for the three dimensional cantilever beam with mesh N×1×1
Order 20×1×1 50×1×1 100×1×1 200×1×1 300×1×1 400×1×1 Reference9) Mode9)
1 27.24 13.77 10.52 9.53 9.34 9.27 8.27 bending(z)
2 31.07 20.28 18.22 17.66 17.56 17.52 16.54 bending(y)
3 171.11 86.28 65.87 59.69 58.48 58.05 51.83 bending(z)
4 194.93 126.95 113.99 110.50 109.84 109.60 103.65 bending(y)
5 481.97 241.66 184.35 167.03 163.62 162.41 145.12 bending(z)
6 547.96 354.85 318.35 308.53 306.67 306.01 284.37 bending(z)
7 639.62 473.85 361.01 326.97 320.27 317.90 290.23 bending(y)
8 953.69 637.16 596.27 539.79 528.70 524.76 470.09 bending(z)
9 1081.13 693.73 621.53 602.14 598.47 597.17 568.74 bending(y)
10 1285.89 784.07 636.57 636.38 636.35 636.33 588.28 torsional

Table 4에서 밑줄은 연속계로 모델링하여 구한 고유 진동수를 기준으로 하여 S-TSCM으로 구한 고유 진동수의 오차가 10% 이내의 범위에 들어온 것을 의미한다. S-TSCM으로 구한 고유 진동수는 길이방향으로 요소의 분할수를 증가할수록 연속계로 모델링하여 구한 고유 진동수에 점차적으로 접근하고 있음을 알 수 있다. 그러나 1차, 3차, 5차, 8차의 고유 진동수는 길이방향으로 최대 400개의 요소로 분할함에도 불구하고 오차백분율은 10%를 상회하였다. 참고로 이들의 공통점은 z축 방향 굽힘진동의 고유 진동수이다.

외팔보를 y축 방향으로 2개의 요소로 분할하고(z축 방향의 요소는 1개) x축 방향으로 20, 50, 100, …, 400으로 직육면체 요소를 이용하여 각각 분할한 후 S-TSCM으로 1차부터 10차까지의 고유 진동수를 계산한 결과를 Table 5에 제시한다.

Table 5 
Computed natural frequencies (Hz) for the three dimensional cantilever beam with mesh N×2×1
Order 20×2×1 50×2×1 100×2×1 200×2×1 300×2×1 400×2×1 Reference9)
1 27.15 13.58 10.27 9.26 9.06 8.98 8.27
2 30.75 19.77 17.64 17.06 16.95 16.91 16.54
3 170.54 85.11 64.32 57.98 56.72 56.28 51.83
4 192.94 123.76 110.38 106.75 106.06 105.82 103.65
5 480.40 238.40 180.03 162.23 158.72 157.47 145.12
6 542.47 345.97 308.32 298.12 296.18 295.49 284.37
7 639.61 467.51 352.58 317.61 310.71 308.26 290.23
8 950.67 637.14 582.41 524.40 512.97 508.91 470.09
9 1070.63 676.55 602.08 581.95 578.12 576.77 568.74
10 1285.88 773.67 636.54 636.34 636.30 636.28 588.28

Table 5에서 외팔보를 300×2×1로 분할하여 10차까지 구한 고유 진동수는 모두 백분율 오차가 10% 이내였다. Table 4에 비해 계산한 고유 진동수의 정도(accuracy)가 상대적으로 우수함을 알 수 있다.

외팔보를 y축 방향으로 4개 그리고 z축 방향으로 2개의 요소로 분할하고 x축 방향으로 20, 50, 100, …, 400으로 직육면체 요소를 이용하여 각각 분할한 후 S-TSCM으로 1차부터 10차까지 고유진동수를 계산한 결과를 Table 6에 제시한다. 외팔보를 200×4×2로 분할하여 10차까지 구한 고유 진동수는 모두 백분율 오차가 10% 이내였다. Table 5에 비해 Table 6의 고유 진동수의 정도가 상대적으로 우수함을 알 수 있다.

Table 6 
Computed natural frequencies (Hz) for the three dimensional cantilever beam with mesh N×4×2
Order 20×4×2 50×4×2 100×4×2 200×4×2 300×4×2 400×4×2 Reference9)
1 27.01 13.29 9.88 8.82 8.61 8.53 8.27
2 30.63 19.57 17.42 16.83 16.72 16.68 16.54
3 169.69 83.31 61.90 55.26 53.94 53.46 51.83
4 192.16 122.51 108.96 105.28 104.58 104.33 103.65
5 478.03 233.39 173.25 154.63 150.92 149.60 145.12
6 540.30 342.44 304.31 293.95 291.98 291.28 284.37
7 616.70 457.74 339.34 302.76 295.48 292.89 290.23
8 946.10 613.93 560.61 499.94 487.88 483.58 470.09
9 1066.38 669.56 594.14 573.69 569.80 568.42 568.74
10 1285.87 757.61 613.23 612.98 612.92 612.90 588.28

따라서 단면의 폭과 높이에 비해 상대적으로 길이가 긴 외팔보형 3차원 구조물의 고유 진동수를 이 연구에서 제시한 방법(S-TSCM)으로 안정적으로 계산할 수 있음을 확인하였다.


4. 결 론

이 연구에서는 3차원 직육면체 요소를 도입한 실베스터-전달강성계수법으로 3차원 고체 구조물의 자유진동 해석을 수행하기 위한 알고리즘을 정식화한 후 이를 바탕으로 MATLAB으로 3차원 고체 구조물의 고유 진동수와 고유 모드를 계산하는 자유진동 해석 프로그램을 개발하였다. 그리고 이를 이용하여 3차원 고체 구조물인 정육면체 구조물, 앤빌, 그리고 3차원 외팔보 구조물의 자유진동 해석을 수행하여 고유 진동수와 고유 모드를 계산하였다.

동일한 모델에 대한 기존의 유한요소법, 실험, 해석해의 결과와 이 연구에서 제안한 방법의 결과를 비교해 보면, 이 연구에서 제안한 방법인 실베스터-전달강성계수법이 직육면체 요소를 도입한 3차원 고체 구조물의 자유진동 해석에서 신뢰성 있는 계산 결과를 제공하고 있음을 확인할 수 있었다.


References
1. M. Choi, 2003, "Free Vibration Analysis of Plate Structures Using Finite Element-Transfer Stiffness Coefficient Method", KSME International Journal, Vol. 17, pp. 805-815.
2. M. S. Choi and D. J. Yeo, 2017, "Free Vibration Analysis of Curved Beams Regarded as Discrete System Using Finite Element-Transfer Stiffness Coefficient Method", Journal of the Korean Society for Power System Engineering, Vol. 21, No. 1, pp. 37-42.
3. M. W. Suh et al., 2000, "Development of Mesh Generation Program for the Primary System of Nuclear Power Plant", Transactions of the Korean Society of Mechanical Engineers - A, Vol. 24, No. 2, pp. 386-393.
4. P. Kattan, 1990, "MATLAB Guide to Finite Elements", Springer, Cambridge, pp. 337-396.
5. M. Petyt, 1990, "Introduction to Finite Element Vibration Analysis", Cambridge University Press, Cambridge, pp. 197-203.
6. M. Choi, T. Kondou and Y. Bonkobara, 2012, "Development of Free Vibration Analysis Algorithm for Beam Structures by Combining Sylvester’s Inertia Theorem and Transfer Stiffness Coefficient Method", Journal of Mechanical Science and Technology, Vol. 26, No. 1, pp. 11-19.
7. Y. Ma, Y. Yang and G. Sun, 2019, "Three Dimensional Vibration Analyses using an eight-node hexahedral element with continuous nodal stress", Computers and Structures, Vol. 212, pp. 58-71.
8. C. W. S. To, 1982, "Application of the Finite Element Method for the Evaluation of Velocity Response of Anvils", Journal of Sound and Vibration, Vol. 84, No. 4, pp. 529-548.
9. S. S. Rao, 2005, "Mechanical Vibrations (Fourth Edition in SI Units)", Prentice Hall, Singapore, pp. 588-656.