Journal Search Engine
Search Advanced Search Adode Reader(link)
Download PDF Export Citaion korean bibliography PMC previewer
ISSN : 1226-525X(Print)
ISSN : 2234-1099(Online)
Journal of the Earthquake Engineering Society of Korea Vol.22 No.7 pp.379-390
DOI : https://doi.org/10.5000/EESK.2018.22.7.379

Practical Numerical Model for Nonlinear Analyses of Wave Propagation and Soil-Structure Interaction in Infinite Poroelastic Media

Lee Jin Ho1)*
1)Assistant Professor, Department of Ocean Engineering, Pukyong National University
Corresponding author: Lee, Jin Ho E-mail: jholee0218@pknu.ac.kr
July 23, 2018 October 4, 2018 October 5, 2018

Abstract


In this study, a numerical approach based on mid-point integrated finite elements and a viscous boundary is proposed for time-domain wave-propagation analyses in infinite poroelastic media. The proposed approach is accurate, efficient, and easy to implement in time-domain analyses. In the approach, an infinite domain is truncated at some distance. The truncated domain is represented by mid-point integrated finite elements with real element-lengths and a viscous boundary is attached to the end of the domain. Given that the dynamic behaviors of the proposed model can be expressed in terms of mass, damping, and stiffness matrices only, it can be implemented easily in the displacement-based finite-element formulation. No convolutional operations are required for time-domain calculations because the coefficient matrices are constant. The proposed numerical approach is applied to typical wave-propagation and soil-structure interaction problems. The model is verified to produce accurate and stable results. It is demonstrated that the numerical approach can be applied successfully to nonlinear soil-structure interaction problems.



무한 다공성 매질에서의 비선형 파전파 해석과 지반- 구조물 상호작용 해석을 위한 실용적 수치 모형

이진호1)*
1)부경대학교 해양공학과 조교수

초록


    Pukyong National University

    1. 서 론

    토목 및 기계공학에서는 다공성 매질에서의 파전파 현상에 대하여 많은 관심을 가져왔는데, 이는 이와 관련된 물리 현상이 다양한 공학적 시스템에 활용되고 있기 때문이다. 이에 대한 전형적인 예로는 다공성 반무한 지반에 서의 지반-구조물 상호작용, 지반의 액상화 해석, 탄성파를 이용한 지반 조 사, 다공성 물질로 이루어진 구조물의 비파괴 검사 등을 들 수 있다. 그러므 로, 다공성 매질의 동역학을 이해하기 위해 많은 연구가 이루어져 왔다[1]. 무한 다공성 매질에서의 파전파 해석 시 가장 중요하게 고려하여야 할 요소 중의 하나가 무한 매질에서의 에너지 방사 현상이다. 만약 이를 적절하게 고 려하지 않는다면 다공성 탄성파(poroelastic wave)가 한정된 영역에 갇히 게 되어 최종 응답이 실제와 많은 차이를 보이게 되기 때문이다. 그러므로, 무한 매질의 효과를 정확하면서도 효율적인 방법으로 고려하는 것은 무한 매질에서의 파전파 해석기법 개발의 가장 중요한 연구 주제이다. 무한 매질 에서의 파전파 문제의 정해는 주파수영역에서 얻을 수 있고, 이를 위해 전달 경계[2], 경계요소[3, 4], 무한요소[5] 등과 같은 수치모형이 주파수영역에 서 개발되어 사용되어 왔다.

    시스템의 비선형 거동은 다공성 매질의 파전파 문제에서 중요한 고려 사 항 중 하나이다. 예를 들면, 강진 지반운동은 지반-구조물 상호작용 시스템의 큰 변형과 이로 인한 지반의 액상화를 유발할 수 있는데, 이는 그 시스템의 안 전성에 큰 영향을 미친다. 재료 비선형성을 포함한 비선형 거동은 시간영역 에서 비선형 유한요소를 사용하여 가장 잘 나타낼 수 있다. 그러므로, 무한 다공성 매질에서의 비선형 파전파 해석의 정확성과 효율성은 무한 영역으로 의 에너지 방사를 시간영역에서 얼마나 정확하고 효율적으로 나타낼 수 있 는지에 의하여 결정된다. 앞에서 언급한 방법들은 주파수 영역에서 개발되 었기 때문에, 이들을 시간영역 해석에 적용하면 convolutional operator의 형태를 가지게 된다. Convolution operator는 현재 시간에서의 응답을 계산 하기 위해 계산을 시작하는 순간부터의 모든 응답 이력을 필요로 하므로, 정 확한 해를 얻을 수는 있지만 효율적인 방법이라고 할 수는 없다. 시간영역에 서 효율적인 지반-구조물 상호작용 해석을 수행하기 위해서는 convolution operator와 같은 global temporal operator보다는 local temporal operator 가 선호되고, 이러한 형태로 표현되는 수치 모형이 필요하게 된다. 고차 흡수 경계조건[6]과 perfectly matched layer(PML)[7, 8]는 시간영역에서 local temporal operator의 형태로 표현할 수 있는 접근법이고, 최근 무한영역에서 의 파전파 문제에 다양하게 활용되고 있다. 특히, 고차 흡수경계조건과 PML 의 장점을 모두 가지고 있는 perfectly matched discrete layer(PMDL)를 이 용하여 반무한 지반에서의 비선형 파전파 해석과 지반-구조물 상호작용 해 석을 위한 수치적 방법이 개발되기도 하였다[9-11].

    고차 흡수경계조건과 PML이 무한 매질에서의 시간영역 파전파 해석에 정확하고 효율적인 방법이지만, 이 방법들을 구현할 때는 특별히 고려해야 할 사항들이 있다. 고차 흡수경계조건은 무한영역의 분산방정식을 유리식 으로 근사하거나 간단한 미분연산자의 조합으로 표현하는 방법이다. 하지 만, 이를 구현하기 위해서는 물리적인 의미가 없는 추가 변수를 사용하여야 한다[6]. PML은 복소좌표변환의 개념을 사용하여 매질에 인위적인 감쇠 를 부여하여 진행하는 파의 진폭을 감소시키는 방법이다. PML은 주로 응 력-속도 정식화를 사용하여 구현되었는데[12], 대부분의 엔지니어는 변위 기반 정식화에 익숙하다. 변위 기반 정식화를 사용하여 PML을 구현할 수 있지만, 이때는 변위를 시간에 대하여 적분한 항이 발생한다[8]. 참고로 변 위를 시간에 대하여 적분한 항은 고차 흡수경계조건에서도 관찰되는 현상 이기도 하다[9-11, 13-15]. 즉, 고차 흡수경계조건과 PML을 구현하기 위 해서는 특별히 고려해야 할 사항들이 있고, 이는 이 방법들의 사용성을 저하 시키는 요인이다. 또한, 이 방법들의 수치 안정성에 대해서는 아직 명확한 결론이 나지 않았고[16, 17], 이는 고차 흡수경계조건과 PML의 활용성을 저해하는 또 다른 요인이기도 하다.

    이 연구에서는 무한 다공성 매질에서의 시간영역 파전파 해석에 활용할 수 있는 정확하고 효율적이면서도 구현하기 쉬운 수치적 방법을 제안하고 자 한다. 이를 위해 다공성 탄성 매질에 대한 mid-point integrated finite element와 점성경계에 근거한 방법을 제안할 것이다. 기존 연구에서 mid-point integrated finite element와 PMDL의 조합을 사용하면 층상 반 무한체를 정확하고 효율적으로 모사할 수 있음을 이미 보였다[18]. Midpoint integrated finite element는 통상적인 유한요소와 유사하게 질량, 감 쇠, 강성 행렬을 사용하여 정식화할 수 있다. 하지만, PMDL은 시간영역 정 식화를 위해서는 변위를 시간에 대하여 적분한 항을 필요로 한다[9-11, 13-15]. 그러므로, 이 연구에서는 무한 매질에 대한 PMDL을 시간영역에 서 간단한 감쇠기로 표현할 수 있는 점성경계로 대체한 수치 모형을 제안하 고자 한다. 점성경계는 기존 상용 유한요소 해석코드에 이미 구현되어 있고, mid-point integrated finite element는 유한요소의 개념을 사용하여 기존 해석코드에 사용자 요소로 쉽게 구현이 가능하기 때문에 제안된 수치 모형 은 무한 다공성 매질에서의 파전파 문제에 대하여 아주 “실용적”인 접근법 이라고 할 수 있다. 이 논문에서는 무한 다공성 탄성 매질에 대한 mid-point integrated finite element와 점성경계에 대해 간략히 살펴보고, 이를 사용 하여 다공성 반무한체를 모사하는 방법을 제안할 것이다. 그 후 제안된 수치 모형의 정확성을 검토하고 다공성 지반에서의 비선형 동적 해석에 적용할 것이다.

    2. 무한 다공성 탄성 매질의 수치 모형

    이 절에서는 무한 다공성 매질에서의 파전파 문제에 적용할 수 있는 “실 용적”인 수치 모형을 제안하고자 한다. 이 모형은 mid-point integrated finite element와 점성경계에 근거하고 있다. 그러므로, 두가지 구성 요소 에 대해 살펴본 후, 무한 다공성 탄성 매질을 나타내는 방법에 대해 제안하 고자 한다. 이 연구에서는 다공성 탄성 매질의 2차원 평면 변형 문제만 다룰 것이지만, 일반적인 3차원 문제로의 확장을 어렵지 않게 수행할 수 있다.

    2.1. 다공성 탄성 매질에 대한 Mid-point Integrated Finite Element

    직교좌표계에서 다공성 매질의 지배방정식은 다음과 같다[19-21]:

    σ z + A ˜ 2 u x 2 + B ˜ 2 u x z M ˜ 2 u t 2 C ˜ u t = 0
    (1)

    σ ( x , z , t ) = [ τ x z ( x , z , t ) σ z ( x , z , t ) 0 p ( x , z , t ) ] T
    (2a)

    u ( x , z , t ) = [ u x ( x , z , t ) u z ( x , z , t ) w x ( x , z , t ) w z ( x , z , t ) ] T
    (2b)

    A ˜ = [ λ + 2 μ + Q α 2 0 Q α 0 0 μ 0 0 Q α 0 Q 0 0 0 0 0 ]
    (3a)

    B ˜ = [ 0 λ + Q a 2 0 Q a μ 0 0 0 0 Q a 0 Q 0 0 0 0 ]
    (3b)

    M ˜ = [ ρ 0 ρ w 0 0 ρ 0 ρ w ρ w 0 ρ w n 0 0 ρ w 0 ρ w n ]
    (3c)

    C ˜ = [ 0 0 0 0 0 0 0 0 0 0 f 0 0 0 0 f ]
    (3d)

    여기서, ux (x, z, t) 와 uz (x, z, t) 는 고체 골격의 변위, wx (x, z, t) 와 wz (x, z, t) 는 공극율을 곱한 고체 골격에 대한 유체의 상대 변위이다. σz (x, z, t) 와 τxz (x, z, t) 는 혼합체의 전응력, p (x, z, t) 는 유체의 간극수압이다. λ와 μ 는 고체 골격의 Lamé 상수, Qα는 다공성 매질의 압축성을 고려하기 위 한 Biot 상수이다. ρw 는 유체의 밀도, ρ = ( 1 n ) ρ s + n p w 는 혼합체의 평균 밀도인데, 여기서 ρs 는 고체의 밀도이다. n은 공극율, f = 1/κ 인데, 여기 서 κ 는 투수계수이다. 식 (2a)의 응력과 간극수압은 식 (2b)의 변위와 관련 된다.(5)

    σ = B ˜ T u x + G ˜ u x
    (4)

    G ˜ = [ μ 0 0 0 0 λ + 2 μ + Q α 2 0 Q α 0 0 0 0 0 Q α 0 0 ]
    (5)

    Fig. 1a에 보인 바와 같이 z = z0 = 0 부터 z = zN = Lz 까지의 영역을 고 려해보자. 이 영역에서의 지배방정식을 근사하기 위해, 이 영역을 N개의 층 으로 나눌 것이다. z j 1 z z j j번째 층에서 z 방향으로의 변위는 선형으 로 변동한다고 가정하고 식 (1)의 지배방정식에 mid-point integration 기 법을 적용하면, 그 층에서의 지배방정식의 근사식을 얻을 수 있다.

    L j 4 [ A ˜ A ˜ A ˜ A ˜ ] { 2 u j 1 x 2 2 u j x 2 } + 1 2 { B ˜ B ˜ T B ˜ B ˜ T B ˜ + B ˜ T B ˜ + B ˜ T } { 2 u j 1 x 2 2 u j x 2 } + 1 L j [ G ˜ -G ˜ G ˜ G ˜ ] { u j 1 u j } + L j 4 [ M ˜ M ˜ M ˜ M ˜ ] { 2 u j 1 t 2 2 u j t 2 } + L j 4 [ C ˜ C ˜ C ˜ C ˜ ] { u j 1 t u j t } = { σ j 1 σ j }
    (6)

    여기서, u j ( x , t ) = u ( x , z j , t ) 이고 σ j ( x , t ) = σ ( x , z j , t ) 이다. 식 (6)의 두 열 을 더한 후, 이를 Lj 로 나누면 다음의 식을 얻을 수 있다.

    σ j σ j 1 L j = A ˜ 1 2 ( 2 u j x 2 + 2 u j 1 x 2 ) B ˜ 1 L j ( u j x - u j 1 x ) + M ˜ 1 2 ( 2 u j t 2 + 2 u j 1 t 2 ) + C ˜ 1 2 ( u j t + u j 1 t )
    (7)

    식 (7)은 z 방향으로의 식 (1)의 Crank-Nicolson 이산화임을 알 수 있다.

    식 (4)도 다음과 같이 이산화할 수 있다. 식 (6)의 2열로부터 1열을 뺀 후 그 결과를 2로 나누면 다음과 같이 이산화된 식을 얻을 수 있다.

    σ j + σ j 1 2 = B ˜ T 1 2 ( u j x + u j 1 x ) + G ˜ ( u j u j 1 L j )
    (8)

    식 (8)은 z 방향으로의 식 (4)의 Crank-Nicolson 이산화이다. 길이 Lj 는 Padé 근사에 의하여 scalar wave 전파로부터 얻은 것과 동일하게 얻을 수 있다[18]. 얻어진 길이는 z 방향으로의 파수에 독립적이기 때문에, 다공 성 매질에서의 vector wave 전파문제의 mid-point integrated finite element에 대해서도 동일하게 적용할 수 있다.

    2.2. 무한 다공성 탄성 매질의 점성경계

    무한 다공성 탄성 매질의 점성경계는 Lee[11]에 유도되어 있고, 여기서 는 그 결과만을 간략히 서술하고자 한다. Fig. 1b에 보인 다공성 반무한체에 대하여 x 방향으로 일정한 운동을 가정하고 무한의 투수계수(즉, f = 0)를 가정하면, 매질은 다음과 같이 감쇠기로 표현할 수 있다.(10a)(10b)

    { t x ( t ) t z ( t ) p ( t ) } = { τ x z ( t ) σ z ( t ) p ( t ) } = [ C v S 0 0 C v P ] { u ˙ x ( t ) u ˙ z ( t ) w ˙ z ( t ) }
    (9)

    C v S = ( 1 n ) ρ s C S 0
    (10a)

    C v P [ ρ ρ w ρ w ρ w n ] [ 1 1 C w 1 0 C w 2 0 ] [ C P 1 0 0 0 C P 2 0 ] [ 1 1 C w 1 0 C w 2 0 ] 1
    (10b)

    여기서 u ˙ x ( t ) , u ˙ z ( t ) , w ˙ z ( t ) 는 무한 매질 표면에서의 속도, σz(t), τxz(t) 는 무한 매질 표면에서의 응력, p(t)는 무한 매질 표면에서의 간극수압, tx(t), tz(t), 는 무한 매질 표면에서의 표면력이다. C S 0 = μ / ( 1 n ) ρ s f = 0일 때 다공성 탄성 매질의 S파 속도이고, 이 매질의 P1과 P2파 속도 C P 1 0 C P 2 0 는 다음의 고유값 문제의 고유값이다.

    ( [ λ + 2 μ + Q α 2 Q α Q α Q ] ( C P 0 ) 2 [ ρ ρ w ρ w ρ w n ] ) { 1 C w 0 } = { 0 0 }
    (11)

    식 (11)의 고유값 중 더 큰 값이 C P 1 0 이고, 작은 다른 하나는 C P 2 0 이다. 각 각의 고유값에 대응하는 고유벡터는 { 1 C w 1 0 } T { 1 C w 2 0 } T 이다. 점성경 계가 수직으로 입사하는 파에 대해서만 정확한 모형이지만, 경사지게 입사 하는 파에 대해서도 아주 효율적인 역학적 모형임이 이미 증명되었다[22]. 식 (9)의 점성경계도 경사지게 입사하는 다공성 탄성파를 흡수하는 데 효율 적임을 3절의 예제에서 보일 것이다.

    2.3. 무한 다공성 탄성 매질의 모사

    무한 다공성 탄성 매질을 mid-point integrated finite element와 점성 경계를 사용하여 나타내는 방법을 제안하고자 한다. Fig. 1c는 이 방법을 적 용할 수 있는 전형적인 예인 다공성 반무한체를 보여주고 있다. 반무한체의 0≤zLz의 영역을 나타내기 위해 일련의 mid-point integrate finite elements를 사용한다. 이 요소들의 길이는 Astaneh and Guddati[18]에 주 어져 있는데, 이 길이는 복소수이기 때문에 실제 구현할 때는 주의하여야 한 다. 특히, 범용 유한요소 해석 코드에 이 요소를 구현할 때는 이 점을 염두에 두어야 할 것이다. 이 연구에서는 0≤zLz의 영역을 N 개의 세부 영역으 로 나누고, Lz / N의 길이를 가지는 mid-point integrate finite element를 사용하는 것을 제안한다. 이 요소들은 실수의 길이를 가지기 때문에, 통상 적인 유한요소와 같은 방법으로 구현할 수 있다. 그리고, mid-point integrated finite elements로 나타낸 0≤zLz의 유한 영역에 매질의 무한 영역을 나타내기 위하여 식 (9)의 점성경계를 부착한다.

    식 (6)과 (9)로부터 무한 다공성 탄성 매질에 대해 제안된 모형이 질량, 감쇠, 강성의 항으로만 표현이 가능하다는 것을 관찰할 수 있다. 시간에 대 하여 변위를 적분한 항은 나타나지 않는다. 그러므로, 고차 흡수경계조건과 PML의 단점을 극복할 수 있다. 또한, 각 항의 계수가 상수이기 때문에, 시 간영역 운동방정식에서 convolution이 나타나지 않게 되어, 주파수영역 정 식화에 기반한 방법의 단점 또한 극복할 수 있게 된다. 그러므로, mid-point integrated finite element와 점성경계에 기반하여 제안한 모형은 토목 및 기계공학 분야의 적용에 아주 유용한 방법이라고 할 수 있다. Mid-point integrated finite element와 점성경계를 구현할 수 있도록 기존의 해석 코 드를 조금만 수정하면, 무한 다공성 매질에서의 파전파 문제를 풀 수 있다. 그러므로, 제안된 모형은 무한 다공성 탄성 매질에서의 동적 거동의 수치 해 석에 “실용적”인 방법이라고 할 수 있을 것이다.

    제안된 방법에 근거하여, Fig. 2에 보인 다공성 매질에서의 비선형 파전 파 문제와 지반-구조물 상호작용 문제의 지배방정식을 다음과 같이 얻을 수 있다.

    [ M i i M i b 0 M b i M b b + M b b f M b e f 0 M e b f M e e f ] { U ¨ i U ¨ b U ¨ e } + [ C i i C i b 0 C b i C b b + C b b f C b e f 0 C e b f C e e f ] { U ˙ i U ˙ b U ˙ e } + { P i i n t ( U , U ˙ ) P b i n t ( U , U ˙ ) 0 } + [ 0 0 0 0 K b b f K b e f 0 K e b f K e e f ] { U i U b U e } = { P i e x t P b e x t 0 }
    (12)

    여기서 M, C, P int ( U , U ˙ ) 는 각각 일반적인 유한요소에 의해 표현되는 구조 물과 지반 근역의 질량 행렬, 감쇠 행렬, 비선형 내력 벡터를 의미한다. M f , C f , K f 는 mid-point integrated finite element와 점성경계에 의해 나타 내어 지는 지반 원역의 질량, 감쇠, 강성 행렬을 의미한다. U , U ˙ , U ¨ 는 각각 변위, 속도, 가속도 벡터이고, Pext 는 외력 벡터이다. 식 (12)에서 하첨자 i, b, e는 각각 구조물과 지반 근역의 유한요소에만 위치하는 절점, 지반 근역 과 원역의 경계면에 위치하는 절점, 지반 원역에만 위치하는 절점을 의미 한다.

    지진파만 입사하고 구조물과 지반 근역에 외력이 가해지지 않을 때, 지 반-구조물 상호작용계의 운동방정식은 다음과 같다[9], [10], [23].

    [ M i i M i b 0 M b i M b b + M b b f M b e f 0 M e b f M e e f ] { U ¨ i U ¨ b U ¨ e } + [ C i i C i b 0 C b i C b b + C b b f C b e f 0 C e b f C e e f ] { U ˙ i U ˙ b U ˙ e } + { P i i n t ( U , U ˙ ) P b i n t ( U , U ˙ ) 0 } + [ 0 0 0 0 K b b f K b e f 0 K e b f K e e f ] { U i U b U e } = { 0 M b b f U ¨ b * M e b f U ¨ b * } + { 0 C b b f U ˙ b * C e b f U ˙ b * } + { 0 K b b f K e b f U b * } + { 0 P b * 0 }
    (13)

    여기서 U ¨ b * , U ˙ b * , U b * , P b * 는 입사 지진파에 의한 자유장 운동을 나타낸다. 이 물리량들은 1차원 파전파 해석을 통해 쉽게 얻을 수 있다. 식 (13)에서 M e b f U ¨ b * , C e b f U ˙ b * , K e b f U b * 는 절점 b와 이웃한 절점에만 0이 아닌 값을 가진 다[9], [10], [23].

    식 (12)와 (13)을 풀어서, 지반-구조물 상호작용계의 비선형 거동을 얻 을 수 있다. 이 연구에서는 일정 가속도법(constant acceleration method) 을 사용할 것이다[9-11].

    3. 검증과 적용

    Mid-point integrated finite element와 점성경계로 구성된 수치모형의 정확성과 적용성을 검토하고자 한다. 제안된 모형의 반사계수를 조사한 후, 다공성 탄성 매질에서의 파전파 문제에 제안된 수치모형을 적용하고자 한 다. 마지막으로 다공성 지반에서의 비선형 동적 문제에 제안된 방법을 적용 하여 그 활용성을 검토하고자 한다.

    3.1. 반사계수

    다공성 탄성 반무한체에 제안된 수치모형을 적용하였을 때의 반사계수 를 계산하여, 제안된 모형의 오차를 조사하고자 한다. 입사파와 반사파가 exp[i(ωt - kxx)] 의 의존성을 가지고 각각의 진폭이 A = { A 1 A 2 A S } T B = { B 1 B 2 B S } T 로 주어질 때, 반사계수 R = BA-1 를 Lee[11]에서 유 도하였다. 여기서, 하첨자 1, 2, S는 각각 P1, P2, S파의 진폭을 나타낸다.

    고려한 매질의 재료 성질은 다음과 같다. Lamé 상수는 λ = 85.4 MPa이 고 = μ 42.7 MPa이다. Biot 상수는 Q = 1488.56 MPa이고 α = 00.9이다. 밀도는 ρs = 2640 kg/m3이고, ρw = 1000 kg/m3 이다. 공극율은 n = 0.3 이 고, 투수계수는 κ = 1.019368 × 10-8 m3 m3 ․ s/kg이다. 이 경우, f = 0 일 때의 P1, P2, S파의 속도는 각각 C P 1 0 = 880.8 m/s, C P 2 0 = 230.6 m/s, C S 0 = 152 m/s 0 = S C 이 다. 이 예제에서 입사파의 가진진동수 ω는 50 rad/s이다.

    반무한체를 mid-point integrated finite element와 점성경계를 사용하 여 나타낸다. 요소의 길이(또는 두께)는 0.1 λ R 0 인데, 여기서 λ R 0 f = 0일 때의 Rayleigh 파의 파장이다. 반무한체의 동적 강성 행렬의 판별식을 0로 두어, Rayleigh 파의 속도 C R 0 와 파장 λ R 0 를 계산하면 이 예제에서는 각각 141.56 m/s와 17.79 m로 계산된다.

    Fig. 3af = 0 이고 Lz/ λ R 0 = 0, 1, 2, 3, 4 일 때의 반사계수를 보여주 고 있다. 반무한체의 표면에 평행하게 입사하는 P1, P2, S파의 정규화된 파 수 k x C S 0 / ω 는 그림에 보인 바와 같이 각각 C S 0 / C P 1 0 = 0.173 , C S 0 / C P 2 0 = 0.659 , 1.0 이다. 평행하게 입사하는 파에 대해서는 흡수할 성분이 없으므 로, 이 파수에서 반사계수는 정확히 1이고 이는 Fig. 3a에서 관찰할 수 있다. Fig. 3a에서 k x C S 0 / ω 가 이 값들보다 작을 때는 mid-point integrated finite element와 점성경계를 사용하는 경우와 점성경계만 사용한 경우의 반사계수의 차이가 거의 없음을 알 수 있다. 이는 진행파(propagating waves) 는 점성경계에 의해서 흡수됨을 의미한다. 반면에, k x C S 0 / ω 가 0.173, 0.659, 1.0보다 커지게 되면, 사용된 mid-point integrated finite elements의 수가 증가할수록 반사계수는 감소함을 확인할 수 있는데, 이는 소멸파(evanescent waves)의 흡수는 mid-point integrated finite element에 의해서 이루어짐 을 의미한다. 즉, 점성경계와 mid-point integrated finite element는 각각 진행파와 소멸파의 흡수에 기여한다. 또한, Fig. 3a에서 제안된 수치모형의 효율은 사용되는 mid-point integrated finite element의 수가 증가할수록 증가함을 관찰할 수 있다. 이러한 관찰은 mid-point integrated finite element와 PMDL을 사용했던 기존 연구에서도 관찰된 사항이다[11].

    Fig. 3bf = 1/κ = 9.81 × 107 kg(m3⋅s) 일 때의 반사계수를 나타낸다. Fig. 3a와 비교하여 Ri2R2i (i = 1, 2, S)에서 큰 변화를 관찰할 수 있는데, 이는 f ≠ 0일 때 P2파의 속도가 크게 변하기 때문이다[11]. 그럼에도 불구 하고, 제안된 수치 모형은 여전히 입사파를 잘 흡수할 수 있음을 확인할 수 있다.

    3.2 다공성 탄성 지반에서의 파전파

    제안된 수치모형을 사용하여 다공성 탄성 지반에서의 파전파 문제(Fig. 4)를 풀고자 한다. 지반의 깊이 H = 40 m이다. 지표면은 자유 표면이고 불 투수성 강체 기반암을 가진다. 지반은 수평방향으로 무한하다. 지반의 재료 성질은 다음과 같다. Lamé 상수는 λ =1439 Mpa이고 μ = 719.8 Mpa이 다. Biot 상수는 Q = 5678 Mpa이고 α = 0.9이다. 밀도는 = ρs 2000 kg/m3 이고 =1000 ρw kg/m3이다. 공극율은 n = 0.3이고 투수계수는 무한대, 즉, f = 0이다.

    이 예제에서는 식 (13)의 Gaussian 펄스를 초기 조건으로 사용한다.

    u x ( x , z , 0 ) = x 10 exp [ x 2 ( z 30 ) 2 ]
    (13a)
    u z ( x , z , 0 ) = z 10 exp [ x 2 ( z 30 ) 2 ]
    (13b)
    w x ( x , z , 0 ) = w z ( x , z , 0 ) = 0
    (13c)
    u ˙ x ( x , z , 0 ) = u ˙ z ( x , z , 0 ) = w ˙ x ( x , z , 0 ) = w ˙ z ( x , z , 0 ) = 0
    (13d)

    식 (13)의 조건은 x = 0 에 대하여 대칭이기 때문에, x = 0에서 대칭 조건 을 사용하여 x ≥ 0인 지반만 다루고자 한다. 대칭 조건은 식 (14)에 주어져 있다.

    u x ( 0 , z , t ) = 0
    (14a)
    τ x z ( 0 , z , t ) = 0
    (14b)
    w x ( 0 , z , t ) = 0
    (14c)

    Fig. 4에 보인 바와 같이 지반의 근역은 길이가 1인 유한요소를 사용하 여 이산화한다. 지반의 원역은 mid-point integrated finite element와 점 성경게를 사용하여 나타낸다. 고려한 원역의 수평 길이는 Lx이고 midpoint integrated finite elements의 수는 N 이다.

    식 (13)의 초기 조건에 의한 동적 응답을 계산하고, 참조해와 비교한다. 참조해는 유한요소망을 x = 640 m 까지 확장하여 얻었다. Fig. 5는 식 (15) 에 주어진 상대 오차의 변화를 보여준다.

    e ( t ) = 0 40 0 40 | u ( x , z , t ) u r e f ( x , z , t ) | 2 d x d z 0 40 0 40 | u r e f ( x , z , t ) | 2 d x d z
    (15)

    여기서 u r e f ( x , z , t ) 는 참조해의 변위를 의미한다. Fig. 5에서 N = 0 인 경우 는 mid-point integrated finite elements를 사용하지 않고 점성경계만 지 반 근역에 직접 부착하여 계산한 경우이다. Fig. 5에서 길이 Lx와 요소 수 N 가 증가할수록 오차는 줄어드는 것을 확인할 수 있다. 현재의 수치 모형에서 는 Lx = 160 m이고 N = 80일 때 최적의 결과를 얻을 수 있고, 이를 통상적인 유한요소만 사용했을 때의 N = 160과 비교하여보면 충분한 이득이 있음을 확인할 수 있다. Fig. 6Lx = 320 m이고 N = 10일 때 ux (x, z,t) 의 분포를 보여주고 있다. 결과로부터 아주 적은 수의 mid-point integrated finite element를 사용하더라도 제안된 수치 모형을 사용하면 좋은 결과를 얻을 수 있음을 확인할 수 있다.

    이 예제와 같은 경계조건을 가지는 waveguide에서의 파전파 문제에 대 하여 PML은 무한히 증가하는 결과를 주는 것이 이미 증명되었다[16]. Root-Finding Absorbing Boundary Condition[24]을 제외한 현재까지 개발된 탄성파에 대한 고차 흡수경계조건은 안정한 결과를 얻기 위하여 주 기적인 경계조건을 가정하기 때문에, 이 예제와 같은 경계조건을 가지는 waveguide에서의 파전파 문제에 적용할 수 없다[17]. 그러므로, 제안된 접 근법의 안정성을 수치적으로 조사하기 위하여, 동적 응답을 500,000 시간 단계동안 계산하였다. 이 논문에 계산 결과를 포함하지는 않았지만, 계산 결과에서 어떠한 불안정성도 관찰되지 않았다. 점성경계의 안정성에 대해 서는 이미 증명이 된 바이므로[25], 이 연구에서 제안한 모형을 사용하면 파 전파 문제에 대해 안정한 결과를 얻을 수 있다.

    3.3 다공성 지반에 설치된 케이슨 방파제의 지진응답해석

    제안된 수치 모형을 사용하여 다공성 지반에 설치된 케이슨 방파제(Fig. 7)의 지진응답을 얻고자 한다. 구조적 변형율은 아주 작고 이 해석에서 주요 관심사는 다공성 지반에서의 소성 거동이므로, 구조물은 간단히 강체 질량 으로 모사한다. 구조물의 질량은 2500 kg/m3이다. 일반적으로 케이슨 방 파제는 사석 위에 설치되지만, 이 해석에서는 사석은 별도로 고려하지 않는 다. 해수의 깊이는 11 m이고 밀도는 1000 kg/m3이다. 구조물 양면에 작용 하는 동수압력의 효과는 Westergaard의 부가질량[26]으로 모사한다. 이 예제에서 양면에 위치한 질량의 크기는 7.0583×104 kg로 계산되었다.

    다공성 지반은 강체 기반암을 가지는 것으로 가정한다. 지반의 두께는 40 m이다. 지반의 재료 성질은 다음과 같다. Lamé 상수는 λ =1439 MPa 이고 μ = 719.8 MPa이다. Biot 상수는 Q = 5678 MPa 이고 α = 0.9이다. 밀도는 = ρs 2000 kg/m3이고 = ρw 1000 kg/m3이다. 공극율은 n = 0.3 이고 투수계수는 κ = 1.019368×10-8 m3 ․ s/kg이다. 지반의 비선형 재료 거 동은 Drucker-Prager 소성 모형을 사용하여 모사한다. 재료의 마찰각과 점 착력은 각각 20° 와 0.219 MPa이다. 지반의 감쇠는 2ξs/ωs = 0.003551 sec의 상수를 가지는 강성 비례 감쇠를 가정한다. 지반의 기본진동수는 ωs = πCS / 2H = 28.16 rad/sec이므로, 전술한 감쇠는 이 기본진동모드에 대하여 5% 감쇠(ξs = 0.05)를 부여한다.

    Fig. 7은 이 문제의 유한요소망을 보여주고 있다. 지반 근역은 통상적인 유한요소에 의해서 나타내어 진다. 지반 원역은 mid-point integrated finite element와 점성경계에 의해서 나타내어진다. 수치 모형에서 길이 Lx 와 고려한 mid-point integrated finite elements의 수 N 은 각각 320 m와 10이다.

    동적 하중이 시스템에 가해지기 전에 정적 조건을 먼저 고려한다. 정적 하중은 구조물과 지반의 자중과 유체의 정수압을 포함한다. 점성경계 없이 mid-point integrated finite element의 강성만을 고려하여 정적 상태를 결 정할 수 있다[27]. 그 후, El Centro 지반운동(Fig. 8)을 자유장 지반의 지표 면 운동으로 가하고 지배 방정식 (13)을 풀어서 시스템의 지진 응답을 얻는 다. 해를 얻기 위하여 일정 가속도법(constant acceleration method)에 기 반한 해법을 사용한다.

    Fig. 9는 계산이 종료한 후의 등가 소성 변형율[28]의 분포를 보여준다. Fig. 10은 구조물의 회전과 Point A와 B에서의 평균 유효 응력의 시간 이력 을 보여준다. 이전의 연구[9-11]와 유사한 비선형 지반-구조물 상호작용의 결과를 확인할 수 있다. 응력 집중으로 인하여 구조물 모서리 근방에서 소성 응답이 발생한다. 지반의 부등 소성 변형으로 인하여 비선형 해석이 종료된 후 구조물의 영구 회전이 존재하게 된다. 이러한 영구 회전은 선형 해석에서 는 관찰할 수 없다. 영구 회전은 모서리와 가까운 좌우 영역에 서로 다른 소 성 응답을 초래하게 된다. 구조물의 시계방향 회전으로 인하여, Point A에 서의 평균 유효 응력이 증가하고 Point B에서의 압축이 심화되며, 결과적 으로 인장에 취약한 지반 재료의 특성으로 인해 Point A에서의 응력의 변화 가 Point B에서의 변화보다 더 크게 발생한다. 비선형 해석으로부터 얻은 지진응답은 선형 해석의 결과와 다른 것을 관찰할 수 있다. Mid-point integrated finite element와 점성경계를 사용하여 다공성 지반에서의 비선 형 지반-구조물 상호작용 해석을 성공적으로 수행할 수 있음을 이 예제에서 확인할 수 있다.

    4. 결 론

    이 연구에서는 무한 다공성 매질의 시간영역 파전파 해석에 사용될 수 있는 정확하고 효율적이며 구현하기 쉬운 수치적 접근법을 제안하였다. 제 안된 접근법은 다공성 탄성 매질의 mid-point integrated finite element와 점성경계에 근거하고 있다. 제안된 모형의 동적 거동은 질량, 감쇠, 강성 행 렬에 의해서만 표현되므로, 변위 기반 유한요소 정식화를 사용하여 쉽게 구 현할 수 있다. 또한 계수 행렬이 상수이므로, 시간 영역 해석을 위해 convolution 연산을 필요로 하지 않는다. 그러므로, 기존 수치해석 코드를 조금만 수정하여 제안된 모형을 쉽게 구현할 수 있다.

    제안된 수치 모형은 다공성 지반에서의 파전파 해석에 적용하여 제안된 수치 접근법의 정확성과 안정성을 검증하였다. 또한, 다공성 지반에 놓인 케이슨 방파제의 지진응답해석에 적용하여, 비선형 지반-구조물 상호작용 해석에의 활용 가능성을 검토하였다. 적용 예제로부터 변위기반 유한요소 정식화에 의해 쉽게 구현할 수 있는 제안된 수치 모형을 사용하면 다공성 지 반에서의 비선형 파전파 문제와 지반-구조물 상호작용 문제에 대하여 정확 하고 안정한 결과를 얻을 수 있음을 결론지을 수 있다.

    / 감사의 글 /

    이 논문은 2016학년도 부경대학교의 지원을 받아 수행된 연구임 (C-D- 2017-0052).

    Figure

    EESK-22-379_F1.gif

    Representation of a poroelastic half-space

    EESK-22-379_F2.gif

    Soil-structure interaction system in a poroelastic half-space

    EESK-22-379_F3.gif

    Reflection coefficients

    EESK-22-379_F4.gif

    Wave-propagation problem in poroelastic soil

    EESK-22-379_F5.gif

    Variations of L2 norms of errors

    EESK-22-379_F6.gif

    Wave fields when Lx = 320 m and N = 10

    EESK-22-379_F7.gif

    Caisson breakwater on poroelastic soil

    EESK-22-379_F8.gif

    Input free-field motion (1940 El Centro earthquake, N-S component, PGA=0.319 g)

    EESK-22-379_F9.gif

    Equivalent plastic strains in the soil

    EESK-22-379_F10.gif

    Earthquake responses in the soil

    Table

    Reference

    1. SchanzM. Poroelastodynamics: Linear models, analytical solutions, and numerical methods . Applied Mechanics Review.2009;62.
    2. KauselE. Forced vibrations of circular foundations on layered media, Research Report R74-11. Department of Civil Engineering, Massachusetts Institute of Technology, Cambridge, Massachusetts,c1974.
    3. BeskosDE . Boundary element methods in dynamic analysis . Applied Mechanics Review.1987;40:1-23.
    4. BeskosDE . Boundary element methods in dynamic analysis: Part II (1986-1996) . Applied Mechanics Review.1998;50:149-197.
    5. AstleyRJ . Infinite elements for wave propagation: a review of current formulations and an assessment of accuracy . International Journal for Numerical Methods in Engineering.2000;49:951-976.
    6. GivoliD. High-order local non-reflecting boundary conditions: a review . Wave Motion.2004;39:319-326.
    7. BasuU , ChopraAK . Perfectly matched layers for time-harmonic elastodynamics of unbounded domains: theory and finite-element implementation . Computer Methods in Applied Mechanics in Engineering.2003;192:1337-1375.
    8. BasuU , ChopraAK . Perfectly matched layers for transient elastodynamics of unbounded domains . International Journal for Numerical Methods in Engineering.2004;59:1039-1074.
    9. LeeJH , KimJK , KimJH . Nonlinear analysis of soil structure interaction using perfectly matched discrete layers. Computers andStructures. 2014;142:28-44.
    10. LeeJH , KimJH , KimJK . Perfectly Matched Discrete Layers for Three-Dimensional Nonlinear Soil-Structure Interaction Analyses .Computers and Structures.2016;165:34-47.
    11. LeeJH . Nonlinear Soil-Structure Interaction Analysis in Poroelastic Soil Using Mid-Point Integrated Finite Elements and Perfectly MatchedDiscrete Layers . Soil Dynamics and Earthquake Engineering.2018; 108:160-176.
    12. Meza-FajardoKC , PapageorgiouAS . Study of the accuracy of the multiaxial perfectly matched layer for the elastic-wave equation . Bulletin of the Seismological Society of America.2012;102:24582467.
    13. GuddatiMN , TassoulasJL . Continued-fraction absorbing boundary conditions for the wave equation . Journal of Computational Acoustics.2000;8:139-156.
    14. GuddatiMN , LimK-W . Continued fraction absorbing boundary conditions for convex polygon domains . International Journal of Numerical Methods in Engineering.2006;66:949-977.
    15. SavadattiS , GuddatiMN . Absorbing boundary conditions for scalar waves in anisotropic media. Part 2: Time-dependent modeling . Journal of Computational Physics.2010;229:6644-6662.
    16. DuruK , KreissG . Numerical interaction of boundary waves with perfectly matched layers in two space dimensional elastic waveguides . Wave Motion.2014;51:445-465.
    17. BaffetD , BielakJ , GivoliD , HagstromT , RabinovichD . Long-time stable high-order absorbing boundary conditions for elastodynamics . Computer Methods in Applied Mechanics and Engineering.2012; 241-244:20-37.
    18. AstanehAV , GuddatiMN . Efficient computation of dispersion curves for multilayered waveguides and half-spaces . Computer Methods in Applied Mechanics and Engineering.2016; 200: 27-46.
    19. BiotMA . Theory of propagation of elastic waves in a fluidsaturated porous solid. I. Low-frequency range . Journal of Acoustical Society of America.1956;28:168178.
    20. ZienkiewiczOC , ChanAHC , PastorM , SchreflerBA , ShiomiT . Computational Geomechanics with Special Reference to Earthquake Engineering. Wiley. c1999.
    21. LewisRW , SchreflerBA . The Finite Element Method in the Static and Dynamic Deformation and Consolidation of Porous Media. Wiley. c1998.
    22. LysmerJ , KuhlemeyerRL . Finite dynamic model for infinite media . ASCE. Journal of the Engineering Mechanics Division.1969;95: 859-877.
    23. GuddatiM , SavadattiS . Efficient and accurate domain-truncation techniques for seismic soil-structure interaction . Earthquakes and Structures.2012; 3: 563-580.
    24. LeeJH , TassoulasJL . Root-Finding Absorbing Boundary Conditions for Problems of Wave Propagation in Infinite Media. Computer Methods in Applied Mechanics and Engineering. c2018.
    25. RabinovichD , GivoliD , BielakJ , HagstromT , The Double Absorbing Boundary method for a class of anisotropic elastic media . ComputerMethods in Applied Mechanics and Engineering.2017; 315:190-221.
    26. WestergaardHM . Water Pressures on Dams during Earthquakes. Transaction ASCE. 1933;98:418-433.
    27. SavadattiS , GuddatiMN . A finite element alternative to infinite elements . Computer Methods in Applied Mechanics and Engineering.2010;199:2204-2223.
    28. ChenW-F . Constitutive Equations for Engineering Materials, Volume 2: Plasticity and Modeling. Elsevier. c1994.