각도 0에서 수직 칸은 0이고, 벽의 4자유도를 남기면 16×16은 특이행렬로 남는다.

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

1956년, 보잉의 M. J. Turner는 후퇴익 응력을 손으로 밀다 막혔다. 부재가 수백 개면 평형식을 하나씩 못 푼다. 동료들과 내놓은 답이 직접강성법(direct stiffness method)이다. 요소 강성을 만들고, 자유도가 겹치는 자리에 더하고, $Ku = f$ 를 한 번 푼다.

각도부터 보자. 봉이 수평으로 누우면 수직 칸이 죽는다. 각도 슬라이더를 0으로 내리면, 4×4 강성에서 수직 자유도의 행과 열이 통째로 0이 된다.

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일 때 비는 칸을 보자. 봉은 자기 축으로만 뻣뻣하다. 각도를 올리면 그 0이 $x$ 와 $y$ 로 나뉜다.

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

봉(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$ 를 잇는다. $\theta = 0$ 이면 $s = 0$ 이고 $c = 1$ 이라, 위 행렬은 $EA/L$ 을 곱한 채 첫 행이 1, 0, −1, 0이 되고 그다음 행은 전부 0이다. 세 번째 행은 첫 행의 부호를 바꾼 것이고, 마지막 행도 0이다. 2×2 블록이 랭크 1인 이유가 이 0들이다. 좌표가 휘면 강성 텐서도 같이 돌아간다. 그 변환은 메트릭으로 구성 텐서를 돌리는 경우에 따로 있다.

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

조립은 겹치는 자유도에서 강성을 더하는 일이다. 노드 $i$ 의 자유도는 전역 인덱스 $2i, 2i+1$ 에 둔다. $x$ 다음이 $y$ 다. 요소의 국소 자유도 4개를 이 인덱스로 보낸 뒤 그 자리에 요소 강성을 누적한다.

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

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

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

조립만 한 $K$ 는 특이행렬이다. 구조가 아직 떠 있어서 강체 이동이 남는다. 지지점에서 변위를 0으로 묶어야 풀린다. 자유도를 자유 집합 $f$ 와 고정 집합 $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 = K_{cf}\,u_f$ 로 되찾는다. $R_c$ 는 지지점이 구조에 되돌려 주는 힘이다.

하중 슬라이더와 EA를 같이 올리면, 끝단 처짐이 어떻게 줄어드는지 바로 보인다.

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를 키우면 같은 하중에도 처짐이 줄어든다. 강성행렬이 하중에서 변위로 가는 그 사상이다. 하중을 고정한 채 EA만 두 배로 올려, 끝단 처짐이 반으로 가는지 보면 된다.

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

벽에 물린 3패널 캔틸레버다. 노드 8개, 부재 12개. 좌표와 연결, 하중을 넣고 변위와 축력을 낸다.

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_가 요소 강성을 전역 칸에 더한다. 이 루프가 상용 FEA 조립과 같은 구조다.

노드 간격은 1 m다. 수평·수직 부재 9개의 길이는 $L = 1$ 이고, 코드의 EA가 $8.0\times 10^6$ N이면 $k = EA/L = 8.0\times 10^6$ N/m다. 끝단 노드 하나의 수직 하중은 $12\times 10^3$ N이므로, 이 봉 하나만 놓고 밀면 변위는 $1.5\times 10^{-3}$ m, 곧 1.5 mm다. 트러스 끝단은 이 나눗셈이 아니다. 부재가 하중을 나눠 갖고, 실제로 푸는 행렬은 축소된 $K_{ff}$ 다.

대각 부재는 세 개다. (0, 1)과 (1, 0)을 잇는 부재의 길이는 $\sqrt{2}$ 라 강성은 수평 부재의 $1/\sqrt{2}$ 배다. 같은 EA라도 대각선은 더 소프트하다. 그 부재의 각은 45°다. $c^2$ 와 $s^2$ 는 $1/2$ 이고, $cs$ 의 절대값도 $1/2$ 다. 4×4에서 0이 아닌 칸은 전부 $EA/L$ 의 절반, 길이가 $\sqrt{2}$ 이니 칸의 크기는 $EA/(2\sqrt{2})$ 다. 수평 부재는 수직 칸이 0으로 비고, 대각 부재는 네 방향이 그 크기로 채워진다. 같은 EA가 각도만으로 이렇게 갈린다. 노드 8개의 자유도는 16개다. 벽의 노드 0과 1이 4개를 묶으므로 푸는 시스템은 12×12다. 2차원에서 최소 구속은 3개인데 이 예는 4개를 묶는다. 벽 노드 하나의 $y$ 만 놓쳐도 개수는 최소에 걸치고, 둘 다 놓치면 $K_{ff}$ 가 다시 특이가 된다.

강성행렬이 특이해질 때#

현장에서 제일 많이 보는 실패는 특이행렬이다. 갈래는 세 가지다.

경계조건이 빠지면 강체 이동이 남고 $K_{ff}$ 가 특이하다. 2차원은 자유도 3개, 3차원은 6개를 적어도 묶어야 한다.

삼각형이 아닌 사각형 패널은 부재가 모자란 메커니즘이다. 트러스는 삼각형으로 채워야 흔들리지 않는다.

길이가 0이거나 노드가 겹치면 $L = 0$ 에서 나눗셈이 터진다. 좌표가 같은 점을 합치기 전에 걸러야 한다.

벽 노드 하나를 놓치면 조립은 거기서 멈춘다#

수평 봉의 1.5 mm는 부재 하나의 이야기다. 12부재 캔틸레버의 끝단 처짐은 $K_{ff}\,u_f = F_f$ 를 풀어야 나오고, 벽의 4자유도가 빠지지 않으면 그 식이 성립하지 않는다. 지지가 모자라 강체 모드가 남으면 변위는 숫자가 아니라 영공간이다.

봉을 면 요소로 바꾸는 순간 같은 조립이 다른 이유로 굳는다. 얇은 쉘에서 그 잠김이 어떻게 나타나는지는 전단 잠김과 MITC가 강성을 되돌리는 과정에 있다. 트러스에서 지지 4개를 먼저 확인하고, 두께를 줄인 모델에서만 그쪽으로 넘어가면 된다.

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