무요소법을 이용한 임의 형상 평판의 고정확도 고유치 추출 기법

Extraction of Accurate Eigenvalues of Plates Using a Meshless Method

Abstract

The Kansa method, which is used for various free vibration problems of arbitrarily shaped plates including membranes, discretizes the domain of a plate using only nodes without elements unlike FEM. The method requires a small amount of computation relative to FEM thanks to this discretization scheme but has limit in the accuracy of its solution. This paper reveals the reason of the limit and, to overcome the limit, proposes the practical method of calculating the singularity of a system matrix and extracting accurate natural frequencies. Case studies for a rectangular plate and an arbitrarily shaped plate validate the proposed method.

1. 서 론

이 논문의 저자는 임의 형상 평판과 박막(membrane)의 고유 진동수와 모드 형상을 정확히 구하기 위한 방법으로 NDIF법(non-dimensional dynamic influence function method)을 1999년 처음으로 개발 하였다(1). NDIF법은 무요소(meshless method)의 일종으로 관심 영역의 경계를 요소를 사용하지 않고 노드만으로 이산화하며, 유한요소법(2)과 경계요소법(3)에 비해 해의 정확도 면에서 탁월한 우수성을 가지고 있음이 저자의 논문들(4~10)을 통해 증명되었다. 그러나 NDIF법은 많은 노드가 사용될 경우 시스템행렬의 특이성(singularity)이 커져서 해가 발산하는 단점이 존재한다. 반면에, NDIF 법과는 달리 관심 영역의 내부와 경계 모두를 노드로 이산화하는 무요소법인 Kansa법(Kansa method)(11~15)은 노드가 증가하더라도 수렴된 해를 제공한다. 그러나 Kansa법은 NDIF법에 비해 해의 정확도가 떨어지는 단점이 존재한다. 뿐만 아니라, Kansa법에서 사용되는 기저함수인 RBF(radial basis function)의 형상 파라미터 값에 따라 해의 정확도가 크게 변하기 때문에, Kansa법은 형상 파라미터 값을 경험적으로 결정해야만 하는 단점도 가진다.

이 논문에서는 Kansa법에서 해의 정확도와 수렴성 저하가 발생하는 원인을 분석하고, 해의 정확도를 올릴 수 있는 방법을 제안한다. 그리고 RBF 형상 파라미터 값을 결정하는 데에 도움이 되는 유용한 방안도 제안한다.

 

2. 고정확도 고유 진동수 추출 기법

2.1 Kansa법(Kansa method)

이 절에서는 Kansa법(11)을 이용한 단순지지 평판의 고유 진동수 추출 이론 정식화에 대해 설명한다. 참고로 고정단 또는 자유단 평판은 향후 논문에서 다루어진다. 단순지지 경계조건을 가진 임의 형상 평판의 운동 지배방정식은 다음과 같다(16).

여기서 W는 평판의 진동 변위, Λ는 주파수 파라미터, ρ는 평판의 밀도, ω는 각주파수, DE는 판강성(flexural rigidity of plate)을 나타낸다. 그리고 단순지지 경계 조건은 변위와 모멘트가 0인 두 개의 식으로 다음과 같이 정의된다.

여기서 Γ는 평판의 경계를 의미한다. 한편, 단순지지 평판의 형상이 다각형일 경우, 지배방정식과 경계조건은 다음과 같이 간단히 나타낼 수 있다.

Fig. 1에서 평판 상의 임의의 점 P에서의 진동 변위는 식 (7)과 같이 RBF의 선형 결합으로 가정한다.

Fig. 1Arbitrarily shaped plate with boundary Γ

여기서 은 임의 점 P에 대한 위치벡터, 는 k번째 노드 Pk에 대한 위치벡터, 는 P와 Pk 사이의 거리, 는 노드 Pk가 함수의 중심인 RBF이다. 그리고 Ak는 RBF의 기여도를 나타내며, N은 평판의 내부와 경계를 이산화하기 위해 사용된 전체 노드수를 의미한다. 참고로, 여러 가지 종류의 RBF가 존재하나 이 논문에서는 ϕ(r) = e−cr2이 사용되었으며, c는 RBF의 형상을 결정짓는 형상 파라미터이다.

해석 대상 평판의 내부는 ND개의 노드로, 경계는 NB개의 노드로 이산화 되었고, 전체 노드수는 N = ND+NB이라고 가정하자. 그리고 ND개의 내부 영역 노드들에 대한 위치벡터는 이고, NB개의 경계 노드들에 대한 위치벡터는 , 이라고 가정한다. 내부 영역 노드들에서는 진동 변위가 지배방정식 식 (5)를 만족해야 하므로, 식 (7)을 식 (5)에 대입하면 다음과 같이 정리된다.

또한, 경계 노드에서의 진동 변위가 경계조건 식 (6)을 만족해야하므로, 식 (7)을 식 (6)에 대입하면 다음과 같이 된다.

이제 식 (8)과 식 (9)는 다음과 같이 하나의 행렬 식으로 나타낼 수 있다.

여기서 기여도 벡터 는 다음과 같으며

크기가 N × N인 정사각 시스템 행렬 과 는 다음과 같다.

여기서 시스템 행렬들의 성분은 다음과 같이 주어진다.

시스템 행렬식 식 (10)은 다음과 같이 대수 고유치 문제로 정식화가 될 수 있다.

여기서 λ = 1/Λ2이다. 대수 고유치 문제 상용 솔버(solver)를 이용하면, 식 (18)로부터 고유 진동수와 고유 모드를 용이하게 구할 수 있다.

2.2 Kansa법 오차 원인 규명

Kansa법의 정확도를 확인하기 위해 엄밀해가 존재하는 단순지지 사각 평판(가로 1.2 m, 세로 0.9 m)을 고려하였다. Fig. 2와 같이 사각 평판은 가로 방향 9개 노드와 세로 방향 7개 노드로 총 64개의 노드로 이산화되었다.

Fig. 2Rectangular plate discretized with 63 nodes (unit : m)

Table 1은 형상 파라미터 c = 5를 사용하여 Kansa법에 의해 구해진 고유 진동수 추출 결과를 보여준다. 엄밀해와 비교해보면 최대 0.48 %의 오차가 발생됨이 확인된다. 반면에 NDIF법(7)은 오차가 없이 엄밀해(17)와 일치한 결과를 보여준다. 참고로, Kansa법에 대한 많은 연구 결과들이 존재하나 이 방법의 오차 원인에 대한 자세한 연구 결과는 없는 것으로 파악된다.

Table 1Natural frequencies of the rectangular plate using c = 5; parenthesized values denote percentage errors

Kansa법의 정확성은 형상 파리미터(c) 값의 크기에 따라 민감하게 변한다. Li(15)는 c의 크기를 크게하는 것이 해의 정확성을 올릴 수 있다고 설명하였으나, 이 연구에서는 Table 2와 같이 Li의 연구와 일치하지 않은 결과가 얻어졌다. Table 2는 c의 값을 5.0부터 0.5단위로 감소시켰을 때 첫 번째 고유진동수 추출 결과를 보여준다. c의 크기가 작아짐에 따라 고유 진동수 오차가 점점 작아짐을 확인할 수 있으며, c = 3일 때는 고유 진동수가 추출되지 않았다. 이러한 사실은 Li(15)의 연구 결과와 상반된 결과이며, 가장 정확한 고유치를 제공할 수 있는 c의 크기는 3 < c < 3.5 범위에 있을 것이라고 추정할 수 있다.

Table 2First natural frequency of the rectangular plate by the Kansa method when c is decreased by 0.5; parenthesized values denote percentage errors

보다 정확한 c의 값을 구하기 위해, c의 값을 0.1 단위로 감소시키면서, 1.5 ≤ c ≤ 6.5인 범위에서 첫 번째 고유 진동수를 추출하였다. 그 결과는 Fig. 3의 그래프에 제시되었다. 그래프를 살펴보면, 4.5 ≤ c < 6.5인 구간에서는 c값이 작아짐에 따라 엄밀해(exact)와의 오차가 점점 작아지면서 수렴함을 확인할 수 있다. 반면에 3.4 ≤ c < 4.5 구간에서는 고유 진동수가 더 이상 수렴하지 않고 심하게 진동함을 확인할 수 있다. 그리고 c < 4.5인 구간에서는 고유 진동수가 추출되지 않았다.

Fig. 3First natural frequency of the rectangular plate by the Kansa method when c is decreased by 0.1 from 6.5 to 1.5

Fig. 3의 그래프에서 해의 수렴성이 유지되는 상태에서 가장 정확한 고유 진동수를 제공하는 c의 값은 4.5인 것으로 확인되며, 이때의 고유 진동수 값은 23.86 Hz로 엄밀해와의 오차가 0.29 %이다.

한편, Fig. 3의 그래프를 살펴보면, c < 4.5인 구간에서도 고유 진동수 값이 진동하지 않고 수렴성이 보장된다면 더 정확한 고유 진동수가 추출될 수 있을 것으로 생각된다. 그래서 c < 4.5 구간에서 고유 진동수가 진동하는 원인에 대해 분석하였다. Kansa법의 최종 결과식 식 (18)을 살펴보면, 시스템 행렬 의 역행렬을 구하는 과정이 포함되어 있다. 일반적으로 특이성이 큰 행렬의 경우 그 역행렬은 많은 오차를 동반한다. 그래서 c값의 크기에 따른 의 특이성의 크기를 계산해보았다. 그 결과는 Fig. 4와 같다. 참고로, 행렬의 특이성 크기를 구하기 위해 SVD(singular value decomposition)을 수행하였다.

Fig. 4Singularity of the system matrix of the rectangular plate when c increases

Fig. 4에서 수직축은 특이성의 역수에 해당하며, 대략 c값이 4보다 커야만 시스템 행렬 의 특이성이 급격히 작아져서, 역행렬을 구하는 과정에서 오차가 발생하지 않을 것이라는 것을 예측할 수 있다. 이러한 사실로부터, Fig. 3에서 고유 진동수 값이 진동하는 원인은, c값이 작아짐에 따라 의 특이성이 커지기 때문이라고 결론 내릴 수 있다.

2.3 오차 최소화 고유 진동수 추출 기법

Kansa법에서 보다 정확한 고유 진동수를 추출할 수 있는 새로운 기법이 이 절에서 제안된다. 2.2절에서 규명한 바와 같이 Kansa법의 오차는 식 (18)에서 시스템 행렬 의 역행렬을 구하는 과정에서 발생한다. 그래서 고유 진동수를 구하기 위하여 식 (18) 대신에 식 (10)이 사용된다. 식 (10)을 다시 적으면 다음과 같다.

여기서

식 (19)에서 기여도 벡터 가 존재하기 위해서는 행렬 의 역행렬이 존재하지 않아야 한다. 결과적으로, 식 (21)과 같이 의 판별식(determinant)이 0이 되게 하는 주파수 파라미터 Λ을 찾는 방법에 의해 고유 진동수가 구해질 수 있다. 참고로, 고유 진동수는 식 (2)를 이용하면 위에서 구한 Λ로부터 계산될 수 있다.

위에서 제안한 방법의 타당성을 검증하기 위하여, Fig. 2의 직사각형 평판에 제안된 방법을 적용하였다. 형상 파라미터의 값을 6.5부터 1.5까지 0.1씩 줄여가면서 첫 번째 고유 진동수를 제안된 방법에 의해 추출하였다. 그 결과는 Fig. 5와 같다.

Fig. 5First natural frequency of the rectangular plate by the proposed method when c is decreased by 0.1 from 6.5 to 1.5

Fig. 5를 살펴보면, 제안된 방법(proposed)은 c < 4.5인 구간에서도 고유 진동수가 진동하지 않고 엄밀해(exact)에 정확히 수렴하는 결과를 보여준다. 반면에 Kansa법(Kansa)은 고유 진동수가 진동하다가 c값이 3.5보다 더 작아지면 고유 진동수가 추출되지 않음을 재확인할 수 있다. 이상의 결과와 같이, Kansa법의 오차 원인은 고유 진동수 추출 과정에서의 시스템 행렬의 역행렬을 구하는 과정에서 발생하는 것으로 확인되었으며, 제안된 방법은 그 오차를 최소화할 수 있는 방법임이 검증되었다.

 

3. 예제 연구

이 절에서는 엄밀해가 존재하는 직사각형 평판과 엄밀해가 존재하지 않는 임의 형상 사각형에 대한 예제 연구가 수행되었다. 예제에 사용된 평판의 두께는 0.005 m, 밀도는 7 800 kg/m3, 푸아송 비는 0.3, 영률은 210 GPa이다.

3.1 직사각형 평판

Fig. 2의 직사각형 평판에 제안된 방법을 적용하여 저차의 6개 고유 진동수를 추출하고자 한다. 먼저 c값을 결정해야 하므로, Fig. 4의 특이성 곡선을 참조하여 c = 2.0로 결정하였다. c = 2.0로 정한 이유는 특이성이 큰 c구간에서 제안된 방법이 고정확도 고유 진동수를 제공함이 앞에서 확인되었기 때문이다. 다음으로, 식 (21)의 의 대수 값(logarithm value)을 4 ≤ Λ ≤ 11 인 구간에서 그래프로 그렸다. 그 결과는 Fig. 6과 같다. 참고로, Λ의 범위는 0부터 20까지 충분히 넓은 범위에서, Λ의 값을 (충분히 큰 값인) 0.01씩 증가하면서 그래프를 그리는 과정을 통해 관심 주파수 대역 4 ≤ Λ ≤ 11을 찾아내었다. 그리고 관심 주파수 대역에서 Λ의 값을 (충분히 작은 값인) 0.0001씩 증가시키면서 Fig. 6과 같이 그래프를 그렸다.

Fig. 6Logarithm graph for det(SM) of the rectangular plate plotted by the proposed method using c = 2.0

Fig. 6의 곡선에서 6개의 극솟값에 해당되는 주파수 파라미터 값들(Λ1 ∼ Λ6)은 각각 4.3631, 6.2914, 7.4538, 8.6015, 8.7238, 10.5123이며, 이들 값들을 식 (2)를 이용하여 고유 진동수로 환산한 값들이 Table 3에 제시되었다. 제안된 방법에 의해 구해진 고유 진동수를 엄밀해(17)와 비교해보면, 제안된 방법이 Kansa법에 비해 확연히 정확한 결과를 제공함이 확인된다. 참고로, Kansa법에서 c = 4.5를 사용한 이유는 Fig. 5에서 Kansa법이 수렴된 결과를 제공하는 c의 최솟값이 4.5인 것으로 확인되었기 때문이다.

Table 3Natural frequencies of the rectangular plate obtained by the Kansa method, the exact solution, and the proposed method; parenthesized values denote percentage errors

3.2 임의 형상 사각형 평판

Fig. 7과 같이 72개의 노드들로 이산화된 임의 형상 사각형 평판이 이 예제에서 고려된다. 먼저, Kansa법을 이용하여 고유 진동수를 구하기 위하여 적절한 c값을 구하고자 한다.

Fig. 7Dimensions of the arbitrarily shaped quadrilateral plate discretized with 72 nodes (unit : m)

이를 위해, 의 특이성 곡선을 Fig. 8과 같이 그렸다. 특이성 곡선을 살펴보면, c > 0.8 구간에서 특이성이 확연히 낮아지는 것을 확인할 수 있다. 그래서, c값을 0.8부터 10까지 증가시키면서 첫 번째 고유 진동수 추출하였으며 그 결과는 Fig. 9와 같다. Fig. 9에서 0.8 < c < 1 구간에서는 고유 진동수 값이 상하로 진동하면서 수렴하지 못하는 것을 확인할 수 있다. 반면에 c > 1.0인 구간에서는 고유 진동수가 수렴되고 있음이 확인된다. 앞의 Fig. 5에서 확인한 바와 같이, Kansa법에서 가장 오차가 적은 고유 진동수를 추출하기 위해선 수렴된 고유 진동수를 제공하는 가장 작은 c값을 찾으면 된다. 그래서 c = 1.2 정도로 선택해서 Kansa법에 의한 고유 진동수를 추출하였다. 그 결과는 Table 4(Kansa method)에 제시되었다.

Fig. 8Singularity of the system matrix of the quadrilateral plate when c increases

Fig. 9First natural frequency of the quadrilateral plate by the Kansa method when c is decreased from 10 to 0.9

Table 4Natural frequencies of the arbitrarily shaped quadrilateral plate obtained by the proposed method and other methods; parenthesized values denote percentage errors

Table 4에서 Kansa법에 의해 구한 고유 진동수의 정확도를 확인하기 위하여 NDIF법과의 오차를 계산하였다. 오차 범위가 0.17 % ~ 1.18 %인 것으로 확인된다. 참고로 현재 가장 정확한 고유 진동수 추출 방법이 NDIF법이기 때문에, NDIF법에 대한 오차를 계산하였다.

다음으로, 제안된 방법을 이용하여 고유 진동수가 구해진다. 이를 위해, Fig. 10과 같이 의 대수 값(logarithm value)을 1.5 ≤ Λ ≤ 4.5인 구간에서 그래프로 그렸다. 그래프에서 가 최소가 되게 하는 고유치 파라미터 값(Λ1 ∼ Λ6)은 각각 1.8984, 2.8816, 3.0817, 3.8007, 4.0015, 4.3343이다. Λ1 ∼ Λ6을 식 (2)를 이용하여 고유 진동수로 환산한 결과는 Table 4(proposed method)에 제시되었다. 제안된 방법을 이용하여 고유 진동수를 구할 때 c값은 0.8로 사용되었다. 그 이유는 Fig. 7에서 의 특이성이 c = 0.8 지점에서 크기 때문이다. 참고로, 0.8 보다 더 작은 c를 적용할 경우에도 Kansa법에 비해 더 정확한 고유 진동수를 추출됨이 확인되었다.

Fig. 10Logarithm graph for det(SM) of the quadrilateral plate plotted by the proposed method using c = 0.8

Table 4에서 제안된 방법은 Kansa법 보다 훨씬 더 정확한 고유 진동수를 제공함이 확인된다. 또한, 단지 72개의 노드를 사용한 제안된 방법은 1812개의 노드를 사용한 FEM(ANSYS)(7) 보다 훨씬 더 정확한 결과를 제공함도 확인된다.

 

4. 결 론

이 논문에서는 Kansa법의 정확도를 올릴 수 있는 새로운 고유 진동수 추출 기법이 제안되었으며, 예제 연구를 통해 그 타당성이 검증되었다. 보다 세부적인 연구 결과로, 기존 논문에서는 Kansa법의 형상 파라미터의 값을 크게 하는 것이 해석 정확도를 올리는데 유리하다고 기술하였으나, 제안된 방법에서는 형상 파라미터 값을 최소화하는 것이 해석 정확도 향상에 유리함이 확인되었다. 또한, Kansa법을 이용해서 정확한 고유 진동수를 추출하기 위해 형상 파라미터 값의 크기를 결정해야 하는데, 이 값을 결정하는 데에 기준이 될 수 있는 방안으로, 의 특이성 곡선을 이용하는 방법을 제안하였다.