김경종
19 KiB
c _ {d} = \sqrt {\frac {M}{\rho}} = \sqrt {\frac {K + 4 \mu / 3}{\rho}} \tag {5.7.8}
M : 팽창탄성계수
K , \mu : 체적탄성계수, 전단강성
임계 시간스텝은 물리적인 시간에 비해 굉장히 작은 경우가 많기 때문에 외연적시간적분법의 효과적인 수행을 위해서는 각 요소의 내력 계산 시간이 매우짧아야 한다. 이로 인하여 고차 형상함수를 가지는 요소는 대부분 사용하지않으며, 저차 형상함수를 가지는 요소 역시 가장 계산량이 적은 요소 기법을사용하는 것이 좋다. 표 5.7.1 은 midas NFX 에서 외연적 시간적분법에 사용할수 있는 요소 기법과 절점 수를 정리한 것이다.
표 5.7.1 외연적 동해석에서 사용 가능한 요소 종류
| 요소 종류 | 요소 기법, 절점 수 |
| Rod | 2 절점 |
| Bar | 2 절점 |
| Pipe | 2 절점 |
| Cable | 2 절점 |
| Membrane, Shell | 3 절점, 변위가정법/ANS4 절점, 감차적분(안정화 기법) |
| Surface | 3 절점4 절점, 감차적분6 절점 |
| Solid | 4 절점5 절점, 감차적분6 절점, 감차적분(안정화 기법)8 절점, 감차적분(안정화 기법)10 절점, 비적합요소 |
- 인공 체적점성(Artificial bulk viscosity)
외연적 동해석에서, 특히 고속거동을 하는 구조물의 해석을 위해서는 인공 체적 점성을 이용하여 중앙차분법의 해를 안정화 시켜야 한다. 체적 점성 효과는 요소의 내력(internal force)을 변화시킴으로써 반영할 수 있다.
체적 점성에 의한 압력은 체적 변형률 속도(volumetric strain rate)가 음수인 경우, 즉 압축되는 형태로 속도가 발생하는 경우와 양수인 경우에 따라 다르게 적용되며, 체적 변형률 속도에 선형으로 비례하는 부분과 제곱에 비례하는 부분으로 나눈다. 예를 들어 3차원 요소에 대해 다음과 같이 압력을 계산한다.
p = \left\{ \begin{array}{l l} \rho L _ {e} (C _ {0} L _ {e} \dot {\varepsilon} _ {k k} ^ {2} - C _ {1} c _ {d} \dot {\varepsilon} _ {k k}) & \text { if } \dot {\varepsilon} _ {k k} < 0 \\ \rho L _ {e} (- C _ {1} c _ {d} \dot {\varepsilon} _ {k k}) & \text { if } \dot {\varepsilon} _ {k k} \geq 0 \end{array} \right. \tag {5.7.9}
위식에 의해 계산된 압력은 요소의 응력에 반영되어 내력을 변화시키게 된다. 단, 점성 효과에 의한 압력은 요소의 변형률과 구성방정식으로부터 발생한 응력이 아니기 때문에 내력 계산에만 반영함에 주의해야 한다. C_{0}, C_{1} 는 인공 체적점성 계산을 위한 계수이며 기본값은 각각 1.5, 0.06 이다.
굽힘(bending)을 받는 요소, 예를 들어 bar 또는 shell 요소의 경우에는 굽힘 변형에 대한 인공 체적점성이 (5.7.9)에 추가적으로 반영된다.
- 감쇠(Damping) 효과
외연적 시간적분법에서 고려하는 감쇠의 종류로는 질량비례감쇠(mass-proportional damping), 강성비례감쇠(stiffness-proportional damping) 그리고 구조감쇠(structural damping)가 있다. 외연적 동해석 평형방정식 (5.7.4)에서 감쇠 효과를 고려하면 다음과 같이 표현된다.
\mathbf {M} \mathbf {a} ^ {n} = \mathbf {f} ^ {\text { ext }} (\mathbf {u} ^ {n}, t ^ {n}) - \mathbf {f} ^ {\text { int }} (\mathbf {u} ^ {n}, t ^ {n}) - \mathbf {f} _ {\text { damp }} ^ {n - 1 / 2} \tag {5.7.10}
\mathbf {f} _ {d a m p} ^ {n - 1 / 2} = \sum_ {e} \alpha_ {e} \mathbf {M} _ {e} \mathbf {v} _ {e} ^ {n - 1 / 2} + \sum_ {e} \left(\beta_ {e} + \frac {\eta}{\omega_ {d} ^ {e}}\right) \int_ {\Omega_ {e}} \mathbf {C}: \dot {\varepsilon} d \Omega \tag {5.7.11}
\alpha_{e} : 질량비례 감쇠계수
\beta_{e} : 강성비례 감쇠계수
\eta, \omega_{d}^{e} : 구조감쇠 계수, 지배 주파수
C : 재료의 탄성계수
여기서 질량비례감쇠, 강성비례감쇠와 구조감쇠는 인공체적점성과 유사하게 각요소 차원에서 계산된다.
중앙차분법에서 감쇠는 임계시간스텝을 작게 만드는 경향을 보인다. 일반적으로감쇠가 시간스텝을 변화시키는 효과는 임계 감쇠비를 통하여 반영된다.
\Delta t _ {e} = \frac {2}{\omega_ {\max} ^ {e}} (\sqrt {\xi^ {2} + 1} - \xi) \tag {5.7.12}
임계 감쇠비 는 질량비례 감쇠와 강성비례 및 구조감쇠를 반영하여 다음과 같이 계산한다.
\xi = \frac {\alpha_ {e}}{2 \omega} + \frac {(\beta_ {e} + \eta / \omega_ {d} ^ {e}) \omega}{2} \tag {5.7.13}
감쇠에 의한 임계시간스텝의 변화뿐만 아니라, 앞서 설명한 인공체적점성 또한임계시간스텝의 변화를 초래하게 된다. 일반적으로 감쇠와 체적점성은 시간스텝을 작게 만들기 때문에 불필요하게 큰 값을 사용하는 것을 권장하지 않는다.
• 질량 스케일(Mass scaling)
질량 스케일은 외연적 동해석의 계산비용을 줄이기 위하여 도입된 방법이다. 앞서 설명한 바와 같이 임계 시간스텝이 요소의 크기에 의해 결정되기 때문에, 일부 요소의 크기가 매우 작거나 변형에 의해 작아지는 경우, 이들 요소에 의해전체 해석의 안정시간 스텝이 매우 작게 결정되며, 결과적으로 해석 시간과 비용이 매우 증가하는 경우가 발생한다. 질량 스케일 방법은 이러한 요소의 질량을 인위적으로 증가시킴으로써 팽창파 속도를 작게 하고, 결과적으로 안정시간스텝이 커지도록 하는 역할을 한다.
주의할 부분은 과도한 질량 스케일의 적용으로 인한 시스템 전체 질량 및 관성력의 변화가 동적 특성을 변화시킬 정도 반영되는 경우이다. midas NFX에서는질량 스케일 방법을 적용하는데 있어서 다음과 같은 방법들을 제공한다.
► 전체 모델의 질량을 일관되게 스케일 하는 방법
► 시간스텝을 특정 값으로 맞추는 방법
또한, 질량 스케일에 의한 계산비용의 증가를 감안하여, 질량 스케일의 빈도를조절할 수 있는 기능을 포함한다.
• 조인트 조건 부가법
조인트(Joint) 요소는 요소에 속하는 두 점 간의 상대적 운동을 미리 정의된 형태로 구속하는 요소이다. 외연적 시간적분법을 사용하는 비선형 과도응답해석에서는 벌칙기법(penalty method)을 이용하여 두 절점이 지정된 상대적 거동을 하도록 구속한다. 비선형 과도응답해석에서 조인트 요소는 절점 수에 따라 1절점(그라운드점i , 절점 ), 또는 2절점(절점 )을 갖는 형태로 구별되며, 각각 6개와 12개의 자유도를 갖는다. 조인트 구속 방정식을 C 라고 했을 때, 구속방정식은 조인트 요소 자유도에 의해 표현된다. 이때, 절점 와 j 에 대한 구속력은 개략적으로 다음과 같이 표현된다.
\mathbf {f} _ {i} = - k C \frac {\partial C}{\partial \mathbf {x} _ {i}}, \quad \mathbf {f} _ {j} = - k C \frac {\partial C}{\partial \mathbf {x} _ {j}} \tag {5.7.14}
k : 벌칙계수
C : 조인트 구속 방정식
여기서 벌칙계수 는 외연적 시간적분법의 안정성을 확보하는 동시에 지정된상대적 거동을 정확히 모사할 수 있도록 적절한 값을 사용해야 한다. midas-NFX에서는 현 단계에서의 임계 시간스텝(critical time step)값과 선형화된 조인트운동방정식의 고유치를 이용해서 벌칙계수를 자동으로 계산한다.
5.7.2 내연적 시간 적분법
midas NFX의 비선형 과도응답해석에서는 5.4.1절에서 소개된 HHT- 방법을 내연적 시간적분법으로 사용한다. 비선형 과도응답해석에서의 동적 평형방정식은(5.4.1)과 동일하며, 비선형 유한요소 해법을 이용한 반복계산을 통해서 누적 증분해를 수렴시키는 방법을 사용하여 각 시간스텝에서의 해를 구한다. (5.4.1)의동적 평형방정식으로부터 구한 불평형력은 다음과 같다.
\mathbf {g} _ {n + 1} = \mathbf {M} \ddot {\mathbf {u}} _ {n + 1} + (1 + \alpha_ {H}) \left(\mathbf {C} _ {n + 1} \dot {\mathbf {u}} _ {n + 1} + \mathbf {f} _ {\text { int }, n + 1} - \mathbf {f} _ {\text { ext }, n + 1}\right) - \alpha_ {H} \left(\mathbf {C} _ {n} \dot {\mathbf {u}} _ {n} + \mathbf {f} _ {\text { int }, n} - \mathbf {f} _ {\text { ext }, n}\right) \tag {5.7.15}
내연적 시간 적분법 내의 비선형 해법은 기본적으로 5.5절에 소개된 방법과 동일하며, (5.7.15)의 불평형력이 0이 되는 방향으로 진행된다. 선형 과도응답해석과 마찬가지로 \alpha _ { H } = - 0 . 0 5 \overset { \equiv } { \equiv } 기본값으로 사용한다.
• 자동 시간스텝 제어
내연적 비선형 과도응답해석에서는 중간스텝 불평형력을 기반으로 한 자동 시간스텝 제어방법을 지원한다. 사용자는 모델에 가해지는 외력의 크기에 대한 중간스텝 불평형력의 크기의 비를 입력하여 허용가능한 중간스텝 불평형력의 크기를정할 수 있고, 내연적 비선형 과도응답해석에서는 이를 기반으로하여 시간스텝의 크기를 제어한다.
• 감쇠(Damping) 효과
내연적 시간 적분법에서는 기본적으로 5.4.1절의 선형 동적응답 해석에서 모드감쇠를 제외한 질량비례감쇠, 강성비례감쇠 및 구조감쇠를 고려한다. 이 경우 감쇠 행렬은 (5.4.5)와 같이 구성되지만, 비선형 과도응답해석에서는 대변형에 의한회전효과가 고려된 질량행렬과 재료 비선형성이 고려된 강성행렬을 사용하여 감쇠 행렬을 구성한다.
5.8 접촉 조건
접촉해석(contact analysis)은 공간상의 두 물체가 서로 맞닿을 수는 있으나, 관통할 수 없다는 조건 (non-penetration condition)을 기본 가정으로 하며, 물리적인관점에서 비선형 거동 또는 조건에 해당한다. 접촉의 종류로는 그림 5.8.1과 같이 물체간의 충돌 및 충돌 시의 마찰을 고려해서 해석하는 일반접촉(generalcontact), 미끄러짐을 반영하지 않는 거친접촉(rough contact), 그림 5.8.2와 같이해석 초기에 두 물체가 접착되어 해석이 진행되는 접착접촉(welded contact) 그리고 접선방향의 미끄러짐만을 고려하는 미끄러짐 접촉(sliding contact) 등이 있다. 이 중에서 접착접촉과 미끄러짐 접촉은 두 물체 간의 해석 초기 위치에 따라 조건이 부여되는 것으로서 선형 조건으로 볼 수 있다.
text_image
General contact Rough contact
그림 5.8.1 일반접촉과 거친접촉의 개념
text_image
Welded contact
그림 5.8.2 접착접촉의 개념
• 접촉조건과 해석종류의 관계
midas NFX에서는 선형 구조해석 및 열전달 해석 시에 초기 인접한 물체 간의접착조건 및 미끄러짐 접촉조건을 사용할 수 있으며 접촉조건과 유사한 보간 요소를 사용할 수 있다. 보간요소는 종속절점(slave node)과 주접촉면(mastersegment)의 절점을 연결하여 이들 간의 상대 운동을 구속하는 방법이다.
비선형 해석(정적, 동적)에서는 접착조건, 미끄러짐 접촉조건 그리고 보간요소법이외에 일반접촉과 거친접촉 조건을 사용할 수 있다. 일반접촉과 거친접촉은 비선형 조건에 해당하며 해석 기법상으로 기하학적 비선형성의 고려 여부에 따라그 거동을 달리한다. 기하학적 비선형을 고려하는 경우에는 대변위가 발생할 수있다는 가정에 의해 모든 주접촉면에 대해 접촉 가능성을 고려하는 반면, 기하학적 비선형을 고려하지 않게 되면 초기에 사용자가 지정한 거리 이내에 인접한주접촉면과 종속절점간의 접촉만을 고려하여 계산한다.
접촉의 정의는 점-면 접촉(node to surface contact) 또는 면-면 접촉(surface tosurface contact), 그리고 단일면 접촉(single surface contact)을 선택할 수 있다.이 중 단일면 접촉은 일반접촉과 거친접촉에 대해서만 적용이 가능하다. 점-면접촉은 계산 시간이 적게 든다는 장점이 있으나 주물체의 절점이 종속물체를 관통하는 경향이 크므로 해의 정확성이 상대적으로 떨어진다. 반면에 면-면 접촉은 점-면 접촉에 비해 계산 시간이 많이 소요되나 비관통 조건을 상대적으로 정확하게 만족하므로 구조물의 거동을 비교적 정확하게 모사할 수 있다.
삼차원 문제의 경우 위와 같은 접촉 정의를 따르며, 이차원 또는 축대칭 문제의
우, 해석 상황에 따라 접촉을 다르게 정의한다.
먼저 이차원 외연적 비선형 해석의 경우, 해석 시간의 이점을 위해 삼차원과 동일한 접촉 방식을 따른다. 하지만 그 외 해석의 경우, 모르타 방법(mortarmethod) 이라고 알려진 접촉 해석방식을 사용한다. 이 방법은 접촉 제한조건을절점(node)에 부여하는 방식이 아닌 접촉면(interface)에 부과하는 방식으로, 절점-면 접촉 방식(node to surface contact) 보다 더 정확한 해로의 수렴성을 가지는 것으로 알려져 있다.
접촉면에서의 접촉 제한 조건을 계산하는 방식은 아래 그림 5.8.3와 같다. 이는종속면(slave segment) 또는 비-모르타 면(non mortar segment) 에서의 적분점을 주 종속면(master segment) 또는 모르타 면(mortar segment)에 직각 투영시켜 접촉면 에서의 적분을 수행한다. 적분점의 개수는 저차요소의 경우 4개 고차요소의 경우 5개를 사용하며, 이는 패치테스트(patch-test) 를 통과하기 위한 최소한의 적분점 개수이다.
text_image
Non-mortar segment :Integration point Mortar segment
그림 5.8.3 모르타 접촉의 개념
• 접촉면의 검색
접촉 검색은 종속절점(slave-node)/주접촉면(master segment) 알고리즘을사용한다. 모르타 접촉의 경우, 종속 절점은 종속 요소의 (slave segment)내부의 적분점(integration point)이다.
이 알고리즘은 종속절점과 주접촉면 간의 인접한 정도 또는 종속절점이
주접촉면을 관통한 정도에 따라 접촉 여부를 판정한다. 일반적으로 종속절점을 포함한 물체와 주접촉면을 포함한 물체는 서로 뒤바뀌어도 관계 없으나, 수치적인 관점에서 볼 때 상대적으로 강성이 큰 물체, 또는 상대적으로 요소가 조밀하지 않은 물체에 주접촉면을 정의해야 좀 더 정확한 해석 결과를 얻을 수 있다.
실제로 종속절점과 주접촉면의 접촉여부 판정을 하기 위해서는 전역 검색(global search) 과정을 거친다. 전역 검색이란 공간 상에서 물체 혹은 접촉면이 충돌 할 수 있는 예비 종속 절점/요소들을 결정하는 과정을 의미한다. 전역 검색에 의해 결정된 종속 절점/요소와 주접촉면 집합에 대해서만 접촉검색을 수행하게 된다.
실제로 절점또는 적분점이 주접촉면에 접촉하는 지를 파악하기 위해서는 그림 5.8.3과 같이 종속절점을 주접촉면 상에 직각 투영(orthogonal projection)해야 한다. 삼차원 접촉의 경우, 벡터 r 을 원점에서 투영된 점(A)까지의 벡터이라하고 x_{s} 를 원점에서 종속절점까지의 벡터라 하면 다음과 같은 식을 얻을 수 있다.
\frac {\partial \mathbf {r} \left(\xi_ {c} , \eta_ {c}\right)}{\partial \xi} \cdot \left[ \mathbf {x} _ {s} - \mathbf {r} \left(\xi_ {c}, \eta_ {c}\right) \right] = 0 \tag {5.8.1}
\frac {\partial \mathbf {r} \left(\xi_ {c} , \eta_ {c}\right)}{\partial \eta} \cdot \left[ \mathbf {x} _ {s} - \mathbf {r} \left(\xi_ {c}, \eta_ {c}\right) \right] = 0
여기서 (\xi_{c},\eta_{c}) 는 접촉점(A)의 위치를 주접촉면 상의 자연좌표계로 표현한 것이다. 위식을 만족하는 (\xi_{c},\eta_{c}) 는 뉴튼랩슨법을 이용하여 수치적으로 계산할 수 있다. 뉴튼랩슨법을 적용하기 위한 좌표 증분 (\Delta\xi_{c},\Delta\eta_{c}) 은 다음과 같다.
\left[ \begin{array}{l} \frac {\partial \mathbf {r}}{\partial \xi} \\ \frac {\partial \mathbf {r}}{\partial \eta} \end{array} \right] \left[ \begin{array}{l l} \frac {\partial \mathbf {r}}{\partial \xi} & \frac {\partial \mathbf {r}}{\partial \eta} \end{array} \right] \left\{ \begin{array}{l} \Delta \xi \\ \Delta \eta \end{array} \right\} = \left[ \begin{array}{l} \frac {\partial \mathbf {r}}{\partial \xi} \\ \frac {\partial \mathbf {r}}{\partial \eta} \end{array} \right] \left\{\mathbf {r} \left(\xi_ {c}, \eta_ {c}\right) - \mathbf {x} _ {s} \right\} \tag {5.8.2}
위 식은 초기조건 \left(\xi_{c},\eta_{c}\right)=(0,0) 을 이용한 경우, 접촉점(A) 또는 종속절점(B)의 위치가 주접촉면으로부터 멀리 떨어지지 않다면 쉽게 수렴한다. 이제 종속절점이 접촉 면을 관통했는지를 검사하고 관통한 경우에 대해서는 관통한
깊이에 비례하는 힘(접촉력)을 종속 절점과 접촉 면에 부가한다.
text_image
B n_A A η ξ
그림 5.8.4 종속절점과 주접촉면 사이의 수직 관계
이차원, 축대칭 문제의 경우, 해당되는 자연좌표계의 개수는 2개가 아닌 이며,계산 방식은 위와 동일하다.
• 접촉력의 계산
접촉으로 판단된 종속절점과 주접촉면 사이의 변위 관계는 벌칙기법(penaltymethod)을 통하여 구속된다. midas NFX에서는 절점 사이의 갭(gap)과 접촉력을다음 (5.8.3), (5.8.4)와 같이 정의한다.
\mathbf {g} _ {N} = \left(\mathbf {x} ^ {B} - \mathbf {x} ^ {A}\right) \cdot \mathbf {n} ^ {A} \tag {5.8.3}
f ^ {C} = - k _ {N} g _ {N} \quad \text { if } g _ {N} < 0 \tag {5.8.4}
k _ { n } : 벌칙 계수
\mathbf { x } ^ { A } , \mathbf { x } ^ { B } : 주접촉면 상의 점(A)와 종속절점(B)의 위치벡터
\mathbf { n } ^ { A } : 주접촉면 상의 점(A)의 수직벡터