연성 재료의 파단 변형률을 지배하는 변수는 식 (1)에 나타낸 응력 삼축비 (𝜂)와 식 (2)에 나타낸 로드각 파라메타 (θ)로 알려져 있다. 여기서 사용된 I1, J2, J3는 각각 응력 텐서 (σ)의 1차 불변량, 편차 응력 텐서 (s)의 2차 불변량과 3차 불변량이다 (식 (3),(4),(5) 참조). 평면 응력 상태에서 응력 삼축비와 로드각의 관계는 식 (6)과 같이 변경된다.

Fig. 1은 하중의 유형 또는 응력의 유형에 따른 응력 삼축비 및 로드각 파라미터의 변동성을 보여 준다. 응력 삼축비와 로드각 파라미터 조합 (η, θ)이 (0.0, 0.0)일 경우 순수 전단, (-1/3, -1.0)일 경우 일축 압축, (1/3, 1.0)일 경우 일축 인장, (-1/3, 0.0)일 경우 평면 변형률 압축, (-1/3, 0.0)일 경우 평면 변형률 인장 등이 대표적인 하중 케이스이다.

Fig. 1 
Stress triaxialities and Load angle parameters according to loading types

Mohr & Marcadet (2015)는 비례 하중 (proportional loading)을 받는 연성 재료의 파단 예측을 위하여 응력 삼축비와 로드각 파라미터를 변수로 가지는 Hosford-Coulomb (HC) 파단 모델을 제시했다. 식 (7)에 보인바와 같이 HC 파단 변형률 ε¯HCpr는 네가지 미지수 (b, c, g, nf)를 가진다. 식 (8)은 미지수 g를 결정한다. 또한 식(9),(10),(11),(12)는 식 (8)에 사용된 미지수를 결정한다. σHF는 HC 등가 응력이며, f1, f2, f3는 로드각 파라미터로부터 결정되므로, 재료 상수는 4가지 (a, b, c, nf)로 압축된다. Roth & Mohr (2016)는 일반 강재에 대하여 nf=0.1을 제안하였으며, 이 경우 HC 파단 변형률 모델은 3가지 재료 상수 (a, b, c)로 구성 된다. 여기서 c는 Mohr-Coulomb 항복 함수의 내부 마찰각 계수 (internal friction coefficient)와 동일하다. 식 (9)에서 σ1, σ2, σ3는 각각 주응력을 의미한다.

b는 HC 파단 변형률의 평균 높이를 결정하며, 등가 이축 인장 실험 (equi-biaxial tension) 또는 일축 인장 실험 (uniaxial tension)을 통해 얻어진다. Fig. 1는 평면 응력 상태에서 a와 c의 변동성에 따른 HC 파단 변형률의 차이를 나타낸다. 일축 인장 응력 상태의 η는 0.33이므로, a와 c는 이를 기점으로 파단 곡선의 포물선을 결정한다.

최종 파단에 이르는 응력 경로 또는 하중 경로의 효과를 고려하기 위하여 식 (13)과 같이 파단 변형률의 누적으로 표현되는 누적 손상 (D)을 적용하였다. 여기서 m은 경로를 결정하는 지수로서 재료 상수이다. m=1.0이면 선형 누적을 의미한다. 본 연구에서는 문제의 간략화를 위하여 m=1.0으로 가정하고자 한다. 손상의 초기 값은 항상 0.0이며, 재료에서 응력과 소성 변형률 이 발달함에 따라 D=1.0일 때 파단이 발생하는 것으로 간주한다.

Fig. 2 
Effect of the material constants a and c of the HC fracture strain model for plane stress state

평판 성형 해석에 적용되는 성형 한계도 (forming limit diagram, FLD)의 개발에 Marciniak-Kuczynski (MK) 분석이 널리 활용되어 왔다 (Banabic, 2010). Pack & Mohr (2017)는 단위 크기의 셸 요소 모델에 대한 MK 분석을 통해 2축 인장을 받는 박판 구조물의 파단 변형률 모델 (DSSE)를 제시하였다.

DSSE는 박판 구조물에서 두께 방향 국부 네킹이 발생하는 소성 변형률과 파단 변형률의 차이가 미소할 것으로 가정하고, 두께 방향 네킹 소성 변형률을 파단 변형률로 적용한 모델이다. DSSE는 2축 응력비가 유지 (1/3≤η≤2/3)되고 비례 멤브레인 하중 (proportional membrane loading)일 경우에만 적용될 수 있다. 본 연구에서는 응력 삼축비가 이 구간을 벗어날 경우 HC 파단 변형률을 적용하였다.

DSSE 파단 변형률 (ε¯DSSEpr)은 식 (14)과 같으며, 여기서 g1 및 g2는 응력 삼축비의 함수로 식 (15),식(16)과 같다. 따라서 확인해야할 재료 상수는 3가지 (b, d, p)이다. 재료 상수 b는 HC 파단 변형률 모델에서 사용되었으므로, 재료 상수는 총 2가지 (d, p)로 압축된다. Pack & Mohr (2017)은 MK 해석을 통하여 일반 강재의 p=0.01을 제시하였다. 따라서 DSSE 파단 변형률 모델에서 확인되어야할 재료 상수는 1가지 (d)가 된다.

Considère 가정에 따르면 소성 불안정은 국부 네킹 (localized necking)이 발생하는 시점에서 시작한다. 이점에 착안하여 평면 변형률 인장 (plane strain tension, PST)을 받을 경우 국부 네킹 소성 변형률 (ε¯DSSEPST)에서의 유동 응력 (kε¯DSSEPST)은 식 (17)과 같다. 평면 변형률 인장 하중을 경험할 경우, 응력 삼축비는 η=1/3이고 이를 식 (14)에 대입하면 식 (14)는 식 (18)과 같이 반복 수치 해석을 통하여 해를 도출할 수 있는 형식으로 변경된다. 선형 누적 손상을 가정할 경우, 하중 경로 효과를 고려한 파단 조건은 식 (19)와 같다.

선행 연구에서 수행된 파단 실험(Choung et al., 2011; Choung et al., 2015a; Choung & Nam, 2013)의 시편 형상은 Fig. 3에 나타나 있다. 시편은 25mm 두께의 EH36 모재의 두께 방향으로 중층 (middle layer), 가공 방향으로 직교하게 시편이 제작되었다. 시편의 형상과 명칭은 Table 1과 같다.

Fig. 3 
Fig. 3 Geometry of specimens (unit: mm)

FB 및 RB는 재료의 유동 응력 (flow stress)를 얻기 위하여 제작되었으며, FNT, RNT, ST 시편은 파단 모델의 재료 상수를 교정하기 위하여 제작되었다. ST의 경우 전단 및 인장 하중이 복합적으로 작용하도록 총 3가지 회전각 (15°, 30°, 45°)을 가지도록 시편이 제작 되었다.

RB, FB, FNT, RNT 시편에는 50mm 신률계 (extensometer)를 부착하여 신률계 변위와 로드셀 하중을 계측하였다. ST 시편 실험에서는 실린더 변위와 로드셀 하중을 계측하였다.

Fig. 4에 시편별 유한 요소 모델을 제시하였다. FNT, RNT 시편의 경우 신률계가 부착된 지점까지 유한 요소를 생성했다. FNT 시편의 경우 시편의 길이, 두께 및 폭 방향 대칭을 감안하여 시편의 1/8을 감차 적분 8절점 요소 (C3D8R)로 모델링한 후 대칭 경계 조건을 부여하였다. RNT 시편의 경우 시편 길이 방향 대칭을 고려하여 시편의 1/2을 감차 적분 4절점 축대칭 요소 (CAX4R)로 모델링 한 후 대칭면에 대칭 경계 조건을 부여하였다. ST 시편의 경우 시편 두께 방향 대칭을 고려하여 시편의 1/2을 감차 적분 20절점 고체 요소 (C3D20R)로 모델링 한 후 대칭면에 대칭 경계 조건을 부여하였다. 모든 시편에 인장 방향 강제 변위를 모델 경계부에 부여하여 인장 하중을 구현하였다. 시편의 길이 방향 중앙부에는 요소 크기에 따른 수렴도 테스트를 통하여 노치 반지름과 요소 길이의 비가 0.1 이하가 유지되도록 하였다. 수치 해석에 사용된 재료의 탄성 계수 (E), 포아송비 (ν) 및 밀도 (ρ)는 206GPa, 0.3, 7.85ton/m³이다. 수치 해석에 사용된 재료의 유동 응력은 3.3절의 과정을 통하여 도출하였다.

Fig. 4 
Numerical simulation models for smooth and notched specimens

상용 유한 요소 해석 코드 Abaqus/Standard (Simulia, 2018)을 이용해 비선형 유한 요소 해석을 수행하였다.

등방성 재료의 von Mises 항복 함수는 재료에서 발달한 von Mises 등가 응력 (σ¯)과 재료의 유동 응력 (k)으로 구성된다 (식 (20) 참고). Fig. 5은 FB 시편을 이용하여 얻어진 균일 진 응력 – 균일 진 소성 변형률 선도이다. 초기 항복 응력 (σ¯0, Initial yield stress)이 ε¯plat까지 유지되는 항복 평탄부 (yield plateau)가 존재하는 것을 확인할 수 있다. 따라서 본 연구에서는 식 (21)와 같이 항복 평탄부 전후의 유동 응력을 나누어 고려하였다. 평탄부 이후의 유동 응력은 소성 불안정 이후 응력 삼축비 효과를 고려하기 위하여 식 (22)의 Swift 구성 방정식과 식 (23)의 Voce 구성 방정식의 선형 조합인 Swift-Voce 구성 방정식으로 표현되었다. Swift-Voce 구성 방정식은 Swift 구성 방정식의 재료 상수 {A, ε0, n}, Voce 구성 방정식의 재료 상수 {k0, Q, β}로 구성된다. α는 구성 방정식의 가중치이다.

Fig. 5 
Flow stress of FB

평활재 (FB 및 RB) 실험으로부터 얻어진 균일 진 응력-균일 진 변형률 곡선 (네킹 이전 유동 응력)에 Swift 구성 방정식과 Voce 구성 방정식을 적합함으로써 각각의 재료 상수를 결정하였다. α를 변경하면서 노치재 (FNT 및 RNT)에 대한 유한 요소 해석을 수행하고, 노치재의 힘-변위 곡선이 실험과 최소 오차를 가지도록 α를 도출하였다. Table 2는 실험 데이터에 가장 잘 적합하도록 결정된 재료 상수를 정리한 표이다.

각각의 구성 방정식을 적용했을 때 유동 응력을 Fig. 6에 제시하였으며, Swift-Voce 모델을 적용할 경우 소성 변형률이 증가함에 따라 응력 삼축비로 인한 강재의 가공 연화가 구현되었음을 확인 가능하다. 수치 해석으로부터 얻은 노치재의 하중-변위 곡선을 실험과 비교하여 Fig. 7에 나타내었다. 이로부터 재료 상수 α가 적절하게 도출되었음을 확인 가능하다.

Fig. 6 
Comparison of flow stresses

Fig. 7 
Comparison of force-displacement curves for FNT and RNT specimens

실험을 통하여 재료 내부의 응력 및 변형률을 계측할 수 없기 때문에, 본 연구에서 수행한 수치 해석은 재료 내부, 즉 파단이 발생하는 지점에서의 응력 및 변형률을 도출하기 위하여 수행된 것이다. 이러한 이유로 파단 조건을 부여하지 않은 비선형 수치 해석을 실시하였다.

Fig. 8에는 실험 및 수치 해석에서 얻어진 신률계 변위-힘 곡선을 제시 하였다. 대부분 시편의 시뮬레이션으로부터 얻은 신률계 변위-힘 곡선이 실험 결과와 상당히 잘 일치하는 것을 확인할 수 있다. 여기서, 실험으로부터 얻어진 하중-변위 곡선의 마지막이 파단으로 간주될 수 있으며, 하중이 급격하게 감소하는 지점을 파단 발생 시점으로 간주하였다.

Fig. 8 
Force – displacement curves based on Swift-Voce constitutive model

Fig. 8에는 파단 시점에서 유한 요소 해석의 변형 형상과 파단 등가 소성 변형률 (ε¯p,f)을 동시에 나타내고 있다. 이를 바탕으로 파단의 발생 지점을 확인하였으며, 이 파단 발생 지점에 대한 응력 삼축비 및 로드각 파라미터의 발달 과정이 Fig. 9에 제시되었다.

Fig. 9 
Evolution of the equivalent plastic strain against two fracture variables

대부분의 시편에서 등가 소성 변형률의 증가에 따라 η도 선형적으로 증가하는 것을 확인 가능하다. 반면, 등가 소성 변형률의 증가에 따라 θ는 서로 다른 양상을 보인다. 예를 들어 RNT 시편의 로드각 파라메타는 1.0을 유지하고 있는 반면, FNT 시편의 θ는 지속적으로 감소하고 있다. ST 시편의 θ는 변곡점을 가지는 것으로 나타났다. η와 θ가 일정하게 유지되는 경우에만 비례 하중으로 간주될 수 있기 때문에 RNT 시편만 약간의 비례 하중을 경험하였으며 (엄밀하게는 비비례 하중), 나머지 시편의 경우 모두 비비례 하중 경로를 나타내었다.

HC 파단 변형률 모델의 재료 상수 집합 αa, b, c를 결정하기 위해서, 최소 자승법 (least square method)를 적용하였다. 즉, 식 (24)와 같은 오차의 제곱 (R)을 최소화하는 α를 탐색하기 위하여 MATLAB의 FMINUNC 함수를 이용하였다 (식 (25) 참조). 여기서 i는 재료 상수 교정에 사용된 실험의 횟수이며, 본 연구에서는 5회의 실험 데이터 (RNT, FNT, ST15, ST30, ST45)가 사용되었다.

최적화를 통해 결정된 파단 모델 재료 상수는 a=1.8648, b=1.5734, c=0.0278이다. 이를 이용하여 가시화하면 Fig. 10과 같은 파단 변형률 평면을 얻을 수 있다. EH36 강재의 파단 변형률 평면은 로드각 파라미터 보다는 응력 삼축비 변동성에 민감한 것을 확인할 수 있다. 식 (17) 및 (18)을 이용하여 결정된 DSSE 파단 변형률 모델의 재료 상수는 ε¯DSSEPST=0.1944 및 d=1.6688이었다. Fig. 11에는 1/3≤η≤2/3이 만족되는 평면 응력 (plane stress) 조건에서 HC와 DSSE를 조합한 파단 변형률 곡선이 제시되었다. 이 구간에서 DSSE 파단 변형률은 HC 파단 변형률에 비하여 현저하게 낮은 것을 확인할 수 있다.

Fig. 10 
HC fracture strain surface for EH36 steel

Fig. 11 
HC and DSSE combined fracture loci for plane stress condition

HC-DSSE 파단 조건을 수치 해석에서 구현하기 위하여 상용 유한 요소 해석 코드 Abaqus/Explicit의 재료 사용자 서브루틴 (VUMAT)을 개발하였다. 사용된 컴파일러는 포트란 (Fortran) 이었다. 유한 요소의 적분점에서 계산된 응력 성분과 변형률 성분을 사용자 서브루틴을 통하여 HC-DSSE 파단 변형률 조건을 구현하였다. 요소의 적분점에서의 누적 손상 D가 1.0에 도달하면 요소를 제거함으로써 파단을 구현하였다.

본 연구에서는 Choung et al. (2015b)에 의하여 수행된 비대칭 노치재 시편의 인장 실험을 검증 모델로 결정하였다 (Fig. 12 참조). 고체 요소 (C3D8R)와 셸 요소 (S4R) 두 종류의 유한 요소 해석 모델이 준비되었다. 신률계가 부착된 50mm 지점까지 모델링되었다. 고체 요소 모델의 경우 요소 크기와 노치 반지름의 비가 0.1 이하가 되도록 유지했다. 고체 요소 모델의 경우 두께 및 길이 방향으로 대칭성을 고려하여 1/4 대칭 모델이 사용되었다. 여기서 고체 요소는 매우 조밀한 크기인 0.1mm를 적용하였으며, HC 파단 변형률을 적용한 케이스를 SD-HC로 명명하였다 (Table 3 참조). 셸 요소 모델의 경우 두께 방향 형상이 없으므로 1/2 대칭 모델이 사용되었다. 단, 셸 요소 모델의 경우 5가지 요소 크기 (0.1mm 0.2mm, 0.4mm, 0.6mm, 0.8mm)를 가지는 해석 모델이 준비되었다 (Fig. 13 참조).

Fig. 12 
Design and photo of asymmetric specimen (Choung et al. 2015b)

Fig. 13 
Numerical simulation models for an asymmetrically notched specimen

하중 지지 능력의 감소가 파단 발생 여부를 결정하는 지표로 예측되었기 때문에 본 연구에서는 하중 지지 능력의 급격한 감소, 즉 하중 – 변위 곡선으로부터 하중이 급격히 감소하는 변위를 파단의 발생 시점으로 판단하였다. EH36 강재의 경우 손상의 시작 (damage initiation) 이후에 손상의 발전 (damage evolution)까지의 과정이 포함되었기 때문에, 실험으로부터 얻은 하중 – 변위 곡선에서 하중의 급격한 변화 시점을 정확하게 파악하는 것은 매우 어려운 일이다. 따라서 본 연구에서는 비교적 파단 발생 변위를 확장 예측하여 Fig. 14 및 Fig. 15에 제시하였다. 실제 파단은 좀 더 작은 변위에서 시작되었을 가능성을 배제할 수 없다는 뜻이다.

Fig. 14 
Fracture simulation results for solid elements with HC fracture strain model

Fig. 15 
Fracture simulation results for shell elements

이렇게 정확한 파단의 발생 변위 또는 시점을 정확히 파악하기 어렵기 때문에 수치 해석을 통한 파단 예측의 정확도를 논의하는 것도 매우 어려운 문제가 된다. 따라서 본 연구에서는 수치 해석에서 얻은 하중이 급격하게 감소하는 시점 (변위)과 이때 누적 손상이 1.0에 근접한지의 관점에서 파단 예측의 정확도를 논의하고자 한다.

실험에서 얻은 인장 하중 – 변위 선도 및 HC 파단 변형률과 0.1mm 요소 크기를 가지는 고체 요소 수치 해석 (SD-HC)에서 얻은 하중 – 변위 선도 및 누적 손상 – 변위 선도가 Fig. 14에 제시되었다. 하중의 관점에서 실험과 수치 해석이 전 변위에 걸쳐 비교적 잘 일치하는 것을 확인 가능하다. 누적 손상 1.0에 상응하는 변위보다 약간 더 증가한 변위에서 파단이 발생한 것으로 관찰된다.

HC 파단 변형률 모델 및 HC-DSSE 파단 변형률 모델을 다양한 요소 크기를 가지는 셸 요소 수치 해석에 적용하여 얻은 하중 – 변위 선도 및 누적 손상 – 변위 선도를 Fig. 15에 제시하였다. 여기서 Fig. 15(a) 및 Fig. 15(b)는 각각 HC 파단 변형률 모델 및 HC-DSSE 파단 변형률 모델을 적용한 결과이다. 요소 크기가 감소함에 따라서 하중이 감소하는 시점 (파단에 이르는 시점)이 빨라지는 것을 확인할 수 있다. HC 파단 변형률을 셸 요소에 적용할 경우 파단 발생 시점에서 누적 손상이 0.4 내외에 불과하다는 것을 Fig. 15(a)로부터 확인 가능하다. 누적 손상이 1.0에 도달하도록 시뮬레이션을 실시하였다면, 매우 과도한 파단 변위를 유발할 것이 예측된다. 이러한 이유로 HC 모델을 셸 요소에 직접 적용할 경우 너무 늦은 파단을 또는 너무 큰 변위에서 파단을 예측하게 되는 것이다. 반면 HC-DSSE 파단 변형률 모델을 적용할 경우, 누적 손상이 1.0에 도달했을 때 비로소 급격한 하중의 감소와 파단이 발생하는 것을 Fig. 15(b)로부터 확인 가능하다. 특히 SH04-DS 케이스의 파단 발생 변위는 SD-HC의 그것과 거의 유사한 것을 확인할 수 있다. 셸 요소 크기에 따른 파단 시점의 변동성 측면에서도 HC-DSSE를 적용한다면 셸 요소 크기에 비교적 덜 민감한 파단 예측이 가능함을 Fig. 15(b)로부터 확인 가능하다. 그러나 요소 크기 민감도에 대하여 좀 더 상세한 연구가 요구된다.

Fig. 16은 HC-DSSE 파단 변형률을 적용해서 누적 손상이 1.0에 도달하였을 변위에서 누적 손상을 보여준다. Fig. 15에서도 보인 것처럼 HC 파단 변형률 모델을 적용시에는 대략 0.4 정도의 누적 손상만을 보여 주고 있다. 그러나 손상의 최대값이 발생하는 위치는 거의 파단 변형률 모델에 따라 거의 동일한 것으로 관찰되었다.

Fig. 16 
Damage distributions when maximum accumulated damage reaches 1.0