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.29 No.6 pp.583-590

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

A Study of Rayleigh Damping Effect on Dynamic Crack Propagation Analysis using MLS Difference Method

Kyeong-Hwan Kim1 , Sang-Ho Lee1 , Young-Cheol Yoon2
1Department of Civil and Environment Engineering, Yonsei University, Seoul, 06285, Korea
2Department of Civil Engineering, Myongji College, Seoul, 03656, Korea
Corresponding author: +82-2-300-1135 (ycyoon@mjc.ac.kr
November 8, 2016 November 21, 2016 November 22, 2016

Abstract

This paper presents a dynamic crack propagation algorithm with Rayleigh damping effect based on the MLS(Moving Least Squares) Difference Method. Dynamic equilibrium equation and constitutive equation are derived by considering Rayliegh damping and governing equations are discretized by the MLS derivative approximation; the proportional damping, which has not been properly treated in the conventional strong formulations, was implemented in both the equilibrium equation and constitutive equation. Dynamic equilibrium equation including time relevant terms is integrated by the Central Difference Method and the discrete equations are simplified by lagging the velocity one step behind. A geometrical feature of crack is modeled by imposing the traction-free condition onto the nodes placed at crack surfaces and the effect of movement and addition of the nodes at every time step due to crack growth is appropriately reflected on the construction of total system. The robustness of the proposed numerical algorithm was proved by simulating single and multiple crack growth problems and the effect of proportional damping on the dynamic crack propagation analysis was effectively demonstrated.


MLS 차분법을 활용한 동적 균열전파해석의 Rayleigh 감쇠영향 분석

김 경 환1 , 이 상 호1 , 윤 영 철2
1연세대학교 토목환경공학과,
2명지전문대학 토목과

초록

본 논문은 강형식 기반의 MLS 차분법에 Rayleigh 감쇠효과를 적용한 동적균열진전 해석기법을 제시한다. Rayleigh 감쇠 효과가 반영된 동적 평형방정식과 구성방정식을 도출하고, MLS 미분근사식을 이용하여 지배방정식들을 이산화하였다. 평형 방정식뿐만 아니라 구성방정식에서도 감쇠효과를 적절하게 고려하여 기존의 무요소 강정식화 기법에서 고려하지 못했던 비 례감쇠 알고리즘을 구현하였다. 시간관련 항을 포함한 동적 평형방정식은 중앙차분법(central difference method)을 이용하 여 시간적분 하였고, 속도에 대한 차분식을 lagging시켜 이산화 방정식을 간소화시켰다. 균열의 기하학적 특성은 표면력 ‘0’ 인 자연경계 조건을 균열면에 놓인 절점들에 부과하여 묘사하였으며, 균열성장으로 인해 해석단계마다 변하는 절점의 생성 및 이동 효과를 계방정식 구성에 반영하였다. 단일균열과 다중균열을 갖는 수치예제를 통해서 제안된 수치기법의 정확성을 검증하였으며, 비례감쇠 효과의 고려가 동적균열진전 해석결과에 미치는 영향을 보였다.


    Ministry of Science, ICT and Future Planning

    National Research Foundation
    No. 2014R1A1A1002000

    Ministry of Land, Infrastructure and Transport
    16RTRP-B104237-02

    1.서 론

    수치해석은 컴퓨터를 활용해 주어진 지배방정식에 대한 수치해을 구하기 때문에 오차를 수반한다. 고체의 동적문제를 기술하는 미분방정식에 대한 시간해석 역시 해석단계마다 오차가 발생되고 또 중첩된다. 이와 같은 현상은 동적 균열 문제와 같이 hyperbolic 방정식의 특성을 나타내는 문제에서 더욱 두드러지게 나타나는데, 감쇠 알고리즘의 도입은 이와 같은 오차 발생을 누그러뜨리고 시간해석의 안정성 확보에 도움을 준다. 감쇠 알고리즘은 재료 자체가 갖고 있는 감쇠 특성을 모사하는 역할도 하기 때문에 재료 및 시스템의 감쇠 특성에 대한 고려가 필요하다. 또한, 외부 하중에 의한 동적 효과가 큰 내진해석에서도 감쇠 알고리즘이 중요한 역할을 한다(Omidi et al., 2013; Hall, 2006).

    1990년대 중반 Belytschko 등(1994)Liu 등(1995)에 의해서 Element-Free Galerkin(EFG)법 및 재생커널 점별법 (reproducing kernel particle method)과 같은 무요소법이 등장한 이래 다른 한편으로 수치적분을 탈피하고 강형식 형태의 지배방정식을 직접 이산화하는 수치해석 방법들도 나타나기 시작했다. 초기에 Atluri와 Zhu(1998)Gu와 Liu(2001)가 제안한 수치방법과 같이 부분적으로 약형식 지배방정식을 사용하는 경우에는 기존의 유한요소법과 유사한 방법으로 감쇠 알고리즘을 적용할 수 있었으나, 이 경우 약형식에 따른 요소망 구성과 적분에 제약으로 무요소법의 장점을 효과적으로 활용할 수 없었다. 이러한 문제를 극복하기 위하여 약형식을 배제하고 강형식 지배방정식을 바로 이산화한 무요소 점별 콜로케이션법 (meshfree point collocation method) 형태의 수치방법 들이 등장하였다(Kim and Kim, 2003; Lee and Yoon, 2004). 콜로케이션 기반의 무요소 기법들은 절점만을 사용 하여 고체문제를 해석하기 때문에 해석 형상이 복잡하거나 시간에 따라 해석 대상이 변화하는 다양한 문제에 적용할 수 있었다(Yoon et al., 2006; Wen and Aliabadi, 2007). 그러나 온전히 강형식만 사용하면서 Rayleigh 감쇠와 같은 비례감쇠 알고리즘을 적용하기 어려워졌다.

    기존의 유한요소법에서 Rayleigh 감쇠 알고리즘은 행렬식 형태로 표현된 운동방정식에 존재하는 강성행렬과 질량행렬에 대한 적절한 비율로 감쇠행렬을 정의하기 때문에 감쇠효과를 손쉽게 반영할 수 있고, 감쇠상수나 감쇠비율의 물리적 의미가 비교적 명확하여 매우 편리하고 효과적인 알고리즘이다. 그러나 강정식화의 경우 약정식화와 같은 강성행렬이나 질량행렬이 별도로 정의되지 않기 때문에 기존의 방식으로 Rayleigh 감쇠 알고리즘을 사용할 수 없다. 그래서 대부분의 강정식화 기반의 동적해석 기법 연구들은 운동방정식에 속도항과 감쇠상수를 추가하는 방법으로 감쇠효과를 고려해 왔다(Sladek et al., 2003; Gu and Liu, 2005; Long et al., 2006). 그러나 이런 방법으로는 강정식화에서 해석대상의 질량에 따른 감쇠영향 밖에 고려할 수 없기 때문에 Rayleigh 감쇠와 같은 비례감쇠 효과를 제대로 반영할 수 없고, 비교적 감쇠상수와 감쇠비의 관계가 명확한 Rayleigh 감쇠 알고리즘을 활용하기 힘들어 진다.

    본 논문에서는 무요소 강정식화 기반의 MLS 차분법에 Rayleigh 감쇠 알고리즘을 적용하여 동적 균열 전파해석을 수행하였다. MLS 차분법은 Taylor 전개를 기반으로 이동 최소제곱법(moving least squares method)을 통해 절점 모델의 미분근사식을 도출하고 지배방정식을 직접 이산화한다. MLS 차분법은 해석과정 중에 해석대상의 기하학적 형상이 변해도 요소망의 재구성없이 절점들의 관계를 통해 해석을 수행할 수 있기 때문에 균열이 진전하면서 수반되는 기하학적 변화를 손쉽게 모사할 수 있다. 또한 Rayleigh 감쇠 알고리즘의 적용을 통해 기하학적 형상변화를 유발하는 동적균열 해석의 안정성 확보와 함께 재료의 감쇠특성도 반영할 수 있다. 기존의 동적 균열전파 해석의 경우 수치기법 자체의 topology 특성 이나 균열의 불연속성 모델링에 의해 떨림현상(oscillation)이 심각하게 발생하여 해석의 안정성이 심각하게 훼손되는 경우가 많았는데, 본 연구에서는 이와 같은 현상이 효과적으로 억제되었다. 또한, 기존의 무요소 강정식화 방법들에서 사용한 감쇠 알고리즘 보다 정확성 높은 해석이 가능하며, 기존의 유한요소법에서 사용했던 감쇠 관련된 기법들을 무요소 강정식화에서도 그대로 사용할 수 있도록 해준다. 동적하중을 받는 균열전파 문제에 대한 수치예제를 통해 새롭게 제안된 해석기법의 정확성을 검증한다.

    2.MLS 차분법을 위한 감쇠 알고리즘

    본 장에서는 강정식화 과정에서 Rayleigh 감쇠를 반영하기 위해 동적 평형방정식과 구성방정식을 수정하는 방법을 소개한 후, MLS 차분법을 이용한 지배방정식의 이산화 과정을 소개 한다. 이산화 과정에서 동적해석을 위한 시간적분은 중앙차분법 (central difference method)을 사용하였다.

    2.1.비례감쇠 효과를 포함한 지배방정식의 이산화

    기존의 무요소 강정식화 기반의 고체역학 문제해석에서는 지배방정식에 감쇠효과를 고려하기 위해 미소변형을 가정한 동적 평형방정식에 다음과 같이 속도항을 추가하여 사용했다.

    σ + b = ρ a + c vin  Ω
    (1)

    여기서, σ는 응력텐서, b는 체적력, ρ는 밀도, a는 가속도, c는 감쇠상수, v는 속도이며, Ω는 내부영역이다. 식 (1)과 같이 감쇠효과를 고려하면 비례감쇠를 기준으로 볼 때 질량에 비례하는 감쇠효과는 고려할 수 있으나 강성에 대한 감쇠 효과는 고려할 수 없다. 본 연구에서는 강정식화 과정에 Rayleigh 감쇠 알고리즘을 적용하여 질량과 강성에 대한 감쇠 효과를 모두 고려하고자 한다. Hughes(1987)의 지적에서도 알 수 있듯이, 아래와 같은 두 식을 동시에 고려해야 강정식화 에서 Rayleigh 감쇠를 올바로 적용할 수 있다.

    σ + b = ρ a + ρ η 1 vin  Ω
    (2)
    σ = 2 μ s u + λ tr s u 1 + η 2 2 μ s v + λ tr s v 1
    (3)

    여기서, η1η2는 Rayleigh 감쇠상수이고, λμ는 Lamé 상수이다. u는 미지의 변위벡터를 나타내며, ∇ su = (∇ u+∇ uT)/2, ∇ sv = (∇v+∇vT)이다. η1η2는 각각 질량과 강성에 대한 감쇠 영향을 반영한다.

    자연경계(Γt)와 필수경계(Γu)에 대한 지배방정식은 각각 다음과 같다.

    σ n = t ¯ on   Γ t
    (4)
    u = u ¯ on   Γ u
    (5)

    여기서, n과 t는 단위수직벡터와 미리 정의된 표면력 벡터를, u는 미리 규정된 변위값을 나타낸다. 자연경계조건을 나타내는 식 (4)는 응력텐서를 포함하고 있기 때문에 식 (3)을 대입하면 변위에 대해 표현할 수 있다. 단, 약정식화와 달리 강정식화 에서는 내부영역과 경계영역에 대해 서로 다른 근사함수를 적용하는 것도 가능하다(Yoon and Song, 2014).

    지배방정식들을 이산화한 후 하나의 시스템으로 조립하기 위하여 편의상 내부영역 및 자연경계에 대한 지배방정식을 변위에 대하여 기술하는 것이 유리하다. 식 (2)와 (3)을 바탕으로 체적력을 무시하면, 아래와 같이 수정된 Navier 방정식이 유도된다.

    μ 2 u + λ + μ u + η 2 μ 2 v + λ + μ v = ρ a + ρ η 1 v
    (6)

    자연경계조건에 대한 식 (4)도 변위에 대한 다음과 같이 표현할 수 있다.

    μ n s u + λ n 1 u = t ¯
    (7)

    수치적인 편의성을 위해 미분연산자 L과 B를 정의하면, 식 (6)과 (7)을 아래와 같이 쓸 수 있다(Yoon et al., 2014).

    Lu + η 2 Lv = ρ Ia + ρ η 1 Iv  in  Ω
    (8)
    Bu + t ¯ on  Γ t
    (9)

    여기서, L = [μ2+(λ+μ)∇∇·], B = [μu·∇ s + λn·1(∇·)]이며, I는 MLS 차분법에서 0차 미분근사에 대한 형상함수로 구성된 행렬이다. 2차원 문제의 경우, α = (α1,α2) 이므로 Yoon 등(2014)이 제시한 미분근사식에 대한 형상 함수 ϕαI를 이용하여 절점 xJ에서 구성한 L(xJ), B(xJ), I(xJ)를 다음과 같이 정의할 수 있다.

    L x J = I = 1 N L 11 I L 12 I L 21 I L 22 I
    (10a)
    B x J = I = 1 N B 11 I B 12 I B 21 I B 22 I
    (10b)
    I x J = I = 1 N I 11 I I 12 I I 21 I I 22 I
    (10c)

    위의 미분연산자에 대한 행렬식의 구성요소는 다음과 같이 주어진다.

    L 11 I x J = λ + 2 μ φ I 2 , 0 x J + μ φ I 0 , 2 x J
    (11a)
    L 22 I x J = μ φ I 2 , 0 x J + λ + 2 μ φ I 0 , 2 x J
    (11b)
    L 12 I x J = L 21 I x J = λ + μ φ I 1 , 1 x J
    (11c)
    B 11 I x J = λ + 2 μ φ I 1 , 0 x J n 1 + μ φ I 0 , 1 x J n 2
    (12a)
    B 22 I x J = μ φ I 1 , 0 x J n 1 + λ + 2 μ φ I 0 , 1 x J n 2
    (12b)
    B 12 I x J = B 21 I x J = μ φ I 1 , 0 x J n 2 + λ φ I 0 , 1 x J n 1
    (12c)
    I 11 I x J = I 22 I x J = φ I 0 , 0 x J
    (13a)
    I 12 I x J = I 21 I x J = 0
    (13b)

    미분근사식의 형상함수와 미분연산자에 대한 좀 더 자세한 내용은 Lee 등(2016), Yoon 등(2014), Yoon과 Song (2014)을 참고할 수 있다.

    2.2.중앙차분법에 의한 시간적분

    앞 절에서는 전체 계 방정식을 구성하기 위해 지배방정식을 변위로 표현하고, MLS 차분법의 미분근사식을 이용하여 지배방정식을 이산화하는 과정을 설명하였다. 본 절은 이산화 방정식을 시간적분하는 과정을 제시한다. 식 (8)에서 보듯이 동적 평형방정식은 속도와 가속도 항을 갖고 있으므로 시간 적분을 통해 변위만의 방정식으로 변환할 필요가 있다. 중앙 차분법을 적용하여 시간적분을 하면, n 시간단계의 속도와 가속도는 다음과 같이 주어진다.

    v n = u n + 1 u n 1 2 Δ t
    (14a)
    a n = u n + 1 2 u n + u n 1 Δ t 2
    (14b)

    여기서, 위첨자 n은 시간단계이고, Δt는 시간간격이다. 식 (14a)~(14b)를 식 (8)에 대입하면, 시간항을 포함한 동적 평형방정식이 다음과 같이 정리된다.

    ρ 1 + η 1 Δ t 2 I n + 1 η 2 Δ t 2 L n + 1 u n + 1 = 2 ρ I n + Δ t 2 L n u n + ρ η 1 Δ t 2 1 I n 1 + η 2 Δ t 2 L n 1 u n 1
    (15)

    식 (15)에서 보듯이 속도와 가속도를 변위로 치환하는 과정에서 미분연산자 L과 I가 시간단계에 따라 구분되었다. 이것은 균열이 성장함에 따라 미분연산자를 구성하는 미분 근사식이 시간단계 별로 변하는 것을 반영한다. 시간해석을 수행하더라도 해석대상의 기하형상이 변하지 않으면 미분 근사식과 미분연산자가 변하지 않지만, 균열이 진전되면 균열 주변의 절점 배치가 바뀌기 때문에 이에 대한 고려가 뒤따라야 한다. MLS 차분법은 Fig. 1과 같이 진전하는 균열을 묘사하기 위해 절점을 이동시키거나 추가할 수 있다. 이런 절점의 변화를 반영하기 위해 시간단계가 구분된 미분연산자를 매 시간단계 마다 갱신(update)해야 한다.

    기존의 유한요소법에서 중앙차분법은 explicit 시간적분법 으로서 적용이 간편하고 n+1 단계 해석시 역행렬 계산이 필요 없다. 그러나 식 (15)는 중앙차분법을 사용하였음에도 불구하고 n+1 단계의 변위를 구하기 위해 역행렬 계산이 필요하다. 이 점을 개선하기 위해 식 (14a)의 속도 점화식을 다음과 같이 한 단계 lagging 시킬 수 있다.

    v n = u n u n 1 Δ t
    (16)

    식 (16)을 식 (8)에 대입하여 정리하면 다음과 같이 좀 더 간략화된 수식이 얻어진다.

    ρ I n + 1 u n + 1 = ρ 2 Δ t η 1 I n + Δ t Δ t + η 2 L n u n + ρ Δ t η 1 1 I n 1 + Δ t η 2 L n 1 u n 1
    (17)

    식 (17)을 식 (15)와 비교하면, 밀도는 상수이고 미분 연산자 I는 대각행렬에 가깝기 때문에 식 (17)의 좌변에서 n+1 단계의 변위에 곱해져 있는 행렬이 훨씬 단순해 진다. 그러나 MLS 차분법의 특성상 영향영역에 포함되는 절점이 유한요소처럼 결정론적으로 정해지지 않고 또 계 방정식 구성시 내부영역뿐만 아니라 경계조건식도 함께 조합되기 때문에 미지변위를 구하기 위한 전체 강성행렬이 대각선행렬로 정리될 수는 없다. 한편, 감쇠상수 값이 모두 0이 되면, 식 (15)와 식 (17)은 같아진다는 것으로부터 감쇠상수가 동적 평형방정식의 속도항에만 영향을 준다는 점을 확인할 수 있다.

    이제 내부영역과 경계에 대한 지배방정식을 모두 조합하여 다음과 같은 전체 계 방정식으로 구성할 수 있다.

    K n + 1 U n + 1 = F n + 1
    (18)

    여기서, Kn + 1는 내부절점, 경계절점에 대한 이산화 방정식 (또는 차분방정식)을 모두 포함한 일반화된 강성행렬, Un + 1는 변위행렬, Fn + 1는 일반화된 외력행렬을 나타낸다. 2차원 해석의 경우 강성행렬의 크기는 전체 절점의 수가 Nt 일 때 (2 Nt×2 Nt)이다. 식 (18)을 임의의 절점 xJ에 대한 차분 방정식으로 표현하면 다음과 같다.

    I = 1 N K 11 I n + 1 x J K 12 I n + 1 x J K 21 I n + 1 x J K 22 I n + 1 x J u 1 I n + 1 u 2 I n + 1 = F 1 n + 1 x J F 2 n + 1 x J
    (19)

    모든 절점(xJ)에 대한 차분방정식을 조립하면 식 (18)의 계 방정식을 구성할 수 있다. Kn + 1와 Fn + 1의 구성성분들은 절점 xJ가 속한 영역에 따라 달라진다. 즉, xJ가 필수경계에 있다면 Kn + 1(xJ)=In + 1(xJ), Fn + 1(xJ)=u(xJ), 그리고 xJ가 자연경계에 있다면 Kn + 1(xJ)=B(xJ), Fn + 1(xJ)=t(xJ)이다. xJ가 내부영역에 존재하면 Kn + 1(xJ)와 Fn + 1(xJ)는 식 (15) 또는 식 (17)에서 정의된 것과 같이 현재 또는 이전 시간단계 미분연산자 및 변위값들로 구성된다.

    3.수치예제

    이 장에서는 단일균열과 다중균열을 갖는 시편에 대한 동적해석 결과를 이론값과 비교하여 MLS 차분법의 동적균열 해석에 대한 정확성과 감쇠 알고리즘의 효과를 살펴본다.

    3.1.단일균열 전파문제

    본 절에서는 단일균열이 있는 문제에 대한 MLS 차분법 균열진전해석 결과를 이론값(Freund, 1990)과 비교한다. 해석 대상은 Fig. 2와 같은 형상을 가지며 윗면에 인장응력을 받고 있다. 평면변형률 상태를 가정하였고, 재료의 물성치는 밀도 ρ=7,800kg/m3, 탄성계수 E=211×109N/m2, 포아 송비 ν=0.30를 적용하였다. 초기균열 길이는 5m이며, 초기 응력 σ0=1000N/m2는 시간에 관계없이 일정하게 적용하였다. 해석모델은 1071개(51×21)의 절점을 배치한 후 초기 균열을 묘사하기 위해 25개의 절점을 추가하여 총 1096개의 절점을 포함한다. 시간간격은 Δt=0.5×10-5s이고, 200단계 해석을 수행하였다. 또한, 균열모드는 mode 1만 고려하였다.

    Fig. 3은 감쇠효과 반영의 유무에 따른 MLS 차분법 해석결과와 이론값을 비교하여 보여준다. 균열속도 υc가 0인 경우는 균열이 진전하지 않으며, υc가 0.4cs인 경우는 동적 응력확대계수(또는 동적 에너지 해방률) 값이 재료의 동적 파괴인성치를 초과하는 시점부터 균열이 진전되기 시작한다. cs는 전단파(shear wave) 속도이다. 식 (17)과 같이 간소화 된 수식을 사용한 결과와 원래의 식 (15)를 사용한 경우를 함께 도시하였다. MLS 차분법으로 계산한 결과가 이론값과 대체로 잘 일치하고 있으며, 감쇠효과를 적용한 경우, 적용 하지 않은 경우 보다 수치해의 떨림현상이 감소된 것을 볼 수 있다. 감쇠효과를 고려한 경우에는 식 (17)를 사용한 경우와 식 (15)를 사용한 경우의 해석결과가 큰 차이를 보이지 않고 비슷하게 나타났다.

    3.2.다중균열 전파문제

    다중균열의 전파해석에 대한 MLS 차분법의 성능 및 감쇠 알고리즘 영향을 조사하기 위해 Fig. 4와 같이 두 개의 균열이 포함된 문제(John and Shah, 1990)를 해석하였다. mode 1과 mode 2가 결합된 다중 균열모드를 고려하였으며, 초기 균열 길이는 c1=0.01905m, c2=0.003175m로 가정하였다. 속도 형태로 주어진 외력은 196μs까지 일정한 비율로 0.06 m/s에 도달하게 한 뒤, 일정하게 유지시켰다. 해석대상의 재료 물성치는 밀도 ρ=32,400kg/m3, 탄성계수 E=31.37×109 N/m2, 포아송비 ν=0.20, 동적 파괴인성치 KIc=0.8Mpa m 가정하였다. 균열속도는 υc=0.05cs로 가정하였고, 시간간격은 Δt=0.7×10-6s로 설정하였다. 해석모델 구성시 1825(73×25)개의 절점을 배치한 후, 균열 묘사를 위해 16개의 절점을 추가하여 총 1841개의 절점을 배치하였다. Crack 1과 crack 2의 간격L에 따라 균열진전 양상이 다르게 나타나는데, L을 0.0778m와 0.0714m로 가정한 후 해석을 수행하였다.

    Fig. 5L이 0.0778m인 경우 얻어진 균열진전 경로를 보여준다. Crack 1은 성장하지 않는 반면 crack 2는 외력이 가해지는 지점을 향하여 진전하고 있는 것을 볼 수 있다. 대체로 감쇠효과의 유무와 상관없이 일관된 균열진전 경로를 나타내고 있다. 감쇠 알고리즘의 영향을 좀 더 살펴보기 위해 Fig. 6과 Fig. 7에 동적 응력확대계수 값을 도시하였다. 이에 대한 정해가 없기 때문에 감쇠효과를 고려하지 않은 경우를 함께 도시하였다. 동적 응력확대계수는 동적 파괴인성치로 정규화하였고, 모드 1과 모드 2 동적 응력확대계수를 조합한 값이 1을 넘으면 균열이 진전되기 시작한다. 이 경우 L이 상대적으로 크기 때문에 응력이 crack 2의 선단으로 집중되고 응력확대계수가 급격하게 증가하면서 균열진전이 시작되어 떨림현상도 함께 발생하는 것을 볼 수 있다. 결국 감쇠 알고리 즘을 적용한 경우 떨림현상을 효과적으로 감소시킬 수 있으며, 해석의 안정성과 정확성을 향상시킬 수 있음을 확인할 수 있다.

    초기 균열사이의 간격 L이 상대적으로 작은 0.0714m일 경우 Fig. 8과 같이 crack 1만 진전한다. crack 1은 휨 응력의 영향으로 약간 휘어지는 경향이 있지만 외력이 가해지는 지점으로 계산된 초기방향을 대체로 유지하면서 성장한다. 균열진전 방향을 정확하게 예측하기 위해 동적 응력확대 계수를 정확하게 계산하는 것이 중요한데, Fig. 9를 보면 감쇠효과를 반영한 경우 떨림현상이 효과적으로 감소된 것을 확인할 수 있다. 특히, Fig. 9(a)의 일부분을 확대한 Fig. 9(b)는 떨림현상에 대한 감쇠영향을 확실하게 보여주고 있다. 동적 균열전파해석 시 응력확대계수의 떨림현상이 심하게 발생하면 오차가 누적되고 해의 정확성이 떨어지는데, 감쇠 알고리즘은 이 상황을 효과적으로 개선하고 있다.

    마찬가지로 Fig. 10(a), (b)에서도 감쇠효과의 도입이 동적 응력확대계수 값의 떨림현상을 크게 억제하고 있음을 확인할 수 있다. 다중균열 문제에서도 단일균열 문제와 비슷하게 식 (15)를 사용한 경우와 식 (17)을 사용한 경우가 큰 차이를 보이지 않았다. 물론 MLS 차분법을 중앙차분법을 이용해 정식화하였으므로 시간간격 조정을 통해 해석의 정확도를 향상시키는 것이 가능하다.

    4.결 론

    본 논문에서는 Rayleigh 감쇠 알고리즘을 적용한 MLS 차분법을 활용하여 동적 균열전파 해석을 수행하였다. MLS 차분법은 Taylor 다항식과 이동최소제곱법을 이용하여 미분 근사식을 도출하기 때문에 절점만으로 복잡한 형상을 갖는 고체역학 동적문제를 효과적으로 해석할 수 있으며, 절점의 이동 및 추가가 자유로워 진전되는 균열을 절점만으로 간편 하게 모사할 수 있다. 기존의 강정식화 기반의 무요소 해석법은 비례감쇠 알고리즘을 제대로 반영하지 못했고 또 균열과 같이 특이성을 갖는 시간해석 문제의 경우 해의 떨림현상을 효과적 으로 억제하지 못해 해석의 정확성과 안정성을 확보하지 못했다. 그러나 본 연구에서는 MLS 차분법에 Rayleigh 감쇠 알고리즘을 적용하고 강정식화에서도 질량과 강성에 대한 비례감쇠 효과를 적절하게 반영하였으며, 중앙차분법을 속도항의 lagging과 함께 적용하여 계산효율성을 높였다.

    감쇠 알고리즘을 포함한 MLS 차분법의 동적 균열문제에 대한 성능을 확인하기 위해 단일균열 및 다중균열의 전파문제를 해석하였다. 강정식화 기반의 무요소 해석방법들은 요소 단위의 수치적분을 사용하지 않기 때문에 근사함수를 구성하는 영향영역 크기에 영향을 받고 동적해석의 안정성이 떨어지는 경우도 있다. 특히, 형상이 복잡하거나 균열과 같은 특이성을 갖는 문제의 경우 이러한 현상이 두드러질 수 있는데 본 연구에서 개발한 수치기법은 동적 균열전파 문제에 있어서 충분한 해석 안정성이 확보되는 것을 확인하였다. 단일균열 문제에 대한 해석결과를 이론해와의 비교를 통해 개발된 해석기법의 정확성을 확인하였다. 또한 다중균열 문제의 경우에는 정해와 비교할 수는 없었지만, 복잡한 기하형상이나 균열의 특이성으로 인해 발생하는 응력의 떨림현상을 효과적 으로 감소시킬 수 있음을 보였고, 동적 응력확대계수나 균열 진전 방향 예측에 대한 정확도를 향상시킬 수 있음을 확인 하였다. 또한, 유한요소법과 같은 약정식화에서 사용하는 감쇠 알고리즘을 강정식화에서도 약간의 수정을 통해 동일하게 적용할 수 있음을 보였으며, 향후 대변형을 동반한 동적문제의 해석에 적용하는 것과 같이 강정식화 기반의 수치기법에 대한 적용성을 높이는 계기가 될 수 있을 것으로 기대된다.

    감사의 글

    이 논문은 2014년도 정부(미래창조과학부)의 재원으로 한 국연구재단의 지원을 받아 수행된 기초연구사업임(No. 2014 R1A1A1002000). 또한, 본 연구는 국토교통부 철도기술연 구사업의 연구비지원(16RTRP-B104237-02)에 의해 수행 되었습니다.

    Figure

    COSEIK-29-6-583_F1.gif

    Crack growth modeling using nodal operation

    COSEIK-29-6-583_F2.gif

    Mode I single crack propagation problem

    COSEIK-29-6-583_F3.gif

    Normalized energy release rates with damping 3% and without damping

    COSEIK-29-6-583_F4.gif

    Multiple crack propagation problem

    COSEIK-29-6-583_F5.gif

    Computed crack growing path(L=0.0778m)

    COSEIK-29-6-583_F6.gif

    Relative dynamic stress intensity factors of crack 1(L=0.0778m)

    COSEIK-29-6-583_F7.gif

    Relative dynamic stress intensity factors of crack 2(L=0.0778m)

    COSEIK-29-6-583_F8.gif

    Computed crack growing path(L=0.0714m)

    COSEIK-29-6-583_F9.gif

    Relative dynamic stress intensity factors of crack 1(L=0.0714m) (a) a entire view (b) a close up shot for rectangular area

    COSEIK-29-6-583_F10.gif

    Relative dynamic stress intensity factors of crack 2(L=0.0714m) (a) a entire view (b) a close up shot for rectangular area

    Table

    Reference

    1. Atluri SN , Zhu T (1998) A New Meshless Local Petrov-Galerkin (MLPG) Approach in Computational Mechanics , Comput. Mech, Vol.22 ; pp.117-127
    2. Belytschko T , Lu YY , Gu L (1994) Element-Free Galerkin Methods , Int. J. Numer. Methods Eng, Vol.37 ; pp.229-256
    3. Freund LB (1990) Dynamic Fracture Mechanics, Cambridge University Press,
    4. Gu YT , Liu GR (2001) A Meshless Local Petrov- Galerkin (MLPG) Method for Free and Forced Vibration Analyses for Solids , Comput. Mech, Vol.27 ; pp.188-198
    5. Gu YT , Liu GR (2005) A Meshfree Weak-Strong (MWS) Form Method for Time Dependent Problems , Comput Mech, Vol.35 ; pp.134-145
    6. Hall JF (2006) Problems Encountered from the Use (or Misuse) of Rayleigh Damping , Earthq. Eng. & Struct. Dyn, Vol.35 ; pp.525-545
    7. Hughes T (1987) The Finite Element Method: Linear Static and Dynamic Finite Element Analysis, Prentice-Hall Inc,
    8. John R , Shah SP (1990) Mixed Mode Fracture of Concrete Subjected to Impact Loading , ASCE J. Struct. Eng, Vol.166 ; pp.585-602
    9. Kim DW , Kim YS (2003) Point Collocation Methods Using the Fast Moving Least-Square Reproducing Kernel Approximation , Int. J. Numer. Methods Eng, Vol.56 ; pp.1445-1464
    10. Lee SH , Kim KH , Yoon YC (2016) Particle Difference Method for Dynamic Crack Propagation , Int. J. Impact Eng, Vol.87 ; pp.132-145
    11. Lee SH , Yoon YC (2004) Meshfree Point Collocation Method for Elasticity and Crack Problems , Int. J. Numer. Methods Eng, Vol.61 ; pp.22-48
    12. Liu WK , Jun S , Zhang Y (1995) Reproducing Kernel Particle Methods , Comput Methods Appl. Mech. & Eng, Vol.191 ; pp.1421-1438
    13. Long SY , Liu KY , Hu DA (2006) A New Meshless Method Based on MLPG for Elastic Dynamic Problems , Eng. Anal. Bound. Elem, Vol.30 ; pp.43-48
    14. Sladek J , Sladek V , Mang HA (2003) Meshless LBIE Formulations for Simply Supported and Clamped Plates under Dynamic Load , Comput. & Struct, Vol.81 ; pp.1643-1651
    15. Omidi O , Valliappan S , Lotfi V (2013) Seismic Cracking of Concrete Gravity Dams by Plastic- Damage Model Using Different Damping Mechanisms , Finite Elem. Anal. & Des, Vol.63 ; pp.80-97
    16. Yoon YC , Kim KH , Lee SH (2014) Analysis of Dynamic Crack Propagation Using MLS Difference Method , J. Comput. Struct. Eng. Inst. Korea, Vol.27 ; pp.17-26
    17. Yoon YC , Lee SH , Belytschko T (2006) Enriched Meshfree Collocation Method with Diffuse Derivatives for Elastic Fracture , Comput. & Math. Appl, Vol.51 (8) ; pp.1349-1366
    18. Yoon YC , Song J-H (2014) Extended Particle Difference Method for Weak and Strong Discontinuity Problems: Part I , Derivation of the Extended Particle Derivative Approximation for the Representation of Weak and Strong Discontinuities. Comput. Mech, Vol.53 (6) ; pp.1087-1103
    19. Wen PH , Aliabadi MH (2007) An Improved Meshless Collocation Method for Elastostatic and Elastodynamic Problems , Commun. Num. Methods Eng, Vol.24 ; pp.635-651