8 변형 문제의 유한요소 수식화

8.1 개론 및 경계치 문제의 정의

그림 8.1은 일반적인 소성가공 시뮬레이션의 흐름도 및 개념을 나타낸 것이다. 그림 8.1에서 보는 바와 같이 소성가공 시뮬레이션은 전처리 프로그램(preprocessor), 솔버(solver), 후처리 프로그램(post-processor)으로 구성되어 있다.

그림 8.1 소성가공 시뮬레이션

처리 프로그램은 기본적으로 CAD 시스템과 연결되며, 재료, 프레스, 마찰 등의 데이타베이스(DB)의 지원을 받는다. 솔버는 속도, 압력, 온도 등의 미지수를 사용자가 요구한 공정 조건에 맞도록 구해주는 역할을 한다. 후처리 프로그램은 전술한 미지수와 이로부터 파생된 변형률속도, 변형, 응력, 성형하중 등의 정보를 사용자가 용이하게 정보화하도록 돕는다. 전후처리기는 사용의 편리성으로 평가되며, 솔버는 결과의 정확도, 계산 속도, 적용의 다양성으로 평가된다.

재료, 프레스, 마찰 등의 정보는 기본적으로 사용자의 영역이다. 물론 소프트웨어가 기본적인 정보를 제공하지만, 궁극적으로는 사용자의 환경에 맞도록 지속적으로 DB를 확충해 가는 자세가 필요하다. 특히 재료의 유동응력은 결과를 좌우하는 것으로 국가마다 재료의 제조사마다 그 물성치, 즉 재료의 거시적 거동특성이 다르고, 열처리 및 인발 등의 전처리 및 단조 조건에 따라 그 성질이 비교적 크게 변한다는 점을 인식할 필요가 있다.

소성가공 시뮬레이션은 자연법칙으로부터 편미분방정식으로 유도된 평형방정식과 열전도방정식을 유한요소기교를 이용하여 이산화시켜서 비선형방정식을 만들고 이를 풂으로써 일차적으로 속도, 압력(라그랑지변수법), 온도변화율 등을 구한 후 이로부터 변형률속도, 변형형상, 온도 분포 등의 각종 파생변수를 구하는 것으로 요약된다.

해석 정보는 크게 공정기하 정보, 공정조건 정보, 재료 물성치 정보 등으로 나누어진다. 공정기하 정보는 공정도, 금형의 형상, 소재의 형상 등을 포함한다. 공정조건 정보 또는 공정조건은 마찰, 금형속도, 열전달계수 등의 데이타를 포함한다. 그리고 재료 물성치 정보에는 재료의 유동응력, 열전도계수, 열용량 등을 포함한다. 그 이외에도 해석 기술자가 결정하거나 정의해야 할 것으로 대칭면 정보, 행정(stroke)정보, 요소의 수, 방정식의 해법 등이 있다. 이 장은 참고문헌 [1.1]의 내용을 요약하고 필요에 따라 새로운 내용을 추가한 것이다.

8.2 강소성 유한요소법

8.2.1 강소성변형 문제의 정의

소성가공 재료는 비압축성, 강소성, 등방성 경화(isotropic hardening) 등의 성질을 가지며 Huber-von Mises 항복이론과 직교 유동법칙(associated flow rule)을 따르고, 접촉면에서의 쿨롱마찰법칙, 일정전단마찰법칙, 하이브리드마찰법칙 등을 따른다고 가정한다. 그리고 자중과 가속도의 영향은 무시할 수 있을 정도로 작다고 가정한다.

그림 8.2에서 보는 바와 같이 소성가공 공정에서 금형과 소재는 해석영역 \(V\) 와 경계 \(S\) 로 나누어지며, 경계 \(S\) 는 표면력(traction)이 \(t_i^{(n)} = \bar{t}_i^{(n)}\) 로 주어진 \(S_{t_i}\), 속도가 \(v_i = \bar{v}_i\) 로 주어진 \(S_{v_i}\), 금형과 소재의 접촉면을 표현하는 \(S_c\) 등으로 세분화될 수 있다. 압축력이 작용하는 금형과 소재의 접촉에서의 법선속도는 동일하다. 즉, \(S_c\) 에서 \(v_n = \bar{v}_n\) 이다.

소성가공 공정의 유동해석 문제는 소재의 속도장 \(v_i\) 와 정수압 \(p\) (라그랑지변수법) 또는 속도장 \(v_i\) (벌칙기법)를 구하고자 하는 것이며, 다음과 같은 경계치 문제로 수식화된다.

그림 8.2 소성역학 문제의 개념도

(1) 평형방정식

\[ \sigma_{ij,j} + f_i = 0 \tag{8.1} \]

(2) 응력-변형률속도 관계식

\[ \sigma_{ij} = -p \delta_{ij} + \sigma'_{ij} \tag{8.2} \]

\[ \sigma'_{ij} = \frac{2\bar{\sigma}}{3\dot{\bar{\varepsilon}}} \dot{\varepsilon}_{ij} \tag{8.3} \]

(3) 속도-변형률속도 관계식

\[\dot{\varepsilon}_{ij} = \frac{1}{2}(v_{i,j} + v_{j,i}) \tag{8.4}\]

(4) 비압축성조건

\[v_{i,i} = 0 \tag{8.5}\]

(5) 경계조건

\[\sigma_{ij}n_j = \bar{t}_i^{(n)} \quad \text{on} \quad S_{t_i} \tag{8.6}\]

\[v_i = \bar{v}_i \quad \text{on} \quad S_{v_i} \tag{8.7}\]

\[v_n = \bar{v}_n \quad \text{on} \quad S_c \tag{8.8}\]

\[\sigma_t = -\mu\sigma_n g(v_t) \quad \text{on} \quad S_c \quad \text{(Coulomb friction)} \tag{8.9}\]

\[\sigma_t = mkg(v_t) \quad \text{on} \quad S_c \quad \text{(constant shear friction)} \tag{8.10}\]

\[\sigma = \min(-\mu\sigma_n, m'k)g(v_t) \quad \text{on} \quad S_c \quad \text{(hybrid friction)} \tag{8.11}\]

여기서 \(\min(*, *)\) 는 두 값 중에서 작은 값을 취함을 의미하며, \(m'\) 은 마찰응력의 최대치를 구속하기 위한 상수, 즉 한계마찰상수이다. 평형방정식에서 \(f_i\) 는 체적력이다. 구성방정식에서 \(\bar{\sigma}\) 와 \(\dot{\bar{\varepsilon}}\) 는 각각 유동응력과 유효변형률속도를 의미한다. 일반적으로 식 (8.3)에서 유동응력

\[ \bar{\sigma} = \bar{\sigma}(\bar{\varepsilon}, \dot{\bar{\varepsilon}}, T) \tag{8.12} \]

는 변형률 \(\bar{\varepsilon}\), 변형률속도 \(\dot{\bar{\varepsilon}}\), 온도 \(T\) 등의 함수이며, 재료의 기계적 성질 또는 상태를 반영하는 함수이다. 식 (8.6)의 벡터 \(n_i\) 는 해석경계에서 정의되며 외향단위법선벡터(outward unit normal vector)이다. 텐서량에서 상첨자 프라임(\('\))은 편차 성분(deviatoric component)을 나타내며, 하첨자 \(t\) 와 \(n\) 은 각각 접선과 법선 성분을 의미한다. 그리고 하첨자가 두 번 반복되는 것은 덧셈에 대한 부호규약을 따른 것이며, 콤마(,) 뒤의 첨자는 편미분을 의미한다. 이러한 부호규약은 모든 수식에 적용된다. 식 (8.9)와 (8.10)에서 \(g(v_t)\)는 다음과 같다.

\[g(v_t) = -\frac{2}{\pi} \tan^{-1} \left( \frac{v_t - \bar{v}_t}{a} \right) \tag{8.13}\]

여기서 \(a\)는 \(|v_t|\)에 비하여 매우 작은 양의 상수이다.

비압축조건을 처리하기 위하여 벌칙기법(penalty method)을 사용할 경우와 라그랑지변수법(Lagrange multiplier method)을 사용할 경우, 위의 경계치 문제의 약형(weak form)은 각각 다음과 같다[1.1].

\[\int_V \sigma'_{ij} \omega'_{ij} dV + \int_V K \dot{\varepsilon}_{ii} \omega_{jj} dV - \int_V f_i \omega_i dV - \sum \int_{S_{t_i}} \bar{t}_i^{(n)} \omega_i dS - \int_{S_c} \sigma_t \omega_t dS = 0 \tag{8.14}\]

\[\int_V \sigma'_{ij} \omega'_{ij} dV - \int_V p \omega_{jj} dV - \int_V f_i \omega_i dV - \int_V v_{i,i} q dV - \sum \int_{S_{t_i}} \bar{t}_i^{(n)} \omega_i dS - \int_{S_c} \sigma_t \omega_t dS = 0 \tag{8.15}\]

여기서 가중함수 \(\omega_i\)는 경계 \(S_{v_i}\)에서 영의 값을 갖는 임의의 함수이고, 가중함수 \(q\)는 전 해석영역에서 임의의 함수이다. 약형에서

\[\omega_{ij} = \left( \frac{\omega_{i,j} + \omega_{j,i}}{2} \right) \tag{8.16}\]

이며, \(\omega_n\)과 \(\omega_t\)는 각각 \(\omega_i\)의 법선성분과 접선성분을 나타낸다. 벌칙기법에서, 벌칙상수 \(K\)는 다음과 같은 물리적 의미를 갖는다.

\[K \dot{\varepsilon}_{ii} = \sigma_{jj} / 3 = -p \tag{8.17}\]

따라서 벌칙기법에서 압력분포는 속도장 해석결과로부터 구해진다. 결론적으로 전술한 경계치 문제를 푸는 것은 라그랑지변수법을 사용할 때는 임의의 가중함수 \(\omega_i\) 와 \(q\) 에 대하여 약형, 즉 식 (8.15)를 항상 만족시키는 속도장 \(v_i\) 와 압력 \(p\) 를 구하는 문제이며, 벌칙기법을 사용할 때는 임의의 가중함수 \(\omega_i\) 에 대하여 약형, 즉 식 (8.14)를 항상 만족시키는 속도장 \(v_i\) 를 구하는 문제로 수식화된다. 물론, 속도장은 \(S_{v_i}\) 와 \(S_c\) 에서 필수경계조건을 만족해야 한다.

라그랑지변수법에서 약형에는 속도장 \(v_i\) 와 압력 \(p\), 가중함수 \(\omega_i\) 와 \(q\) 등이 바탕을 형성하고 있다. 그 이외의 것은 이들로부터 파생된 것이다. 따라서 이 네 개의 함수를 보간하여 절점의 속도와 요소의 압력이 미지수인 비선형연립방정식을 얻을 수 있다. 벌칙기법에서는 속도장 \(v_i\) 와 가중함수 \(\omega_i\) 를 보간하여 절점의 속도가 미지수인 비선형연립방정식을 얻을 수 있다. 유한요소법에 관한 기초지식을 부록 B에 수록하였다. 이 책에서는 강소성 유한요소법의 유한요소방정식의 유도과정은 참고문헌[1.1]로 대신하는 반면, 탄소성 유한요소법에 관한 그 과정을 제8.3절에서 상술한다.

8.2.2 최소유효변형률속도

식 (8.3)의 분모에는 유효변형률속도 \(\dot{\bar{\varepsilon}}\) 가 존재한다. 이 값은 이론적으로는 탄성영역에서 0.0이며, 이 때 수치적으로 계산이 불가능한 문제에 빠지게 된다. 이러한 문제를 해결하기 위하여 유효변형률속도가 허용최소유효변형률속도 이하가 되면 그 값 대신 이 허용최소유효변형률속도를 사용한다. 이 수치는 해석결과의 정확도와 수렴특성에 큰 영향을 미칠 수 있다. 만약 이 값을 크게 잡으면, 해는 비교적 쉽게 구해지겠지만, 특성영역의 과도한 수치적 변형을 야기시키므로 대부분 해의 정확도를 떨어뜨리는 부정적 역할을 한다. 반대로 이 값을 매우 작게 잡으면, 해는 보다 이론해에 가깝게 될 가능성은 높지만 해의 수렴특성이 나빠지게 된다. 따라서 이 값을 결정하는 데도 비교적 세심한 주의가 필요하다.

8.2.3 금형침투허용 및 업데이팅

일반적으로 소성가공 시뮬레이션 중 속도장 정보를 바탕으로 절점의 위치 결정 시에 절점의 금형 침범을 허용하고 있다. 절점이 금형을 침범하지 못하게 하면, 해석의 전진이 늦어지는 문제가 있다. 물론 침투한 절점의 위치는 인위적으로 가까운 금형 상의 점으로 수정되거나 아니면 침투한 상태를 유지하되 금형표면으로 밀어낸다. 전자는 인위적 부피변화 측면에서 유리하지 않다.

절점의 금형침투허용깊이는 예측결과의 정확도와 계산시간 등에 무시할 수 없는 정도의 영향을 주고 있다. 경험에 의하면, 원하는 예측결과의 정확도의 50% 이하로 절점의 금형침투허용깊이를 정하는 것이 바람직하다. 가령 0.02mm의 정확도로 해석하고자 할 때는 0.01mm 이하를 허용하는 것이 바람직하다. 침투절점 밀어내기 기법을 사용한다면, 이 값을 3배까지 증가시켜도 무방하다. 대개의 프로그램들은 사용자의 큰 실수를 방지하기 위하여 기본값(default value)을 정해 두고 있다. 소성가공 시뮬레이터로 0.001mm의 정확도를 원하는 경우는 일반적으로 흔하지 않기 때문에 세심한 주의를 기울이지 않으면, 마이크로 포밍공정(부록 D.2.1 참조)과 같이 치수 자체가 매우 작은 공정의 해석 시에 문제가 될 수가 있다. 프로그램의 지능화가 많이 되어 있을수록 이런 문제는 더욱 심하게 된다.

8.3 탄소성 유한요소법

강소성 유한요소해석은 금속재료가 소성가공 중에 소성변형만 발생한다고 전제한다. 즉, 탄성변형이 무시된다. 전체적인 재료의 유동이나 하중 등의 관찰에는 강소성 유한요소해석으로 충분한 경우가 많고 짧은 계산 시간 및 수렴성의 장점이 있어 많은 경우에 유용하게 활용되고 있다. 소성가공 과정에서 재료는 소성변형과 동시에 탄성변형을 겪는 반면 강소성 유한요소해석은 탄성변형을 고려하지 않기 때문에 응력 및 변형의 계산에서 어느 정도의 오차가 불가피하며 하중이 제거될 때 탄성회복을 계산하지 못하게 된다. 강소성 유한요소해석에 비해 탄소성 유한요소해석은 이러한 문제점을 해결할 수 있는 반면에 해석시간이 많이 소요되고 해의 안정적 획득 측면에서 단점이 있다. 따라서 강소성 유한요소해석과 탄소성 유한요소해석을 상황에 따라 적절하게 사용하는 것이 바람직하다.

8.3.1 탄소성변형 문제의 정의

그림 8.3은 탄소성 재료의 유한요소법, 즉 탄소성 유한요소법을 전개하기 위한 소성가공 공정의 개념도이다. 소성가공 공정에서 금형과 소재는 해석영역 \(V\)와 경계 \(S\)로 나눌 수 있으며, 경계 \(S\) 는 표면력이 \(t_i^{(n)} = \bar{t}_i\)로 이루어진 \(S_t\), 변위가 \(u_i = \bar{u}_i\)로 주어진 \(S_u\), 금형과 소재의 접촉면을 표현하는 \(S_c\) 등으로 표현할 수 있다. 다음의 이론적 전개에서 소재는, 등방탄소성(isotropic elasto-plastic), 속도무의존성(rate independence), 등방성경화(isotropic hardening) 등의 성질을 가지고 있으며 von Mises 항복이론과 직교유동법칙(associated flow rule)을 따르는 경우에 국한한다. 또한, 가속도의 영향은 무시할 수 있을 정도의 준정적(quasi-static)현상이 공정 중에 발생한다고 가정하여 정적 내연적 유한요소법(static implicit finite element method)에 대하여 기술한다.

그림 8.3 탄소성 유동해석 문제의 정의

탄소성 재료의 소성가공 공정 문제는 소재의 변위 \(u_i\)를 구하는 것이며, 다음과 같은 경계치 문제로 수식화된다.

(1) 평형방정식

\[\sigma_{ij,j} + f_i = 0 \quad \text{in} \quad V \tag{8.18}\]

(2) 변형률증분 분해

\[d\varepsilon_{ij} = d\varepsilon_{ij}^e + d\varepsilon_{ij}^p \tag{8.19}\]

(3) 구성방정식

\[\text{탄성} : d\varepsilon_{ij}^e = \frac{d\sigma'_{ij}}{2G} + \frac{1-2\nu}{E} d\sigma_m \delta_{ij} \tag{8.20}\]

\[\text{소성} : d\varepsilon_{ij}^p = d\lambda \frac{\partial f}{\partial \sigma_{ij}} = \frac{3}{2} d\lambda \frac{\sigma'_{ij}}{\bar{\sigma}} \tag{8.21}\]

(4) 변위-변형률증분 관계식

\[d\varepsilon_{ij} = \frac{1}{2} \left( \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i} \right) \tag{8.22}\]

(5) 항복식

\[f = \sigma_e - \sigma_y = \left( \frac{3}{2} \sigma'_{ij} \sigma'_{ij} \right)^{\frac{1}{2}} - \sigma_y = 0 \tag{8.23}\]

(6) 비압축성 조건

\[d\varepsilon_{ii}^p = 0 \tag{8.24}\]

(7) 경계조건

\[\sigma_{ij}n_j = \bar{t}_i \quad \text{on} \quad S_t \tag{8.25}\]

\[u_i = \bar{u}_i \quad \text{on} \quad S_u \tag{8.26}\]

\[u_n = \bar{u}_n \quad \text{on} \quad S_c \tag{8.27}\]

\[\sigma_t = -\mu\sigma_n g(v_t) \quad \text{on} \quad S_c \tag{8.28}\]

여기서 상첨자 \(p\) 와 \(e\) 는 각각 소성성분과 탄성성분을 나타내는 지수이다. 또, 하첨자 \(n\) 과 \(t\) 는 각각 법선성분과 접선성분을 나타내는 지수이다. 변형률증분 \(d\varepsilon_{ij}\) 는 재료의 회전을 고려하여 구해지며 한 계산증분의 크기가 작아 선형화하여 계산할 수 있는 경우로 국한하여 설명한다. 식 (8.28)의 \(g(v_t)\) 는 마찰과 속도와의 관계를 반영하는 함수이며 앞 절에서 설명하고 있다. 또한 식 (8.22)의 적용은 소변형률증분(small strain increment)의 경우로 제한되며 대변형률증분(large strain increment) 및 대회전(large rotation)을 수반하는 대변형(large deformation)을 고려하는 방법에 대해서는 8.3.7에서 기술한다.

8.3.2 von Mises 항복식

식 (8.21)의 \(d\lambda\) 를 von Mises 항복식의 경우에는 간단하게 구할 수 있다. 식 (8.21)에서

\[d\varepsilon_{ij}^p = \frac{3}{2} d\lambda \frac{\sigma'_{ij}}{\bar{\sigma}} \tag{8.29}\]

이며, 소성유효변형률 증분 \(d\bar{\varepsilon}^p\) 는 다음 식으로 주어진다.

\[d\bar{\varepsilon}^p = \left( \frac{2}{3} d\varepsilon_{ij}^p d\varepsilon_{ij}^p \right)^{\frac{1}{2}} \tag{8.30}\]

식 (8.29)를 식 (8.30)에 대입하면,

\[d\bar{\varepsilon}^p = \left( \frac{2}{3} \frac{3}{2} d\lambda \frac{\sigma'_{ij}}{\bar{\sigma}} \frac{3}{2} d\lambda \frac{\sigma'_{ij}}{\bar{\sigma}} \right)^{\frac{1}{2}} = d\lambda \frac{\left( \frac{3}{2} \sigma'_{ij} \sigma'_{ij} \right)^{\frac{1}{2}}}{\bar{\sigma}} = d\lambda \tag{8.31}\]

즉, von Mises항복식의 경우에는 \(d\lambda\) 가 유효소성변형률증분 \(d\bar{\varepsilon}^p\) 와 동일함을 알 수 있다. 따라서 식 (8.21)의 소성유동법칙은 다음과 같이 다시 쓸 수 있다.

\[d\varepsilon_{ij}^p = \frac{3}{2} d\bar{\varepsilon}^p \frac{\sigma'_{ij}}{\bar{\sigma}} \tag{8.32}\]

8.3.3 등방성 가공경화 및 일관성조건

대부분의 금속은 소성변형 과정을 거치면서 경화가 된다. 즉, 소성가공을 더 발생시키기 위해서는 더 큰 응력이 필요하게 된다. 유효소성변형률, \(\bar{\varepsilon}^p\) 는 다음 식과 같이 변형경로를 따라 적분된다.

\[\bar{\varepsilon}^p = \int d\bar{\varepsilon}^p \tag{8.33}\]

이 때 \(d\bar{\varepsilon}^p\) 는 식 (8.30)으로 주어진다. 그림 8.4는 평면응력조건에서 가공경화가 발생하는 단축인장에서의 변형률-응력 곡선과 이에 따른 초기 항복곡면과 가공경화에 따라 점점 팽창하는 항복곡면을 나타내고 있다. 이 때 부호가 다른 방향을 포함한 다른 모든 방향에서도 같은 크기의 응력에서 항복이 발생한다면 등방성이 성립하며 재료가 등방성가공경화의 특성을 가진다고 한다. 그림 8.4에서 초기에 2 축 방향으로 응력이 증가하여 초기항복곡면에 이르기까지는 탄성변형만이 발생한다. 점차 응력이 증가하여 \(\sigma_2 = \sigma_y\) 가 되면 항복이 발생한다. 변형이 더욱 진행함에 따라 소재는 가공경화하게 되며 가공경화량을 유효소성변형률 \(\bar{\varepsilon}^p\) 로 표현할 수 있으므로 식 (8.23)의 항복식을 다음 식과 같이 표현할 수 있다.

\[f(\sigma, \bar{\varepsilon}^p) = \sigma_e - \sigma_y(\bar{\varepsilon}^p) = 0 \tag{8.34}\]

그림 8.4 등방성가공경화

여기서 \({\sigma}\) 는 응력 텐서이다. 한편, 단축인장에서는 위 식은 다음 식과 동일하게 된다.

\[\sigma_e - \sigma_y(\bar{\varepsilon}^p) = \sigma_2 - \sigma_y(\varepsilon_2^p) = 0 \tag{8.35}\]

따라서 그림 8.4에서 \(\varepsilon_2^p (=\bar{\varepsilon}^p)\) 가 증가함에 따라 유동응력, \(\sigma_2 (=\sigma_e)\) 는 증가시켜야 유효응력에 의해 크기가 결정되는 항복곡면 또한 항복식을 만족하도록 증가하여야 항복이 지속된다. 그림 8.4에서 \(r\) 은 가공경화량을 나타낸다. 이와 같이 가공경화가 발생할 때 지속적으로 소성변형이 발생하기 위해서는 응력이 항복곡면 위에 존재해야 하며 이를 일관성 조건(consistency condition)이라고 한다.

8.3.4 탄소성 유한요소법의 수식화

변위장이 주어진 경계에서 0이 되는 임의의 변위의 변분 \(\delta{\omega}\) 에 대해서 변위장 \(\mathbf{u}\) 를 구하기 위한 약형(weak form)은 다음과 같다.

\[\int_V \delta\varepsilon_{ij}\sigma_{ij} dV - \int_V \delta\omega_i f_i dV - \int_{S_t} \delta\omega_i t_i dS - \int_{S_c} \delta\omega_t \sigma_t dS = 0 \tag{8.36}\]

여기서 \(f_i\) 와 \(t_i\) 는 각각 체적력 성분과 표면력 성분을 나타낸다.

유한요소 내의 변위 및 가중함수(또는 가상변위)를 보간함수(또는 형상함수)를 이용하여 다음과 같이 근사화한다.

\[\mathbf{u} = \mathbf{N}\mathbf{U} \tag{8.37}\]

\[{\omega} = \mathbf{N}\mathbf{W} \tag{8.38}\]

여기서 \(\mathbf{N}\) 은 요소 내에서 변위를 근사화하는 보간함수를 포함하는 행렬이다. \(\mathbf{U}\) 와 \(\mathbf{W}\) 는 변위와 가중함수의 절점 값을 나타내는 벡터이다.

응력 텐서 \(l{\sigma} = [\sigma_{ij}]\) 와 변형률 텐서 \({\varepsilon} = [\varepsilon_{ij}]\) 를 다음 식과 같이 벡터 형태로 나타내어 유한요소 수식화에 활용한다.

\[\tilde{{\sigma}} = \begin{bmatrix} \sigma_{xx} \\ \sigma_{yy} \\ \sigma_{zz} \\ \sigma_{xy} \\ \sigma_{yz} \\ \sigma_{zx} \end{bmatrix}, \quad \tilde{{\varepsilon}} = \begin{bmatrix} \varepsilon_{xx} \\ \varepsilon_{yy} \\ \varepsilon_{zz} \\ \varepsilon_{xy} \\ \varepsilon_{yz} \\ \varepsilon_{zx} \end{bmatrix} \tag{8.39}\]

변형률 증분은 식 (8.39)과 식 (8.37)을 이용하여 다음 식으로 나타낼 수 있다.

\[d\tilde{{\varepsilon}} = \begin{bmatrix} d\varepsilon_{xx} \\ d\varepsilon_{yy} \\ d\varepsilon_{zz} \\ 2d\varepsilon_{xy} \\ 2d\varepsilon_{yz} \\ 2d\varepsilon_{zx} \end{bmatrix} = \frac{\partial \mathbf{N}}{\partial \mathbf{x}} \mathbf{U} = \mathbf{B}\mathbf{U} \tag{8.40}\]

탄소성 해석에서 응력증분 \(\Delta{\sigma}\) 는 다음 형태로 주어진다.

\[d\tilde{{\sigma}} = \tilde{\mathbf{C}}^{ep} d\tilde{{\varepsilon}} \tag{8.41}\]

여기서 \(\tilde{\mathbf{C}}^{ep}\) 는 식 (8.61)에서 4차텐서인 \(\mathbf{C}^{ep}\) 의 차수를 2차로 낮추어 사용하였다. 식 (8.36)은 다음과 같이 정리된다.

\[\int_V \delta\tilde{{\varepsilon}}^T \tilde{{\sigma}} dV - \int_V \delta{\omega}^T \mathbf{f} dV - \int_{S_t} \delta{\omega}^T \mathbf{t}^{(n)} dS - \int_{S_c} \delta\omega_t \tilde{\sigma}_t dS = 0 \tag{8.42}\]

탄소성해석에서는 비선형성을 고려하기 위하여 여러 개의 계산증분으로 나누어서 순차적으로 계산해야 한다. 시간 \(t_n\) 에서 해가 구해졌을 때 시간 \(t_{n+1}\) 에서 해를 구하기 위하여 시간증분 \(\Delta t = t_{n+1} - t_n\) 구간에서의 변형에 대하여 식 (8.42)에 식 (8.41)를 적용하여 다음 식을 얻을 수 있다.

\[\int_V \delta \tilde{\varepsilon}^T (\tilde{\sigma}_n + \tilde{\mathbf{C}}^{ep} d\tilde{\varepsilon}) dV - \int_V \delta \omega^T \mathbf{f} dV - \int_{S_t} \delta \omega^T \mathbf{t}^{(n)} dS - \int_{S_c} \delta \omega_t \tilde{\sigma}_t dS = 0 \tag{8.43}\]

식 (8.43)에 식 (8.37), (8.38), (8.40)을 적용하면 다음과 같이 정리된다.

\[\int_V \delta(\mathbf{B}\mathbf{W})^T (\tilde{\sigma}_n + \tilde{\mathbf{C}}^{ep}\mathbf{B}\mathbf{U}) dV - \int_V \delta(\mathbf{N}\mathbf{W})^T \mathbf{f} dV - \int_{S_t} \delta(\mathbf{N}\mathbf{W})^T \mathbf{t}^{(n)} dS - \int_{S_c} \delta(\mathbf{N}\mathbf{W})^T \mathbf{t} \tilde{\sigma}_t dS = 0 \tag{8.44}\]

\[\delta \mathbf{W}^T \left[ \int_V \mathbf{B}^T (\tilde{\sigma}_n + \tilde{\mathbf{C}}^{ep} \mathbf{B}\mathbf{U}) dV - \int_V \mathbf{N}^T \mathbf{f} dV - \int_{S_t} \mathbf{N}^T \mathbf{t}^{(n)} dS - \int_{S_c} \mathbf{N}^T \mathbf{t} \tilde{\sigma}_t dS \right] = 0 \tag{8.45}\]

임의의 \(\delta \mathbf{W}\) 에 대하여 성립해야 하므로 다음 식을 얻는다.

\[\int_V \mathbf{B}^T (\tilde{\sigma}_n + \tilde{\mathbf{C}}^{ep} \mathbf{B}\mathbf{U}) dV - \int_V \mathbf{N}^T \mathbf{f} dV - \int_{S_t} \mathbf{N}^T \mathbf{t}^{(n)} dS - \int_{S_c} \mathbf{N}^T \mathbf{t} \tilde{\sigma}_t dS = 0 \tag{8.46}\]

식 (8.46)을 정리하면 다음 유한요소방정식을 얻는다.

\[[\mathbf{K}]\{\mathbf{U}\} + \{\mathbf{F}^n\} - \{\mathbf{F}^{ext}\} = 0 \tag{8.47}\]

여기서

\[[\mathbf{K}] = \int_V \mathbf{B}^T \tilde{\mathbf{C}}^{ep} \mathbf{B} dV \tag{8.48}\]

\[\mathbf{F}^{ext} = \int_V \mathbf{N}^T \mathbf{f} dV + \int_{S_t} \mathbf{N}^T \tilde{\mathbf{t}} dS + \int_{S_c} \mathbf{N}^T \mathbf{t} \tilde{\sigma}_t dS \tag{8.49}\]

\[\mathbf{F}^n = \int_V \mathbf{B}^T \tilde{\sigma}_n dV \tag{8.50}\]

이다. 소성가공 공정은 여러 가지의 비선형성으로 인해 식 (8.47)을 만족시키는 해를 Newton-Raphson 반복법을 적용하여 구한다. 식 (8.47)에서 발생하는 잔류힘(residual force)은 다음 식으로 표현된다.

\[\Psi(\mathbf{U}) = [\mathbf{K}]\{\mathbf{U}\} + \{\mathbf{F}^n\} - \{\mathbf{F}^{ext}\} \tag{8.51}\]

Taylor 방법으로 식 (8.51)을 근사하면 다음 식을 얻는다.

\[\Psi(\mathbf{U}) + \frac{\partial \Psi(\mathbf{U})}{\partial \mathbf{U}} \Delta\mathbf{U} + O(\Delta\mathbf{U}^2) = 0 \tag{8.52}\]

여기서 행렬 \(\mathbf{J} = \partial \Psi / \partial \mathbf{U}\) 는 자코비안(Jacobian) 또는 유효접선강성행렬로 호칭한다. 식 (8.52)에서 1차항만 사용하여 선형화하면 다음 식을 얻게 된다.

\[\frac{\partial \Psi(\mathbf{U})}{\partial \mathbf{U}} \Delta\mathbf{U} = \mathbf{J}(\mathbf{U}) \Delta\mathbf{U} = -\Psi(\mathbf{U}) \tag{8.53}\]

식 (8.53)에서의 반복계산은 식 (8.51)의 잔류힘이 미리 정해 놓은 한계치 이하로 감소할 때까지 계속한다. 반복계산 과정에서 변위 \(\mathbf{U}\) 는 다음과 같이 갱신된다. 여기서 \(m\) 은 반복계산횟수를 나타낸다.

\[\mathbf{J}(\mathbf{U}_m) \Delta\mathbf{U}_m = -\Psi(\mathbf{U}_m) \tag{8.54}\]

\[\mathbf{U}_{m+1} = \mathbf{U}_m + \Delta\mathbf{U}_m \tag{8.55}\]

8.3.5 응력 및 상태변수의 적분

수치적 해를 구할 때 소성변형이 진전됨에 따라 관련 식들을 만족하도록 응력, 변형률, 상태변수 등이 변형시간에 따라 적분되어야 한다. 이 때 적분은 외연적(explicit)적분방법과 내연적(implicit)적분방법으로 수행될 수 있다.

응력증분과 이에 따른 유효소성변형률 증분이 발생하였을 때, 일관성을 유지하기 위한 조건은 다음과 같이 기술할 수 있다.

\[f(\sigma + d\sigma, \bar{\varepsilon}^p + d\bar{\varepsilon}^p) = 0 \tag{8.56}\]

식 (8.56)를 전개하면

\[f(\sigma + d\sigma, \bar{\varepsilon}^p + d\bar{\varepsilon}^p) = f(\sigma, \bar{\varepsilon}^p) + \frac{\partial f}{\partial \sigma} : d\sigma + \frac{\partial f}{\partial \bar{\varepsilon}^p} d\bar{\varepsilon}^p = 0 \tag{8.57}\]

이다. 여기서 연산자 ( : )은 텐서축약(tensor contraction)을 수행하여 스칼라 값을 계산한다. 식 (8.57)이 성립하기 위해서는 다음 식을 만족해야 한다.

\[\mathbf{n} : d\sigma + \frac{\partial f}{\partial \bar{\varepsilon}^p} d\bar{\varepsilon}^p = 0 \tag{8.58}\]

여기서

\[\mathbf{n} = \frac{\partial f}{\partial \sigma} \tag{8.59}\]

이다. Hooke의 법칙을 이용하여 응력증분과 탄성변형률증분 사이의 관계식을 구할 수 있다.

\[d\sigma = \mathbf{C} d\varepsilon^e = \mathbf{C}(d\varepsilon - d\varepsilon^p) \tag{8.60}\]

여기서 \(\mathbf{C}\) 는 응력텐서와 변형률텐서의 탄성변형관계를 나타내는 4차 텐서이다. 식 (8.21)의 소성변형률증분, \(d\varepsilon^p\) 의 표현을 이용하면 다음과 같이 구할 수 있다.

\[d\sigma = \mathbf{C} d\varepsilon^e = \mathbf{C}(d\varepsilon - d\lambda \mathbf{n}) \tag{8.61}\]

그리고 식 (8.58)에 식 (8.31)과 식 (8.61)를 적용하면 다음 식을 얻을 수 있다.

\[\mathbf{n} : \{\mathbf{C}(d\varepsilon - d\lambda \mathbf{n})\} + \frac{\partial f}{\partial \bar{\varepsilon}^p} d\lambda = 0 \tag{8.62}\]

이 식으로부터 \(d\lambda\) 에 대한 식을 다음과 같이 구할 수 있다.

\[d\lambda = \frac{\mathbf{n} : \mathbf{C} d\varepsilon}{\mathbf{n} : \mathbf{C}\mathbf{n} - \left( \frac{\partial f}{\partial \bar{\varepsilon}^p} \right)} \tag{8.63}\]

식 (8.61)에 식 (8.63)을 대입하여 정리하면 응력증분에 대한 다음 식을 얻는다.

\[d\sigma = \mathbf{C} \left( d\varepsilon - \mathbf{n} \frac{\mathbf{n} : \mathbf{C} d\varepsilon}{\mathbf{n} : \mathbf{C}\mathbf{n} - \left( \frac{\partial f}{\partial \bar{\varepsilon}^p} \right)} \right) = \left( \mathbf{C} - \frac{(\mathbf{C} : \mathbf{n}) \otimes (\mathbf{C} : \mathbf{n})}{\mathbf{n} : \mathbf{C}\mathbf{n} - \left( \frac{\partial f}{\partial \bar{\varepsilon}^p} \right)} \right) d\varepsilon \tag{8.64}\]

또는

\[d\sigma = \mathbf{C}^{ep} d\varepsilon \tag{8.65}\]

이다. 여기서 \(\mathbf{C}^{ep}\) 는 탄소성접선강성을 나타낸다. 연산자 ( \(\otimes\) )는 텐서곱(tensor product or dyadic)을 나타낸다. 소성변형이 없을 경우에는 \(d\lambda = 0\) 이 되며 이 때는 \(\mathbf{C}^{ep} = \mathbf{C}\) 가 되어 탄성변형의 경우로 복귀한다.

식 (8.65)에서 전체 변형증분 \(d\varepsilon\) 이 주어지고, 전 단계시간 \(t_n\) 에서의 응력 \(\sigma_{t_n}\) 및 가공경화상태변수 \(\bar{\varepsilon}_{t_n}^p\) 이 주어졌을 때 다음 시간 단계 \(t_{n+1}\) 에서의 물리량들을 시간 \(t_n\) 에서의 값들만으로 계산하는 방법을 외연적 적분방법이라고 한다. 위의 식들을 이용하여 필요한 값들을 단계시간 \(t_n\) 에서의 값들을 사용하여 계산한다.

\[d\lambda_{t_n} = \frac{\mathbf{n}_{t_n} : \mathbf{C} d\varepsilon_{t_n}}{\mathbf{n}_{t_n} : \mathbf{C}\mathbf{n}_{t_n} - \left( \frac{\partial f}{\partial \bar{\varepsilon}^p} \right)_{t_n}} \tag{8.66}\]

\[d\sigma_{t_n} = \mathbf{C} \left( d\varepsilon_{t_n} - d\lambda_{t_n} \mathbf{n}_{t_n} \right) \tag{8.67}\]

\[dr_{t_n} = \frac{\partial f}{\partial \bar{\varepsilon}_{t_n}^p} d\bar{\varepsilon}_{t_n}^p \tag{8.68}\]

시간적분에 의해 시간 \(t_{n+1}\) 에서의 값들을 다음과 같이 구한다.

\[\sigma_{t_{n+1}} = \sigma_{t_n} + d\sigma_{t_n} \tag{8.69}\]

\[\varepsilon_{t_{n+1}}^p = \varepsilon_{t_n}^p + d\varepsilon_{t_n}^p \tag{8.70}\]

\[r_{t_{n+1}} = r_{t_n} + dr_{t_n} \tag{8.71}\]

이러한 적분방법은 1차 전진오일러(forward Euler) 외연적적분방법이며 간단하게 구현할 수 있는 반면에 여러 가지 단점들을 가지고 있다. 시간증분이 충분히 작아야 해가 불안정해지지 않으며 정확성도 시간증분의 크기에 많이 영향을 받는다. 또한, 최종 응력 상태가 항복곡면에 존재해야하는 일관성 조건을 만족한다는 보장이 없다.

따라서 시간 \(t_{n+1}\) 에서 일관성조건이 만족되도록 하기 위해서는 내연적(implicit) 적분방법을 사용하여야 한다. 그림 8.5는 내연적(implicit)적분방법을 사용하여 시간 \(t_{n+1}\) 에서 응력상태가 von Mises 항복곡면 위에 있도록 하는 과정을 보여 준다. 3차원 응력공간에서 von Mises 항복곡면은 정수압(hydrostatic stress)을 나타내는 직선에 중심축이 일치하는 원통면으로 나타나게 된다. 그림 8.5는 응력공간을 그 직선에 평행한 방향으로 투영하여 도시하였으며 이 때 von Mises 항복곡면은 원으로 나타난다. 내연적 적분 과정에서는 먼저 탄성적으로 응력, \(\mathbf{\sigma}_{n+1}^{tr} (= \mathbf{\sigma}_{n+1}^{trial})\) 을 예측한다. \(\mathbf{\sigma}_{n+1}^{tr}\) 은 일반적으로 가공경화에 의하여 크기가 변화한 시간에서의 항복곡면과 일치하지 않는다. 이 때 von Mises 항복곡면의 경우는 그림 8.5와 같이 원으로 나타난 항복곡면의 중심을 향하여 복귀시켜 간단히 새로운 항복곡면과 일치하게 할 수 있다. 이 때 응력수정성분은 소성수정성분이며, 그림 8.5에서 \(2G\Delta\varepsilon^p\) 에 해당한다. 이러한 방법은 반경방향복귀법(radial return method) [8.1] 라고 하며 von Mises 항복곡면의 경우 매우 간단하게 수행될 수 있다. 아래에 반경방향복귀법에 대하여 설명하며 인식편의성을 위하여 다음 수식들에서 시간 \(t_{n+1}\) 에 해당하는 항의 경우 시간첨자를 생략한다.

그림 8.5 von Mises 항복식의 t_n에서 t_n+1 동안의 내연적 응력적분 과정

Hooke의 법칙을 이용하여 탄성응력에 대한 식을 다음과 같이 구할 수 있다.

\[\sigma = 2G\varepsilon^e + \lambda Tr(\varepsilon^e)\mathbf{I} \tag{8.72}\]

여기서 \(\mathbf{I}\) 는 단위행렬을 의미한다. 시간 \(t_{n+1}\) 에서 탄성변형률을 다음과 같이 쓸 수 있다.

\[\varepsilon^e = \varepsilon_{t_n}^e + \Delta\varepsilon^e = \varepsilon_{t_n}^e + \Delta\varepsilon - \Delta\varepsilon^p \tag{8.73}\]

식 (8.73)을 식 (8.72)에 대입하면

\[\sigma = 2G(\varepsilon_{t_n}^e + \Delta\varepsilon) + \lambda Tr(\varepsilon_{t_n}^e + \Delta\varepsilon)\mathbf{I} - 2G\Delta\varepsilon^p = \sigma^{tr} - 2G\Delta\bar{\varepsilon}^p \mathbf{n} \tag{8.74}\]

이다. 식 (8.74)는 정수압응력성분과 편차응력성분, \(\sigma'\) 으로 나누어 다음과 같이 쓸 수 있다.

\[\sigma' + \frac{1}{3}(\sigma:\mathbf{I}) = \sigma'^{tr} - 2G\Delta\bar{\varepsilon}^p \frac{3}{2}\frac{\sigma'}{\sigma_e} \tag{8.75}\]

식 (8.75)를 정리하면

\[\left( 1 + 3G\frac{\Delta\bar{\varepsilon}^p}{\sigma_e} \right)\sigma' = \sigma'^{tr} \tag{8.76}\]

이다. 위 식의 양변을 제곱하고, 식 (8.23)을 이용하여 정리하면 다음 식을 얻을 수 있다.

\[\sigma_e + 3G\Delta\bar{\varepsilon}^p = \sigma_e^{tr} \tag{8.77}\]

이 식으로부터 항복식은 다음과 같이 유도된다.

\[f = \sigma_e^{tr} - 3G\Delta\bar{\varepsilon}^p - \sigma_y(\bar{\varepsilon}^p + \Delta\bar{\varepsilon}^p) = 0 \tag{8.78}\]

식 (8.78)을 선형화하면

\[\sigma_e^{tr} - 3G\Delta\bar{\varepsilon}^p - \sigma_y(\bar{\varepsilon}^p + \Delta\bar{\varepsilon}^p) + \left( -3G - \left( \frac{\partial \sigma_y}{\partial \bar{\varepsilon}^p} \right) \right) d\Delta\bar{\varepsilon}^p = 0 \tag{8.79}\]

이다. 이 식을 \(d\Delta\bar{\varepsilon}^p\) 에 대하여 정리하면

\[d\Delta\bar{\varepsilon}^p = \frac{\sigma_e^{tr} - 3G\Delta\bar{\varepsilon}^p - \sigma_y(\bar{\varepsilon}^p + \Delta\bar{\varepsilon}^p)}{3G + \left( \frac{\partial \sigma_y}{\partial \bar{\varepsilon}^p} \right)} \tag{8.80}\]

이다. 선형화된 식들을 반복계산하여 식 (8.78)을 만족하는 \(d\Delta\bar{\varepsilon}^p\) 를 구한다. 다음에 식 (8.76)과 식 (8.77)을 이용하여 소성변형률증분을 다음 식과 같이 구할 수 있다.

\[\Delta\varepsilon^p = \Delta\bar{\varepsilon}^p \frac{3}{2} \frac{\sigma'}{\sigma_e} = \Delta\bar{\varepsilon}^p \frac{3}{2} \frac{\sigma'^{tr}}{\sigma_e^{tr}} \tag{8.81}\]

탄성변형률증분 \(\Delta\varepsilon^e\) 를 계산한 뒤 응력증분을 Hooke의 법칙을 이용하여 계산한다.

\[\Delta\varepsilon^e = \Delta\varepsilon - \Delta\varepsilon^p \tag{8.82}\]

\[\Delta\sigma = 2G\Delta\varepsilon^e + \lambda \mathbf{I} \Delta\varepsilon^e : \mathbf{I} \tag{8.83}\]

시간 단계 \(t_{n+1}\) 에서 응력 및 소성변형률은 다음 식에 의하여 적분하여 구한다.

\[ \begin{aligned} \sigma_{t_{n+1}} &= \sigma_{t_n} + \Delta\sigma \\ \bar{\varepsilon}^p &= \bar{\varepsilon}_{t_n}^p + \Delta\bar{\varepsilon}^p \end{aligned} \tag{8.84} \]

8.3.6 탄소성재료특성의 프로그램 구현

내연적탄소성유한요소해석(static implicit elasto-plastic finite element) 프로그램에서는 기본적으로 식 (8.54)과 식 (8.55)의 식을 이용하여 시간열 \(\mathbf{t} = \{t_0(=0), t_1, t_2, \cdots, t_{n+1}, \cdots, t_{end}\}\) 에 따라 순차적으로 탄소성문제를 해석해 나간다. \(t_n\) 까지 해석이 수행되었을 때 다음 시간단계 \(t_{n+1}\) 까지 해석을 수행하기 위해서는 식 (8.54)과 식 (8.55)의 식들을 구성해서 해석을 수행해야 한다. 유한요소법은 각 요소행렬들을 계산하여 총합하여 전체 문제를 풀어 나간다. 이 때 요소수준에서 요소행렬들을 계산해야 하며 요소행렬들의 계산은 사용하는 요소의 특성에 따라 다르게 구현된다. 한편, 요소행렬들의 계산은 적분점을 이용한 수치적분을 사용하여 수행되며 재료의 탄소성 특성의 반영은 적분점 수준에서 구현된다. 따라서 요소의 종류가 상이하여도 적분점 수준에서 탄소성재료특성의 구현은 동일하다. 따라서 이 절에서는 적분점 수준에서 구현되는 von Mises 항복식을 따르는 재료의 탄소성 특성이 구현되는 과정에 대하여 설명한다.

시간단계 \(t_n\) 에서 응력 \(\sigma_{t_n}\) 및 가공경화상태변수인 \(\bar{\varepsilon}_n^p\) 이 주어졌을 때 다음 시간단계 \(t_{n+1}\) 까지 요소의 적분점 수준에서 내연적 적분법을 사용하여 von Mises 항복식을 따르는 재료의 탄소성 프로그램에서 처리하는 과정을 순차적으로 정리하였다.

(1) 먼저 변형률증분 \(\Delta\varepsilon\) 를 계산한다.

\[\Delta\tilde{\varepsilon} = \mathbf{B}\mathbf{U} \tag{8.85}\]

(2) 탄성 예측응력(elastic trial stress)를 계산한다.

\[\sigma^{tr} = \sigma_{t_n} + 2G\Delta\varepsilon + \lambda\mathbf{I}\Delta\varepsilon : \mathbf{I} \tag{8.86}\]

(3) 항복식을 계산한다.

\[f = \sigma_e^{tr} - \sigma_y = \left( \frac{3}{2} \sigma'^{tr} : \sigma'^{tr} \right)^{\frac{1}{2}} - \sigma_y \left( \bar{\varepsilon}^p \right) \tag{8.87}\]

(4) 항복 여부 판정에 따라 다음 절차를 따라 계산한다.

(i) \(f > 0\) 이면, 유효소성변형률증분, \(d\Delta\bar{\varepsilon}^p\) 를 다음 식들의 반복계산으로 구한다. 이 때 \(d\Delta\bar{\varepsilon}^p\) 가 정해진 오차 이하로 감소할 때까지 반복한다.

\[\bar{\varepsilon}^{p(k)} = \bar{\varepsilon}_{t_n}^p + \Delta\bar{\varepsilon}^{p(k)}\]

\[\sigma_y^{(k)} = \sigma_y \left( \bar{\varepsilon}^{p(k)} \right)\]

\[d\Delta\bar{\varepsilon}^p = \frac{\sigma_e^{tr} - 3G\Delta\bar{\varepsilon}^{p(k)} - \sigma_y^{(k)}}{3G + \frac{\partial \sigma_y}{\partial \bar{\varepsilon}^{p(k)}}} \tag{8.88}\]

\[\Delta\bar{\varepsilon}^{p(k+1)} = \Delta\bar{\varepsilon}^{p(k)} + d\Delta\bar{\varepsilon}^p\]

(ii) 그 이외의 경우에는 소성변형이 발생하지 않으므로

\[\Delta\bar{\varepsilon}^p = 0 \tag{8.89}\]

(5) 소성변형률증분, 탄성변형률증분, 응력증분 등을 계산한다.

\[ \begin{aligned} \Delta\varepsilon^p &= \frac{3}{2} \Delta\bar{\varepsilon}^p \frac{\sigma'^{tr}}{\sigma_e^{tr}} \\ \Delta\varepsilon^e &= \Delta\varepsilon - \Delta\varepsilon^p \\ \Delta\sigma &= 2G\Delta\varepsilon^e + \lambda\mathbf{I}\Delta\varepsilon^e : \mathbf{I} \end{aligned} \tag{8.90} \]

시간 \(t_{n+1}\) 에서의 응력과 유효소성변형률을 적분하여 구한다.

\[ \begin{aligned} \sigma_{t_{n+1}} &= \sigma_{t_n} + \Delta\sigma \\ \bar{\varepsilon}^p &= \bar{\varepsilon}_{t_n}^p + \Delta\bar{\varepsilon}^p \end{aligned} \tag{8.91} \]

유한요소방정식 (8.54)를 구성하는 데 필요한 적분점에서의 유효접선강성행렬에의 기여항, \(\mathbf{J}\) 를 계산한다.

\[\mathbf{J} = \frac{\partial \mathbf{\Psi}}{\partial \mathbf{U}} \tag{8.92}\]

8.3.7 대변형의 고려

대변형이 수반된 변형은 대회전 및 대변형률을 적절하게 고려해야 한다. 이를 위한 방법은 여러 가지가 있으나 변형구배텐서(deformation gradient tensor) \(\mathbf{F}\) 를 이용한 방법을 소개한다. 매 해석스텝에서 변형구배텐서 \(\mathbf{F}\) 를 구한 뒤, 이것에 극분해(polar decomposition)를 적용하여 회전텐서(rotation tensor) \(\mathbf{R}\) 와 우신장텐서(right stretch tensor) \(\mathbf{U}\) 를 다음 식과 같이 구한다.

\[\mathbf{F} = \mathbf{R}\mathbf{U} \tag{8.93}\]

식 (8.93)의 우신장텐서 \(\mathbf{U}\) 로부터 대변형률을 고려하기 위하여 변형률증분, \(\Delta\varepsilon_{ij}\) 을 다음 식으로 구한다.

\[\Delta\varepsilon_{ij} = \ln(U_{ij}) \tag{8.94}\]

한 단계의 해석스텝 동안의 von Mises 항복식을 따르는 응력증분 \(\Delta\mathbf{\sigma}\) 의 계산은 8.3.5절에 자세히 설명되어 있다. 회전성분을 포함하지 않는 변형률 증분, \(\Delta\varepsilon_{ij}\) 을 사용하여 시간 \(t_{n+1}\) 에서 회전을 고려하지 않은 응력 \(\mathbf{\sigma}^*_{n+1}\) 을 다음 식으로 구한다.

\[\sigma^*_{n+1} = \sigma_n + \Delta\sigma = \sigma_n + \mathbf{C}(\Delta\varepsilon - \Delta\varepsilon^p) = \mathbf{C}^{ep}\Delta\varepsilon \tag{8.95}\]

여기서 \(\sigma_n\) 은 시간 \(t_n\) 에서의 응력을 나타내며, 직전 해석스텝 계산단계에서 이미 구해진 값이다. \(\mathbf{C}\) 는 탄성강성이며 \(\Delta\varepsilon^p\) 는 소성변형률증분이다. 식 (8.95)에서 해를 구하기 위한 반복계산 과정에서의 수렴성 향상을 위하여 탄소성일관접선강성(elastoplastic consistent tangent stiffness) \(\mathbf{C}^{ep}\) [8.2]을 사용하였다. 탄소성일관접선강성과 응력증분 \(\Delta\sigma\) 를 구하는 과정은 항복식의 종류 및 \(\Delta\varepsilon^p\) 를 구하는 방법에 따라 달라진다.

식 (8.95)에 의하여 구해진 응력 \(\sigma^*_{n+1}\) 을 식 (8.93)에서 구한 회전텐서 \(\mathbf{R}\) 을 사용하여 변환하여 회전을 고려한 해석스텝 \(t_{n+1}\) 에서의 응력, \(\sigma_{n+1}\) 을 구한다.

\[\sigma_{n+1} = \mathbf{R}\sigma^*_{n+1}\mathbf{R}^T \tag{8.96}\]

8.3.8 3차원 탄소성해석에 MINI-요소의 활용

육면체요소를 사용하는 것이 수치적 안정성과 해석결과의 정확도에서 유리하다는 사실은 일반적으로 알려져 있는 사실이나, 소성가공 공정의 해석을 위한 육면체요소망의 자동생성은 용이하지 않다. 사면체요소는 요소망생성 및 재구성이 용이한 관계로 금속성형해석에 많이 사용되고 있다. 그러나 사면체요소를 사용할 때는 비압축성 조건이 요소에 너무 과하게 작용하여 수치적 로킹 현상(locking phenomenon)이 발생하게 된다. 이러한 사면체요소 사용 시에 발생하는 문제점을 해결하기 위해서 절점에 정수압 자유도를 추가하여 수치적 로킹 문제를 해결할 수 있다[8.3].

일반적으로 속도장이 주어진 경계 \(S_u\), 표면력(traction)이 주어진 경계 \(S_t\) 를 가지는 변형되는 물체 \(V\) 에 대한 탄소성 유한요소해석을 위한 2 변수(변위, 정수압) 방법에 의한 약형을 유도한다. 주어진 변위경계조건을 만족하는 임의의 변위의 변분 \(\delta\mathbf{\omega}\) 와 임의의 정수압의 변분 \(\delta p\) 에 대해서 변위장 \(\mathbf{u}\) 와 정수압 \(p\) 를 구하기 위한 약형은 다음과 같다.

\[\int_V \delta\mathbf{\varepsilon}^T \mathbf{\sigma} dV - \int_V \delta\mathbf{\omega}^T \mathbf{f} dV - \int_{S_t} \delta\mathbf{\omega}^T \mathbf{t}^{(n)} dS = 0 \tag{8.97}\]

\[\int_V \delta q \left[ \varepsilon_v + \frac{p}{K} \right] dV = 0 \tag{8.98}\]

여기서 \(\sigma\) 과 \(\varepsilon\) 은 각각 응력 텐서와 변형률 텐서를 나타내고, \(\mathbf{f}\) 와 \(\mathbf{t}^{(n)}\) 는 각각 체적력과 표면력 벡터를 나타낸다. 그리고 \(\varepsilon_v\) 는 부피변화율을 나타내고, \(K\) 는 체적탄성계수(bulk modulus of elasticity)를 나타낸다. 한편, \(\mathbf{m}^T = [1 \ 1 \ 1 \ 0 \ 0 \ 0]\) 로 정의하면, 부피변화율은 다음과 같이 표현된다.

\[\varepsilon_v = \mathbf{m}^T \varepsilon \tag{8.99}\]

탄소성 유한요소해석에서 비선형성을 고려하기 위하여 여러 개의 계산증분으로 나누어서 순차적으로 계산한다. 시간 \(t_n\) 에서 해가 구해졌을 때, 시간 \(t_{n+1}\) 에서 해는 식 (8.97)와 식 (8.98)에 식 (8.95)과 식 (8.99)을 적용함으로써 구해진다. 즉, 다음의 수식으로부터 시간 \(t_{n+1}\) 에서 해를 구할 수 있다.

\[\int_V \mathbf{B}^T (\sigma'_n + \mathbf{D}_d \mathbf{B}\Delta\mathbf{U}) dV - \int_V \mathbf{B}^T \mathbf{m} \mathbf{N}_p \mathbf{P} dV - \int_V \mathbf{N}_u^T \mathbf{f} dV - \int_{S_t} \mathbf{N}_u^T \mathbf{t}^{(n)} dS = 0 \tag{8.100}\]

\[-\int_V \mathbf{N}_p^T \mathbf{m}^T \varepsilon_n dV - \int_V \mathbf{N}_p^T \mathbf{m}^T \mathbf{B}\Delta\mathbf{U} dV - \int_V \mathbf{N}_p^T \frac{1}{K} \mathbf{N}_p \mathbf{P} dV = 0 \tag{8.101}\]

여기서 \(\sigma'\) 는 편차응력을, \(\mathbf{D}_d\) 는 편차응력-변형률 관계 행렬을, \(\mathbf{B}\) 는 변형률-변위 관계행렬을 나타낸다. 사면체요소에서 변위장 근사 형상함수 \(\mathbf{N}_u\) 와 정수압 근사 형상함수 \(\mathbf{N}_p\) 를 동일한 차수로 보간하면 정수압이 진동하는 현상이 발생하게 된다. 이를 해결하기 위하여 MINI-요소[8.4]에서는 버블 형상함수, \(\mathbf{N}_B\) 와 버블 변위 \(\Delta\mathbf{U}^B\) 를 추가하여 변위장의 차수를 높여 주어 이러한 현상을 해결한다.

그리고 식 (8.100)와 식 (8.101)을 Newton-Raphson반복법과 버블 형상함수를 적용하면 다음 식을 얻는다.

\[ \begin{bmatrix} \mathbf{K}_0 & \mathbf{h}_{ub} & \mathbf{h}_{up} \\ \mathbf{h}_{ub}^T & \mathbf{h}_{bb} & \mathbf{h}_{bp} \\ \mathbf{h}_{up}^T & \mathbf{h}_{bp}^T & \mathbf{V} \end{bmatrix} \begin{pmatrix} \Delta\mathbf{U} \\ \Delta\mathbf{U}^B \\ \Delta\mathbf{P} \end{pmatrix} = \begin{pmatrix} \mathbf{f}^{ext} \\ 0 \\ 0 \end{pmatrix} - \begin{pmatrix} \mathbf{f}^{int} \\ 0 \\ \mathbf{f}^p \end{pmatrix} \tag{8.102} \]

여기서 첨자 0은 표준 형상함수와 관련된 항을 나타내고 첨자 \(b\) 는 버블 형상함수와 관련된 항을 나타내며 첨자 \(p\) 는 정수압과 관련된 항을 나타낸다. 식 (10)의 강성행렬의 부행렬은 다음과 같다.

\[ \begin{aligned} \mathbf{K}_0 &= \int_V \mathbf{B}^{0^T} \mathbf{D}_d \mathbf{B}^0 dV \\ \mathbf{h}_{ub} &= \int_V \mathbf{B}^{0^T} \mathbf{D}_d \mathbf{B}^b dV \\ \mathbf{h}_{bu} &= \mathbf{h}_{ub}^T = \int_V \mathbf{B}^{b^T} \mathbf{D}_d \mathbf{B}^0 dV \\ \mathbf{h}_{bb} &= \int_V \mathbf{B}^{b^T} \mathbf{D}_d \mathbf{B}^b dV \\ \mathbf{h}_{up} &= -\int_V \mathbf{B}^{0^T} \mathbf{m} \mathbf{N} dV \\ \mathbf{h}_{bp} &= -\int_V \mathbf{B}^{b^T} \mathbf{m} \mathbf{N} dV \\ \mathbf{V} &= -\int_V \mathbf{N}^T \frac{1}{K} \mathbf{N} dV \\ \mathbf{f}^{ext} &= \int_V \mathbf{N}^T \mathbf{f} dV + \int_{S_t} \mathbf{N}^T \mathbf{t}^{(n)} dS \\ \mathbf{f}^{int} &= \int_V \mathbf{B}^{0^T} \{ \sigma'_{n+1} + \mathbf{m} \mathbf{N}_p \mathbf{P} \} dV \\ \mathbf{f}^p &= -\int_V \mathbf{N}^T \left( \mathbf{m}^T \varepsilon + \frac{1}{K} \mathbf{N} \mathbf{P} \right) dV \end{aligned} \tag{8.103} \]

식 (8.102)에서 버블 변위 \(\Delta\mathbf{U}^B\) 를 소거하면 다음 유한요소방정식을 얻는다.

\[ \begin{bmatrix} \mathbf{K}_0 - \mathbf{h}_{ub} \mathbf{h}_{bb}^{-1} \mathbf{h}_{ub}^T & \mathbf{h}_{up} - \mathbf{h}_{ub} \mathbf{h}_{bb}^{-1} \mathbf{h}_{bp} \\ \mathbf{h}_{up}^T - \mathbf{h}_{bp}^T \mathbf{h}_{bb}^{-1} \mathbf{h}_{ub}^T & \mathbf{V} - \mathbf{h}_{bp}^T \mathbf{h}_{bb}^{-1} \mathbf{h}_{bp} \end{bmatrix} \begin{pmatrix} \Delta\mathbf{U} \\ \Delta\mathbf{P} \end{pmatrix} = \begin{pmatrix} \mathbf{f}^{ext} \\ 0 \end{pmatrix} - \begin{pmatrix} \mathbf{f}^{int} \\ \mathbf{f}^p \end{pmatrix} \tag{8.104} \]

소성가공 공정해석에서는 여러 가지 비선형성이 존재하므로 식 (8.104)를 만족시키는 해를 8.3.4와 동일한 방법으로 Newton-Raphson 반복법을 적용하여 구한다.