Newton-Raphson 법¶

선형화와 반복 점화식¶

외력의 가상 일과 전체 방정식의 조립에서 얻은, 시각 \(t_{n+1}\)에서의 절점 변위 \(\boldsymbol{u}_{n+1}\)에 관한 비선형 방정식을 Newton-Raphson 법으로 푼다. 시각 \(t_n\)까지의 절점 변위 \(\boldsymbol{u}_n\)은 알려져 있다고 하고, 변위 증분 \(\Delta\boldsymbol{u}\)를 미지수로 하여

\[ \boldsymbol{u}_{n+1} = \boldsymbol{u}_n + \Delta\boldsymbol{u} \]

를 구한다. 이후에는 외력 벡터의 절점 변위 의존성을 무시하고, \(\boldsymbol{F}(\boldsymbol{u}_{n+1}) = \boldsymbol{F}_{n+1}\)로 두어

\[ \boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u}) - \boldsymbol{F}_{n+1} = \boldsymbol{0} \]

를 푼다.

현재 해 \(\Delta\boldsymbol{u}\)에서의 접선 강성

\[ \boldsymbol{K} = \left. \frac{\partial \boldsymbol{Q}}{\partial \boldsymbol{u}} \right|_{\boldsymbol{u}_n + \Delta\boldsymbol{u}} \]

을 사용하여 비선형 방정식을 선형화하면

\[ \boldsymbol{K}\, d\boldsymbol{u} + \boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u}) - \boldsymbol{F}_{n+1} = \boldsymbol{0} \]

이 된다. 제 \(i\) 반복의 수정량을 \(d\boldsymbol{u}_i\), 반복 시작 시의 잔차 벡터를

\[ \boldsymbol{R}_{i-1} = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u}) \]

라고 쓰면, 반복 점화식은

\[ \boldsymbol{K}_i\, d\boldsymbol{u}_i = \boldsymbol{R}_{i-1}, \qquad \Delta\boldsymbol{u} \leftarrow \Delta\boldsymbol{u} + d\boldsymbol{u}_i \]

의 형태가 된다. 잔차 \(\boldsymbol{R}_i\)는 평형 상태에서 벗어난 힘의 불균형에 해당하는 양이다.

접선 강성 행렬의 구성¶

접선 강성 \(\boldsymbol{K} = \partial\boldsymbol{Q}/\partial\boldsymbol{u}\)는 내력의 가상 일의 이산화에서 얻은 요소 내력 벡터를 절점 변위로 편미분하고, 요소별 피적분항을 요소 영역에서 적분한 것을 모아 구성한다. 요소 수준의 피적분항을 \(\boldsymbol{K}^e_X\)(기준 배치 표기, TL 법) 또는 \(\boldsymbol{K}^e_x\)(현재 배치 표기, UL 법)라고 쓰면, 요소 접선 강성은

\[ \boldsymbol{K}^e = \int_{\Omega^e_0} \boldsymbol{K}^e_X\, dV \quad (\text{TL}), \qquad \boldsymbol{K}^e = \int_{\Omega^e} \boldsymbol{K}^e_x\, dv \quad (\text{UL}) \]

로 주어진다. 다음에 TL/UL 피적분항의 최종 형태를 나타낸다. 어느 경우나 재료 강성 항(초기 변위 항)과 기하 강성 항(초기 응력 항)의 합으로 분해된다.

Total Lagrange 법¶

Total Lagrange 법에서는 제2 Piola-Kirchhoff 응력률 \(\dot{\boldsymbol{S}}\)와 Green-Lagrange 변형률률 \(\dot{\boldsymbol{E}}\)의 선형 관계 \(\dot{\boldsymbol{S}} = \boldsymbol{\mathsf{C}}:\dot{\boldsymbol{E}}\)를 가정한다. 이는 선형 탄성체(St.Venant-Kirchhoff 체)와 초탄성체의 구성 법칙에 대응하며, FrontISTR에서는 Total Lagrange 법을 이러한 재료에 사용한다. 이때 요소 접선 강성의 피적분항은 텐서 형식으로

\[ \delta\boldsymbol{u}^{eT}\, \boldsymbol{K}^e_X\, \dot{\boldsymbol{u}}^e = \dot{\boldsymbol{S}}:\delta\boldsymbol{E} + \boldsymbol{S}:(\delta\boldsymbol{F}^T \dot{\boldsymbol{F}}) \]

라고 쓸 수 있다. 우변의 제1항이 재료 강성 항(초기 변위 항), 제2항이 기하 강성 항(초기 응력 항)이다.

FrontISTR 구현에서는 이 피적분항을 Voigt 표기를 사용한 행렬 형식

\[ \boldsymbol{K}^e_X = (\boldsymbol{B}_L + \boldsymbol{B}_{NL})^T\, \tilde{\boldsymbol{C}}\, (\boldsymbol{B}_L + \boldsymbol{B}_{NL}) + \boldsymbol{F}_9^T\, \boldsymbol{S}_9\, \boldsymbol{F}_9 \]

으로 계산한다. 각 행렬은 다음과 같다. \(\boldsymbol{B}_L, \boldsymbol{B}_{NL}\)은 내력의 가상 일의 이산화에서 도입한 B 행렬이며, \(\tilde{\boldsymbol{C}}\)는 구성 텐서 \(\boldsymbol{\mathsf{C}}\)를 Voigt 표기로 나타낸 \(6\times 6\) 재료 강성 행렬(텐서 표기와 수학적 기초)이다. \(\boldsymbol{S}_9, \boldsymbol{F}_9\)는 기하 강성 항을 행렬곱으로 나타내기 위해 사용하는 다음의 재배열 행렬이다. 먼저 \(3\times 3\)의 2계 텐서 \(\boldsymbol{A}\)를 9성분 벡터로 재배열하는 기호 \([\,\cdot\,]\)를

\[ [\boldsymbol{A}] = (A_{11}, A_{21}, A_{31}, A_{12}, A_{22}, A_{32}, A_{13}, A_{23}, A_{33})^T \]

로 정의하면, \(\boldsymbol{F}_9\)는 변형 구배의 변분을 \([\delta\boldsymbol{F}] = \boldsymbol{F}_9\, \delta\boldsymbol{u}^e\)의 형태로 나타내는 \(9\times d n_e\) 행렬이며, 요소 절점 \(\alpha = 1, \ldots, n_e\)에 대응하는 \(9\times d\) 블록

\[ [\boldsymbol{F}_9]_\alpha = \begin{bmatrix} (\partial N_\alpha^e/\partial X_1)\, \boldsymbol{I} \\ (\partial N_\alpha^e/\partial X_2)\, \boldsymbol{I} \\ (\partial N_\alpha^e/\partial X_3)\, \boldsymbol{I} \end{bmatrix} \qquad (\boldsymbol{I} \text{는 } 3\times 3 \text{ 단위행렬}) \]

을 요소 절점 순서대로 가로로 배열한 \(\boldsymbol{F}_9 = [[\boldsymbol{F}_9]_1, \ldots, [\boldsymbol{F}_9]_{n_e}]\)로 주어진다. \(\boldsymbol{S}_9\)는 이것과 조합하여 기하 강성 항이 \(\delta\boldsymbol{u}^{eT}\, \boldsymbol{F}_9^T \boldsymbol{S}_9 \boldsymbol{F}_9\, \dot{\boldsymbol{u}}^e\)로 표현되도록 선택된 \(9\times 9\) 행렬

\[ \boldsymbol{S}_9 = \begin{bmatrix} S_{11} \boldsymbol{I} & S_{12} \boldsymbol{I} & S_{13} \boldsymbol{I} \\ S_{21} \boldsymbol{I} & S_{22} \boldsymbol{I} & S_{23} \boldsymbol{I} \\ S_{31} \boldsymbol{I} & S_{32} \boldsymbol{I} & S_{33} \boldsymbol{I} \end{bmatrix} \]

이다.

Updated Lagrange 법¶

Updated Lagrange 법에서는 상대 Kirchhoff 응력 텐서의 Jaumann 응력률 \(\hat{\boldsymbol{\sigma}}^{\nabla J}\)와 변형률 속도 텐서 \(\boldsymbol{D}\)의 선형 관계 \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\)를 가정한다. 이는 선형 탄성체·탄소성체·크리프 재료에 공통되는 hypoelastic 구성 법칙의 형식이며, FrontISTR에서는 Updated Lagrange 법을 이러한 재료에 사용한다. 이때 현재 배치로 나타낸 요소 접선 강성의 피적분항은 텐서 형식으로

\[ \delta\boldsymbol{u}^{eT}\, \boldsymbol{K}^e_x\, \dot{\boldsymbol{u}}^e = \boldsymbol{\sigma}^{\nabla T}:\delta\boldsymbol{A}_{(L)} + \boldsymbol{\sigma}:(\delta\boldsymbol{F}_t^T\, \boldsymbol{L}) \]

라고 쓸 수 있다(\(\boldsymbol{\sigma}^{\nabla T}\)는 Truesdell 응력률, \(\boldsymbol{A}_{(L)}\)은 Almansi 변형률의 선형부, \(\boldsymbol{F}_t = \partial\boldsymbol{u}/\partial\boldsymbol{x}\)는 현재 배치를 기준으로 한 변위 구배, \(\boldsymbol{L}\)은 속도 구배 텐서). 우변의 제1항이 재료 강성 항, 제2항이 기하 강성 항이다.

FrontISTR 구현에서는 이 피적분항을 Voigt 표기를 사용한 행렬 형식

\[ \boldsymbol{K}^e_x = \boldsymbol{b}^T\, (\tilde{\boldsymbol{C}} - \boldsymbol{G})\, \boldsymbol{b} + \boldsymbol{f}_9^T\, \boldsymbol{\sigma}_9\, \boldsymbol{f}_9 \]

으로 계산한다. \(\boldsymbol{b}\)는 현재 배치에서 구성한 B 행렬(내력의 가상 일의 이산화)이다. \(\boldsymbol{\sigma}_9, \boldsymbol{f}_9\)는 TL 법에서 정의한 \(\boldsymbol{S}_9, \boldsymbol{F}_9\)에서 제2 PK 응력 \(\boldsymbol{S}\)를 Cauchy 응력 \(\boldsymbol{\sigma}\)로, 기준 배치 구배 \(\partial N_\alpha^e/\partial X_i\)를 현재 배치 구배 \(\partial N_\alpha^e/\partial x_i\)로 치환한 것이다.

\(\boldsymbol{G}\)는 hypoelastic 구성 법칙 \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\)를 Truesdell 응력률 기반의 구성 법칙으로서 접선 강성의 틀과 일치시키는 데 필요한 Cauchy 응력 의존 보정 행렬이며, 4계 텐서 성분 \(G_{ijkl} = \delta_{il}\sigma_{kj} + \delta_{jl}\sigma_{ik}\)를 \(6\times 6\) Voigt 표기로 배열한

\[ \boldsymbol{G} = \begin{bmatrix} 2\sigma_{11} & 0 & 0 & \sigma_{12} & 0 & \sigma_{31} \\ 0 & 2\sigma_{22} & 0 & \sigma_{12} & \sigma_{23} & 0 \\ 0 & 0 & 2\sigma_{33} & 0 & \sigma_{23} & \sigma_{31} \\ \sigma_{12} & \sigma_{12} & 0 & \tfrac{\sigma_{11}+\sigma_{22}}{2} & \tfrac{\sigma_{12}}{2} & \tfrac{\sigma_{23}}{2} \\ 0 & \sigma_{23} & \sigma_{23} & \tfrac{\sigma_{12}}{2} & \tfrac{\sigma_{22}+\sigma_{33}}{2} & \tfrac{\sigma_{12}}{2} \\ \sigma_{31} & 0 & \sigma_{31} & \tfrac{\sigma_{23}}{2} & \tfrac{\sigma_{12}}{2} & \tfrac{\sigma_{33}+\sigma_{11}}{2} \end{bmatrix} \]

이다.

전체 강성 행렬의 조립¶

전체 접선 강성 \(\boldsymbol{K}\)는 요소 강성 \(\boldsymbol{K}^e\)를 절점별 \(d\times d\) 블록 \(\boldsymbol{K}^e_{\alpha\beta}\)로 나누고, 요소 절점 물리량의 조립에서 도입한 2계 텐서용 조립 집합 \(\mathcal{E}^2(i_g, i_h)\)를 사용하여

\[ \boldsymbol{K}_{i_gi_h} = \sum_{(e,\alpha,\beta) \in \mathcal{E}^2(i_g, i_h)} \boldsymbol{K}^e_{\alpha\beta} \]

를 \(i_g\)행 \(i_h\)열에 배치한 행렬로 얻어진다. 구현에서는 집합 \(\mathcal{E}^2\)를 명시적으로 만들지 않고, 요소 루프 안에서 대응하는 블록에 직접 더한다. 절점당 자유도 \(\times\) 전체 절점 수 \(n_g\)차의 정방행렬이 되지만, 요소를 통해 연결된 절점 사이 이외의 성분은 \(0\)이므로 희소 행렬 형식으로 저장된다.

TL 법과 UL 법의 요소 강성은 기준 배치의 전환(절점 좌표와 B 행렬의 구성 원천)과 \(\boldsymbol{G}\) 행렬의 유무를 제외하면 같은 형태이므로, FrontISTR에서는 양자를 공통 서브루틴으로 구현한다.

반복 알고리즘¶

이상을 정리하면, 반복 초기에는 \(\Delta\boldsymbol{u} = \boldsymbol{0}\)으로 두고 초기 잔차 \(\boldsymbol{R}_0 = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n)\)를 구한 다음, 제 \(i\) 반복에서 다음 절차를 수행한다.

현재 변위 \(\boldsymbol{u}_n + \Delta\boldsymbol{u}\)에서의 접선 강성 \(\boldsymbol{K}_i = \partial\boldsymbol{Q}/\partial\boldsymbol{u}|_{\boldsymbol{u}_n + \Delta\boldsymbol{u}}\)를 접선 강성 행렬의 구성의 절차로 계산한다.
기하학적 경계 조건을 반영하기 위해 변위 구속이 부과된 자유도에 대해 접선 강성 행렬과 잔차 벡터를 변형하여 \(\tilde{\boldsymbol{K}}_i, \tilde{\boldsymbol{R}}_{i-1}\)를 얻는다(기하학적 경계 조건의 처리).
선형 방정식 \(\tilde{\boldsymbol{K}}_i\, d\boldsymbol{u}_i = \tilde{\boldsymbol{R}}_{i-1}\)를 풀어 수정량 \(d\boldsymbol{u}_i\)를 구한다. 이 절차는 반복 계산에서 계산 부하의 대부분을 차지하는 경우가 많다.
변위 증분을 \(\Delta\boldsymbol{u} \leftarrow \Delta\boldsymbol{u} + d\boldsymbol{u}_i\)로 갱신하고, 이에 따라 내력 벡터 \(\boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u})\) 및 잔차 \(\boldsymbol{R}_i = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u})\)를 계산한다.
수렴 판정을 수행하고, 수렴하면 반복을 종료한다. 잔차 \(\boldsymbol{R}_i\) 중 기하학적 경계 조건이 부과된 자유도에는 구속 반력에 해당하는 성분이 나타나므로, 이들을 제외한 성분 \(\tilde{\boldsymbol{R}}_i\)로 판정 지표를 구성한다. 구체적인 판정 지표와 임계값은 수렴 판정에서 다룬다. 수렴하지 않고 반복 상한에 도달한 경우에는 반복 실패로 한다.

반복이 수렴한 시점의 \(\Delta\boldsymbol{u}\)를 \(\boldsymbol{u}_n\)에 더하여 시각 \(t_{n+1}\)의 누적 변위 \(\boldsymbol{u}_{n+1} = \boldsymbol{u}_n + \Delta\boldsymbol{u}\)로 하고, 다음 시간 단계로 진행한다.