(5) 12-16.pdf1.11MB

• 1. 서 론
• 2. 곡면 쉘요소의 정식화
• 2.1 기본가정 및 좌표계
• 2.1.1 전체좌표계(x,y,z)
• 2.1.2 자연좌표계(ξ,η,ζ)
• 2.1.3 절점좌표계 (Ve1I,Ve2I,Ve3I)
• 2.1.4 국지좌표계(x′,y′,z′)
• 2.2 쉘요소의 정식화를 위한 변분법
• 2.3 쉘요소의 기하학적 형상
• 2.4 변위장의 보간
• 2.5 공변 변형률 (Covariant Strain)
• 2.6 공변 변형률의 국지좌표계로의 변환
• 2.7 회전장의 보간과 스핀텐서
• 2.8 강성행렬
• 2.9 질량행렬
• 2.10 응력계산과 부재력
• 3. 수치시험
• 3.1 고유치 시험
• 3.2 조각보 시험 (Patch Test)
• 3.2.1 평면응력 조건에서의 조각보 시험
• 3.2.2 평판 조각보 시험
• 3.3 Cook의 평면 문제
• 3.4 원통형 지붕쉘(Scordelis-Lo Roof Shell)
• 3.5 집중하중을 받는 구멍이 있는 반구체 쉘
• 3.6 집중하중을 받는 원통형 쉘
• 3.7 개발된 쉘요소의 수렴 성능
• 3.8 응용예제
• 4. 결 론
• 감사의 글

1. 서 론

종래에는 원전구조물의 내진설계를 위한 설계 층응답스펙트럼을 산정하기 위하여 Frame 요소를 이용하는 방법이 많이 사용되어 왔으나, 내진설계 과정이 고도화됨에 따라 그림 1과 같은 쉘요소(Shell Element)를 이용한 3차원 해석 모형의 사용이 증가하고 있다. 쉘요소는 다른 유한요소에 비해 거동특성이 복잡하므로 최근까지도 지속적으로 개선이 이루어졌으며, 면내회전자유도가 포함된 쉘요소의 경우 80년대 후반 이후에 본격적으로 연구되기 시작하였다.^(1-17)

한편 2007년 개정된 SRP 기준⁽¹⁸⁾에 따라 거의 모든 부지에 대한 지반-구조물 상호작용 (SSI, Soil-Structure Interaction) 해석이 일반화됨에 따라 원전시설물의 내진설계분야에서 이에 대한 해석 수요가 크게 증가하고 있다. 하지만, 기존의 지명도 있는 SSI용 해석프로그램들이 1980년대에 개발된 점을 볼 때 이들 프로그램에서 최신의 쉘요소 연구 결과들을 반영하지 못할 수 있다고 판단된다. 따라서, 최근의 쉘요소를 이용한 내진해석 추세를 반영하기 위해서는 SSI해석용 특수프로그램에서도 정확성 향상을 위하여 최신의 쉘요소 기법을 적용할 필요성이 있다. 이러한 이유로 이 연구가 시작되었으며, 이 논문에서는 국내에서 SSI해석용 특수프로그램으로 개발된 KIESSI-3D⁽¹⁹⁾에 추가될 쉘요소의 정식화와 검증과정을 기술하였다.

<그림1> 쉘요소를 이용한 원전구조물의 3차원 모델링 예

 이 연구에서는 감절점쉘요소(Degenerated Shell Element)의 개념에 근거하여 쉘요소를 정식화하였으며, 회전장이 독립변수로 도입된 범함수(Functional)에 의하여 면내회전자유도(Drilling Rotational Degrees of Freedom)를 도입함으로써 절점당 6자유도를 갖는 쉘요소를 개발하였다. 쉘요소의 면내거동(In-plane Behavior) 개선을 위하여 4개의 비적합변위형(Nonconforming Displacement Modes)에 의한 비적합변위를 면내방향의 변위성분에 추가하였고, 면외거동(Out-of-plane Behavior) 개선을 위하여 대체전단변형률장(Substitute Shear Strain Field)을 적용하였다. 그리고 비적합변위형을 구현하기 위하여 비적합변위형의 직접수정법(Direct Modification Method)을 적용하였으며, 쉘요소의 부피적분을 위한 수치적분에는 9점 적분법(9-points Integration Rule)이 사용되었다.

 이 연구에서의 쉘요소와 가장 관련성이 큰 기존의 연구로는 최창근 등^(1),(7),(8)의 회전자유도를 갖는 8절점 입체요소(Solid Element)가 있다. 회전자유도를 갖는 입체요소의 정식화에 사용된 비적합변위형의 직접수정법과 회전자유도의 도입방법이 쉘의 면내거동 개선을 위해 거의 동일하게 적용되었다.

백종균⁽²⁾의 면내회전자유도를 갖지 않는 감절점쉘요소의 경우 면외거동을 개선하기 위해서 대체전단변형률장이 사용된 것은 이 연구에서의 접근방법과 동일하나, 백종균은 면내거동을 개선하기 위해서는 대체막변형률장(Assumed Membrane Shear Strain)을 사용하였고, 이 연구에서는 면내거동의 개선을 위해서 비적합변위형만 적용하였다.

이태열⁽³⁾ 쉘요소의 경우 최창근⁽¹⁾의 입체요소 정식화에서의 회전자유도 도입방법과 비적합변위형의 직접수정법을 적용하고, 강성행렬의 부행렬(Sub-matrix)별로 다른 수치적분을 적용함으로써 편평한 쉘요소(Flat Shell Element)를 구성하였다. 이태열의 쉘요소는 평판요소의 회전장에 비적합변위형을 적용하고 추가적으로 대체전단변형률장을 적용하였으나, 이 연구에서는 쉘요소의 면외거동개선을 위해서는 대체전단변형률장만 사용하였다.

 위에서 언급한 비적합변위장을 적용한 기존의 연구들에서는 전체좌표계 방향의 변위장에 비적합변위를 추가 적용하던 방식이었으나, 이 연구에서는 곡선좌표계인 자연좌표계 방향의 변위장에 비적합변위를 적용하였다. 이러한 자연 좌표계방향 변위장에 비적합변위를 추가함에 따라 곡면을 갖는 쉘요소에서 면내방향으로만 비적합변위형을 적용하는 것이 가능하게 되었다.

2. 곡면 쉘요소의 정식화

2.1 기본가정 및 좌표계

 이 연구에서 사용한 쉘요소의 정식화에 있어서 변형전 중립면(Neutral Surface)에 수직인 직선면은 변형 후에도 직선을 유지하는 것으로 가정하고, 중립면에 수직방향으로 작용하는 응력에 의해 생성된 변형에너지는 무시하였다.

 이 연구에서는 감절점 쉘요소의 기본 개념을 이용하여 편평하지 않은 4절점 쉘요소를 개발하였으나, 쉘요소의 형상을 정의함에 있어서 쉘요소의 중립면에 위치한 4개 절점에서의 좌표와 두께를 이용하여 그림 2와 같이 쉘요소의 기하학적 형상을 정의하였다.

 감절점쉘요소의 경우 쉘요소의 상면과 하면에서의 좌표차이를 이용하여 쉘 두께방향을 정의하는 것이 보통이다. 하지만 이러한 접근방법은 일반적으로 널리 사용되는 상용 소프트웨어에서의 입력방식과 달라 상용 유한요소프로그램의 모델링 기능을 이용하기 곤란한 단점을 갖는다. 따라서 이 연구에서는 쉘의 기본적인 정식화과정은 감절점쉘요소의 접근법을 이용하지만, 각 절점에서의 두께방향은 절점좌표로부터 정의되는 중립면에 수직한 방향으로 가정하고, 가정된 쉘두께 방향벡터와 두께값을 이용하여 쉘의 윗면 및 아랫면의 형상을 정의하였다.

 이 연구에서의 쉘요소 정식화에 있어서 다음과 같은 4 종류의 좌표계를 사용하였다.

2.1.1 전체좌표계(x,y,z)

 구조물의 기하학적 형상을 정의하기 위해 사용된 직각좌표계(Cartesian Coordinate System)이다. 전체좌표계를 이용하여 구조물의 강성행렬 및 하중벡터의 조립이 이루어지며, 이 좌표계에서 운동방정식의 해를 구한다.

<그림2> 쉘요소의 기하학적 형상 및 자연좌표계

<그림3> 절점 I에서의 절점좌표계 개념도

2.1.2 자연좌표계(ξ,η,ζ)

 자연좌표계(Natural Coordinate System)는 그림 2와 같이 쉘요소의 중립면에 따른 곡선좌표계(Curvilinear Coordinate System) ξ및 η와 쉘요소의 두께방향에 대한 좌표축 ζ로 이루어진다. 이 자연좌표계에서는 ξ,η및 ζ값이 -1에서 +1 사이의 값을 갖는다.

2.1.3 절점좌표계 (Ve1I,Ve2I,Ve3I)

 절점좌표계는 쉘 요소의 중립면에 위치한 절점에서 정의되는 좌표계이다(그림 3). I절점의 중립면에 수직한 단위벡터 는 자연좌표계 ξ축 및 η축에 접하는 벡터의 벡터곱(Cross Product)으로부터 정의한다.

 전체좌표계상의 좌표 x, y및 z를 x_i(즉, x₁ , x₂및 x₃)로 표현하고, 전체좌표계의 각 방향의 기저벡터(Base Vector)를 각각 e₁, e₂및 e₃(혹은 e_x, e_y 및 e_z)라고하면 전체좌표계상의 임의 점에서의 위치벡터 x는 x=∑x_ie_i로 표현할 수 있다. 전체좌표계로 표현된 I절점의 좌표를 x^e_I 라고하면 I절점에서 중립면에 수직한 단위벡터는 다음과 같이 정의할 수 있다.

여기서 연산자 ×는 벡터곱이며, ∥⋅∥은 벡터의 Norm으로서 벡터의 크기이다. 이렇게 계산된 벡터 는 I절점의 절점좌표계를 표현하는 추가적인 단위벡터및 를 정의하는데 사용된다. 만약 가 전체좌표계의 y방향의 기저벡터인 e_y와 평행하다면 = e_z로 정의하고, 그 외의 경우는 다음의 식에 의해 벡터를 정의한다.

 절점좌표계를 표현하는 벡터는 위에서 정의한 와 를 이용하여 다음과 같이 정의한다.

따라서,,로 표현되는 절점좌표계는 직교좌표계(Orthogonal Coordinate System)이다.

2.1.4 국지좌표계(x′,y′,z′)

국지좌표계(Local Coordinate System)는 각 적분점에서 응력과 변형률을 계산하는데 사용하는 직교 좌표계이다. 쉘요소의 경우 쉘의 두께방향인 국지좌표계 z′방향으로의 수직응력 σ_z′가 영(즉,  σ_z′=0)이라는 가정에 근거하여 요소의 정식화가 이루어지므로 국지좌표계에서 응력과 변형률의 관계를 정의하는 것이 보통이다. 이 연구에서는 절점위치에서는 국지좌표계와 절점좌표계가 일치하도록 정의하였으며, 절점위치가 아닌 요소 내 임의의 중립면 위치에서 절점좌표계는 절점좌표계와 유사하게 정의된다.

국지좌표계 z′방향의 기저벡터 e_z ′는 자연좌표계 ξ축 및 η축에 접하는 벡터의 벡터곱으로 정의한다.
 

위의 식에서 ξ_i 는 자연좌표계를 첨자로 표시한 것이며,  ξ₁ = ξ,  ξ₂ = η,  ξ₃= ζ의 의미로 사용되었다.

만약 e_z′ 가 전체좌표계 y방향의 기저벡터인 e_y와 평행하다면 e_x′= e_z로 정의하고, 그 외의 경우는 다음 식에 의해  e_x′벡터를 정의한다.

국지좌표계와 전체좌표계 사이의 좌표변환을 수행하는 방향여현행렬(Direction Cosine Matrix) [Q]는 국지좌표계의 단위 기저벡터를 이용하여 다음과 같이 정의할 수 있다.
 

2.2 쉘요소의 정식화를 위한 변분법

 이 연구에서는 쉘요소의 면내회전자유도를 도입하기 위하여 독립적인 회전장(Rotation Field)이 고려된 범함수를 사용하였다.^(4),(5) Hughes와 Brezzi⁽⁴⁾에 의해 제안된 범함수는 쉘 문제에서 다음과 같이 표현될 수 있다.

 여기서, {E′}는 변위장(Displacement Field)의 미분을 통하여 계산되는 국지좌표계에서의 변형률 벡터,

 [D′] 는 국지좌표계에서의 재료구성방정식과 관련된 탄성 행렬(Constitutive Matrix), θ_z′ 는 국지좌표계 z′축 방향의 회전장(Rotation Field), µ는 전단탄성계수, W는 외력이 한 일(Work)이다. 그리고 a는 문제종속적인 상수로서 이 연구에서는 0.01의 값을 사용하였다. [D′] 행렬은 등방성 재료인 경우 탄성계수 E와 포아송비 ν에 의해 다음과 같이 표현된다.

2.3 쉘요소의 기하학적 형상

 일반적으로 감절점쉘요소에서의 쉘요소의 형상은 쉘요소의 중립면에서의 절점좌표와 쉘 두께를 이용하여 표시하면 다음과 같다.⁽²⁾

 여기서, t^e_I 는 절점 I에서의 쉘요소의 두께이고, x^e_I, y^e_I및 z^e_I는 절점 I에서 x, y및 z방향 절점좌표이다. 또한, V^e_3Ix,V^e_3Iy및 V^e_3Iz는 각각 절점 I에서 쉘두께방향의 절점좌표계 방향벡터의 전체좌표계 x, y및 z방향 벡터성분이다. 또한 N_I 는 절점 I에서의 형상함수로서 다음과 같다.

여기서 ξ_I 및 η_I 는 각각 절점 I에서의 국지좌표계 ξ 및 η의 값이다.

 연쇄법칙(Chain Rule)에 의해 자연좌표계로의 편미분과 전체좌표계로의 편미분간의 관계는 다음과 같이 자코비안행렬(Jacobian Matrix) [J]에 의해 표현되며, 위에서 정의한 기하학적인 형상은 자코비안행렬의 계산에 이용된다.

2.4 변위장의 보간

 쉘요소의 정식화(Formulation)에서 요소의 적합변위장은 중립면에 위치한 각 절점에서 변위와 회전변위를 이용하여 아래와 같이 표현한다.⁽⁶⁾

여기서, u^e_I, v^e_I및 w^e_I는 각각 절점 I에서의 전체좌표계 x,y및 z방향의 변위이고, β^e_1I, β^e_2I및 β^e_3I는 절점 I에서의 절점좌표계 방향벡터 V^e_1I, V^e_2I및 V^e_3I방향 절점회전각(Nodal Rotation Angle)이다(그림 4).
이 적합변위장의 보간을 간단히 정리하면 다음과 같다.

 여기서 {u^h }는 전체좌표계에서의 변위장, {d_I }는 절점 I에서의 전체좌표계에서의 절점변위와 절점좌표계에서의 절점 회전, {d}는 쉘요소 절점변위벡터, 형상함수행렬 [N] 은

<그림4> 절점 I에서의 절점변위 및 절점회전 개념도

 이다. V^e_1Ix, V^e_1Iy및 V^e_1Iz는 각각 벡터 V^e_1I의 x, y및 z방향의 성분이며, V^e_2Ix, V^e_2Iy및 V^e_2Iz는 각각 벡터 V^e_2I의 x, y 및 z방향의 성분이다.

 이 연구에서는 쉘요소의 면내거동을 개선하기 위한 방법으로 다음과 같은 4개의 비적합변위형을 사용하였다.⁽⁸⁾

 이 연구에서 다루는 4절점 쉘요소는 편평한 쉘(Flat Shell)이 아니므로, 전체좌표계상에서 면내방향으로 비적합변위를 직접 정의하기 보다는 곡선좌표계인 자연좌표계에서 비적합변위를 정의하는 것이 편리하다. 이 연구에서의 비적합변위형은 쉘의 면내거동을 개선하기 위해서 자연좌표계 ξ및 η방향의 변위성분에 적용되었으며, 자연좌표계 ξ및 η방향으로 추가되는 비적합변위장 및 는 다음과 같다

 여기서, 및 는 각각 I 번째 비적합변위형의 진폭과 관련된 매개변수이다.

 이 연구에서는 비적합변위장의 수치적 구현에 있어서 일정한 변형률 상태를 표현할 수 있도록 하기위해 비적합변위형의 직접수정법을 적용하였다.^{(1),(3),(7-9)}따라서, 이 연구에서의 비적합변위형 의 편미분과 관련된 항은 다음과 같이 수정된 값이 수치구현에 사용되었으며, 비적합변위형의 자연자표계로 편미분한 값  를 전체좌표계로 편미분한 값 로 변환하는 관계식에서는 자연좌표계의 원점에서 구한 자코비안행렬의 역행렬이 사용되었다.

 여기서, J는 자코비안(즉, J= det([J])이고, [J]는 자코비안 행렬)이고, J₀ 은 자연좌표계 원점에서의 자코비안 값이다.

2.5 공변 변형률 (Covariant Strain)

이 연구에서는 쉘요소의 면내거동 개선을 위해서 곡선좌표계인 자연좌표계 ξ및 η방향으로 비적합변위형을 적용하고, 면외거동 개선을 위해서는 대체전단변형률장을 적용하였다. 이러한 이유로 요소의 정식화 과정에서 우선적으로 자연좌표계상의 변형률을 정의하고 필요에 따라 다른 좌표계로 좌표변환을 하는 것이 편리하다.

전체좌표계의 변위성분을 자연좌표계로 미분한 값으로부터 자연좌표계의 변위성분을 자연좌표계로 미분한 값을 구하는 것은 다음의 관계식에 의해 산정할 수 있다.⁽²⁾ 

 여기서, 는 자연좌표계 ξ, η, Ϛ 방향 변위함수이며, [J] 은 자코비안행렬이다.

위의 관계식으로부터 자연좌표계의 변위를 자연좌표계로 미분한 {}벡터와 전체좌표계의 변위를 자연좌표계로 미분한 {E}벡터와의 관계는 다음과 정리할 수 있다.
 

여기서 

이고, [I]는 3X3단위행렬(Unit Matrix), J_ij 는 자코비안 행렬 [J]의 i행 j열의 성분을 의미한다. 변위장을 편미분하여 구할 수 있는 변위미분벡터와 절점변위와의 관계식은 다음과 같이 간단히 표현될 수 있다.
 

여기서
 

이고, [G]는 적합변위의 미분과 절점변위와의 관계행렬로서 변위장의 보간에 사용된 형상함수행렬을 미분하여 결정할 수 있고, []는 비적합변위의 자연좌표계로의 편미분과 비적합변위형의 내부변위간의 관계행렬로서 비적합변위장의 미분을 통하여 결정할 수 있다.

 자연좌표계상의 변형률인 공변변형률은 자연좌표계상의 ξ, η 및  Ϛ방향 변위성분 및 의 편미분에 의해 구할 수 있다.

자연좌표계상의 공변변형률 {}은 자연좌표계상의 변위  및 를 자연좌표계  ξ, η 및  Ϛ로 미분한 변위의 미분 {}으로 다음과 같이 나타낼 수 있다. 

여기서
 

이제 식 (29)과 {} = [J_E ]{E}관계에 의해 자연좌표계에서 공변변형률과 절점변위와의 관계식(Strain-Displacement Relationship Matrix)은 다음과 같이 표현할 수 있다.

 

여기서
 

한편, 대체전단변형률장은 쉘요소의 4개의 변(ξ=±1, η =±1)에서 결정되는 적합변형률을 선형적으로 보간하여 요소의 내부에 위치한 임의 지점의 적합공변변형률을 대체한 것이다. 이 연구에서는 쉘요소의 4개 변상의 적합변형률을 대체전단변형률장 산정에 사용하였으며, 요소의 변상에서 구한 전단변형률과 대체전단변형률 사이의 관계는 다음과 같다.^(2),(10)

2.6 공변 변형률의 국지좌표계로의 변환

자연좌표계에서 정의된 변형률을 전체좌표계 및 국지좌표계로 변환하는 것은 텐서변환 관계에 의해 이루어진다. 곡선좌표계인 자연좌표계의 전체좌표계로의 변환관계식은 다음과 같다.
 

여기서, [J^-1 ]는 자코비안행렬의 역행렬이며, 변형률텐서행렬 은 3X크기의 대칭행렬이다. 또한, 직교좌표계간의 좌표변환인 전체좌표계의 국지좌표계로의 변환은 다음과 같은 변형률텐서의 변환관계가 성립한다.
 

여기서, [Q]은 국지좌표계의 기저벡터를 이용해 정의한 방향여현행렬이다. 따라서, 국지좌표계에서 정의된 변형률을 국지좌표계로 변환하는 관계식은 다음과 같다.
 

이 연구에서는 [ ϵ′] = [Q][ϵ][Q]^T  변환관계의 수치구현에 있어서 2계 텐서인 변형률텐서를 Voigt 표기법으로 벡터로 표현한 다음의 변환 관계를 좌표변환에 이용하였다.

여기서
 

 이고, [T_Q ]는 다음과 같다.

위의 [T_Q]변환행렬에서 [Q]대신 [J^-1 ]가 사용된 변환행렬을 [T_J^-1]라고하면, 다음과 같은 자연좌표계에서의 전체좌표계로의 변형률 변환식이 성립한다.

 따라서 자연좌표계에서 정의된 변형률을 국지좌표계로 변환하는 관계식은 다음과 같다.

즉
 

여기서
 

수치적구현에 있어서 비적합변위형의 직접수정법 개념에 따라 비적합변위와 관련된 행렬계산과 관련된 [T_J^-1  ]는 원점에서의 자코비안행렬의 역행렬을 이용하여 구성되었다.
 

2.7 회전장의 보간과 스핀텐서

전체좌표계 x, y, z축 방향의 회전 θ_x , θ_y , θ_z 는 형상함수와 절점 I에서의 전체좌표계의 절점회전 θ^e_lx , θ^e_ly , θ^e_lz를 이용하여 다음과 같이 보간할 수 있다.
 

 절점 I에서의 절점좌표계에서의 절점회전 B^e_il를 전체좌표계에서의 절점회전 θ^e_il로 변환하기 위해서는 다음과 같은 변환관계가 사용된다.

여기서
 

따라서, 절점좌표계로 표현된 절점변위/회전과 전체좌표계에서의 변위장과의 관계는 다음과 같다.
 

여기서

국지좌표계에서의 회전 θ_x' , θ_y', θ_z'는 국지좌표계의 방향여현행렬 [Q]를 이용하여 다음과 같이 구할 수 있다.

 즉

따라서 앞서 언급한 범함수에서 독립변수로 도입된 회전장 θ^h_z'는 다음과 같이 표현된다.
 
즉

여기서

한편, 3차원 연속체 문제에서의 회전은 전체좌표계의 x_k축을 둘레로 회전하는 회전 Ψ_k은 변위장의 미분으로 구해지는 스핀텐서(Spin Tensor)와 순환기호(Permutation Symbol) ϵ_ijk 를 이용하여 다음과 같이 정의할 수 있다.

즉

 여기서

이며, 행렬 [h]은 다음과 같다.

전체좌표계의 변위를 전체좌표계로 미분한 벡터 {E^*}와 전체좌표계의 변위를 자연좌표계로 미분한 벡터 {E}는 자 코비안 행렬의 역행렬 [J^-1 ]을 이용하여 다음의 변환관계로 계산할 수 있다.
 
여기서

쉘요소의 경우 면에 수직한 두께방향이 국지좌표계 z′방향으로 표현되므로 위의 회전 Ψ_k를 국지좌표계의 회전 Ψ'_k로 좌표변환하는 것이 필요하다. 이 좌표변환은 앞에서 정의된 방향여현행렬 [Q]를 이용하여 다음과 같이 표현된다.
 


여기서

 이상의 관계식들과 {E} = [J_E ]^-1 {}을 이용하여 국지좌표계에서의 회전을 구하면 다음과 같다.

따라서 국지좌표계 z′축 회전 Ψ_z'는 다음과 같이 계산된다.

최종적으로

여기서

수치적구현에 있어서 비적합변위형의 직접수정법 개념에 따라 비적합변위와 관련된계산과 관련된 [J^*_E ]는 원점에서의 자코비안행렬의 역행렬을 이용하여 구성되었다.
  

2.8 강성행렬

앞서 정의된 변형률 및 회전장과 관련된 식들을 범함수에 대입하여 정리하면 다음과 같은 범함수를 얻는다.
 

 여기서 일 W는 절점변위 {d}와 절점하중 {f}의 스칼라곱인 W = {f}^T {d}으로 나타낼 수 있다. 위의 범함수에 정류 조건을 적용하면 다음과 같은 방정식을 얻는다.

여기서

 이 연구에서는 쉘요소의 부피적분에 9점적분법을 사용하였다.⁽⁸⁾ 


여기서, W₀  = 0.001, W_α  = 1- W₀ /8, α = 1/이다.

요소종속적인 내부변위로 볼 수 있는 비적합변위와 관련된 자유도를 정적 응축을 통하여 소거하면, 다음과 같은 쉘요소의 강성행렬을 구할 수 있다.
 

식 (66)의 쉘요소 자유도는 병진변위는 전체좌표계이고, 회전은 절점좌표계에 대한 것이다. 따라서 절점좌표계에서의 절점회전 β^e_il를 전체좌표계에서의 회전변위 θ^e_il로 변환하기 위해서는 절점변위의 변환관계가 사용된다. 전체좌표계에서의 최종적인 쉘요소의 강성행렬 [K]는 다음과 같다.
 

2.9 질량행렬

 체적력(Body Force)은 단위 체적당의 하중으로서 중력과 같은 정역학적인 체적력과 관성력과 같은 동역학적인 체적력이 있다. 만약 체적력이 관성력이고, 물체의 밀도를 ρ,물체의 임의 지점에서의 가속도를 ϋ라고하면, D'Alembert원리에 의하여 관성력에 의한 체적력이-ρϋ로 표현된다(ρ=질량밀도). 체적력에 의한 등가절점하중을 {f}, 가상절점변위를 δ{d}, 임의 점에서의 가상변위를 δ{u}라고하면, 가상일의 원리에 의하여 다음과 같은 내적가상일과 외적가상일에 관한 방정식이 정의된다.

 변위장 {u}를 절점변위벡터 {d}를 이용하여 {u} = [N]{d}로 보간하면, 가속도장 {}은 절점가속도벡터 {}를 이용하여 {} = [N]{}로 표시된다. 따라서 위의 가상일의 방정식은 다음과 같이 정리된다.


즉

 여기서 질량행렬 [m] 은 변위장의 보간에 사용되는 형상함수행렬을 이용해서 다음과 같이 정의된다.

따라서 전체좌표계에서의 최종적인 쉘요소의 질량행렬 [M]는 다음과 같다.
 

2.10 응력계산과 부재력

 비적합변위형이 사용된 유한요소에서의 응력은 절점변위 {d′}와 비적합변위형의 진폭 {}를 이용하여 다음과 같이 산정한다.

 여기서, {E′₀ }와 {σ′₀}는 각각 임의 초기 변형률과 초기응력을 표시한다. 비적합변위형의 진폭 {}은 절점변위와 강성 행렬 산정시에 사용되는 부행렬 [k_NN ] 및 [k_NC ]을 이용하여로 산정할 수 있다. ⁽⁶⁾ 

쉘요소의 단위길이당 부재력/모멘트는 쉘요소의 국지좌표계의 응력을 적분하여 다음과 같이 산정하였다(그림 5).
 

<그림5> 쉘요소의 단위길이당 부재력/모멘트 표기법

 ① 면내(In-plane) 단위길이당 축력 및 전단력


② 두께방향 단위길이당 전단력


③ 단위길이당 휨모멘트 및 비틀림

3. 수치시험

3.1 고유치 시험

바람직하지 못한 영에너지모드(Spurious Zero Energy Mode)가 존재하는지의 여부를 확인하기 위해 찌그러지지 않은 정사각형 형상의 단일 쉘요소의 강성행렬에 대한 고유치해석을 수행하였다. 경계조건이 없는 상태에서 3차원 쉘요소는 총 6개의 강체운동과 관련된 영에너지 모드들을 갖는 것으로 나타났다. 따라서 바람직하지 못한 영에너지모드는 없는 것으로 판단된다.

3.2 조각보 시험 (Patch Test)

개발된 쉘요소가 일정한 변형률 상태를 표현할 수 있는가의 여부를 확인하기위하여 조각보시험을 수행하였다.

3.2.1 평면응력 조건에서의 조각보 시험

 쉘요소의 면내거동과 관련된 평면응력 상태에서 일정한 변형률을 표현할 수 있는지의 여부를 확인하기 위하여5개의 유한요소로 이루어진 해석모델에 u= 10^-3(x + y / 2) , v= 10^-3(x + y / 2)의 경계조건을 가하는 조각시험을 수행하였다.⁽¹¹⁾ 탄성계수 E= 1.0 X 10⁶, 포아송비 ν = 0.25, 두께 t=0.001의 조건에서 응력과 변형률은 각각 E_x = E_y = γ_xy = 10^-3 , σ_x =  σ_y = 1333, T_xy = 400 의 이론 해를 갖는다(그림 6).

 수치 시험 결과 이 연구에서 개발된 쉘 요소는 이론해와 일치하는 결과를 보였으며, 조각보시험을 통과함을 알 수있었다.

<그림6> 평면 응력 조각보 시험에 사용된 요소체눈

3.2.2 평판 조각보 시험

 쉘요소의 면외거동과 관련된 조각시험을 위하여 Hinton과 Huang⁽¹²⁾이 제안한 휨모멘트, 전단력, 비틀림하중의 세 가지 하중 경우에 대하여 각각 평판 조각보 시험을 수행하였으며, 이 연구에서의 쉘 요소는 모든 경우의 조각보 시험을 통과하였다.

3.3 Cook의 평면 문제

 이 평면 문제는 쉘 요소의 면내 방향 거동을 평가하는 문제이다. 이 문제에서는 해석대상의 두께는 1.0이며, 한쪽 면은 고정되어 있고 다른 면은 크기가 1.0인 분포하중이 작용한다(그림 7). 재료의 탄성계수 1.0이며 포아송비는 0.33이다. 끝단의 중앙점 A에서의 수직변위에 대한 이론치는 23.91이다. 해석결과는 표 1에 나타내었다.

이 연구에서의 쉘 요소의 경우 면내거동개선을 위해 사용한 비적합변위형의 영향으로 양호한 거동을 보여주었다.

<그림7> Cook의 평면 문제

<표1> Cook의 평면 문제의 끝단 중앙에서의 수직변위 비교

<표2> 원통형 지붕쉘의 변위비교

<그림8> 원통형 지붕쉘과 재료특성

<그림9> 구멍이 있는 반구체 쉘

3.4 원통형 지붕쉘(Scordelis-Lo Roof Shell)

 곡면의 양단은 격판(Diaphragm)으로 지지되고 원통형 길이방향의 가장자리는 자유지지된 원통형 지붕쉘에 자중이 작용하는 경우에 대해 해석을 수행하였다(그림 8). 지붕의 길이는 50, 반경은 25, 두께는 0.25이며, 지붕에는 단위 면적당 90인 자중이 작용한다. 재료의 탄성계수는 4.32×10⁸이며 포아송비는 0.0이다. 이 문제에서의 양단은 격판으로 지지되어 있다. 이 문제에서의 거동특성은 면내응력상태가 대부분이고 비신장휨(Inextensible Bending)이 작게 작용한다. 이 지붕의 변의 중앙점 (B점)에서의 수직변위의 이론값은 0.3024로 알려져 있다.⁽¹¹⁾ 이 문제는 대칭성을 이용하여 1/4부분만 모델링하여 해석하였으며, 자중에 의한 힘은 질량행렬 구성후 중력가속도를 곱하여 산정하였다. 수치시험 결과는 표 2에 나타내었다.

 이 문제의 해석결과들을 비교해볼 때 면내거동개선을 위해 사용한 비적합변위형이 면내응력성분이 대부분인 원통형 지붕쉘 문제의 거동개선에 효과가 있음을 알 수 있다.

3.5 집중하중을 받는 구멍이 있는 반구체 쉘

이 MacNeal과 Harder⁽¹¹⁾에 의해 제안된 이 문제는 비신장 휨과 큰 강체회전의 거동을 시험하는 예제이다(그림 9). 이 문제에서의 재료의 탄성계수는 6.825×10⁷, 포아송비는 0.3이고, 해석대상 구조물의 쉘의 반경은 10.0, 두께는 0.04이다. x축방향의 하중재하점에서의 하중방향 변위를 참고값 0.093과 비교하여 표 3에 나타내었다.

<표3> 구멍이 있는 반구체 쉘 변위비교

<그림10> 자유단을 갖는 원통형 쉘 모델과 재료특성

<그림11> 격판을 갖는 원통형 쉘 모델과 재료특성

<표4> 자유단을 갖는 원통형 쉘의 변위 비교

이 문제의 경우 면외거동이 지배적인 문제이므로 제안된 요소에서 채택한 전단대체변형률장의 영향으로 양호한 결과를 보여주는 것으로 판단된다.

3.6 집중하중을 받는 원통형 쉘

이 원통형 쉘 문제는 비신장휨과 복합적인 면내응력상태에 대해 요소의 거동을 테스트하는 예제이다. 여기서는 양단의 지지형태가 자유단인 경우(그림 10)와 격판으로 지지된 경우(그림 11)의 두 가지 경우에 대하여 해석을 수행하였다. 하중과 기하학적형상이 대칭인 문제이므로 1/8모델이 사용된다.

 먼저 양단이 자유단인 경우에는 원통의 길이가 10.35, 반경은 5.0, 두께는 0.094이다. 사용된 재료의 탄성계수는 1.05×10⁷ , 포아송비는 0.3125이다. 이 경우의 하중재하점의 수직변위의 이론치는 0.1139로 알려져 있으며(2), 해석결과는 표 4에 나타내었다.

양단이 격판으로 지지된 경우에는 원통의 길이가 600, 반경은 300, 두께는 3이고, 재료의 탄성계수 3.0×10⁶이고 포아송비는 0.3이다. 이 문제에 대한 이론치는 하중재하점의 수직변위가 0.18248×10^-4인 것으로 알려져 있으며⁽²⁾, 해석 결과는 표 5에 나타내었다.

<표5> 격판을 갖는 원통형 쉘의 변위 비교

<그림12> 개발된 쉘요소의 체눈세분화에 따른 해의 수렴곡선

양단이 자유단인 경우에 비해 양단이 격판으로 지지된 경우에는 다른 예제들에 비해 모든 요소들의 수렴속도가 느린 것을 알 수 있으며, 이 예제로부터 구속조건이 많은 문제에서는 보다 조밀한 유한요소체눈(Finite Element Mesh)이 필수적인 것을 알 수 있다. 

3.7 개발된 쉘요소의 수렴 성능

 앞에서 수행한 검증예제들에 대한 개발된 쉘요소의 수렴 특성을 평가하기 위하여, 앞의 해석결과를 그림 12와 같은 수렴곡선으로 정리하였다. 이 그림으로부터 이 연구에서 제안한 쉘요소는 요소의 세분화됨에 따라 해가 수렴해가는 것을 확인할 수 있다.

3.8 응용예제

KIESSI-3D에 이식된 개발된 쉘요소의 적용성을 보이기 위하여 그림 13과 같은 원자로건물 해석모델⁽²⁰⁾ 에 대한 x-방향 수평지진응답해석을 수행하였다. 원자로 격납건물의 반경은 18.9m, 높이 63.1m, 원자로 외벽두께는 1.07m, 원자로 돔 두께는 0.76m이고, 지반조건은 경암 위에 두께 21.5m인 수평층상 연암이 놓인 조건이며, 수평지진입력의 통제점은 자유장지반에서 경암과 연암의 경계면이다.

<그림13> SASSI NPP 예제의 KIESSI-3D 해석 모델

해석모델의 차이(쉘요소 해석모델과 Beam요소 해석모델)가 지진응답에 미치는 영향을 평가하기 위하여, 두 해석 모델에 의한 격납건물 응답의 통제운동에 대한 수평변위의 전달함수를 격납건물 상단(그림 13에서 A점)과 원자로 돔 하단(그림 13에서 B, B1, B2점)에서 계산하고 이를 그림14와 같이 비교하였다. 비교결과, 격납건물 상단에서 수평 응답은 해석모델에 따른 차이가 거의 없는 것으로 나타났으나, 돔 하단에서는 수평응답은 쉘요소 모델의 3차원 변형효과로 인하여 두 번째 Peak에서 진폭의 차이가 첫 번째 Peak에 비해 상대적으로 큼을 알 수 있었다. 아울러 이 예제에서는 Beam 해석모델에 의한 전달함수의 진폭이 가장 커서 Beam 해석모델은 쉘요소를 이용한 3D 해석모델에 비해 보수적인 응답을 나타냄도 알 수 있었다.
 

<그림14> 전달함수 비교

4. 결 론

 이 연구에서 개발된 4절점 곡면 쉘요소는 단일요소의 고유치해석결과 강체운동과 관련된 6개의 영에너지모드만 가지며, 바람직하지 못한 영에너지모드는 발견되지 않았다. 또한 조각보시험을 수행한 결과 이 연구의 쉘요소는 모든 조각보시험을 통과함을 확인하였다. 따라서 이 연구에서 개발된 쉘요소는 일정변형률상태를 표현할 수 있음을 확인하였으며, 요소의 크기가 작아짐에 따라 참값에 수렴할 것으로 기대된다.

 이 연구에서는 쉘요소의 수치모형화의 편의성을 위하여 두께방향 절점좌표계 벡터를 근사적으로 결정하였다. 즉, 감절점쉘요소의 상하면 절점좌표의 차이에 의해 두께방향 절점좌표계를 결정하는 것이 아니라, 곡선좌표계상의 기저벡터의 벡터곱에 의해 요소별로 독립적인 두께방향 절점벡터를 구하였다. 이러한 방법은 매우 작은 정도의 정확도의 저하를 가져오지만, 수치모형화의 편의성을 증진시킨다.

이상의 결과들을 살펴볼 때 개발된 요소는 전반적으로 양호한 결과를 보여주었으며, 실제문제에 효과적으로 사용할수 있는 유한요소라고 판단된다. 

감사의 글

이 연구는 지식경제부 기술혁신사업 중 원자력융합원천기술과제(2011T100200078)의 연구비 지원으로 수행하였습니다.