김경종
17 KiB
좌굴(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}{d t} \left(\frac {\partial L}{\partial \dot {\mathbf {q}}}\right) = 0, \quad \text { where } L = U + V - T \tag {2.80}
\mathbf {M} \{\ddot {\mathbf {q}} \} + \left(\mathbf {K} + \mathbf {K} _ {g}\right) \{\mathbf {q} \} = \{\mathbf {0} \} \tag {2.81}
여기서 M은 질량 행렬(mass matrix), K는 강성 행렬(stiffness matrix), K_{g} 는 기하 강성 행렬(geometric stiffness matrix)를 의미한다. 지배 방정식으로부터 정적 좌굴 해석은 식(2.82)의 형태로 나타난다.
\left(\mathbf {K} + \lambda \mathbf {K} _ {g}\right) \{\mathbf {q} \} = \{\mathbf {0} \} \tag {2.82}
이 식은 고유치 문제(eigenvalue problem)와 유사한 형태를 가지고 있으므로 고유치 해석 솔버를 사용하여 고유치 \lambda_{i} 와 고유치
벡터(eigenvector)를 구할 수 있다. 이 때 고유치 해석을 위해서 BlockLanczos 방법을 사용한 “BLZPACK”이라는 오픈 소스(open source)를활용하였다. 일반적인 고유치 해석의 경우 기하 강성 행렬 K 의 자리에질량 행렬 M을 사용하며, 질량 행렬의 특성 상 집중 행렬(lumped matrix)또는 일관 행렬(consistent matrix)을 이용하여 고유치를 계산한다. 하지만정적 좌굴 해석의 경우 기하 강성 행렬 K 는 집중 행렬 형태를 구현하기어려우므로 일관 행렬 형태로 계산을 수행해야 한다. 따라서 정적 좌굴해석의 경우 식(2.82)에 대한 고유치 해석을 통해 식(2.83)과 같이 임계좌굴 압력(critical buckling pressure) 또는 임계 좌굴 하중(critical bucklingload)를 계산하게 된다. [6]
P _ {c r _ {i}} = \lambda_ {i} P \quad \text {or} \quad F _ {c r _ {i}} = \lambda_ {i} F \tag {2.83}
2.3.2 Dynamic Buckling
동적 좌굴 해석의 경우 정적 좌굴 해석과는 달리 식(2.84)와 같이시간에 따라 변화하는 하중이 가해진다. 여기서 P 는 정적 압축하중(static compressive loading), P 는 축 방향의 동적 하중(dynamic axialloading), 는 축 방향의 동적 하중에 대한 가진 주파수(excitationfrequency), 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) \{\mathbf {q} \} = \{\mathbf {0} \} \tag {2.85}
여기에서 K 는 하중이 가해지지 않은 상태에서의 강성 행렬(stiffnessmatrix)이고, (s) K 는 정적 압축 하중에 대한 기하 강성 행렬(geometric
stiffness matrix)이며, \mathbf{K}_{g}^{(d)} 는 축 방향의 동적 하중에 대한 기하 강성 행렬이다. M 은 질량 행렬(mass matrix)이고, \beta 는 동적 하중의 척도 인자(dynamic load scale factor)를 의미한다. 이 때 주요 동적 불안정 경계 영역(principal region of dynamic instability)을 찾기 위해서 식(2.86)과 같이 주기가 2T(=2\times2\pi/\theta) 인 해를 가정하였다. 여기서 a 와 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)과 같이 나타낼 수 있다.
\begin{array}{l} \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 {\dots \dots} _ {\text { high order term }} = 0 \tag {2.87} \\ \end{array}
식(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 \ \textdegree ] 된다. 그리고 이로부터구조물이 불안정해지는 가진 주파수(excitation frequency) \theta _ { i } 를 구할 수있으며, 동적 하중 변화에 따른 불안정 경계 영역(instability region)을 구할수 있다.
2.3.3 Dynamic Buckling Theory of Beam
text_image
P₀ + Pₜ cos(θt) L v(x,t) 1954 y x
Fig. 6 Dynamic Buckling Model of Beam
Fig. 6과 같이 보(beam) 구조물에 축 방향으로 정적 압축하중(compressive static load)과 동적 압축 하중(compressive dynamic load)이동시에 가해질 경우 보의 동적 좌굴 현상을 이론적으로 전개할 수 있다.
우선 보의 정적 휨(static bending)에 관한 식은 식(2.90)과 같다. [7]
E I \frac {d ^ {2} v}{d x ^ {2}} + P v = 0 \tag {2.90}
식(2.90)을 두 번 미분하면 식(2.91)과 같이 쓸 수 있다.
E I \frac {d ^ {4} v}{d x ^ {4}} + P \frac {d ^ {2} v}{d x ^ {2}} = 0 \tag {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)과 같이 표현할 수 있다.
E I \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[ E I \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)}{d t ^ {2}} \right] \sin \frac {k \pi x}{L} = 0 \tag {2.95}
식(2.95)를 만족하기 위해서는 모든 k에 대해 sin 함수의 계수가 영(雫)이 되어야 하며, 따라서 모든 k에 대해 식(2.96)이 만족해야 한다.
E I \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)}{d t ^ {2}} = 0 \quad (k = 1, 2, 3, \dots) \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)}{d t ^ {2}} + \omega_ {k} ^ {2} \left(1 - \frac {P _ {0} + P _ {t} \cos (\theta t)}{P _ {k} ^ {c r}}\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 {E I}{\rho}} \tag {2.98}
P _ {k} ^ {c r} = \frac {k ^ {2} \pi^ {2}}{L ^ {2}} E I \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)}{d t ^ {2}} + \Omega_ {k} ^ {2} \left(1 - 2 \mu_ {k} \cos (\theta t)\right) f _ {k} (t) = 0 \quad (k = 1, 2, 3, \dots) \tag {2.100}
\Omega_ {k} ^ {2} \equiv \omega_ {k} ^ {2} \frac {\left(P _ {k} ^ {c r} - P _ {0}\right)}{P _ {k} ^ {c r}} \tag {2.101}
\mu_ {k} \equiv \frac {P _ {t}}{2 \left(P _ {k} ^ {c r} - P _ {0}\right)} \tag {2.102}
그리고 여기서 k=1일 때의 기본모드(fundamental mode)만 고려하면 식(2.100)은 식(2.103)과 같이 나타낼 수 있고, 마찬가지로 식(2.101)과 식(2.102)도 식(2.104)와 식(2.105)로 나타낼 수 있다.
\frac {d ^ {2} f (t)}{d t ^ {2}} + \Omega^ {2} \big (1 - 2 \mu \cos (\theta t) \big) f (t) = 0 \tag {2.103}
\Omega^ {2} = \left(\frac {\pi^ {2}}{L ^ {2}} \sqrt {\frac {E I}{\rho}}\right) ^ {2} \left[ 1 - P _ {0} \frac {L ^ {2}}{\pi^ {2} E I} \right] \tag {2.104}
\mu = \frac {P _ {t}}{2 \left[ \left(\pi^ {2} E I\right) / L ^ {2} - P _ {0} \right]} \tag {2.105}
이 때 식(2.103)은 식(2.106)과 같은 Mathieu-Hill 방정식의 형태를 갖게 된다. [8]
f ^ {\prime \prime} + \Omega^ {2} [ 1 - 2 \mu \Phi (t) ] f = 0 \tag {2.106}
여기서 \Phi(t) 는 가진 주기가 T=2\pi/\theta 인 임의의 주기 함수이고, 다음과 같이 매개변수의 변화에 따라 해가 안정하거나 불안정해질 수 있다.
⇒ 해의 주기가 T(또는 2T)에서 주기 2T(또는 T)로 바 Parks 경우 : 안정
⇒ 해의 주기가 T(또는 2T)에서 주기 T(또는 2T)로 바 Parks 경우 : 불안정
즉 주기 2T, T인 해가 존재하므로 이를 푸리에 급수(Fourier series)로 가정하고, 이를 통해 매개변수에 따른 해의 안정성을 판별할 수 있다.
불안정 영역(instability region)을 구하기 위해 주기가 2T(=2\times2\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)과 같이 나타낼 수 있다.
\begin{array}{l} - \frac {\theta^ {2}}{4} a \sin \frac {\theta t}{2} - \frac {\theta^ {2}}{4} b \cos \frac {\theta t}{2} \tag {2.109} \\ + \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 \\ \end{array}
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} _ {\text { high order term }} = 0 \tag {2.110}
그리고 식(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
3.1.1 Patch Test
scatter
| Point | X1 | X2 | Label |
|---|---|---|---|
| 1 | 0,0 | 0 | 1 (0,0) |
| 2 | 2,2 | 0 | 2 (2,2) |
| 3 | 10,0 | 0 | 3 (10,0) |
| 4 | 8,3 | 0 | 4 (8,3) |
| 6 | 4,7 | 0 | 6 (4,7) |
| 7 | 10,10 | 0 | 7 (10,10) |
| 8 | 8,7 | 0 | 8 (8,7) |
| The label 'INHE UNIVERSITY' is not present in the image. The data points are explicitly labeled with numbers and coordinates. |
Fig. 7 Patch Test Mesh
패치 테스트 모델(patch test model)은 Fig. 1과 같다. 이 모델에 대해서constant curvature, constant shear, constant twist 패치 테스트를 수행하였으며,재료의 탄성계수(Young‟s modulus)는 2.1E+06이고, 푸아송의 비(Poisson‟sratio)는 0.3이다. 그리고 모델의 두께는 1.0과 0.001 두 모델에 대해 패치테스트를 수행하였다.
3.1.1.1 Constant Curvature Patch Test
text_image
U_{1-2-3} = 0 β = 0 BENDING U_{1-3} = 0 β = 0 M M
Fig. 8 Constant Curvature Patch Test Model
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 패치테스트 해석 결과 각 절점의 곡률는 식(3.1)과 같이 구할 수 있다.
\frac {1}{\rho} = \frac {d \theta}{d x} = \frac {d ^ {2} w}{d x ^ {2}} \tag {3.1}