강성행렬

1. 개요

강성행렬 Stiffness Matrix
분야	유한요소법 × 구조해석
지배 방정식	K u = f
주요 성질	대칭성, 희소성, 양정치성
기호	K (전역), ke (요소)

구조물을 풀면 결국 딱 한 줄이 남는다. Ku = f. 나머지는 이 행렬을 어떻게 만드느냐의 문제다.

강성행렬(Stiffness Matrix, $\mathbf{K}$)은 유한요소법에서 구조물의 변위와 그에 필요한 힘 사이의 관계를 담은 행렬로, 선형 정적 구조해석의 지배 방정식 $\mathbf{K}\mathbf{u} = \mathbf{f}$의 핵심이다. 여기서 $\mathbf{u}$는 절점 변위 벡터, $\mathbf{f}$는 절점 하중 벡터다. 스칼라 스프링의 후크 법칙 $f = ku$를 다자유도로 일반화한 것이라 보면 정확하다. 강성행렬은 “이 절점을 이만큼 밀면, 저 절점에 얼마의 힘이 걸리는가”라는 구조물의 저항 특성을 통째로 인코딩한 물건이다.

유한요소 해석의 실체는 사실상 이 거대한 행렬을 조립하고, $\mathbf{K}\mathbf{u}=\mathbf{f}$라는 선형 연립방정식을 푸는 일이다. 자유도가 수백만~수억 개에 이르는 대규모 문제에서, 이 행렬의 성질을 잘 이해하고 활용하느냐가 해석 시간을 하루로 만들지 한 달로 만들지를 가른다.¹

2. 요소 강성행렬과 조립

전역 강성행렬은 하늘에서 뚝 떨어지지 않는다. 먼저 구조물을 잘게 나눈 각 요소마다 요소 강성행렬(element stiffness matrix) $\mathbf{k}^e$를 구한다. 이는 변형률-변위 행렬 $\mathbf{B}$(형상함수 $N_i$를 미분해 만든 것으로, $\boldsymbol{\varepsilon} = \mathbf{B}\mathbf{u}^e$)와 재료 물성 행렬 $\mathbf{D}$(구성방정식 $\boldsymbol{\sigma} = \mathbf{D}\boldsymbol{\varepsilon}$, 응력과 변형률 참고)로부터 다음 적분으로 얻어진다. 변형에너지 $\tfrac12\int \boldsymbol{\varepsilon}^T\mathbf{D}\boldsymbol{\varepsilon}\,d\Omega$에 $\boldsymbol{\varepsilon} = \mathbf{B}\mathbf{u}^e$를 넣으면 곧바로 $\tfrac12 \mathbf{u}^{eT}\mathbf{k}^e\mathbf{u}^e$가 나오므로, 아래 식은 에너지의 2차 형식을 행렬로 적은 것에 불과하다.

$$\mathbf{k}^e = \int_{\Omega^e} \mathbf{B}^T \mathbf{D} \, \mathbf{B} \; d\Omega$$

이 적분은 대개 가우스 구적으로 수치적분하며, 뒤틀린 요소는 등매개변수 요소 사상의 자코비안 행렬식을 곱해 기준요소 위에서 적분한다. 이 식 자체가 어디서 나오는지는 갤러킨 방법이 약형식에서 유도한다. 가장 간단한 예로, 길이 $h$·단면적 $A$·영률 $E$인 1차원 봉 요소는 $\mathbf{k}^e = \frac{EA}{h}\begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}$이다. 그다음, 각 요소가 공유하는 절점을 매개로 전역 자유도 위치에 요소 행렬을 더해 넣는 과정이 조립(assembly)이다.

$$\mathbf{K} = \sum_{e} \mathbf{L}^{e^T} \, \mathbf{k}^e \, \mathbf{L}^e$$

여기서 $\mathbf{L}^e$는 요소의 지역 자유도를 전역 자유도로 매핑하는 위치 행렬이다. 실제 구현에서는 이 행렬 곱을 명시적으로 하지 않고, 인접 요소들이 공유 절점에서 기여를 겹쳐 더하는 “스캐터-애드(scatter-add)” 방식으로 처리한다. 이 겹침이 바로 서로 연결된 절점끼리만 행렬 성분이 생기는 구조를 만들며, 뒤에 나올 희소성의 원천이 된다.

3. 세 가지 핵심 성질

강성행렬이 사랑받는(그리고 활용되는) 이유는 세 가지 아름다운 성질 덕분이다.

3.1. 대칭성 (Symmetry)

$\mathbf{K} = \mathbf{K}^T$가 성립한다. 이는 우연이 아니라 물리 법칙의 결과다. 선형 탄성체는 보존계라 변형에너지가 존재하고, 강성행렬은 이 에너지의 2차 형식의 헤시안이기 때문이다. 맥스웰-베티 상반정리(절점 $i$를 밀 때 $j$에서 느끼는 힘 = 그 반대)가 대칭성의 물리적 표현이다.² 거꾸로 말하면 대칭이 깨지는 곳에는 이유가 있다 — 종동력(follower load)이나 비연관 소성, 마찰 접촉의 접선강성은 비대칭이 되어 대칭 솔버를 못 쓴다(비선형 구조해석). 실용적으로는 저장 공간과 계산량이 절반으로 줄어든다는 뜻이라, 매우 반가운 성질이다.

3.2. 희소성 (Sparsity)

한 절점은 자신과 요소를 공유하는 이웃 절점하고만 직접 연결된다. 따라서 강성행렬의 대부분 성분은 0이고, 0이 아닌 성분은 대각선 주변에 몰려 있다. 100만 자유도 문제라면 행렬은 $10^{12}$개 성분을 가지지만, 실제 0이 아닌 것은 그 중 수천만 개 남짓이다. 3차원 8절점 육면체 격자로 세면, 내부 절점 하나는 이웃 27개 절점(자기 포함)과 연결되고 절점당 3자유도라 한 행에 $27 \times 3 = 81$개, 전체 약 8100만 개다(전체의 0.008%). 이 **희소성**을 무시하고 꽉 찬(dense) 행렬로 다뤘다간 배정밀도 $10^{12} \times 8$바이트 = 8TB가 필요하다. CSR(값 8바이트 + 열 인덱스 4바이트)로 저장하면 약 1GB, 대칭성으로 절반만 저장하면 그 절반이다.

3.3. 양정치성 (Positive Definiteness)

강체 운동(rigid body motion)을 경계 조건으로 제대로 구속한 강성행렬은 대칭 양정치(SPD)다. 즉 0이 아닌 모든 변위에 대해 변형에너지 $\frac{1}{2}\mathbf{u}^T\mathbf{K}\mathbf{u} > 0$이다. 물리적으로 “0이 아닌 변형은 항상 양의 에너지를 요구한다”는 당연한 사실이다. 이 SPD 성질 덕분에 피벗팅 없이 촐레스키 분해가 가능하고, SPD 전용 반복법인 켤레기울기법(CG)을 쓸 수 있다. 또 $\mathbf{K}\mathbf{u}=\mathbf{f}$의 해는 총 퍼텐셜 에너지 $\tfrac12\mathbf{u}^T\mathbf{K}\mathbf{u} - \mathbf{f}^T\mathbf{u}$의 유일한 최솟점이기도 하다 — 구조해석이 곧 2차 함수 최소화라는 뜻이고, CG가 잘 듣는 이유가 여기 있다.

4. 경계 조건 처리

조립 직후의 강성행렬은 사실 특이(singular)하다. 아무 데도 고정하지 않은 구조물은 허공에서 자유롭게 떠다닐 수 있어(강체 운동), $\mathbf{K}$가 특이해 방정식이 풀리지 않는다. 강체 운동은 변형에너지가 0인 변위라 $\mathbf{K}$의 영공간을 이루고, 그 차원은 1차원 봉이 1(병진), 2차원 평면 문제가 3(병진 2 + 회전 1), 3차원 고체가 6(병진 3 + 회전 3)이다. 따라서 3차원에서는 최소 6개의 독립 구속이 필요하다. 구속 없는 구조물의 모드 해석에서 고유진동수 0인 모드가 정확히 6개 나오는 것이 같은 사실이고, 그 개수가 6을 넘으면 요소 연결이 빠진 곳(메커니즘)이 있다는 신호다. 그래서 경계 조건으로 일부 자유도를 구속해줘야 비로소 $\mathbf{K}$가 양정치가 되고 유일해가 존재한다. 지정 변위(디리클레 조건)를 부과하는 방법은 크게 셋이다.

• 행/열 소거: 구속된 자유도에 해당하는 행과 열을 아예 시스템에서 제거한다. 가장 깔끔하지만 인덱스 관리가 번거롭다.
• 페널티법: 해당 대각 성분에 아주 큰 수를 더해 그 자유도를 사실상 고정한다. 구현은 쉽지만 조건수를 악화시킨다.³
• 라그랑주 승수법: 구속을 별도 미지수로 도입한다. 정확하지만 행렬 크기가 커지고 양정치성이 깨진다(대각 블록이 0인 안장점 형태라 부정치가 된다).

5. 검산 — 조건수와 페널티의 값

1차원 봉($E = A = 1$, 길이 1)을 요소 $N$개로 나누고 왼쪽 끝을 고정하면 $h = 1/N$이고 요소 강성은 $1/h = N$이다. 고정 전 행렬의 최소 고유값은 반올림 수준($\sim 10^{-16}$)으로 강체 병진 모드 하나가 보이고, 고정 후의 조건수는 다음과 같다. 해석해는 $\kappa = \sin^2\!\frac{(2N-1)\pi}{4N+2} \big/ \sin^2\!\frac{\pi}{4N+2}$이다.

$N$	고정 후 조건수 (실측 = 해석해)	$\kappa / N^2$
4	29.3	1.83
8	113.5	1.77
16	437.7	1.71
32	1708.7	1.67

$\kappa/N^2$이 $16/\pi^2 \approx 1.62$로 다가간다. 즉 격자를 두 배로 세밀하게 하면 조건수는 네 배($\kappa \propto h^{-2}$)가 되고, 전처리 없는 CG의 반복 수는 $\sqrt\kappa$에 비례해 두 배가 된다. 2·3차원 2차 타원형 문제도 같은 $h^{-2}$ 법칙을 따르며, 좋은 전처리기가 해야 하는 일이 바로 이 $h$ 의존성을 끊는 것이다.

같은 봉($N = 8$)에서 고정을 행 소거 대신 페널티 $\alpha \cdot k_{ii}$(대각 성분의 $\alpha$배)로 걸고 끝에 하중 1을 주면, 정답 끝변위 1에 대해 이렇게 된다.

$\alpha$	고정점이 미끄러진 양 = 끝변위 오차	조건수
$10^{2}$	$1.25\times10^{-3}$	$3.0\times10^{3}$
$10^{4}$	$1.25\times10^{-5}$	$2.9\times10^{5}$
$10^{6}$	$1.25\times10^{-7}$	$2.9\times10^{7}$
$10^{8}$	$1.25\times10^{-9}$	$2.9\times10^{9}$

구속 오차는 $1/\alpha$로 줄고 조건수는 $\alpha$로 는다. 배정밀도($\varepsilon \approx 10^{-16}$)에서 직접법이라면 유효숫자의 절반을 구속 오차와 반올림에 나눠 주는 $\alpha \approx \varepsilon^{-1/2} \approx 10^{8}$ 언저리가 흔히 쓰는 경험칙이고, 반복법은 조건수에 바로 반복 수가 묶이므로 훨씬 작은 값을 쓰거나 아예 행 소거를 택한다.

6. 대역폭과 절점 번호 매기기

희소 행렬이라도 어떻게 저장하고 푸느냐가 성능을 좌우한다. 직접법(가우스 소거, 촐레스키)에서는 절점 번호를 어떻게 매기느냐에 따라 0이 아닌 성분들이 대각선에서 얼마나 멀리 퍼지는지, 즉 대역폭(bandwidth)이 달라진다. 대역폭 $b$인 밴드 행렬을 촐레스키 분해하는 비용은 대략 $O(n b^2)$이므로, 대역폭을 줄이는 것이 곧 속도다. 절점 번호를 잘 재배열해 대역폭을 줄이는 알고리즘으로 커스힐-맥키(Cuthill-McKee, 1969)와 그 역순 버전(RCM)이 고전이다.⁴

더 근본적으로, 분해 과정에서 원래 0이었던 자리에 0이 아닌 값이 새로 채워지는 필-인(fill-in)을 최소화하는 재배열(예: 최소차수 순서화, 중첩절단)이 대규모 직접법의 핵심 기술이다. 차이는 규모로 보인다. $m \times m$ 2차원 격자($n = m^2$)를 행 순서로 번호 매기면 대역폭이 $m$이라 밴드 촐레스키가 $n b^2 = m^4 = n^2$ 연산이지만, 중첩절단은 $O(n^{3/2})$로 끝난다(George, 1973). $n = 10^6$이면 $10^{12}$ 대 $10^{9}$ 규모, 상수를 빼고도 세 자릿수다. 이런 재배열과 분해를 묶은 것이 희소 직접 솔버다. 반면 문제가 너무 커지면 직접법을 포기하고 크리로프 부분공간법 같은 반복법으로 넘어가는데, 이때는 전처리기의 품질이 수렴 속도를 결정한다. SPD 행렬에 대한 불완전 촐레스키 전처리와 켤레기울기법의 조합이 구조해석 반복법의 단골 메뉴다.

7. 관련 문서

• 유한요소법
• 구조해석
• 크리로프 부분공간법
• 전처리기
• 경계 조건
• 모드 해석
• 응력과 변형률 · 등매개변수 요소 · 갤러킨 방법
• 촐레스키 분해 · 켤레기울기법 · 희소 직접 솔버 · 희소행렬
• 조건수 · 비선형 구조해석 · 접촉 해석

8. Footnotes

1. 유한요소 코드의 겉모습은 화려한 3D 형상과 컬러 컨투어지만, 내부에서 벌어지는 일의 90%는 이 희소 행렬을 조립하고 분해하는 지루한 선형대수다. 그래픽은 결과 발표용, 강성행렬은 밥벌이용. ↩

2. “A점에 힘을 줘 B점에서 생긴 처짐 = B점에 같은 힘을 줘 A점에서 생긴 처짐”은 맥스웰(1864)의 상반정리이고, 베티(1872)가 이를 임의의 두 하중계 사이의 일의 상반성으로 일반화했다. 엄밀히 이 문장이 직접 말하는 것은 유연도 행렬의 대칭 $F_{ij}=F_{ji}$이고, 구속된 강성행렬은 그 역행렬($\mathbf{K}=\mathbf{F}^{-1}$)이라 대칭이 그대로 넘어온다. 160년 전 정리가 오늘날 메모리 절반을 아껴준다. 페널티 쪽 단골 무대는 접촉 해석이다. ↩

3. 페널티 값이 너무 작으면 구속이 헐거워 물체가 조금 움직이고, 너무 크면 조건수가 나빠져 반복법이 죽는다. 접촉 해석의 페널티 강성 튜닝과 정확히 같은 딜레마다. 세상 모든 페널티법은 결국 “얼마나 큰 수를 넣을 것인가”로 귀결된다. ↩

4. 커스힐-맥키 순서는 1969년에 나왔고, 그 순서를 뒤집기만 해도 필-인이 줄어든다는 것(RCM)은 앨런 조지가 1971년 학위논문에서 보였다. 지금도 대부분의 희소 행렬 라이브러리(MATLAB, SciPy 등)에 기본 탑재되어 있다. 반세기 넘게 현역인 걸 보면, 좋은 알고리즘엔 은퇴가 없다. ↩