1. 서 론

유한요소해석은 금속재료의 소성변형을 예측하고 해석 하기 위하여 널리 사용되고 있다. 유한요소해석에서 재료 의 소성 변형률 분포와 응력분포를 계산하기 위해서는 재 료의 비선형 방정식을 하중 경로를 따라 풀어야 한다. 변 위에 대한 비선형 방정식의 해를 구하기 위하여 회귀착점 방법(return-mapping method)이 널리 사용되고 있으며, 이 를 구현하기 위하여 뉴튼-랩슨 축차(Newton-Raphson iteration) 방법이 널리 사용되고 있다.

3-D 적층 제조 공정으로 제작된 재료는 하중 방향에 따 라 강한 이방성을 나타낸다고 보고되고 있다[1, 2]. 이러한 이방성 금속재료의 소성변형을 해석하기 위해 이방성 항 복 함수가 사용되며, Hill의 이방성 2차 항복함수가 그 중 한 예이다. 그러나, 금속재료가 강한 이방성을 나타내거나 항복함수가 모서리를 가진다면, 뉴튼-랩슨 축차는 수렴하 지 않는 경우가 발생한다[3]. 이러한 경우를 해결하기 위 하여 Scherzinger[4]는 선분탐색 기법(line search method) 을 뉴튼-랩슨 방법과 연동하여 사용할 것을 제안하였다. 그러나 선분탐색 기법을 뉴튼-랩슨 축차에 추가할 경우, 계산 식과 알고리즘이 매우 복잡해지므로 이를 코드로 구 현하는 것은 매우 까다로운 작업이다.

본 논문에서는 소성 변형의 유한요소 해석에서 뉴튼-랩 슨 방법을 사용하지 않고 단순한 반복계산 만을 사용하여 소성변형의 해를 구하는 고정점 축차(Fixed-point iteration) 기법을 제안하였다. 고정점 축차는 고전적인 수치해 석 분야에서 비선형 방정식의 해를 구하기 위하여 널리 사용되어 왔다. Akpama 등[5]은 고정점 축차 방법을 결정 소성학의 활주 량을 구하는데 사용하였으며, 계산 속도를 뉴튼-랩슨 기법과 비교하였다. 고전 소성변형의 해석에서 고정점 축차를 적용하여 문제를 해결한 연구는 시도된 바 가 없다. 본 논문에서 시도한 고정점 축차 기법은 항복함 수의 2계 도함수를 계산할 필요가 없으며, 반복계산의 수 렴범위가 명확히 제시되므로 소성변형의 해를 구하는데 매우 유용하리라 판단된다.

다음 장에서는 Hil l의 이방성 2차 항복함수에 대하여 간 략히 기술하였고, 소성 변형에서의 고정점 축차 기법에 대 하여 기술하였다. 그 후 고정점 축차를 이용하여 3차원 유 한요소 해석을 수행하였고, 계산시간을 뉴튼-랩슨 방법과 비교하여 고정점 축차의 우월성을 증명하였다.

2. Hill의 이방성 2차 항복함수

금속분말을 3D Printing 기법으로 제작한 구조물은 적층 방향에 따라 이방성(anisotropy) 을 나타낸다고 보고되어 있다[2]. 이러한 이방성이 강한 금속 소재의 소성변형을 해석하기 위하여 본 논문에서는 Hil l의 2차 항복함수를 사 용하였다. 항복 함수의 형태는 다음 식과 같다[6].

(1)

2ϕ2(σ)=(G+H)(σxx−σzz)2−2H(σxx−σzz)(σyy−σzz)+(F+H)(σyy−σzz)2+2Nσxy2+2Lσyz2+2Mσzx2=2σf2

F, G, H, L, M, N은 물질의 이방성에 의존하는 재료 상 수이고 σ_f는 현재의 항복강도를 나타낸다. 본 논문에서는 등방성 변형경화를 고려하였다. σ_ij는 실험실 좌표계에서 의 응력의 성분이다. 만약 F = G = H = 1, L = M = N = 3이 라면, 항복함수는 von Mises의 등방성 항복함수로 귀결된 다. 본 논문에서는 강한 이방성을 나타내기 위하여 참고문 헌[2]의 이방성 재료를 모델링한 참고문헌[7]에서 사용되 었던 재료 상수를 그대로 사용하였다. 이 경우 H il l의 항 복함수는 다음 식과 같다.

(2)

(σxx−σzz)2−74(σxx−σzz)(σyy−σzz)+(σyy−σzz)2+3τxy2+3τyz2+3τzx2=σf2

따라서 F = G = 1, H = 7 / 4에 해당한다. (2)식에 제시된 이방성 2차 항복함수의 초기 항복면과 변형 경화된 이후 의 항복면을 편차응력 공간에 도시하였다(Fig. 1 참조). σ_f = 1은 항복강도가 1MPa 인 경우이고 σ_f = 2는 항복강 도가 2MPa 로 증가된, 즉 변형경화 된, 경우이다. π-평면 상에서, von Mises 항복조건일 경우는 잘 알려진 바와 같 이 원형의 항복곡선을 나타내고 이방성일 경우 타원을 나 타낸다. 본 논문에서는 Voigt 표기법을 사용하여 응력과 변형률은 벡터 V=(σxx,σyy,σzz,τxy,τxz,τyz) 등으로 나타낼 수 있고 강성은 강성행렬 C로 나타낼 수 있다.

Fig. 1

Yield locus of Eq. (2) on the π-plane (before (σ_f = 1) and after (σ_f = 2) the strain hardening).

3. 소성변형에서의 고정점 축차 및 수렴조건

소성변형 하는 금속 구조물에서 변형률 분포와 응력 분 포를 계산하기 위해서는 항복조건(yield condition), 유동법 칙(flow rule), 응력과 탄성변형률 간의 탄성구성방정식을 적용하여야 한다. 보통 회귀착점 방법(return-mapping method)을 사용하여 항복면과 이에 대응하는 응력, 소성변 형률 상태 등을 구하게 된다. 본 논문에서는 최근에 제시 된 고정점 축차(fixed-point iteration) 기법을 사용하여 소 성변형 해석을 수행하고자 한다[8].

수학에서 고정점이라 함은 함수에 의해 그 값이 변하는 않는 점이다. 즉 고정점에 대하여 다음 식이 성립한다.

x=f(x)

위 식에서 x를 고정점이라 부른다. 비선형 방정식의 해 를 구하는 과정도 일종의 고정점을 찾는 과정으로 간주할 수 있으며, 이를 고정점 축차라고 부른다.

소성변형에서 어떠한 시간 구간 동안의 등가소성변형률 증분은 다음의 비선형 변형경화 식을 만족시킨다.

nCnΔε¯p=nCΔε−hΔε¯p

n은 소성유동의 방향을 나타내고 C는 탄성 계수를 나타 내고, Δε는 주어진 시간 스텝 n과 n + 1 사이의 변형률 증 분을 나타낸다. 이 식은 소성변형 하는 재료의 항복함수의 정합조건(consistency condition)을 나타낸다. 위 식의 해는 Δε¯p 가 된다. 위의 비선형 방정식에 고정점 축차기법 (fixed-point iteration method)을 적용한다면, 등가 소성 변 형률의 증분에 대하여 다음과 같은 축차 방정식(iteration equation)을 적용하게 된다.

(3)

nCnΔε¯p(k+1)=nCΔε−hΔε¯p(k)

여기서 (k)는 현재 축차 회수를 나타내며, 통상적으로 변 형율 증분은 주어진 값이라고 가정한다. (3)식의 기본원리 는 우변에 주어진 이전 축차의 Δε¯p 값으로부터 좌변에 주어진 현재 축차의 미지수 값 Δε¯p 을 구하는 것이다. 종 속 변수의 값들은 축차 과정 중의 임시 상태에서의 값을 나타낸다. 위의 고정점 축차가 수렴하기 위한 조건은 수축 사상정리(contraction mapping theorem)에 의해 다음과 같 이 주어진다.

(4)

‖Δε‖<12‖C‖‖∂2ϕ∂σ2‖

여기서 ‖∂2ϕ∂σ2‖ 를 ‖H‖ 로 표기하고, 축차 과정 중의 항복함 수 ϕ의 2계 도함수 행렬(또는 Hessian 행렬)의 크기(norm) 를 나타낸다. 위 식은 주어진 변형률 증분의 크기 ∥∆ε∥ 가 우변에서 정의된 값보다 작아야만 고정점 축차의 수렴이 보 장될 수 있음을 말하고 있다. H 행렬의 원소는 다음과 같다.

H11=∂2ϕ∂σxx∂σxx, H12=∂2ϕ∂σxx∂σyy, …,

4. 평면 응력 조건에서의 수렴 테스트

4.1. 서브-인크리먼트 방법(sub-increment method)

본 논문에서 가정한 Hil l의 이방성 2차 항복함수에 대하 여 2차 도함수 행렬(또는 Hessian 행렬)의 크기를 구하면 ‖∂2ϕ∂σ2‖=‖∂2ϕ∂σ2‖∞=0.006×10−6 (1/Pa)이다. 여기서 벡터 ‖v‖ 와 행렬 ‖H‖ 의 크기는 다음 식을 이용하여 정의하였다[9].

(5)

‖v‖=‖v‖∞=maxi|vi|‖H‖=‖H‖∞=maxi∑j|Hij|

Fig. 2를 참고하면 ‖∂2ϕ∂σ2‖ 의 값은 모든 평면 응력 조건 에서 같은 값을 가짐을 알 수 있다. 탄성계수 E = 160 GPa, 푸아송 비 v = 0.3, 항복응력 σ_f= 500 MPa이고 변형 경화가 없는 경우를 고려하였다. 고정점 축차가 수렴하기 위한 조건은 식 (4)에 의해 주어지며, 그 값은 0.000208 이 된다. 즉 고정점 축차 계산에 사용되는 변형률 증분의 크 기가 ‖Δε‖ < 0.000208을 만족한다면 고정점 축차가 수렴하 게 된다. 이를 검증하기 위하여 평면 응력 조건을 고려하 였다. 즉, 초기에 응력이 없는 질점에 ε_xx와 ε_yy가 가해지는 조건이며, 평면 응력조건을 구현하기 위하여 ε_zz를 조정해 가며 nested loop 기법을 이용하였다[10]. Fig. 3은 주어진 변형률 증분의 크기가 ‖Δε‖ = 0.0002 인 경우 이고, Fig. 4 는 주어진 변형률 증분의 크기가 ‖Δε‖ = 0.0008인 경우이 다. 그림 상의 최종 변형률 상태에 도달하기 위하여, 전체 변형률을 ‖Δε‖ 으로 나눈 후, 각각의 미소 변형률 증분에 대하여 해를 구하고, 그 해를 누적해 나가는 방식으로 수 치계산을 진행하였다. 이러한 기법을 서브-인크리먼트 방 법으로 칭하기로 한다. Fig. 3의 경우는 변형률 증분이 임 계 변형률 증분 0.000208보다 작으므로 모든 축차에서 수 렴함을 볼 수 있다. 반면에 Fig. 4는 변형률 증분이 임계 변형률 증분을 크게 넘어서므로, 어떠한 하중 조건에서는 축차가 발산함을 보여준다(본 논문에서는 축차의 횟수가 20을 넘어서면 발산하는 것으로 간주하였다). 따라서 식 (4)의 임계 변형률 증분이 고정점 축차의 수렴 여부를 판 단하는 기준으로 활용될 수 있음을 알 수 있다. Fig. 3와 4에 나타낸 고정점 축차의 횟수는 주어진 변형률에 도달 하기 위한 마지막 서브-인크리먼트에서의 축차 횟수이다.

Fig. 2

‖∂2ϕ∂σ2‖ has a constant value in the case of the quadratic yield function without strain hardening. The plane stress condition was imposed. The interior ellipse indicates the yield locus in the plane stress condition.

Fig. 3

The number of fixed-point iterations with ‖Δε‖ = 0.0002.

Fig. 4

The number of fixed-point iterations with ‖Δε‖ = 0.0008.

5. 3차원 유한요소 해석

일반적인 응력조건에서 고정점 축차의 수렴 거동을 검 증하기 위하여 3차원 유한요소 해석을 수해하였다. Fig. 5 에 보여진 것과 같은 길이가 10 mm 인 인장시편을 고려 하였으며, 상용 유한요소해석 프로그램 ABAQUS/Explicit [11]를 이용하여 해석을 수행하였다. 고정점 축차 알고리 즘은 사용자 정의 재료 서브루틴(user material subroutine) 을 이용하여 유한요소 해석에 적용하였다. 해석에 사용된 변형률 증분의 크기 ‖Δε‖ 는 현재 질점의 항복강도에서의 ‖∂2ϕ∂σ2‖ 값을 이용하여 계산하였다. 즉, Hessian 행렬은 응력 상태에 따라 달라지지만, 현재 항복면에서의 Hessian 행렬 의 크기는 하나의 값으로 정의 될 수 있다. 변형경화 후의 Hessian 행렬의 크기는 ‖H‖=σ0/σf‖H0‖ 를 이용하여 초기 상태의 크기를 알면 비례식을 이용하여 구할 수 있다. 여 기서 σ₀ 와 σ_f 는 초기와 현재의 항복강도를 나타내고, ‖H0‖ 와 ‖H‖ 는 초기와 현재의 Hessian 행렬의 크기를 나타 낸다. 따라서 유한요소 계산에서 처음에 한 번만 항복함수 의 2계 도함수를 구하면 그 이후에는 값을 계산할 필요가 없다. 만약 ABAQUS 주 프로그램으로부터 주어진 시간 스텝 동안의 변형률 증분이 임계변형률 증분 ‖Δε‖ 보다 크 다면, 주어진 변형률 증분을 ‖Δε‖ 으로 나누어 해를 누적 하였다(이는 앞 절에서 설명한 서브-인크리먼트와 같은 방 법이다).

Fig. 5

Finite element mesh of tensile specimen.

Fig. 6은 계산에 사용된 서브 인크리먼트의 수를 나타낸 다. ABAQUS의 변형률 증분은 ABAQUS 주 프로그램에 서 결정되어 그 값이 user material subroutine 으로 주어진 다. 만약 그 크기가 고정점 축차가 수렴하기 위한 임계 변 형률 증분의 크기보다 크다면, 서브 인크리먼트가 1개 이 상 필요하다. ABAQUS/Explicit의 시간 증분은 매우 작으 므로 소성변형 하는 영역에서는 서브 인크리먼트의 수가 1 임을 Fig. 6으로부터 확인할 수 있다. Fig. 7은 수렴을 하 기 위해 필요한 고정점 축차의 횟수를 나타낸다. 소성변형 하는 영역에서 2번 이하의 축차 후에 수렴하는 것을 확인 할 수 있다.

Fig. 6

The number of sub-increment in ABAQUS/Explicit computation.

Fig. 7

The number of fixed-point iterations in ABAQUS/Explicit computation.

Table 1은 뉴튼-랩슨 축차를 사용하는 경우와 고정점 축 차를 사용하는 경우의 ABAQUS/Explicit 계산 시간을 비 교한 표이다. 일반적으로 UMAT에 의한 계산이 ABAQUS 내장 재료함수보다 많은 시간이 소요된다는 것을 확인할 수 있으며, 고정점 축차가 축차식이 단순하고 축차 과정 동안 항복함수의 2계 도함수를 계산하지 않으므로, 전체 적인 계산시간이 뉴튼-랩슨 기법보다 상당히 단축됨을 확 인할 수 있다. 계산의 정확도를 검증하기 위하여 고정점 축차 해석에서 얻어진 인장시편의 하중-변위 선도를 ABAQUS 내장 재료함수를 사용한 경우와 비교하였다. Fig. 8에서 볼 수 있듯이 3가지 해석 방법 모두 거의 동일 한 하중-변위 선도를 나타내었다. 시편의 넥킹 부위도 세 가지 해석 모두 동일하였으며, 따라서 고정점 축차에 의한 해석 결과는 신뢰할 수 있다고 판단된다.

Table 1

Comparison between computational times in ABAQUS/ Explicit

Fig. 8

Comparison of force-displacement curves of the tensile specimen.

6. 결 론

본 논문에서는 이방성 금속재료의 소성변형을 해석하기 위한 고정점 축차 방법을 제시하였다. 이방성 소성변형을 기술하기 위하여 Hil l의 2차 항복함수를 사용하였다. 제안 한 고정점 축차 기법이 수렴하기 위한 충분조건을 제시하 였으며, 수치적으로 수렴성을 검증하기 위하여 2차원 평 면응력 조건을 고려하였다. 또한 제시한 방법의 수치적 효 율성을 증명하기 위하여 3차원 인장시편의 유한요소 해석 을 수행하였고, 계산 시간을 뉴튼-랩슨 방법을 이용한 회 귀 착점 방법과 비교하였다. ABAQUS/Explicit를 이용한 3차원 유한요소 해석을 통하여, 제시한 고정점 축차 방법 이 뉴튼-랩슨 방법보다 효율적임을 증명하였다. 단, implicit FEM과 같이 시간 증분이 매우 큰 경우에는, 고정 점 축차는 많은 수의 서브 인크리먼트가 필요하다. 이 경 우 고정점 축차에 의한 CPU 시간은 뉴튼-랩슨 기법에 의 한 CUP 시간 보다 증가할 것으로 예상된다.

Acknowledgements

감사의 글

이 논문은 한국기술교육대학교 교수 교육연구진흥과제 지원에 의하여 연구되었음.

• 1. R. Wauthle, B. Vrancken, B. Beynaerts, K. Jorissen, J. Schrooten, J.-P. Kruth and J. V. Humbeeck: Addit. Manuf., 5 (2015) 77.Article
• 2. L.-C. Shen, X.-H. Yang, J.-R. Ho, P.-C. Tung and C.-K. Lin: Materials, 13 (2020) 5142.ArticlePubMedPMC
• 3. S.-Y. Yang and W. Tong: NUMISHEET 2022, Proceedings of the 12^th International Conference and Workshop on Numerical Simulation of 3D Sheet Metal Forming Processes, (2022) 3.Article
• 4. W. M. Scherzinger: Comput. Methods Appl. Mech. Eng., 317 (2017) 526.Article
• 5. H. Akpama, M. B. Bettaieb and F. A. Abed-Meraim: Int. J. Numer. Methods Eng., 108 (2016) 363.Article
• 6. R. Hill: The Mathematical Theory Of Plasticity, Oxford Univ. Press, (1998).
• 7. S. Y. Yang, D. H. Jin and J. H. Kim: J. Powder Mater., 29 (2022) 303.Article
• 8. S. Y. Yang, J. H. Kim and D. H. Jin: The 17^th International Symposium on Novel and Nano Materials Abstract Book, (2022) 279.
• 9. B. Noble and J. W. Daniel: Applied Linear Algebra, Prentice- Hall Inc., (1977).
• 10. E. A. de Souza Neto, D. Peric and D. R. J. Owen: Computational Methods for Plasticity: Theory and Applications, Wiley (2008).ArticlePubMed
• 11. ABAQUS 6.13 Documentation-Dassault systemes.

Figure & Data

References

Citations

Citations to this article as recorded by