1. 서 론

페리다이나믹 모델(peridynamic model)은 재료 및 거동의 불연속성을 해석하기 위해 운동방정식에서 나타나는 응력 미분을 대신하여 유한한 범위 내에 흩어져 있는 절점들의 비국부 상관관계(nonlocal interactions)를 적분하여 사용한다. 비교적 최근에 결합 기반(bond-based) 페리다이나믹 모델(Silling, 2000)이 개발된 이후 유리 및 복합재 등 취성 재료에 대한 다양한 동적파괴 현상(분기균열, 균열전파, 균열각도, 이차균열, 계단균열, 파쇄 등)이 조사되었다(Ha and Bobaru, 2010; Ha and Bobaru, 2011; Ha and Cho, 2011; Hu et al., 2011; Bobaru et al., 2012; Ha and Cho, 2012; Hu et al., 2012). 그러나 결합 기반 모델은 개별 결합들의 상호작용이 서로 독립적으로 구현되기 때문에 간단한 형태의 탄성 모델 외에 일반적이고 엄밀한 재료 모델의 구현이 어렵다. 특히 포아송비가 1/3(2차원 모델)로 한정되고 전단 변형이 표현되지 못하기 때문에 엄밀한 균열패턴 분석에 한계를 가진다.

본 연구에서는 일반화된 페리다이나믹 모델인 상태 기반(state-based) 페리다이나믹 모델(Silling et al., 2007)을 통해 동적파괴 해석을 수행한다. 상태 기반 페리다이나믹 모델은 연속체 기반의 재료 구성 모델로부터 직접 페리다이나믹 재료 모델을 구성하기 때문에 다양한 재료 모델을 엄밀하게 구현할 수 있고, 체적 및 전단 변형이 모두 표현되기 때문에 보다 정확한 모델 구성 및 파괴 해석이 가능하다. 선형 탄성체 모델의 평면응력 문제에 대한 상태 기반 페리다이나믹 모델(Le et al., 2014)을 구성하고 상태 기반 모델에 적합한 파괴 모델(Ha, 2015)을 적용하여 동적취성파괴 해석을 수행한다. 실험을 통해 조사된 균열면에 평행한 압축응력과 균열패턴의 상관관계(Ravi-Chandar and Knauss, 1984)를 페리다이나믹 동적파괴 해석 모델을 통해 검증하였다. 상태 기반 페리다이나믹 모델은 일반적인 재료 구성이 가능하기 때문에 재료 특성에 따른 균열패턴을 또한 분석할 수 있었다. 특히 재료의 포아송 비를 변경하여 압축 변형과 균열패턴의 상관관계가 보다 엄밀하게 분석되었다.

2. 평면응력 탄성체의 상태 기반 페리다이나믹 동적파괴 모델

페리다이나믹 모델은 Fig. 1과 같이 유한한 길이 δ를 반경으로 하는 비국부 영역(H_x)으로 이산화된다.

시간 t에서 절점 x에 대한 비국부 영역 내에 포함되는 모든 절점들의 상호작용력을 적분하여 페리다이나믹 운동방정식을 표현하면 다음과 같다.

는 절점 x'으로부터 x에 가해지는 힘 성분으로 단위는 힘/(체적)²이 된다. 여기서 는 힘 벡터 상태(force vector state)라고 하며, 절점 x와 연결되는 모든 절점에 작용한다. 결합 벡터 ξ=x'-x에 대해 적용시킨 는 힘 벡터를 구성한다. 벡터 상태는 벡터를 변환한다는 점에서 이차텐서와 유사하지만 상태는 선형이거나 연속일 필요가 없다. 본 연구에서는 Fig. 1과 같이 두 절점 결합의 방향에 평행하게 과 이 작용하는 Ordinary 모델을 사용한다. 이 때 두 힘 벡터의 크기는 서로 같을 필요가 없다. 변형된 절점 y와 y'로부터 변형 벡터 상태는 로 정의된다. 힘 벡터 및 변형 벡터 상태를 변형기울기(deformation gradient)로 구성되는 탄성체의 구성방정식과 연관지어 페리다이나믹 탄성체 모델을 구성할 수 있다. 페리다이나믹 상태 및 3차원 탄성체 정식화에 대한 자세한 설명은 Silling 등(2007)을 참조하기 바란다.

Figure 1

Force vector state in an ordinary, state-based peridynamic model

본 연구에서는 취성 재료로 만들어진 평면 시편에 대한 동적 균열전파 특성을 조사하기 위해 평면응력 탄성체에 대한 상태 기반 페리다이나믹 모델을 사용한다. 2차원 모델의 상태 기반 페리다이나믹 정식화(Silling et al., 2007)는 3차원 정식화와 유사하지만, 2차원 모델에서는 반경 의 비국부 영역(H_x)이 구가 아닌 둥근 디스크가 되며, 절점들은 단층의 디스크에 분포한다는 차이점이 있다. 평면응력 모델 구성을 위한 상세한 과정은 Le 등(2014)에 기술되어 있으므로 본 장에서는 간략하게 소개하기로 한다. 평면 탄성 모델의 변형에너지 밀도와 페리다이나믹 변형에너지 밀도를 대응시켜서 평면응력 페리다이나믹 모델의 힘 벡터 상태를 유도하면 다음과 같다.

여기서, 힘 스칼라 상태는 다음과 같다.

θ는 체적팽창율(volume dilatation)을 의미하며 는 편향변형율(deviatoric strain) 텐서를 의미한다. 는 결합의 변형 전 크기 를 나타내는 스칼라 상태이고, 는 결합의 길이 변화 를 나타내는 스칼라 상태이고, 편향 길이 변화 이다. 는 비국부 영역 적분 시 사용되는 영향 함수(influence function)로서 일반적으로 비국부 영역 내에서 상수 1 혹은 삼각형 형태 등으로 정의된다. 페리다이나믹 상태에 대한 내적은 다음과 같이 정의된다.

두 변형에너지 밀도를 대응시켜 페리다이나믹 계수 k'과 α를 계산하면 다음과 같다(Le et al., 2014).

k와 μ는 각각 체적 및 전단 탄성계수이다. 힘 벡터 상태를 식 (2)에 대입하면 평면응력 선형탄성체의 페리다이나믹 모델을 구성할 수 있다.

Figure 2

Schematic of fracture surface(2D)

페리다이나믹 해석에서 파괴 단면(fracture surface)은 여러 결합들이 응집되어 연쇄적으로 끊어지면서 생성된다. 이때 새로운 단위 파괴단면을 생성하는데 필요한 에너지해방율(energy release rate) G₀은 Fig. 2와 같이 0 < z < δ인 점선을 따라 존재하는 모든 절점 A이 파괴단면을 통과하여 연결되는 모든 점 B에 대해 수용할 수 있는 최대 에너지 포텐셜 ω₀을 빗금친 부분에 대해 적분하여 계산한다.

각 결합의 변형에너지 포텐셜 ω_ε이 보다 커지면 해당 결합을 끊어서 파괴 모델을 구성할 수 있다. 이때 평면응력 탄성체 상태 기반 페리다이나믹 모델에서 각 결합의 에너지 포텐셜은 다음과 같이 계산할 수 있다(Ha, 2015).

$w_{\text{ξ}}=\frac{d\mathit{\text{W}}}{d{\mathit{\text{V}}}_{\text{ξ}}}=k\text{'}\mathit{\text{θ}}\frac{2(2v-1)}{m(v-1)}\underset{-}{w}\text{<ξ>}\left|\text{ξ}\right|\underset{-}{e}\text{<ξ>}$

G₀ 값은 재료파괴 시험을 통해 실험적으로 찾을 수 있다. 평면응력 문제의 상태 기반 페리다이나믹 모델의 수치 모델 구성과정 및 파괴 모델에 대한 상세한 설명 및 수렴성 검증 등은 Ha(2015)에 수록되어 있다.

3. 균열면에 평행한 압축응력과 균열패턴의 상관관계

평면응력 탄성체의 상태 기반 페리다이나믹 동적파괴 해석 모델을 통해 균열면에 평행한 압축응력과 균열패턴의 상관관계에 대해 조사하였다. Ravi-Chandar와 Knauss(1984)은 Fig. 3과 같이 취성 재료 평판에 초기 균열면에 인장 하중 P₁을 가하면서 오른쪽 경계에 균열면에 압축 하중 를 동시에 가해서 압축응력이 균열분기 및 패턴에 미치는 상관관계를 실험을 통해 분석하였다. 인장 파괴만 발생하는 초기에는 균열이 수평 방향으로 진전하다가 일정 이상의 응력 집중(stress intensity factor)이 발생하면 분기되어 진전한다. 그러나 균열이 분기되기 전에 균열 끝단에 압축응력이 가해지면 응력 집중이 감소되어 분기 시점이 늦어지는 현상(delay of branching)이 발생한다. 압축 하중의 크기를 변화시켜가면서 실험을 통해 현상을 확인하였다(Fig. 4).

Figure 3

Experimental problem description

Figure 4

Experimental results for the influence of crack-parallel compressive stress waves(cited and modified from Ravi-Chandar et al.(1984))

Figure 5

Numerical problem description: crack-face tensile and crack-parallel compression

전산 시뮬레이션을 통해 실험결과를 검증하기 위해 페리다이나믹 해석 모델을 Fig. 5와 같이 가로, 세로가 각각 0.1m, 0.04m인 직사각형 평판으로 구성하였다. 평판의 좌측면에는 가로 방향으로 길이 0.05m의 초기 균열이 있고, 초기 균열면을 따라 P₁=24MPa의 인장 하중을 가한다. 균열면에 평행한 압축응력을 생성하기 위해 오른쪽 면에 균일한 하중 P₂를 작용시킨다. 해석 모델의 비국부 영역은 δ=0.002m인 원으로 설정하고 16,646개 균일 격자 절점으로 이산화하였다. 비국부 영역 당 최대 절점 수는 49개로 설정하였다(Ha and Bobaru, 2010).

Figure 6

Crack propagation simulations with increasing crack-parallel compression

페리다이나믹 동적취성파괴 해석에 사용된 재료는 Soda- lime 유리로 기본 물성치 및 에너지해방율 등은 Table 1과 같이 실험(Bowden et al., 1967)으로 측정된 값을 사용하였다. Fig. 6에서 압축 하중의 크기를 0, 6MPa, 12MPa로 증가시키면서 균열패턴을 비교하였다. 페리다이나믹 균열은 각 비국부 영역 내의 초기 결합의 수에 대한 끊어진 결합의 수의 비율을 손상지수로 사용하여 표현하였다. 즉, 손상지수가 0이면 초기 결합 대비 끊어진 결합이 없는 비손상 상태이고, 손상지수가 1이면 모든 결합이 끊어진 완전 손상상태를 의미한다. 따라서 초기 균열면 주변에서 손상지수는 0.4~0.5가 된다. 압축 하중을 가하지 않은 경우(Fig. 6-상)와 비교하여 6MPa의 압축 하중을 가한 경우(Fig. 6-중)는 균열분기 시점이 확연히 늦추어진 것을 확인할 수 있다. 또한 과도한 압축하중을 가하는 경우(Fig. 6-하)에는 분기균열이 발생하지 못하였다. 이는 분기균열이 응력 집중 정도에 따라 나타나는 것을 잘 표현하고 있는 결과이며, 페리다이나믹 시뮬레이션을 통해 분기 시점이 늦어지는 현상을 검증할 수 있었다. Fig. 4의 실험 결과와 Fig. 6의 페리다이나믹 해석 결과에서 압축 하중이 증가하는 경우 균열분기 시점이 늦어지는 현상과 분기된 균열이 전파하는 각도가 줄어드는 현상을 공통적으로 확인할 수 있었다.

Figure 7

Pre-notched specimen: edge tensile only

Figure 8

Crack path for change of Poisson’s ratio

균열분기 당시의 압축응력의 크기 뿐만 아니라 압축응력이 균열 끝단에 도달하는 시기 및 변형 형상의 변화 등이 함께 작용하여 이와 같이 현상이 일어난다고 판단하였고, 이를 검증하기 위해 Fig. 7과 같이 인장 하중만을 부과하는 모델에 대해 포아송 비를 0.11, 0.22, 0.44로 변경하면서 동적파괴 시뮬레이션을 수행(t_final=60μm)하였다. 더 정밀한 균열패턴 분석을 위해 이후의 예제들은 δ=0.001m인 원형 비국부 영역 내에 최대 49개의 절점이 분포하도록 설정하고 65,448개의 균일 격자 절점으로 이산화하였다. Fig. 8에 도시한 바와 같이 인장력이 수직 방향으로 균일하고 고정된 인장변형 a를 발생시킨다고 가정했을 때, 포아송 비가 작을수록 수평 방향 압축변형 b가 크고 압축응력이 빠르게 전파될 것을 예측할 수 있다. 포아송 비가 가장 작은 0.11의 경우(Fig. 8-상)와 포아송 비가 0.22의 경우(Fig. 8-중)를 비교하면 균열분기 시점이 눈에 띄게 늦어지거나 분기 각도가 크게 줄어드는 경향은 없었다. 그러나 포아송 비가 가장 큰 0.44의 경우(Fig. 8-하)와 비교하면 분기균열 시점이 0.005m 이상 늦어지는 것을 확인할 수 있었다(수직 실선과 점선 비교). 분기 각도도 포아송 비가 0.11인 경우에 약간 작아졌지만 이후 분기균열이 진전되는 각도에는 크게 변화가 없었다. 이를 통해 균열분기 시점은 압축응력파가 도달하는 시점 및 변형 형상이 상대적으로 주요하게 작용한다고 예상할 수 있다. 분기된 균열이 진전하는 각도와 압축응력의 크기의 관계를 조사하기 위해 Fig. 9와 같이 균열 면에서 떨어진 위치에 고정된 수직 인장력 f_y와 변화하는 수평 압축력 f_x를 동시에 가하면서 균열패턴을 Fig. 10에 비교하였다. 인장과 압축 응력파가 간섭하는 시점이 비슷하기 때문에 균열분기 시점은 4가지 경우 모두 비슷하였다(수직 실선). 그러나 지속적으로 높은 압축력을 받을수록 분기된 균열이 진전하는 각도가 좁아지고 동일 시간(60μm)동안 진전하는 거리도 다소 짧아지는 것을 확인할 수 있었다.

Figure 9

Pre-notched specimen: edge tensiles and compressions

Figure 10

Crack path for increasing crack-parallel compression under edge tensile

4. 복합적인 계단균열 및 2차 균열분석

본 논문에서 구성하는 균열진전 문제는 인장모드(K-I)가 균열진전에 지배적으로 작용하는 인장 균열 모델이다. 이 때 인장력의 크기가 균열패턴에 미치는 영향을 조사하기 위해 Fig. 7과 같이 인장 하중만을 받는 평판 모델에 대해서 하중의 크기를 12, 24, 36, 48MPa로 변화시키면서 균열패턴을 분석하였다. Fig. 11에서 확인할 수 있듯이, 12MPa 인장하중 모델의 결과를 기준으로 하중의 크기가 2배(24MPa)로 증가하면 균열분기가 매우 일찍 발생하며 분기각도가 절반 이하로 작아진다. 추가적인 분기균열도 발생하지만 선행 분기균열이 지배적인 것을 알 수 있다. 그러나 하중이 3배(36MPa)로 증가하면 추가적인 분기균열이 두드러질 뿐만 아니라 다수의 2차 균열(secondary crack)이 발생하게 된다. 여러 가지 이유에 의해 기존 균열면 일부에 응력집중이 발생하면 집중된 에너지를 해방하기 위해 기존의 균열면에 거의 수직으로 2차 균열이 발생하게 된다. 2차 균열은 인장균열(tensile crack)만으로는 에너지를 적절히 해방할 수 없는 경우에 전단모드(K-II)가 지배적으로 작용하면서 발생한다. 2차 균열과 탄성에너지의 관계는 Ha와 Bobaru(2011)에서 연구된 바 있다. 인장하중이 매우 큰 경우(48MPa)에는 초기 균열로부터 분기된 균열이 진전하면서 동시에 오른쪽 면 중앙에서 새로운 균열이 발생하여 서로 반대 방향으로 전파한다. 새로운 균열은 선행 분기균열에 가까워지면서 둘로 갈라져서 유착(coalescence)된다. 연쇄적이고 동시다발적으로 나타나는 균열들이 최종적으로 복잡한 계단균열(cascade cracks) 패턴을 형성한다. 결합 기반 페리다이나믹 모델을 암석 재료에 적용하여 외력의 증가에 의한 균열전파, 분기, 2차 균열, 계단균열, 유착 등의 현상을 조사하기도 하였다(Ha et al., 2014).

Figure 11

Crack path for change of loading magnitude (Soda-lime glass)

5. 재료 특성에 따른 균열패턴 분석

재료 특성에 따른 균열패턴을 분석하기 위해 Duran 50 유리와 Soda-lime 유리에 대한 동적파괴 해석 결과를 비교하였다. Duran 50 유리의 물성치들은 실험(Döll, 1975)을 통해 규명된 사용하였다. Duran 50 유리는 열 및 화학적 내성이 증대된 강화유리로서 Table 2에서 확인할 수 있듯이 구조강성(Young’s modulus)은 Soda-lime 유리에 비해 다소 작으나(83%) 에너지해방율이 매우 높다(151%). Fig. 7과 동일한 모델에 대해 Duran 50 유리 재료를 사용하고 인장 하중의 크기를 12Mpa, 24MPa, 36MPa로 조정하여 동일 시간(60μm)에서 균열패턴을 조사하였다(Fig. 12).

Figure 12

Crack path for change of loading magnitude (Duran 50 glass)

Soda- lime 유리를 사용한 균열패턴(Fig. 11)과 비교하여 12MPa 인장 하중의 경우 두 유리 시편 모두 1번만 균열이 분기되었으나 Duran 50 유리가 균열이 더 느리게 진전한 것을 볼 수 있다. 또한 24MPa과 36MPa의 두 인장하중의 경우 Duran 50 유리 시편은 Soda-lime 유리 시편에서 균열이 분기된 시점을 한 차례 견디고 0.015m 이상 더 진전한 후 분기균열이 발생하였다. 전체적인 균열 양상도 Duran 50 유리 시편에서 완화되는 것을 확인할 수 있다. 이와 같은 분석을 통해 Soda-lime 유리에 비교하여 강성(Young’s modulus)은 작지만 에너지해방율이 큰 Duran 50 유리가 파괴 내성이 더욱 우수한 것을 확인할 수 있었다.

Fig. 13에서 두 유리 재료에 대한 2차 균열패턴을 또한 검토하였다. 문제 구성은 Fig. 7과 동일하고 재료만 각각 Soda-lime과 Duran 50 유리를 사용하였다. 12MPa 인장 하중에서 선행 분기된 균열이 충분히 진전한 후 85μm에서 균열패턴을 비교하였다. Soda-lime 유리 시편은 선행 균열 분기 직후 추가로 분기되려다가 구속(arrested crack)된 지점(Fig. 11-상 참조)에서 2차 균열이 발생하였다. Duran 50 유리 시편은 선행 균열분기 이후 균열구속 등 별도의 취약지점이 발생하지 않았고(Fig. 12-상 참조), 선행 균열전파 각도가 변경되는 균열 선의 중앙 부근에서 에너지가 집중된 후 2차 균열이 발생하였다.

Figure 13

Secondary crack patterns for different glass materials

6. 결 론

본 논문에서는 평면 시편에 대한 균열전파, 분기균열 패턴과 응력파 및 변형 형상에 대한 영향을 조사하였다. 동적파괴 전산 시뮬레이션을 위해 평면응력 탄성체에 대한 상태 기반 페리다이나믹 모델을 적용하고, 상태 기반 모델에 적합한 손상 모델을 적용하였다. 실험을 통해 관측된 균열면에 평행한 압축응력과 균열패턴의 상관관계(Ravi-Chandar and Knauss, 1984)를 전산 시뮬레이션을 통해 분석하였다. 재료 특성을 정확히 표현할 수 있는 상태 기반 페리다이나믹 모델을 활용하여 재료 물성치를 적절히 변화시켜가면서 실험 결과에 대한 심층적인 분석이 가능하였다. 하중의 크기에 따라 복합적인 계단균열 및 2차 균열이 발생하는 것을 또한 확인할 수 있었고 재료 특성에 따른 균열패턴이 타당하게 모사되는 것을 볼 수 있었다. 본 연구를 통해 상태 기반 페리다이나믹 모델이 동적파괴 거동 분석 및 엄밀한 손상 해석에 유용한 것을 확인하였다.