Journal Search Engine
Search Advanced Search Adode Reader(link)
Download PDF Export Citaion korean bibliography PMC previewer
ISSN : 1225-7591(Print)
ISSN : 2287-8173(Online)
Journal of Korean Powder Metallurgy Institute Vol.27 No.4 pp.318-324
DOI : https://doi.org/10.4150/KPMI.2020.27.4.318

Computation of Stress Field During Additive Manufacturing by Explicit Finite Element Method

Seung-Yong Yanga*, Jeoung Han Kimb,*
aKorea University of Technology and Education, School of Mechanical Engineering, Cheonan, 31253, Republic of Korea
bHanbat National University, Department of Materials Science and Engineering, Daejeon, 34158, Republic of Korea
-

양승용·김정한: 교수


*Corresponding Author: Seung-Yong Yang, TEL: +82-41-560-1459, FAX: +82-41-560-1459, E-mail: ysy@kut.ac.kr
Jeoung-Han Kim, TEL: +82-42-821-1240, FAX: +82-42-821-1592, E-mail: jh.kim@hanbat.ac.kr
August 7, 2020 August 24, 2020 August 26, 2020

Abstract


In the present work, an explicit finite element analysis technique is introduced to analyze the thermal stress fields present in the additive manufacturing process. To this purpose, a finite element matrix formulation is derived from the equations of motion and continuity. The developed code, NET3D, is then applied to various sample problems including thermal stress development. The application of heat to an inclusion from an external source establishes an initial temperature from which heat flows to the surrounding body in the sample problems. The development of thermal stress due to the mismatch between the thermal strains is analyzed. As mass scaling can be used to shorten the computation time of explicit analysis, a mass scaling of 108 is employed here, which yields almost identical results to the quasi-static results.



외연적 유한요소법을 이용한 적층제조 공정 중 응력 장 변화 계산

양 승용a*, 김 정한b,*
a한국기술교육대학교 기계공학부
b한밭대학교 신소재공학과

초록


    Korea University of Technology and Education

    1. 서 론

    3차원 적층제조(Additive Manufacturing) 기법은 통상적 으로 레이저 또는 전자빔 등의 열원을 이동시켜가며 금속 등의 분말을 용융온도까지 가열하고 응고시켜 원하는 형 태를 만드는 과정으로 구성되어 있다[1]. 급속한 가열과 냉각으로 인하여 금속재료는 팽창 및 수축하는 과정에서 잔류응력을 발생시키게 되고 재료에 이방성 등의 미세조 직을 초래할 수 있다[2]. 잔류응력은 구조물의 피로수명 및 파괴거동에 영향을 미칠 수 있으므로 제조 과정에서 잔류응력을 조절할 수 있는 방법이 요구된다. 적층과정 에 서의 온도 분포와 잔류응력 분포를 시험을 통하여 측정하 는 것은 복잡한 과정을 거치고 비용이 많이 들게 되므로 유한요소법 등을 이용하여 계산하는 방법이 많이 사용되 고 있다[3].

    이러한 열응력을 해석하기 위하여 유한요소법, 유한체 적법 등의 수치해석 기법이 사용되고 있다. Mukherjee 등 [4, 5]은 용융 금속내의 대류 운동을 고려하기 위하여 유 체역학에서 사용되는 유한체적법을 사용하였다. Abott 등 [6]은 3차원 적층제조 공법으로 제조된 시편에 대하여 유 한요소 해석을 수행하였고 결과를 실험결과와 비교하였다.

    본 논문에서는 외연적(Explicit) 유한요소법을 이용하여 금속 재료의 가열 및 냉각과정에서의 응력 변화를 해석하 고자 한다. 외연적 유한요소법(explicit method)은 구조물 의 충돌 및 파괴와 같이 짧은 시간에 급격한 변화가 발생 하는 문제에 주로 사용되는 해석 기법으로서 구조물의 운 동방정식을 시간에 대하여 풀어서 물체내의 변위와 응력 을 계산하는 기법이다. 질량 행렬을 대각화 시키므로 질량 행렬의 역행렬은 행렬 연산 없이 바로 대각 원소의 역수 로 구할 수 있고, 시간 증분이 작으므로 변위를 시간에 대 하여 근사적으로 적분함으로써, 뉴튼 반복(Newton iteration) 계산을 하지 않더라도 쉽게 변위를 구할 수 있으므로 해 석 코드를 단순화 할 수 있는 장정이 있다. 본문에서는 평 형 및 연속 방정식을 기술한 후 유한요소 정식화 과정을 기술하였고 간단한 예제에 적용한 결과를 제시하였다.

    2. 열팽창을 고려한 유한요소 방정식

    물질의 푸아송 비가 0.5 이거나 전위에 의해 소성변형이 발생하는 경우 변형이 비압축성을 나타내게 된다. 이러한 비압축성 변형의 경우를 일반적인 유한요소법으로 해석하 면 체적 잠김(volumetric locking) 현상이 발생하고 압력의 계산값이 불안정해 지는 문제가 있다. 이를 해결하기 위하 여 Zienkiewicz는 연산자분리(operator split) 법을 개발하 여 운동방정식의 적분에 이용하였다[7]. 연산자 분리 기법 에서는 운동량 보존 방정식을 다음과 같이 변형한다. 밀도 변화가 크지 않고 체적력을 무시할 수 있는 경우 운동량 보존 방정식은 다음과 같다[7].

    ρ 0 v t D i v τ p = ρ 0 v * t p
    (1)

    여기서 ρ0 는 물체의 밀도이고 τ 는 편차응력 이고 p는 압 력이다. v*Div τ 에 해당하는 중간 속도이다. 편차응력 에 의한 속도 v* 를 분리해서 구한 후, 연속 방정식에 대 입하여 압력의 시간적분에 이용한다. 연속 방정식은 다음과 같다. 연속 방정식과 압력에 대한 구성방정식을 결합하면

    ρ 0 K p t = ρ t = D i v ( ρ v ) ρ 0 D i v v

    K는 체적 탄성계수이다.

    온도 변화에 의한 체적 변형률 ε ˙ v t 이 존재하는 열응력 문제라면, 압력 변화와 체적 변화 사이의 관계식은 다음과 같다.

    1 c 2 p t = ρ 0 ( ε ˙ i i ε ˙ v t ) = ρ 0 D i v v + ρ 0 ε ˙ v t
    (2)

    (1)과 (2)는 속도 v와 압력 p을 구하기 위한 미분방정식 이다. 각 절점에서의 속도 v와 압력 P을 이용하여 유한요 소 내의 임의의 점에서의 속도와 압력 근사값을 표현하면

    v = N v V p = N p P

    형상함수 N은 다음과 같은 벡터이다.

    N v ( x ) = N p ( x ) = ( N 1 ( x ) , N 2 ( x ) , , N M + 1 ( x ) )

    갈레아킨 방법[8]을 적용하기 위하여 운동방정식에 형 상함수를 곱하고 적분한 값을 고려하면 다음과 같은 식을 얻을 수 있다.

    N v T ρ 0 v t d V = B v T τ d V + ( N v T T ) m p d V + N v T t d A
    (3)

    = N v T ρ 0 v * t d V + B v T m N p d V P N v T p n d A
    (4)

    여기서 B v = N v 이고 m = (1, 1, 1, 0, 0, 0)T이다[9]. 운동방정식 (4)식을 행렬식으로 나타내면,

    M Δ V ¯ Δ t = M Δ V ¯ * Δ t + Q P n + f p
    (5)

    시간 단계 tn 에서의 압력과 (3)식으로부터 구한 ΔV를 이용하여 (5)식에서 ΔV*를 구하게 된다. 구해진 ΔV* 값 은 다음 시간에서의 압력 Pn+1을 구하는 데 이용한다. 즉 연속방정식 (2)식을 유한요소 행렬식으로 나타내면

    N p T 1 ρ 0 c 2 p t d V = N p T D i v ( v + Δ v * Δ t ρ 0 p ) d V N p T ε ˙ v t d V M ˜ P ¯ n + 1 P ¯ n Δ t = Q T V ¯ n Q ¯ Δ V ¯ * + Δ t H P ¯ n + f v E ˙ t
    (6)

    여기서

    Q ( 4 × 12 ) T = N p T m T N v d V Q ¯ ( 4 × 12 ) = ( ¯ N p ) T N v d V f v = N p T n T ( Δ v * Δ t ρ 0 Δ p ) d A = N p T n T Δ v d A E ˙ t = N p T ε ˙ v t d V H ( 4 × 4 ) = ( ¯ N p ) T 1 ρ 0 ¯ N p d V

    선형 4면체 요소에서는 열 팽창 변형률 ε ˙ v t 가 각 유한요 소 체적 내에서 상수이다.

    (6)식은 압력의 구배가 있으므로 확산 방정식의 형태를 띄고 있으며, 이 식을 풀면 Vn 과 ΔV*Vn을 알 때 Pn+1 을 구할 수 있다. 외연적 시간 적분법에서는 시간 증분이 다음의 CFL 안정 조건을 만족하여야 한다.

    K ρ Δ t Δ x < 1

    보통 이 시간 증분은 매우 작은 값이므로, 해석 시간이 긴 문제에서는 밀도를 실제 값보다 큰 값을 사용하여 시 간 증분을 증가시키게 된다. 이를 질량 스케일링(mass scaling)이라고 부른다. 질량 스케일링을 적당히 사용하면 계산 시간을 줄이면서도 비교적 정확한 값을 얻을 수 있다.

    3. 구성방정식

    열응력 문제에서 (3) 식에 의하여 ΔV 를 구하기 위해서 는 현재 시간에서의 편차응력 τ를 알아야 한다. 물체의 온 도가 변화하는 경우 열변형률이 발생하므로, 주어진 변형 률 변화량과 온도 변화량에 대하여 응력 증가량은 다음 식으로 표현된다.

    σ ˙ = C ( ε ˙ ε ˙ p ε ˙ t )
    (7)

    ε ˙ p 는 소성변형률 속도이고, ε ˙ t 는 열변형률 속도이다. 소 성변형률 속도는 등방성 소성변형을 가정하여 Prandtl-Reuss 방정식으로부터 구하게 되고, 열변형률은 등방성 열팽창 률 α(T) 를 이용하여 다음과 같이 나타낼 수 있다 [10].

    ε t = α ( T ) ( T T r e f ) δ α ( T ) ( T 1 T r e f ) δ
    (8)

    여기서 δ 는 단위 행렬이다.

    Tref는 기준온도이고 T1은 물체의 초기 온도이다. 위 식 은 물체가 초기 온도에서 열변형률이 없음을 나타낸 식이 다. 온도가 변함에 따라 열변형률이 변화하므로 응력의 변 화가 발생하게 된다. 탄성계수 C도 온도의 함수이고 따라 서 시간의 함수일 수 있으나, C ˙ 의 영향은 계산 결과에서 크지 않으므로 무시하였다.

    (1)식 또는 (3)식에서 편차응력이 필요하므로, 편차응력 τ ˙ 은 (7)식에서 편차 변형률 만을 고려하여 구하게 된다.

    τ ˙ = 2 μ ( ε ˙ t ε ˙ p ε ˙ t )

    (9) 식에는 전단 탄성계수 μ 만 들어가고, 비압축성 재 료의 경우도 μ 가 무한대가 되는 경우가 발생하지 않는다.

    4. 열응력 수치해석 결과

    앞 절에서 기술한 외연적 유한요소 행렬식을 이용하여 해석 코드 NET3D 를 개발하였으며 코드의 검증을 위하 여 간단한 해석 모델에 대하여 열응력 장을 계산하였다. 적층 제조공정의 기초적인 해석을 위하여 용접과정[11]과 의 유사성에 착안하여 간단한 열응력 문제를 고려하고자 한다. 먼저 그림 1과 같은 10 × 10 × 10 mm 육면체에 대하 여 바닥면을 고정시키고 물체의 온도를 200초 동안 1500 K에서 100 K로 선형적으로 감소시켰다. 바닥면이 고 정되어 있고 물체가 온도 감소에 의해 수축하려 하므로 물체에 비틀림이 발생하고 수축 과정에서 소성 변형이 발 생하게 된다. 물체의 재료 상수 값은 Table 1에 나타나 있 는 Ti-6Al-4V 의 값을 이용하였다.

    재료 상수 값은 ABAQUS 해석 시에는 테이블로 주어진 참고문헌의 값을 직접 사용하였고 저자의 해석 코드에서 는 온도의 2차 함수로 근사한 함수 식을 사용하였다. 그림 2에 시간에 따른 물체의 평균적인 von Mises 응력을 나타 내었다. 또한 시간에 따른 항복응력의 변화를 같이 나타내 었다. 시간이 지남에 따라 물체의 온도가 감소하므로 물체 가 경화되어 항복응력은 증가하게 되고 열변형에 따라 von Mises 응력이 증가하다가 항복응력과 같아지는 시간 에 소성변형이 발생하게 됨을 그림 2가 보여주고 있다. 해 석 코드의 결과는 ABAQUS의 결과와 같은 경향을 보이 며 약간의 차이를 보이는 이유는 저자의 해석 코드는 물 성 데이터로써 2차함수 근사식을 사용하였기 때문이라고 판단된다.

    다음 예제로 그림 3과 같은 400 × 400 × 100 mm 육면체 에 대하여 초기 조건으로 중앙 실린더 블록은 1500 K로 지정하고 주변 블록은 300 K로 지정하여 시간에 따라 자 연적으로 열전도가 발생하도록 한 경우를 해석하였다. 테 두리의 x면과 y면은 모든 방향의 변위가 발생하지 않도록 구속하였다. 물성값은 Table 1의 값을 사용하였다. 시간이 경과하면서 중앙 코어 블록은 냉각되므로 수축하게 되고 주변부는 가열되므로 팽창하게 되어 그림 4와 같은 열팽 창 변형률의 분포를 보이게 된다. 열에 의한 비틀림에 의 하여 냉각과정에서 소성변형이 발생하게 되고 그림 5와 같은 등가 소성변형률의 분포를 보이게 된다. 그림 6은 단 면에서의 압력 분포를 ABAQUS 결과와 비교한 것이다. 두 결과가 동일한 경향을 나타내고 있으나, 본 논문의 해 석 코드가 조금 더 완만한 압력 변화를 나타내고 있다. 이 는 압력을 확산방정식을 풀어서 구하게 되므로 변화가 완 만하게 구해지기 때문이라고 판단된다. 질량 스케일링은 108을 사용하였다. 외연적 ABAQUS Explicit와 Standard 를 이용하여 질량 스케일링을 108 사용할 경우와 사용하 지 않을 경우의 준정적(quasi static) 해석결과가 거의 동일 함을 확인하였고, 따라서 질량 스케일링 108은 해석 결과 를 크게 왜곡하지 않음을 알 수 있었다. 외연적(explicit) 유한요소 코드를 사용하여 열응력을 해석할 경우, 질량 스 케일링을 통해 계산 시간을 단축 할 수 있음을 알 수 있 다. Figure 78은 저자의 해석코드와 ABAQUS에서 계 산된 von Mises 응력을 비교한 것이다. 두 코드의 해석결 과가 거의 동일함을 확인할 수 있었다.

    다음 예로서, 적층 제조 공정에서 나타나는 표면 가열에 의한 열응력 발생 문제를 고려하였다. 문제의 구조물은 그 림 3과 같은 형태이고 중앙 실린더 블록의 윗면에 레이저 가 조사되어 가열된다고 가정하였다. 표면에 투입되는 복 사 에너지의 양은 20초 동안 5 × 106 J/sm2이다. 투입된 열 은 자연 전도에 의해 주변으로 전달된다. 물체의 둘레는 전과 같이 모든 방향으로 고정된 것으로 가정하였다. 물성 값은 Table 2의 C r-Mo 값을 사용하였다. 그림 910은 저자의 프로그램과 ABAQUS에서 계산된 등가 소성 변형 률의 분포를 단면에 대하여 비교한 것이다. 그리고 그림 11은 압력의 분포를 대각선 경로를 따라 비교한 것이다. 두 결과가 거의 비슷한 값을 나타낸다. 그림 6에서와 같이 압력이 약간의 부정확 성을 보이는 것은 (6)식에 시간 증 분 Δt가 있으며 열전도 해석에서 일반적으로 시간 증분이 크기 때문이 것으로 생각된다. 연산자 분리 방법에서 계산 오차는 시간 증분에 비례함을 Aref와 Balachandar[12]는 제시하였다.

    5. 결 론

    적층제조 공정에서 재료는 급속한 가열 및 냉각 과정을 거치게 되고, 가열 및 냉각 과정 모두에서 주변부와의 온 도 차이로 인하여 뒤틀림을 받게 되고, 항복현상이 발생하 게 된다. 이러한 역학적 과정을 해석하기 위하여 외연적 (explicit) 유한요소법을 적용하였으며, 간단한 모델에 대하 여 적용한 결과, 상용 해석 코드인 ABAQUS와 동일한 경 향을 보임을 확인할 수 있었다. 질량 스케일링(mass scaling) 을 통하여 시간 증분을 증가시켰으며 계산시간을 단축할 수 있었다. 질량 스케일링 한 결과는 준 정적 해석 결과와 거의 동일함을 확인하였다.

    앞으로 다양한 구속조건에서 과정을 해석함으로써, 적 층제조 공정에서 발생하는 열응력에 의한 미세조직 발생 연구 등에 활용될 수 있을 것으로 판단된다.

    감사의 글

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

    Figure

    KPMI-27-4-318_F1.gif
    Finite element model of asymmetric contraction. Bottom surface was pinned and temperature changed from 1500 K to 100K.
    KPMI-27-4-318_F2.gif
    Average von Mises stress and yielding stress. Dotted line is yielding stress. Red line is NET3D result and blue line is ABAQUS results.
    KPMI-27-4-318_F3.gif
    Initial temperature distribution in Kelvin unit. The core block was assigned 1500 K initially.
    KPMI-27-4-318_F4.gif
    Volumetric thermal strain at time = 50 s along a diagonal path.
    KPMI-27-4-318_F5.gif
    Equivalent plastic strain distribution in the cross section of the block at time = 50 s.
    KPMI-27-4-318_F6.gif
    Pressure distribution at time = 50 along a diagonal path. Red line is ABAQUS and green line is NET3D.
    KPMI-27-4-318_F7.gif
    Computed von Mises stress distribution by NET3D at time = 50 s.
    KPMI-27-4-318_F8.gif
    Von Mises stress distribution by ABAQUS at time = 50s.
    KPMI-27-4-318_F9.gif
    Equivalent plastic strain by NET3D at time = 110 s. The cylindrical top surface was heated by surface heating.
    KPMI-27-4-318_F10.gif
    Equivalent plastic strain by ABAQUS at time = 110 s.
    KPMI-27-4-318_F11.gif
    Pressure distribution at time = 110 s along a diagonal path.

    Table

    Ti-6Al-4V material values [4]
    Cr-Mo material values [5]

    Reference

    1. N. K. Adomako, S. Noh, C. Oh and S. Yang: Mater. Res. Lett., 7 (2019) 259.
    2. T. DebRoy, H. L. Wei, J. S. Zuback, T. Mukherjee, J. W. Elmer, J. O. Milewski, A. M. Beese, A. Wilson-Heid, A. De and W. Zhang: Prog. Mater. Sci., 92 (2018) 112.
    3. M. M. Francois, A. Sun, W. E. King, N. J. Henson, D. Tourret, C. A. Bronkhorst, N. N. Carlson, C. K. Newman, T. Haut, J. Bakosi, J. W. Gibbs, V. Livescu, S. A. Vander Wiel, A. J. Clarke, M. W. Schraad, T. Blacker, H. Lim, T. Rodgers, S. Owen, F. Abdeljawad, J. Madison, A. T. Anderson, J. L. Fattebert, R. M. Ferencz, N. E. Hodge, S. A. Khairallah and O. Walton: Curr. Opin. Solid State Mater. Sci., 21 (2017) 198.
    4. T. Mukherjee, W. Zhang and T. DebRoy: Comp. Mater.Sci., 126 (2017) 360.
    5. T. Mukherjee, J. S. Zuback, W. Zhang and T. DebRoy: Comp. Mater. Sci., 143 (2018) 325.
    6. D. W. Abbot, D. V. V. Kallon, C. Anghel and P. Dube: Proce. Manuf., 35 (2019) 164.
    7. O. C. Zienkiewicz, J. Rojek, R. L. Taylor and M. Pastor: Int. J. Num. Methods Engrg., 43 (1998) 565.
    8. O. C. Zienkiewicz and K. Morgan: Finite elements and approximations, Dover Publications, (1983).
    9. O. C. Zienkiewcz, R. L. Taylor and P. Nithiarasu: The finite element method for fluid dynamics, 7th ed. Butterworth-Heinemann, (2014).
    10. ABAQUS 6.14 User’s Guide, Dassault Systemes.
    11. S. Y. Yang, B. C. Koo and H. C. Jeong: J. Korean Soc. Railway, 7 (2004) 121.
    12. H. Aref and S. Balachandar: A first course in computational fluid dynamics, Cambridge Univ. Press, (2018).