Journal Search Engine
Search Advanced Search Adode Reader(link)
Download PDF Export Citaion korean bibliography PMC previewer
ISSN : 1229-3059(Print)
ISSN : 2287-2302(Online)
Journal of the Computational Structural Engineering Institute of Korea
Vol.30 No.1 pp.1-6

DOI : https://doi.org/10.7734/COSEIK.2017.30.1.1

Development of Nine-node Co-rotational Planar Element for Elastoplastic/Contact Analysis

Hae-Seong Cho1, Hyun-Shig Joo1, Sang Joon Shin1
1Department of Aerospace and Mechanical Engineering, Seoul National University, Seoul, 08826, Korea
Corresponding author : +82-2-880-1642; E-mail: ssjoon@snu.ac.kr
September 19, 2016 December 1, 2016 December 19, 2016

Abstract

This paper presents development of the nine-node co-rotational(CR) planar element applicable for elastoplastic and contact analysis. The CR formulation is one of the efficient geometrically nonlinear formulations. It is based on the assumptions of small strain and large displacement. Further, it is extended to both elastoplastic analysis and contact analysis in this paper. For accurate plastic analysis, nine-node quadrilateral element, which can provide accurate stress prediction, is chosen. Bi-linear hardening rule based on Newton- Raphson return-mapping is employed. Also, Lagrange multiplier is used in order for constraints regarding the contact analysis. The present development is validated via the time transient problems. The present results are compared with those obtained by the other existing software.


탄소성/접촉 해석을 위한 Co-rotational 정식화 기반의 9절점 평면 요소 개발

조 해 성1, 주 현 식1, 신 상 준1
1서울대학교 기계항공공학부

초록

본 논문에서는 비교적 최근 정립된 co-rotational 이론을 기반으로 한 4절점 평면요소 정식화를 확장하여 9절점 평면 요소에 적합한 CR 정식화를 제시하고자 한다. 그리고 등방성 재료의 소성 해석을 위해, 선형 경과 규칙(bi-linear hardening rule)을 바탕으로 하는 Newton-Raphson return-mapping 알고리즘을 적용하였다. 이때, von Mises 기준을 적용하여 소성 변형 상태를 예측하였다. Lagrange 승수를 도입하여 2차원 접촉에 대한 구속조건을 부여하였다. 개발한 요소는 상용프로그램인 ABAQUS 해석결과와 비교 검증하였다.


    HighSpeed

    1.서 론

    다양한 분야에서 구조물의 기하학적 비선형성 및 재료적 비선형성을 포함한 정밀구조해석의 요구는 점점 높아지고 있다. 특히 고온, 고압의 환경에 노출되어 있는 구조물은 기하비선형 및 재료 비선형성이 크게 발생할 수 있다. 항공우주 분야에서는 높은 압력 및 온도에 의해 다양한 물성을 가진 구조물에 대한 기하학적 탄성 및 소성 변형이 일어난다. 예를 들어 고속비행체 발사 시, 추진기관으로부터 발생하는 고온, 고압의 동적 하중은 수직 발사장치 내부 구조물들과의 유동-구조 연계현상을 일으켜 초기 발사 안정성에 큰 영향을 미치게 된다. 다중발사관의 형상을 가진 수직 발사장치의 경우, 인접한 발사관 내 고속비 행체의 추진화염으로부터 내부 고속비행체를 보호하기 위해 후방덮개 구조물을 갖는다. 이 후방덮개 구조물은 길이에 비해 얇은 두께를 가진 단열재/금속재로 구성되며, 특정한 각도 및 위치를 가진 추가 구조물에 부딪혀 추진기관 화염의 방향을 편향하게 해주는 장치이다. 추진기관의 화염에서 발생하는 동적 하중들에 의해 길이에 비해 얇은 두께를 갖고 금속 및 복합재로 구성된 후방덮개는 기하학적 비선형 및 소성영역에서의 거동을 나타낸다.

    비선형 구조 해석으로의 확장을 위해 대표적으로 total Lagrangian, updated Lagrangian 기법이 널리 사용되었다 (Bathe et al., 1975). 이러한 전통적인 해석 기법은 요소 자 체 가정에 따라 비선형 정식화 가정에 따라 상당한 어려움이 따를 수 있다(Simo et al., 1993). 즉, 요소에 따라 기하학적 비선형성을 고려하기 위한 추가의 수학적 모델링이 요구된다.

    최근 기하학적 비선형 해석을 위해 co-rotational(CR) 기법이 많이 사용되고 있으며 작은 변형률과 대 변위를 갖는 구조물의 해석에 적합하다(Rankin et al., 1985). CR 기법은 요소에 독립적으로 적용 가능하여 요소 절점 수 및 절점 자유도에 따라 동일하게 적용하여 기존의 선형 요소를 기하비선형 해석으로 확장시킬 수 있다. 평면 요소에 대한 CR 기법 정식화는 Battini(2008)에 의해 처음 제시되었다. 본 논문의 저자는 면내 회전자유도를 갖는 평면 요소에 적용 가능한 CR 정식화를 제시하였다(Cho et al., 2015; 2016).

    본 논문에서는 Battini에 의해 제시된 정식화를 확장하여 면내 잠김 현상에 자유롭고 응력 예측 정확성이 비교적 높은 9 절점 사각 평면 요소에 적합한 CR 정식화를 제시하고자 한다. 그리고 등방성 재료의 소성 해석을 위해, 선형 경과 규칙 (bi-linear hardening rule)을 바탕으로 하는 Newton-Raphson return-Mapping 알고리즘을 적용하였다. 이때, von Mises 기준을 적용하여 소성 변형 상태를 예측하였다. Lagrange 승수를 도입하여 2차원 접촉에 대한 구속조건을 부여 하였다. 개발한 요소는 상용프로그램인 ABAQUS 해석결과 와 비교 검증하였다.

    2.본 론

    2.1.CR 이론 기반 9절점 평면 요소

    Fig. 1은 Co-rotational(CR) 기법의 좌표계를 나타낸다. CR 기법은 기하학적 비선형적 거동을 예측하기 위해 전체 거동을 CR 좌표계를 기준으로 강체 거동과 순수 구조 변형으로 분리하는 두 단계로 나누어 고려한다. 순수 구조 변형과 강체 회전거동을 분리하여 고려함에 따라 CR 좌표계가 도입되며 기존의 존재하는 선형의 요소를 적용하면 탄성 역역에서의 기하학적 비선형성을 예측할 수 있는 구조해석을 도모할 수 있다. 본 논문에서 좌표계는 무게중심에 놓여지고, 국부 좌표계 내 변형량을 고려하기 위해 CR 좌표계를 설정하였다. 전역좌표계 에서 CR 좌표계로의 변환에 따른 좌표계의 회전은 식 (1)과 같이 도출할 수 있다(Battini, 2008).

    tan  θ = i = 1 9 [ x i ( Y i + V i Y c V c ) y i ( X i + U i X c U c ) ] i = 1 9 [ y i ( Y i + V i Y c V c ) + x i ( X i + U i X c U c ) ]
    (1)

    이때, θ 는 요소 강체 회전 변환 행렬 [R ] 을 구성한다. 회전 변환 행렬을 사용하여 각각 식 (2)와 같이 순수 변형의 병진 자유도를 나타낼 수 있다.

    { u i ¯ υ i ¯ } = [ R ] T { X i + U i X c U c Y i + V i Y c V c } { x i y i }
    (2)

    최종적으로 회전 변환 행렬, [R ] 은 요소에 대한 회전 변환 행렬 [B ] 로 식 (3)과 같이 정의할 수 있고 이를 도입하여 전역좌표계에서 정의된 강성행렬 [KG ] 과 내력 벡터 {fG}는 식 (4)와 같다.

    [ B ] = [ P ] [ d i a g ( [ R ] ) ] T [ K g ] = [ B ] T [ K l ] [ B ] + [ K h ] ,    [ f g ] = [ B ] T { f l }
    (3)

    where,

    where, [ K h ] = δ [ d i a g ( [ R i ] ) ] [ P ] T { f l } + [ d i a g ( [ R i ] ) ] δ [ P ] T { f l }
    (4)

    행렬 [P ] 는 전체 요소 거동에서 순수 변형에 대한 물리량을 투영하는 행렬을 나타내며 최종적으로 식 (5)와 같이 구성할 수 있다.(6)

    [ P i ] = I 2 { A i } { G i } T
    (5)

    { A i } = { { r l } i } T , { G i } = 1 k = 1 9 { r } i T { r G } i { { r } i } T
    (6)

    이때 { r l } i = [ R ] { r } i = [ R ] { X i + U i X c U c Y i + V i Y c V c } . .

    2.2동적해석 정식화

    9절점 요소의 관성행렬은 요소의 형상함수를 이용하여 구할 수 있다. 이렇게 정의된 관성행렬과 요소의 감쇠 행렬은 앞서 정의한 요소의 강체 거동에 따른 회전변환 행렬 [ diag ([R ]) ] 를 통해 식 (8)과 같이 전역좌표계에서 정의된 물리량으로 표현 가능하다. 최종 요소의 운동에 대한 비선형 지배방정식은 식 (9)와 같다.(7)

    [ M G ] = [ d i a g ( [ R ] ) ] T [ M L ] [ d i a g ( [ R ] ) ]
    (7)

    M G ( s ) x ¨ ( t ) + K G ( s ) x ( t ) = F ( s , t )
    (8)

    본 논문에서는 비선형의 지배방정식을 수치적으로 해석하기 위해 Newton-Raphson 반복 해법과 Hilber-Hughes-Taylor -α(HHT-α) 내재적 시간적분 기법을 적용하였다(Crisfield, 1997).

    2.3.Newton-Raphson return-Mapping 기법 - 재료 비선형

    CR 정식화는 대변위, 작은 변형률의 가정을 바탕으로 한다. 그리고 정식화에 도입되는 요소의 변형은 국부좌표계에서 계산된다. 하지만 소성 변형의 경우, 일반적인 탄성 기하비선형 해석에 비해 국부 변형의 크기가 커 그에 따른 적절한 좌표계 도입이 요구된다. 따라서 본 논문에서는 Lagragian 좌표를 도입하여 소성 변형에 따라 발생하는 변형률과 응력을 고려 하고자 한다. 소성을 고려한 국부 요소의 강성행렬은 최종적으로 아래의 식과 같다(Bathe, 1996).

    [ K l ] = [ K C ] + [ K T ] = Ω [ B ] [ F ] T [ C p l ] [ F ] [ B ] T d Ω + Ω [ B ] [ S ] [ B ] T d Ω
    (9)

    이때, 행렬 [B ] , [F ] 그리고 [S ] 는 변형률-변위 관계 행렬, 변형구배행렬 그리고 소성변화에 따라 규정되는 응력행렬을 나타낸다.

    재료의 소성해석을 위해 일정한 기울기를 갖는 선형 경화 규칙을 따르는 Newton-Raphson return-mapping 알고 리즘을 적용하였다(Neto et al., 2008). von-Mises 응력을 판단 기준으로 도입, 계산된 내부 응력과 비교함으로써 소성 영역에 대해 해석한다. 구조물의 변형 후 상태를 도출하기 위해 변형률 증분은 아래와 같이 구할 수 있다.(11)(12)

    Δ = i + 1 i
    (10)

    그리고 Fig. 2와 같이 응력-변형률 관계에서 정의되는 탄성 변형률과 소성 변형률은 변형률 증분을 사용하여 계산할 수 있다.

    n + 1 e t r i a l = n e + Δ
    (11)

    n + 1 p t r i a l = n p
    (12)

    이렇게 정의된 변형률 성분들을 사용하여 항복함수에 의거한 응력을 도출할 수 있고 von Mises 응력을 판단 기준에 따라 소성변형/탄성변형 상태를 판단할 수 있다.(2)(13)

    Φ = S σ y ( n + 1 p t r i a l )
    (13)

    구조물이 항복점을 지나 소성상태가 될 경우, 구조물의 경화 흐름을 만족하는 변형률과 응력을 반복계산을 통해 도출할 수 있다.(14)(15)(16)

    n + 1 p = n + 1 p t r i a l + γ
    (14)

    σ n + 1 = A ( γ ) σ d n + 1 t r i a l
    (15)

    n + 1 e = [ C p l ] 1 σ n + 1
    (16)

    2.4.2차원 접촉 해석

    본 논문에서는 마찰에 의한 효과는 무시하고 강체와 탄성체의 접촉에 의한 효과를 고려하기 위해 Lagrange 승수를 도입하여 탄성체에 미치는 반력을 계산하고자 하였다(Zienkiewicz, et al., 2005). Fig. 3은 본 연구에서 도입한 접촉 해석 기법의 개념도를 나타낸다.

    A는 강체로 가정된 경계면을 나타내고, B는 탄성체를 나타낸다. 이때 탄성체의 거동에 따라 강체의 경계면을 넘어 가는 절점에 대해 Lagrange 승수를 활용하여 탄성체의 절점과 강체경계면 사이에 구속조건을 부여한다.

    탄성체의 절점이 강체 경계면을 통과한 경우, 경계면에 수직한 방향으로 통과한 변형량을 lp로 정의하면 아래와 같이 통과한 변형량에 대한 구속조건을 정의할 수 있다.

    l p = 0
    (17)

    따라서 탄성체와 강체 경계면에 접촉이 일어났을 경우, 전체 계의 지배 방정식은 강체 경계면에 대한 penalty 항(p)을 포함하여 아래와 같이 정립할 수 있다.

    [ [ K g ] 0 [ C 2 ] T 0 p [ I ] [ C 1 ] T [ C 2 ] [ C 1 ] 0 ] { { u } 0 { λ } } = { { f } p l p 0 }
    (18)

    이 때, 식 (18)의 지배방정식은 식 (17)의 구속조건에 의거한 접촉 조건을 만족할 경우 규정된다. 그리고 비선형 시간응답해석 알고리듬에서 접촉 조건을 만족한 경우, Fig. 4의 corrector-contact 단계에서 계산된다.

    3.수치해석결과

    본 장에서는 기 개발한 9절점 평면 요소의 동적 해석을 수행 하고 상용프로그램인 ABAQUS의 해석결과와 비교 검증하였다. 해석에 사용된 재료의 물성은 아래의 Table 1과 같다.

    3.1.외팔보 동적 응답 해석

    첫 번째 검증 예제로, 동적 끝단 하중을 받는 외팔보의 탄성과 소성해석을 수행하였다(Fig. 5). 본 논문에서 개발한 평면 요소는 총 50개의 요소를 사용하여 606개의 총 절점 자유도가 고려되었으나, ABAQUS 해석에는 총 500개의 요소를 사용 하여 총 3,422개의 절점 자유도가 고려되었다. 비교 결과, 현 해석을 통해 도출된 끝단 응답결과가 ABAQUS 해석 결과와 잘 일치하는 것을 확인하였다. 또한 소성 변형으로 인해 발생 하는 거동적 특징이 잘 예측되는 것을 확인하였다. 해석 결과는 Fig. 6에 나타내었다.

    3.2.링 구조물-강체 접촉 해석

    두 번째 검증 예제로, 집중하중을 받는 링 구조물을 고려 하였다. 링 구조물은 구속조건이 없고 아래로 향하는 집중하중에 의해 아래로 이동하다가 강체의 벽면에 접촉이 발생하는 예제를 Fig. 7과 같이 설정하였다. 본 논문에서 개발한 평면요소는 총 120개의 요소를 사용하여 720개의 총 절점 자유도가 고려 되었으나, ABAQUS 해석에는 총 540개의 요소를 사용하여 총 3,960개의 절점 자유도가 고려되었다. 검증을 위해 집중 하중을 받는 절점의 변위를 비교하였고 Fig. 8에 나타내었다, 현 해석을 통해 도출된 끝단 응답결과가 ABAQUS 해석 결과와 잘 일치하는 것을 확인하였다. 또한 본 논문에서 제시한 소성 해석이 탄성해석과 구분되는 거동 특성이 잘 예측되는 것을 확인하였다(Fig. 9).

    3.3.외팔보-강체 접촉 해석

    마지막 검증 예제로, 동적 끝단 하중을 받는 외팔보의 탄성과 소성해석을 수행하였다. 이때 외발보는 기울어진 강체 벽면과 보의 바로 위에 강체 벽면이 위치하도록 설정하였다(Fig. 10). 본 논문에서 개발한 평면요소는 총 50개의 요소를 사용하여 606개의 총 절점 자유도가 고려되었으나, ABAQUS 해석에는 총 500개의 요소를 사용하여 총 3,422개의 절점 자유도가 고려되었다. 비교 결과, 현 해석을 통해 도출된 끝단 응답결과가 ABAQUS 해석 결과와 비교적 잘 일치하는 것을 확인하였으나 시간이 지남에 따라 오차가 발생하는 것을 확인하였다. 이는 ABAQUS에서 시간응답 해석 시, 기법간의 차이로 인해 발생 한 것으로 사료된다. ABAQUS의 경우, 기 개발 요소의 해석 에서 적용한 시간간격에서 수렴하지 않는 문제가 발생하여 auto-time-stepping 기법을 적용하였다. 본 예제에서 역시 기 개발 정식화를 적용함에 따라 소성 변형으로 인해 발생하는 거동적 특징이 잘 예측되는 것을 확인하였다. 해석 결과는 Fig. 11에 나타내었다.

    4.결 론

    본 논문에서는 CR기법을 적용하여 9절점 비선형 사각 평면 요소를 개발하였다. 또한 동적해석을 위해 Newton-Raphson 기법과 HHT-α 내재적 시간적분 기법을 적용하여 지배 방정식을 해석하였다. 소성 해석을 위해 CR 국부 좌표계 내, Lagrangian 좌표를 도입하여 국부에서 발생하는 변형을 정밀하게 예측 하고자 하였다. 소성영역에서의 응력해석을 위해 일정한 기울기를 갖는 선형 경화 규칙을 따르는 Newton-Raphson return- Mapping 알고리즘을 적용하였고 Lagrange 승수를 도입한 접촉해석으로 추가 확장하였다. 수치해석은 동적해석의 몇 가지 예제를 사용하였고 해석결과는 상용프로그램 해석결과와 비교 검증하였다. 검증 결과, 변위 응답이 동일조건의 조밀한 구조 격자를 사용한 ABAQUS 해석결과와 잘 일치하는 것을 확인 하였다. 향후, 본 논문을 통해 개발한 평면 요소를 미사일 발사 장치 후방 덮개의 구조해석에 적용하여 유체-구조 결합해석을 수행할 것이다.

    감사의 글

    This research was supported by a grant to High- Speed Vehicle Research Center of KAIST with the support of Defense Acquisition Program Administration (DAPA) and Agency for Defense Development(ADD).

    Figure

    COSEIK-30-1_F1.gif

    Geometry and coordinate of a planar element

    COSEIK-30-1_F2.gif

    Strain-Stress curve

    COSEIK-30-1_F3.gif

    Conceptual schematic of two-dimensional contact analysis

    COSEIK-30-1_F4.gif

    Computational algorithm for contact analysis

    COSEIK-30-1_F5.gif

    Analysis condition of a cantilevered beam under harmonic tip load

    COSEIK-30-1_F6.gif

    Comparison of tip displacement history in Example 1

    COSEIK-30-1_F7.gif

    Analysis condition of a ring structure under contact with rigid wall boundary

    COSEIK-30-1_F8.gif

    Comparison of tip displacement history in Example 2

    COSEIK-30-1_F9.gif

    Deformed configurations

    COSEIK-30-1_F10.gif

    Analysis condition of a cantilevered beam under two rigid wall boundaries

    COSEIK-30-1_F11.gif

    Comparison of tip displacement history in Example 3

    Table

    Material properties

    Reference

    1. Bathe K-J , Ramm E , Wilson EL (1975) Finite Element Formulations for Large Deformation Dynamic Analysis , Int. J. Nummer. Meth. Eng, Vol.9 ; pp.353-386
    2. Bathe K-J (1996) Finite Element Procedures, Englewood Cliffs, Trentice Hall, Inc,
    3. Battini J-M (2008) A Non-linear Corotational 4-node Plane Element , Mech. Res. Communications, Vol.35 ; pp.408-413
    4. Cho H , Shin SJ (2015) Triangular Planar Element Based on Co-rotational Framework , J. Comput. Struct. Eng. Inst. Korea, Vol.28 ; pp.485-491
    5. Cho H , Kwak JY , Shin SJ , Lee N , Lee S (2016) Flapping Wing Fluid-Structureal Interaction Analysis using Co-rotatioonal Triangular Planar Structural Element , AIAA J, Vol.54 ; pp.2265-2276
    6. Crisfield MA (1997) Non-linear Finite Element Analysis of Solid and Structures, Advanced Topics, Wiley, Vol.2
    7. Rankin CC , Nour-Omid B (1986) An Element Independent Co-rotational Procedure for the Treatment of Large Rotations , ASME J. Pressure Vessel Tech, Vol.108 ; pp.165-174
    8. Simo J , Armero F (1993) Improved Version of Assumed Enhanced Strain Tri-linear Element for 3d Finite Deformation Problems , Comput. Meth. Appl. Mech. Eng, Vol.110 ; pp.359-386
    9. Neto EAS , Peric D , Owen DRJ (2008) Computational Methods for Plasticity, John Wiley & Sons Ltd Publication,
    10. Zienkiewicz O , Taylor R , Zhu J (2005) The Finite Element Mehod: Solid Mechanics, Elsevier,