1. 서 론
재료 파괴 및 손상 예측은 공학 및 과학 분야에서 매우 도전적인 문제이며 전산역학 해석 분야에서의 중요한 연구 분야이다. 페리다이나믹 이론(Silling, 2000)은 원거리 상관관계(Long-range interaction)에 기반을 둔 연속체 기반의 전산해석 방법론으로써 기존 해석법들이 재료파괴를 표현하는 한계를 극복하기 위해 제안되었다. 기존 방법론들과 달리 페리다이나믹에서는 변위장의 공간미분성을 가정하지 않기 때문에 균열 및 파괴 해석 모델과 같은 불연속 변위장을 가지는 문제를 해석하기에 적합하다.
페리다이나믹 모델은 비국부 모델로서 인접한 영역에 속한 절점들의 상호작용을 해석한 적분방정식을 사용하기 때문에 많은 계산메모리와 계산시간이 필요하다. 본 연구에서는 탄성체 모델과 연성으로 효율적이고 엄밀한 페리다이나믹 해석법을 소개한다. 탄성체 모델과 페리다이나믹 모델의 연성기법은 여러 연구자에 의해 제시되어 왔다(Han et al., 2012; Liu et al., 2012; Seleson et al., 2013). 그 중 Seleson 등(Seleson et al., 2013)은 1차원 페리다이나믹 모델과 파동방정식을 연성하는 방법론을 제안하였다. 이 방법론은 원자-연속체 연성(Fish et al., 2007)을 위해 사용되는 힘-기반 블랜딩(force-based blending)과 유사하며, 단일 모델에서 연성을 구현할 수 있다는 강점이 있다. 본 논문에서는 힘-기반 블랜딩 방법론을 통해 2차원 페리다이나믹 모델과 탄성체 모델을 연성하는 방법론(Seleson et al., 2014)을 간략히 소개한다. 또한 블랜딩 영역의 크기가 해석결과에 미치는 영향을 조사하고 정적파괴 모델에 적용하여 기존 방법론들과 비교한다.
2. 페리다이나믹 이론
페리다이나믹 모델은 운동방정식에서 응력 미분을 대체하는 절점력들의 적분을 사용(Silling, 2000)하기 때문에 재료 불연속성을 해석하는데 적합하다. 또한 페리다이나믹 모델은 연속체 기반 비국부 모델로서 유한한 범위 내에 흩어져 있는 모든 절점들이 상호 작용하도록 모델을 구성한다. 절점 
와 시간 
에서 페리다이나믹 운동방정식은 다음과 같다.
(1)
여기서, 
는 범위 
에 속한 두 절점 
과 
로 구성되는 페리다이나믹 결합(bond)에 대한 짝힘(pairwise force) 함수이며, 
는 변위 벡터이다. 
는 체적력, 
는 밀도이다. 본 연구에서 절점 에 대한 비국부 영역 는 Fig. 1과 같이 반경 
의 원형 영역으로 가정하였다.
(2)
미소탄성 재료(Silling, 2000)는 다음과 같이 짝힘과 선형 미소탄성 포텐셜 ω의 관계로부터 다음과 같이 정의된다.
(3)
여기서, 
, 
이다. 
에 대해 짝힘을 선형화하면 미소재료상수는 다음과 같이 정의된다.
(4)
본 연구에서는 미소탄성 취성 페리다이나믹 모델(Silling et al., 2005)을 사용하며, 이에 대한 지배방정식과 등방성 재료상수는 각각 다음과 같다.
(5)
(6)
여기서, 
이다.
3. 힘-기반 페리다이나믹-탄성체 연성 모델
이 장에서는 페리다이나믹 모델과 탄성체 모델을 힘-기반으로 연성하는 기법(Seleson et al., 2014)에 대해 소개한다. 도메인은 3개의 서로 겹치지 않는 영역으로 분할된다.
(7)
여기서, 
는 각각 페리다이나믹 모델, 탄성체 모델 영역과 두 영역이 연성되는 영역이다. 아래첨자 PD와 CE는 각각 페리다이나믹(PeriDynamic)과 탄성체(Classical Elasticity) 모델의 약자이고, 는 블랜딩 함수에 의해 연성되는 것(Blended)을 의미한다. 연성 영역은 페리다이나믹과 탄성체 영역을 다음의 블랜딩 함수를 사용하여 적절히 구성된다.
(8)
여기서, 
는 0에서 1 사이의 값을 가지는 연속다항식이다. 블랜딩 함수는 서로 다른 모델들의 기여도에 대한 가중함수 역할을 또한 수행한다. 본 연구에서 사용되는 힘-기반 연성 기법에서는 운동방정식 레벨에서 서로 다른 모델들을 결합한다. 또한 연성 영역 
주변에 비국부 영역의 반경 로 조정되는 연결 영역이 각각 
와 
에 걸쳐 존재하게 된다. 에서 연결 영역을 제외한 영역은 페리다이나믹 모델로만 구성이 된다. 반면 연결 영역이 제외된 영역에서는 변형이 부드럽다고 가정하여 탄성체 네이비어 방정식으로 표현한다. 이러한 두 영역은 연결 영역을 중간에 두어 서로 직접적으로 연관되지 않는다.
블랜딩 모델은 다음의 3단계를 거쳐 구성된다.
Step 1: 블랜딩 함수를 통한 분할
식 (8)에 소개된 블랜딩 함수를 사용하여 식 (1)의 우항을 다음과 같이 두 부분으로 분할한다.
(9)
식 (9)의 우항의 두번째 항은 식 (8)에 의해 연결 영역이 배제된 에서 사라진다. 반대로 우항의 첫번째 항은 연결 영역이 배제된 에서 사라진다.
Step 2: 선형화
에서 변형이 부드럽다고 가정했기 때문에 본 연구에서는 식 (9)의 우항의 두번째 항을 다음과 같이 선형화할 수 있다.
(10)
Step 3: Gradient expansion
부드러운 변형을 가지는 도메인 영역 내의 질점 에 대해서 역시 도메인 내의 질점 
에서의 변형을 표현하면 다음과 같다.
(11)
따라서 식 (10)은 다음과 같이 근사화될 수 있다.
(12)
항들을 정리하면 페리다이나믹 방정식과 네이비어 방정식의 블랜딩 기반 연성방정식은 다음과 같다.
(13)
4. 이산화된 2차원 모델방정식
이 장에서는 2차원 정적 문제에 대한 연성방정식을 유도하고 경계조건에 대해 논한다. 여기서 
이다.
다음과 같은 텐서를 정의한다.
(14)
여기서, 
이다. 이를 이용하여 2차원 정적 문제에 대한 블랜딩방정식은 다음과 같다.
(15a)
(15b)
식 (15a)와 (15b)의 이산화를 위해 다음과 같은 중앙차분법을 사용하였다. 균일 정규 그리드를 가정하였고, 그리드 크기는 
이며, 
이다. 정밀도는 
이다.
(16a)
(16b)
(16c)
(16d)
페리다이나믹 방정식이 작동하는 비국부 영역에 대한 경계조건 부과방식은 아직 연구가 필요하다. 본 연구에서는 네이비어 방정식이 작동하는 국부 영역에 대한 경계조건만을 고려한다. 국부 영역에 대한 Dirichlet 경계조건은 명백하므로 여기에서는 Neumann 경계조건에 대해서만 논하기로 한다. Neumann 혹은 traction 경계조건은 다음과 같은 traction 벡터의 정의에 기반한다.
(17)
Piola-Kirchoff 응력텐서 는 다음과 같다.
(18)
여기서, I는 단위텐서이고, 
은 가 포함된 경계의 단위법선벡터이다. 식 (17)과 (18)을 조합한 조건을 경계에 가한다. Neumann 경계조건을 위해 필요에 따라 의 정밀도를 가지는 단방향 유한차분(One-sided Finite Difference)법을 사용할 수 있다.
5. 수치예제
5.1 집중하중 해석
이 장에서는 2차원 수치예제를 통해 연성 모델을 검증한다. 해석 모델은 Fig. 3과 같이 가로, 세로가 동일하게 1m인 정사각영역이다. E=1Pa이고 2차원 평면응력 모델을 사용한다. 그리드 크기는 h=0.01m로 총 10,201개의 절점이 사용되었다. Homogeneous Dirichlet 경계조건이 모든 경계에 작용하였다. 이러한 문제는 국부적인 응력집중 및 심각한 변형을 초래하기 때문에 일반적인 탄성 모델로는 엄밀한 결과를 얻기 어렵다. 반면 페리다이나믹 모델은 더 엄밀한 결과를 주지만 계산비용이 탄성 모델보다 많이 비싸다는 문제가 있다. 본 연구에서 제안하는 연성으로는 집중하중이 가해지는 영역을 포함하는 주변은 페리다이나믹 모델로 구성하고 그 외의 영역은 일반 탄성 모델로 구성하므로써 계산상의 이점과 해석의 엄밀함을 동시에 추구한다. 집중하중은 다음 식과 같이 영역 중앙에서 
방향으로 부과된다.
(19)
여기서, 
는 Fig. 3에서와 같이 영역 중심이다.
계산에 사용되는 연성 모델은 Fig. 3과 같다. 연성 영역 
구성을 위해 식 (8)을 사용하고 
0.5인 상수로 정하였다. 영역 분할을 위한 파라메터는 L1=0.2이고 L2=0.8이다(Fig. 2 참고). 비국부 영역은 식 (2)와 같이 구성되며, 여기서 =4h로 정하였다.
연성 모델을 해석한 변위 결과를 Fig. 6에 나타내었다. Fig. 7의 좌측열은 CE 변위 결과를 우측열은 PD 변위 결과를 각각 비교하였다. 집중하중 영역 근처에서 CE 결과가 다른 두 변위에 비해 상당히 다른 것을 알 수 있다. 영역 중심의 하중 지점에서 최대 변위가 나타난다. 하중 지점의 수평변위 v를 비교해 보면, PD 모델 9.45×10-3, 연성 모델 8.44×10-3, CE 모델 7.33×10-3으로 나타난다. 또한 하중 지점의 수직변위 w는 PD 모델 3.27×10-1, 연성 모델 2.59×10-1, CE 모델 1.56×10-1으로 나타난다. PD 모델과 비교하여 연성 모델의 수평 변위는 약 89.5%, 수직 변위는 약 77.6%로 나타나는 반면 CE 모델의 수평 변위는 약 86.6%, 변위는 약 47.7%로 나타났다. Fig. 6과 7을 통한 해의 분포를 함께 비교하면, 연성 모델이 변위 값 뿐만 아니라 해의 거동 또한 PD 모델에 많이 근접함을 알 수 있다.
Figure 7
Displacement profiles for a point load, for the CE model(left column) and the PD model(right column)(case 1)
연성 영역 구성에 대한 영향을 검토하기 위해 Fig. 5와 같이 PD 영역이 넓은 연성 모델(L1=0.5, L2=0.8)을 구성해서 동일한 조건에서 해석을 수행하였다. 하중 지점에서 연성 모델의 수평 변위는 8.46×10-3, 수직 변위는 2.59×10-1으로 case 1과 매우 유사함을 확인하였다. 다양한 연성 영역을 사용해서 검토해 본 결과 case 1에서 사용한 좁은 PD 영역만으로도 하중 지점에 국부적으로 집중된 변형 및 응력을 엄밀하게 구현할 수 있음을 알 수 있었다. 또한 연성 모델이 PD 모델의 특성을 더 잘 표현하기 위해서는 식 (15)에서 사용하는 
값들이 PD 영역에서 정확하게 구현되어야 할 필요가 있다. 비국부 모델인 PD 모델은 경계 주변에서의 엄밀한 계산이 국부 모델보다 어렵다. 이를 해결하기 위해서는 보다 엄밀한 수치적분이 필요하다. PD 영역의 수치적분이 값에 미치는 영향은 Seleson 등(Seleson et al., 2013)이 검토한 바 있지만 이에 대해서는 추가적인 연구가 필요하다.
5.2 정적파괴 해석
Fig. 8과 같이 중앙에 0.01m의 수평 균열이 있고 인장하중 상태에 놓인 모델을 고려하자. 평판의 너비와 높이는 0.05m로 동일하며, E=192GPa이다. 그리드 크기 h=0.000355m로 총 20,164개의 절점이 사용되었으며, =0.0016m이다. 판의 상하 수평경계에 일정한 변위 경계 조건 
0, 
0.003m를 가정하였고 좌우 수직경계에는 zero traction 조건을 가정하였다. L1=0.024m이고 L2=0.038m(Fig. 2 참고)인 블랜딩 모델에 대해서 정적파괴 해석을 수행한 결과를 CE 모델 해석 결과와 비교하여 Fig. 9에 나타내었다. 균열선단 주변에서 많은 차이를 보이는 것을 확인할 수 있다. 참고로 PD 모델 해석결과는 연성
모델 결과와 거의 동일하다. 두 결과의 차이를 보다 정확하게 확인하기 위해 Fig. 10에서 수평 변위의 수평 방향 기울기를 비교하였다. 균열선단 뿐만 아니라 균열면 주변에서 CE 모델이 연성 모델(혹은 PD 모델)과 상당히 다른 경향을 나타내는 것을 확인할 수 있었다.
