Journal of the Computational Structural Engineering Institute of Korea. 2014. 17-26
https://doi.org/10.7734/COSEIK.2014.27.1.17

ABSTRACT


MAIN

1. 서 론

해석영역 내부에 존재하는 경계면의 처리시 요소(mesh)의 제약이 발생하는 유한요소법을 보완할 수 있는 대안으로 제시된 무요소법은 1990년대 중반 Element-Free Galerkin (EFG; Belytschko et al., 1994)법과 재생커널 점별법(Reproducing Kernel Particle Method; Liu et al., 1995) 등의 형태로 제시되었다. 개발 초기에 무요소법은 수치해석시 절점만을 이용하여 해석할 수 있다는 장점이 부각되었으나, 실제로는 무요소법이 Galerkin법을 기반으로 지배방정식에 대한 약형식(weak form)을 유지하고 있어 수치적분을 피할 수 없을 뿐만 아니라, 필수경계조건 처리가 어려운 단점이 제기되었다. 이러한 단점을 극복하기 위하여 필수경계조건 처리에 대한 연구들이 진행되었으나, 대부분의 연구들이 유한요소와 결합하거나 Lagrange Multiplier Method 혹은 Penalty Method와 같은 구속조건식을 도입하여 필수경계를 처리했기 때문에 전체 알고리즘이 유한요소법과 큰 차이 없이 정식화되는 경향을 보였다(Krongauz et al., 1996; Gosz et al., 1996; Fernandez-Mendez et al., 2004). 결국, 이러한 접근은 요소망의 제약으로부터 자유롭고자 했던 무요소법의 장점을 희석시키는 결과를 초래했다.

한편, 무요소법의 장점을 극대화 시키고자 Galerkin법에 기초한 약정식화를 탈피하여 Collocation 방법을 시도한 접근법도 있었다(Zhu et al., 1998). Zhang 등(2001)에 의해 최소제곱 콜로케이션법(Least-Squares Collocation Method)의 형태로 무요소법에 Collocation법이 합성되었으나, 해석의 정확도가 기대만큼 좋지 않았고 최소제곱법에 의한 계 방정식 구성이 기존 약형식 기반의 무요소법보다 더 많은 계산량을 필요로 하여 처음 기대만큼 좋은 결과를 얻지 못했다. 그 후, 요소망의 의존성을 완전히 탈피하기 위하여 강형식(strong form)기반의 무요소 콜로케이션 방법(Meshfree Point Collocation Method)이 제안되었다(Kim et al., 2003; Lee et al., 2004). 이 방법은 고속 미분근사를 바탕으로 지배 미분방정식을 그대로 이산화하기 때문에 수치적분과 같이 약형식으로부터 야기되는 수치적인 계산부담을 효과적으로 피할 수 있으며, 필수경계조건 처리도 매우 수월했다. 이러한 장점을 기반으로 정적균열문제 해석(Yoon et al., 2007) 및 계면경계문제 해석(Kim et al., 2007, Yoon et al., 2009) 등에 관한 연구가 수행되었고, 고체의 동적해석으로 확장되었으나(Yoon et al., 2012), 충격하중(impact loading)을 받는 동적균열의 전파해석까지 연결되지는 못했다.

동적균열전파 문제는 시간이 경과할 때 해석모델 내의 균열이 진전하면서 기하학적 형상이 계속 변화하기 때문에 약형식을 기반으로 하는 경우 요소망 구성의 제약이 발생한다. 이런 이유로 무요소법이 등장할 때부터 동적균열전파 해석으로의 적용이 많은 관심을 받아 왔다. 특히, EFG법을 기반으로 동적균열전파 해석에 대한 연구들이 많이 이루어졌는데(Belytschko et al., 1996; Krysl et al., 1999; Belytschko et al., 2000), 이런 연구들도 필수경계조건 처리와 같이 EFG법이 갖고 있는 태생적인 단점을 극복하기 위한 구속조건식들을 포함하고 있어서 무요소법의 장점을 생각만큼 극대화시키지 못했다. 이러한 단점을 극복하고자 유한요소법과 무요소법을 결합하는 기법들이 제안되었고(Belytschko et al., 1995), 이를 동적균열전파 문제에 적용한 연구들도 있었다(Gu et al., 2008; Rajesh et al., 2010). 그러나 이러한 방법들 역시 무요소법이 지향하는 요소망 제약을 완벽하게 극복하지는 못했다.

본 논문에서 제시하는 MLS(Moving Least Squares) 차분법은 Taylor 전개에 기초한 고속 미분근사와 미분방정식의 강정식화(strong formulation)를 기반으로 한 수치기법이다. 절점모델을 기반으로 이동최소제곱법을 통해 도출된 미분근사식이 고차 미분까지 간편하게 계산해주기 때문에 지배 미분방정식을 강형식 그대로 이산화하는 것이 가능하다. 결과적으로 상당한 계산량을 요구하는 수치적분 과정이 불필요하게 되고, 전 해석과정에서 요소망 정보를 전혀 필요로 하지 않는다. 이러한 특징들은 무요소법이 애초에 추구했던 장점을 극대화 시켜주는 동시에 계산 효율성도 향상시켜 주는 효과가 있다. 본 논문에서는 MLS 차분법을 이용하여 고체역학문제에서의 동적균열전파 해석 알고리즘을 제안하고, 수치예제를 통하여 제안된 알고리즘의 정확성을 검증하고자 한다.

2. MLS 차분법의 동적해석 알고리즘

이 장에서는 본 논문에서 사용하고 있는 MLS 차분법의 특징을 설명하기 위하여 고속 미분근사식의 유도과정을 설명한다. 그리고 동적해석 수행을 위한 지배방정식의 이산화 알고리즘을 바탕으로 전체 계 방정식을 구성하는 과정에 대한 설명이 뒤따른다.

2.1 고속 미분근사식의 유도

MLS 차분법은 Taylor 전개를 바탕으로 이동최소제곱법을 이용하여 형상함수를 도출하기 때문에 미분근사식의 유도과정이 간편하다. 미분근사식의 유도시 다중지수 표기법(multi-index notation)을 사용한다. https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F101.pngn차원 좌표를 갖는 벡터를 의미하며, https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F103.png의 α-차 지수와 편미분은 다음과 같이 정의된다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F105.png      (1a)

https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F106.png      (1b)

여기서, https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F107.png을 의미한다. 예를 들어, 2차원의 경우 https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F103.png에 대한 α=(1,1)차 지수는 식 (1a)에 의하여 https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F110.png와 같이 표현되며, https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F103.png1https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F206.png방향, https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F103.png2는 방향 y성분을 나타낸다. 이 때, 임의의 기준점 https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F103.png0에 대한 https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F103.png점에서의 https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F116.png번 미분 가능한 함수 https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F117.png를 고려하고 m차 이상의 고차항을 무시하면, https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F117.png에 대한 m차 Taylor 다항식을 다음과 같이 정의할 수 있다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F121.png      (2)

식 (2)를 다항식 기저벡터와 미분계수에 대한 계수벡터의 스칼라곱으로 분해하면 다음과 같다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F122.pnghttps://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F123.pnghttps://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F124.png

https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F125.png      (3)

여기서, https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F126.png, https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F127.png, https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F128.png이고, https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F129.png는 m차 다항식 기저벡터이며, https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F131.pnghttps://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F132.png의 각 성분들을 순서대로 정렬한 미지벡터이다. 이동최소제곱법을 이용하여 식 (3)을 식 (4)와 같이 x0의 영향영역 내에 포함된 N개의 절점에 대해서 가중 함수식 형태의 잔차식(residual)을 정리하고 그것을 최소화하는 과정을 통해 미분차수별 도함수를 포함하는 https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F131.png를 얻는다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F136.png      (4)

여기서, https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F137.png는 compact support를 갖는 가중함수이고, https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F138.pnghttps://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F139.png에서의 해(solution) 값이다. 가중함수에서 r은 기준점 x0에 대한 영향반경을 나타낸다. 식 (4)를 최소화하는 과정에서 기준점 x0를 실제로 함수값을 계산하려는 위치 x로 치환하여 정리하면 식 (5)와 같이 c(x)가 유도된다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F145.png

https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F146.png      (5)

여기서, https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F147.png는 x의 영향반경 안에 속하는 절점변위(nodal solution) 값들을 모아 놓은 벡터이다. 식 (3)에서 x0를 x로 치환한 후, 식 (5)를 대입하면 최종적인 Taylor 다항식이 얻어진다. 결과적으로, 가중함수, 기저함수벡터, 절점변위를 이용하여 미분근사식을 유도할 수 있음을 알 수 있으며 일반화된 형상함수(shape function)를 사용하여 다음과 같이 정리할 수 있다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F151.png      (6)

여기서, 일반화된 형상함수 https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F152.png는 절점 I에 대한 형상함수의 α차 미분에 대한 근사값을 나타낸다. MLS 차분법에서는 형상함수를 구하는 과정에서 미분근사식이 함께 구해지기 때문에 추가의 계산 없이도 원함수(original function)에 대한 고속 미분근사가 가능하다. 또한, 미분근사식의 유도과정이 Taylor 전개에 기초를 두고 있기 때문에 수학적 기초가 탄탄하고, 기저벡터의 차수 조정을 통해 미분근사식의 정확도를 손쉽게 조절할 수 있다. 고체역학문제의 지배방정식은 2차 미분이 필요하기 때문에 2차(quadratic) 기저벡터를 기본으로 사용하며, 필요하다면 3차(cubic), 4차(quartic) 등의 고차 근사식으로 확장하는 것도 매우 용이하다.

2.2 지배방정식의 이산화

미소변형을 가정한 고체역학 문제는 다음과 같이 내부영역과 경계에 대한 3개의 지배방정식으로 정의할 수 있다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F155.png in       (7a)

https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F157.png on https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F158.png      (7b)

https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F159.png on https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F160.png      (7c)

여기서, 는 내부영역, https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F158.png는 자연경계, https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F160.png는 필수경계를 의미하며, σ는 응력텐서, b는 체적력, ρ는 밀도, a는 가속도이다. n과 https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F169.png는 각각 자연경계에서의 단위수직벡터와 미리 규정된 표면력을, u와 https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F171.png는 각각 미지변위와 필수경계에서 미리 규정된 변위값을 나타낸다. 각 재배방정식을 결합하여 하나의 시스템 방정식을 만들기 위해서는 지배방정식의 모든 변수를 변위로 변환하는 과정이 필요하다. 식 (7a)의 운동방정식은 체적력을 무시하고 선형탄성재료에 대한 구성방정식을 적용하여 변위에 대해 표현하면 아래와 같은 Navier 방정식으로 유도된다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F172.png in       (8)

여기서, https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F174.png는 Laplace 연산자이고, λ와 µ는 Lame 상수이다. 같은 방법으로, 식 (7b)의 자연경계 조건식을 변위에 대해 표현하면 다음과 같다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F177.png on https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F158.png      (9)

식 (7c)의 필수경계방정식은 변위로 정의되어 있어 따로 변환할 필요가 없다. 수치적인 편의를 위해 식 (8)과 식 (9)를 미분연산자 과 를 사용하여 다음과 같이 정의한다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F181.png      (10)

https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F182.png      (11)

응력항이 변위로 변환된 운동방정식인 식 (8)은 가속도 항을 포함하고 있기 때문에 시간적분을 이용하여 변위로 변환하는 작업이 필요하다. 본 논문에서는 아래와 같은 Newmark 시간적분법(Newmark, 1959)을 적용한다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F183.png      (12a)

https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F184.png      (12b)

여기서, https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F185.png는 속도, https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F186.png는 시간단계(time step)의 크기를 의미하고, 상첨자는 시간단계를 가리킨다. 식 (12a)와 식 (12b)를 변위와 가속도에 관한 식으로 정리하고, 미분연산자 L을 사용하여 식 (8)을 다시 쓰면 다음과 같이 변위에 관한 운동방정식을 얻을 수 있다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F188.png

https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F189.png      (13)

이 때, 2차원 문제의 경우 절점 I에 대해 미분연산자 L과 단위행렬 I는 크기가 2×2인 행렬로 표현되며, 각각의 구성성분은 아래와 같다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F193.png      (14a)

https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F194.png      (14b)

https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F195.png      (14c)

https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F196.png      (15a)

https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F197.png      (15b)

여기서, 아래첨자 Ihttps://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F199.png의 영향영역 내에 포함된 절점들의 번호이다.

자연경계조건에 대한 미분연산자인 B행렬의 구성성분은 다음과 같다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F201.png      (16b)

https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F202.png      (16b)

https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F203.png      (16b)

여기서, https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F204.pnghttps://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F205.png는 각각 단위수직벡터의 https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F206.png, https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F207.png 방향성분이다. 이제 https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F208.png1 시간단계에서 내부, 자연경계, 필수경계에 놓인 절점들에 대한 지배방정식을 다음과 같은 일반화된 행렬식의 형태로 나타낼 수 있다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F209.png      (17)

여기서, https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F199.png가 포함된 영역에 따라서 https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F211.pnghttps://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F212.png의 구성이 달라진다. 이렇게 각 절점 위에서 지배방정식을 모두 조립(assemble)하면 전체 계 방정식을 다음과 같이 얻는다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F213.png      (18)

여기서, U는 전제 절점변위 해를 모두 조립한 벡터를 의미한다. 일반화된 계수행렬 K는 (NXN), F는 (NX1)의 크기를 갖는데, 이 때 N은 전체 절점수와 해석의 차원 n을 곱한 값이다. K의 구성성분은 https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F199.png가 속해 있는 영역에 따라 다음과 같이 주어진다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F224.png, https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F199.png in       (19a)

https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F225.png, https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F199.png in https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F158.png      (19b)

https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F226.png, https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F199.png in https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F160.png      (19c)

한편, 일반화된 힘벡터 는 가 포함된 영역에 따라 아래와 같이 정의된다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F235.png

https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F236.png

      ,https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F199.png in       (20a)

https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F239.pnghttps://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F199.png in https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F158.png      (20b)

https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F240.pnghttps://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F199.png in https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F160.png      (20c)

결국, 식 (18)의 전체 계 방정식을 풀어서 n+1 시간단계에서의 모든 절점변위를 구할 수 있다. 고체역학문제의 동적해석을 위한 2차원 시스템 방정식에 대한 보다 자세한 내용은 Yoon 등(2012)을 참고할 수 있다.

3. 균열전파 메커니즘

본 절에서는 앞서 유도된 계 방정식을 바탕으로 균열전파 해석을 수행하기 위해 절점모델 상에서 균열을 묘사하는 방법을 설명하고, 균열진전 여부와 그 방향을 결정하는 방법에 대해 서술한다.

3.1 가시한계법(Visibility criterion)

MLS 차분법의 미분근사식(또는 형상함수)은 이동최소제곱법에 근거한 Taylor 전개를 기반으로 기준점과 주변 점들 간의 거리관계에 의해서 정의되기 때문에 균열의 불연속성 묘사시 균열 주변점들 사이의 거리관계를 정립하는 것이 중요하다. 본 논문에서는 균열 반대쪽 절점의 영향을 제거하는 가시한계법(Visibility criterion; Fleming et al., 1997)을 사용하여 불연속면을 묘사하였다. 균열이 없는 경우, 영향영역은 반경 r을 반지름으로 하는 원의 형태이나, 균열선단 근처에 위치한 기준점의 영향반경은 Fig. 1과 같이 기준점과 선단을 통과하는 직선을 이용하여 균열 반대편의 절점들을 구분하고 그 기여도(contribution)를 무시하였다. 즉, Fig. 1에서 검은색으로 채워진 절점만 계산에 포함되고, A점에 대해 균열 반대편에 있는 속이 빈 원으로 표시된 절점들은 근사함수 계산시 제외된다. 또한, 기준점이 영향영역 내에 이웃절점들을 일정한 개수만큼 포함시키도록 영향반경 r을 설정하는 알고리즘(Kim et al., 2004)을 적용하여, 균열전파에 따라 해석모델의 절점배치에 부분적인 변화가 생겨도 일정한 조건하에서 형상함수가 구성될 수 있도록 하였다.

가시한계법에 대한 실제적인 구현은 형상함수 구성과정 중 가중함수를 정의할 때 이루어진다. 식 (21)에는 본 연구에서 사용한 가중함수를 제시했다. 이동최소제곱법 적용시 기준점으로부터의 거리가 영향반경보다 먼 이웃절점은 자연스럽게 가중함수값이 0이 되어 형상함수 구성에서 제외되는 것처럼 균열 반대편에 위치한 이웃절점들은 영향영역 내에 포함되었을지라도 가시한계법에 의해 불연속면 반대편에 위치한 것으로 인식하여 그 영향을 무시한다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F249.png (21)

https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F248.png
Figure

Influence domain for shape function construction based on the Visibility criterion

3.2 동적 응력확대계수(Stress Intensity Factor)

동적균열전파 해석의 가장 중요한 부분 중 하나는 시간에 따른 균열의 성장여부와 그 방향을 결정하는 것이다. 이를 위하여 본 논문에서는 동적 응력확대계수를 사용하였다. 동적 응력확대계수는 에너지해방률(energy release rate)를 이용하여 구할 수 있는데, 우선, 모드 I과 모드 II 응력확대계수가 포함된 에너지해방률 G는 평면변형률(plane strain) 상태에서 다음과 같이 정의할 수 있다(Freund, 1990).

https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F251.png      (22)

여기서, E는 탄성계수, https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F253.png는 포아송비이고 KK는 모드Ⅰ, 모드Ⅱ 응력확대계수이며, https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F256.png는 균열전파 속도이다. 나머지 변수들은 아래와 같이 주어진다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F257.png, https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F258.png      (23a)

https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F263.png,https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F262.png      (23b)

https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F261.png      (23c)

여기서, Cd는 팽창파(dilatational wave) 속도이고, Cs는 전단파(shear wave) 속도이다.

본 논문에서는 J-적분을 면적분 형태로 변환시켜 에너지해방률을 계산하였으며(Moran et al., 1987), 이를 바탕으로 식 (22)로부터 모드 I와 모드 II의 응력확대계수를 다음과 같이 결정하였다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F264.png      (24)

위 식에서 https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F265.pnghttps://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F266.png는 탄성재료를 가정할 때 각각 모드 I과 모드 II의 가상 변위장을 결합한 상호작용(interaction) 에너지해방률을 의미한다(Yau et al., 1980). Fig. 2는 균열에 대한 J-적분을 수행하기 위한 영역을 나타내는데, 적분영역에 새로운 균열방향으로 설정되는 가상균열이 추가되어, 기존의 균열 경계면 https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F267.pnghttps://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F268.png외에도 가상균열 경계면 https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F269.pnghttps://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F270.png가 계산에 포함된 것을 볼 수 있다. 이러한 적분과정을 통해 https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F265.png는 다음과 같이 유도된다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F280.png
Figure 2

The domain of interaction J-integral

https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F272.png

https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F273.png

https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F274.png

https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F275.png      (25)

여기서, 상첨자 https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F276.png는 시간미분을 가리킨다. 임의의 가중함수 https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F277.png는 식 (26)과 같이 정의되고, m은 적분영역의 단위수직벡터이다. 모드 II에 대한 https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F266.png도 동일하게 구할 수 있다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F281.png      (26)

이와 같이 계산된 동적 응력확대계수를 사용하여 모드 I 등가 응력확대계수(equivalent Stress Intensity Factor)를 다음과 같이 정의할 수 있다(Freund, 1990).

https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F282.png      (27)

여기서, θ는 균열진전방향이다. 계산된 등가 응력확대계수를 해석대상 재료의 동적파괴인성치(dynamic facture toughness) KIC와 비교하여 다음과 같은 기준을 만족할 때 균열이 성장하도록 하였다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F285.png      (28)

즉, 등가응력확대계수가 동적파괴인성치를 넘게 되면 자연스럽게 균열이 성장하기 시작하고, 이 때, θ값이 균열의 진전방향이 된다.

4. 수치예제

4.1 모드 I 균열전파문제

MLS 차분법의 모드 I 동적균열전파 해석에 대한 정확성을 검증하기 위하여 Fig. 3과 같이 인장응력을 받는 편측 균열을 포함하고 있는 무한체를 해석한다(Freund, 1990). 초기 균열길이 a0는 5m이고, 무한체 상단에 작용하는 응력 σ0는 시간에 따라 일정하게 1000N/m2을 적용하였으며, 평면변형률 상태를 가정하였다. 밀도 ρ=7800kg/m3, 탄성계수 E=211×109N/m2, 포아송 비 https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F291.png=0.30으로 가정하였으며, 시간단계 크기 https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F292.png는 0.5×10-6로 200단계의 해석을 수행하였다. 해석대상에 기본적으로 1000개(50×20)의 절점을 배치하고, 균열선단 주변에 149개의 절점을 추가하여, 총 1149개의 절점을 배치하였다. 본 해석기법은 절점의 추가배치가 자유롭다는 것을 다시 한번 주목할 필요가 있다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F294.png
Figure 3

Mode 1 crack propagation problem (the crack propagates in a prescribed direction)

Fig. 4에서는 본 예제에 대한 이론해(analytic solution)와 MLS 차분법으로 계산한 에너지해방률을 비교했다. 위쪽에 그려진 선은 균열이 전파하지 않는 경우의 시간에 따른 에너지해방률이고, 아래쪽 선은 외부하중에 의한 응력파가 균열면에 도달하는 즉시 균열이 진전되기 시작하는 경우를 나타낸다. 전자는 균열속도 https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F295.png를 0으로, 후자는 https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F295.png를 0.4cs로 설정했다. 검증을 위해 MLS 차분법과 EFG법으로 계산된 에너지해방률 값을 함께 도시했다. EFG법의 경우 총 절점수를 비슷하게 맞추기 위하여 1150개(50×23)의 절점을 사용하였다. 균열이 진전하지 않는 경우(https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F295.png=0), 본 연구의 해석결과는 이론해와 매우 잘 일치하며, 균열이 진전하는 경우(https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F295.png=0.4cs)에는 계산된 에너지해방률이 약간의 떨림현상(oscillation)을 보인다. 추가적인 비교를 위해 총 4100개의 절점을 사용한 EFG법의 해석결과를 함께 도시했다. MLS 차분법은 균열선단 주변에 절점을 추가하여 EFG법보다 훨씬 적은 수의 절점으로도 상당히 높은 정확도의 해석을 할 수 있음을 알 수 있다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F301.png
Figure 4

Energy release rates for a stationary crack and a crack propagation at t=1.0H/cd

Fig. 5는 응력파가 균열면에 도달한 후에도 균열이 일정 시간동안 진전되지 않다가 특정한 시간이 되었을 때 갑자기 진전을 시작하는 경우를 모사한 결과이다(https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F295.png를 0으로 유지하다가, t=1.5H/cd일 때 0.4cs로 바꾸어 주었다.). 4100개 절점으로 모델링한 EFG법은 균열진전이 시작되는 순간에 에너지해방률이 급격이 감소하는 현상을 대체로 잘 모사하고 있으나, MLS 차분법은 약간의 overshooting 하고 그 후에 이론값을 대체로 잘 따라가는 모습을 보인다. 또한, MLS 차분법이 EFG법에 비교하여 절점수는 상대적으로 작지만 오히려 떨림현상(oscillation)이 적게 나타나는 것을 확인할 수 있다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F308.png
Figure 5

Energy release rates of MLSDM and EFGM for a crack propagation at t=1.5H/cd

4.2 혼합모드 균열전파문제

본 절에서는 Fig. 6과 같은 Kalthoff와 Winkler(1987)의 실험에 대한 모사를 통해 충격하중에 의해 균열이 순간적으로 혼합모드 상태로 성장하는 문제를 해석한다. 대칭조건을 이용하여 Fig. 7과 같이 위쪽 절반만 모델링했다. 해석모델 구성시 2704(52×52)개의 절점을 기본적으로 배치하고 균열주변에 131개를 추가하여, 총 2835개의 절점을 사용했다. 평면병형률로 가정하였고, 외부에서 가해지는 충격하중의 속도 https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F256.png는 16.5m/s로 일정하게 적용하였다. 탄성계수 E=190GPa, 포아송 비 https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F291.png=0.30, 동적파괴인성치 KIC=68 MPahttps://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F315.png, 밀도 ρ=8000kg/m3로 가정했고, 시간단계 크기 https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F292.png는 2.0×10-7s로 설정하고 300단계를 해석했다. 초기 균열길이는 0.05m, 균열 속도 https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F295.png는 0.15cs로 가정하였으며, 기타의 다른 조건들은 앞 예제와 동일하게 적용하였다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F321.png
Figure 6

Experiment of dynamic crack propagation done by Kalthoff and Winkler(1987)

https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F322.png
Figure 7

Dynamic crack growth simulation for an edge- cracked plate under mixed-mode loading

Fig. 8에는 MLS 차분법으로 계산한 균열의 전파경로를 도시했다. 균열이 단계적으로 전파될 때, 균열의 전파방향이 일관성 있게 계산되고 있으며, 초기 균열선단에서 마지막 단계에서 해석된 균열선단까지 연결한 선분의 각도가 약 72도로 계산됐다. Fig. 9에는 EFG법으로 계산한 균열전파 경로를 도시했다. MLS 차분법과 비교하기 위해 2916(54×54)개의 절점을 사용했다. EFG법으로 계산된 균열성장 각도는 약 76도로 나타났다. Table 1에 정리된 것과 같이 실험으로부터 측정된 균열 성장각은 70도인데, MLS 차분법이 EFG법보다 실험값에 더 근접한 결과를 주는 것을 알 수 있다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F323.png
Figure 8

Crack path computed by the MLSDM (The average propagation angle is 72°)

https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F324.png
Figure 9

Crack path computed by the EFGM (The average propagation angle is 76°)

Fig. 10은 MLS 차분법과 EFG법으로 계산된 동적 응력확대계수값의 시간에 따른 변화를 보여준다. 이 때 동적 응력확대계수는 재료의 동적 파괴인성치로 나누어 정규화하였다. 두 방법 모두 초기에 균열이 발생하지 않는 상태에서는 모드 II의 응력확대계수가 음의 값을 가지며 어느 정도 발달한 후, 균열이 진전하기 시작한 이후로는 모드 I 응력확대계수가 재료의 동적 파괴인성치의 2배 가까이 증가하고, 모드 II 응력확대계수는 0 주변으로 떨리면서 그 값을 유지하는 것으로 나타났다. 결국, 처음 균열이 정지상태에서 진전하기 전까지는 모드 II의 영향이 상대적으로 컸지만, 균열이 진전을 시작한 이후로는 모드 I의 영향이 지배적으로 되기 때문에 시간이 경과하면서 균열이 일정한 방향을 유지하면서 진전하게 되는 것을 미루어 알 수 있다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2014-027-01/03TK022014270104/images/10.7734.27.1.17.F325.png
Figure 10

Dynamic stress intensity factors during crack growth computed by the MLSDM and the EFGM

Table 1

Comparison of the average propagation angles

MethodsExperimentMLSDMEFGM
degree(°)707276

5. 결 론

본 논문에서는 수치적분이 없고 구속조건식의 도입에 의한 필수경계조건 처리가 필요없어 무요소법의 장점을 극대화 할 수 있는 절점기반의 MLS 차분법을 이용하여 동적균열전파 해석을 수행하였다. MLS 차분법은 이동최소제곱법과 Taylor 다항식을 이용하여 절점만으로 구성된 수치모델을 바탕으로 해석을 수행하기 때문에 요소 및 요소망에 전혀 의존적이지 않다. MLS 차분법의 고속 미분근사식이 고차의 미분까지 손쉽게 계산해주기 때문에 지배방정식을 강형식으로 이산화하는데에 매우 유리하다. 운동방정식을 변위만의 식으로 정리한 후 Newmark 법을 사용하여 시간적분을 수행하였다. 또한, 가시한계법을 사용하여 균열의 불연속성을 간단한 절차를 통해 묘사하였고, 균열전파 여부를 판단하기 위해 동적 응력확대계수를 산정하여 활용하였다. 동적 응력확대계수 산정을 위해 J-적분을 면적적분으로 변환하고 상호작용 적분법을 적용하여 동적 에너지해방률을 계산하였다.

MLS 차분법을 바탕으로 개발된 동적균열전파 해석 알고리즘의 정확성을 검증하기 위해 모드 I상태와 혼합모드 상태의 균열전파 문제를 해석하였다. 이를 통해 절점만 사용하여 균열의 거동을 묘사하고, 운동방정식을 풀어 시간에 따른 균열의 전파과정을 안정적으로 모사할 수 있음을 보였다. MLS 차분법으로 얻어진 해석결과를 이론값 또는 실험값과 비교하였으며, EFG법으로 계산한 결과와도 비교하였다. 대체로 본 논문에서 제시한 알고리즘이 EFG법보다 떨림현상이 적게 나타났으며, 균열성장 각도에 대한 계산 정확성도 상대적으로 높다는 것을 확인할 수 있었다. 매 해석단계마다 절점을 재배치할 필요없이 균열선단 주변에서 절점 한두 개만 추가 또는 삭제하는 국부적인 수정만으로도 균열성장과정이 모사되기 때문에 계산 효율성도 우수하다. 더욱이 완벽하게 절점만 이용하는 해석이 이루어지기 때문에 해석모델 안에서 절점 추가배치가 자유로워 균열선단 주변에 절점을 손쉽게 추가할 수 있다. 향후에는 여러 개의 균열이 동시에 성장하는 문제 그리고 비선형 재료모델을 포함한 동적균열전파 해석으로 연구가 확장될 수 있을 것으로 기대된다.

Acknowledgements

이 논문은 2012년도 정부(교육과학기술부)의 재원으로 한국연구재단의 지원을 받아 수행된 기초연구사업임(2012R1 A1A2008248).

References

1
Belytschko, T., Lu, Y.Y., Gu, L.1994Element- free Galerkin MethodsInternational Journal for Numerical Methods in Engineering37229256610.1002/nme.1620370205
2
Belytschko, T., Organ, D., Gerlach, C.2000Element-free Galerkin Methods for Dynamic Fracture in ConcreteComputer Methods in Applied Mechanics and Engineering187385399610.1016/S0045-7825(00)80002-X
3
Belytschko, T., Organ, D., Gerlach, C.1995A Coupled Finite Element-Element-Free Galerkin MethodComputational Mechanics17186195610.1007/BF00364080
4
Belytschko, T., Tabbara M.1996Dynamic Fracture using Element-free Galerkin MethodsInternational Journal for Numerical Methods in Engineering39923938610.1002/(SICI)1097-0207(19960330)39:6%3C923::AID-NME887%3E3.0.CO;2-W
5
Fernandez-Mendez, S., Huerta, A. 2004Imposing Essential Boundary Conditions in Mesh-Free MethodsComputer Methods in Applied Mechanics and Engineering19312571275610.1016/j.cma.2003.12.019
6
Fleming, M., Moran, B., Belytschko, T.1997Enriched Element-Free Galerkin Methods for Crack Tip FieldsInternational Journal for Numerical Methods in Engineering4014831504610.1002/(SICI)1097-0207(19970430)40:8%3C1483::AID-NME123%3E3.0.CO;2-6
7
Freund, L.B.1990Dynamic Fracture MechanicsCambridge: Cambridge University Press
8
Gosz, J., Liu, W.K.1996Admissible Approximations for Essential Boundary Conditions in the Reproducing Kernel Particle MethodComputational Mechanics19120135610.1007/BF02824850
9
Gu, Y.T., Zhang, L.C.2008Coupling of the Meshfree and Finite Element Methods for Determination of the Crack Tip FieldsEnginnering Fracture Mechanics759861004610.1016/j.engfracmech.2007.05.003
10
Kalthoff, J.F., Winkler, S.1987Failure Mode Transition at High Rates of Shear LoadingInternational Conference on Impact Loading and Dynamic Behavior of Materials1185195
11
Kim, D.W., Kim, H.K.2004Point Collocation Method Based on The FMLSRK Approximation for Electromagnetic Field AnalysisIEEE Transactions on Magnetics40210291032610.1109/TMAG.2004.824612
12
Kim, D.W., Kim, Y.S.2003Point Collocation Methods using The Fast Moving Least-Square Reproducing Kernel ApproximationInternational Journal for Numerical Methods in Engineering5614451464610.1002/nme.618
13
Kim, D.W., Liu, W.K., Yoon, Y.C., Belytschko, T., Lee, S.H.2007Meshfree Point Collocation Method with Intrinsic Enrichment for Interface ProblemsComputational Mechanics4010371052610.1007/s00466-007-0162-1
14
Krongauz, Y., Belytschko, T.1996Enforcement of Essential Boundary Conditions in Meshless Approximations using Finite ElementsComputer Methods in Applied Mechanics and Engineering131133145610.1016/0045-7825(95)00954-X
15
Krysl, P., Belytschko, T.1999The Element Free Galerkin Method for Dynamic Propagation of Arbitrary 3-D CracksInternational Journal for Numerical Methods in Engineering44767800610.1002/(SICI)1097-0207(19990228)44:6%3C767::AID-NME524%3E3.0.CO;2-G
10.1002/(SICI)1097-0207(19990228)44:6<767::AID-NME524>3.0.CO;2-G
16
Lee, S.H., Yoon, Y.C.2004Meshfree Point Collocation Method for Elasticity and Crack ProblemsInternational Journal for Numerical Methods in Engineering612248610.1002/nme.1053
17
Liu, W.K., Jun, S., Zhang, Y.1995Reproducing Kernel Particle MethodsComputer Methods in Applied Mechanics and Engineering19114211438
18
Moran, B., Shih, C.F.1987Crack Tip and Associated Domain Integrals From Momentum and Energy BalanceEngineering Fracture Mechanics27615642
19
Newmark, N.M.1959A Method of Computation for Structural DynamicsJournal of the Engineering Mechanics Division, ASCE856794
20
Rajesh, K.N., Rao, B.N.2010Coupled Meshfree and Fractal Finite Element Method for Mixed Mode Two-dimensional Crack ProblemsInternational Journal for Numerical Methods in Engineering84572609
10.1002/nme.2910
21
Yau, J.F., Wang, S.S., Corten, H.T.1980A Mixed-mode Crack Analysis of Isotropic Solids using Conservation Laws of ElasticityJournal of Applied Mechanics47335341610.1115/1.3153665
22
Yoon, Y.C., Kim, D.J., Lee, S.H.2007A Gridless Finite Difference Method for Elastic Crack AnalysisJ. Comput Struct. Eng.203321327
23
Yoon, Y.C., Lee, S.H.2009Intrinsically Extended Moving Least Squares Finite Difference Method for Potential Problems with Interfacial BoundaryJ. Comput Struct. Eng.225411420
24
Yoon, Y.C., Kim, K.H., Lee, S.H.2012Dynamic Algorithm for Solid Problems using MLS Difference MethodJ. Comput. Struct. Eng.2521391486
25
Zhang, X., Liu, X.H., Song, K.Z., Lu, M.W.2001Least-squares Collocation Meshless MethodInternational Journal for Numerical Methods in Engineering5110891100610.1002/nme.200
26
Zhu, T., Atluri, S.N.1998A Modified Collocation Method and a Penalty Formulation for Enforcing the Essential Boundary Conditions in the Element Free Galerkin MethodComputational Mechanics21211222610.1007/s004660050296
페이지 상단으로 이동하기