Journal of the Computational Structural Engineering Institute of Korea. 2015. 485-490
https://doi.org/10.7734/COSEIK.2015.28.5.485

ABSTRACT


MAIN

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)를 사용하여 잠김 현상에서 자유로운 비선형해석을 수행하였으나 궁극적으로 회전자유도를 포함하고 있지 않아 플래핑 날개와 같이 구조물의 회전 거동을 포함한 구조해석에 적용하기에는 분명한 한계가 있다. 회전자유도를 포함한 평면요소는 회전 거동 및 토오크를 직접 고려할 수 있을 뿐만 아니라, 잠김 현상에서 자유로울 수 있어 평면거동을 하는 구조물에 효과적인 해석을 수행할 수 있다.

따라서 본 논문에서 회전자유도를 갖는 평면요소에 대한 CR 정식화를 수행하였으며 이를 동적해석으로 확장하였다. 개발한 요소는 선행 연구에서 개발된 평면요소들과 비교 검증하였고 상용프로그램인 ANSYS 해석결과와 비교 검증하였다.

2. 본 론

이 장에서는 CR 기법과 이를 적용한 회전자유도를 갖는 평면 요소의 비선형해석 정식화 과정을 소개하고자 한다.

2.1 CR 비선형 구조해석 기법

https://cdn.apub.kr/journalsite/sites/jcoseik/2015-028-05/07TK102015280502/images/10.7734.28.5.485.F100.png
Figure 1

Geometry and coordinate of a planar element

CR 기법은 국부요소를 독립적으로 사용하여 기존의 기하학적으로 선형이며 정확한 예측이 가능한 요소를 기하학적 비선형해석으로 확장할 수 있다(Rankin et al., 1986, Felippa et al., 2005). Fig. 1은 CR 기법의 좌표계를 나타낸다. CR 기법은 기학적 비선형적 거동을 예측하기 위해 전체 거동을 CR 좌표계를 기준으로 강체 거동과 순수 구조 변형으로 분리하는 두 단계로 나누어 고려한다. 이때 기존의 선형요소는 순수 구조변형을 나타내는 국부요소의 형태로 적용되고, 회전변환행렬을 도입하여 강체 거동을 모사할 수 있다. 따라서 기존 선형요소의 강성행렬과 내력 벡터에 회전변환행렬을 도입하여 전역좌표계에서의 강성행렬과 내력벡터를 구성할 수 있다.

2.2 회전자유도를 갖는 CR 삼각평면요소

Fig. 1에서 나타낸 O는 전역좌표계이고, ri(X,Y), ui(U,V)는 각각 전역좌표계에서의 위치와 변위를 나타낸다. 하첨자 J는 요소 절점을 나타내고 N은 전체 요소의 절점 개수이다. 각 좌표계는 요소 형상의 무게 중심을 기준으로 고려되었고, 무게중심에 해당하는 물리량은 하첨자 c로 나타내었다. 또한 하첨자 gl은 각각 전역좌표계와 국부좌표계에서 정의된 물리량을 나타낸다. 그리고 X, Y, U, V는 각각 절점 좌표와 변위를 나타내며 x, y는 무게중심 기준으로 나타낸 상대적 절점 좌표를 나타낸다.

본 논문에서 좌표계는 무게중심에 놓여지고, 국부 요소의 회전자유도를 고려하기 위해 CR 좌표계를 설정하였다. 전역좌표계에서 CR 좌표계로의 변환에 따른 좌표계의 회전은 식 (1)과 같이 도출할 수 있다(Battini, 2008).

tanθ= i=13 [ xi(Yi+Vi-Yi-Vc)-yi(Xi+Ui-Xi-Uc) ] i=13 [ yi(Yi+Vi-Yi-Vc)+xi(Xi+Ui-Xi-Uc) ]       (1)

이때, θ는 회전변환행렬 [R]을 구성하고, 변형 전 좌표계에서 정의되는 θg는 회전변환행렬 [Rg]을 구성한다. 회전변환행렬을 사용하여 각각 식 (2), (3)과 같이 순수 변형의 병진자유도와 회전자유도를 나타낼 수 있다.

ui- υ i=RTXi+Ui-Xc-UcYi+Vi-Yc-Vc-xiyi      (2)

R-=RT [ Rg ]       (3)

최종적으로 회전변환행렬, [R]은 요소에 대한 회전변환행렬 [B]로 식 (4)와 같이 정의할 수 있고, 이를 도입하여 전역좌표계에서 정의된 강성행렬 https://cdn.apub.kr/journalsite/sites/jcoseik/2015-028-05/07TK102015280502/images/10.7734.28.5.485.F131.png과 내력 벡터 https://cdn.apub.kr/journalsite/sites/jcoseik/2015-028-05/07TK102015280502/images/10.7734.28.5.485.F132.png는 식 (5)와 같다.

B=P [ diag( [ R ] g) ]  T ,i=1:N      (4)

Kg=BT [ Kl ] [ B ] + [ Kh ]  , { fg } = [ B ]  T { fl }

where

Kh=δ [ diag( [ Ri ] ) ] [ P ] T { fl } + [ diag( [ Ri ] ) ] δ [ P ] T { fl }       (5)

행렬 [P]는 식 (6)과 같이 전역좌표계에서 정의된 변위를 국부 좌표계에 대한 미분으로 유도되며 전체 요소 거동에서 순수 변형에 대한 물리량을 추출한다. 이때, 벡터 {Ai}와 {Gi}는 절점 병진자유도에 의해 유도되며, 최종적으로 행렬 [P]는 식 (7)과 같이 구성할 수 있다(Battini, 2007).

[ Pi ] =δui-δui,gδui-δθi,gδθi-δui,gδθi-δθi,g  , i  =  1 :N      (6)

Pi=I3- { Ai } { Gi } T      (7)

국부요소는 OPtimal Triangular(OPT) 평면 요소(Felippa, 2003)를 적용하였다. 각 절점에 회전자유도를 포함하는 OPT 평면요소는 기존의 유한요소 정식화와 달리 자유 정식화(free formulation) 과정을 바탕으로 하여, 요소의 세장비에 관계없이 면 내, 굽힘 거동을 정확하게 예측 가능한 선형 요소이다. 또한 OPT 평면요소를 정의하는데 있어 평면 응력(plane stress) 연속 방정식을 적용하였다.

2.3 관성력에 의한 가상일

Fig. 1과 같이 정의된 요소의 좌표계에서 운동에너지는 아래와 같이 나타낼 수 있다.

K=12ρAx.Tx.dA=12ρhAx.Tx.dA+ρh312A θ  .Tθ.dA      (8)

식 (8)에서 ρ, A, h는 각각 밀도, 면적, 두께를 나타낸다. 운동에너지와 관성력에 의한 가상일의 관계에 따라 관성력 벡터는 다음과 같이 나타낼 수 있다.

δW=- δK  =δqTfk      (9)

ρhAdT'd.d+ρh312δθA [ θ..+θ.2 ] dA

q는 요소 절점자유도의 변분 항이며 병진자유도와 회전자유도로 구성된다. Fig. 1에서 정의한 좌표계에 따라 요소 절점자유도는 국부 좌표계에서 정의된 물리량으로 고려할 수 있다. 이 때, 요소의 절점자유도는 미소 물리량이므로 요소의 형상함수를 통해 요소 변위로 나타낼 수 있다.

δd=i=1NNiδdi-,   δθ =i=1NNiδθi-      (10)

최종적으로 요소의 관성력 벡터는 아래와 같이 나타낼 수 있다.

fk,l,i=ANiρhd-..dAρh36ANi [ (θ-..)2+θ-.. ] dA      (11)

국부좌표계에서 정의된 전단관성행렬은 국부 자유도에 대한 관성력의 미분으로 나타낼 수 있다. 따라서 관성력 벡터와 동일한 방법으로 요소자유도를 도출하여 적용하면 아래와 같이 관성행렬, 감쇠행렬을 표현할 수 있다.

Ml,i,j= [ ANiNjρhdA00ρh36ANiNj(2θ-.)dA ]       (12)

Cl,i,j= [ 000ρh36ANiNj(2θ-.)dA ]  , i,j   =1,2,3

이렇게 정의된 각각의 행렬과 벡터는 회전변환 행렬 [B]를 통해 전역좌표계에서 정의된 물리량으로 표현 가능하며 요소의 운동에 대한 비선형 지배방정식은 식 (14)와 같다.

Mg= [ B ] TMl [ B ] ,Cg= [ B ] TCl [ B ]       (13)

Mg(s)x..(t)+Cg(s)x..(t)+Kg(s)x..(t)=F(s,t)      (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와 같다. 해석결과를 살펴보면 기하학적 비선형적 거동특성이 잘 나타남을 확인할 수 있으며, 질 좋은 격자를 사용하여 도출된 결과와 비교하였을 때, 상대오차는 0.9%이다. 또한 기존의 선행 결과와 비교하였을 때, 본 연구를 통해 개발된 평면요소가 비 정렬 격자에서 더 좋은 성능을 나타냄을 확인하였다. 하중 인자가 1일 때, 각각의 결과 비교를 Table 1에 정리하였다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2015-028-05/07TK102015280502/images/10.7734.28.5.485.F123.png
Figure 2

Geometry and properties of an angle frame

https://cdn.apub.kr/journalsite/sites/jcoseik/2015-028-05/07TK102015280502/images/10.7734.28.5.485.F124.png
Figure 3

Comparison of the horizontal displacement at Point A

https://cdn.apub.kr/journalsite/sites/jcoseik/2015-028-05/07TK102015280502/images/10.7734.28.5.485.F125.png
Figure 4

Deformed configuration by the present analysis (λ=1.0)

이어서, Bathe 등(1975)에 의해 수행된 선행연구 사례를 고려하였다. 한 쪽 면이 지지된 평판에 윗면과 아랫면에 압력이 가해지는 조건이 적용되었다. 해석에 사용된 구조물의 형상과 물성은 Fig. 5와 같다. 해석결과, 본 연구를 통해 도출된 해석결과와 선행 연구결과가 잘 일치하는 것을 확인하였다. 해석결과는 Fig. 6에 나타내었다.

Table 1

Comparison of the horizontal displacement at Point A(λ=1.0)

Reference Displacement Difference(%)
Present6.710.95
Qm6(Battini, 2008)6.395.51
Qnew(Battini, 2008)6.494.03
Fine grid6.77-
https://cdn.apub.kr/journalsite/sites/jcoseik/2015-028-05/07TK102015280502/images/10.7734.28.5.485.F126.png
Figure 5

Cantilevered flat plate loaded by the surface pressure

https://cdn.apub.kr/journalsite/sites/jcoseik/2015-028-05/07TK102015280502/images/10.7734.28.5.485.F127.png
Figure 6

Comparison of the tip deflection

3.2 동적해석 검증

동적해석 검증을 위해 앞서 정적해석의 검증 예제로 수행되었던 Bathe 등(1975)의 예제를 고려하였다(Fig. 5). 시간에 따라 일정한 압력이 가해짐에 따라 발생하는 끝단의 변위 응답을 비교하였다. 본 연구에서 개발된 평면 요소는 총 20개가 사용되어 총 66개의 절점 자유도를 갖는 반면, 선행 연구에서는 8절점 사각 평면 요소 10개가 적용되어 총 84개의 절점자유도가 적용되었다. 상대적으로 적은 수의 자유도가 적용되었음에도 불구하고 해석결과, 선행 연구결과와 잘 일치하는 것을 확인하였다.

또한, 고차 요소를 사용한 결과와 비교 검증을 통해 본 요소의 효율성과 도출되는 해의 정확성을 확인하기 위해 3차원 고체 요소를 적용한 ANSYS 해석과의 비교를 시도하였다. 해석에 적용된 형상과 물성은 Fig. 8과 같고 시간에 따라 변화하는 끝단 집중하중이 가해짐에 따라 발생하는 끝단의 변형 응답을 비교 검증하였다. 가진 주파수는 10Hz로 선정하였다. 본 논문에서 개발한 평면요소는 총 20개의 요소를 사용하여 66개의 총 절점 자유도가 고려되었으나, 3차원 ANSYS 해석에는 총 80개의 요소를 사용하여 총 594개의 절점 자유도가 고려되었다. 비교 결과, 현 해석을 통해 도출된 끝단 응답결과가 ANSYS에 비하여 매우 적은 개수의 요소를 적용했음에도 불구하고 ANSYS 해석결과와 잘 일치하는 것을 확인하였다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2015-028-05/07TK102015280502/images/10.7734.28.5.485.F128.png
Figure 7

Transient response of the cantilevered flat plate under constant surface pressur

https://cdn.apub.kr/journalsite/sites/jcoseik/2015-028-05/07TK102015280502/images/10.7734.28.5.485.F129.png
Figure 8

Cantilevered plate under the harmonic tip load

https://cdn.apub.kr/journalsite/sites/jcoseik/2015-028-05/07TK102015280502/images/10.7734.28.5.485.F130.png
Figure 9

Transient response of the cantilevered plate under the harmonic tip load

4. 결 론

본 논문에서는 CR 기법을 적용하여 회전자유도를 갖는 3절점 비선형 삼각평면요소를 개발하였다. 국부 요소의 회전을 고려하기 위해 요소 형상의 무게중심에 CR 좌표계를 정의하였다. 또한 동적해석을 위해 운동에너지의 변분으로부터 전단 질량행렬, 감쇠행렬 그리고 관성력 벡터를 유도하였고 Newton-Raphson 기법과 HHT-α 내재적 시간적분 기법을 적용하여 지배방정식을 해석하였다. 수치해석은 정적해석과 동적해석의 몇 가지 예제를 사용하였고 해석결과는 선행 연구결과 및 상용프로그램 해석결과와 비교 검증하였다. 검증 결과, 정적해석에서의 처짐은 0.9%의 오차를 나타내었으며, 동적해석결과, 변위 응답이 3차원 고체 요소를 사용한 ANSYS 해석결과와 잘 일치하는 것을 확인하였다. 향후, 본 논문을 통해 개발한 평면 요소를 플래핑 날개의 구조해석에 적용하여 유체-구조 결합해석을 수행하고 나아가 국부 Lagrangian 승수를 적용한 다물체 동역학 해석으로 확장을 수행할 것이다.

Acknowledgements

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.

References

1
D.J., Allman, 1984. A Compatible Triangular Element Including Vertex Rotations for Plane elasticity Problems. Compos. Struct., 19, pp.1-8.
10.1016/0045-7949(84)90197-4
2
D.J., Allman, 1988. A Quadrilateral Finite Element Including Vertex rotations for Plane Elasticity Problems. Int. J. Nummer. Meth. Eng., 18, pp.717-739.
10.1002/nme.1620260314
3
K-J., Bathe, E., Ramm and E.L., Wilson, 1975. Finite Element Formulations for Large Deformation Dynamic Analysis. Int. J. Nummer. Meth. Eng., 9, pp.353-386.
10.1002/nme.1620090207
4
J-M., Battini, 2007. A Modified Corotational Framework for Triangular Shell Elements. Comput. Meth. Appl. Mech. Eng., 196, pp.1905-1914.
10.1016/j.cma.2006.10.006
5
J-M., Battini, 2008. A Non-linear Corotational 4-node Plane Element. Mech. Res. Commun., 35, pp.408-413.
10.1016/j.mechrescom.2008.03.002
6
P.G., Bergan and C.A., Fellipa, 1985. A Triangular Membrane Element with Rotational Degrees of Freedom. Comput. Meth. Appl. Mech. Eng., 50, pp.25-60.
10.1016/0045-7825(85)90113-6
7
N., Carpenter, H., Stolarski and T., Belytschko, 1985. A Flat Triangular Shell Element with Improved Membrane Interpolation. Comput. Meth. Appl. Mech. Eng., 1, pp.161-168.
10.1002/cnm.1630010405
8
M.A., Crisfield, 1990. A Consistent Co-rotational Formulation for Non-linear, Three-dimenstional, Beam-elements. Comput. Meth. Appl. Mech. Eng., 81, pp.2969-2992.
10.1016/0045-7825(90)90106-V
9
Crisfield M.A. 1997Non-linear Finite Element Analysis of Solid and StructuresAdvanced Topics2Wiley, London
10
C.A., Felippa, 2003. A Study of Optimal Membrane Triangles with Drilling Freedoms. Comput. Meth. Appl. Mech. Eng.,, 192, pp.2125-2168.
10.1016/S0045-7825(03)00253-6
11
C.A., Felippa and B., Haugen, 2005. A Unified Formulation of Small-strain Coroational Finite Elements: I. Theory. Comput. Meth. Appl. Mech. Eng., 194, pp.2285-2335.
10.1016/j.cma.2004.07.035
12
Pai P.F. 2007 Highly Flexible Structures: ModelingComputation, and Experimentation AIAA, IncVA
10.2514/4.861925
13
C.C., Rankin and B., Nour-Omid, 1986. An Element Independent Co-rotational Procedure for the Treatment of Large Rotations. ASME J. Pressure Vessel Techn., 108, pp.165-174.
10.1115/1.3264765
14
J., Simo and F., Armero, 1993. Improved Version of Assumed Enhanced Strain Tri-linear Element for 3d Finite Deformation Problems. Comput. Meth. Appl. Mech. Eng., 110, pp.359-386.
10.1016/0045-7825(93)90215-J
15
R., Taylor, P., Beresford and E., Wilson, 1976. A Non-conforming Element for Stress Analysis. Int. J. Nummer. Meth. Eng., 10, pp.1211-1219.
10.1002/nme.1620100602
페이지 상단으로 이동하기