OpenSees 좌굴과 한계점: 요소 하나짜리 기둥은 좌굴하지 못하고, 선형 변환도 마찬가지다

요소 하나로 모델링한 기둥은 좌굴할 수 없습니다. "부정확하게 좌굴한다"가 아닙니다. 어떤 하중에서도 접선강성이 특이해지지 않습니다. 두 지점 사이의 elasticBeamColumn 하나에는 내부 절점이 없고, 반정현파 좌굴 모드가 존재할 자리가 없기 때문입니다. geomTransf Linear로 모델링한 기둥도 좌굴할 수 없습니다. 어떤 메시에서도 그렇고, 같은 이유가 다른 옷을 입은 것입니다. 기하강성이 없으면 축력이 횡방향 강성을 줄일 방법이 없습니다.

두 모델 모두 깔끔하게 돌아갑니다. 둘 다 좌굴하중의 열 배에서 수렴한 결과를 보고합니다. 어느 쪽도 한마디도 하지 않습니다. 이 글은 좌굴 명령이 없는 프레임워크에서 좌굴하중을 찾는 법, 초기 불완전이 질문 자체를 어떻게 바꾸는지, 그리고 하중제어뿐 아니라 변위제어까지 무너뜨리는 유일한 경로 형태에 관한 것입니다. 모든 것을 OpenSeesPy 3.8.0에서 닫힌해와 대조해 측정했습니다. 이 글은 13부입니다. 1–12부에서 객체 모델, 탄성 골조, 파이버 단면, 동적해석, 소스 읽기, 철근콘크리트, 지진기록, 3차원 건물, 수렴 실패, Tcl, 연속체 요소, 후처리를 다뤘습니다.

1. 좌굴 명령 없이 좌굴하중 찾기

OpenSees에는 선형 좌굴 해석이 없습니다. 필요하지도 않습니다. 좌굴하중의 정의를 그대로 쓸 수 있기 때문입니다. Pcr는 접선강성이 양정치이기를 그만두는 하중입니다. 네 줄이면 도달합니다.

def min_eig_of_tangent(P): build_model() # 공회전 좌표계 변환 ops.system('FullGeneral') # printA 가 조밀 K 를 주도록 ops.integrator('LoadControl', P) ops.analysis('Static'); ops.analyze(1) n = ops.systemSize() K = np.array(ops.printA('-ret')).reshape((n, n)) return np.min(np.linalg.eigvalsh(0.5 * (K + K.T))) # 그다음 min_eig_of_tangent(P) 가 0을 지날 때까지 P 를 이분한다

세부 둘이 중요합니다. system('FullGeneral')이 필요한 것은 printA가 솔버가 들고 있는 것을 그대로 돌려주기 때문입니다. 밴드나 희소 시스템이면 직접 풀어야 하는 압축 배열이 나옵니다. 그리고 eigvalsh 전에 0.5(K + Kᵀ)로 대칭화하면 공회전 좌표계 요소가 남기는 반올림 비대칭이 제거됩니다.

길이 3 m, EI = 2000 kN·m²인 기둥에 요소 20개로 돌린 결과입니다.

단부조건kPcr = π²EI/(kL)²OpenSees오차

고정–자유2548.3114 kN548.7438+0.08%

핀–핀12193.2454 kN2200.1792+0.32%

고정–핀0.69924486.2595 kN4515.9240+0.66%

고정–고정0.58772.9817 kN8884.8417+1.28%

네 오차가 모두 양수인 것이 변위기반 이산화의 특징입니다. 가정된 3차 형상함수가 좌굴 형상을 구속하므로 모델이 연속체보다 뻣뻣하고 좌굴하중이 높게 나옵니다. 단부가 뻣뻣할수록 오차가 커지는데, 유효길이가 짧을수록 요소당 표현해야 할 곡률이 크기 때문입니다.

2. 기둥이 좌굴이라도 하려면 요소가 몇 개 필요한가

요소 수Pcr오차

1어떤 하중에서도 교차 없음—

22670.23 kN+21.75%

42310.81 kN+5.36%

82224.05 kN+1.40%

162202.73 kN+0.43%

322197.42 kN+0.19%

오차는 대략 1/n²로 떨어지므로 메시를 두 배로 하면 오차가 4분의 1이 됩니다. 소성힌지 모델링에서 넘어온 습관인 부재당 요소 4개는 좌굴하중에 5% 오차를 남깁니다. 8개면 1.4%가 됩니다.

변환은 같은 이야기의 나머지 절반입니다.

geomTransf요소 20개에서의 PcrPcr / 정확값

Linear결코 특이해지지 않음—

PDelta2197.76 kN1.0021

Corotational2200.18 kN1.0032

PDelta와 Corotational 둘 다 분기점을 0.3% 이내로 찾고, 이 문제에서는 PDelta가 근소하게 더 가깝게 나왔습니다. 이산화 오차 둘 사이의 천분의 일 차이일 뿐, 그쪽을 선호할 이유는 아닙니다. 둘은 나중에, 즉 PDelta의 미소회전 가정이 깨지는 좌굴 후 경로에서 갈라집니다. 8부의 leaning column이 PDelta를 쓴 이유이고, 이 글이 Corotational을 쓰는 이유입니다.

함정은 Linear입니다. 대부분의 튜토리얼의 기본 선택이고, 모든 정적해석을 아름답게 수렴시키며, 현상 자체를 조용히 지웁니다. 모델에 세장한 압축부재가 있는데 geomTransf Linear가 있다면, 그 모델의 안정성은 아예 검토되고 있지 않습니다.

3. 불완전한 기둥은 좌굴하지 않는다. 증폭한다

위의 네 수치는 분기 하중입니다. 완전히 곧은 기둥이 두 번째 평형경로를 얻는 하중이죠. 실제 기둥은 곧지 않고, 불완전한 기둥의 비선형 해석은 결코 실패하지 않고, 정점도 없으며, Pcr를 알려주지도 않습니다. 쌍곡선을 따라갈 뿐입니다.

진폭 e₀인 정현파 초기 휨에 대한 고전적 결과는 이렇습니다.

delta(P) = e0 / (1 - P/Pcr)

모델 자신의 좌굴하중 2200.18 kN을 기준으로 측정한 값입니다.

e₀P/Pcrδ OpenSeese₀/(1−P/Pcr)오차

L/10000.49340.005922 m0.0059220.000%

L/10000.80740.015578 m0.0155800.015%

L/5000.80740.031142 m0.0311600.059%

L/2000.80740.077617 m0.0779010.364%

좌굴하중의 절반에서 정확하고, 변형이 2차항이 문제 될 만큼 커질 때만 벌어집니다. 그래서 오차가 e₀에 따라 커지고 다른 무엇에도 따르지 않습니다. 놓치기 쉬운 공학적 요점은 이것입니다. 불완전한 탄성 기둥의 비선형 해석에는 찾아낼 파괴 사건이 없습니다. 0.99 Pcr까지 밀어도 수렴합니다. 변위가 클 뿐입니다. 수렴하지 않는 지점을 지켜보며 좌굴하중을 찾고 있다면 찾지 못할 것이고, 만약 찾았다면 그것은 안정 한계가 아니라 재료 파괴이거나 수치적 실패입니다.

4. Southwell:도달한 적 없는 좌굴하중을 해석에서 뽑기

증폭식을 정리하면 δ/P = δ/Pcr + e₀/Pcr이므로, δ/P를 δ에 대해 그리면 기울기가 1/Pcr인 직선입니다. 좌굴하중 근처의 측정값이 전혀 필요 없습니다.

e₀기울기에서 얻은 Pcr모델의 Pcr오차R²

L/10002200.25 kN2200.18 kN+0.003%1.000000

L/5002200.47 kN+0.013%1.000000

L/2002202.02 kN+0.084%1.000000

추정하려는 하중의 80%에서 멈춘 데이터로 0.003% 오차입니다. 접선 고유값 방법이 현실적이지 않을 때 — 큰 모델, 비탄성 재료, 임계 모드를 분리할 수 없는 구조 — 꺼내 쓸 기법입니다. 실물 실험에서 같은 양을 뽑아내는 방법이기도 해서, 모델과 실험을 비교하는 자연스러운 접점이 됩니다.

5. 스냅백: 변위제어도 따라갈 수 없는 경로

9부에서 한계점에서 하중제어는 그냥 지나쳐 가고 변위제어는 하강 분기를 정확히 따라간다는 것을 보였습니다. 제어하는 변위가 경로를 따라 단조 증가할 때 참인 이야기입니다. 항상 그렇지는 않습니다.

9부의 얕은 von Mises 트러스를 가져와 정점과 재하점 사이에 800 kN/m 탄성 스프링을 직렬로 넣습니다. 이제 재하점 변위는 u₄ = u₃ + P/k이고, 하강 분기에서는 하중이 충분히 빨리 떨어져 P/k가 u₃가 커지는 속도보다 빨리 줄어듭니다. 제어하는 변위가 거꾸로 갑니다.

정확 경로는 u₄ = 0.143925 m에서 0.056075 m까지 87.85 mm 되돌아가고, 그동안 트러스는 계속 변형합니다. 그 한 문제에 적분기 넷을 걸었습니다.

적분기결과근거

LoadControl실패 없음, 전부 건너뜀u₃가 0.227 m로 점프

DisplacementControl회귀점에서 실패회귀점 0.143925에 대해 u₄ = 0.1436 m에서 정지

ArcLength(0.0015)경로 전체를 추적u₄에서 역방향 87스텝, u₃가 0.292 m 도달

ArcLength(0.0010)왔던 경로로 되돌아감역방향 254스텝, u₃가 0.043 m에서 정체

가져갈 것이 셋입니다.

변위제어는 정직하고 정확하게 실패했습니다. 참 회귀점 0.143925 m에 대해 0.1436 m에서 멈췄습니다. 실패한 위치가 곧 진단입니다. 특정 변위에서 멈추고 더 나아가기를 거부하는데 잔차는 여전히 줄고 있다면, 그 해석은 경로가 거기서 꺾인다고 말하고 있는 것입니다.

호장법은 자동이 아닙니다. 호장을 0.0015에서 0.0010으로 줄였더니 좋아진 게 아니라 나빠졌습니다. 회귀점에서 부호가 뒤집혀 방금 올라온 분기를 되짚어 내려갔고, 254스텝을 거꾸로 쓰면서 u₃ = 0.043 m를 끝내 넘지 못했습니다. 호장법은 직전 스텝에서 방향을 고르는데, 급한 회귀점 근처에서는 그 선택이 틀릴 수 있습니다. 증상은 상태변수가 왔던 곳으로 돌아가는데도 해석이 "성공"하는 것이고, 그래서 여기서의 진단 도구는 analyze()의 반환값이 아니라 전진해야 할 양의 증분 부호입니다.

그리고 또 한 번, 아무 불평도 하지 않은 유일한 것은 하중제어였습니다.

6. 분기점, 한계점, 스냅백

무슨 일이 일어나는가무엇으로 찾는가

분기점두 번째 경로가 생긴다. 접선이 특이하지만 주경로는 계속된다접선 최소 고유값이 0을 지나는 지점

한계점하중이 최대를 가진다. 경로는 하중이 줄면서 계속된다변위제어 또는 호장법

스냅백제어하는 변위도 최대를 가진다호장법만, 그것도 호장을 맞춰야

불완전 기둥위의 어느 것도 아니다. 매끄러운 쌍곡선Southwell, 또는 완전 모델의 접선

마지막 줄의 실무적 귀결은 따로 말할 값어치가 있습니다. "좌굴시키려고" 모델에 불완전을 넣는 것은 정반대의 일을 합니다. 분기점을 없애고, 그 자리에 특별한 점이 하나도 없는 곡선을 놓습니다. 불완전은 사용하중에서의 응답이 필요할 때, 또는 구조물에 진짜 한계점이 있어 그 정점이 필요할 때 모델에 넣는 것입니다. 좌굴하중을 찾는 데는 도움이 되지 않습니다.

7. 안정 문제 점검표

변환이 PDelta나 Corotational인가? Linear는 좌굴하지 못한다.

압축부재당 요소가 최소 4개, 가급적 8개인가?

찾고 있는 것이 분기점인가, 한계점인가, 스냅백인가?

분기점이라면 수렴 거동이 아니라 접선의 최소 고유값을 확인했는가?

모델에 불완전이 있다면, 그것을 넣음으로써 분기점을 없앴다는 것을 알고 있는가?

제어하는 자유도가 예상하는 경로를 따라 단조 증가하는가?

호장법을 쓴다면 두 번째 호장으로도 돌려봤고, 둘이 일치하는가?

analyze()가 성공을 보고하는데 거꾸로 가고 있는 상태변수가 있는가?

가장 가까운 고전 사례의 π²EI/(kL)²와 대략적으로라도 대조했는가?

8. 같은 문제를 다섯 단계 깊이로

입문자는 세장한 기둥에 좌굴하중이 있다는 것을 압니다. 주니어는 공회전 좌표계 변환을 쓰고 부재당 요소를 하나보다 많이 씁니다. 실무 해석자는 해석이 수렴을 멈추는 지점이 아니라 접선강성에서 Pcr를 찾습니다. 시니어는 분기점과 한계점을 구분해 적분기를 거기 맞추고, 불완전이 분기점을 지운다는 것을 압니다. 전문가는 반환값이 아니라 상태변수를 지켜봅니다. 회귀점에서는 성공한 스텝과 역주행한 스텝이 밖에서 보기에 똑같기 때문입니다.

다음 글은 Windows에서 OpenSees를 돌리는 이야기입니다. 설치, 거기서만 나오는 오류, 경로, 그리고 툴체인 중 다르게 동작하는 부분들.

전체 스크립트

좌굴하중을 구하는 접선 고유값 방법을 고전 단부조건 넷에 적용하고, 이어서 아예 좌굴할 수 없는 두 가지 모델 선택에 적용합니다.

# Part 13 - a buckling load with no buckling command # units: kN, m import math import numpy as np import openseespy.opensees as ops L, E, I, A = 3.0, 2.0e8, 1.0e-5, 0.01 EI = E * I CASES = {'pinned-pinned': 1.0, 'fixed-free': 2.0, 'fixed-fixed': 0.5, 'fixed-pinned': 0.6992} def build(case, nele, transf='Corotational'): ops.wipe(); ops.model('basic', '-ndm', 2, '-ndf', 3) for i in range(nele + 1): ops.node(i + 1, 0.0, i * L / nele) if case == 'pinned-pinned': ops.fix(1, 1, 1, 0); ops.fix(nele + 1, 1, 0, 0) elif case == 'fixed-free': ops.fix(1, 1, 1, 1) elif case == 'fixed-fixed': ops.fix(1, 1, 1, 1); ops.fix(nele + 1, 1, 0, 1) else: ops.fix(1, 1, 1, 1); ops.fix(nele + 1, 1, 0, 0) ops.geomTransf(transf, 1) for i in range(nele): ops.element('elasticBeamColumn', i + 1, i + 1, i + 2, A, E, I, 1) ops.timeSeries('Linear', 1); ops.pattern('Plain', 1, 1) ops.load(nele + 1, 0.0, -1.0, 0.0) def min_eig(case, nele, P, transf='Corotational'): """The smallest eigenvalue of the tangent stiffness under an axial load P.""" build(case, nele, transf) ops.constraints('Transformation'); ops.numberer('Plain') ops.system('FullGeneral') # so printA returns a dense matrix ops.test('NormDispIncr', 1e-10, 30); ops.algorithm('Newton') ops.integrator('LoadControl', P); ops.analysis('Static'); ops.analyze(1) n = ops.systemSize() K = np.array(ops.printA('-ret')).reshape((n, n)) return float(np.min(np.linalg.eigvalsh(0.5 * (K + K.T)))) def critical(case, nele=20, transf='Corotational'): lo, hi = 1e-3, 3.0 * math.pi**2 * EI / (CASES[case] * L)**2 f_lo = min_eig(case, nele, lo, transf) if f_lo * min_eig(case, nele, hi, transf) > 0: return None # never goes singular for _ in range(50): mid = 0.5 * (lo + hi) if f_lo * min_eig(case, nele, mid, transf) <= 0: hi = mid else: lo = mid return 0.5 * (lo + hi) print(f'EI = {EI:.1f} kN m2, L = {L} m, 20 elements\n') print(f"{'end conditions':16s}{'k':>8s}{'exact (kN)':>13s}{'OpenSees':>13s}{'error':>9s}") worst = 0.0 for case, k in CASES.items(): ex = math.pi**2 * EI / (k * L)**2 got = critical(case) err = abs(got - ex) / ex; worst = max(worst, err) print(f'{case:16s}{k:8.4f}{ex:13.4f}{got:13.4f}{100*err:8.2f}%') assert worst < 0.02, f'worst error {100*worst:.2f}% is larger than expected' print('\nwith geomTransf Linear, and with a single element:') for label, kw in (('Linear transformation', dict(transf='Linear')), ('one element', dict(nele=1))): got = critical('pinned-pinned', **kw) print(f' {label:24s} -> ' + ('no crossing at any load' if got is None else f'{got:.4f} kN')) assert got is None, f'{label} was not expected to buckle' print('\nOK four classical loads recovered; neither shortcut can buckle at all')

실행하면 이렇게 나옵니다.

EI = 2000.0 kN m2, L = 3.0 m, 20 elements end conditions k exact (kN) OpenSees error pinned-pinned 1.0000 2193.2454 2200.1792 0.32% fixed-free 2.0000 548.3114 548.7438 0.08% fixed-fixed 0.5000 8772.9817 8884.8417 1.28% fixed-pinned 0.6992 4486.2595 4515.9240 0.66% with geomTransf Linear, and with a single element: Linear transformation -> no crossing at any load one element -> no crossing at any load OK four classical loads recovered; neither shortcut can buckle at all

OpenSees 시리즈 전체

각 편은 닫힌해·독립적인 풀이·OpenSees 소스 중 하나와 대조해 검증했다. 순서대로 읽도록 썼지만 각 편이 자기 전제를 스스로 밝힌다.

1부 — 도메인 모델과 해석 조립

2부 — 단면·기하변환·분포하중

3부 — 파이버 단면·모멘트–곡률·푸시오버

4부 — 질량·감쇠·시간적분

5부 — 라이브러리 읽기와 재현 가능한 워크플로

6부 — 철근콘크리트 파이버 단면과 구속

7부 — 실제 지진기록·층간변위·증분동적해석

8부 — 3차원: 강막·비틀림·leaning column

9부 — 수렴하지 않을 때: 실패·알고리즘·한계점

10부 — Tcl 읽기와 공개 스크립트 이식

11부 — 쉘·솔리드·지반: 락킹·메시·부지응답

12부 — 후처리: recorder·질의, 그리고 그림이 증명하는 것

13부 — 좌굴과 한계점: 분기점·스냅백·호장법 — 지금 보는 글

14부 — Windows: 설치·파일 형식·조용한 실패

15부 — 검증된 정답이 딸린 연습문제와 한영 용어집

16부 — 어떻게 동작하는가: 시행상태·확정상태·해석 루프

17부 — 커뮤니티가 말하는 것, 측정해보면

18부 — 절점 변위에서 파이버 응력까지

19부 — 면진과 감쇠: 속도 의존 요소와 베어링

참고 자료

Timoshenko, S. P., Gere, J. M., Theory of Elastic Stability, McGraw-Hill — 고전 네 사례와 증폭식.

Southwell, R. V., "On the analysis of experimental observations in problems of elastic stability", Proceedings of the Royal Society A 135, 1932.

Crisfield, M. A., Non-linear Finite Element Analysis of Solids and Structures — 호장법, 스냅백, 회귀점에서의 방향 제어.

OpenSeesPy, printA, systemSize — 1절의 접선강성 되읽기.

OpenSeesPy, Corotational 변환, PDelta 변환.

OpenSeesPy, ArcLength, DisplacementControl.

Scott, M. H., OpenSees Digital — 기하변환과 안정성에 관한 반복적인 실무 해설.

출처 확인 2026-08-31. 모든 결과는 OpenSeesPy 3.8.0에서 산출했으며 본문의 스크립트로 재현할 수 있습니다. 좌굴하중은 탄성 기둥 하나와 트러스 하나에 대한 것이며, 옮겨 쓸 수 있는 결론은 방법, 각 적분기의 실패 양상, 그리고 선형 변환이나 요소 하나짜리 부재가 모델에서 안정성을 통째로 제거한다는 사실입니다. 결과는 프레임워크의 거동을 보여주는 것으로, 프로젝트별 해석이나 설계기준 검토를 대체하지 않습니다.