Fig. 1은 고유치 해석 대상 임의 형상 음향 공동을 보여준다. 일반적인 고유치 해석 이론에서 음향 공동의 내부 음압은 식 (1)과 같은 헬름홀츠 방정식(Helmholtz equation) 형태의 지배방정식을 만족한다⁽²⁰⁾.

Fig. 1 
Arbitrarily shaped acoustic cavity with rigid-wall boundary Γ, which is discretized with boundary nodes P₁,P₂,…,P_N

여기서 k는 파수(wavenumber)이며 각주파수 ω와 음속 c로 구성된 식 k-ω/c에 의해 구해진다. 그리고 p(r)은 Fig. 1에서 음향 공동 내부의 한 점 P에서의 음압(sound pressure)을 나타내며, r은 점 P에 대한 위치 벡터를 의미한다.

음향 공동의 경계 조건은 강체벽(rigid-wall) 경계조건과 개방(pressure-released) 경계 조건이 존재하며, 이 논문에서는 가장 일반적인 강체벽 경계 조건만을 고려한다. 강체벽 경계 조건은 식 (2)와 같이 1차 미분 방정식의 형태를 가진다.

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

NDIF법을 적용하기 위하여 Fig. 1과 같이 음향 공동의 경계 Γ는 N개의 노드들 P₁,P₂,…,P_N으로 이산화된다. 다음으로, s번째 경계 노드 P_s에 대한 무차원 동영향 함수를 다음과 같이 정의한다.

여기서 J₀(⋯)는 제1종 0차 베셀 함수를 나타내며, 벡터 r_s는 경계 Γ에 위치한 s번째 노드에 대한 위치 벡터를 의미한다. 참고로, 식 (3)의 무차원 동영향 함수는 물리적으로 노드 P_s에서 단위 음압이 ω의 주파수로 발생되었을 때 음향 공동의 내부 점 P에서의 음압의 크기를 의미한다. 참고로 식 (3)은 지배방정식 식(1)을 정확히 만족하며 음향 공동의 음압을 구하기 위한 기저 함수로 사용된다.

음향 공동 내부 음압 p(r)은 경계 노드에서 정의된 무차원 동영향 함수 식 (3)의 선형 결합으로 식 (4)와 같이 근사화된다.

여기서 A_s는 s번째 노드에서 정의된 무차원 동영향 함수가 내부 음압에 미치는 기여도를 나타내는 계수이다. 참고로, 식 (4) 또한 지배방정식 식 (1)을 만족한다.

음향 공동의 경계 Γ에 주어진 강체벽 경계 조건 식 (2)는 N개의 경계 노드 P₁,P₂,…,P_N에 대해 식 (5)와 같이 이산화된다.

여기서 r_i와 n_i는 Fig. 1과 같이 경계 노드 P_i에 대한 위치 벡터와 법선 벡터를 각각 의미한다.

근사화된 음압 식 (4)를 이산화된 경계조건 식 (5)에 대입하면 다음의 식을 얻을 수 있다.

식 (6)에서 법선 방향 n_j에 대한 미분을 수행하면, 식 (7)이 얻어진다.

여기서 J₁(⋯)은 제1종 1차 베셀 함수를 나타낸다.

식 (7)을 행렬식의 형태로 정리하면 식 (8)과 같이 된다.

여기서 SM(k)는 파수 k를 변수로 가지는 시스템 행렬이고, A는 기여도 벡터를 의미하며 벡터 A의 s번째 성분은 A_s이다. 시스템 행렬 SM(k)의 성분은 식 (9)에 의해 주어진다.

음향 공동의 고유치는 식 (8)에서 주어진 시스템 행렬 SM(k)의 판별식을 계산함으로써 추출될 수 있다. 음향 공동의 고유치들은 식 (10)을 만족하는 파수 값들과 일치한다.

하지만, 식 (10)에 의해 구해진 고유치들 속에는 음향 공동고유치 이외의 허위 고유치들(spurious eigenvalues)이 포함되어 있다. 이들 허위 고유치들은 음향 공동과 동일한 형상을 가진 멤브레인의 고유치들에 해당됨이 과거 연구^(2,3)에서 밝혀졌다. 허위 고유치가 발생하는 이유는 시스템 행렬 SM(k)이 멤브레인에 대한 시스템 행렬 SM_mem(k)와 식 (11)과 같이 종속되어 있기 때문이다.

여기서 행렬 SM_net(k)은 허위 고유치를 제외한 음향 공동 고유치만을 제공하는 시스템 행렬을 의미한다. 식 (11)의 양변에 대한 판별식 값을 구하면 다음과 같이 된다.

식 (12)에서 det(SM_net(k))=0을 만족하는 파수 값들이 음향 공동의 고유치와 정확히 일치하며, 식 (12)를 이용하면 det(SM_net(k))=0은 식 (13)으로 대체 가능하다.

식 (13)으로부터 허위고유치를 제외한 음향 공동의 고유치만을 구할 수 있다. 아울러, 음향 공동의 j번째 고유모드 형상은 j번째 고유치를 식 (8)에 대입하여 j번째 기여도 벡터를 구하고, 이 기여도 벡터의 성분값들과 j번째 고유치를 식 (4)에 대입하는 방법에 의해 구해진다.

NDIF법에서 해석 대상 음향 공동의 고유치를 추출하기 위해선 관심 주파수 범위를 설정해야 하는데, 이 논문에서는 관심 주파수 범위의 변화에 따라 추출되는 고유치의 정밀도 저하가 발생하는 문제점을 해결하고자 한다. 이러한 정밀도 저하 현상을 면밀히 분석하기 위해 Fig. 2와 같은 엄밀해⁽²¹⁾가 존재하는 원형 음향 공동이 고려된다. NDIF법을 적용하기 위하여 Fig. 2와 같이 음향 공동의 경계는 32개의 노드로 이산화된다. 식 (13)을 이용하여 원형 음향 공동에 대한 판별식 곡선을 그리면 Fig. 3의 실선과 같다.

Fig. 2 
Circular acoustic cavity of unit radius discretized with 32 boundary nodes

Fig. 3의 실선을 살펴보면, 저주파수 영역에서 판별식 곡선이 발산함을 확인할 수 있다. 이러한 발산 원인은 식 (14)⁽²²⁾를 이용하여 설명할 수 있다.

Fig. 3 
Determinant curve for the circular acoustic cavity (solid line: original NDIF method, dotted line: improved NDIF method)

여기서 λi(c)와 λi(m)은 각각 시스템 행렬 SM과 SM_(mem)의 i번째 고유치, S는 시스템 행렬의 크기를 의미한다. 식 (14)는 Fig. 3의 실선을 그리기 위해 NDIF법에서 사용되는 판별식이며, 맨 우측 식의 분자와 분모 값을 계산해서 판별식 곡선이 그려진다. 저주파수 영역에서 시스템 행렬 SM또는 SM_(mem)의 랭크가 풀 랭크(full rank)가 아니라면, λi(c)와 λi(m)중 일부가 무의미한 값을 가지게 되어 식 (14)가 발산하게 된다는 것을 최근 연구^(18,19)에서 검증하였다. 그리고 이러한 발산 현상을 방지하기 위해 식 (14)를 보완한 식 (15)가 최근 연구에서 제안되었다.

여기서 R_c(k)와 R_m(k)은 파수 k에 따라 변하는 SM와 SM_(mem)의 랭크 함수(rank function)를 각각 의미하며 Fig. 4에 제시되었다. 그리고 min(R_c(k)), min(R_m(k))은 관심 파수 구간 k_min≤k≤k_max에서 R_c와 R_m의 최소값을 각각 의미한다. 일반적으로, 파수와 랭크는 비례하므로, 관심 파수 구간에서 랭크의 최소값은 식 (16)과 식 (17)과 같이 최소 파수일때의 랭크 값이다.

Fig. 4 
Rank functions R_c(k) and R_m(k) of the acoustic cavity and membrane discretized using 32 nodes, respectively

예를 들어, 현재의 예제인 Fig. 3의 경우는 k_min=1.5이므로 Fig. 4로부터 식 (18)과 식 (19)와 같이

임을 알 수 있다. 또한 k_b는 SM이 풀 랭크 값을 가지기 시작하는 파수를 의미하며, Fig. 4로부터 k_b=4.8임을 알 수 있다. 식 (15)를 사용하여 판별식 곡선을 그리면 Fig. 3의 점선과 같이 저주파수 범위에서 판별식 값이 발산하는 현상이 사라지면서 저주파수 고유치들(C₁~C₄)이 엄밀해와 정확히 일치된 값으로 성공적으로 추출됨이 확인되었다.

그런데, 만약 관심 파수의 최소값 k_min을 변화시킨다면 식 (15) 속에 포함된 min(R_c(k))와 min(R_m(k))의 값도 변하게 되어서, 식 (15)가 k_min에 종속되는 문제점이 존재한다. Table 1은 k_min을 변화시켰을 때 식 (15)에 의해 구해진 저주파수 고유치들의 정밀도 저하 현상을 보여준다. k_min=1.5일 때는 구해진 고유치 4개 모두가 엄밀해와 정확히 일치하나, k_min이 작아 짐에 따라 고유치들의 정밀도 저하가 발생함을 확인할 수 있다. 이러한 정밀도 저하의 원인은, k_min이 작아지면 min(R_c(k))와 min(R_m(k))도 동시에 작아지기 때문에 λi(c)와 λi(m)값들 중 유의미한 값들이 식 (15)의 계산에서 제외되기 때문인 것으로 분석된다.

관심 파수의 최소값 k_min의 변화에 의한 고유치 정밀도 저하 문제점을 해결하기 위해, 이전 연구에서 제안된 식 (15)를 개선하고자 λi(c)와 λi(m)값들이 식 (15)를 계산할 때 제외되지 않도록 하기 위해, 이 연구에서는 k_min에 종속되지 않는 식 (20)을 제안한다.

여기서 랭크 함수 R_c(k)와 R_m(k)은 Fig. 4의 의해 주어진다. 참고로, 식 (20)의 추가적인 장점은 파수 값 k_b를 기준으로 양쪽 구간에 대해 두 개의 식으로 분리되는 식 (15)와 다르게 하나의 식으로 시스템 행렬의 판별식을 계산할 수 있다는 것이다.

새로이 제안된 식 (20)을 사용하여 구한 원형 음향 공동에 대한 판별식 곡선은 Fig. 5에 제시되었다. 이 판별식 곡선으로부터 구한 고유치들(C₁~C₈)은 곡선에 직접 표시하였다. 이들 고유치들은 엄밀해(21)와 정확히 일치하는 것으로 확인되며, 새로이 제안된 식 (20)은 타당하고 정확한 것으로 판단된다. 하지만, 이 판별식 곡선은 시스템 행렬이 풀 랭크가 아닌 저주파수 구간(즉, 0≤k≤4.8)에서 많은 불연속 지점들을 가진다. 참고로, 이러한 불연속 지점들은 Fig. 4의 두 랭크 곡선에서 불연속 지점들과 정확히 일치한다. 판별식 곡선에서의 이러한 불연속 현상은 고유치를 추출하는 데에 혼선을 줄 수 있기에, 멤브레인에 대한 이전 연구⁽¹⁷⁾에서도 개발된 불연속 지점 제거 방법을 식 (20)의 분자와 분모에 동시에 응용해서 적용한다. 그러면, 식 (20)은 식 (21)과 같이 된다.

Fig. 5 
Determinant curve for the circular acoustic cavity obtained by the proposed NDIF method using Eq. (20)

여기서 ∏λadd c(k)은 음향 공동 시스템 행렬(SM)의 랭크 변화에 의해 추가적으로 곱해지는 고유치들의 곱을 의미하며, ∏λadd m(k)은 멤브레인 시스템 행렬(SM_(mem))의 랭크 변화에 의해 추가적으로 곱해지는 고유치들의 곱을 의미한다.

식 (21)를 사용하여 구한 판별식 곡선은 Fig. 6과 같다. Fig. 5와 Fig. 6을 비교해보면, 새로 제안된 식 (21)에 의해 판별식 곡선의 불연속성이 성공적으로 제거 되어 졌음을 확인할 수 있으며, 구해진 고유치들(C₁~C₈)도 엄밀해⁽²¹⁾와 정확히 일치함이 확인되었다.

Fig. 6 
Determinant curve for the circular acoustic cavity obtained by the proposed NDIF method using Eq. (21)