J. Comput. Struct. Eng. Inst. Korea, 28(5) pp., October, 2015 http://dx.doi.org/10.7734/coseik.2015.28.5. pissn 1229-3059 eissn 2287-2302 Computational Structural Engineering Institute of Korea Co-rotational 이론 기반 비선형 삼각평면 유한요소의 개발 조 해 성 1 신 상 준 1 1 서울대학교 기계항공공학부 Development of Nonlinear Triangular Planar Element Based on Co-rotational Framework Hae-Seong Cho 1 and Sang-Joon Shin 1 1 Department of Aerospace and Mechanical Engineering, Seoul National University, Seoul, 08826, Korea Abstract This paper presents development of a geometrically nonlinear triangular planar element including rotational degrees of freedom, based on the co-rotational(cr) formulation. The CR formulation is one of the efficient geometrically nonlinear formulations and it is based on the assumptions on small strain and large rotation. In this paper, modified CR formulation is suggested for the developemnt of a triangular planar element. The present development is validated regarding the static and time transient problems. The present results are compared with the results predicted by the previous researchers and those obtained by the existing commercial software. Keywords : geometrically nonlinear structural analysis, co-rotational formulation, triangular planar element 1. 서 론 최근까지 많은 선행 연구들을 통해 다양한 형태로 정밀한 구 조해석 모델이 개발되어 왔고, 상당부분 유한요소를 사용하고 있다. 그 중, 저차의 유한요소는 비교적 간단한 정식화를 포함 하고 효율적으로 구조물을 해석할 수 있다. 이러한 이유로 많 은 연구자들에 의해 정확하고 효율적인 4절점 사각평면요소 와 3절점 삼각평면 요소들이 개발되어 왔다(Allman, 1984; 1988; Bergan et al., 1985; Carpenter et al., 1985). 하 지만 대부분의 요소는 선형의 거동을 예측할 수 있다. 비선형 구조해석으로의 확장을 위해 대표적으로 total Lagrangian, updated Lagrangian 기법이 널리 사용되었다(Bathe et al., 1975). 이러한 전통적인 해석기법은 요소의 형태와 요소 절점의 자유도에 따라 비선형성을 고려하는데 추가의 수학적 모델링이 요구되며 요소 자체의 가정에 따라 상당한 어려움이 따를 수 있다(Simo et al., 1993). 비교적 최근 Rankin에 의해 co-rotational(cr) 기법이 정립되었으며 작은 변형률과 큰 회전거동을 갖는 구조물의 기 하학적 비선형해석에 적용 가능하다(Rankin et al., 1985). CR 기법의 가장 큰 장점은 요소에 독립적으로 적용가능 함이 다. 따라서 요소의 절점 수, 절점의 자유도에 따라 동일하게 기법을 적용하여 기존 요소를 비선형해석으로 확장 가능하다. 하지만, CR 기법은 대부분 보, 쉘 요소에 적용되었고 (Crisfield, 1990; Battini, 2007), 최근에 들어서 Battini (2008)에 의해 CR 기법을 적용한 비선형 평면요소의 정식화 가 수행되었다. 하지만 Battini에 의해 정립된 비선형 평면 요소는 평면 상 회전자유도가 고려되지 않았고 병진 거동에 대한 자유도만을 갖는 사각 평면요소를 고려하였다. 국부요소 로 비 적합 4절점 사각평면 요소(Taylor et al., 1976)를 사 용하여 잠김 현상에서 자유로운 비선형해석을 수행하였으나 궁극적으로 회전자유도를 포함하고 있지 않아 플래핑 날개와 같이 구조물의 회전 거동을 포함한 구조해석에 적용하기에는 분명한 한계가 있다. 회전자유도를 포함한 평면요소는 회전 거동 및 토오크를 직접 고려할 수 있을 뿐만 아니라, 잠김 현 Corresponding author: Tel: +82-2-880-1642; E-mail: ssjoon@snu.ac.kr Received June 16 2015; Revised July 21 2015; Accepted July 22 2015 C 2015 by Computational Structural Engineering Institute of Korea This is an Open-Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License(http://creativecommons. org/licenses/by-nc/3.0) which permits unrestricted non-commercial use, distribution, and reproduction in any medium, provided the original work is properly cited. 한국전산구조공학회 논문집 제28권 제5호(2015.10) 11
Co-rotational 이론 기반 비선형 삼각평면 유한요소의 개발 상에서 자유로울 수 있어 평면거동을 하는 구조물에 효과적인 해석을 수행할 수 있다. 따라서 본 논문에서 회전자유도를 갖는 평면요소에 대한 CR 정식화를 수행하였으며 이를 동적해석으로 확장하였다. 개 발한 요소는 선행 연구에서 개발된 평면요소들과 비교 검증하 였고 상용프로그램인 ANSYS 해석결과와 비교 검증하였다. 다. 각 좌표계는 요소 형상의 무게 중심을 기준으로 고려되었 고, 무게중심에 해당하는 물리량은 하첨자 로 나타내었다. 또 한 하첨자 와 은 각각 전역좌표계와 국부좌표계에서 정의된 물리량을 나타낸다. 그리고,,, 는 각각 절점 좌표와 변위를 나타내며, 는 무게중심 기준으로 나타낸 상대적 절 점 좌표를 나타낸다. 본 논문에서 좌표계는 무게중심에 놓여지고, 국부 요소의 2. 본 론 회전자유도를 고려하기 위해 CR 좌표계를 설정하였다. 전역 좌표계에서 CR 좌표계로의 변환에 따른 좌표계의 회전은 식 이 장에서는 CR 기법과 이를 적용한 회전자유도를 갖는 평면 요소의 비선형해석 정식화 과정을 소개하고자 한다. (1)과 같이 도출할 수 있다(Battini, 2008). 2.1 CR 비선형 구조해석 기법 (1) 이때, 는 회전변환행렬 을 구성하고, 변형 전 좌표계 에서 정의되는 는 회전변환행렬 을 구성한다. 회전변 환행렬을 사용하여 각각 식 (2), (3)과 같이 순수 변형의 병 진 자유도와 회전자유도를 나타낼 수 있다. Fig. 1 Geometry and coordinate of a planar element (2) (3) CR 기법은 국부요소를 독립적으로 사용하여 기존의 기하 학적으로 선형이며 정확한 예측이 가능한 요소를 기하학적 비선형해석으로 확장할 수 있다(Rankin et al., 1986, Felippa et al., 2005). Fig. 1은 CR 기법의 좌표계를 나 타낸다. CR 기법은 기학적 비선형적 거동을 예측하기 위해 전체 거동을 CR 좌표계를 기준으로 강체 거동과 순수 구조 변형으로 분리하는 두 단계로 나누어 고려한다. 이때 기존의 선형요소는 순수 구조변형을 나타내는 국부요소의 형태로 적 용되고, 회전변환행렬을 도입하여 강체 거동을 모사할 수 있 다. 따라서 기존 선형 요소의 강성행렬과 내력 벡터에 회전 변환행렬을 도입하여 전역좌표계에서의 강성행렬과 내력벡터 를 구성할 수 있다. 최종적으로 회전변환행렬, 은 요소에 대한 회전변환행 렬 로 식 (4)와 같이 정의할 수 있고, 이를 도입하여 전 역좌표계에서 정의된 강성행렬 과 내력 벡터 는 식 (5)와 같다. (4) where (5) 2.2 회전자유도를 갖는 CR 삼각평면요소 Fig. 1에서 나타낸 는 전역좌표계이고,, 는 각각 전역좌표계에서의 위치와 변위를 나타낸다. 하 첨자 는 요소 절점을 나타내고 은 전체 요소의 절점 개수이 행렬 는 식 (6)과 같이 전역좌표계에서 정의된 변위를 국부 좌표계에 대한 미분으로 유도되며 전체 요소 거동에서 순수 변형에 대한 물리량을 추출한다. 이때, 벡터 와 는 절점 병진 자유도에 의해 유도되며, 최종적으로 행렬 는 식 (7)과 같이 구성할 수 있다(Battini, 2007). 12 한국전산구조공학회 논문집 제28권 제5호(2015.10)
조해성 신상준 (6) (7) 국부요소는 OPtimal Triangular(OPT) 평면 요소(Felippa, 2003)를 적용하였다. 각 절점에 회전자유도를 포함하는 OPT 평면요소는 기존의 유한요소 정식화와 달리 자유 정식화(free formulation) 과정을 바탕으로 하여, 요소의 세장비에 관계없이 면 내, 굽힘 거동을 정확하게 예측 가능한 선형 요소이다. 또한 OPT 평면요소를 정의하는데 있어 평면 응력(plane stress) 연 속 방정식을 적용하였다. 2.3 관성력에 의한 가상일 국부좌표계에서 정의된 전단 관성 행렬은 국부 자유도에 대한 관성력의 미분으로 나타낼 수 있다. 따라서 관성력 벡 터와 동일한 방법으로 요소 자유도를 도출하여 적용하면 아 래와 같이 관성행렬, 감쇠행렬을 표현할 수 있다. (12) 이렇게 정의된 각각의 행렬과 벡터는 회전변환 행렬 를 통해 전역좌표계에서 정의된 물리량으로 표현 가능하며 요소의 운동에 대한 비선형 지배방정식은 식 (14)와 같다. Fig. 1과 같이 정의된 요소의 좌표계에서 운동에너지는 아래와 같이 나타낼 수 있다. (8) 식 (8)에서,, 는 각각 밀도, 면적, 두께를 나타낸 다. 운동에너지와 관성력에 의한 가상일의 관계에 따라 관성 력 벡터는 다음과 같이 나타낼 수 있다. (9) 위 는 요소 절점 자유도의 변분 항이며 병진 자유도와 회 전자유도로 구성된다. Fig. 1에서 정의한 좌표계에 따라 요 소 절점 자유도는 국부 좌표계에서 정의된 물리량으로 고려 할 수 있다. 이 때, 요소의 절점 자유도는 미소 물리량이므로 요소의 형상함수를 통해 요소 변위로 나타낼 수 있다. (10) 최종적으로 요소의 관성력 벡터는 아래와 같이 나타낼 수 있다. (11) (13) (14) 본 논문에서는 지배방정식을 수치적으로 해석하기 위해 Newton-Raphson 기법과 Hilber-Hughes-Taylor-α(HHTα) 내재적 시간적분 기법(Crisfield, 1997)을 적용하였다. 3. 수치해석결과 본 장에서는 앞서 개발된 요소의 정적/동적 해석을 수행하 고 선행 연구결과 및 상용프로그램인 ANSYS의 해석결과와 비교 검증하였다. 3.1 정적해석 검증 정적해석 검증을 위해 Battini(2008)에 의해 수행되었던 예제를 고려하였다. 해석에 적용된 형상 및 물성은 Fig. 2와 같다. 본 예제에서는 비 정렬 격자를 사용함에 따라 도출되 는 해의 정확성을 질 좋은 격자를 사용하여 도출된 해와 비 교하였다. 비 정렬 격자는 절점의 위치를 정렬된 위치 기준 20% 뒤틀린 형상을 고려하였다. 점 A 위치에 축방향의 집 중하중이 가해졌고 하중의 크기에 따라 발생하는 수평방향 변위를 선행연구에서 제시된 결과들과 비교 검증하였다. 점 A 위치에 가해진 하중의 크기는 40000이고, 하중인자()를 변경함에 따라 발생하는 변위를 고려하였다. 해석결과는 Fig. 3에 나타내었으며 하중 인자가 1일 때의 형상은 Fig. 4와 같다. 해석결과를 살펴보면 기하학적 비선형적 거동특성 한국전산구조공학회 논문집 제28권 제5호(2015.10) 13
Co-rotational 이론 기반 비선형 삼각평면 유한요소의 개발 Table 1 Comparison of the horizontal displacement at Point A(λ=1.0) Reference Displacement Difference(%) Present 6.71 0.95 Qm6(Battini, 2008) 6.39 5.51 Qnew(Battini, 2008) 6.49 4.03 Fine grid 6.77 - Fig. 5 Cantilevered flat plate loaded by the surface pressure Fig. 2 Geometry and properties of an angle frame Fig. 6 Comparison of the tip deflection Fig. 3 Comparison of the horizontal displacement at Point A 이어서, Bathe 등(1975)에 의해 수행된 선행연구 사례를 고려하였다. 한 쪽 면이 지지된 평판에 윗면과 아랫면에 압 력이 가해지는 조건이 적용되었다. 해석에 사용된 구조물의 형상과 물성은 Fig. 5와 같다. 해석결과, 본 연구를 통해 도 출된 해석결과와 선행 연구결과가 잘 일치하는 것을 확인하 였다. 해석결과는 Fig. 6에 나타내었다. 3.2 동적해석 검증 Fig. 4 Deformed configuration by the present analysis (λ=1.0) 이 잘 나타남을 확인할 수 있으며, 질 좋은 격자를 사용하여 도출된 결과와 비교하였을 때, 상대오차는 0.9%이다. 또한 기존의 선행 결과와 비교하였을 때, 본 연구를 통해 개발된 평면요소가 비 정렬 격자에서 더 좋은 성능을 나타냄을 확인 하였다. 하중 인자가 1일 때, 각각의 결과 비교를 Table 1 에 정리하였다. 동적해석 검증을 위해 앞서 정적해석의 검증 예제로 수행 되었던 Bathe 등(1975)의 예제를 고려하였다(Fig. 5). 시 간에 따라 일정한 압력이 가해짐에 따라 발생하는 끝단의 변 위 응답을 비교하였다. 본 연구에서 개발된 평면 요소는 총 20개가 사용되어 총 66개의 절점 자유도를 갖는 반면, 선행 연구에서는 8절점 사각 평면 요소 10개가 적용되어 총 84개 의 절점 자유도가 적용되었다. 상대적으로 적은 수의 자유도 가 적용되었음에도 불구하고 해석결과, 선행 연구결과와 잘 일치하는 것을 확인하였다. 또한, 고차 요소를 사용한 결과와 비교 검증을 통해 본 요 소의 효율성과 도출되는 해의 정확성을 확인하기 위해 3차원 고체 요소를 적용한 ANSYS 해석과의 비교를 시도하였다. 14 한국전산구조공학회 논문집 제28권 제5호(2015.10)
조해성 신상준 Fig. 7 Transient response of the cantilevered flat plate under constant surface pressure 적용하여 지배방정식을 해석하였다. 수치해석은 정적해석과 동적해석의 몇 가지 예제를 사용하였고 해석결과는 선행 연 구결과 및 상용프로그램 해석결과와 비교 검증하였다. 검증 결과, 정적해석에서의 처짐은 0.9%의 오차를 나타내었으며, 동 적해석결과, 변위 응답이 3차원 고체 요소를 사용한 ANSYS 해석결과와 잘 일치하는 것을 확인하였다. 향후, 본 논문을 통해 개발한 평면 요소를 플래핑 날개의 구조해석에 적용하여 유체-구조 결합해석을 수행하고 나아가 국부 Lagrangian 승 수를 적용한 다물체 동역학 해석으로 확장을 수행할 것이다. 감사의 글 Fig. 8 Cantilevered plate under the harmonic tip load This work was supported by Advanced Research Center Program(No. 2013073861) through the National Research Foundation of Korea(NRF) grant funded by the Korea government(msip) contracted through Next Generation Space Propulsion Research Center at Seoul National University and also be a grant to Bio-Mimetic Robot Research Center funded by Defense Acquisition Program Administration. Fig. 9 Transient response of the cantilevered plate under the harmonic tip load 해석에 적용된 형상과 물성은 Fig. 8과 같고 시간에 따라 변 화하는 끝단 집중하중이 가해짐에 따라 발생하는 끝단의 변 형 응답을 비교 검증하였다. 가진 주파수는 10Hz로 선정하 였다. 본 논문에서 개발한 평면요소는 총 20개의 요소를 사용 하여 66개의 총 절점 자유도가 고려되었으나, 3차원 ANSYS 해석에는 총 80개의 요소를 사용하여 총 594개의 절점 자유 도가 고려되었다. 비교 결과, 현 해석을 통해 도출된 끝단 응 답결과가 ANSYS에 비하여 매우 적은 개수의 요소를 적용했 음에도 불구하고 ANSYS 해석결과와 잘 일치하는 것을 확인 하였다. 4. 결 론 본 논문에서는 CR 기법을 적용하여 회전자유도를 갖는 3 절점 비선형 삼각평면요소를 개발하였다. 국부 요소의 회전 을 고려하기 위해 요소 형상의 무게중심에 CR 좌표계를 정 의하였다. 또한 동적해석을 위해 운동에너지의 변분으로부터 전단 질량행렬, 감쇠행렬 그리고 관성력 벡터를 유도하였고 Newton-Raphson 기법과 HHT-α 내재적 시간적분 기법을 References Allman, D.J. (1984) A Compatible Triangular Element Including Vertex Rotations for Plane elasticity Problems, Compos. Struct., 19, pp.1 8. Allman, D.J. (1988) A Quadrilateral Finite Element Including Vertex rotations for Plane Elasticity Problems, Int. J. Nummer. Meth. Eng., 18, pp.717 739. Bathe, K-J., Ramm, E., Wilson, E.L. (1975) Finite Element Formulations for Large Deformation Dynamic Analysis, Int. J. Nummer. Meth. Eng., 9, pp.353~ 386. Battini, J-M. (2007) A Modified Corotational Framework for Triangular Shell Elements, Comput. Meth. Appl. Mech. Eng., 196, pp.1905~1914. Battini, J-M. (2008) A Non-linear Corotational 4-node Plane Element, Mech. Res. Commun., 35, pp.408~413. Bergan, P.G., Fellipa, C.A. (1985) A Triangular Membrane Element with Rotational Degrees of Freedom, Comput. Meth. Appl. Mech. Eng., 50, pp.25~60. Carpenter, N., Stolarski, H., Belytschko, T. 한국전산구조공학회 논문집 제28권 제5호(2015.10) 15
Co-rotational 이론 기반 비선형 삼각평면 유한요소의 개발 (1985) A Flat Triangular Shell Element with Improved Membrane Interpolation, Comput. Meth. Appl. Mech. Eng., 1, pp.161~168. Crisfield, M.A. (1990) A Consistent Co-rotational Formulation for Non-linear, Three-dimenstional, Beam-elements, Comput. Meth. Appl. Mech. Eng., 81, pp.2969~2992. Crisfield, M.A. (1997) Non-linear Finite Element Analysis of Solid and Structures, Advanced Topics, 2, Wiley, London. Felippa, C.A. (2003) A Study of Optimal Membrane Triangles with Drilling Freedoms, Comput. Meth. Appl. Mech. Eng., 192, pp.2125~2168. Felippa, C.A., Haugen, B. (2005) A Unified Formulation of Small-strain Coroational Finite Elements: I. Theory, Comput. Meth. Appl. Mech. Eng., 194, pp.2285~2335. Pai, P.F. (2007) Highly Flexible Structures: Modeling, Computation, and Experimentation, AIAA, Inc. VA. Rankin, C.C., Nour-Omid, B. (1986) An Element Independent Co-rotational Procedure for the Treatment of Large Rotations, ASME J. Pressure Vessel Techn., 108, pp.165~174. Simo, J., Armero, F. (1993) Improved Version of Assumed Enhanced Strain Tri-linear Element for 3d Finite Deformation Problems, Comput. Meth. Appl. Mech. Eng., 110, pp.359~386. Taylor, R., Beresford, P., Wilson, E. (1976) A Non-conforming Element for Stress Analysis, Int. J. Nummer. Meth. Eng., 10, pp.1211~1219. 요 지 구조의 기하학적 비선형해석을 위해 대표적으로 Total Lagrangian, Updated Lagrangian 정식화 기법이 있다. 이러한 고 전적인 정식화 과정은 요소의 변형률을 가정하는 방법에 따라 그리고 요소의 절점 수에 따라 추가의 수학적 정식화 과정이 요구된다. 하지만 비교적 최근에 정립된 Co-roational(CR) 이론은 기 존재하는 보, 판, 쉘 요소에 독립적으로, 요소 절점 자 유도에 따라 일정하게 적용하여 대변위, 작은 변형률을 갖는 구조의 기하비선형 해석을 가능케 한다. 본 논문에서는 회전 자 유도를 갖는 삼각평면요소에 대한 CR 기법을 정식화하였고 동적해석으로 확장하여 이를 상용프로그램과 검증하였다. 해석 에 사용한 삼각평면요소는 OPtimal Triangular(OPT) 평면요소이다. 핵심용어 : 기하학적 비선형 구조해석, Co-rotational 이론, 삼각평면요소 16 한국전산구조공학회 논문집 제28권 제5호(2015.10)