봉 요소 강성행렬을 조립해 트러스 변위를 푸는 FEM 직접강성법

다리는 왜 그렇게 처지는가 — 트러스 유한요소의 직접강성법#

1956년, 보잉의 엔지니어 M. J. Turner는 후퇴익의 응력을 손으로 계산하다 벽에 부딪혔다. 부재가 수백 개인 구조를 평형식으로 하나씩 풀 수는 없었다. 그와 동료들이 내놓은 답이 직접강성법(direct stiffness method)이다. 요소 하나의 강성행렬을 만들고, 겹치는 자유도에서 더해 하나의 큰 행렬로 조립한다. 그리고 $Ku = f$ 한 번을 푼다.

이 글은 2차원 트러스를 예로, 봉 요소 하나의 강성행렬을 유도하고 → 국소 좌표를 전역으로 회전시키고 → 전역 행렬로 조립하고 → 경계조건을 넣어 변위를 푸는 전 과정을 Python으로 직접 짠다. 유한요소법의 골격이 실은 선형대수 한 줄이라는 걸 보고 나면, FSI 솔버의 구조 파트나 상용 FEA의 조립 코드도 같은 눈으로 읽힌다.

스프링 하나에서 시작하는 강성#

봉(bar) 요소는 축 방향으로만 힘을 받는 스프링이다. 길이 $L$ , 단면적 $A$ , 탄성계수 $E$ 인 봉의 축 강성은 $k = EA/L$ 이다(늘어난 길이당 필요한 힘). 양 끝 노드의 축 변위를 $u_1, u_2$ 라 하면 노드 힘은 다음과 같다.

\begin{pmatrix} f_1 \\ f_2 \end{pmatrix} = \frac{EA}{L} \begin{pmatrix} 1 & -1 \\ -1 & 1 \end{pmatrix} \begin{pmatrix} u_1 \\ u_2 \end{pmatrix}

$f_1, f_2$ 는 양 끝 노드 힘, 대괄호 안이 국소 강성행렬 $k_e$ 다. 행의 합이 0인 건 강체 이동(두 노드가 같이 움직임)에는 힘이 들지 않는다는 뜻이다. 이 특이성이 나중에 왜 경계조건을 반드시 넣어야 하는지를 설명한다.

국소에서 전역으로 — 회전 변환#

트러스 부재는 사방으로 기울어져 있다. 국소 축 변위를 전역 $x, y$ 성분으로 바꿔야 조립이 된다. 부재가 $x$ 축과 이루는 각을 $\theta$ , $c = \cos\theta$ , $s = \sin\theta$ 라 하자. 축 변위는 전역 변위의 사영이다: $u_{\text{axial}} = c\,u_x + s\,u_y$ . 이 변환을 강성행렬 양쪽에 곱하면 4×4 전역 요소 강성이 나온다.

k_e = \frac{EA}{L} \begin{pmatrix} c^2 & cs & -c^2 & -cs \\ cs & s^2 & -cs & -s^2 \\ -c^2 & -cs & c^2 & cs \\ -cs & -s^2 & cs & s^2 \end{pmatrix}

각 항은 노드1의 $x,y$ 와 노드2의 $x,y$ 자유도를 잇는 강성이다. 아래 시뮬레이션에서 각도를 직접 돌려보자.

Orientation θ: 30°Auto-rotate

Teal cells are positive, pink negative. At θ = 0° only the horizontal DOFs carry stiffness; rotate the bar and the entries redistribute as c², s², and cs — every 2×2 block is rank one.

$\theta = 0$ 이면 $s = 0$ 이라 수직 자유도의 행·열이 통째로 0이 된다. 봉은 자기 축 방향으로만 뻣뻣하다. 각 2×2 블록이 랭크 1인 이유가 여기 있다.

겹치는 자유도를 더한다 — 전역 조립#

조립(assembly)은 이름만 거창하지, 겹치는 자유도에서 강성을 더하는 일이다. 노드 $i$ 의 자유도는 전역 인덱스 $2i, 2i+1$ ( $x, y$ )에 놓는다. 요소의 4개 국소 자유도를 이 전역 인덱스로 매핑(scatter)한 뒤 해당 위치에 요소 강성을 누적한다.

K = \sum_{e} \mathbf{L}_e^\top \, k_e \, \mathbf{L}_e

$\mathbf{L}_e$ 는 요소 자유도를 전역 자유도로 보내는 선택행렬, $K$ 는 전역 강성행렬이다. 실제 코드에선 선택행렬을 곱하지 않고 인덱스 배열로 바로 더한다. 한 노드에 여러 부재가 만나면 그 대각 블록에 여러 요소의 기여가 쌓인다.

경계조건: 고정단을 지운다#

조립만 끝낸 $K$ 는 특이행렬(singular)이다. 구조가 아직 공간에 떠 있어 강체 이동이 자유롭기 때문이다. 지지점(support)에서 변위를 0으로 고정해야 풀린다. 가장 깔끔한 방법은 자유도를 자유(free) 집합 $f$ 와 고정(constrained) 집합 $c$ 로 나누는 것이다.

\begin{pmatrix} K_{ff} & K_{fc} \\ K_{cf} & K_{cc} \end{pmatrix} \begin{pmatrix} u_f \\ u_c \end{pmatrix} = \begin{pmatrix} F_f \\ R_c \end{pmatrix}

$u_c = 0$ (고정)이므로 위쪽 행만 떼면 $K_{ff}\,u_f = F_f$ . 이 축소된 시스템은 비특이라 풀린다. 반력 $R_c$ 는 사후에 $R_c = K_{cf}\,u_f$ 로 회수한다. 아래 트러스에서 하중과 강성을 바꿔보자.

Load P: 24 kNAngle: 0°EA: 8 MNPulse load

Red members are in tension, blue in compression; thickness scales with axial force. Raise EA and the same load bends the truss far less — stiffness is literally the matrix that maps load to displacement.

EA를 키우면 같은 하중에도 처짐이 급감한다. 강성행렬이 곧 하중→변위 사상 그 자체이기 때문이다.

Python — 12부재 트러스를 푼다#

벽에 물린 3패널 캔틸레버 트러스(노드 8개, 부재 12개)를 numpy로 조립하고 푼다. 입력은 노드 좌표·부재 연결·하중, 출력은 노드 변위와 부재 축력이다.

import numpy as np
 
nodes = np.array([[0,0],[0,1],[1,0],[1,1],[2,0],[2,1],[3,0],[3,1]], float)
members = [(0,2),(2,4),(4,6),(1,3),(3,5),(5,7),
           (2,3),(4,5),(6,7),(1,2),(3,4),(5,6)]
EA = 8.0e6           # 축강성 EA [N]
fixed = [0, 1]       # 벽에 고정된 노드
 
def bar_stiffness(p1, p2, EA):
    d = p2 - p1
    L = np.hypot(*d)
    c, s = d / L
    k = EA / L * np.array([[ c*c,  c*s, -c*c, -c*s],
                           [ c*s,  s*s, -c*s, -s*s],
                           [-c*c, -c*s,  c*c,  c*s],
                           [-c*s, -s*s,  c*s,  s*s]])
    return k, L, (c, s)
 
ndof = nodes.shape[0] * 2
K = np.zeros((ndof, ndof))
geom = []
for a, b in members:                       # 전역 강성 조립
    k, L, cs = bar_stiffness(nodes[a], nodes[b], EA)
    dof = [2*a, 2*a+1, 2*b, 2*b+1]
    K[np.ix_(dof, dof)] += k
    geom.append((L, cs))
 
F = np.zeros(ndof)                          # 자유단(노드 6,7)에 24 kN 하중
for n in (6, 7):
    F[2*n+1] -= 12.0e3
 
fixed_dof = [d for n in fixed for d in (2*n, 2*n+1)]
free_dof = [d for d in range(ndof) if d not in fixed_dof]
 
u = np.zeros(ndof)                          # 축소 시스템 K_ff u_f = F_f
u[free_dof] = np.linalg.solve(K[np.ix_(free_dof, free_dof)],
                              F[free_dof])
 
for m, (a, b) in enumerate(members):        # 부재 축력 N = (EA/L)[-c,-s,c,s]·u
    L, (c, s) = geom[m]
    dof = [2*a, 2*a+1, 2*b, 2*b+1]
    N = EA / L * np.array([-c, -s, c, s]) @ u[dof]
    print(f"member {a}-{b}: N = {N/1e3:+7.2f} kN "
          f"({'tension' if N > 0 else 'compression'})")
 
print(f"free-end drop = {u[2*6+1]*1e3:.3f} mm")

np.ix_로 인덱스 격자를 만들어 요소 강성을 정확한 전역 위치에 더한다. 이 20줄이 상용 FEA 솔버 심장부와 구조적으로 동일하다.

강성행렬이 특이해질 때#

현장에서 가장 흔한 실패는 "행렬이 특이하다"는 에러다. 원인은 대개 셋 중 하나다.

첫째, 경계조건 누락. 강체 이동을 막는 지지가 부족하면 $K_{ff}$ 가 여전히 특이하다. 2D는 최소 3개, 3D는 6개 자유도를 구속해야 한다.

둘째, 메커니즘(mechanism). 삼각형화가 안 된 사각형 패널은 부재가 모자라 흐물거린다. 트러스는 항상 삼각형으로 채워야 한다.

셋째, 0 길이·중복 노드. 두 노드가 같은 좌표면 $L = 0$ 에서 나눗셈이 폭발한다. 메시 병합 전에 좌표 중복을 걸러라.

다시 읽지 않을 사람을 위한 요약#

유한요소법의 뼈대는 "요소 강성 만들기 → 전역으로 회전 → 겹치는 자유도 더하기 → 경계조건 넣고 $Ku = f$ 풀기"다.
조립 전 $K$ 는 항상 특이하다. 강체 이동을 막는 지지가 비특이성을 만든다.
봉 요소의 2×2 블록은 랭크 1 — 축 방향으로만 뻣뻣하다. 각도 변환이 그것을 전역 좌표로 흩뿌린다.