I. 서 론

음향 방사를 모사하는 Kirchhoff 적분식을 풀 때 시간에 대해 순차적으로 응답을 계산하는 방식인 시간 진행 기법 (time marching scheme)을 사용하는 시간 영역 음향 경계요소법 (Transient acoustic Boundary ElementMethod, TBEM)은 다양한 음향 문제를 시간 영역에서 다룰 수 있는 방법이다 ^[1-4]. TBEM 방법은 중/저 주파수를 포함하는 과도적 음향 문제들을 다룰 수 있는 장점이 있지만, 해가 지수함수적으로 발산하는 수치적 불안정성으로 인해 널리 사용되지 못하여 왔다. 특히 외부 방사 문제를 다루는 경우, 비유일성 (non-uniqueness) 문제에 따라 내부 영역의 음향 모드의 영향을 받는 수치적 불안정성이 발생하는 문제가 있다 ^[5-7].

주파수 영역에서 Helmholtz 적분식을 이용할 때도 마찬가지이지만, Kirchhoff 적분식을 이용하여 외부 방사 문제를 계산할 때에는 대상 방사체의 내부 영역에서의 공명 주파수와 그 근방에서 정확한 해를 구할 수 없는 비유일성 문제가 발생한다. 이는 외부 문제와 Dirichlet 경계 조건이 주어진 내부 문제가 동일한 고유 방정식을 소유하기 때문에 발생하며, Kirchhoff-Helmholtz 적분식의 일반화된 형태인 Fredholm 적분식의 고유한 수학적인 문제이다 ^[8]. 따라서 TBEM 계산에서는 가상적인 내부 공명에 기인한 비유일성 문제와 수치적 불안정성 문제를 동시에 해결하여야 한다. 주파수 영역 경계요소법의 경우에는 비유일성 문제의 원인과 해결 방안에 대해 이미 많은 연구가 수행되었다. 특히, 경계 적분식과 그 공간 미분된 적분식을 선형 조합한 Burton-Miller 방법을 이용하는 경우에는 비유일성 문제를 해결할 수 있음이 수학적으로 증명된 바 있다 ^[9]. TBEM 방법에서도 Burton- Miller 방법을 적용하여 모든 가상적인 음향 모드들에 대한 영향을 제거함으로써 비유일성 문제를 해결할 수 있었고, 수치적 안정성을 크게 향상시킬 수 있었다 ^[10]. 그러나 이 방법에서는 공간 미분된 적분 방정식을 사용하기 때문에 수식이 복잡하여 추가적인 계산량이 발생하며, 수음점이 경계 표면 절점과 가까워짐에 따라 피 적분 함수가 거리 제곱의 역수에 비례하는 초강 특이성 (hyper-singularity) 문제를 해결해야 했고 ^[11], 수치 오차가 주 원인인 수치적 불안정성을 완벽히 제거할 수 없었다 ^[12]. 수식의 복잡함이나 초강 특이성 문제를 피하면서 비유일성 문제를 해결하는 간단한 방법으로서, 내부의 가상적인 수음점에서의 응답을 0으로 구속하는 조건을 이용한 CHIEF (Combined Helmholtz Integral Equation Formulation) 방법이 주파수 영역 해석에서 잘 사용되어 왔으나 ^[13], 이 방법이 시간 영역 해석에 적용된 예는 아직 없었다.

본 연구에서는 내부의 가상적인 수음점과 경계 표면의 절점들간의 최소 거리에 대한 음파의 지연 시간을 고려하여 미지수인 계산하고자 하는 현재 시간의 경계 표면 음장을 구속함으로써, 순차적인 계산을 하는 시간 영역 음향 경계요소법에 적합한 CHIEF 방법을 수식화하였다. 예제로서 반지름 방향으로 진동하는 구의 음향 방사 문제를 다루었으며, CHIEF 방법의 효과를 높이기 위해 유효 주파수 범위 내에서 발생하는 가상적인 저차 음향 모드들의 비절점 (anti-node)에 근접하도록 가상 수음점들을 내부에 위치시켰다 ^[14]. 단일 등비 행렬식을 이용한 고유 모드 해석 기법 ^[5]을 이용하여 가상적인 내부 모드들의 특성을 관찰함으로써, 비유일성 문제 해결 여부와 수치적 안정화의 가능성을 검토하고자 하였으며, 저차의 가상 음향 모드 성분들이 크게 감쇠되는 것을 관찰하였다.

II. CHIEF 방법이 적용된 TBEM 수식

2.1 이산화된 Kirchhoff 적분식

시간 영역 음향 경계요소법은 다음의 Kirchhoff 적분식을 기초로 한다 ^[3].

(1)

여기서, c(r₀)는 입체각 (solid angle),  p(r, t)는 위치 벡터 r과 시간 t에서의 음압, R_s는 경계 표면 위의 지점인 r_s와 수음점 r₀간의 거리 (R_s=|r₀–r_s|), t₀는 계산하고자 하는 현재 시간, t_ret는 거리 R_s에 대한 음향 전파의 지연 시간을 고려한 경계 표면 지점 r_s에서의 과거 시간 (t_ret=t₀–R_s/c)을 나타낸다. c₀는 균일한 매질 내 자유 음장에서의 음속, S는 방사체 표면, n은 경계 표면에서 해석 영역 쪽으로 향하는 단위 법선 벡터, (∂/∂n)는 법선 방향으로의 공간 미분을 나타낸다.

Kirchhoff 적분식을 수치적으로 계산하기 위해서는 이산화된 시간의 경계 절점에서의 응답을 이용하여 임의의 위치와 시간에서의 응답을 나타낼 수 있는 보간법과 경계 요소 상에서 정밀하고 효율적으로 계산할 수 있는 적분 방법이 필요하다. 본 연구에서는 2차의 삼각 요소와 7점의 Gaussian 구적법을 이용하여 경계 적분을 계산하였으며, 3차의 Lagrange 시간 형상 함수를 이용하여 시간 이산화 오차를 줄이고자 하였다. 이와 같은 수치적 처리 기법을 사용하고, 모든 경계 표면의 절점들이 수음점으로 고려된 경우, 다음과 같이 이산화된 Kirchhoff 적분식을 얻을 수 있다.

       (2)

위 식에서, C는 입체각을 나타내는 대각 행렬(diagonal matrix)이며, P_k와 [∂P/∂n]_k는 각각 시간 k∆t에서의 음압과 그의 법선 방향 공간 미분 벡터, α_i와β_i는(k–i)∆t에서의 경계 표면 변수들에 의한 쌍극 (dipole) 성분을 나타내는 상수 행렬, γ_i는 단극 (monopole) 성분을 나타내는 상수 행렬, W는 경계 모델 형상에 따른 최대 지연 시간 간격 수이다.

2.2 CHIEF 방법의 적용

이미 언급된 바와 같이 CHIEF 방법은 내부의 임의 수음점에서의 응답을 0으로 만드는 구속 조건을 지배적분식에 부가하여 이를 과결정 (over-determination)시킴으로써 비유일성 문제를 해결하는 방법이다. 본 TBEM 해석의 경우에는 아래의 구속 조건 식을 이용하였다.

(3)

여기서, r_C는 내부에 존재하는 임의의 CHIEF 점을 나타내며, m∆t는 CHIEF 점과 경계 절점간의 최소 거리에 대한 지연 시간이다. 이 식에서 (k+m)∆t에서의 음압을 계산하는 이유는 제한 조건 식에 미지수인 계산하고자 하는 현재 시간 k∆t에서의 응답을 구속시키기 위함이다. 앞선 절과 동일한 이산화 방법을 적용하면 식 (3)은 다음과 같이 변환된다.

        (4)

위 식에서, α_i^C, β_i^C는 CHIEF 점들에 대한 (k–i)∆t에서의 경계 표면 변수들의 쌍극 성분을 내포하는 상수 행렬, γ_i^C는 단극 성분에 대한 상수행렬이다. 식 (4)와 식 (2)를 결합시키면 다음과 같이 과결정된 수식을 얻을 수 있다:

(5)

식 (2)와 식 (5)의 시간 미분항은 효율성과 정확성 그리고 수치적 안정성 측면에서 다양한 방법으로 이산화할 수 있으며 ^[15,16], 본 연구에서는 3차의 중앙 차분 기법 (central difference scheme)을 이용하였다. 또한, 인과 관계를 성립시키기 위해서 k∆t, (k−1)∆t, (k−2)∆t에서의 시간 미분항들은 각각 1차 후방 차분 기법 (backward difference scheme), 1차와 2차 중앙 차분 기법을 사용하였다. 이에 따라 식 (5)는 다음과 같이 다시 표현된다.

                                          (6)

여기서, ψ_i와 φ_i는 재배열된 상수 행렬로서, CHIEF 제한 조건 식들로 인해 비정방 (non-square) 구조를 가진다.

본 연구에서는 외부 음향 문제 해석의 대부분을 이루는 과도적 음향 방사나 산란 문제를 다룬다. 이 경우 식 (6)을 이용하여 미지수인 경계 표면에서의 음압 벡터 P_k를 시간의 흐름에 따라 순차적으로 계산할 수 있다. 실제 계산에 있어서, 현재 시간에서의 음압 벡터 P_k는 최대 지연 시간 이전의 과거 시간에서의 음압 벡터들과 현재 및 과거 시간의 음압의 공간 미분 벡터들로부터 다음과 같이 계산된다.

                     (7)

여기서, (  )⁺는 Moore-Penrose 유사 역행렬 (pseudo- inverse)을 나타낸다 ^[17]. 현재 시간 k∆t에서의 계산이 완료되면, 최대 지연 시간 이전의 음압 벡터들과 이의 공간 미분 벡터들이 시간 (k+1)∆t를 기준으로 갱신되어 P_k+1가 계산되게 된다. 식 (6) 및 (7)의 순환 계산 구조는 다중 입력 다중 출력 (multi-input, multi- output; MIMO) 무한 임펄스 응답함수 (infinite impulse response; IIR) 구조의 이산화된 컨볼루션 (convolution) 과정을 표현하고 있다. 그림 1은 이러한 MIMO IIR 계산 구조를 도식화 해주고 있다. 임피던스 (impedance) 경계 조건을 포함한 혼합 경계 조건이 주어지는 경우, 경계 표면의 변수 벡터들과 그와 관련된 상수 행렬들은 경계 조건에 따라 재배열되어야 하며, 이 때 시간의 함수로 표현된 임피던스 데이터가 필요하다 ^[18]. 어떤 경우에도, 최종적인 TBEM 수식은 식 (7)과 마찬가지로 MIMO의 IIR 형식의 순환 연산 구조를 가진다.

CHIEF 방법의 제한 조건 식들로 인해 식 (6)의 상수 행렬들인 ψ_i와 φ_i가 비정방 구조를 갖는다. 그러나 최종 식 (7)에서 ψ₀⁺가 곱해짐에 따라 최종 상수 행렬들인 (ψ₀⁺ψ_i)와 (ψ₀⁺φ_i)는 다시 정방 구조의 형태를 취한다. 이와 같은 변환은 임피던스가 주어진 혼합 경계 조건의 TBEM 계산에 있어서도 마찬가지이다. 결국, 최종 수식의 상수 행렬들이 정방 행렬로 변환이 되므로 CHIEF 방법 적용 전과 후의 TBEM 계산 방식에는 차이가 없게 된다. 따라서, 이러한 간편한 수식 전개를 통해 저차의 내부 공명 모드들에 관련된 비유일성 문제를 효율적으로 해결할 수 있게 된다. 그러나 주파수 영역 계산과 마찬가지로, CHIEF 방법을 효과적으로 적용하기 위해서는 내부에 두는 가상 수음점들의 위치가 관련된 내부 음향 모드의 비절점에 근접해야 하며 ^[14], 빠른 시간 변화를 보이는 고주파수 대역의 음향 모드들에 관련된 비유일성 문제 해결에는 효과적이지 않다.

III. 고유 모드 해석의 통한 내부 공명 모드 특성 분석

3.1 고유 모드 해석

식 (7)은 계산의 편의를 위해 다음과 같은 표준형 (canonical form)으로 다시 표현할 수 있다 ^[5].

(8)

또는,

                                        (9)

여기서, 〈∂  P/∂ n〉 _k는 경계 조건으로부터 주어지며, 순차 (sequential) 음압 벡터〈P〉_k는 음압 벡터 P_k의 최대 지연 시간 동안의 시간 변화를 나타낸다. 행렬 M은 TBEM 계산의 특성을 나타내는 단일 반복 행렬 (single iterative matrix)이다. 고유 방정식 (eigen equation)의 형식으로 수식화하기 위해, 아래와 같이 출력 변수인 음압 벡터 P_k가 시간에 대한 지수함수 λ^k와 임의의 공간 벡터 u의 곱으로 주어진다고 가정한다.

                      (10a,b)

위 식의 λ는 단위 시간 간격에 대한 응답의 변화율을 표현하며, α와 ω는 각각 감쇠 계수 (attenuation constant)와 주파수를 나타내고, u는 응답의 방사 모드 형상이다. 위의 가정에 따라 순차 음압 벡터〈P〉_k도 동일한 방법으로 아래와 같이 표현할 수 있다.

(11a,b)

여기서, 위 첨자 T는 전치 (transpose) 행렬을 나타내기 위한 표현이다. 식 (11a)를 식 (9)에 대입하고 입력 변수〈∂  P/∂ n〉 _k를 0으로 가정하면, 다음과 같은 고유 방정식을 얻을 수 있다.

                                                           (12)

식 (12)의 고유 해는 TBEM 계산에서 의미 있는 해 (non-trivial solution)로서, 해당 음향 문제의 내부 공명 모드와 관련된다.

고유해는 TBEM 계산의 모드 변수들인 모드 형상, 감쇠율, 모드 주파수에 대한 물리적인 정보를 가지고 있다. 고유 벡터의 실수 및 허수부는 고유 모드 형상의 시간 영역에서의 변화를 나타내며, 이를 시간 영역 파동 벡터 (time-domain wave vector)라고 정의한다 ^[18]. 고유치의 크기는 해당 모드의 단위 시간당 변화율을 나타내므로, 감쇠율 D_i는 다음 식과 같이 계산할 수 있다.

                                         (13)

또, 고유치의 위상을 이용하면 모드 주파수 f_i를 다음과 같이 계산할 수 있다 ^[7].

                                         (14)

여기서, arg(λ_i)는 고유치 λ_i의 위상을 나타낸다.

3.2 비유일성 문제에 의한 수치적 불안정성

TBEM 계산의 수치적 불안정성은 고유치의 크기가 1보다 큰 특정 고유 모드에 의해 결정되며, 지수 함수적으로 발산하는 특성을 갖는다 ^[5]. TBEM 계산 응답은 그와 같은 파동 벡터 성분을 포함하고 있을 때불안정해지게 되며, 다음과 같이 수식화 할 수 있다.

                                            (15)

위 식의 q_i는 시간 영역 파동 벡터이다.

수치적 불안정성의 직접적인 원인은 수치오차이며, 이로 인해 TBEM 계산의 극점 (pole)인 특정 모드의 고유치의 크기가 단위 원 바깥에 위치하게 되고, 응답은 지수 함수적으로 발산하게 된다. 해당 모드의 선형 변화율을 나타내는 고유치의 크기는 이론적으로 1에 가까운 것이 일반적이기 때문에, 매우 작은 수치 오차에도 불안정해질 수 있다. 즉, 해당 모드의 감쇠가 거의 없거나 혹은 수치 오차가 증가한 경우 불안정성의 발생확률이 높아지게 된다. 외부 음향 문제를 다루는 경우, TBEM 계산의 극점은 비유일성 문제와 관련된 내부 공명 모드들에 의해 결정이 되며, 이들에 의해 수치적 불안정성이 발생하게 된다. 특히, 음향 산란이나 방사 문제의 경우, 가상적인 내부 모드의 물리적인 감쇠가 없기 때문에, 이론적인 고유치의 크기는 1이며 수치적 불안정성의 확률이 매우 높아지게 된다. 본 연구에서는 고유 모드 해석을 통해 CHIEF 방법의 적용 전과 후의 가상적인 내부 모드들의 거동을 관찰하고, 이를 통해 비유일성 해결과 수치적 안정화의 가능성을 검토하였다.

IV. 수치 예제

수치 예제로서 반지름 방향으로 진동하는 구에서의 음향 방사 문제를 해석하였다. 구의 반지름은 0.25 m이었으며, 주어진 경계 모델은 136개의 2차 삼각 요소와 274개의 절점으로 구성되었다. 요소의 최대 크기는 0.155 m이어서 λ/3 조건에 의한 유효 주파수 범위가 740 Hz 미만이었다. 시간 간격 크기는 0.08 ms으로 결정하였고, 그에 따라 다음 수식과 같이 최소 절점간 거리 ∆h_min에 대한 시간 간격 크기의 비인 무차원수 ν는 0.87이었다.

                                                       (16)

위 수식의 무차원수 ν는 시간 영역 유한 차분 계산의 수치적 안정성과 관련하여 잘 알려진 CFL (Courant- Friedrichs-Lewy) 조건에서 사용된 바 있다 ^[19]. 또, 모델의 형상에 의한 최대 지연 시간은 18･∆t이었다. 구의 표면 진동 속도는 모든 경계 절점에 대해 동일하게 옥타브 대역의 충격 응답 형태로 주어지며, 이는 아래 식과 같다.

                      (17)

여기서, h(ω)는 옥타브 대역 필터이며, 진동 속도 u₀는 0.001 m/s이다. 반지름 방향으로 조화 진동하는 구의 방사 음압에 대한 주파수 영역의 해석해가 존재하기 때문에 ^[20], 이를 역 Fourier 변환함으로써 옥타브 대역의 충격 가진에 의한 경계 표면의 방사 음압의 시간 이력에 대한 해석해는 다음과 같이 계산할 수 있다.

(18)

위 수식에서, k와 a는 각각 파수 (wave number)와 구의 반지름이다. 본 연구에서는 500 Hz의 옥타브 대역 충격 응답을 계산하였다.

음향 산란이나 방사 문제와 같이 Neumann 경계 조건이 주어지는 외부 음향 문제를 해석하는 경우에는 Dirichlet 경계 조건이 주어진 내부 음향 문제에 대한 얇은 물체의 간접적 경계 적분식 (thin-body indirect boundary integral)의 해가 실제 해와 공유되면서 비유일성 문제가 발생하게 된다 ^[8]. 얇은 물체라 함은 대상체의 외부 경계면과 내부 경계면의 속도는 같지만, 음압은 서로 구분되는 경우를 의미한다. 가상적인 내부 음향 모드에 대한 해석해는 구면 좌표계 (r, θ, ψ)에서의 파동 방정식의 해를 이용하여 다음과 같이 표현할 수 있다 ^[20].

(19)

여기서, j_n과 P_n^m은 각각 구면 Bessel 함수 (spherical Bessel function)와 Legendre 연관 다항 함수 (associated Legendre polynomials)를 나타내며, Legendre 연관 다항 함수는 m의 절대값이 n 이하일 경우에만 의미 있다고 정의된다. 가상적인 내부 음향 모드들은 Dirichlet 경계 조건이 주어진 문제의 의미 있는 해이기 때문에, 대상체의 내부 경계면에서의 음압이 0이 되는 j_n(ka)=0의 조건을 만족해야 한다. 한편 얇은 경계면의 내부와 외부 사이의 음압 변화 정리 (Jump theorem)에 따라 외부 경계면은 비절면을 형성하게 되며 ^[8], 3.1절의 고유 모드 해석을 통해 얻은 방사 모드는 이러한 외부 경계면의 음압 분포를 나타낸다.

표 1은 1.5 kHz 이내에 있는 저차의 가상적인 내부 음향 모드들에 대한 해석해를 요약한 것이며, 본 연구에서는 이들 모드의 비절점 주위에 총 77개의 CHIEF 수음점들을 위치시켰다. 그림 2는 이들 CHIEF 점들의 위치와 경계 요소 모델을 보여준다. n=0인 모드에 대응하여 |r|<0.050 m 이내에 18개의 CHIEF 점들을 고르게 분포시켰고, n=1과 n=2인 모드들에 대응하여서는 |r|=0.116 m와 |r|=0.145 m의 거리에 각각 20개와 39개의 CHIEF 점들을 위치시켰다.

본 시험 예제의 TBEM 계산 모델은 6302개의 고유 모드를 가지며, 이는 경계 절점의 총 개수 n_p와 최대 지연 시간 간격 수 (W+5)∆t의 곱으로 주어진다. 그림 3는 CHIEF 방법의 적용 전과 후에 있어서 고유 모드들의 고유치 및 모드 주파수를 서로 비교해서 보여준다. 이 그림에서는 CHIEF 방법의 적용 여부에 따른 변화를 뚜렷하게 보여주기 위해서 TBEM 응답에 대한 기여도를 나타내는 고유치의 크기가 0.9보다 큰 모드들을 도시하였다. 그림 3 (a)와 같은 CHIEF 방법이 적용되기 전의 경우에는 j_n(ka)=0의 조건을 만족하는 주파수에서 고유치의 크기가 거의 1인 내부 모드들이 관찰되었으며, 이들은 식 (17)를 통해 얻은 해석해와 거의 일치하였다. 또한 주파수가 높아짐에 따라 이산화에 따른 수치 감쇠로 인해 고유치의 크기가 점차 작아지는 경향을 보였다. 총 20개의 불안정 모드들은 1.3 ~ 2.4 kHz 범위에서 존재하였으며, 최대 고유치의 크기는 1.0001이고, 이에 따른 감쇠율은 –9.40 dB/s이었다.

그림 3 (b)에서 보듯이 CHIEF 방법을 적용한 경우, 저차의 가상적인 내부 음향 모드들에 대한 감쇠율이 크게 증가하였다. 특히, 표 1에서 정리한 10개의 저차의 가상 음향 모드들의 고유치 크기는 0.90336 ~ 0.97188의 범위에서 관찰되었으며, 이에 따른 감쇠율의 범위는 3096 ~ 11035 dB/s이었다. 이는 제안된 방법에 의해 저차 모드들에 의한 비유일성 문제를 해결할 수 있음을 분명히 보여준다. 한편 가상 모드들의 주파수는 해석해와 비교하여 최소 70 Hz 이상 저주파수 측으로 이동하였고, 총 64개의 불안정 고유 모드들이 경계 요소 크기에 의한 고주파수 한계인 740 Hz를 넘어선 1.7 ~ 2.6 kHz에서 관찰되었다. 최대 고유치의 크기는 1.0187이며 이에 해당하는 감쇠율은 –2008 dB/s이었다.

그림 4는 CHIEF 방법 적용 전과 후의 (1,0,0) 모드의 시간에 따른 모드 형상 변화를 나타내는 파동 벡터를 보여주며, 모드 형상이 왜곡되었음을 알 수 있다. 모드 형상의 변화 정도를 파악하기 위해 CHIEF 방법 적용 전과 후의 모드 형상 u_a와 u_b에 대한 MAC (Modal Assurance Criterion) 값을 아래 식과 같이 계산하였는데, MAC값은 모드 해석 과정에서 모드 형상의 상관 관계를 파악하기 위해 널리 쓰이는 인자이다 ^[21].

                             (20)

예로서, (1,0,0) 모드의 경우에는 MAC 값이 0.470이었다. 이와 같이 모드 형상이 왜곡됨에도 해당 모드는 해석해의 고유치 크기 및 모드 주파수와 근접한 범위에서 해석해의 모드 형상과 가장 높은 상관 관계를 보였고, 표 1에 보인 다른 저차 모드들에서도 이와 비슷한 경향이 관찰되었다.

비유일성 문제를 해결하는 또 다른 방법인 Burton- Miller 기법을 적용하는 경우 공간 미분된 경계 적분식에서 기인한 정적인 내부 모드와 관련하여 단조 변화하는 비요동 (non-oscillatory) 모드들이 단조 증가하는 불안정성을 초래하게 된다 ^[22]. 대부분의 불안정 모드들은 유효 주파수 범위를 넘어서는 영역에서 발생하는데 반해, 비요동 모드에 의한 불안정성은 유효 주파수 범위 내에서 발생하기 때문에 해의 정확도에 크게 영향을 주게 된다. CHIEF 방법의 적용 전과 후에 있어서의 고유모드들에서는 비요동 모드가 관찰되지 않았는데, 이는 경계면의 음장이 0으로 제한된 조건하에서 내부 음장이 일정하게 분포된 정적인 내부 모드가 존재할 수 없었기 때문이다.

그림 5에서는 500 Hz 옥타브 대역 충격 응답으로 가진된 구의 표면 음압에 대한 계산 결과를 해석해와 비교하였다. TBEM 계산의 정확도를 파악하기 위해서 해석해에 대한 상대 오차 (relative error norm)를 다음 식과 같이 계산하였다.

                 (21)

여기서, p_A와 p_T는 방사 음압에 대한 해석해와 TBEM 계산 결과를 나타낸다. CHIEF 방법이 적용되기 전의 대한 계산 결과인 그림 5 (a)에서는 (1,0,0)의 가상적인 내부 음향 모드의 성분으로 인한 비유일성 문제가 발생하였음을 알 수 있다. 또한 비록 불안정한 고유 모드들을 포함하고 있지만, 최대 고유치에 의한 감쇠율이 –9.4 dB/s로 매우 작기 때문에 수치적 불안정성은 관찰되지 않았다. 그림 5 (b)에서는 CHIEF 방법이 적용됨에 따라 (1,0,0) 모드 성분이 크게 감쇠되었고, 해석해와 TBEM 계산 결과가 거의 일치함을 볼 수 있다. 수치적 불안정성이 관찰되기 되기 시작한 0.05 s 이전에 대한 TBEM 계산 결과의 상대오차는 CHIEF 방법이 적용되기 전의 경우에는 14.8 %이었지만, 적용 후에 1.7 %로 크게 줄었다. 그러나 유효 주파수 범위를 넘어서 존재하는 불안정 고유 모드들의 개수와 최대 고유치의 크기가 증가하였기 때문에, 지수함수적인 발산이 관찰되었다.

이와 같이, CHIEF 방법을 사용함으로써 가상의 내부 저차 음향 모드들에 기인한 TBEM 계산의 비유일성 문제를 해결할 수 있었고, Burton-Miller 방법을 사용함에 따른 비요동 모드에 의한 불안정성을 피할 수 있음을 관찰하였다. 그러나, 유효 주파수보다 높은 고주파수 대역의 고차 내부 음향 모드들에 의한 비유일성 문제는 여전히 존재하였고, 이들에 의한 수치적 불안정성이 증가하였다.

V. 결 론

본 연구에서는 비유일성 문제를 해결하기 위해 쉽게 사용할 수 있는 CHIEF 방법을 TBEM 해석에 적합하게 수식화하였고, 비유일성 문제 해결 여부 및 수치적 안정성을 검토하였다. 본 연구에서 수식화된 CHIEF 방법을 통해 저차의 내부 음향 모드들에 기인한 비유일성 문제를 해결할 수 있었으며, Burton- Miller 기법의 고질인 비요동 모드에 의한 수치적 불안정성 문제를 피할 수 있었다. 그러나 유효 주파수를 넘어선 고주파수 대역에서의 고차 내부 음향 모드들에 관련된 비유일성 문제가 남았으며, 수치적 불안정성은 CHIEF 방법의 적용 전과 비교하여 오히려 증가하였다. 따라서, 본 연구에서 수식화된 CHIEF 방법과 함께, 모든 고유 모드에 인위적인 감쇠를 적용하는 방법 ^[12]이나 불안정한 모드들만을 필터링하는 방법 ^[18] 등을 이용하여 고주파수 대역의 불안정성을 해결할 수 있다면, TBEM 계산에 있어서의 비유일성 문제를 완벽히 해결할 수 있으리라 생각된다.