Acoustic Wave Propagation Modeling using High-order Cell-based Finite-difference Methods

Jun Hyeon Jo; Wansoo Ha

doi:10.32390/ksmer.2021.58.4.281

Preview

Research Paper

Journal of the Korean Society of Mineral and Energy Resources Engineers. 31 August 2021. 281-289
https://doi.org/10.32390/ksmer.2021.58.4.281

Acoustic Wave Propagation Modeling using High-order Cell-based Finite-difference Methods

고차 셀 기반 유한 차분법을 이용한 음향파 파동 전파 모델링

Jun Hyeon Jo¹

Wansoo Ha²^*

조 준현¹

하 완수²^*

¹Postgraduate Student, Department of Energy Resources Engineering, Pukyong National University, Busan, Korea

²Associate Professor, Department of Energy Resources Engineering, Pukyong National University, Busan, Korea

¹부경대학교 에너지자원공학과 대학원생

²부경대학교 에너지자원공학과 부교수

^{*Corresponding Author}

ABSTRACT

We propose high-order finite-difference methods for acoustic wave propagation modeling considering the P-wave velocity and density. We can simulate wave propagation without directly differentiating discontinuous medium parameters using the cell-based finite-difference method. However, a limitation of the conventional cell-based finite-difference method using the second-order scheme is that it requires a large number of grids per wavelength to obtain an accurate result. We improve the conventional cell-based method to exploit arbitrary high-order schemes. In a numerical example, we compare results of the proposed methods with that of the analytic solution, and show that we can obtain accurate results with small number of grids per wavelength using the high-order cell-based methods.

Keywords

Wave propagation modeling

Acoustic waves

Cell-based finite-difference method

이 연구에서는 P파 속도와 밀도를 고려한 음향파 파동 전파 모델링을 위한 고차 셀 기반 유한 차분법을 제안하였다. 셀 기반 유한 차분법을 사용하면 불연속적인 매질 특성을 직접 미분하지 않고도 파동 전파를 모델링할 수 있다. 그러나 기존의 2차 차분식을 사용하는 셀 기반 유한 차분법에서는 정확한 결과를 얻기 위해 파장당 격자수가 커야 한다는 한계가 있었다. 이 연구에서는 임의의 고차 차분식을 사용할 수 있도록 기존 기법을 개선하였다. 수치 예제에서 해석해와의 비교를 통해 고차 셀 기반 유한 차분법을 이용하면 파장당 격자수가 작아도 정확한 결과를 얻을 수 있음을 보였다.

키워드

파동 전파 모델링

음향파

셀 기반 유한 차분법

MAIN

서 론
이 론
셀 기반 유한 차분법
고차 셀 기반 유한 차분법
2차원 및 3차원 확장
알고리듬의 효율적인 구현
차수에 따른 알고리듬의 정밀도 검증
분산 분석 및 안정성 분석
적용성 분석
결 론

서 론

파동 전파 모델링은 컴퓨터를 이용해 매질 내에서의 파동 전파를 구현해내는 기술이다. 탄성파 탐사시 현장에서 얻은 자료를 처리하기 위해 파동 전파 모델링을 사용한다. 실제 지하 매질은 탄성 매질이지만 자료 처리시 P파만 사용하기 위해 또는 컴퓨터 연산량 제한으로 인해 문제를 단순화시켜 탄성파 대신 음향파 모델링을 사용하는 경우도 많다(Schuster, 2017). 파동 전파 모델링에는 일반적으로 테일러 전개로부터 유도할 수 있는 유한 차분법을 가장 많이 이용하며(Alford et al., 1974), 유한 요소법(Marfurt, 1984)이나 스펙트럴 요소법(Komatitsch and Tromp, 1999)도 사용할 수 있다.

일반적인 유한 차분법 사용시 P파 속도 외에 밀도 변화를 고려할 경우 밀도에 대한 직접 미분항이 발생하게 된다. 연속적인 파동장과 달리 밀도와 같은 매질 특성은 지층 경계에서 불연속적일 수 있으므로 매질 특성의 공간에 대한 직접 미분은 피하는 것이 바람직하다(Ikelle and Amundsen, 2018). 매질 특성의 직접 미분을 피하는 방법의 하나로 엇격자 기법(staggered grid)을 사용할 수 있다(Virieux, 1986). 엇격자 기법은 응력과 입자의 속도로 표현된 1계 미분 파동 방정식에서 파동장과 매질의 특성을 다른 격자에 할당하는 기법으로, 변위로 파동장을 표현하는 2계 미분 파동 방정식을 사용할 경우에는 셀 기반 유한 차분법을 사용해 밀도의 직접 미분을 피할 수 있다(Min et al., 2004).

셀 기반 유한 차분법은 유한 차분법을 사용하여 격자에서 파동장을 계산하되, 모델 파라미터를 격자가 아닌 격자 사이의 공간인 셀에 할당하는 기법이다(Min et al., 2004; Min and Kim, 2006). 셀 기반 유한 차분법은 주파수 영역 탄성파 파동 방정식 사용시 자유 경계조건을 묘사하기 위해 처음 제안되었지만(Min et al., 2004) 시간 영역 모델링이나(Min and Kim, 2006; Lee et al., 2009) 음향파 모델링에도 사용할 수 있다(Lee et al., 2009). 그러나 기존의 셀 기반 유한 차분법은 공간 2계 미분시 2차 유한 차분식을 사용하므로 수치 분산을 피하기 위해 파장당 격자수가 20개 이상으로 커진다는 한계가 있다(Lee et al., 2009).

이 연구에서는 P파 속도와 밀도를 고려한 음향파 파동 방정식에 대해 고차 셀 기반 유한 차분식을 제안하였다. 균질 매질에서 일반적인 유한 차분식과 동일해지도록 수식을 유도하였으며, 수치 예제를 통해 고차 셀 기반 유한 차분식 사용시 2차 차분식을 이용한 방법에 비해 수치 분산이 감소하는 것을 보였다.

이 론

셀 기반 유한 차분법

이 논문에서는 간편한 수식을 이용하여 개념을 보다 정확히 설명하기 위해 먼저 1차원 파동 방정식을 이용해 기존 셀 기반 유한 차분법과 고차 셀 기반 유한 차분법을 설명한 후, 2차원과 3차원으로 확장하는 방법에 대해 설명하였다. P파 속도와 밀도를 고려한 1차원 음향파 파동 방정식은 다음과 같다.

(1)

\frac{1}{K} \frac{\partial^{2} P}{\partial t^{2}} = \frac{\partial}{\partial x} (\frac{1}{ρ} \frac{\partial P}{\partial x}) + f

위 식에서 P는 파동장, $ρ$ 는 밀도, $K = ρ V_{P}^{2}$ 는 체적탄성률, $V_{P}$ 는 P파 속도, $f$ 는 송신원이다.

셀 기반 유한 차분법에서 특정 격자 위치에서의 매질 특성은 인접 셀에 할당된 값의 평균값을 사용한다(Min et al., 2004). Fig. 1의 격자와 셀을 이용하여 셀 기반 유한 차분법을 적용할 경우 식 (1)의 좌변은 아래와 같이 표현할 수 있다.

(2)

\frac{1}{K} \frac{\partial^{2} P}{\partial t^{2}} \approx \frac{1}{2} (\frac{1}{K_{i - h}} + \frac{1}{K_{i + h}}) \frac{P_{i}^{n + 1} - 2 P_{i}^{n} + P_{i}^{n - 1}}{Δ t^{2}} = β_{i} \frac{P_{i}^{n + 1} - 2 P_{i}^{n} + P_{i}^{n - 1}}{Δ t^{2}}

위 식에서 체적탄성계수의 역수(압축률)는

(3)

β_{i} = \frac{1}{2} (\frac{1}{K_{i - h}} + \frac{1}{K_{i + h}})

이고, $h = 1 / 2$ 은 수식을 간편하게 만들기 위해 도입하였다. 정수는 격자 번호에 해당하고, $i + h$ 는 셀 번호에 해당한다(Fig. 1).

https://static.apub.kr/journalsite/sites/ksmer/2021-058-04/N0330580402/images/ksmer_58_04_02_F1.jpg

Fig. 1.

One-dimensional grids and cells.

한편, 식 (1)의 우변의 공간에 대한 미분은 다음과 같다.

(4)

\frac{\partial}{\partial x} (\frac{1}{ρ} \frac{\partial P}{\partial x}) \approx \frac{1}{Δ x} [\frac{1}{ρ_{i + h}} \frac{(P_{i + 1} - P_{i})}{Δ x} - \frac{1}{ρ_{i - h}} \frac{(P_{i} - P_{i - 1})}{Δ x}] = \frac{1}{Δ x^{2}} [ν_{i + h} P_{i + 1} + ν_{i - h} P_{i - 1} - (ν_{i + h} + ν_{i - h}) P_{i}]

위 식에서 밀도의 역수(비체적)는 다음과 같이 정의한다.

(5)

ν_{i + h} = \frac{1}{ρ_{i + h}} ν_{i - h} = \frac{1}{ρ_{i - h}}

식 (2)와 (4)를 $P_{i}^{n + 1}$ 에 대해 정리하면 다음과 같다.

(6)

P_{i}^{n + 1} = 2 P_{i}^{n} - P_{i}^{n - 1} + \frac{Δ t^{2}}{β_{i}} [\frac{1}{Δ x^{2}} (ν_{i + h} P_{i + 1} + ν_{i - h} P_{i - 1} - (ν_{i + h} + ν_{i - h}) P_{i}) + f_{i}^{n}]

이 식을 이용하면 시간 전진 기법으로 파동 전파 모델링을 수행할 수 있다.

고차 셀 기반 유한 차분법

공간에 대한 고차 셀 기반 유한 차분식은 공간 차분시 반격자에 대해 1계 미분의 고차 차분식을 적용하여 유도할 수 있다. 식 (4)의 공간 미분식에서 바깥쪽 미분에 임의의 $N = 2 M$ 차 중앙 차분식을 적용하고 안쪽 미분에는 반격자에 대한 2차 중앙 차분식을 적용하면 다음과 같다.

(7)

\frac{\partial}{\partial x} (\frac{1}{ρ} \frac{\partial P}{\partial x}) \approx \frac{1}{(\frac{Δ x}{2})} \sum_{\begin{matrix} m = - M \\ m \neq 0 \end{matrix}}^{M} C_{m}^{1, N} [v_{i + m h} \frac{P_{i + m} - P_{i}}{m Δ x}] = \frac{1}{Δ x^{2}} \sum_{\begin{matrix} m = - M \\ m \neq 0 \end{matrix}}^{M} \frac{2}{m} C_{m}^{1, N} v_{i + m h} (P_{i + m} - P_{i})

위 식에서 $m > 0$ 일 때

(8)

v_{i + m h} = \frac{1}{m} \sum_{p = 1}^{m} \frac{1}{ρ_{i + (p - h)}}

이고, $m < 0$ 일 때

(9)

v_{i + m h} = \frac{1}{| m |} \sum_{p = 1}^{| m |} \frac{1}{ρ_{i - (p - h)}}

이다. 1계 미분의 $N$ 차 유한 차분식 계수 $C_{m}^{1, N}$ 은 Table 1에 제시하였다.

Table 1.

Central finite-difference coefficients of the first derivative ( $C_{m}^{1, N}$ ), where $N$ is the order of a finite-difference scheme, and $m$ is the index of a finite-difference coefficient

Order ( $N$ )	$m$
Order ( $N$ )	-5	-4	-3	-2	-1	0	1	2	3	4	5
2					-1/2	0	1/2
4				1/12	-2/3	0	2/3	-1/12
6			-1/60	3/20	-3/4	0	3/4	-3/20	1/60
8		1/280	-4/105	1/5	-4/5	0	4/5	-1/5	4/105	-1/280
10	-1/1260	5/504	-5/84	5/21	-5/6	0	5/6	-5/21	5/84	-5/504	1/1260

식 (7)은 1계 미분 유한 차분 계수( $C_{m}^{1, N}$ )를 사용한 식으로, 수식을 간략하게 만들고 안정성 및 분산 분석시 기존의 유한 차분법 수식과 비교하기 위해 2계 미분 유한 차분 계수( $C_{m}^{2, N}$ )를 사용하도록 수식을 수정할 수 있다. Table 1의 1계 미분 유한 차분 계수와 Table 2의 2계 미분 유한 차분 계수 사이에 존재하는 다음 관계식을 이용하면

(10)

\frac{2}{m} C_{m}^{1, N} = C_{m}^{2, N} (m \neq 0)

식 (7)의 1계 미분 유한 차분 계수는 2계 미분 유한 차분 계수를 이용하여 다음과 같이 바꿔 쓸 수 있다.

(11)

\frac{\partial}{\partial x} (\frac{1}{ρ} \frac{\partial P}{\partial x}) \approx \frac{1}{Δ x^{2}} \sum_{\begin{matrix} m = - M \\ m \neq 0 \end{matrix}}^{M} C_{m}^{2, N} v_{i + m h} (P_{i + m} - P_{i})

위 식은 $N$ 차 셀 기반 유한 차분식으로 다음과 같이 $P_{i}$ 항을 분리하여 표현할 수도 있다.

(12)

\frac{\partial}{\partial x} (\frac{1}{ρ} \frac{\partial P}{\partial x}) \approx \frac{1}{Δ x^{2}} [\sum_{\begin{matrix} m = - M \\ m \neq 0 \end{matrix}}^{M} (C_{m}^{2, N}, v_{i + m h}, P_{i + m}) - (\sum_{\begin{matrix} m = - M \\ m \neq 0 \end{matrix}}^{M} C_{m}^{2, N} v_{i + m h}) P_{i}]

Table 2.

Central finite-difference coefficients of the second derivative ( $C_{m}^{2, N}$ )

Order ( $N$ )	$m$
Order ( $N$ )	-5	-4	-3	-2	-1	0	1	2	3	4	5
2					1	-2	1
4				-1/12	4/3	-5/2	4/3	-1/12
6			1/90	-3/20	3/2	-49/18	3/2	-3/20	1/90
8		-1/560	8/315	-1/5	8/5	-205/72	8/5	-1/5	8/315	-1/560
10	1/3150	-5/1008	5/126	-5/21	5/3	-5269/1800	5/3	-5/21	5/126	-5/1008	1/3150

2차원 및 3차원 확장

식 (12)를 2차원과 3차원으로 확장할 경우 매질 특성의 평균값 계산시 $y, z$ 방향 인접 셀을 고려해주면 된다. $y$ 방향 격자 변호를 $j, z$ 방향 격자 번호를 $k$ 라고 하면 2차원 수식에서 좌변의 평균 압축률은 다음과 같고,

(13)

β_{i, j} = \frac{1}{4} (\frac{1}{K_{i - h, j - h}} + \frac{1}{K_{i - h, j + h}} + \frac{1}{K_{i + h, j - h}} + \frac{1}{K_{i + h, j + h}})

우변의 $x$ 방향 공간 미분식과 평균 비체적은 Fig. 2의 격자상에서 다음과 같이 쓸 수 있다.

(14)

\frac{\partial}{\partial x} (\frac{1}{ρ} \frac{\partial P}{\partial x}) \approx \frac{1}{Δ x^{2}} [\sum_{\begin{matrix} m = - M \\ m \neq 0 \end{matrix}}^{M} (C_{m}^{2, N} v_{i + m h, j} P_{i + m, j}) - (\sum_{\begin{matrix} m = - M \\ m \neq 0 \end{matrix}}^{M} C_{m}^{2, N} v_{i + m h, j}) P_{i, j}] v_{i \pm | m | h, j} = \frac{1}{2 | m |} \sum_{p = 1}^{| m |} (\frac{1}{ρ_{i \pm (p - h), j - h}} + \frac{1}{ρ_{i \pm (p - h), j + h}})

$y$ 방향 미분도 위 식과 같은 방식으로 계산할 수 있다.

(15)

\frac{\partial}{\partial x} (\frac{1}{ρ} \frac{\partial P}{\partial y}) \approx \frac{1}{Δ y^{2}} [\sum_{\begin{matrix} m = - M \\ m \neq 0 \end{matrix}}^{M} (C_{m}^{2, N} v_{i, j + m h} P_{i, j + m}) - (\sum_{\begin{matrix} m = - M \\ m \neq 0 \end{matrix}}^{M} C_{m}^{2, N} v_{i, j + m h}) P_{i, j}] v_{i, j \pm | m | h} = \frac{1}{2 | m |} \sum_{p = 1}^{| m |} (\frac{1}{ρ_{i - h, j \pm (p - h)}} + \frac{1}{ρ_{i + h, j \pm (p - h)}})

3차원의 경우도 2차원과 같은 방식으로 좌변의 평균 압축률은

(16)

β_{i, j, k} = \frac{1}{8} (\frac{1}{K_{i - h, j - h, k - h}} + \frac{1}{K_{i - h, j - h, k + h}} + \frac{1}{K_{i - h, j + h, k - h}} + \frac{1}{K_{i - h, j + h, k + h}} + \frac{1}{K_{i + h, j - h, k - h}} + \frac{1}{K_{i + h, j - h, k + h}} + \frac{1}{K_{i + h, j + h, k - h}} + \frac{1}{K_{i + h, j + h, k + h}})

우변의 방향 공간 미분과 평균 비체적은

(17)

\frac{\partial}{\partial x} (\frac{1}{ρ} \frac{\partial P}{\partial x}) \approx \frac{1}{Δ x^{2}} [\sum_{\begin{matrix} m = - M \\ m \neq 0 \end{matrix}}^{M} (C_{m}^{2, N} v_{i + m h, j, k} P_{i + m, j, k}) - (\sum_{\begin{matrix} m = - M \\ m \neq 0 \end{matrix}}^{M} C_{m}^{2, N} v_{i + m h, j, k}) P_{i, j, k}] v_{i \pm | m | h, j} = \frac{1}{4 | m |} \sum_{p = 1}^{| m |} (\frac{1}{ρ_{i \pm (p - h), j - h, k - h}} + \frac{1}{ρ_{i \pm (p - h), j - h, k + h}} + \frac{1}{ρ_{i \pm (p - h), j + h, k - h}} + \frac{1}{ρ_{i \pm (p - h), j + h, k + h}})

와 같이 계산할 수 있고, $y$ 방향, $z$ 방향 미분과 평균 비체적도 같은 방식으로 계산할 수 있다.

https://static.apub.kr/journalsite/sites/ksmer/2021-058-04/N0330580402/images/ksmer_58_04_02_F2.jpg

Fig. 2.

Two-dimensional grids and cells for differentiation along the x-axis.

알고리듬의 효율적인 구현

식 (6)에서 시간 전진 기법 사용시 평균 압축률은 시간에 따라 변화하지 않으므로 처음 한 번만 계산한 후 평균값을 저장하여 사용하면 모델링 계산 시간을 줄일 수 있다(Park and Ha, 2019).

식 (12)에서 비체적 $v_{i + m h}$ 계산시 8차 유한 차분식을 사용할 경우

(18)

v_{i + h} = \frac{1}{ρ_{i + h}} v_{i + 2 h} = \frac{1}{2} (\frac{1}{ρ_{i + h}} + \frac{1}{ρ_{i + 3 h}}) v_{i + 3 h} = \frac{1}{3} (\frac{1}{ρ_{i + h}} + \frac{1}{ρ_{i + 3 h}} + \frac{1}{ρ_{i + 5 h}}) v_{i + 4 h} = \frac{1}{4} (\frac{1}{ρ_{i + h}} + \frac{1}{ρ_{i + 3 h}} + \frac{1}{ρ_{i + 5 h}} + \frac{1}{ρ_{i + 7 h}})

와 같은데, 고차 차분식 계산을 위해 동일한 연산을 반복하게 되므로 비효율적이다. 실제 계산시에는 다음과 같이 점화식을 이용하여 연산 횟수를 줄일 수 있다.

(19)

v_{i \pm | m | h} = \frac{1}{| m |} ((| m | - 1) v_{i \pm (| m | - 1) h} + \frac{1}{ρ_{i \pm (2 | m | - h)}}) v_{i + 2 h} = \frac{1}{2} (1 v_{i + h} + \frac{1}{ρ_{i + 3 h}}) v_{i + 3 h} = \frac{1}{3} (2 v_{i + 2 h} + \frac{1}{ρ_{i + 5 h}}) v_{i + 4 h} = \frac{1}{4} (3 v_{i + 3 h} + \frac{1}{ρ_{i + 7 h}})

비체적 또한 시간에 따라 변화하지 않으므로 미리 계산한 후 저장해서 사용할 수 있지만, 유한 차분 차수가 올라갈수록 더 많은 메모리가 필요하므로 계산 환경에 따라 비체적의 평균값을 다시 계산하는 것이 나을 수도 있다. 이 연구의 수치 예제에서는 평균 비체적을 시간 격자마다 새로 계산하였다.

차수에 따른 알고리듬의 정밀도 검증

개발된 알고리듬을 검증하기 위해 2차원 2층 무한 매질에서 파동 전파 해석해와 셀 기반 유한 차분법을 이용한 모델링 결과를 비교하였다. 상부층의 P파 속도는 2.0 km/s, 밀도는 1.8 g/cm³이고, 하부층의 P파 속도는 3.0 km/s, 밀도는 2.5 g/cm³이다. Fig. 3에 매질의 P파 속도와 밀도, 송수신기 위치를 표시하였다. 모델링을 위해 Fig. 4에 제시한 최대 주파수 80 Hz의 송신 파형을 사용하였다. 공간 격자 크기는 10 m로, 파장당 격자수는 최대 주파수 기준 2.5개, 최대 진폭 주파수 기준 10개에 해당한다. 샘플링 간격 1 ms로 1초간 모델링을 수행하였다. 이 연구에서는 고정단 경계 조건을 사용하였으나, 필요할 경우 다양한 흡수 경계 조건을 적용할 수 있다.

https://static.apub.kr/journalsite/sites/ksmer/2021-058-04/N0330580402/images/ksmer_58_04_02_F3.jpg

Fig. 3.

A source (star) and a receiver (triangle) in an infinite two-layer model.

https://static.apub.kr/journalsite/sites/ksmer/2021-058-04/N0330580402/images/ksmer_58_04_02_F4.jpg

Fig. 4.

The source wavelet and its frequency spectrum used in the numerical example.

2차원 무한 매질에서 음향파 파동 방정식의 해석해는 다음과 같다(Arfken et al., 2012).

(20)

\frac{i}{4} H_{0}^{(1)} (\frac{ω}{c} | X - X' |)

위 식에서 $i = \sqrt{1}$ , $H_{0}^{(1)}$ 는 1종 한켈함수, $ω$ 는 각주파수, $c$ 는 전파 매질의 상속도, $X$ 와 $X'$ 는 각각 수신기와 송신원 위치이다. 위의 해석해를 이용해 주파수 영역에서 직접파와 반사파의 그린 함수를 구하고 송신 파형을 곱한 후 역푸리에 변환하여 시간 영역 해석해를 구했다. 셀 기반 유한 차분법에서 공간 유한 차분 차수를 다르게 하여 모델링을 수행하고 Fig. 3의 수신기 위치에서 얻은 트레이스를 해석해와 비교하였다(Fig. 5). 공간 유한 차분 차수가 낮을 때는 분산이 크게 나타나지만 차수가 높아질수록 분산이 줄어들며 모델링 결과가 해석해로 수렴하는 것을 확인할 수 있다. Fig. 6에 모델링 시작 후 0.8초 시점에서의 전체 파동장을 제시하였다. 이 결과에서도 차수가 높아질수록 분산이 감소하는 것을 볼 수 있다.

https://static.apub.kr/journalsite/sites/ksmer/2021-058-04/N0330580402/images/ksmer_58_04_02_F5.jpg

Fig. 5.

The analytic solution and results of cell-based finite-difference methods with different orders.

https://static.apub.kr/journalsite/sites/ksmer/2021-058-04/N0330580402/images/ksmer_58_04_02_F6.jpg

Fig. 6.

Snapshots at 0.8 s from the cell-based finite-difference modeling in the order of (a) 2, (b) 4, (c) 6, and (d) 8.

모델링시 차수 증가에 따른 계산 시간은 Table 3에 제시하였다. 위 예제에서 격자 크기를 절반 크기인 5 m로 줄인 후 2차 셀 기반 차분식을 이용해 얻은 결과는 Fig. 7에 제시하였다. 이때 계산 시간은 3.8010 s로, 10 m 격자에서 8차 차분식을 사용한 결과(Table 3, Fig. 5)와 비교하면 작은 격자에서 2차 차분식을 사용하는 것보다 큰 격자에서 고차 차분식을 사용하는 것이 더 빠르고 정확함을 알 수 있다.

Table 3.

Order of the cell-based finite-difference schemes and its calculation time

Order	Time (s)
2	0.9817
4	1.8049
6	2.7748
8	3.0769
10	4.0976

https://static.apub.kr/journalsite/sites/ksmer/2021-058-04/N0330580402/images/ksmer_58_04_02_F7.jpg

Fig. 7.

The analytic al solution and the result of the second order cell-based finite-difference method with the grid size of 5 m.

분산 분석 및 안정성 분석

균질 매질에서 1차원 파동 방정식에 일반적인 고차 유한 차분식을 사용할 경우 $x$ 방향 2차 미분항은 다음과 같이 정리할 수 있다.

(21)

\frac{\partial}{\partial x} (\frac{1}{ρ} \frac{\partial P}{\partial x}) \approx \frac{1}{ρ} \frac{1}{Δ x^{2}} \sum_{m = - M}^{M} C_{m}^{2, N} P_{i + m}

고차 셀 기반 유한 차분식 식 (12)는 균질 매질의 경우 다음과 같이 정리할 수 있다.

(22)

\frac{\partial}{\partial x} (\frac{1}{ρ} \frac{\partial P}{\partial x}) \approx \frac{1}{ρ} \frac{1}{Δ x^{2}} [\sum_{\begin{matrix} m = - M \\ m \neq 0 \end{matrix}}^{M} (C_{m}^{2, N} P_{i + m}) - (\sum_{\begin{matrix} m = - M \\ m \neq 0 \end{matrix}}^{M} C_{m}^{2, N}) P_{i}]

위 식과 Table 2에서

(23)

C_{0}^{2, N} = - \sum_{\begin{matrix} m = - M \\ m \neq 0 \end{matrix}}^{M} C_{m}^{2, N}

이므로,

(24)

\frac{\partial}{\partial x} (\frac{1}{ρ} \frac{\partial P}{\partial x}) \approx \frac{1}{ρ} \frac{1}{Δ x^{2}} [\sum_{\begin{matrix} m = - M \\ m \neq 0 \end{matrix}}^{M} (C_{m}^{2, N} P_{i + m}) + C_{0}^{2, N} P_{i}] = \frac{1}{ρ} \frac{1}{Δ x^{2}} \sum_{m = - M}^{M} C_{m}^{2, N} P_{i + m}

이 되고, 위 식은 식 (21)에 제시된 기존의 고차 유한 차분식과 동일하다. 일반적으로 분산 및 안정성 분석시에는 균질 매질에 대해 분석하므로 고차 셀 기반 유한 차분식의 균질 매질에서의 분산 및 안정성 특성은 일반적인 유한 차분식의 분산 및 안정성 특성과 같아지게 된다. 일반적인 유한 차분식의 분산 및 안정성은 잘 알려져 있으므로(Alford et al., 1974; Cohen, 2002; Zauderer, 2006; Min et al., 2016; Ikelle and Amundsen, 2018) 안정성 결과만 제시하면 Table 4와 같다(Lines et al., 1999). Table 4에서 $p$ 는 Courant- Friedrichs-Lewy (CFL)수로, $Δ x = Δ y = Δ z$ 일 때 다음과 같이 정의한다.

(25)

p = \frac{v Δ t}{Δ x} \leq p_{m a x}

Table 4.

Stability limits $p_{m a x}$ for finite-difference schemes

Dimension	2nd order	4th order	6th order	8th order	10th order
1D	1.0000	0.8660	0.8134	0.7843	0.7654
2D	0.7071	0.6123	0.5752	0.5546	0.5412
3D	0.5773	0.5000	0.4696	0.4528	0.4419

적용성 분석

이 연구에서는 시간 영역 파동 전파 모델링에 고차 셀 기반 유한 차분법을 사용하였다. 그러나 셀 기반 유한 차분법은 공간 미분에만 사용하였으므로 이 연구에서 사용한 고차식 유도 방법은 주파수 영역 모델링에도 동일하게 적용 가능하다. 탄성파 파동 방정식 모델링에도 고차 셀 기반 차분법을 적용할 수 있으나, 탄성파에 적용하기 위해서는 공간에 대한 $x, y$ 방향, $y, z$ 방향 및 $z, x$ 방향 교차 미분항을 위해 별도의 차분식을 사용할 필요가 있다. 파형 역산에도 적용이 가능하나 이 경우 일반적인 고차 유한 차분법에 비해 그래디언트 계산 복잡도가 증가하게 된다. 일반적인 고차 유한 차분법을 사용하면 특정 격자에서의 파동장 계산에 해당 격자에서의 매질 특성만 사용하지만, 고차 셀 기반 유한 차분식을 사용하면 한 셀에서의 매질 특성이 주변 여러 격자상에서의 파동장 계산에 사용되므로 그래디언트 계산 식이 상대적으로 복잡해지게 된다.

결 론

이 연구에서는 P파 속도와 밀도를 고려한 음향파 파동 전파 모델링을 위한 고차 셀 기반 유한 차분법을 제안하였다. 셀 기반 유한 차분법을 이용하면 불연속적인 매질 특성에 대한 직접 미분을 피할 수 있으나 기존 2차 셀 기반 차분식의 경우 다른 일반적인 2차 차분식의 경우와 마찬가지로 파장당 격자수를 많이 사용해야했다. 이 연구에서는 2차 셀 기반 차분식을 고차 셀 기반 유한 차분식으로 개선하여 수치 모델링 시 격자 크기를 증가시킬 수 있게 함으로써 파장당 격자수를 줄일 수 있게 하였다. 수치 예제에서는 파동 전파 해석해와의 비교를 통해 고차 셀 기반 유한 차분식을 이용하면 파장당 격자수가 작아도 차수를 증가시킴으로써 정확한 결과를 얻을 수 있음을 증명하였다. 이러한 파장당 격자수의 제한을 완화시키는 것은 격자수에 따라 모델링 시간이 급격히 증가하는 3차원 모델링에서 더욱 효율적인 파동 전파 모델링을 가능하게 한다. 또한, 고차 셀 기반 유한 차분법은 시간 영역 음향파 파동 전파 모델링 뿐 아니라 주파수 영역 모델링과 탄성파 파동 전파 모델링에도 응용될 수 있다.

Acknowledgements

본 연구는 산업통상자원부(MOTIE)와 한국에너지기술평가원(KETEP)의 지원을 받아 수행한 연구 과제입니다(No. 20182510102470).

References

Alford, R.M., Kelly, K.R., and Boore, D.M., 1974. Accuracy of finite-difference modeling of the acoustic wave equation, Geophysics, 39(6), p.834-842. 10.1190/1.1440470

Arfken, G.B., Weber, H.J., and Harris, F.E., 2012. Mathematical Methods for Physicists: A Comprehensive Guide, Academic Press, Elsevier, p.447-468.

Cohen, G., 2002. Higher-order numerical methods for transient wave equations, Springer, p.35-164. 10.1007/978-3-662-04823-8_4

Ikelle, L. and Amundsen, L., 2018. Introduction to petroleum seismology, Society of Exploration Geophysicists, Tulsa, OK USA, p.828-844. 10.1190/1.978156080344729721953

Komatitsch, D. and Tromp, J., 1999. Introduction to the spectral element method for three-dimensional seismic wave propagation, Geophysical Journal International, 139(3), p.806-822. 10.1046/j.1365-246x.1999.00967.x

Lee, H.Y., Lim, S.C., Min, D.J., Kwon, B.D., and Park, M., 2009. 2D time-domain acoustic-elastic coupled modeling: a cell-based finite-difference method, Geosciences Journal, 13(4), p.407-414. 10.1007/s12303-009-0037-x

Lines, L.R., Slawinski, R., and Bording, R.P., 1999. A recipe for stability of finite‐difference wave‐equation computations, Geophysics, 64(3), p.967-969. 10.1190/1.1444605

Marfurt, K.J., 1984. Accuracy of finite-difference and finite- element modeling of the scalar and elastic wave equations, Geophysics, 49(5), p.533-549. 10.1190/1.1441689

Min, D.J. and Kim, H.S., 2006. Feasibility of the surface-wave method for the assessment of physical properties of a dam using numerical analysis, Journal of Applied Geophysics, 59, p.236-243. 10.1016/j.jappgeo.2005.09.004

Min, D.J., Pyun, S., Ha, W., Kwak, S., Chung, W., and Shin, C., 2016. Numerical analysis for Geophysics, CIR, Seoul, Korea, p.44-51.

Min, D.J., Shin, C., and Yoo, H., 2004. Free-surface boundary condition in finite-difference elastic wave modeling, Bulletin of the Seismological Society of America, 94(1), p.237-250. 10.1785/0120020116

Park, B. and Ha, W., 2019. Efficient 3D Acoustic Wave Propagation Modeling using a Cell-based Finite Difference Method, Geophysics and Geophysical Exploration, 22(2), p.56-61.

Schuster, G.T., 2017. Seismic Inversion, Society of Exploration Geophysics, Tulsa, OK USA, p.225-234.

Virieux, J., 1986. P-SV wave propagation in heterogeneous media: Velocity‐stress finite‐difference method, Geophysics, 51(4), p.889-901. 10.1190/1.1442147

Zauderer, E., 2006. Partial differential equations of applied mathematics, John Wiley & Sons, New Jersey, USA, p.768- 821. 10.1002/9781118033302

Journal of the Korean Society of Mineral and Energy Resources Engineers ISSN:2288-0291(Print) 2288-2790(Online) 한국자원공학회지

Preview

Acoustic Wave Propagation Modeling using High-order Cell-based Finite-difference Methods

ABSTRACT

MAIN

(1)

(2)

(3)

Fig. 1.

One-dimensional grids and cells.

(4)

(5)

(6)

(7)

(8)

(9)

Table 1.

Central finite-difference coefficients of the first derivative (Cm1,N), where N is the order of a finite-difference scheme, and m is the index of a finite-difference coefficient

(10)

(11)

(12)

Table 2.

Central finite-difference coefficients of the second derivative (Cm2,N)

(13)

(14)

(15)

(16)

(17)

Fig. 2.

Two-dimensional grids and cells for differentiation along the x-axis.

(18)

(19)

Fig. 3.

A source (star) and a receiver (triangle) in an infinite two-layer model.

Fig. 4.

The source wavelet and its frequency spectrum used in the numerical example.

(20)

Fig. 5.

The analytic solution and results of cell-based finite-difference methods with different orders.

Fig. 6.

Snapshots at 0.8 s from the cell-based finite-difference modeling in the order of (a) 2, (b) 4, (c) 6, and (d) 8.

Table 3.

Order of the cell-based finite-difference schemes and its calculation time

Fig. 7.

The analytic al solution and the result of the second order cell-based finite-difference method with the grid size of 5 m.

(21)

(22)

(23)

(24)

(25)

Table 4.

Stability limits pmax for finite-difference schemes

Acknowledgements

References

Central finite-difference coefficients of the first derivative ( $C_{m}^{1, N}$ ), where $N$ is the order of a finite-difference scheme, and $m$ is the index of a finite-difference coefficient

Central finite-difference coefficients of the second derivative ( $C_{m}^{2, N}$ )

Stability limits $p_{m a x}$ for finite-difference schemes