• 초기 강성법(Initial stiffness), 뉴튼 랩슨법(Newton Raphson), 수정(Modified)뉴튼 랩슨법 비선형해석에서의 반복법은 접선강성의 계산 시점에 따라 초기 강성법, 뉴튼 랩슨법, 수정 뉴튼 랩슨법 등으로 분류할 수 있다. 초기 강성법에서는 해석 시작시점에서 계산된 접선강성을 계속 유지하고, 뉴튼 랩슨법에서는 매 반복 계산마다 접선강성을 재계산하며, 수정 뉴튼 랩슨법은 외력의 변화가 발생하는 시점에서 접선강성을 계산한다. 접선강성행렬의 계산 및 행렬분해는 많은 계산시간을요구하므로 초기 강성법과 수정 뉴튼 랩슨 법을 활용할 경우, 수렴 과정에 문제가 발생하지 않는다면 뉴튼 랩슨 법에 비해 빠른 계산 속도를 얻을 수 있다.midas NFX에서는 초기 강성법 또는 뉴튼 랩슨법 등을 명시적으로 구분하지 않고, 접선강성의 재계산 시점을 정의함으로써 모든 반복법의 효과를 얻을 수 있다. • 자동 강성행렬 재계산 해석하고자 하는 모델의 비선형 정도, 수렴해의 평탄함 정도 등 해당 문제 특징에 따라 적절한 접선강성의 계산 시점을 선택하는 것이 중요하다. midas NFX에서는 비선형 유한요소 해법으로 해당 비선형 문제의 특징, 즉 반복계산시의 수렴특징, 발산여부 등을 종합적으로 판단하여 적절한 시점에 강성행렬을 재계산하는 자동 강성행렬 재계산 방법(automatic tangential stiffness update)을 제공한다. 기본적으로 다음과 같은 몇 가지 조건이 만족되는 경우 강성행렬 재계산이수행된다. ► 사용자가 지정한 최대 반복계산 횟수보다 예상되는 반복계산 횟수가 많은 경우 ► 해가 발산하는 것으로 판정된 경우 • 수렴조건, 발산 판정 및 하중 이등분(Load bisection) 반복법에서 해의 수렴여부는 부재력 기준(force norm), 변위기준(displacementnorm), 그리고 에너지 기준(energy norm) 으로 판정한다. $$ \text { Force norm ratio } = \frac {\sqrt {\mathbf {g} _ {i} ^ {T} \mathbf {g} _ {i}}}{\sqrt {\Delta \mathbf {f} _ {\text { int } , i} ^ {T} \Delta \mathbf {f} _ {\text { int } , i}}} \tag {5.5.9} $$ $$ \text { Displacement norm ratio } = \frac {\sqrt {\delta \mathbf {u} _ {i} ^ {T} \delta \mathbf {u} _ {i}}}{\sqrt {\Delta \mathbf {u} _ {i} ^ {T} \Delta \mathbf {u} _ {i}}} \tag {5.5.10} $$ $$ \text { Energy magnitude ratio } = \left| \frac {\delta \mathbf {u} _ {i} ^ {T} \mathbf {g} _ {i}}{\Delta \mathbf {u} _ {i} ^ {T} \Delta \mathbf {f} _ {\text { int } , i}} \right| \tag {5.5.11} $$ midas NFX에서는 이 세 가지 기준 중 하나 또는 다수의 기준에 대하여 사용자가 제공하는 허용치와의 비교를 통해 수렴여부를 판단한다. • 발산 판정 및 하중단계 이등분 해의 발산 여부는 자동 강성행렬 재계산 방법에서 중요한 기준으로 사용되며,발산율(divergence rate) 을 기초로 결정한다. $$ E _ {i} = \frac {\delta \mathbf {u} _ {i} ^ {T} \mathbf {g} _ {i}}{\delta \mathbf {u} _ {i} ^ {T} \mathbf {g} _ {i - 1}} \tag {5.5.12} $$ 발산율의 절대값이 1보다 큰 경우( ), 비선형해석의 해는 발산의 위험이 있다고 판단할 수 있으며, 강성행렬 재계산 또는 하중 이등분 등의 알고리즘 상필요한 조치가 취해진다. 하중 이등분은 해가 발산하거나 수렴을 위한 반복계산 수가 사용자가 지정한 최대 횟수보다 많아지는 등, 현재 하중단계의 증분이 수렴해를 얻기에 너무 크다고 판정되는 경우에 적용된다. 기존 하중 증분을 이등분하여 반복계산을 재시작함으로써, 부적절한 하중 증분 크기에 대해 유연하게 대처할 수 있다. midasNFX에서는 사용자가 지정하는 최대 이등분 단계(maximum bisection level)에 도달할 때까지 필요에 따라 하중 이등분 자동 수행하게 된다. • 시간 증분 자동조절 비선형 해석의 효율성을 증가시키기 위하여, midas NFX에는 시간 증분의 크기를비선형해석의 수렴성을 기초로 자동 조절하는 기능을 포함한다. 기본이 되는 시간 증분 크기와 최대 증분 크기는 사용자의 입력에 의하여 결정된다. 비선형 해석에서 시간증분 자동 조절 기능을 사용할 경우, 이전 증분단계에서 수렴하는데필요로 하는 반복계산 횟수를 기초로 다음 단계의 시간 증분크기는 증가 또는 감소된다. $$ \Delta t ^ {i + 1} = n _ {s} \Delta t ^ {i} \quad (1 \leq n _ {s} \leq n _ {s, \max}) \tag {5.5.13} $$ 여기서 증분 조정계수 ( )는 자연수로 국한시킴으로써 사용자가 의도한 시점또는 하중 크기에서 비선형 해를 최대한 얻을 수 있도록 하였다. 증분 조정계수는 초기 증분량을 의미하는 1을 최소로, 사용자가 제공하는 값 $\mathit { \Pi } _ { . \mathit { n } _ { s , \operatorname* { m a x } } }$ )을 최대로하는 범위를 갖는다. # 5.6 대변형을 고려한 변형률/응력 산출법 기하학적 선형 해석에서는 변형 전과 변형 후의 형상 차이를 고려하지 않고 변형률과 응력을 정의한다. 일반적으로 기하학적 선형 해석에서의 변형률은 다음과 같이 정의된다. $$ \boldsymbol {\varepsilon} = \frac {1}{2} \left[ \frac {\partial \mathbf {u}}{\partial \mathbf {X}} + \left(\frac {\partial \mathbf {u}}{\partial \mathbf {X}}\right) ^ {T} \right] \tag {5.6.1} $$ u : 변위 X : 변형 전 또는 변형 후 좌표 대변형을 고려한 기하학적 비선형 해석에서는 다양한 방법으로 변형률을 정의할수 있으며, 가상일(virtual work)을 정의할 수 있도록 각각의 변형률에 대응하는응력이 존재한다. • 변형률의 정의 midas NFX에서 사용하는 변형률로는 그린(Green) 변형률 또는 그린-라그랑지(Green-Lagrange) 변형률과 변형 속도(rate of deformation) 또는 변형률 속도(strain rate)가 있다. 그린 변형률 텐서 는 다음과 같이 정의된다. $$ d s ^ {2} - d S ^ {2} = 2 d \mathbf {X} \cdot \mathbf {E} \cdot d \mathbf {X} \tag {5.6.2} $$ X : 구조물 특정 위치의 변형 전 좌표 그린 변형률은 변형 전 미소 길이 의 제곱과 변형 후 미소길이 제곱의차이에 관한 값으로 볼 수 있다. 그린 변형률 텐서는 다음과 같이 변형구배(deformation gradient)에 의해 정의할 수 있다. $$ \mathbf {E} = \frac {1}{2} (\mathbf {F} ^ {T} \cdot \mathbf {F} - \mathbf {I}) \tag {5.6.3} $$ 강체운동만이 있는 경우 그린 변형률은 발생하지 않기 때문에, 변형을 측정하는값으로 타당하다. midas NFX에서는 초탄성 재료와 같이 에너지 포텐셜이 존재하 는 경우에 대하여, 계산의 편이를 도모하기 위해 그린 변형률을 사용하지만 사용자에게 제공되는 결과값은 다음과 같이 주 신장률(principal stretch)과 그 방향을 이용하여 계산한 변형률이다. $$ \boldsymbol {\varepsilon} = \ln \left(\lambda_ {1}\right) \mathbf {n} _ {1} \mathbf {n} _ {1} + \ln \left(\lambda_ {2}\right) \mathbf {n} _ {2} \mathbf {n} _ {2} + \ln \left(\lambda_ {3}\right) \mathbf {n} _ {3} \mathbf {n} _ {3} \tag {5.6.4} $$ : 주 신장률 n : 주 신장 방향 벡터 (변형 후) 초탄성 재료 이외의 모든 재료에 대하여 사용되는 변형속도는 속도구배(velocitygradient)에 의해 다음과 같이 정의된다. $$ \mathbf {D} = \frac {1}{2} (\mathbf {L} + \mathbf {L} ^ {T}) = \text { sym } [ \frac {\partial \mathbf {v}}{\partial \mathbf {x}} ] \tag {5.6.5} $$ V : 속도 벡터 X : 구조물 특정 위치의 변형 후 좌표 즉 속도구배 텐서의 대칭 부분에 해당한다. 변형속도는 미소길이 제곱의 변화율에 대한 값으로 볼 수 있다. $$ \frac {\partial d s ^ {2}}{\partial t} = 2 d \mathbf {x} \cdot \mathbf {D} \cdot d \mathbf {x} \tag {5.6.6} $$ 변형속도 역시 강체운동만이 존재하는 경우 발생하지 않으며, 그린 변형률과 다음의 관계를 가진다. $$ \mathbf {D} = \mathbf {F} ^ {- T} \cdot \dot {\mathbf {E}} \cdot \mathbf {F} ^ {- 1} \tag {5.6.7} $$ 변형속도는 시간에 대한 변화율이기 때문에 이를 시간적분하여 변형률로 사용하는 것이 일반적이다. midas NFX에서 초탄성 재료 이외의 재료에 대해 기하학적비선형성을 고려한 해석을 수행하게 되면 변형속도를 시간에 대해 적분하여 변형률로 산출한다. • 응력의 정의 기하학적으로 변형이 커지게 되면 응력 또한 다양한 방법으로 정의할 수 있다.midas NFX에서 사용하는 응력의 종류로는 코시 응력(Cauchy stress)와 2ndPK(Piola-Kirchhoff) 응력이 있다. 코시 응력은 현재 형상에 대해 평형 방정식을 만족하는 값이기 때문에 진응력(true stress)이라 불리기도 하며, 다음과 같이 정의된다. $$ \mathbf {n} \cdot \boldsymbol {\sigma} \cdot d \Gamma = d \mathbf {f} = \mathbf {t} d \Gamma \tag {5.6.8} $$ 변형 전 형상에 대해 코시응력은 다음과 같이 2nd PK 응력으로 변환할 수 있다. $$ \mathbf {S} = \boldsymbol {J} \mathbf {F} ^ {- 1} \cdot \boldsymbol {\sigma} \cdot \mathbf {F} ^ {- T} \tag {5.6.9} $$ 

text_image

n₀ df dΓ₀

Reference configuration 

text_image

n df dΓ F⁻¹df

Current configuration 그림 5.6.1 응력 정의를 위한 변형 전/후 형상과 힘의 방향 2nd PK 응력은 물리적 의미가 명쾌하지 않으나 그린 변형률과 짝을 이루어 운동 방정식을 기술하는 데 유용하기 때문에 에너지 포텐셜을 가지는 고무 재료등의 거동을 정의하는 데 주로 사용된다. 특히 초탄성 재료의 경우 그린 변형률로부터 2nd PK 응력을 계산할 수 있기 때문에 별도의 응력 적분법이 필요하지않다. 그러나, 초탄성 재료의 경우에도 사용자에게 제공되는 응력값은 (5.6.9)의 역변환을 통하여 계산된 코시 응력이다. 그 밖의 모든 재료에 대하여는 다음에설명하는 응력과 변형률의 적분 방법을 이용한다. • 응력 속도와 변형률 속도의 적분 탄소성 재료, 점탄성 재료 등은 에너지 포텐셜이 주어지지 않는 반면 변형률 속도(strain rate)와 객관 응력 속도(objective stress rate)의 관계식으로 표현되는 구성방정식을 이용한다. midas NFX에서 사용하는 Jaumann 응력속도는 다음과 같이 정의할 수 있다. $$ \dot {\boldsymbol {\sigma}} ^ {J} = \dot {\boldsymbol {\sigma}} - \mathbf {w} \cdot \boldsymbol {\sigma} - \boldsymbol {\sigma} \cdot \mathbf {w} ^ {T} \tag {5.6.10} $$ 변형률 속도와 객관 응력속도는 재료의 구성 방정식에 의해 관계를 가진다. $$ \dot {\boldsymbol {\sigma}} ^ {J} = \mathbf {C}: \mathbf {D} \tag {5.6.10} $$ 스텝 에서 계산된 응력과 변형률 증분을 이용하여 스텝 에서의 응력을 계산하기 위해 (5.6.10)에 중앙 차분법(central difference)을 반영하고 구조물의 회전을 고려한다. $$ \boldsymbol {\sigma} _ {n + 1} = \Delta \mathbf {R} \cdot \boldsymbol {\sigma} _ {n} \cdot \Delta \mathbf {R} ^ {T} + \mathbf {C}: \Delta \boldsymbol {\varepsilon} \tag {5.6.11} $$ 회전량 증분 은 다음과 같이8 계산함으로써 객관 응력 증분(incrementallyobjective) 조건을 만족할 수 있다. $$ \Delta \mathbf {R} = (\mathbf {I} - \frac {1}{2} \Delta \mathbf {W}) ^ {- 1} (\mathbf {I} + \frac {1}{2} \Delta \mathbf {W}) \tag {5.6.12} $$ 특히 변형률 증분과 증분 스핀(spin) 계산은 형상에서 수행함에 주의해야 한다. $$ \Delta \boldsymbol {\varepsilon} = \frac {1}{2} \left(\frac {\partial \Delta \mathbf {u}}{\partial \mathbf {x} _ {n + 1 / 2}} + \left[ \frac {\partial \Delta \mathbf {u}}{\partial \mathbf {x} _ {n + 1 / 2}} \right] ^ {T}\right), \Delta \mathbf {w} = \frac {1}{2} \left(\frac {\partial \Delta \mathbf {u}}{\partial \mathbf {x} _ {n + 1 / 2}} - \left[ \frac {\partial \Delta \mathbf {u}}{\partial \mathbf {x} _ {n + 1 / 2}} \right] ^ {T}\right) \tag {5.6.13} $$ 변형률 속도의 적분 역시 (5.6.11)과 동일하게 구조물의 회전량 증분 을 이용하여 계산한다. # 5.7 비선형 동적응답 해법 midas NFX는 내연적 시간적분법(implicit time integration) 또는 외연적 시간적분법(explicit time integration)을 사용하는 비선형 과도응답해석(nonlinear transient response analysis)을 모두 지원한다. # 5.7.1 외연적 시간 적분법 midas NFX의 비선형 과도응답해석에서는 외연적 시간적분법 중 가장 널리 알려진 중앙차분법(central difference)을 이용하며 집중질량(lumped mass)을 사용하고 있다. # - 중앙차분법 시간축에 대해 중앙차분법을 적용하기 위해 시간 스텝을 $n+1/2$ 과 $n, n+1$ 스텝으로 구분한다. $$ \Delta t ^ {n + 1 / 2} = t ^ {n + 1} - t ^ {n}, t ^ {n + 1 / 2} = \frac {1}{2} \left(t ^ {n + 1} + t ^ {n}\right), \Delta t ^ {n} = t ^ {n + 1 / 2} - t ^ {n - 1 / 2} \tag {5.7.1} $$ 스텝 $n+1$ 에서의 변위는 $n+1/2$ 스텍에서의 속도로부터 계산할 수 있다. $$ \dot {\mathbf {u}} ^ {n + 1 / 2} = \mathbf {v} ^ {n + 1 / 2} = \frac {1}{\Delta t ^ {n + 1 / 2}} \left(\mathbf {u} ^ {n + 1} - \mathbf {u} ^ {n}\right), \quad \mathbf {u} ^ {n + 1} = \mathbf {u} ^ {n} + \Delta t ^ {n + 1 / 2} \mathbf {v} ^ {n + 1 / 2} \tag {5.7.2} $$ 스텝 $n+1/2$ 에서의 속도는 n 스텝에서 가속도로부터 다음과 같이 계산한다. $$ \ddot {\mathbf {u}} ^ {n} = \mathbf {a} ^ {n} = \frac {1}{\Delta t ^ {n}} \left(\mathbf {v} ^ {n + 1 / 2} - \mathbf {v} ^ {n - 1 / 2}\right), \quad \mathbf {v} ^ {n + 1 / 2} = \mathbf {v} ^ {n - 1 / 2} + \Delta t ^ {n} \mathbf {a} ^ {n} \tag {5.7.3} $$ 가속도 $a^{n}$ 은 공간상의 이산화(discretization)을 통하여 다음과 같이 계산할 수 있다. $$ \mathbf {M} \mathbf {a} ^ {n} = \mathbf {f} ^ {n} = \mathbf {f} ^ {\text { ext }} (\mathbf {u} ^ {n}, t ^ {n}) - \mathbf {f} ^ {\text { int }} (\mathbf {u} ^ {n}, t ^ {n}), \mathbf {a} ^ {n} = \mathbf {M} ^ {- 1} \left(\mathbf {f} ^ {\text { ext }} (\mathbf {u} ^ {n}, t ^ {n}) - \mathbf {f} ^ {\text { int }} (\mathbf {u} ^ {n}, t ^ {n})\right) \tag {5.7.4} $$ 결국, 중앙차분법을 이용한 외연적 시간 적분은 (5.7.2)에서 (5.7.4)까지의 반복계산을 의미한다. 또한 대각행렬 형태의 질량 행렬을 사용함으로써, (5.7.4)에서 역행렬을 계산하는 비용은 무시할 수 있을 정도가 된다. 이러한 절차를 시간 증분에 대하여 반복함으로써 구조물의 과도 특성, 준정적 특성 등, 일반적인 동적 거동을 효과적으로 모사할 수 있다. # - 임계 시간스텝(Critical time step) 외연적 시간적분법은 알고리즘이 비교적 간단하며 강건한 특성을 갖는데 반해, 조건적 안정성(conditional stability)을 갖는다. 즉, 시간스텝의 크기가 임계값을 초과하였을 경우, 외연적 시간적분법에 의한 해는 발산하게 된다. 여기서 임계시간스텝(critical time step) 또는 최소 안정시간스텝(minimum stable time step) 은해석 모델에 포함된 모든 요소 각각의 안정시간스텝 중 가장 작은 값을 기준으로 계산된다. $$ \Delta t = \alpha \Delta t _ {c r i t}, \Delta t _ {c r i t} = \frac {2}{\omega_ {\max}} \leq \min _ {e} \left\{\Delta t _ {e} \right\} \tag {5.7.5} $$ 여기서 $\omega_{max}$ 는 선형화된 시스템 전체의 최대 주파수이며, $\Delta t_{e}$ 는 각 요소의 최소 안정시간스텝이다. 또한 $\alpha$ 는 외연적 시간적분법의 안정성을 향상시키기 위하여 도입되는 스케일 팩터로 0.85를 기본값으로 갖는다. 각 요소의 최소 안정시간스텝은 일반적으로 다음과 같이 계산할 수 있다. $$ \Delta t _ {e} = \frac {2}{\omega_ {\max} ^ {e}} = \min \left\{\frac {L _ {e}}{c _ {d}} \right\} \tag {5.7.7} $$ 여기서 $L_{e}$ 는 요소의 특성길이(characteristic length)이며, 일반적으로 요소 내부를 통과하는 가장 짧은 직선 길이로 볼 수 있다. 팽창파 속도(dilation wave speed) $c_{d}$ 는 재료 물성치에 의해 계산되며, 예를 들어 3차원 요소에 대해 다음과 같다.