Journal Search Engine
Search Advanced Search Adode Reader(link)
Download PDF Export Citaion korean bibliography PMC previewer
ISSN : 1229-3059(Print)
ISSN : 2287-2302(Online)
Journal of the Computational Structural Engineering Institute of Korea
Vol.31 No.2 pp.87-95

DOI : https://doi.org/10.7734/COSEIK.2018.31.2.87

Explicit Transient Simulation of SH-waves Using a Spectral Element Method

Seungwook Youn1, Jun Won Kang1†
1Department of Civil Engineering, Hongik University, Seoul, 04066, Korea
Corresponding author: Tel: +82-2-320-1601; E-mail: jwkang@hongik.ac.kr
January 18, 2018 March 7, 2018 March 20, 2018

Abstract


This paper introduces a new explicit spectral element method for the simulation of SH-waves in semi-infinite domains. To simulate the wave motion in unbounded domains, it is necessary to reduce the infinite extent to a finite computational domain of interest. To prevent the wave reflection from the trunctated boundaries, perfectly matched layer(PML) wave-absorbing boundary is introduced. The forward problem for simulating SH-waves in PML-truncated domains can be formulated as second-order PDEs. The second-order semi-discrete form of the governing PDEs is constructed by using a mixed spectral elements with Legendregauss-Lobatto quadrature method, which results in a diagonalized mass matrix. Then the second-order semi-discrete form is transformed to a first-order, whose solutions are calculated by the fourth-order Runge-Kutta method. Numerical examples showed that solutions of SH-wave in the two-dimensional analysis domain resulted in stable and accurate, and reflections from truncated boundaries could be reduced by using PML boundaries. Elastic wave propagation analysis using explicit time integration method may be apt for solving larger domain problems such as three-dimensional elastic wave problem more efficiently.



스펙트럴 요소법을 이용한 SH파 전파의 외연적 시간이력해석

윤 승 욱1, 강 준 원1†
1홍익대학교 토목공학과

초록


이 논문에서는 스펙트럴 요소법과 외연적 시간적분법을 이용해 SH파의 전파 거동을 계산하는 수치해석 기법을 제시한다. 2차원 영역에서의 탄성파 해석을 위해 해석영역을 유한 영역으로 한정하고 파동이 반사되지 않도록 수치적 파동흡수 경계조 건인 perfectly matched layer(PML)를 도입하였다. PML이 포함된 시간영역 파동방정식의 유한요소해법을 위해 스펙트럴 요소법을 적용하였고 Legendre- Gauss-Lobatto 수치적분법을 사용하여 질량행렬을 대각화하였다. 2차 미분방정식 시스템 의 파동방정식을 1차 미분방정식 시스템으로 변환하였고 병렬화를 통한 탄성파 해석 성능의 최적화를 위해 외연적 시간적분 법인 4차 Runge-Kutta 방법을 이용해 해석영역에서의 변위응답을 계산하였다. 2차원 해석영역에서 SH파의 전파 거동을 계 산하는 수치예제를 통해 제시한 외연적 스펙트럴 요소법의 정확성을 검증하였고 PML로 인한 반사파의 감쇠효과를 확인하 였다. 외연적 시간적분법을 통한 탄성파 해석 기법은 3차원 영역과 같은 대규모 문제에서의 탄성파 수치해석을 효율적으로 수행할 수 있을 것으로 기대된다.



    National Research Foundation of Korea
    2017R1C1B2004975

    1. 서 론

    탄성파의 수치해석은 지진공학, 지반조사, 비파괴검사, 구조물의 손상평가 등의 다양한 분야에서 지속적으로 발전해 왔다(Epanomeritakis et al., 2008; Kang et al., 2011; 2013). 이러한 연구는 구조체의 미지 물성 분포를 도출하는 탄성파의 역해석(inverse analysis)을 중심으로 지속적으로 발전해 왔는데(Chew et al., 1990; Pratt et al., 1999; Operto et al., 2004; Epanomeritakis et al., 2008), 특히 고차원 해석의 경우 컴퓨터 메모리와 계산 시간 면에서 고도의 연산 기법을 필요로 한다. 따라서 한정된 계산 자원을 효율적으로 이용하여 대규모 문제의 수치해석에 적용될 수 있는 계산 알고리즘이 필요하다.

    반무한 영역에서의 탄성파 전파 수치해석을 위해서는 원래의 영역을 대신할 수 있는 유한 해석영역이 필요하다. 이를 위해 2차원 해석영역의 경계에서 반사파가 발생하지 않도록 수치적 파동흡수 경계조건인 perfectly matched layer(PML)(Kang and Kallivokas, 2010)를 도입하였으며 혼합유한요소법 (mixed finite element method)을 이용해 PML을 경계로 하는 유한영역에서의 SH파의 전파거동을 계산하는 정해석 (forward analysis)을 수행하였다.

    이 논문에서는 2차원 반무한 영역에서의 SH파 전파 거동을 계산하기 위해 PML을 고려한 시간영역 파동방정식의 해를 구하는 외연적 수치적 기법을 제시한다. 우선 파동방정식의 수치해를 구하기 위해 Legendre-Gauss-Lobatto 수치적분 기반 스펙트럴 요소를 사용해 질량행렬을 대각화한다. 그리고 2차 미분방정식 시스템으로 표현되는 semi-discrete form을 1차 미분방정식 시스템으로 변환하고 Runge-Kutta 방법을 이용한 외연적 시간적분을 통해 이 시스템의 수치해를 구한다. 외연적 시간적분법은 내연적 시간적분법에 비해 해의 안정성이 낮지만 계산 속도가 빠르고 효율적이라는 장점이 있다.

    수치예제를 통해 외연적 시간적분법을 이용한 SH파의 정해석 결과를 도출하고 해석영역에서 PML로 인한 반사파의 감쇠 효과를 확인한다. 그리고 외연적 시간적분법과 내연적 시간적 분법의 비교를 위해 Runge-Kutta 방법을 사용한 정해석 결과를 정해 및 Newmark-β 방법을 이용해 구한 수치해와 비교하여 정확도를 검토한다, 이 논문에서 제시하는 기법은 코드의 병렬화를 통한 탄성파의 최적 해석에 적용될 수 있다.

    2. PML을 경계로 하는 SH파의 정해석

    이 논문에서는 Fig. 1(a)와 같은 반무한 탄성영역에서의 SH파 전파 거동을 계산하고자 한다. SH파의 거동을 계산하기 위해 반무한 탄성영역의 양 측면과 아래 부분에 PML을 도입 하여 Fig. 1(b)와 같은 유한 해석영역을 구성하였다. PML을 경계조건으로 하는 2차원 유한영역에서의 SH파 전파거동은 다음과 같은 초기-경계값 문제로 표현될 수 있다(Kang and Kallivokas, 2010).(2d)(2e)(2f)

    f m 2 υ t 2 + c s g c υ t + c s 2 g k υ ( F ˜ e s t + F ˜ p s ) = 0
    (1a)

    F e 2 s t 2 + F p s t c s 2 υ t = 0 in Ω × ( 0 , T ]
    (1b)

    υ ( x , t ) = 0 on Γ fixed × ( 0 , T ]
    (2a)

    s 2 t ( x , t ) = p ( x , t ) on Γ free × ( 0 , T ]
    (2b)

    υ ( x , 0 ) = 0 on Ω
    (2c)

    υ t ( x , 0 ) = 0 on Ω
    (2d)

    s ( x , 0 ) = 0 on Ω
    (2e)

    s t ( x , 0 ) = 0 on Ω
    (2f)

    식 (1)은 변위(υ)와 응력이력(s1,s2)이 혼합된 연립미분 방정식으로서 PML을 포함하는 2차원 영역에서의 SH파 지배 방정식을 나타낸다(Kang and Kallivokas, 2010). 여기서 υ(= ρu)는 매질의 밀도를 포함하는 정규화된 변위이고, s i ( u = 1 , 2 ) 는 응력의 시간에 대한 적분값으로서 다음과 같이 표현된다.(3)

    s i = 0 t σ i ( x , τ ) d τ
    (3)

    c s = ( = μ / ρ ) 는 매질의 탄성파 속도이며, f m , g c , g k , F ˜ e , F ˜ p 는 다음과 같은 식들로 정의된다.(4)(5)(6)(7)(8)

    f m = ( 1 + f 1 e ) ( 1 + f 2 e )
    (4)

    g c = g 2 p ( 1 + f 1 e ) + g 1 p ( 1 + f 2 e )
    (5)

    g k = g 1 p g 2 p
    (6)

    F ˜ e = [ 1 + f 2 e 0 0 1 + f 1 e ]
    (7)

    F ˜ p = [ c s g 2 p 0 0 c s g 1 p ]
    (8)

    여기서, f i e f i p 는 각각 i축 (i =1,2) 방향으로의 소멸파 (evanescent wave)와 진행파(propagating wave)에 대한 감쇠함수이고, g i p 는 이러한 감쇠함수를 영역의 특성길이 b로 나눠 정규화한 함수이다. 감쇠함수 f i e 는 다음과 같은 식으로 표현된다.

    f i e , p ( i ) = { 0 , 0 i L 3 2 L P M L log ( 1 | R | ) ( i L L PML ) 2 L i L t
    (9)

    식 (9)에서 L 은 일반영역의 길이이고, Lt는 PML을 포함 하는 해석영역 전체 길이이며, LPML = (= Lt - L )은 PML의 길이를 나타낸다. R 은 파동이 PML 끝의 고정단에 도달했을 때 발생하는 반사파의 크기를 나타내는 탄성파의 이론적 반사 계수이다. 식 (2a)와 (2b)는 경계조건이고 식 (2c)~(2f)는 초기조건을 나타낸다.

    이러한 초기-경계값 문제에서 변위(υ)와 응력(s1, s2)은 혼합유한요소법을 이용하여 계산할 수 있다. 파동방정식 (1a)와 (1b)에 가중치 함수 w(x), p(x), q(x) 를 곱하고 전체 영역에 대해 적분한 후 식 (2a), (2b)의 경계조건을 적용하면 다음과 같은 변분식을 얻을 수 있다.(10a)(10b)(10c)

    Ω w ( f m 2 υ t 2 + c s g c υ t + c s 2 g k υ ) d Ω + Ω w ( F ˜ e s ˙ + F ˜ p s ) d Ω = Γ w ( F ˜ e s ˙ + F ˜ p s ) n d Γ
    (10a)

    Ω p { ( 1 + f 1 e ) 2 s 1 d t 2 c s 2 2 υ x 1 t + c s g 1 p s 1 t } d Ω = 0
    (10b)

    Ω q { ( 1 + f 2 e ) 2 s 2 d t 2 c s 2 2 υ x 2 t + c s g 2 p s 2 t } d Ω = 0
    (10c)

    식 (10)에서 변위와 응력의 근사함수인 υs1, s2 υ H 1 s 1 L 2 , s 2 L 2 의 함수공간에 속하고 서로 다른 기저함수(basis function) ϕ i ( x ) ψ i ( x ) 를 사용하여 다음과 같이 표현할 수 있다. 또한 가중치 함수 w ( x ) , p ( x ) , q ( x ) 역시 기저함수 ϕ i ( x ) ψ i ( x ) 를 사용하여 다음과 같이 나타 낼 수 있다.

    υ ( x , t ) ϕ ( x ) T v ( t )
    (11a)

    s 1 ( x , t ) ψ ( x ) T s 1 ( t )
    (11b)

    s 2 ( x , t ) ψ ( x ) T s 2 ( t )
    (11c)

    w ( x ) w T ϕ ( x )
    (12a)

    p ( x ) p T ψ ( x )
    (12b)

    q ( x ) q T ψ ( x )
    (12c)

    식 (11)과 (12)를 식 (10)에 대입하면 M u ¨ + C u ˙ + Ku = p ( t ) 형태의 운동방정식을 얻을 수 있고, 이 운동방정식을 계산하여 각 노드의 변위와 응력을 계산할 수 있다. 운동방정 식의 질량행렬, 감쇠행렬, 강성행렬, 힘벡터, 해벡터는 다음과 같이 표현된다.(14)(16)(17)

    M = [ Ω f m ϕ ϕ T d Ω 0 0 0 Ω ( 1 + f 1 e ) ψ ψ T d Ω 0 0 0 Ω ( 1 + f 2 e ) ψ ψ T d Ω ]
    (13)

    C = [ Ω c s g c ϕ ϕ T d Ω Ω ( 1 + f 2 e ) ϕ x 1 ψ T d Ω Ω ( 1 + f 2 e ) ϕ x 1 ψ T d Ω Ω c s 2 ψ ϕ T x 1 d Ω Ω c s g 1 p ψ ψ T d Ω 0 Ω c s 2 ψ ϕ T x 2 d Ω 0 Ω c s g 2 p ψ ψ T d Ω ]
    (14)

    K = [ Ω c s 2 g k ϕ ϕ T d Ω Ω c s g 2 p ϕ x 1 ψ T d Ω Ω c s g 1 p ϕ x 2 ψ T d Ω 0 0 0 0 0 0 ]
    (15)

    p ( t ) = [ Γ free ϕ p d Γ free 0 0 ]
    (16)

    u = [ v s 1 s 2 ]
    (17)

    3. 외연적 시간적분법을 이용한 SH파 정해석

    이 장에서는 4차 Runge-Kutta 방법과 Legendre- Gauss- Lobatto 수치적분법을 이용하여 SH파의 전파 거동을 계산 하는 방법을 소개한다.

    3.1. 4차 Runge-Kutta 방법

    Runge-Kutta 방법은 1차 미분방정식을 푸는 방법으로서 이 수치적분 방법을 이용해 앞서 유도한 semi-discrete 형태의 운동방정식의 해를 구할 수 있다. 이 논문에서는 여러 Runge- Kutta 방법들 중에서 가장 일반적으로 사용되는 4차 Runge- Kutta 방법을 이용하였다. 4차 Runge-Kutta 방법을 이용해 1차 미분방정식의 해를 계산하는 과정은 다음과 같은 식으로 표현된다.(18)(19a)(19b)(19c)

    u ˙ = f ( t , u )
    (18)

    k 1 = f ( t i , u i )
    (19a)

    k 2 = f ( t i + 1 2 Δ t , u i + 1 2 k 1 Δ t )
    (19b)

    k 3 = f ( t i + 1 2 Δ t , u i + 1 2 k 1 Δ t )
    (19c)

    k 4 = f ( t i + Δ t , u i + 1 2 k 3 Δ t )
    (19d)

    u i + 1 = u i + 1 6 ( k 1 + 2 k 2 + 2 k 3 + k 4 ) Δ t
    (20)

    여기서, u ˙ 은 변수 u의 시간에 대한 1차 미분이다. 식 (19)에서 tii번째 시간이고, uii번째 시간에서의 변위이며 Δ t = ( = t i + 1 t i ) 는 시간 증분을 나타낸다. k 1 , k 2 , k 3 , k 4 는 기울 기에 의한 u의 증분을 나타내는데, 식 (20)을 사용하여 각각 의 증분에 가중치를 곱해 i + 1번째 시간에서의 변위 ui+1을 계산할 수 있다.

    앞서 유도된 운동방정식을 Runge-Kutta 방법으로 풀기 위해서는 2차 미분방정식 시스템을 1차 미분방정식 시스템 으로 변형시켜야 한다. 그렇게 변형된 미분방정식 시스템의 1차 미분항에 대한 계수행렬을 단위행렬 I로 만든 후 Runge-Kutta 방법에 의해 해를 구할 수 있다.

    첫째로, SH파 운동방정식을 1차 미분방정식 시스템으로 표현하기 위해 2차 미분방정식 시스템인 M u ¨ + C u ˙ + Ku = p ( t ) 에서 u ˙ y로 치환한다. 그러면 2차 미분방정식 시스템을 식 (21)과 같이 2개의 1차 미분방정식 시스템으로 나타낼 수 있다. 이를 행렬 형태로 표현하면 식 (22)와 같다.(21a)(21b)

    u ˙ y = 0
    (21a)

    M y ˙ + Cy + Ku = p ( t )
    (21b)

    [ I 0 0 M ] [ u ˙ y ˙ ] + [ 0 -I K C ] [ u y ] = [ 0 p ( t ) ]
    (22)

    둘째로, 식 (22)에서 1차 미분항의 계수행렬을 단위행렬로 만들기 위해 질량행렬 M을 단위행렬 I로 변형시킨다. 질량행렬 M을 단위행렬 I로 만들기 위해서는 M을 대각화하는 과정이 필요하다. 질량행렬 M을 대각화하기 위한 가장 간단한 방법은 집중질량을 가정하는 것이다. 하지만 보다 높은 정확도를 확보 하기 위해 이 논문에서는 Legendre-Gauss- Lobatto 수치적 분법을 이용해 질량행렬 M을 대각화하였으며, 이를 위해 9개의 노드를 가진 2차원 사각형 요소를 사용하여 방정식의 변위 및 응력 해를 근사시켰다.

    3.2. Legendre-Gauss-Lobatto 수치적분 기반 스펙트럴 요소법

    Legendre-Gauss-Lobatto 수치적분법은 고차원의 수치 적분을 효율적으로 할 수 있는 방법이다. 이 논문에서는 9절점 사각형 유한요소 국부좌표계의 각 방향에 대해 3개의 적분점을 사용해 수치적분을 수행하였으며, 이에 대한 노드의 위치, 적분점의 위치, 가중치를 Table 1에 나타냈다. Table 1에서 볼 수 있듯이 Legendre-Gauss-Lobatto 수치적분법은 한 요소에 대해 노드의 위치와 적분점의 위치가 같다. Legendre- Gauss-Lobatto 수치적분법에서 적분점의 위치와 가중치를 계산하는 방법은 다음과 같다.

    적분점 ( L e g e n d r e G a u s s L o b a t t o p o i n t s ) : x 0 = 1 , x j = zeros of L N , x N = 1 1 j N 1
    (23)

    가중치 ( L e g e n d r e G a u s s L o b a t t o w e i g h t s ) : w j = 2 N ( N + 1 ) 1 [ L N ( x j ) ] 2 j = 0 , , N
    (24)

    식 (23), (24)에서 N 은 Legendre-Gauss-Lobatto 수치 적분법의 차수를 나타내고, 적분점의 개수+1로 표현할 수 있다. N 차 Legendre-Gauss-Lobatto 수치적분법의 경우 2N -1차 다항식까지 정확하게 적분할 수 있다. LN 은 르장 드르 다항식(legendre polynomials)을 의미하며 다음과 같은 식으로 표현할 수 있다.(25a)(25b)(25c)

    L 0 ( x ) = 1
    (25a)

    L 1 ( x ) = x
    (25b)

    L N + 1 ( x ) = 2 N + 1 N + 1 x L N ( x ) N N + 1 L N 1 ( x )
    (25c)

    Legendre-Gauss-Lobatto 수치적분법에서 사용되는 기저 함수 역시 르장드르 다항식이 포함된 식으로 다음과 같이 표현 할 수 있다.(26)

    ϕ i = 1 N ( N + 1 ) L N ( x i ) ( 1 x 2 ) L N ( x ) x x i
    (26)

    이 논문에서는 Legendre-Gauss-Lobatto 수치적분법을 사용해 식 (13)~(15)의 질량행렬 M, 감쇠행렬 C, 강성행렬 K를 구성하였다. Legendre-Gauss-Lobatto 수치적분법을 사용하여 질량행렬 M을 대각화하면 식 (21)을 다음과 같은 2개의 1차 미분방정식 시스템으로 바꿀 수 있다.(27a)(27b)

    u ˙ = y
    (27a)

    y ˙ = M 1 ( Cy Ku + p ( t ) )
    (27b)

    식 (27)을 식 (28)과 같이 행렬 형태로 표현할 수 있고, 식 (29)의 초기조건을 이용하여 Runge-Kutta 방법에 적용 하면 변위(υ)와 응력(s1, s2)을 구할 수 있다.

    [ u ˙ y ˙ ] = [ 0 -I M -1 K M -1 C ] [ u y ] + [ 0 M 1 p ( t ) ]
    (28)

    u ( 0 ) = 0
    (29a)

    y ( 0 ) = 0
    (29b)

    4. 수치예제

    4.1. SH파 정해석 결과

    이 논문에서 제시한 2차원 영역에서의 SH파 정해석 알고리 즘의 정확도를 검증하기 위해 Fig. 2와 같이 수직방향으로 균질한 물성을 갖는 반무한 고체영역의 전파 거동을 계산하였다. 해석영역의 밀도는 2200kg/m3이고 탄성파 속도는 100m/s 이다.

    Fig. 3과 같이 예제로 사용된 탄성영역의 크기는 60m× 30m이며, 탄성영역의 양 측면과 아래 부분에 6m 길이의 PML을 도입하였다. 해석 영역의 요소 길이는 탄성영역과 PML에서 각각 1m, 0.5m이다.

    Fig. 4는 반무한 영역의 표면에 가해진 하중의 시간이력을 나타낸다. 전체 해석 시간은 2s이고, 시간 단계는 0.0025s이다. 시간에 대해 Gaussian 분포를 가지는 하중을 탄성영역 표면 -2m≤x≤2m 위치에 작용시켰으며 이는 다음과 같은 함수로 표현된다.

    p ( t ) = 10000 exp [ ( t 0.22 ) 2 / 0.0027 ]
    (30)

    Fig. 5는 하중이 가해졌을 때 Runge-Kutta 방법을 이용 하여 계산한 해석영역 표면(x =0m)에서의 변위응답의 시간 이력을 나타낸다. Fig. 5는 각각 A(0,0), B(10,0), C(20,0), D(30,0) 위치에서 변위응답의 시간이력을 나타내며 각각의 위치를 Fig. 3에 표시하였다.

    Fig. 6은 각각 t =0.2s, 0.4s, 0.6s, 0.8s일 때 전체 해석 영역에서의 SH파형을 나타낸다. Runge-Kutta 방법과 PML을 연동한 해석 결과 반무한 영역에서 SH파가 잘 전파되는 것을 확인할 수 있었으며, 파동이 양 측면과 하단의 PML 영역에서 효과적으로 흡수되는 것을 확인할 수 있었다.

    4.2. PML의 반사파 감쇠 효과

    2차원 영역에서 PML에 의한 반사파 감쇠 효과를 검증하기 위해 Fig. 7, 8과 같은 두 영역에서 동일한 방법으로 SH파 전파거동을 계산하였다. Fig. 7은 기존의 해석영역에서 PML을 제거하고 x, y 방향으로 각각 5배씩 확대한 영역으로서 전체 해석시간 2s 동안 기존의 해석영역에 반사파가 도달하지 않도록 설정하였다. Fig. 8은 기존의 해석영역과 크기는 같지만 PML 이 없이 고정 경계조건을 가지는 영역이다.

    Fig. 9는 PML영역과 확장된 영역에서 계산한 변위응답의 시간이력을 A(0,0), B(10,0), C(20,0), D(30,0) 위치에서 비교한 결과를 나타낸다. 두 변위응답을 비교한 결과가 거의 일치하는 것을 확인할 수 있다. 따라서 이 논문에서 제시한 PML을 경계조건으로 하는 정해석 알고리즘의 반사파 감쇠 효과가 뛰어나다고 평가된다.

    Fig. 10은 PML영역과 고정된 경계조건을 갖는 영역에서 계산한 변위응답의 시간이력을 A(0,0), B(10,0), C(20,0), D(30,0) 위치에서 비교한 결과를 나타낸다. 두 변위응답을 비교한 결과 고정된 경계로부터 발생하는 반사파의 존재를 확인할 수 있으며, PML의 도입으로 반사파가 발생하지 않는 것을 확인할 수 있다.

    4.3. 외연적 시간적분 결과와 내연적 시간적분 결과의 비교

    Runge-Kutta 방법을 이용한 외연적 시간적분 방법의 타당 성을 검토하기 위해 이 방법을 사용해 구한 결과를 Newmark- β 방법을 이용한 내연적 시간적분 결과와 비교하였다. Fig. 3과 같은 영역에서 SH파의 전파거동을 계산하였으며, 하중은 식 (30)과 동일한 Gaussian 형태를 사용하였다. 각 해석 영역의 노드 및 요소 정보를 Table 2에 나타냈다.

    Fig. 11은 A(0,0), B(10,0), C(20,0), D(30,0) 위치에서 Runge-Kutta 방법과 Newmark-β 방법을 사용하여 계산한 변위응답을 비교한 결과를 나타낸다. 변위응답을 비교한 결과 두 가지 방법이 유사한 결과를 보여 외연적 시간적분법의 정확 성과 안정성을 확인할 수 있다.

    외연적 시간적분법과 내연적 시간적분법의 정확도를 검증 하기 위해 Runge-Kutta 방법과 Newmark-β 방법을 이용 하여 계산한 근사해와 SH파의 정해를 계산하여 비교하였다. SH파의 정해는 시간영역 SH파의 Green’s function을 사용 하여 식 (31)으로 계산하였고 trapzoidal rule을 이용하여 적분하였다. 하중은 Fig. 12와 같은 단위 충격하중을 사용 하였다.

    u ( x , t ) = x 1 x 2 1 2 π μ H ( t t s ) t 2 t s 2 d x
    (31)

    t s = r β
    (32)

    여기서, μ는 전단탄성계수를 나타내고 H 는 Heaviside step function을 나타낸다. x는 정해를 계산하고 싶은 위치를 나타 내고 x1x2는 분포하중이 시작되는 지점과 끝나는 지점을 나타낸다. t는 계산하고자 하는 시간을 나타내고 ts 는 계산 하고자 하는 위치에 탄성파가 도달하는 시간으로서 식 (32)와 같이 나타낼 수 있다. r은 하중이 가해지는 지점으로부터 계산 하고자 하는 위치까지의 거리를 나타내고 β는 전단파속도를 나타낸다. 이와 같은 과정으로 계산한 정해를 Runge-Kutta 방법과 Newmark-β 방법을 통해 계산한 변위응답과 비교 하여 Fig. 13과 Table 3, 4에 나타냈다.

    Fig. 13은 (10,0), (20,0), (30,0), (30,10), (30,20) 위치에서 정해와 근사해의 시간이력을 비교한 결과를 나타낸다. 각각의 위치에서 정해와 근사해가 유사한 결과를 나타내는 것을 확인할 수 있다. 하지만 Fig. 13(a)에서 Runge-Kutta 방법을 이용한 해석 결과에 떨림 현상이 발생하였다. 이는 외연적 시간적분 방법이 내연적 시간적분 방법에 비해 해의 안정성이 낮아서 발생한 것으로 판단된다.

    Table 3은 (10,0), (20,0), (30,0), (30,10), (30,20) 위치에서 최대 변위가 발생하는 시간에서의 정해와 근사해를 비교한 결과를 나타낸다. 각각의 위치에서 변위를 비교하여 오차를 계산한 결과 모든 위치에서 Runge-Kutta 방법을 사용한 결과가 Newmark-β 방법을 사용한 결과보다 오차가 적은 것을 확인할 수 있었다.

    Table 4는 1.0s s 일 때 (10,0), (20,0), (30,0), (30,10), (30,20) 위치에서 정해와 근사해를 비교한 결과를 나타낸다. 이 경우 해석이 많이 진행된 상태로 해의 크기가 매우 작기 때문에 두 가지 방법 모두 오차가 크게 발생하였다. 동일한 방법으로 오차를 계산한 결과 Runge-Kutta 방법의 경우 약 22%의 오차를 나타낸 반면 Newmark-β 방법의 경우 약 26%의 오차가 발생하여 Runge-Kutta 방법으로 계산한 해의 정확성이 더 높음을 알 수 있었다.

    4. 결 론

    이 논문에서는 균질한 2차원 반무한 탄성체를 대상으로 PML 경계조건이 도입된 영역에서 스펙트럴 요소법과 외연적 시간적분법을 이용한 SH파 정해석 알고리즘을 제시하였다. 그리고 2차원 탄성영역에 대한 SH파의 전파거동을 계산하는 수치예제를 통해 정해석 알고리즘을 검증하였다.

    Legendre-Gauss-Lobatto 수치적분법을 기반으로 한 스펙트럴 요소법을 사용해 질량행렬 M을 대각화했다. 이를 바탕으로 2차 미분방정식 시스템의 파동방정식을 1차 미분방 정식 시스템으로 변환하였으며 외연적 시간적분 방법인 Runge- Kutta 방법에 적용하여 수치해를 계산하였다.

    이 논문에서 제시한 방법으로 구한 변위응답을 확장된 영역과 고정된 경계조건을 가진 영역에서의 변위응답과 비교하여 PML 경계에서 파동이 감쇠하는 것을 확인하였다. 그리고 내연적 시간적분법 중 하나인 Newmark-β 방법을 이용한 변위응답 및 정해와 비교하여 제시한 외연적 스펙트럴 요소법의 정확도 및 안정성을 검증하였다.

    이 논문의 연구는 3차원과 같은 고차원의 탄성파 수치해석을 효율적으로 수행하기 위한 기반 연구로 코드의 병렬화를 통한 계산 성능의 최적화가 가능할 것으로 기대된다.

    감사의 글

    이 연구는 2017년도 미래창조과학부 재원에 의한 한국연구 재단의 지원(NRF-2017R1C1B2004975)에 의하여 수행되 었습니다. 연구비 지원에 깊은 감사를 드립니다.

    Figure

    COSEIK-31-87_F1.gif

    (a) Semi-infinite domain; (b) Two-dimensional PML-truncated computational domain

    COSEIK-31-87_F2.gif

    Schematic of a homogeneous semi-infinite solid medium

    COSEIK-31-87_F3.gif

    A two-dimensional PML-truncated domain

    COSEIK-31-87_F4.gif

    Time history of a surface stress load p(t)

    COSEIK-31-87_F5.gif

    Time history of displacement at sampling points

    COSEIK-31-87_F6.gif

    Snapshots of displacement

    COSEIK-31-87_F7.gif

    Two-dimensional enlarged domain

    COSEIK-31-87_F8.gif

    Two-dimensional fixed boundary domain

    COSEIK-31-87_F9.gif

    Comparison of the displacement u calculated in the PML-truncated and enlarged domains

    COSEIK-31-87_F10.gif

    Comparison of the displacement u calculated in the PML-truncated and fixed domains

    COSEIK-31-87_F11.gif

    Comparison of the displacement u calculated by the developed explicit method and the Newmark-β method

    COSEIK-31-87_F12.gif

    Time history of a unit impulse load

    COSEIK-31-87_F13.gif

    Comparison of exact solution and approximate solutions

    Table

    Legendre-Gauss-Lobatto quadrature rule

    Mesh information in the wave propagation analysis

    Errors obtained by the developed explicit method and the Newmark-β method at peak point

    Error obtained by the developed explicit method and the Newmark-β method in 1.0s

    Reference

    1. Epanomeritakis, I ., Akcelik, V ., Ghattas, O ., Bielak, J . (2008) A Newton-CG Method for Large-scale Three-dimensional Elastic Full- waveform Seismic Inversion, Inverse Problems24(034015),
    2. Fathi, A ., Kallivokas, LF ., Poursartip, B . (2015) Full-waveform Inversion in Three- dimensional PML-truncated Elastic Media, Comput. Methods Appl. Mech. & Eng.296,pp.39~72.
    3. Fathi, A ., Poursartip, B ., Kallivokas, L.F . (2015) Time-domain Hybrid Formulations for Wave Simulations in Three-dimensional PML- truncated Heterogeneous Media, Int. J. Numer. Methods Eng.101(3), pp.165~198.
    4. He, Q ., Jiao, D . (2011) An Explicit Time-domain Finite-element Method that is Unconditionally Stable, Antennas and Propagation (APSURSI), 2011 IEEE Int. Symp.pp.2976~2979.
    5. Kang, J.W ., Kallivokas, LF . (2010) Mixed Unsplit-field Perfectly-matched-layers for Transient Simulations of Scalar Waves in Heterogeneous Domains, Comput. Geosci.14(4), pp.623~648.
    6. Kang, J.W ., Kallivokas, LF . (2011) The Inverse Medium Problem in Heterogeneous PML-truncated Domains Using Scalar Probing Waves, Comput. Methods Appl. Mech. & Eng.200(1-4), pp.265~283.
    7. Kang, J.W . (2013) Time-domain Elastic Full-waveform Inversion Using One-dimensional Mesh Continuation Scheme, J. Comput. Struct. Eng. Inst. Korea26(4), pp.213~221.
    8. Operto, S ., Ravaut, C ., Improta, L ., Virieux, J ., Herrero, A ., Dell'Aversana, P . (2004) Quantitative Imaging of Complex Structures from Dense Wide-aperture Seismic Data by Multiscale Traveltime and Waveform Inversions: a Case Study, Geophys. Prospect.52(6), pp.625~651.
    9. Pratt, R.G ., Shipp, R.M . (1999) Seismic Waveform Inversion in the Frequency Domain, Part 2: Fault Delineation in Sediments Using Crosshole Data, Geophys.64(3), pp.902~914.
    10. Van de Vosse, F.N ., Minev, P.D . (1996) Spectral Elements Methods: Theory and Applications, EUT Report 96-W-001, Eindhoven University of Technology.