두께	NODE	x-coord.	$\beta$	$\rho$
1.0	2	2.0	0.0011429	1750
	3	10.0	0.0057143	1750
	4	8.0	0.0045714	1750
	6	4.0	0.0022857	1750
	7	10.0	0.0057143	1750
	8	8.0	0.0045714	1750
0.001	2	2.0	0.114286	17.50
	3	10.0	0.571429	17.50
	4	8.0	0.457143	17.50
	6	4.0	0.228571	17.50
	7	10.0	0.571429	17.50
	8	8.0	0.457143	17.50

Table. 1 Constant Curvature Patch Test Results Table. 1은 constant curvature 패치 테스트 결과로 두 가지 두께에 대해구속조건이 적용된 절점을 제외한 나머지 절점에서 모두 일정한 곡률이나오는 것을 볼 수 있다. # 3.1.1.2 Constant Shear Patch Test 

text_image

U₃ = 0 SHEAR U₁₋₂ = 0 β = 0 U₃ = 0 Q

Fig. 9 Constant Shear Patch Test Model Constant shear 패치 테스트 모델에 대한 경계조건은 우선 1번과 5번절점에 대해 3 방향의 변위를 구속하였고, 모든 절점에 대해 1,2 방향의변위와 2 방향의 회전을 구속하였다. 또한 외력은 3번, 7번 절점에 3방향의 집중 하중을 가하였으며, 이 때 힘의 크기는 두께 1.0일 때Q=2000, 두께 0.001일 때 Q=2.0을 가하였다. 이러한 경계 조건과 외력은Fig. 9로부터 확인할 수 있으며, constant shear 패치 테스트 해석 결과는Table. 2와 같다.

두께	NODE	x-coord.	z-disp. (w)	dw/dx
1.0	2	2.0	0.000495238	2.47619E-04
	3	10.0	0.00247619	2.47619E-04
	4	8.0	0.00198095	2.47619E-04
	6	4.0	0.000990476	2.47619E-04
	7	10.0	0.00247619	2.47619E-04
	8	8.0	0.00198095	2.47619E-04
0.001	2	2.0	0.000495238	2.47619E-04
	3	10.0	0.00247619	2.47619E-04
	4	8.0	0.00198095	2.47619E-04
	5	4.0	0.000990476	2.47619E-04
	6	10.0	0.00247619	2.47619E-04
	7	10.0	0.00247619	2.47619E-04
	8	8.0	0.00198095	2.47619E-04

Table. 2 Constant Shear Patch Test Results Constant shear 패치 테스트 결과로 두 가지 두께에 대해 구속조건이적용된 절점을 제외한 나머지 절점에서 모두 일정한 기울기가 나타나는것을 확인할 수 있으며, 실제로 3 방향 변위는 식(3.2)를 사용하여 계산할수 있다. $$ \tau = \frac {V}{A} = G \gamma , \quad \gamma = \frac {w}{L}, \quad G = \frac {E}{2 (1 + \nu)} \quad \Rightarrow \quad w = \frac {2 (1 + \nu) V L}{E A} \tag {3.2} $$ # 3.1.1.3 Constant Twist Patch Test 

text_image

U₃ = 0 P TWISTING U₁₋₂ = 0 U₃ = 0 U₃ = 0

Fig. 10 Constant Twist Patch Test Model Constant twist 패치 테스트 모델에 대한 경계조건은 우선 1, 3, 5번절점에 대해 3 방향의 변위를 구속하였고, 모든 절점에 대해 1,2 방향의변위를 구속하였다. 또한 외력은 7번 절점에 3 방향의 집중 하중을가하였으며, 이 때 힘의 크기는 두께 1.0일 때 P=1000, 두께 0.001일 때P=1.0E-06을 가하였다. 이러한 경계 조건과 외력은 Fig. 10으로부터확인할 수 있으며, constant twist 패치 테스트 해석 결과는 Table. 3과 같다.

두께	NODE	x-coord.	$\alpha$	$\rho_1$	y-coord.	$\beta$	$\rho_2$
1.0	2	2.0	0.00751761	266.04	2.0	0.0075345	265.45
	4	8.0	0.03103950	257.74	3.0	0.0114079	262.98
	6	4.0	0.01518240	263.46	7.0	0.0272324	257.05
	7	10.0	0.03767130	265.45	10.0	0.0379748	263.33
	8	8.0	0.03104290	257.71	7.0	0.0272190	257.17
0.001	2	2.0	0.00742857	269.23	2.0	0.00742857	269.23
	4	8.0	0.02971430	269.23	3.0	0.0111429	269.23
	6	4.0	0.01485710	269.23	7.0	0.0260000	269.23
	7	10.0	0.03714290	269.23	10.0	0.0371429	269.23
	8	8.0	0.02971430	269.23	7.0	0.0260000	269.23

Table. 3 Constant Twist Patch Test Results Constant twisting 테스트 결과를 보면 두께가 0.001일 때 x와 y방향에대해 일정한 곡률이 나오는 것을 볼 수 있지만 두께가 1.0일 때는 그렇지않음을 알 수 있다. 이 문제는 두께가 두꺼울 때 횡 전단 변형(transverseshear deformation)이 지배적으로 나타나기 때문에 발생하는 것으로 전단보정 계수(shear correction factor)를 사용하여 전단 변형(shear deformation)을억제하거나 사각형 요소를 사용함으로써 이 문제를 해결할 수 있다.[1] # 3.1.2 Pinched Cylinder 

text_image

P D C R A B End diaphragm End diaphragm P

Fig. 11 Pinched Cylinder Model 

text_image

P/4 Sym. Sym. BC(12) Sym. y z x

Fig. 12 Pinched Cylinder 1/8 Model Fig. 11과 같이 양 끝에 막이 있는 pinched cylinder 쉘은 이론적인 해를비교할 수 있기 때문에 쉘 검증 문제에 적합하다. 이 때 하중은중간면에서 서로 반대 방향으로 집중 하중 P가 가해지고 있으며, 형상 및하중이 대칭이기 때문에 Fig. 12와 같이 1/8 모델을 사용하였다. Pinched cylinder 쉘의 형상은 길이(L) 600, 반경(R) 300, 두께(t) 3이고,재료 특성은 탄성 계수(Young‟s modulus) 3.0E+06, 푸아송 비(Poisson‟s ratio)0.3이며, 외력 P의 크기는 1이다. 격자(mesh)는 20x20, 30x30, 40x40 격자를사용하였으며, 경계 조건은 절단면은 각각 대칭 경계 조건을적용하였으며, pinched cylinder 쉘의 end diaphragm 부분은 X와 Y 방향의변위를 구속하였다. 해석 결과 Table. 4 로부터 ABAQUS와 현재 쉘 모두 이론적인 해와유사하게 나옴을 확인할 수 있으며, ABAQUS결과가 보다 이론적인 해에 가까움을 확인할 수 있다. 이는 쉘 요소가 다르기 때문에 발생하는차이로 ABAQUS 쉘 요소가 현재 쉘 요소에 비해 보다 유연함을 알 수있다. 그리고 Fig. 13으로부터 요소의 수가 증가함에 따라 ABAQUS와현재 쉘의 결과가 모두 이론적인 해에 수렴함을 알 수 있다. 또한 Fig.14는 Y방향의 변위에 대해 해석 결과를 도시한 그림으로 변형 형상이유사하게 나옴을 확인할 수 있다.

Mesh	Exact	Present	$W_{present}$ / $W_{exact}$	ABAQUS	$W_{abaqus}$ / $W_{exact}$
20x20	1.8248E-05	1.74362E-05	0.9555	1.77866E-05	0.9747
30x30		1.79593E-05	0.9842	1.81676E-05	0.9956
40x40		1.81878E-05	0.9967	1.82150E-05	0.9982

Table. 4 Comparison of Linear Static Analysis for Pinched Cylinder with Exact Solution 

line

| Number of elements per side | Present | ABAQUS | | --------------------------- | ------- | ------ | | 20 | 0.955 | 0.975 | | 30 | 0.985 | 0.995 | | 40 | 0.998 | 0.999 |

Fig. 13 Comparison of Convergence for Pinched Cylinder with ABAQUS  Fig. 14 Comparison of Linear Static Analysis for Pinched Cylinder with ABAQUS # 3.1.3 Hemispherical Shell 

text_image

Sym. Sym. P P x y z

Fig. 15 Hemispherical Shell Model Fig. 15와 같은 반구형(hemispherical) 쉘의 1/4 모델을 사용하여이론적인 해와 ABAQUS 결과값에 대해 비교해 보았다. 먼저 반경은 10m,두께는 0.04m이고, 요소는 한 면당 9개, 17개 격자(mesh)를 사용하였다.탄성계수는 68.25MPa이고, 푸아송 비는 0.3의 물성치를 주었다.경계조건은 좌우 면에 대해 대칭 경계조건을 적용하였고, 하중은 Fig.15와 같은 지점에 (+)Z, (-)X방향으로 각각 1씩 가하였다. 해석 결과 Table. 5로부터 ABAQUS는 이론값보다 큰 값이 나오고, 현재 코드 결과는이론값보다 작은 값이 나옴을 알 수 있는데 이는 사용한 쉘 요소가다르므로 요소의 특성에 따른 차이로 볼 수 있다. 하지만 Fig. 16로부터두 결과값 모두 요소 수가 증가함에 따라 이론값에 수렴함을 확인할 수있으며, Fig. 17은 Y방향 변위에 대해 ABAQUS와 현재 코드의 해석결과를 직접 도시한 그림으로 서로 유사한 결과가 나옴을 알 수 있다.

Node/side	Exact	Present	$W_{present}$ / $W_{exact}$	ABAQUS	$W_{abaqus}$ / $W_{exact}$
9	0.0924	0.0888439	0.9615	0.0941153	1.0186
17		0.0919091	0.9947	0.0933083	1.0098
25		0.0921386	0.9972	0.0928799	1.0052

Table. 5 Comparison of Linear Static Analysis for Hemispherical Shell with Exact Solution 

line

| Number of nodes per side | Present | ABAQUS | | ------------------------ | ------- | ------ | | 9 | 0.96 | 1.02 | | 17 | 0.995 | 1.01 | | 25 | 0.998 | 1.005 |

Fig. 16 Comparison of Convergence for Hemispherical Shell with ABAQUS  Fig. 17 Comparison of Linear Static Analysis for Hemisphrical Shell with ABAQUS # 3.2 Geometric Nonlinear Analysis 

text_image

BC(Fixed) y x Moment

Fig. 18 Beam Model for Geometric Nonlinear Analysis Fig. 18과 같은 보(beam) 형상에 대해 기하비선형 정적 해석(geometricnonlinear static analysis)을 수행하여 ABAQUS와 그 결과를 비교해 보았다.먼저 가로 12m, 세로 1m, 두께 0.01m 이고, 요소는 총 12개를사용하였다. 탄성계수는 1.0MPa이고, 푸아송 비는 0.0의 물성치를 주었다.경계조건은 한 면은 XYZ방향의 변위와 회전을 모두 구속하였고, 다른 한면은 Y방향으로 모멘트를 가하였다. 그리고 물체가 Y방향에 대해 회전할때 X방향의 회전으로 인한 형상의 뒤틀림을 방지하기 위해 모든 절점에 대해 X방향 회전을 구속하였다. 해석 결과 Fig. 19로부터 ABAQUS 해석결과와 유사한 결과가 나옴을 확인 할 수 있으며, 이 때 힘의 증가에따른 형상 변화는 Fig. 20과 같다. 

line

| Displacement & Rotation | Present_Disp.X | Present_Disp.Z | Present_Rota.Y | Abaqus_Disp.X | Abaqus_Disp.Z | Abaqus_Rota.Y | | ----------------------- | -------------- | -------------- | -------------- | ------------- | ------------- | ------------- | | 0 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | | 2 | 0.0070 | 0.0020 | 0.0100 | 0.0080 | 0.0040 | 0.0060 | | 4 | 0.0100 | 0.0050 | 0.0200 | 0.0120 | 0.0080 | 0.0140 | | 6 | 0.0150 | 0.0080 | 0.0300 | 0.0160 | 0.0120 | 0.0220 | | 8 | 0.0180 | 0.0120 | 0.0350 | 0.0200 | 0.0160 | 0.0280 | | 10 | 0.0200 | 0.0160 | 0.0380 | 0.0240 | 0.0200 | 0.0320 | | 12 | 0.0220 | 0.0200 | 0.0400 | 0.0280 | 0.0240 | 0.0360 | | 14 | 0.0250 | 0.0240 | 0.0420 | 0.0320 | 0.0280 | 0.0400 | | 16 | 0.0280 | 0.0280 | 0.0440 | 0.0360 | 0.0320 | 0.0440 |

Fig. 19 Comparison of Geometric Nonlinear Analysis for Beam with ABAQUS 

line

| M | Value | |-------|-------| | 0.04 | 1.0 | | 0.03 | 0.8 | | 0.02 | 0.6 | | 0.01 | 0.4 |

Fig. 20 Geometry Change of Beam According to Loads Increase # 3.3 Static Buckling Analysis # 3.3.1 Rectangular Plate Shell 정적 좌굴 해석 프로그램의 검증을 위해 먼저 Fig. 21과 같은 직사각형평판 쉘(rectangular plate shell) 형상에 대해 정적 좌굴 해석을 수행하였다.크기는 가로 20m, 세로 8m이고, 두께는 0.01m로 4노드 쉘 40x16 격자로모델링 하였으며, 재료의 물성치는 탄성계수 29MPa, 푸아송 비는 0.3을적용하였다. 경계 조건은 한 면은 XYZ방향의 변위를 구속하였고, 반대쪽면은 YZ방향의 변위를 구속하였다. 그리고 나머지 두 면은 Z방향의변위를 구속하였다. 마지막으로 하중은 YZ방향의 변위를 구속한 면에X방향으로 총 8N의 압축력을 가하였다. 

text_image

BC(23) F BC(3) BC(3) BC(123) Z X Y

Fig. 21 Rectangular Plate Shell Model 해석 결과 나오는 고유치(eigenvalue) $\lambda _ { i }$ 는 Table. 6과 같고, ABAQUS와유사한 결과가 나옴을 확인할 수 있다. 또한 고유치 벡터(eigenvector)로부터 확인할 수 있는 좌굴 형상 역시 ABAQUS와 유사함을 확인할 수있다.(Table. 7)