I. 서 론
건축, 자동차, HVAC, 덕트 등 여러 분야에서 소음 절감을 위해 흡음재를 사용한다. 일일이 흡음재를 제작해서 성능을 테스트 해 보는 방법도 있겠지만 많은 시간과 비용이 소모될 것이다. 따라서 흡음재의 흡음 성능을 미리 예측하는 일은 산업적으로 중요한 일이다.
흐름 저항으로 흡음재의 특성 임피던스를 예측하는 Delany 와 Bazley[1] 모델 등이 알려져 있다. Selamet et al.[2] 등은 흐름 저항으로 구한 특성 임피던스와 전파 상수를 이용해 섬유계 흡음재가 2개 장착된 원통형 덕트의 투과손실을 구했다.
흐름저항으로 특성임피던스를 예측하는 모델은 글라스울 등의 섬유계 물질의 흡음 성능의 예측에 있어서는 꽤 양호한 성능을 보여주지만 골격의 운동이 소음 감쇠에 영향을 미치는 폼형 물질의 흡음 성능을 예측하는 데는 적절치 않다. 따라서 이런 경우는 Biot 모델을 적용하는 것이 적합하다.
Biot[3,4]는 토양에 관한 연구 중 완전 포화된 다공성 탄성체 내의 유체와 고체의 상호작용을 고려한 파의 전파 이론을 제시했다. Biot 이론에서 요구되는 여러 계수들 가운데 유효 밀도와 유효 체적 탄성 계수를 측정 가능한 변수들로 표현할 수 있도록 해주는 Johnson-Champoux-Allard-Lafarge(JCAL) 모델이 잘 알려져 있다.[5] Nennig et al.[6]는 Biot 모델과 JCAL 모델을 이용해서 한 가지 흡음재가 채워진 원통형 덕트의 투과손실을 구하는 모드 매칭법을 제시했다.
본 연구는 한 개의 흡음재만 다룬 Reference [6]의 모드 매칭법을 다층 흡음재로 확장시켰다. Biot의 이론과 JCAL 모델을 통해 원통형 다층 소음기의 해석에 필요한 연립방정식들을 얻어내었다. 수치적인 테스트를 통해 2층 소음기에서 처음 12개의 모드만으로도 꽤 정확한 투과손실을 얻을 수 있음을 보였다. COMSOL Multiphysics 소프트웨어를 통해 투과손실을 다시 계산했고 이를 본 연구에서 제시한 모드 매칭법의 결과와 비교함으로써 모드 매칭법의 유효성을 입증했다.
II. 지배방정식
Fig. 1과 같이 길이 L의 원통형 덕트에 N개의 다공성 탄성 재료들을 층층이 배열한 N층 소음기를 고려한다. 중앙부에는 평균 유량이 0인 공기가 있다. 중앙부를 기준으로 흡음재 1까지의 거리를 , 흡음재 2까지의 거리를 흡음재 N까지의 거리를 라고 하고, 소음기의 벽까지의 거리를 라고 한다.
공기가 흐르는 중앙부()에서, 공기의 음속을 , 파수를 , 각주파수를 𝜔 라고 하고, z방향 파수를 라고 하자. 음압과 변위는 파동방정식 (1)과 오일러 방정식 (2)을 만족시킨다.
를 만족하는 변위퍼텐셜 를 Eq. (2)에 대입하면 음압이 다음과 같이 표현된다.
여기서 은 공기의 r방향 파수를 의미한다. 이후 나오는 식들에 항은 생략한다.
다공성 탄성체 내부() 에는 팽창과 수축을 반복하면서 전파되는 압력파와 회전운동을 통해 전파되는 전단파가 있다. Biot 모델에서, 다공성 탄성체 내에는 고체 전파음과 유체 전파음이 있으며 각각 두 개의 압력파와 한 개의 전단파의 합으로 모델링 가능하다. 헬름홀츠 분해를 통해 고체의 변위 u와, 유체의 평균 변위 U는 다음과 같이 변위 퍼텐셜에 관한 식으로 표현할 수 있다.
여기서 는 유체 상태와 고체 상태의 압력파들 사이의 진폭의 비율이며 는 전단파들 사이의 진폭 비율이다. 의 정확한 표현식은 이 절의 마지막에서 설명하겠다.
또한 Eqs. (6)과 (7)에서 과 는 두 개의 압력파의 변위 퍼텐셜을, 𝜓는 전단파의 변위 퍼텐셜을 의미하며 아래와 같이 표현된다.[7]
는 원통좌표계에서 𝜓의 성분을 나타내며 원통좌표계에서 헬름홀츠 방정식을 풀면 베셀방정식 형태를 띄는 것을 생각해보면, Eqs. (8), (9), (10), (11), (12)에서 변위 퍼텐셜들이 베셀 함수 형태를 가지는 것은 자연스럽다. 우리는 𝜃 방향에 의존하지 않는 축대칭인 m = 0 번 모드만을 사용할 것 이다. 첨자 n은 소음기에서 n번째 층의 흡음재를 의미한다. Eq. (13)은 Eq. (5)와 사실상 같은 식이다.
고체 상태와 유체 상태의 응력-변형률 관계는 다음과 같다.
전체 응력 (total stress)과 공극압 (pore pressure) 은 다음과 같이 표현된다.
Eqs. (14), (15)의 P, Q, R, N은 Biot의 Gedanken 실험을 통해 유추해 볼 수 있는 계수들이다. N 은 흡음재의 골격의 복소 전단 탄성 계수이다. Reference [5]에 따르면, P, Q, R 은 골격을 이루는 물질이 비압축성이라는 가정하에 다음과 같이 표현됨이 알려져 있다.
여기서 프레임의 체적 탄성 계수 는 Lame 상수와 푸아송 비 𝜈 사이의 관계식에서 다음과 같다.
다공성 탄성체 내부의 공기는 프레임과 상호작용하기 때문에 공기의 체적 탄성 계수 와 밀도 는 진동수에 따라 복잡하게 변하게 된다. 본 연구에서는 JCAL 모델을 적용하여 와 를 표현했다.[5]
여기서 는 굴곡도, 𝜙는 공극률, 𝜎는 흐름 저항, 는 공기의 밀도, 𝜂는 공기의 점성, 𝜅는 공기의 열전도율, 는 정압 비열, 𝛬는 점성 특성 길이, 은 열적 특성 길이, 𝛾는 비열비, 는 대기압이다. 은 정적 열 투과율이며 간소화된 Lafarge 모델에서 다음과 같다.[5]
Biot 이론에서, 다공성 탄성체 내에는 두 개의 압력파(첨자 1, 2)와 한 개의 전단파(첨자 3)가 존재하며 각각의 파수는 다음과 같다.
는 프레임을 이루는 물질의 밀도이며, 는 공극률을 고려한 프레임의 겉보기 밀도이다.
지금까지 설명한 표현식들을 이용하면, Eq. (7)의 를 아래와 같은 식을 이용해 구할 수 있다.
III. 경계조건
공기와 다공성 탄성체 사이 () 의 경계조건은 다음과 같다.[6]
첨자 1은 첫 번째 층의 흡음재이고, 는 법선 방향 벡터다. Eq. (36a)은 응력과 압력의 평형을, Eq. (36b)는 질량 플럭스의 연속성을, Eq. (36c)는 압력의 연속성을 나타낸다.
다공성 탄성체와 다공성 탄성체 사이의 경계면에서 경계조건은 다음과 같다.[8]
Eq. (37a)은 응력의 연속성을, Eq. (37b)는 고체상태의 변위의 연속성을, Eq. (37c)은 질량 플럭스의 연속성을, Eq. (37d)는 공극압의 연속성을 나타낸다.
II 장에서 설명한 Biot 이론의 표현식들을 이용하면, 경계조건 Eqs. (36a), (36b), (36c), (37a), (37b), (37c), (37d)는 총 6N+1 개의 방정식으로 표현된다. N은 흡음재의 개수이다. 이 경계조건들을 행렬 형태로 표현하면 다음과 같다.
C는 크기의 정사각행렬이며, 행렬의 각각의 성분들은 부록에 명시하겠다.
IV. 근을 찾는 방법
Eq. (39)가 비자명 해를 가지려면 을 만족시키는 를 찾아야 한다. 수치적으로 근을 찾을 때 많이 쓰이는 Newton 법을 사용해도 충분하며 본 논문에서는 Muller 법을 사용했다.
본 연구에서는 조건수가 매우 큰 행렬을 다루기 때문에 반올림 오차를 주의 깊게 다뤄야 한다. 층의 개수가 늘어날수록 이런 오차는 더욱 커진다. 모드 매칭법을 사용하려면 반올림 오차를 최소화하는 과정이 반드시 필요하다.
수치적으로 근을 찾는 방법인 Newton법과 할선법, Muller 법 등은 모두 초기값을 필요로 하며 초기값에 따라 어떤 근을 찾게 될지 예측이 힘들기 때문에 빠지는 근이 없도록 충분히 많은 초기값들을 반복해서 대입해야 한다.
Reference [6]에서, Nennig는 편각원리를 이용해 빠지는 근들이 없도록 근을 찾는 방법을 제안했다. 그러나 이 방법은 Newton 법 보다 시간이 더 오래 걸리는 단점이 있으며 반올림 오차가 커지면 잘못된 결과를 주는 경우도 생겨서 본 연구에서는 사용하지 않았다.
V. 모드 매칭
기본적으로 Reference [6]의 표기법을 따랐다. 소음기 내의 공기층에서 압력 p는 모드를 이용하여 다음과 같이 근사적으로 표현 가능하다.
여기서 I, II, III 이다(Fig. 1 참고). K는 사용한 모드의 개수이고, 는 영역에서 찾은 K개의 (모드)들 중에서 허수부의 크기가 m 번째로 작은 원소이다. 위에 붙은 + 와 – 부호는 파의 진행 방향이다. 흡음재가 채워진 II 영역에서 를 구하는 과정은 지금까지 설명했다. 공기 통로 I 또는 III 영역에서는 wall condition 을 통해 을 구한 후 Eqs. (5), (41)을 이용하면 된다.
Eq. (40)에서 는 압력의 m번째 고유함수이며 Eq. (3)과 같다. 는 m번째 고유함수 앞에 붙은 계수이다.
z = 0 일 때( I 과 II 의 경계) Eq. (40)를 간단히 표현하면 아래와 같다.
여기서 과 는 각각의 모드를 원소로 가지는 길이가 K 인 벡터이며 다음과 같이 정의된다.
z = L 의 경우에는( II 와 III의 경계) Eq. (40)은 다음과 같다.
흡입구[Eqs. (42), (43)]와 배출구[Eqs. (46), (47)]에서의 압력 표현식 차이는 단지 길이 L만큼 파동이 이동하면서 생긴 위상차만큼의 대각행렬 를 곱해주면 되므로 여기서는 z = 0 (I, II) 기준으로 설명하고 z = L 부분에 대해서는 같은 과정을 반복해서 기술하지 않겠다.
Eq. (42)와 유사하게, 흡입구에서 공기의 변위 w 그리고 흡음재 내의 유체의 변위 U와 고체의 변위 u 에 대해서도 다음 형태로 나타낼 수 있다.
여기서도 Eq. (44)과 유사하게 다음과 같은 정의를 사용한다.
z = 0에서의 영역 I과 영역 II 사이의 경계조건은 다음과 같다.
Eq. (58a)는 압력의 연속성을, Eq. (58b)는 공기의 변위의 연속성을, Eqs. (58c), (58d), (58e)는 흡음재의 wall condition을 의미한다.
Eqs. (58a), (58b), (58c), (58d), (58e)의 경계조건을 모두 만족하려면 무한개의 모드가 필요하지만 지금은 K개의 모드만을 고려중이기 때문에 방정식을 풀기 위해서는 위의 경계조건을 적분한 형태로 바꿔서 쓰는 것이 관례이다. 적분 형태로 바꾼 경계조건은 다음과 같다.
여기서 는 가중 함수이고 m 은 모드 수(m = 1 ~ K)를 뜻한다. 이 가중 함수에 알맞은 함수를 선택하면 성능이 좋아지고 보통은 각각의 물리량에 맞는 고유 함수를 선택한다. 이 기법의 동기는 적분 계산을 수행할 때 직교성을 이용하여 관심있는 성분을 극대화 하는 것에 있다. 본 연구에서는 흡입구(z = 0) 에서 다음과 같은 가중 함수를 사용했다.
유사하게, 배출구(z = L)에서 사용한 가중 함수는 다음과 같다.
가중 함수의 선택에 관한 더 깊은 논의는 References [6], [9]를 참고하면 좋다.
흡입구(z = 0)에서 적분 형태의 경계조건 Eqs. (59a), (59b), (59c), (59d), (59e)를 간단히 표현하면 다음과 같다.
여기서 는 다음과 같이 정의된다. 첨자 *는 켤레 전치를 의미한다.
경계조건 Eq. (62a)는 z = 0에서 압력의 연속성을, Eq. (62b)는 공기의 변위의 연속성을, Eq. (62c)는 n 번째 흡음재의 wall condition을 의미한다.
경계조건 Eqs. (59c), (59d), (59e)를 각각 계산할 경우 반올림 오차가 발생해 성능이 떨어진다. 따라서 Nennig는 Reference [6]에서 Eq. (62c)처럼 식 3개를 전부 더해서 조건을 완화시켜 사용했으며, 본 연구에서도 그것을 따랐다.
이제 다음과 같은 표기법을 정의한다.
그러면 흡입구(z = 0)에서 경계조건 Eqs. (62a), (62b), (62c)는 다음과 같이 간단한 행렬식으로 표현된다.
앞서 설명했듯이, 배출구(z = L)에서는 흡입구랑 유사하지만 위상차만 곱해주면 되며 다음과 같이 표현된다.
여기서 는 다음과 같다.
Eqs. (73a), (73b)을 푸는 것이 본 연구의 목적이다. 시스템 Eqs. (73a), (73b)의 시스템을 푸는 방법은 일단 평면파 가정 과 배출구 부분에서 들어오는 파가 없다는 가정 으로 시작한다. 그 다음, 으로 두고 Eq. (73a)를 풀어서 를 구한 뒤 그렇게 구한 를 Eq. (73b)에 대입 후 를 구하고 그렇게 구한 를 이용해서 다시 Eq. (73a)를 푸는 과정을 수렴할 때 까지 반복하는 것이다.[10]
반복적으로 해를 찾아가는 이유는 행렬을 역산해서 시스템을 풀 때마다 발생하는 반올림 오차를 반복을 거듭하면서 줄일 수 있기 때문이다.
유의사항으로 가 정사각행렬이 아니기 때문에 역산할 때 유사 역행렬을 곱해야 한다. 또 변위와 압력의 크기를 비슷하게 맞추기 위해 모든 변위에 을 곱해주는 과정이 필요하다.
VI. 결과 및 해석
모드 매칭법의 유효성 검증을 위해 2층 소음기에 대해 수치적인 테스트를 진행했다. 공기의 밀도와 음속은 각각 1.212 kg/m3, 342.208 m/s이다. 테스트에 사용한 소음기의 규격은 Table 1, 흡음재의 물성치들은 Table 2와 같다. 흠음재의 이름과 물성치는 Reference [12]에서 가져왔다.
Table 2.
Sound absorber | 𝜙 |
𝜎 (kNs/m4) |
𝛬 (μm) | (μm) | (kg/m3) |
N (kPa) | 𝜈 | |
FM3 | 0.97 | 87 | 2.52 | 37 | 119 | 31 |
55 (1-0.055j) | 0.3 |
XFM | 0.98 | 13.5 | 1.7 | 80 | 160 | 30 |
200 (1-0.1j) | 0.35 |
RGW2 | 0.99 | 9 | 1 | 192 | 384 | 16.3 |
220 (1-0.1j) | 0 |
Fig. 2에서, 2층 소음기 A 의 두 층에 모두 같은 흡음재 XFM을 채운 경우에 대해 모드 매칭법을 적용해 투과손실을 구했다. 같은 흡음재를 2개 채운 것은 하나의 흡음재를 채운 것과 사실상 동치이므로 결과는 Fig. 2에 나타난 것처럼 단층 소음기의 투과손실과 일치했다. 테스트에는 XFM과 Silencer A를 사용하였으며 단층 소음기에 대한 투과손실 데이터는 Reference [6]에서 가져왔다.
다음으로 2-layered silencer A의 중심에서 가까운 순서대로 XFM과 FM3을 채운 상황을 고려했다. 편의를 위해 이후 내용에서는 이런 소음기를 XFM- FM3 silencer A 라고 칭한다.
Fig. 3에서, 모드의 개수에 따른 XFM-FM3 silencer A의 투과손실의 수렴 여부를 테스트했다. 그 결과, 모드 개수가 6개(+, – 방향을 합치면 12개)만 있어도 투과손실이 피크 구간을 제외하면 수렴값과 1 dB 내외로 같았다. 모드 개수가 12개(+, – 방향을 합치면 24개)가 되면 피크에서도 거의 수렴했다.
Fig. 4에서, XFM-FM3 silencer A에 대해 모드 매칭법과 유한요소법으로 각각 투과손실을 계산한 후 비교했다. 모드 매칭법에서는 25개(+, – 방향을 합치면 50개)의 모드를 사용했다. 유한요소법은 COMSOL Multiphysics 소프트웨어를 사용해서 수행했다. 반올림 오차의 영향 때문에 피크 부분에서 1 dB 내외로 오차가 약간 발생했지만 두 결과가 거의 일치한다고 볼 수 있다.
Fig. 5에서, Fig. 4의 FM3를 글라스울 RGW2로 바꾼 XFM-RGW2 silencer A 에 대해 모드 매칭법과 유한요소법으로 구한 각각의 투과손실을 비교했다. 그 결과, 글라스울 흡음재에 대해서도 모드 매칭법이 유한요소법의 결과와 일치함을 보여주었다. 이로써, 본 연구에서 제시한 모드 매칭법이 유효함을 입증할 수 있었다.
VII. 결 론
본 논문에 제시한 원통형 다층 소음기의 투과손실을 예측하는 모드 매칭법을 2층 소음기에 적용해 본 결과, 처음 6개(+, - 방향 12개)의 모드만 써도 피크 이외에선 수렴값과 1 dB 정도 차이밖에 없었다. 또 본 논문에서 제시한 모드 매칭법을 COMSOL Multiphysics 소프트웨어를 사용한 유한요소법의 결과와 비교한 결과, 폼과 글라스울 양쪽 모두 유한요소법과 모드 매칭법으로 구한 투과손실이 거의 차이가 없었다. 이를 통해, 본 연구에서 제시한 모드 매칭법의 유효성을 입증했다.