오목 멤브레인의 전 주파수 대역 고유치 추출을 위한 3영역 분할법 기반 무차원동영향함수법 개발: 제1부(이론 정립과 검증)

2.1 3영역 분할법

3영역 분할법 기반 무차원동영향함수법(NDIF법)은 오목한 형상의 멤브레인을 고정밀도로 해석하기 위해 개발된 기법으로, 해석 대상의 전체 영역을 세 개의 볼록 영역으로 분할하여 해석하는 방법론을 채택한다. 이 연구에서는 Fig. 1과 같이 오목 형상 멤브레인을 세 개의 볼록 영역 D₁, D₂, D₃로 분할하며, 각각의 영역은 외부 경계(실선)와 접촉 경계(점선)로 둘러싸여 있다 가정한다.

Fig. 1

Arbitrarily shaped concave membrane divided into three convex domains D1, D2, and D3

Fig. 1과 같이 볼록 영역 D₁의 경계는 외부 경계(실선) Γ₁과 접촉 경계(점선) Γ_a+Γ_b로 구성된다. 여기서 외부 경계는 변위 또는 기울기가 경계조건으로 주어지는 경계를 의미하며, 접촉 경계는 분할된 다른 볼록 영역과 접촉되어 있는 경계를 의미한다. 경계 Γ₁, Γ_a, Γ_b는 각각 절점 $$ P_1^{\left(1\right)},P_2^{\left(1\right)},\cdots ,P_{N_1}^{\left(1\right)}$$, $$ P_1^{\left(a\right)},P_2^{\left(a\right)},\cdots ,P_{N_a}^{\left(a\right)}$$, $$ P_1^{\left(b\right)},P_2^{\left(b\right)},\cdots ,P_{N_b}^{\left(b\right)}$$로 이산화된다. 마찬가지 방식으로, 두번째 볼록 영역 D₂의 경계는 외부 경계 Γ₂과 접촉 경계 Γ_a로 구성된다. 경계 Γ₂는 $$ P_1^{\left(2\right)},P_2^{\left(2\right)},\cdots ,P_{N_2}^{\left(2\right)}$$로 이산화된다. 마지막으로 세 번째 볼록 영역 D₃의 경계는 외부 경계 Γ₃과 접촉 경계 Γ_b로 구성된다. Γ₃는 절점 $$ P_1^{\left(3\right)},P_2^{\left(3\right)},\cdots ,P_{N_3}^{\left(3\right)}$$로 이산화된다.

2.2 영역 내부의 변위 함수 설정

분할된 영역 D₁ 내부의 변위 W₁(r)는 해당 영역을 둘러싼 외부 경계 Γ₁와 접촉 경계 Γ_a+Γ_b에 놓인 절점들에서 정의된 무차원동영향함수들의 선형 결합 형태로 식 (1)과 같이 가정된다⁽¹⁾.

여기서 벡터 r은 영역 내부 한점에 대한 위치 벡터, ∧는 주파수 파라미터, J₀(⋯)는 멤브레인의 무차원동영향함수(제1종 0차 베셀함수)이다. $$ {\mathbf{r}}_i^{\left(1\right)}$$은 Γ₁에서 i번째 절점에 대한 위치벡터, $$ {\mathbf{r}}_j^{\left(a\right)}$$은 Γ_a에서 j번째 절점에 대한 위치벡터를 나타낸다. 또한, $$ A_i^{\left(1\right)}$$, $$ A_j^{\left(a\right)}$$, $$ A_j^{\left(b\right)}$$는 해당 절점에서 정의된 무차원동영향함수의 기여도 계수이다.

분할된 영역 D₂와 D₃의 내부 변위도 외부 경계와 접촉 경계에 놓인 절점에 대해 정의된 무차원동영향함수들의 선형 결합으로 식 (2)와 식 (3)과 같이 각각 가정된다.

식 (2)와 식 (3)에 주어진 변수들에 대한 설명은 식 (1)에 대한 설명을 참고하면 알 수 있기에 생략한다.

2.3 영역별 외부 변위와 접촉 변위 가정

영역 D₁의 외부 경계 절점 $$ P_1^{\left(1\right)},P_2^{\left(1\right)},\cdots ,P_{N_1}^{\left(1\right)}$$에 부여된 변위 조건은 식 (4)와 같이 가정한다.

여기서 $$ U_k^{\left(1\right)}$$는 외부 경계 Γ₁에 놓인 k번째 절점 $$ P_k^{\left(1\right)}$$에 부여된 변위를 의미한다. 마찬가지 방식으로 영역 D₂와 D₃의 외부 경계 Γ₂와 Γ₃에 놓인 절점에 부여된 변위 경계 조건은 식 (5), 식 (6)과 같이 각각 가정한다.

여기서 $$ U_k^{\left(2\right)}$$와 $$ U_k^{\left(3\right)}$$는 외부 경계 Γ₂와 Γ₃ 놓인 절점 $$ P_k^{\left(2\right)}$$와 $$ P_k^{\left(3\right)}$$에 부여된 변위를 각각 의미한다.

영역 D₁의 접촉 경계 Γ_a의 절점 $$ P_1^{\left(a\right)},P_2^{\left(a\right)},\cdots ,P_{N_a}^{\left(a\right)}$$과 Γ_b의 절점 $$ P_1^{\left(b\right)},P_2^{\left(b\right)},\cdots ,P_{N_b}^{\left(b\right)}$$에서의 변위는 식 (7), 식 (8)과 같이 각각 가정한다.

여기서 $$ U_m^{\left(a\right)}$$와 $$ U_m^{\left(b\right)}$$는 각각 접촉 경계 Γ_a와 Γ_b에 놓인 m번째 절점에서의 변위를 의미한다.

한편, 영역 D₂의 관점에서 접촉 경계 Γ_a의 절점에서의 변위는 식 (9)와 같이 가정한다.

마찬가지로, 영역 D₃의 관점에서 접촉 경계 Γ_b의 절점에서의 변위는 식 (10)과 같이 가정한다.

2.4 시스템 행렬 구성

영역 D₁의 내부 변위 식 (1)을 변위 경계 조건 식 (4)와 접촉 변위 식 (5), 식 (6)을 대입하면 식 (11)을 각각 얻을 수 있다.

식 (11) ~ 식 (13)은 식 (14) ~ 식 (16)과 같이 행렬식으로 각각 표현할 수 있다.

여기서 시스템 행렬 SM_αβ의 p번째 행, q번째 열 성분은 $$ J_0\left(\wedge \left|{\mathit{r}}_p^{\left(\alpha \right)}-{\mathit{r}}_q^{\left(\beta \right)}\right|\right)$$이고, 기여도 벡터 A₁, A_a, A_b의 성분은 각각 $$ A_i^{\left(1\right)}$$, $$ A_j^{\left(a\right)}$$, $$ A_j^{\left(b\right)}$$로 주어지며, 변위 벡터 U₁, U_a, U_b의 성분은 각각 $$ U_k^{\left(1\right)}$$, $$ U_m^{\left(a\right)}$$, $$ U_m^{\left(b\right)}$$로 주어진다.

영역 D₂의 내부 변위 식 (2)를 변위 경계 조건 식 (5)와 접촉 변위 식 (9)에 대입한 후 행렬식으로 변환하면 식 (17), 식 (18)을 각각 얻을 수 있다.

마찬가지로 영역 D₃의 내부 변위 식 (3)을 변위 경계 조건 식 (6)과 접촉 변위 식 (10)에 대입한 후 행렬식으로 변환하면 식 (19), 식 (20)을 각각 얻을 수 있다.

Fig. 1에서와 같이 접촉 경계 Γ_a는 D₁과 D₂의 공통 접촉 경계이므로, 식 (15)와 식 (18)이 같다는 변위 연속 조건을 고려하면 식 (21)을 얻을 수 있다.

또한, 접촉 경계 Γ_b는 D₁과 D₃의 공통 접촉 경계이므로, 식 (16)과 식 (20)이 같다는 변위 연속 조건을 고려하면 식 (22)를 얻을 수 있다.

다음으로, 식 (1) ~ 식 (3)을 접촉 경계 Γ_a의 법선 방향으로 미분을 한 후, Γ_a에서의 기울기 연속 조건을 고려하면 식 (23)의 행렬식을 얻을 수 있다.

마찬가지 방식으로, Γ_b에서의 기울기 연속 조건을 고려하면 식 (24)의 행렬식을 얻을 수 있다.

여기서 PM_αβ = ∂ SM_αβ/∂ n을 의미하며, 시스템 행렬 PM_αβ의 p번째 행, q번째 열 성분은 ∂$$ J_0\left(\wedge \left|{\mathbf{r}}_p^{\left(\alpha \right)}-{\mathbf{r}}_q^{\left(\beta \right)}\right|\right)$$/ ∂n이다.

2.5 최종 시스템 행렬 구성 및 고유치 추출 판별식

2.4절에서 유도한 7개의 시스템 행렬식 식 (14), 식 (17), 식 (19), 식 (21) ~ 식 (24)로부터 최종 시스템 행렬을 추출하기 위해 이전 연구와 마찬가지로 행렬 축약법을 먼저 이용하였다⁽²⁾. 고정단 경계 조건 U₁ = U₂ = U₃ = 0을 식 (14), 식 (17), 식 (19)에 적용하면 식 (25) ~ 식 (27)을 각각 얻을 수 있다.

식 (25) ~ 식 (27)을 식 (21) ~ 식 (24)에 대입하여 정리하면 식 (28)과 같이 축약된 최종 시스템 행렬식을 얻을 수 있다.

여기서 SM_r(Λ)와 Dr은 해석 대상 오목 멤브레인에 대한 최종 축약 시스템 행렬과 축약 기여도 벡터를 의미하며 식 (29)의 식들에 의해 주어진다.

여기서 SM_r(Λ)의 성분은 식 (30) ~ 식 (41)과 같다.

식 (28)로부터 고유치 추출 판별식 식 (29)를 얻을 수 있으며 식 (42)의 해가 고유치에 해당된다.

여기서 축약된 최종 시스템 행렬 SM_r(Λ)은 정사각 행렬이며 그 크기는 접촉 경계 절점 개수의 두배인 2(n_a + n_b)에 불과하다.

3.1 제안된 이론 검증

식 (42)가 전 주파수 대역에서 정확한 고유치를 제공하는지를 확인하기 위해 엄밀해를 가지는 직사각형 멤브레인이 고려된다. Fig. 2(a), Fig. 2(b)는 의도적으로 3개 영역으로 분할된 가로 1.2 m, 세로 0.9 m 직사각 멤브레인을 보여주며, 이들은 각각 40개와 52개의 절점으로 이산화되었다. 식 (42)를 이용하여 이들 멤브레인에 대한 판별식 곡선을 구하면 Fig. 3과 같다. Fig. 3에서 40 절점 곡선과 52 절점 곡선 둘 다 저주파수 구간에서 시스템 행렬이 발산하여 저차 고유치들이 추출되지 않음을 확인할 수 있다. Fig. 3의 두 판별식 곡선으로부터 추출된 고유치들은 Table 1에 요약되었다⁽⁵⁾.

Fig. 2

Rectangular membranes divided into 3 domains and discretized with 40 and 52 nodes, respectively

Fig. 3

Determinant curves of the rectangular membrane for 40 and 52 nodes plotted by Eq. (42)

Table 1

Eigenvalues of the rectangular membrane obtained by NDIF method (Eq. (42)), the exact solution, and FEM (ANSYS)

Table 1에서, 40 절점의 경우 1차 고유치가, 52 절점의 경 1차 ~ 5차 고유치가 미추출되었음이 확인된다. 한편, 성공적으로 추출된 고유치들은 엄밀해와 정확히 일치했음을 확인할 수 있는데, 1089 절점을 사용한 FEM보다 NDIF법이 훨씬 정확하다고 말할 수 있다.

이 논문에서는 절점의 개수가 증가하면 저차 고유치 미추출 문제가 발생하는 식 (42)의 약점을 개선하는 방안이 제안된다. 그 첫 번째 과정이 SM_r(Λ)곡선의 판별식을 계산하기 위해 이전 논문의 연구 결과를 확장하여 식 (43)을 이용하는 것이다⁽²⁾.

여기서 R(Λ)는 SM_r(Λ)의 랭크(rank)이며 주파수 파라미터의 함수이다. 그리고 Ω_k(Λ)는 SM_r(Λ)에 대한 대수 고유치 문제에서의 k번째 고유치를 의미한다⁽²⁾. 한편, 식 (43)이 유용성을 가지려면 R(Λ)가 상수가 아니어야 한다⁽²⁾. 그래서 Fig. 2의 사각 멤브레인의 SM_r(Λ)에 대한 R(Λ) 곡선을 Fig. 4와 같이 구해보았다. Fig. 4에서 40 절점과 52 절점 두 경우 모두 랭크가 20과 24로 각각 일정하게 유지됨이 확인되며, 이러한 사실로부터 SM_r(Λ)의 판별식 값을 구할 때 식 (43)은 유용하지 않다고 결론지을 수 있다.

Fig. 4

Rank of SMr for the rectangular membrane discretized with 40 and 52 nodes as shown in Fig. 2

3.2 저차 고유치 추출을 위한 새로운 방안

이 연구에서는 식 (43)을 사용하기 위해 최종 시스템 행렬을 새로운 방법으로 구하는 방안을 시도한다. 이 방법의 핵심은 시스템 행렬을 구성할 때 로컬 시스템 행렬을 역행렬화시키는 과정을 배제하는 것이다. 그래서, 식 (14), 식 (17), 식 (19), 식 (21) ~ 식 (24)을 변형하지 않고 그대로 하나의 행렬식으로 나타내는 방법을 채택한다. 그러면 식 (44)를 얻을 수 있다.

여기서 TM(Λ), D, U는 해석 대상 오목 멤브레인에 대한 최종 시스템 행렬, 기여도 벡터, 외각 경계 변위 벡터를 의미하며, 식 (45) ~ 식 (47)에 의해 주어진다.

해석 대상 오목 멤브레인이 고정단 경계를 가진다면 U₁ = U₂ = U₃ = 0을 만족하므로 식 (44)는 식 (48)과 같이 된다.

식(48)에서 D ≠ 0을 만족하기 위해선 시스템 행렬의 판별식이 식 (49)와 같이 0이 되어야 한다.

해석 대상 오목 멤브레인의 고유치는 식 (49)를 만족하는 주파수 파라미터(Λ) 값으로부터 구할 수 있다. 이제 식 (49)의 판별식을 계산하기 위해 식 (43)에 식 (45)의 TM(Λ)을 대입하면 식 (50)과 같이 된다.

식 (50)의 유용성을 확인하기 위해 TM(Λ)의 랭크 곡선을 그려 본 결과 주파수 파라미터에 따라 랭크 값의 변화가 있음이 Fig. 5와 같이 확인되었다. 마지막으로 식 (50)을 이용하여 판별식 곡선을 그린 후 평탄화 과정을 거치면 Fig. 6과 같은 최종 판별식 곡선을 얻을 수 있다^(2~4). 참고로 평탄화 과정에 대한 구체적인 연구 결과는 논문 2부에서 제시된다.

Fig. 5

Rank of TM for the rectangular membrane discretized with 40 and 52 nodes as shown in Fig. 2

Fig. 6

Determinant curves of the rectangular membrane for 40 and 52 nodes plotted by Eq. (50)

Fig. 6의 판별식 곡선으로부터 얻은 고유치들은 판별식 곡선에 직접 표기하였다. 이들 값들을 살펴 보면, 저주파수 영역을 포함한 전 주파수 대역에서 모든 고유치들이 추출되었음이 확인되며, Table 1의 엄밀해와 정확히 일치함도 확인된다.