NINI
322 lines
23 KiB
Markdown
322 lines
23 KiB
Markdown
좌굴(static buckling)과 동적 좌굴(dynamic buckling)로 나눌 수 있으며, 이절에서는 정적 좌굴과 동적 좌굴에 대해 알아보도록 하겠다.
|
|
|
|
#### 2.3.1 Static Buckling
|
|
|
|
이처럼 좌굴은 다양한 힘과 형상, 조건 등에서 발생하며 이러한 좌굴 문제를 해석하기 위한 지배방정식은 가상일에 기반을 둔 운동 에너지(total potential energy) 식(2.78)과 위치 에너지(kinetic energy)에 Hamilton's principle을 적용하여 구할 수 있다.
|
|
|
|
$$
|
|
\Pi = U + V = \frac{1}{2} \mathbf{q}^T \mathbf{K} \mathbf{q} + \frac{1}{2} \mathbf{q}^T \mathbf{K}_g \mathbf{q} \tag{2.78}
|
|
$$
|
|
|
|
$$
|
|
T = \frac{1}{2} \dot{\mathbf{q}}^T \mathbf{M} \dot{\mathbf{q}} \tag{2.79}
|
|
$$
|
|
|
|
식(2.78)과 식(2.79)에 Euler-Lagrange 방정식 식(2.80)을 적용하면, 좌굴해석을 위한 지배 방정식 식(2.81)을 구할 수 있다.
|
|
|
|
$$
|
|
\frac{\partial L}{\partial \mathbf{q}} - \frac{d}{dt} \left( \frac{\partial L}{\partial \dot{\mathbf{q}}} \right) = 0 , \quad where \quad L = U + V - T \tag{2.80}
|
|
$$
|
|
|
|
$$
|
|
\mathbf{M}\{\ddot{\mathbf{q}}\} + \left(\mathbf{K} + \mathbf{K}_{g}\right)\left\{\mathbf{q}\right\} = \left\{\mathbf{0}\right\} \tag{2.81}
|
|
$$
|
|
|
|
여기서 M 은 질량 행렬(mass matrix), K 는 강성 행렬(stiffness matrix), K $_g$ 는 기하 강성 행렬(geometric stiffness matrix)를 의미한다. 지배 방정식으로부터 정적 좌굴 해석은 식(2.82)의 형태로 나타난다.
|
|
|
|
$$
|
|
\left(\mathbf{K} + \lambda \mathbf{K}_{g}\right)\left(\mathbf{q}\right) = \left\{\mathbf{0}\right\} \tag{2.82}
|
|
$$
|
|
|
|
이 식은 고유치 문제(eigenvalue problem)와 유사한 형태를 가지고 있으므로 고유치 해석 솔버를 사용하여 고유치 $\lambda$ , 와 고유치
|
|
|
|
벡터(eigenvector)를 구할 수 있다. 이 때 고유치 해석을 위해서 Block Lanczos 방법을 사용한 "BLZPACK"이라는 오픈 소스(open source)를 활용하였다. 일반적인 고유치 해석의 경우 기하 강성 행렬 $\mathbf{K}_g$ 의 자리에 질량 행렬 $\mathbf{M}$ 을 사용하며, 질량 행렬의 특성 상 집중 행렬(lumped matrix) 또는 일관 행렬(consistent matrix)을 이용하여 고유치를 계산한다. 하지만 정적 좌굴 해석의 경우 기하 강성 행렬 $\mathbf{K}_g$ 는 집중 행렬 형태를 구현하기 어려우므로 일관 행렬 형태로 계산을 수행해야 한다. 따라서 정적 좌굴 해석의 경우 식(2.82)에 대한 고유치 해석을 통해 식(2.83)과 같이 임계 좌굴 압력(critical buckling pressure) 또는 임계 좌굴 하중(critical buckling load)를 계산하게 된다. [6]
|
|
|
|
$$
|
|
P_{cr_i} = \lambda_i P \quad or \quad F_{cr_i} = \lambda_i F \tag{2.83}
|
|
$$
|
|
|
|
### 2.3.2 Dynamic Buckling
|
|
|
|
동적 좌굴 해석의 경우 정적 좌굴 해석과는 달리 식(2.84)와 같이 시간에 따라 변화하는 하중이 가해진다. 여기서 $P_0$ 는 정적 압축 하중(static compressive loading), $P_t$ 는 축 방향의 동적 하중(dynamic axial loading), $\theta$ 는 축 방향의 동적 하중에 대한 가진 주파수(excitation frequency), t는 시간을 의미한다.
|
|
|
|
$$
|
|
P(t) = P_0 + P_t \cos(\theta t) \tag{2.84}
|
|
$$
|
|
|
|
식(2.84)와 같이 시간에 따라 변화하는 하중이 가해질 때 운동 방정식(equation of motion)은 지배 방정식 식(2.81)로부터 식(2.85)와 같이 구할 수 있다.[6]
|
|
|
|
$$
|
|
\mathbf{M}\{\ddot{\mathbf{q}}\} + \left(\mathbf{K} + \mathbf{K}_{g}^{(s)} + \beta \cos(\theta t) \mathbf{K}_{g}^{(d)}\right) \left(\mathbf{q}\right) = \left\{\mathbf{0}\right\} \tag{2.85}
|
|
$$
|
|
|
|
여기에서 $\mathbf{K}$ 는 하중이 가해지지 않은 상태에서의 강성 행렬(stiffness matrix)이고, $\mathbf{K}_{g}^{(s)}$ 는 정적 압축 하중에 대한 기하 강성 행렬(geometric
|
|
|
|
stiffness matrix)이며, $\mathbf{K}_g^{(d)}$ 는 축 방향의 동적 하중에 대한 기하 강성행렬이다. $\mathbf{M}$ 은 질량 행렬(mass matrix)이고, $\boldsymbol{\beta}$ 는 동적 하중의 척도인자(dynamic load scale factor)를 의미한다. 이 때 주요 동적 불안정 경계 영역(principal region of dynamic instability)을 찾기 위해서 식(2.86)과 같이주기가 $2T(=2\times 2\pi/\theta)$ 인 해를 가정하였다. 여기서 $\mathbf{a}$ 와 $\mathbf{b}$ 는 임의의벡터를 의미한다.[6]
|
|
|
|
$$
|
|
\mathbf{q}(t) = \mathbf{a}\sin\frac{\theta t}{2} + \mathbf{b}\cos\frac{\theta t}{2} \tag{2.86}
|
|
$$
|
|
|
|
식(2.86)을 식(2.85)에 대입하여 푸리에 급수 전개(Fourier series expansion)를 하면 식(2.87)과 같이 나타낼 수 있다.
|
|
|
|
$$
|
|
\mathbf{a}\sin\frac{\theta t}{2} \left[ -\frac{\theta^2}{4}\mathbf{M} + \mathbf{K} + \mathbf{K}_g^{(s)} + \frac{\beta}{2}\mathbf{K}_g^{(d)} \right] + \mathbf{b}\cos\frac{\theta t}{2} \left[ -\frac{\theta^2}{4}\mathbf{M} + \mathbf{K} + \mathbf{K}_g^{(s)} - \frac{\beta}{2}\mathbf{K}_g^{(d)} \right] + \underbrace{\qquad \qquad }_{\text{high order term}} = 0 \tag{2.87}
|
|
$$
|
|
|
|
식(2.87)에서 고차항을 무시하면 식(2.88)과 같이 쓸 수 있다.
|
|
|
|
$$
|
|
\mathbf{a} \sin \frac{\theta t}{2} \left[ -\frac{\theta^2}{4} \mathbf{M} + \mathbf{K} + \mathbf{K}_g^{(s)} + \frac{\beta}{2} \mathbf{K}_g^{(d)} \right] + \mathbf{b} \cos \frac{\theta t}{2} \left[ -\frac{\theta^2}{4} \mathbf{M} + \mathbf{K} + \mathbf{K}_g^{(s)} - \frac{\beta}{2} \mathbf{K}_g^{(d)} \right] \approx 0 \tag{2.88}
|
|
$$
|
|
|
|
모든 시간에 대해 식(2.88)를 만족하는 자명하지 않은 해(non-trivial solution)를 얻기 위해서는 식(2.89)와 같은 고유치 문제(eigenvalue problem) 형태를 갖는다.[6]
|
|
|
|
$$
|
|
\left| \mathbf{K} + \left( \mathbf{K}_{g}^{(s)} \pm \frac{\beta}{2} \mathbf{K}_{g}^{(d)} \right) - \frac{\theta^{2}}{4} \mathbf{M} \right| = 0 \tag{2.89}
|
|
$$
|
|
|
|
동적 좌굴 해석의 경우 정적 좌굴 해석 문제와는 달리 질량 행렬을 사용하며, 본 논문에서는 집중 질량 행렬(lumped mass matrix)을 사용하여 고유치를 구했으며, 이 때 고유치 $\lambda_i = \theta_i^2/4$ 이 된다. 그리고 이로부터 구조물이 불안정해지는 가진 주파수(excitation frequency) $\theta_i$ 를 구할 수 있으며, 동적 하중 변화에 따른 불안정 경계 영역(instability region)을 구할 수 있다.
|
|
|
|
#### 2.3.3 Dynamic Buckling Theory of Beam
|
|
|
|
|
|
|
|
Fig. 6과 같이 보(beam) 구조물에 축 방향으로 정적 압축 하중(compressive static load)과 동적 압축 하중(compressive dynamic load)이 동시에 가해질 경우 보의 동적 좌굴 현상을 이론적으로 전개할 수 있다.
|
|
|
|
우선 보의 정적 휨(static bending)에 관한 식은 식(2.90)과 같다. [7]
|
|
|
|
$$
|
|
EI\frac{d^2v}{dx^2} + Pv = 0 (2.90)
|
|
$$
|
|
|
|
식(2.90)을 두 번 미분하면 식(2.91)과 같이 쓸 수 있다.
|
|
|
|
$$
|
|
EI\frac{d^4v}{dx^4} + P\frac{d^2v}{dx^2} = 0 {(2.91)}
|
|
$$
|
|
|
|
이 때 E는 탄성 계수(Young's Modulus)이고, I는 면적 관성 모멘트(area moment of inertia)이며, P는 축 방향의 힘을 의미한다. 그리고 Fig. 6과 같이 정적 압축 하중과 동적 압축 하중이 동시에 가해지는 경우 힘 P는 식(2.92)와 같이 나타낼 수 있다.
|
|
|
|
$$
|
|
P(t) = P_0 + P_t \cos(\theta t) \tag{2.92}
|
|
$$
|
|
|
|
여기서 $P_0$ 는 보에 가해지는 정적 압축 하중이고, $P_t$ 는 축 방향의 동적 하중에 대한 진폭(amplitude)를 의미하며, $\theta$ 는 축 방향의 동적 하중에 대한 가진 주파수(excitation frequency)를 의미한다. 이 때 동적 하중의 주기 T는 $2\pi/\theta$ 이다.
|
|
|
|
횡 방향 관성력만을 고려하여 운동 방정식을 구성하면, 식(2.91)은 식(2.93)과 같이 표현할 수 있다.
|
|
|
|
$$
|
|
EI\frac{\partial^4 v(x,t)}{\partial x^4} + \left(P_0 + P_t \cos(\theta t)\right) \frac{\partial^2 v(x,t)}{\partial x^2} + \rho \frac{\partial^2 v(x,t)}{\partial t^2} = 0 \tag{2.93}
|
|
$$
|
|
|
|
식(2.93)에서 v(x,t)는 시간에 따른 처짐을 의미하며, $\rho$ 는 단위 길이당 밀도를 의미한다. 위의 식(2.93)이 식(2.94)와 같은 형태의 해를 갖는다면 변수 분리(separation of variable)를 통해 식(2.93)은 식(2.95)와 같이 표현할수 있다. 이 때 식(2.94)는 경계 조건(boundary condition)을 만족한다.
|
|
|
|
$$
|
|
v(x,t) = \sum_{k=1}^{\infty} f_k(t) \sin \frac{k\pi x}{L} \tag{2.94}
|
|
$$
|
|
|
|
$$
|
|
\sum_{k=1}^{\infty} \left[ EI \frac{k^4 \pi^4}{L^4} f_k(t) - \left( P_0 + P_t \cos(\theta t) \right) \frac{k^2 \pi^2}{L^2} f_k(t) + \rho \frac{d^2 f_k(t)}{dt^2} \right] \sin \frac{k \pi x}{L} = 0 \tag{2.95}
|
|
$$
|
|
|
|
식(2.95)를 만족하기 위해서는 모든 k에 대해 $\sin$ 함수의 계수가 영(零)이 되어야 하며, 따라서 모든 k에 대해 식(2.96)이 만족해야 한다.
|
|
|
|
$$
|
|
EI\frac{k^{4}\pi^{4}}{L^{4}}f_{k}(t) - (P_{0} + P_{t}\cos(\theta t))\frac{k^{2}\pi^{2}}{L^{2}}f_{k}(t) + \rho\frac{d^{2}f_{k}(t)}{dt^{2}} = 0 \qquad (k = 1, 2, 3, \cdots) \tag{2.96}
|
|
$$
|
|
|
|
그리고 이를 자유진동 상태에서의 고유 주파수(natural frequency) $\omega_k$ 와 오일러 보(Euler beam)의 임계 좌굴 하중(critical buckling load) $P_k^{cr}$ 을 각각식(2.98), (2.99)와 같이 정의하고 이를 사용하여 식(2.96)를 다시 쓰면식(2.97)과 같이 나타낼 수 있다.
|
|
|
|
$$
|
|
\frac{d^2 f_k(t)}{dt^2} + \omega_k^2 \left( 1 - \frac{P_0 + P_t \cos(\theta t)}{P_k^{cr}} \right) f_k(t) = 0 \quad (k = 1, 2, 3, \dots) \tag{2.97}
|
|
$$
|
|
|
|
$$
|
|
\omega_k = \frac{k^2 \pi^2}{L^2} \sqrt{\frac{EI}{\rho}} \tag{2.98}
|
|
$$
|
|
|
|
$$
|
|
P_k^{cr} = \frac{k^2 \pi^2}{L^2} EI \tag{2.99}
|
|
$$
|
|
|
|
또한 일정한 압축 하중 $P_0$ 가 가해지는 보의 고유 주파수 $\Omega_k$ 와 가진 매개변수(excitation parameter) $\mu_k$ 를 식(2.101), (2.102)와 같이 정의하고, 이를 식(2.97)에 적용하면 식(2.100)과 같이 쓸 수 있다.
|
|
|
|
$$
|
|
\frac{d^2 f_k(t)}{dt^2} + \Omega_k^2 (1 - 2\mu_k \cos(\theta t)) f_k(t) = 0 \qquad (k = 1, 2, 3, \dots) \tag{2.100}
|
|
$$
|
|
|
|
$$
|
|
\Omega_k^2 \equiv \omega_k^2 \frac{\left(P_k^{cr} - P_0\right)}{P_k^{cr}} \tag{2.101}
|
|
$$
|
|
|
|
$$
|
|
\mu_k = \frac{P_t}{2(P_k^{cr} - P_0)} \tag{2.102}
|
|
$$
|
|
|
|
그리고 여기서 k=1일 때의 기본모드(fundamental mode)만 고려하면 4(2.100)은 4(2.103)과 같이 나타낼 수 있고, 마찬가지로 4(2.101)과 4(2.102)도 4(2.104)와 4(2.105)로 나타낼 수 있다.
|
|
|
|
$$
|
|
\frac{d^2 f(t)}{dt^2} + \Omega^2 (1 - 2\mu \cos(\theta t)) f(t) = 0 \tag{2.103}
|
|
$$
|
|
|
|
$$
|
|
\Omega^{2} = \left(\frac{\pi^{2}}{L^{2}} \sqrt{\frac{EI}{\rho}}\right)^{2} \left[1 - P_{0} \frac{L^{2}}{\pi^{2} EI}\right] \tag{2.104}
|
|
$$
|
|
|
|
$$
|
|
\mu = \frac{P_t}{2[(\pi^2 EI)/L^2 - P_0]} \tag{2.105}
|
|
$$
|
|
|
|
이 때 식(2.103)은 식(2.106)과 같은 Mathieu-Hill 방정식의 형태를 갖게 된다.[8]
|
|
|
|
$$
|
|
f'' + \Omega^2 [1 - 2\mu \Phi(t)] f = 0 \tag{2.106}
|
|
$$
|
|
|
|
여기서 $\Phi(t)$ 는 가진 주기가 $T=2\pi/\theta$ 인 임의의 주기 함수이고, 다음과 같이 매개변수의 변화에 따라 해가 안정하거나 불안정해질 수 있다.
|
|
|
|
- ⇒ 해의 주기가 T(또는 2T)에서 주기 2T(또는 T)로 바뀔 경우 : 안정
|
|
- ⇒ 해의 주기가 T(또는 2T)에서 주기 T(또는 2T)로 바뀔 경우 : 불안정
|
|
|
|
즉 주기 2T, T인 해가 존재하므로 이를 푸리에 급수(Fourier series)로 가정하고, 이를 통해 매개변수에 따른 해의 안정성을 판별할 수 있다.
|
|
|
|
불안정 영역(instability region)을 구하기 위해 주기가 $2T(=2\times 2\pi/\theta)$ 인해를 식(2.107)과 같이 정의하면 다음과 같다.
|
|
|
|
$$
|
|
f(t) = \sum_{k=1,3,5,\dots}^{\infty} \left( a_k \sin \frac{k\theta t}{2} + b_k \cos \frac{k\theta t}{2} \right) \tag{2.107}
|
|
$$
|
|
|
|
이 중에서 k=1인 경우만 고려하면, 식(2.108)과 같다.
|
|
|
|
$$
|
|
f(t) = a\sin\frac{\theta t}{2} + b\cos\frac{\theta t}{2} \tag{2.108}
|
|
$$
|
|
|
|
이를 식(2.103)에 대입하면 식(2.109)와 같으며, 이 때 주기 2T에 대해 푸리에 급수 전개를 하면 식(2.110)과 같이 나타낼 수 있다.
|
|
|
|
$$
|
|
-\frac{\theta^{2}}{4}a\sin\frac{\theta t}{2} - \frac{\theta^{2}}{4}b\cos\frac{\theta t}{2} + \Omega^{2}(1 - 2\mu\cos(\theta t))a\sin\frac{\theta t}{2} + \Omega^{2}(1 - 2\mu\cos(\theta t))b\cos\frac{\theta t}{2} = 0
|
|
a\sin\frac{\theta t}{2}\left(-\frac{\theta^{2}}{4} + \Omega^{2} + \mu\Omega^{2}\right) + b\cos\frac{\theta t}{2}\left(-\frac{\theta^{2}}{4} + \Omega^{2} - \mu\Omega^{2}\right) + \underbrace{\dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots
|
|
$$
|
|
|
|
그리고 식(2.110)에서 주기가 2T인 항을 취하고 나머지 고차항을 무시하면, 식(2.111)과 같이 나타낼 수 있고, 이 때 모든 시간에 대해 식(2.111)을 만족하는 자명하지 않은 해(non-trivial solution)을 얻기 위해서는 식(2.112)를 만족해야 한다.
|
|
|
|
$$
|
|
a\sin\frac{\theta t}{2}\left(-\frac{\theta^2}{4} + \Omega^2 + \mu\Omega^2\right) + b\cos\frac{\theta t}{2}\left(-\frac{\theta^2}{4} + \Omega^2 - \mu\Omega^2\right) \approx 0 \tag{2.111}
|
|
$$
|
|
|
|
$$
|
|
\Omega^2 \left[ 1 \pm \mu - \frac{1}{4} \left( \frac{\theta}{\Omega} \right)^2 \right] = 0 \tag{2.112}
|
|
$$
|
|
|
|
이는 불안정 영역을 결정짓는 매개변수와 가진 진동수의 근사화된 식으로 식(2.113)과 같이 무차원 가진 진동수 $\theta/(2\Omega)$ 를 매개변수 $\mu$ 에 대한 양함수의 형태로 표현할 수 있으며, 이를 도시하여 가진 매개 변수 $\mu$ 에 따른 불안정 영역의 변화를 살펴 볼 수 있다.
|
|
|
|
$$
|
|
\frac{\theta}{2\Omega} = \sqrt{1 \pm \mu} \tag{2.113}
|
|
$$
|
|
|
|
# 3. Numerical Example
|
|
|
|
본 연구에서 개발한 동적 좌굴 유한요소해석 프로그램을 이용하여 동적 좌굴 해석을 하기 위한 기본적인 해석을 수행하였다. 기본적으로 선형 정적 해석, 비선형 정적 해석 그리고 정적 좌굴 해석을 수행하였으며, 최종적으로 본 연구의 목적인 동적 좌굴 해석을 수행하여 이론값 또는 실험값과의 비교를 통해 동적 좌굴 유한요소해석 프로그램의 타당성과 신뢰성을 검증하였다.
|
|
|
|
## 3.1 Linear Static Analysis
|
|
|
|
|
|
|
|
|
|
|
|
패치 테스트 모델(patch test model)은 Fig. 1과 같다. 이 모델에 대해서 constant curvature, constant shear, constant twist 패치 테스트를 수행하였으며, 재료의 탄성계수(Young"s modulus)는 2.1E+06이고, 푸아송의 비(Poisson"s ratio)는 0.3이다. 그리고 모델의 두께는 1.0과 0.001 두 모델에 대해 패치 테스트를 수행하였다.
|
|
|
|
#### 3.1.1.1 Constant Curvature Patch Test
|
|
|
|
|
|
|
|
Constant curvature 패치 테스트 모델에 대한 경계조건은 우선 1번 절점에 대해 1,3 방향의 변위와 2 방향의 회전을 구속하였고, 5번 절점에 대해 1,2,3 방향의 변위와 2 방향의 회전을 구속하였다. 또한 외력은 3번, 7번 절점에 2 방향의 모멘트를 가하였으며, 이 때 힘의 크기는 두께 1.0일 때 M=1000, 두께 0.001일 때 M=0.0001을 가하였다. 이러한 경계조건과 외력은 Fig. 8로부터 확인할 수 있으며, constant curvature 패치테스트 해석 결과 각 절점의 곡률 $(\rho)$ 는 식(3.1)과 같이 구할 수 있다.
|
|
|
|
$$
|
|
\frac{1}{\rho} = \frac{d\theta}{dx} = \frac{d^2w}{dx^2} \tag{3.1}
|
|
$$
|
|
|
|
[Table HTML](tables/chunk-003-table-001.html)
|
|
|
|
Table. 1은 constant curvature 패치 테스트 결과로 두 가지 두께에 대해 구속조건이 적용된 절점을 제외한 나머지 절점에서 모두 일정한 곡률이 나오는 것을 볼 수 있다.
|
|
|
|
# 3.1.1.2 Constant Shear Patch Test
|
|
|
|
|
|
|
|
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와 같다.
|
|
|
|
[Table HTML](tables/chunk-003-table-002.html)
|
|
|
|
Constant shear 패치 테스트 결과로 두 가지 두께에 대해 구속조건이 적용된 절점을 제외한 나머지 절점에서 모두 일정한 기울기가 나타나는 것을 확인할 수 있으며, 실제로 3 방향 변위는 식(3.2)를 사용하여 계산할 수 있다.
|
|
|
|
$$
|
|
\tau = \frac{V}{A} = G\gamma, \quad \gamma = \frac{w}{L}, \quad G = \frac{E}{2(1+v)} \quad \Rightarrow \quad w = \frac{2(1+v)VL}{EA} \tag{3.2}
|
|
$$
|
|
|
|
#### 3.1.1.3 Constant Twist Patch Test
|
|
|
|
|
|
|
|
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과 같다.
|
|
|
|
[Table HTML](tables/chunk-003-table-003.html)
|
|
|
|
Constant twisting 테스트 결과를 보면 두께가 0.001일 때 x와 y방향에 대해 일정한 곡률이 나오는 것을 볼 수 있지만 두께가 1.0일 때는 그렇지 않음을 알 수 있다. 이 문제는 두께가 두꺼울 때 횡 전단 변형(transverse shear deformation)이 지배적으로 나타나기 때문에 발생하는 것으로 전단 보정 계수(shear correction factor)를 사용하여 전단 변형(shear deformation)을 억제하거나 사각형 요소를 사용함으로써 이 문제를 해결할 수 있다.[1]
|
|
|
|
### 3.1.2 Pinched Cylinder
|
|
|
|
|
|
|
|
Fig. 11 Pinched Cylinder Model 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방향의 변위에 대해 해석 결과를 도시한 그림으로 변형 형상이 유사하게 나옴을 확인할 수 있다.
|
|
|
|
[Table HTML](tables/chunk-003-table-004.html)
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
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와 현재 코드의 해석 결과를 직접 도시한 그림으로 서로 유사한 결과가 나옴을 알 수 있다.
|
|
|
|
[Table HTML](tables/chunk-003-table-005.html)
|
|
|
|
|
|
|
|
|
|
|
|
# 3.2 Geometric Nonlinear Analysis
|
|
|
|
|
|
|
|
Fig. 18과 같은 보(beam) 형상에 대해 기하비선형 정적 해석(geometric nonlinear static analysis)을 수행하여 ABAQUS와 그 결과를 비교해 보았다. 먼저 가로 12m, 세로 1m, 두께 0.01m 이고, 요소는 총 12개를 사용하였다. 탄성계수는 1.0MPa이고, 푸아송 비는 0.0의 물성치를 주었다. 경계조건은 한 면은 XYZ방향의 변위와 회전을 모두 구속하였고, 다른 한 면은 Y방향으로 모멘트를 가하였다. 그리고 물체가 Y방향에 대해 회전할 때 X방향의 회전으로 인한 형상의 뒤틀림을 방지하기 위해 모든 절점에
|
|
|
|
대해 X방향 회전을 구속하였다. 해석 결과 Fig. 19로부터 ABAQUS 해석 결과와 유사한 결과가 나옴을 확인 할 수 있으며, 이 때 힘의 증가에 따른 형상 변화는 Fig. 20과 같다.
|
|
|
|
|
|
|
|
|
|
|
|
## 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의 압축력을 가하였다.
|
|
|
|
|
|
|
|
해석 결과 나오는 고유치(eigenvalue) i 는 Table. 6과 같고, ABAQUS와 유사한 결과가 나옴을 확인할 수 있다. 또한 고유치 벡터(eigenvector)로 부터 확인할 수 있는 좌굴 형상 역시 ABAQUS와 유사함을 확인할 수 있다.(Table. 7)
|