임의 형상 음향 공동의 전 주파수 대역 고유치 추출을 위한 개선된 NDIF법 개발: 제1부(이론 정립)

2.1 지배방정식과 경계조건

Fig. 1과 같은 임의 형상 음향 공동의 고유치 해석에 대한 지배방정식은 식 (1)과 같이 헬름홀츠 방정식(Helmholtz equation)의 형태를 가진다⁽¹²⁾.

Fig. 1

Arbitrarily shaped, 2-D acoustic cavity discretized with boundary nodes P1,P2,...,PN

여기서 r은 음향 공동 내부의 한 점 P에 대한 위치 벡터를 의미(Fig. 1)하고, p(r)은 음향 공동 내부의 한 점 P에서의 음압(sound pressure)을 나타내며, k는 파수(wavenumber)이다. 그리고 강체벽 음향 공동의 경계 조건은 다음과 같이 주어진다⁽¹²⁾.

여기서 r_Γ는 음향 공동 경계 상의 한 점에 대한 위치 벡터이며, n은 경계에서의 법선 방향을 의미한다(Fig. 1).

2.2 음향 공동 내부 음압 가정과 경계 조건 이산화

NDIF법을 적용하기 위하여 Fig. 1과 같이 음향 공동의 경계는 N개의 노드들 P₁,P₂,...,P_N으로 이산화 된다. 그리고 음향 공동 내부 음압 p(r)은 각각 경계 노드들에서 정의된 무차원 동영향 함수들의 선형 결합으로 다음과 같이 가정된다⁽⁷⁾.

여기서 J₀는 제1종 0차 베셀 함수를 나타내며, 벡터 r_s는 경계 Γ에 위치한 s번째 노드에 대한 위치 벡터를 의미하고, A_s는 s번째 노드에서 정의된 무차원 동영향 함수의 기여도를 나타내는 기여도 계수이다.

음향 공동의 경계를 따라 연속적으로 주어진 강체벽 경계 조건 식 (2)는 다음과 같이 경계 노드에 대한 식으로 이산화된다.

여기서, r_i와 n_i는 강체벽 경계에 놓인 노드 P_i에 대한 위치 벡터를 각각 의미한다.

2.3 시스템 행렬 추출

앞에서 가정한 음압 식 (3)을 이산화된 강체벽 경계조건 식 (4)에 대입하면 다음의 식을 얻을 수 있다.

식 (5)에서 법선 방향 n_j에 대한 미분을 수행하면, 식 (5)는 다음과 같이 된다.

여기서 J₁은 제1종 1차 베셀 함수를 나타낸다. 식 (6)을 행렬식의 형태로 정리하면 식 (7)과 같이 나타낼 수 있다.

여기서 SM(k)는 NDIF법에서의 시스템 행렬이고, A는 기여도 벡터이다. 시스템 행렬 SM(k)의 i번째 행과 s번째 열 위치에 있는 성분은 다음의 식에 의해 주어진다.

또한, 기여도 벡터 A의 s번째 성분은 A_s로 주어진다.

음향 공동의 고유치는 식 (9)와 같이 시스템 행렬의 판별식이 0을 만족하는 파수(wavenumber) 값에 의해 구해진다.

하지만, 식 (9)에 의해 구해진 고유치들 속에는 음향 공동과 동일한 형상을 가지는 멤브레인의 고유치들에 해당하는 허위 고유치들(spurious eigenvalues)이 포함되어 있다⁽⁷⁾. 그래서, 허위 고유치를 제거하기 위하여 식 (10)과 같이, 식 (9)의 양변을 멤브레인 시스템 행렬 SM_mem(k)의 판별식 값으로 나눈 방법에 의해, 허위 고유치가 배제된 음향 공동의 고유치를 구할 수 있다⁽⁷⁾.

아울러, 음향 공동의 고유모드 형상은 위에서 구해진 고유치를 식 (7)의 시스템 행렬에 대입하여 얻은 고유벡터를 식 (3)에 대입하는 방법에 의해 그려질 수 있다.

3.1 원형 음향 공동의 고유치 추출 문제점

NDIF법에서 노드 수를 증가시킬 경우 고유치 추출 정확성(수렴성 포함)과 부가적인 문제점을 확인하기 위해, Fig. 2와 같은 엄밀해가 존재하는 원형 음향 공동을 고려하였다. Fig. 2에서 알 수 있듯이, 원형 음향 공동의 경계는 16개와 32개 노드로 각각 이산화되었다.

Fig. 2

Circular acoustic cavity of unit radius discretized with 16 and 32 boundary nodes, respectively

Fig. 2의 이산화 모델에 NDIF법을 적용하여, 식 (10)의 좌변에 해당하는 시스템 행렬의 판별식 값 det(SM_net(k))을 파수 k의 함수로 판별식 곡선을 그리면 Fig. 3과 같다. Fig. 3에서 점선은 Fig. 2 (a)의 16노드로 이산화된 경우이고, 실선은 Fig. 2 (b)의 32 노드로 이산화된 경우이다.

Fig. 3

Determinant curve for the circular acoustic cavity (dotted line: 16 nodes, solid line: 32 nodes)

16노드를 사용한 점선 판별식 곡선에서 총 8개의 고유치 C₁ ~ C₈이 구해졌으며, 이들 고유치들은 Table 1의 두번째 열에 요약되었다. Table 1에서 이들 고유치들 중에 일곱 번째 고유치(C₇)를 제외한 나머지 7개 고유치들은 엄밀해⁽¹³⁾와 정확히 일치하는 것을 확인할 수 있다. 이 일곱 번째 고유치를 엄밀해에 수렴시키기 위해 노드의 개수를 2배로 증가시켜 32개의 노드로 음향 공동을 이산화 시켰다. 이 경우에 대한 판별식 곡선은 Fig. 3의 실선에 해당된다. 실선 판별식 곡선에서 총 4개의 고차 고유치 C₅ ~ C₈만이 구해졌으며, 이들 고유치들은 이들 고유치들은 Table 1의 세번째 열에 요약되었다. 참고로, Fig. 3에서 M₁ ~ M₅는 음향 공동과 같은 형상을 가진 멤브레인의 고유치를 의미하며, 이들 고유치들은 허위고유치(spurious eigenvalue)로 불리어진다.

Table 1

Eigenvalues of the circular acoustic cavity by NDIF method, the exact solution(13), and ANSYS (parenthesized values denote errors (%) with respect to the exact solution)

Fig. 3의 실선 판별식 곡선의 저주파수 영역에서 곡선이 불규칙적으로 발산하는 원인은 다음 절에서 규명되어 진다.

Table 1에서 16개 노드를 사용한 경우 일곱 번째 고유치 값은 6.4160으로 엄밀해 6.4156과 0.0062 %의 오차를 가지나, 32개의 노드를 사용한 경우 일곱 번째 고유치 값은 엄밀해에 정확히 수렴했음을 알 수 있다. 위와 같이, NDIF법은 노드 수를 증가시키면 고유치가 엄밀해에 수렴하는 장점을 가지는 것으로 확인되나, 반면에 저차의 고유치가 구해지지 않는 단점도 가지고 있음을 확인할 수 있다. 이 논문에서는 이러한 단점을 극복하기 위한 새로운 기법을 제안하고자 한다.

한편, Table 1에서 2840개의 많은 노드를 사용한 ANSYS에 의해 구해진 고유치들을 엄밀해와 비교해보면, NDIF법이 적은 개수의 노드를 사용했음에도 불구하고 NDIF법의 정확성이 ANSYS에 비해 우수함을 확인할 수 있다.

3.2 저주파수 영역 발산 문제점 해결 방안

일반적으로, 어떤 정사각 행렬 M의 판별식은 행렬 M의 고유치들의 곱으로 식 (11)과 같이 구해질 수 있다⁽¹⁴⁾.

여기서 S는 행렬 M의 크기, λ_i는 행렬 M의 i번째 고유치이며 다음의 대소 관계를 만족한다.

여기서 R은 행렬 M의 랭크(rank)이며, 행렬 M이 불량 조건 행렬(ill-conditioned matrix)이 아니라면 R=S를 만족한다. 만약 행렬 M이 불량 조건 행렬이라면, λ_R보다 작은 고유치 λ_R+1, λ_R+2, ..., λ_S은 무의미한 값을 가지게 된다. Fig. 3의 실선 판별식 곡선이 저주파수 영역에서 발산하는 이유는, 저주파수 영역에서 시스템 행렬이 불량 조건 행렬이 되어 무의미한 고유치들(λ_R+1, λ_R+2, ..., λ_S)이 시스템 행렬의 판별식 계산 때 곱해졌기 때문인 것으로 추정된다. 이 추정의 타당성을 확인하기 위해, Fig. 2 (b)의 음향 공동에 대한 시스템 행렬의 랭크 곡선을 Fig. 4의 실선과 같이 구하였다. 이 랭크 곡선을 살펴보면, 파수가 약 4.9 이하인 저주파수 영역에서는 랭크 값이 32보다 작아서(즉, R < S), 시스템 행렬이 불량 조건 행렬이 되었음을 확인할 수 있다. 참고로, Fig. 4에서 점선은 음향 공동과 같은 형상을 가진 멤브레인의 시스템 행렬 랭크 곡선이다.

Fig. 4

Rc and Rm for the acoustic cavity and membrane discretized using 32 nodes

이 논문에서는 판별식 곡선이 발산하는 문제점을 해결하기 위해, 위에서 설명한 무의미한 고유치들을 판별식 계산 때 제외시키는 방안을 사용하고자 한다. 이를 위해, 판별식 계산식 식 (10)의 좌변은 식 (11)에 의해 다음과 같이 우선 변형한다.

여기서 $$ {\lambda }_i^{\left(c\right)}$$와 $$ {\lambda }_i^{\left(m\right)}$$은 각각 SM와 SM_(mem)의 i번째 고유치를 의미한다. 상술한 바와 같이, 저주파수 영역에서 시스템 행렬의 판별식 값이 발산하는 이유는, Fig. 4에서 랭크 값이 32보다 작은 파수 영역(대략 k < 4.9)에서 무의미한 고유치들이 식 (13)의 계산에서 사용되기 때문이다. 이러한 문제점을 해결하기 위해, R < S인 저주파수 영역(k < 4.9)에서는 식 (13)을 계산할 때 곱해지는 고유치들의 개수가 최소 랭크 값과 같도록 식 (13)을 식 (14)와 같이 변형한다.

여기서 min(R_c)와 min(R_m)는 각각 SM와 SM_(mem)의 랭크의 최소값을 의미하며, k_b는 SM이 더 이상 발산하지 않는 파수의 최소값을 의미한다. 참고로, 32개 노드로 이산화된 원형 음향 공동의 경우, min(R_c) = 19, min(R_m) = 19, k_b = 4.9이다.

식 (14)를 이용하여 32개 노드로 이산화된 원형 음향 공동(Fig. 2 (b))에 대한 판별식 곡선을 그리면 Fig. 5와 같다. 이 판별식 곡선과 Fig. 3의 실선 판별식 곡선을 비교해보면, 식 (14)를 사용하는 새로운 기법이 저주파수 영역에서 판별식 곡선이 발산하는 문제점을 성공적으로 해결해 줌을 확인할 수 있다. Fig. 5로부터 추출된 고유치 C₁ ~ C₈은 Table 1의 엄밀해와 정확히 일치함이 확인되었으며, 이로부터 제안된 기법은 해의 정확성 측면에서도 우수하다고 말할 수 있다.

Fig. 5

Determinant curve for the circular acoustic cavity discretized using 32 nodes, plotted by Eq. (14)

Fig. 5에서 추출한 원형 음향 공동의 j번째 고유치 C_j를 SM(k)에 대입한 SM(C_j)로부터 j번째 고유벡터 V^(j)를 식 (15)와 같이 구할 수 있다.

이 고유벡터의 성분 $$ {\mathrm{V}}_1^{\left(j\right)}, {\mathrm{V}}_2^{\left(j\right)},\cdots ,{\mathrm{V}}_{32}^{\left(j\right)}$$을 식 (3)의 기여도 계수에 대입하면, 식 (16)과 같이 j번째 고유모드 형상 식이 구해진다.

식 (16)으로부터 1차부터 8차까지 8개의 고유모드 형상이 Fig. 6과 같이 성공적으로 추출되었으며, 이들 형상들은 2840개의 많은 노드를 사용한 Fig. 7의 ANSYS 결과와 정확하게 일치함을 확인할 수 있다.

Fig. 6

Mode shapes obtained by the proposed method using 32 nodes

Fig. 7

Mode shapes obtained by ANSYS using 2840 nodes