제4부 기계 설계·구조·재료

4.8 진동과 구조 동특성: 모드 해석·공진 회피

4.8.1 이 장에서 배우는 것 / 선수 지식 / 특히 중요한 독자

이 장은 로봇 구조(링크·프레임·하우징)의 선형 진동을 다룬다. 구조는 질량·강성·감쇠로 이루어진 계이고, 모드 해석은 이 계가 외력 없이 진동하는 고유 진동수와 모드 형상을 일반화 고유값 문제로 구하는 일이다. 공진 회피는 가진 진동수(모터 회전, 기어 맞물림, 접촉 충격)가 구조의 고유 진동수에 가까워지지 않도록 설계하거나 운전 범위를 제한하는 일이다. 이 장은 에르미트 보 요소의 모드 해석과 수렴, 팁 질량(가반하중)과 뿌리 회전 강성이 1차 고유 진동수에 미치는 영향, 유효 모달 질량과 모드 선별, 레일리 감쇠와 주파수 응답, 반전력 대역폭에 의한 감쇠비 추정, 허용 증폭에서 역산한 진동수 분리 여유, 기저 가진 전달률, 캠벨 선도를 계산으로 확인한다. 관절 유연성과 링크 유연성을 포함한 로봇 동역학 모델은 3.9, 연속체 로봇의 동역학은 3.10, 재료 선정과 비강성은 4.4, 링크·관절의 강성 설계는 7.4, 서보 루프의 공진 억제는 5.9, 입력 성형은 12.11, 대형 팔의 처짐과 진동은 7.6, 유한요소 도구는 4.14, 로봇 소음은 26.26, 진동 증상의 현장 진단은 37.12이 맡는다.

배우는 것:

• 모드 해석의 구조 — $K\phi=\omega^2M\phi$, $M$-직교성, 모드 좌표의 분리, 비례 감쇠 $C=\alpha M+\beta K$에서 모드별 감쇠비.
• 보 요소 모델과 수렴 — 에르미트 3차 요소의 $K$·$M$, 고유값 수렴 차수 4, 특성 방정식의 해석해, 레일리 몫(상한)과 던컬리 합(하한).
• 모드 선별 — 유효 모달 질량과 정적 처짐의 모드 포착: 몇 개의 모드가 응답을 대표하는가.
• 감쇠와 주파수 응답 — 공진 배율 $1/(2\zeta)$, 반전력 대역폭 추정의 오차($\zeta^2$ 정도로 커짐).
• 공진 회피 — 허용 증폭 $a$에서 진동수비 한계 $r\le r_{lo}$ 또는 $r\ge r_{hi}$를 역산, 전달률 $T(\sqrt2)=1$, 캠벨 선도로 운전 속도 범위의 공진 교차 식별.

선수 지식: 3.1(강체 동역학: 관성과 렌치), 3.3(라그랑주 정식화), 1.1(행렬 분해: 촐레스키·고유값 계열).

특히 중요한 독자: 링크·하우징·프레임의 강성과 질량을 정하는 기계 설계자, 서보 대역폭과 구조 공진의 관계를 조정하는 제어 엔지니어(5.9), 진동 시험과 모달 시험을 계획하는 시험 담당자, 진동 증상을 진단하는 현장 엔지니어(37.12).

출처 방침. 기준 시점은 2026-11-05이다. 직접 읽은 자료: COMET-FEniCS 문서의 외팔보 모달 해석 예제(3차원 고체 모델의 고유 진동수와 유효 모달 질량 표, WebFetch 요약으로 읽음); Wikipedia "Dunkerley's method"의 서지 정보(Dunkerley 1894). 서지만 확인: Dunkerley(1894, Phil. Trans. R. Soc. A 185, 279–360, doi:10.1098/rsta.1894.0008), Ewins의 모달 시험 교재·Blevins의 고유 진동수 공식집·Clough & Penzien의 구조 동역학 교재(Crossref에는 서평만 등재되어 있어 책의 존재만 확인했고 판·본문은 읽지 못했다). 읽지 못한 것: 위 교재들의 본문, ISO 20816-1·ISO 7626 시리즈·IEC 60068-2-64·ISO 13373-1·MIL-STD-810의 본문(iso.org 접속 차단·유료 표준), IEC 60068-2-6은 판매처 목록의 소개문만 확인했다. 이 장의 이론(고유값 문제·모드 중첩·유효 질량·레일리 몫)은 교과서의 표준 결과이며 이 장이 직접 유도하거나 구현하고 독립 라이브러리와 해석해로 대조했다. 재료 물성·단면·감쇠비·가진 차수는 이 장이 정한 가정이며 실제 제품 값이 아니다.

4.8.2 개념 — 구조 진동은 질량·강성·감쇠의 고유값 문제이고 공진 회피는 가진 진동수와의 분리 여유를 허용 증폭에서 역산하는 일이다

결론부터 적으면 다음과 같다. 선형 구조의 자유 진동은 서로 $M$-직교하는 모드의 합으로 표현되고, 각 모드는 고유 진동수 $\omega_i$·모드 형상 $\phi_i$·감쇠비 $\zeta_i$를 갖는 1자유도 진동자처럼 거동한다. 외력의 진동수가 어떤 모드의 고유 진동수와 같으면 그 모드의 응답은 정적 응답의 $1/(2\zeta)$배까지 커진다. 로봇에서는 1차 모드가 대개 응답의 대부분을 대표하며(유효 모달 질량 비), 가반하중이 늘면 고유 진동수가 $1/\sqrt{M}$로 낮아져 같은 모터 회전수의 가진이 공진에 접근한다. 공진 회피 여유는 경험 값이 아니라 허용 증폭 $a$에서 닫힌 식으로 정해진다. 구조를 고정단 보로 모델링하면 관절 컴플라이언스를 무시하므로 1차 고유 진동수의 상한을 준다.

이 장의 계산 결과를 미리 적으면 다음과 같다. ① 알루미늄 관 링크(외경 60 mm·두께 3 mm·길이 0.8 m, 질량 $1.160$ kg)를 에르미트 요소로 이산화하면 1차 고유 진동수가 요소 수에 대해 차수 $3.9\sim4.0$으로 해석해 $89.844$ Hz에 수렴한다(40자리 계산의 차수 $3.998$). ② 팁 질량 $5$ kg(질량비 $4.31$)이면 1차 $20.760$ Hz이고 해석해와 같다. 레일리 몫 $20.761$ Hz는 상한, 던컬리 합 $20.745$ Hz는 하한이다. ③ 유효 모달 질량 비는 $0.6131$·$0.1883$·$0.0647$(해석 모드 형상의 적분과 $10^{-5}$ 안에서 같고, 공개된 3차원 고체 모델의 $0.6130$·$0.1885$·$0.0649$와 같은 값)이며, 끝 하중 정적 처짐의 $97.07\,\%$를 1차 모드가 포착한다. ④ 3자유도 모델의 반전력 추정은 감쇠비 $0.0254$·$0.0163$·$0.0135$를 $0.64\,\%$·$-0.59\,\%$·$0.01\,\%$ 안에서 복원한다. ⑤ 허용 배율 $1.25$는 진동수비 $r\le0.447$ 또는 $r\ge1.341$($\zeta=0.02$)을 요구하고, 전달률은 감쇠비와 무관하게 $r=\sqrt2$에서 $1$이다. ⑥ 가반하중 $5$ kg에서 모터 회전수 $1246$ rpm의 1차 성분이 1차 모드와 일치하고, 허용 회전수는 $557$ rpm 이하 또는 $1670$ rpm 이상이다. ⑦ 뿌리 회전 강성이 링크 굽힘 강성 $EI/L$의 $10$배이면 1차 고유 진동수는 고정단 값의 $87.5\,\%$이다. ⑧ scikit-fem·SciPy·mpmath·solve_ivp와 독립 대조했다.

4.8.3 원리

4.8.3.1 1자유도 진동자와 공진

질량 $m$·강성 $k$·점성 감쇠 $c$의 1자유도 계는 $\omega_n=\sqrt{k/m}$, $\zeta=c/(2\sqrt{km})$로 특성이 정해진다. 결론은 다음과 같다. 정현파 힘에 대한 변위 응답의 정적 대비 배율은 진동수비 $r=\omega/\omega_n$에서 다음이다.

$$\frac{\lvert X\rvert}{X_{st}}=\frac{1}{\sqrt{(1-r^2)^2+(2\zeta r)^2}},\qquad \frac{\lvert X\rvert}{X_{st}}\Big\rvert_{r=1}=\frac1{2\zeta}.$$

최대 배율은 $r=\sqrt{1-2\zeta^2}$에서 $1/(2\zeta\sqrt{1-\zeta^2})$이다. $\zeta=0.02$이면 공진에서 배율이 $25.0$이다(전달률은 $\sqrt{1+4\zeta^2}/(2\zeta)=25.02$, 표 ⑦). 자유 진동은 $e^{-\zeta\omega_nt}$로 감쇠하므로 $2\,\%$ 정착 시간은 $\ln(50)/(\zeta\omega_n)$이다. 이 장의 링크에 $5$ kg을 달면 $f_1=20.76$ Hz이고 $\zeta=0.01$·$0.02$·$0.05$에서 정착 시간이 $3.00$·$1.50$·$0.60$ s이다. 말단 진동의 정착 시간은 위치 결정 사이클 타임에 직접 더해지므로 고유 진동수와 감쇠비는 사이클 타임의 요소이다.

기저가 $y(t)=Y\sin\omega t$로 움직일 때 질량의 절대 변위 진폭비(전달률)는 다음이다.

$$T(r,\zeta)=\sqrt{\frac{1+(2\zeta r)^2}{(1-r^2)^2+(2\zeta r)^2}}.$$

$r=\sqrt2$에서 분자와 분모가 같아져 $T=1$이고 $r>\sqrt2$에서만 $T<1$(절연)이 된다. 감쇠를 키우면 공진 배율은 줄지만 $r>\sqrt2$에서의 절연은 나빠진다(표 ⑦: $r=3$에서 $\zeta=0.02$는 $0.126$, $\zeta=0.2$는 $0.193$).

4.8.3.2 다자유도 계의 모드 해석

$n$자유도 계 $M\ddot x+C\dot x+Kx=f$에서 $C=0$, $f=0$이면 $x=\phi\,e^{i\omega t}$를 대입해 일반화 고유값 문제 $K\phi=\omega^2M\phi$를 얻는다. $K$가 반정부호이고 $M$이 정부호이면 고유값 $\omega_i^2\ge0$이고 모드가 $\phi_i^TM\phi_j=\delta_{ij}$로 정규화된다(질량 정규화). 이때 $\phi_i^TK\phi_j=\omega_i^2\delta_{ij}$이다. $x=\Phi q$로 좌표를 바꾸면 방정식이 모드마다 분리된다.

$$\ddot q_i+2\zeta_i\omega_i\dot q_i+\omega_i^2q_i=\phi_i^Tf,\qquad \zeta_i=\frac{\alpha}{2\omega_i}+\frac{\beta\omega_i}{2}\ \ (C=\alpha M+\beta K).$$

비례 감쇠가 아닌 $C$는 모드를 분리하지 못하며 복소 모드가 필요하다. 주파수 응답은 모드 중첩으로 $H_{jk}(\omega)=\sum_i\phi_{ji}\phi_{ki}/(\omega_i^2-\omega^2+2i\zeta_i\omega_i\omega)$이다. 풀이 예제에서 이 식은 직접 역행렬과 $10^{-15}$ 수준에서 같다. 수치 해법은 $M=CC^T$(촐레스키)로 표준형 $C^{-1}KC^{-T}v=\omega^2v$로 변환한 뒤 대칭 행렬의 고유값 알고리즘(야코비·QR·란초스)을 적용하고 $\phi=C^{-T}v$로 돌려받는다. 큰 모델에서는 부분공간 반복이나 란초스로 낮은 모드 몇 개만 구한다(4.14, 3.10).

유효 모달 질량은 기저 가진에서 각 모드가 얼마의 질량을 움직이는지의 척도이다. 기저 방향 영향 벡터 $r$(강체 병진)과 고정 자유도와의 질량 결합 $c$에 대해 참여 계수 $\Gamma_i=\phi_i^T(Mr+c)$, 유효 질량 $m_{eff,i}=\Gamma_i^2$이다. 모드를 몇 개 남겨야 하는지는 누적 유효 질량의 목표값(해석 목적이 정한다)으로 정한다. 끝 하중의 정적 처짐은 $\sum_i\phi_i(L)^2/\omega_i^2$로 모드 중첩되며, 고차 모드를 자르면 이 합이 정적 처짐보다 작아지므로 정적 포착 비율을 함께 확인한다. 표 ③에 두 수치를 적었다.

4.8.3.3 오일러-베르누이 보와 해석해

균일 외팔보의 굽힘 진동 $EIw''''+\rho A\ddot w=0$의 모드는 $\lambda=\beta L$이 $\cos\lambda\cosh\lambda=-1$을 만족할 때 존재하며 $\lambda_n=1.8751,\ 4.6941,\ 7.8548$이다(이 장의 이분법 계산값과 SciPy brentq가 $9\times10^{-16}$ 안에서 같다). 고유 진동수는 다음이다.

$$f_n=\frac{\lambda_n^2}{2\pi}\sqrt{\frac{EI}{\rho AL^4}}.$$

끝에 집중 질량 $M$이 있으면 특성 방정식이 $1+\cos\lambda\cosh\lambda+\lambda\mu(\cos\lambda\sinh\lambda-\sin\lambda\cosh\lambda)=0$($\mu=M/(\rho AL)$)로 바뀐다. 이 식이 맞음은 FE 32개 모델과 6개 팁 질량에서 $10^{-8}$ 안에서 일치하는 것으로 확인했다(표 ②). 에르미트 3차 요소(요소당 두 절점의 처짐과 기울기)는 요소 길이 $h$에서 $K_e=\frac{EI}{h^3}[\,\cdot\,]$, $M_e=\frac{\rho Ah}{420}[\,\cdot\,]$의 $4\times4$ 행렬을 갖는다(구현은 beam_mk, 스크립트의 행렬 값은 scikit-fem 조립과 $2.3\times10^{-12}$에서 같다). 고유값의 오차는 요소 길이의 4제곱에 비례하며(표 ①의 관측 차수 $3.9\sim4.0$) 상대 오차가 $10^{-8}$에 이르면 부동소수 반올림이 섞인다. 요소 32개에서 이 장의 야코비가 낸 상대 오차 $1.03\times10^{-8}$은 SciPy의 $8.6\times10^{-9}$, 40자리 계산의 $8.18\times10^{-9}$와 약 $2\times10^{-9}$ 차이가 나는데 이것이 반올림 오차이며, 40자리 계산의 차수 $\log_2(e_{16}/e_{32})$는 $3.998$이다.

손으로 하는 추정에는 두 근사를 사용한다. 레일리 몫 $R(\psi)=\psi^TK\psi/\psi^TM\psi$는 임의의 시험 형상 $\psi$에서 $\omega_1^2$ 이상이며 시험 형상이 참 모드에 가까울수록 정확하다. 끝 하중의 정적 처짐 형상 $w(x)\propto3x^2L-x^3$을 사용하면 $\omega_1^2\le3EI/\big(L^3(M+\tfrac{33}{140}m_b)\big)$이다($m_b=\rho AL$; 계수 $33/140=0.235714$는 $\int w^2dx/w(L)^2$의 적분이며 quad와 같다). 던컬리 합 $1/\omega_1^2\approx1/\omega_{beam}^2+1/\omega_{tip}^2$($\omega_{tip}^2=3EI/(ML^3)$)는 $\omega_1$의 하한을 준다(Dunkerley 1894). 이 장의 6개 팁 질량에서 레일리 몫은 항상 정확한 값 이상, 던컬리 합은 이하였다(표 ②). 두 값 사이에 참값이 있으므로 설계 초기에는 두 식으로 $f_1$의 범위를 즉시 얻는다.

4.8.3.4 경계 강성과 접합부 — 고정단 가정은 상한이다

실제 로봇 링크는 관절 컴플라이언스(감속기·베어링·체결부)로 지지되므로 뿌리 회전 강성 $k_\theta$가 유한하다. 이 장은 고정단을 회전 스프링 $k_\theta=\kappa\,EI/L$로 바꾼 모델을 풀었다(표 ⑨). 끝 하중이 지배적인 경우($M\gg m_b$) 정적 컴플라이언스가 $L^3/(3EI)+L^2/k_\theta$의 합이므로 $f_1\approx f_c/\sqrt{1+3/\kappa}$($f_c$는 고정단 값)가 성립한다. $\kappa=10$이면 $0.877$, $\kappa=20$이면 $0.933$이며 FE 해와 $\kappa=5\sim20$에서 $0.4\,\%$ 이내, $\kappa=0.5$에서 $0.9\,\%$ 이내로 같다. 결론은 다음이다. 근사식에서 $f_1/f_c\ge0.9$이려면 $\kappa\ge12.8$이므로 관절 강성이 링크 굽힘 강성 $EI/L$의 약 13배 이상이어야 고정단 모델의 오차가 $10\,\%$ 아래로 내려가며, 그렇지 않으면 관절 모델(3.9)을 포함해야 한다. 볼트 체결부의 강성과 감쇠는 해석이 아니라 시험으로 보정하는 것이 일반적이다. 접합부의 마찰과 미세 슬립은 비선형이며 이 장은 선형 모델에 한정한다.

4.8.3.5 감쇠비 추정과 모달 시험

시험에서 감쇠비는 주파수 응답 함수(FRF)의 반전력 대역폭으로 추정한다. 봉우리 진폭의 $1/\sqrt2$가 되는 두 진동수 $\omega_1<\omega_2$에서 $\hat\zeta=(\omega_2-\omega_1)/(2\omega_{peak})$이다. 이 식은 $\zeta$가 작을 때만 정확하다. 1자유도 계의 표 ⑤에서 오차는 $\zeta=0.01$에서 $0.04\,\%$, $0.05$에서 $0.51\,\%$, $0.1$에서 $2.1\,\%$, $0.2$에서 $9.1\,\%$, $0.3$에서 $25.6\,\%$로 $\zeta^2$ 정도로 증가한다. 다자유도 계에서는 모드가 가까우면 이웃 모드의 꼬리가 봉우리를 변형한다. 3자유도 예에서 모드 1·2(고유 진동수 비 $0.57$)의 추정 오차가 $+0.64\,\%$·$-0.59\,\%$이고 유효 질량이 작은 모드 3($0.6\,\%$)은 봉우리가 뚜렷해 $+0.01\,\%$였다. 모드 형상의 시험·해석 상관은 모달 보증 기준(MAC)으로 하며 정의는 $\mathrm{MAC}=(\phi_a^T\phi_b)^2/\big((\phi_a^T\phi_a)(\phi_b^T\phi_b)\big)$이다. 이 장은 MAC의 계산 예와 시험 장비(임팩트 해머·가진기·가속도계) 설정을 다루지 않았다 (확인 필요).

4.8.3.6 공진 회피 — 허용 증폭에서 진동수 여유 역산

공진 회피의 첫째 기준은 가진 진동수와 고유 진동수의 비이다. 허용 배율 $a$에서 $1/\sqrt{(1-r^2)^2+(2\zeta r)^2}\le a$인 영역은 $r\le r_{lo}$ 또는 $r\ge r_{hi}$이다. 감쇠가 작으면($\zeta\to0$) $r_{lo}=\sqrt{1-1/a}$, $r_{hi}=\sqrt{1+1/a}$이다. $a=1.25$이면 $0.447$·$1.342$, $a=1.5$이면 $0.577$·$1.291$, $a=2$이면 $0.707$·$1.225$이다. 표 ⑥은 $\zeta$가 $0.01$에서 $0.1$까지 커져도 한계가 $1\sim2\,\%$만 이동함을 보인다. 이 여유는 변위 배율의 정현파 정상 응답에 대한 것이며 과도 응답(공진 통과)과 비정현 가진은 포함하지 않는다.

둘째 기준은 가진 진동수가 운전 조건에 따라 움직인다는 점이다. 모터 회전수 $N$[rpm]의 $k$차 성분은 $f=kN/60$ Hz이고, 가진원은 이 선이 고유 진동수를 지나는 회전수에서 공진한다. 이 선들을 한 평면에 그린 것이 캠벨 선도(도해 3)이다. 가반하중이 바뀌면 고유 진동수가 이동하므로 가반하중 범위 전체를 띠로 표시한다. 가진 성분의 차수는 부품 형식이 정한다. 모터의 회전 주파수(1차), 극수·슬롯 수에 따른 토크 리플, 기어 맞물림(잇수 $\times$ 회전수), 파동기어 파형 발생기의 2차 성분 같은 값을 감속기와 모터 해설(5.10, 5.5)에서 확인해 넣어야 하며 이 장의 1·2·6차는 가정이다 (확인 필요).

대응 수단은 네 가지로 나뉜다. 구조 변경(강성 증가·질량 감소로 고유 진동수 이동), 감쇠 추가(구속층·점탄성 재료), 가진원 변경(운전 속도 범위 제한·속도 건너뛰기·가속 프로파일), 제어(서보 노치 필터와 입력 성형, 5.9·12.11). 금속 재료의 비강성 $E/\rho$는 알루미늄($70\,\mathrm{GPa}/2700\,\mathrm{kg/m^3}$)과 강($200\,\mathrm{GPa}/7850\,\mathrm{kg/m^3}$) 모두 약 $26\times10^{6}\ \mathrm{m^2/s^2}$로 비슷하다(이 장의 가정값에서 계산, 물성은 (확인 필요)). 재료 교체보다 단면 형상(관성 모멘트 증가)이 고유 진동수를 높이는 주된 수단이라는 뜻이다(4.4).

4.8.4 알고리즘

입력: 구조 모델(요소·재료·경계 조건), 접합부 강성, 가진원 목록(차수·회전수 범위), 허용 배율 a, 감쇠비 추정 zeta
1. 모델: 요소로 K, M을 조립한다. 관절·체결부는 스프링으로 넣고 시험으로 보정한다.
2. 모드 해석: M = C C^T, 표준형 A = C^-1 K C^-T, 대칭 고유값 해(야코비·QR·란초스), phi = C^-T v
3. 모드 선별: 유효 모달 질량 누적과 정적 처짐 포착 비율이 해석 목적의 목표값에 이르는 모드까지 남긴다.
4. 여유: a 와 zeta 에서 r_lo, r_hi 를 구해 f_exc / f_n 이 허용 영역 안인지 점검한다. 가반하중 범위 전체의 f_n 띠로 점검한다.
5. 캠벨: 각 가진 성분 f = k N / 60 과 f_n 띠의 교차 회전수를 구하고 운전 회전수를 r 한계 밖으로 둔다.
6. 위반 시: 구조 변경, 감쇠 추가, 가진 변경, 제어(5.9, 12.11) 중 하나를 고르고 2로 돌아간다.
7. 시험: 임팩트 해머나 가진기로 FRF를 얻어 고유 진동수·감쇠비·모드 형상을 식별하고 모델을 보정한다.
출력: 고유 진동수 표, 유효 질량 표, 허용 운전 회전수 구간, 보정된 모델

4.8.5 풀이 예제 — 링크와 팁 질량의 모드 해석, 감쇠 추정, 공진 회피

결론부터 적으면 다음과 같다. ① 에르미트 요소 모델의 고유 진동수는 해석해로 차수 4로 수렴하고 팁 질량 $0\sim10$ kg에서 해석해와 $10^{-8}$ 안에서 같다. ② 1차 모드는 보 단독에서 유효 질량의 $61.3\,\%$, $5$ kg 팁 질량에서 $90.9\,\%$를 대표하고 끝 하중 정적 처짐의 $97.1\,\%$ 이상을 포착한다. ③ 1차 고유 진동수는 가반하중 $0\to10$ kg에서 $89.84\to14.87$ Hz로 낮아진다. ④ 같은 링크가 뿌리 회전 강성 $EI/L$의 $10$배로 지지되면 고정단 대비 $12.5\,\%$ 낮다. ⑤ 허용 배율 $1.25$에서 $5$ kg 팁 질량의 허용 회전수는 $557$ rpm 이하 또는 $1670$ rpm 이상이다.

구현. beam_mk(보 요소 $K$·$M$·질량 결합), chol·jacobi·geig(일반화 고유값), eff_mass(유효 모달 질량), lam_roots(특성 방정식의 근), rayleigh_w2·dunkerley_w2(근사). 주파수 응답·반전력·진동수비 한계는 코어 밖의 보조 함수이다. 표 생성·도해·캐시는 스크립트에 있다.

from math import sqrt, pi, cos, sin, cosh, sinh, log

def beam_mk(n, L, EI, rA, M=0.0, kth=None):
    """외팔보 에르미트 요소 n개. 고정단 자유도 제거 후
    (K, M, c): c는 고정단 처짐 자유도와 질량 결합.
    kth가 있으면 고정단 대신 회전 스프링 kth로 지지"""
    h = L / n
    N = 2 * n + 2
    K = [[0.0] * N for _ in range(N)]
    Mm = [[0.0] * N for _ in range(N)]
    a, b = EI / h ** 3, rA * h / 420
    h2 = h * h
    k4 = [[12, 6 * h, -12, 6 * h],
          [6 * h, 4 * h2, -6 * h, 2 * h2],
          [-12, -6 * h, 12, -6 * h],
          [6 * h, 2 * h2, -6 * h, 4 * h2]]
    m4 = [[156, 22 * h, 54, -13 * h],
          [22 * h, 4 * h2, 13 * h, -3 * h2],
          [54, 13 * h, 156, -22 * h],
          [-13 * h, -3 * h2, -22 * h, 4 * h2]]
    for e in range(n):
        for i in range(4):
            for j in range(4):
                K[2 * e + i][2 * e + j] += a * k4[i][j]
                Mm[2 * e + i][2 * e + j] += b * m4[i][j]
    Mm[N - 2][N - 2] += M
    d = 2
    if kth is not None:
        K[1][1] += kth
        d = 1
    cut = lambda A: [r[d:] for r in A[d:]]
    return cut(K), cut(Mm), [r[0] for r in Mm[d:]]

def chol(A):
    n = len(A)
    C = [[0.0] * n for _ in range(n)]
    for i in range(n):
        for j in range(i + 1):
            s = A[i][j]
            s -= sum(C[i][k] * C[j][k] for k in range(j))
            C[i][j] = sqrt(s) if i == j else s / C[j][j]
    return C

def fsub(C, b):
    y = []
    for i, row in enumerate(C):
        s = b[i] - sum(row[k] * y[k] for k in range(i))
        y.append(s / row[i])
    return y

def bsub(C, b):
    n = len(C)
    x = [0.0] * n
    for i in range(n - 1, -1, -1):
        s = b[i]
        s -= sum(C[k][i] * x[k] for k in range(i + 1, n))
        x[i] = s / C[i][i]
    return x

def jacobi(A, sweeps=60):
    """대칭 행렬의 (고유값, 고유벡터 열) 오름차순"""
    n = len(A)
    a = [r[:] for r in A]
    V = [[float(i == j) for j in range(n)] for i in range(n)]
    for _ in range(sweeps):
        off = sum(a[i][j] ** 2 for i in range(n)
                  for j in range(i))
        dg = sum(a[i][i] ** 2 for i in range(n))
        if off < 1e-28 * dg:
            break
        for p in range(n - 1):
            for q in range(p + 1, n):
                if a[p][q] == 0.0:
                    continue
                th = (a[q][q] - a[p][p]) / (2 * a[p][q])
                t = 1.0 if th >= 0 else -1.0
                t /= abs(th) + sqrt(th * th + 1)
                c = 1 / sqrt(t * t + 1)
                s = t * c
                for k in range(n):
                    x, y = a[k][p], a[k][q]
                    a[k][p] = c * x - s * y
                    a[k][q] = s * x + c * y
                for k in range(n):
                    x, y = a[p][k], a[q][k]
                    a[p][k] = c * x - s * y
                    a[q][k] = s * x + c * y
                for k in range(n):
                    x, y = V[k][p], V[k][q]
                    V[k][p] = c * x - s * y
                    V[k][q] = s * x + c * y
    od = sorted(range(n), key=lambda i: a[i][i])
    return ([a[i][i] for i in od],
            [[V[k][i] for i in od] for k in range(n)])

def geig(K, M):
    """K phi = w2 M phi. 질량 정규화 모드(행)를 반환"""
    n = len(K)
    C = chol(M)
    R = range(n)
    cols = [fsub(C, [K[i][j] for i in R]) for j in R]
    A = [fsub(C, [cols[j][i] for j in R]) for i in R]
    A = [[(A[i][j] + A[j][i]) / 2 for j in range(n)]
         for i in range(n)]
    w2, V = jacobi(A)
    Phi = [bsub(C, [V[k][j] for k in R]) for j in R]
    return w2, Phi

def eff_mass(Phi, M, r, c=None):
    """모드별 유효 모달 질량. 영향 벡터는 M r + c"""
    n = len(r)
    Mr = [sum(M[i][k] * r[k] for k in range(n))
          for i in range(n)]
    if c is not None:
        Mr = [u + v for u, v in zip(Mr, c)]
    return [sum(p * q for p, q in zip(f, Mr)) ** 2
            for f in Phi]

def lam_roots(mu, nroot=3):
    """특성식 1 + cos l cosh l + l mu g(l)의 근,
    g = cos l sinh l - sin l cosh l"""
    def f(l):
        g = cos(l) * sinh(l) - sin(l) * cosh(l)
        return 1 + cos(l) * cosh(l) + l * mu * g
    out, x = [], 0.01
    while len(out) < nroot:
        y = x + 0.005
        if f(x) * f(y) < 0:
            lo, hi = x, y
            for _ in range(100):
                m = (lo + hi) / 2
                if f(lo) * f(m) > 0:
                    lo = m
                else:
                    hi = m
            out.append((lo + hi) / 2)
        x = y
    return out

def rayleigh_w2(EI, L, mb, M):
    """끝 하중 처짐 형상의 레일리 몫 3EI/(L^3(M+33mb/140))"""
    return 3 * EI / (L ** 3 * (M + 33 / 140 * mb))

def dunkerley_w2(EI, L, rA, M):
    """던컬리 합 1/w2 = 1/w2(보) + 1/w2(팁 질량), 하한"""
    wb2 = lam_roots(0.0, 1)[0] ** 4 * EI / (rA * L ** 4)
    wm = M * L ** 3 / (3 * EI)
    return 1 / (1 / wb2 + wm)

표 ① 요소 수에 따른 외팔보 고유 진동수와 수렴.

요소 $2$개에서 4개, 8개, 16개로 갈 때 1차 고유 진동수의 관측 차수가 $3.30$·$3.89$·$3.97$·$3.99$로 4에 가까워진다. 요소 32개의 $3.67$은 상대 오차가 $10^{-8}$에 이르러 반올림 오차가 섞인 것이며 40자리 계산은 $3.998$이다(라이브러리 대조). 높은 모드일수록 같은 요소 수에서 오차가 크다(요소 $4$개의 2차 $1.2\times10^{-3}$, 3차 $7.7\times10^{-3}$).

표 ② 팁 질량(가반하중)에 따른 1차 고유 진동수.

팁 질량이 $0\to10$ kg으로 늘면 $f_1$이 $89.84\to14.87$ Hz로 낮아지고 질량비 $M/m_b\ge4$에서 레일리 몫의 오차가 $4\times10^{-5}$ 이하이다. 보 단독($M=0$)에서는 레일리 몫이 $+1.47\,\%$로 가장 크다. 던컬리 합은 팁 질량이 있을 때 $-0.04\sim-0.34\,\%$이다. 도해 1은 처음 3개 모드의 형상이며 2차 모드는 $x/L=0.783$에서 절점을 갖는다(32개 요소 모드 형상의 선형 보간).

표 ③ 유효 모달 질량과 끝 하중 정적 처짐의 포착(요소 32개).

보 단독에서 모드 4개의 누적 유효 질량이 $0.899$이고 팁 질량 5 kg이면 모드 1이 $0.9095$, 모드 4개가 $0.980$이다. 끝 하중의 정적 처짐은 모드 1개로 $97.07\,\%$, 2개로 $99.54\,\%$를 포착한다. 응답이 아닌 질량 관점의 누적($0.80$)은 정적 처짐 포착($0.995$)보다 느리게 오른다. 목적이 정하는 모드 수가 다르다는 뜻이며, 가반하중이 클수록 모드 1 하나가 대표하는 비중이 커진다. 공개된 3차원 고체 모델(COMET-FEniCS, $200\times6\times11$ 격자의 1차 라그랑주 요소)의 굽힘 모드 유효 질량 비 $0.6130$·$0.1885$·$0.0649$와 이 장의 보 요소 $0.6131$·$0.1883$·$0.0647$이 같은 값이다.

표 ④ 3자유도 모델: 레일리 감쇠와 반전력 추정.

표 ⑤ 1자유도 계의 반전력 추정 오차.

레일리 감쇠 $\alpha=15$, $\beta=10^{-5}$에서 모드 1~3의 감쇠비가 $0.0254$·$0.0163$·$0.0135$이다. 참 감쇠비는 상태공간 고유값에서 구한 값과 $1.2\times10^{-14}$ 안에서 같다. 반전력 추정은 모드 1·2에서 $+0.64\,\%$·$-0.59\,\%$이다. 1자유도에서 $\zeta\ge0.1$이면 오차가 $2\,\%$를 넘으므로 감쇠가 큰 구조의 FRF에는 곡선 맞춤(모드 중첩 식으로 $\omega_i$·$\zeta_i$·잔류 항을 함께 맞추기)이 필요하다. 모드 3의 유효 질량이 $0.6\,\%$로 작아도 말단 응답 $H_{33}$에서는 봉우리가 뚜렷하게 나타난다. 모드 중첩 FRF와 직접 해는 $1.7\times10^{-15}$ 안에서 같다.

표 ⑥ 허용 배율에서 역산한 진동수비 한계.

표 ⑦ 기저 가진 전달률 $T(r,\zeta)$.

$\zeta=0.02$에서 허용 배율 $1.25$이면 $r\le0.4474$ 또는 $r\ge1.3410$이다. 감쇠를 $0.1$로 키워도 $0.4529$·$1.3247$로 이동폭이 작다. 전달률은 $r=\sqrt2$에서 모든 감쇠비에 대해 $1.000$이고 $r=2$에서 $0.334\,(\zeta=0.02)\to0.412\,(\zeta=0.2)$로 감쇠가 클수록 절연이 나빠진다. 기저 가진의 정상 응답을 solve_ivp로 적분한 진폭은 이 식과 $9.5\times10^{-7}$ 안에서 같다.

표 ⑧ 캠벨 선도와 허용 회전수.

가반하중 $2\to10$ kg에서 1차 고유 진동수가 $31.61\to14.87$ Hz로 낮아지면 1차 가진 성분의 공진 회전수가 $1897\to892$ rpm으로 이동한다. 도해 3의 띠가 가반하중 범위의 고유 진동수 대역이고, 1차 성분 직선이 이 띠를 지나는 회전수 $892\sim1897$ rpm이 공진 위험 구간이다. 2차 성분은 같은 대역을 $446\sim948$ rpm에서, 6차 성분은 $149\sim316$ rpm에서 지난다. 1차 성분만 보면 허용 회전수는 가반하중 $5$ kg에서 $557$ rpm 이하 또는 $1670$ rpm 이상이다. 가반하중이 $2\sim10$ kg으로 변하는 제품이라면 1차 성분의 공진 회피가 어려워지고(허용 구간 $399$ 이하 또는 $2543$ 이상) 운전 범위를 제한하거나 가반하중별 속도 한계를 두어야 한다.

표 ⑨ 뿌리 회전 강성에 따른 1차 고유 진동수(팁 질량 5 kg).

라이브러리 독립 대조

대조 결과 이상 징후는 없었다. 이 장의 야코비는 SciPy eigh와 $2.4\times10^{-9}$ 안에서 같고, 두 구현의 차이는 반올림 오차 수준이다. 1자유도 반전력 추정의 격자 보간은 brentq와 $2.2\times10^{-4}$ 안에서 같다(격자 간격에 의한 오차). 처음에는 반전력 식의 분모를 고유 진동수 $\omega_n$으로 두고 brentq와 비교해 $9\times10^{-2}$ 어긋났다. 이 장은 시험에서 알 수 있는 봉우리 진동수 $\omega_{peak}=\omega_n\sqrt{1-2\zeta^2}$로 나누므로 정의를 맞춘 뒤 같아졌다.

한계. 선형 모델이며 비례 감쇠만 다루었다. 접합부의 마찰·슬립·간극, 비선형 강성, 회전체의 자이로 효과, 공진 통과의 과도 응답, 진동 피로는 계산하지 않았다. 링크 단면은 균일하며 관절 모듈의 질량 분포·감속기 강성·케이블 영향은 넣지 않았다. 재료 물성·감쇠비·가진 차수는 이 장의 가정이다. 모달 시험(임팩트 해머·가진기)과 MAC 상관은 계산하지 않았다 (확인 필요).

견해: 설계 초기에는 레일리 몫과 던컬리 합으로 $f_1$의 범위를 구하고, 가반하중 범위 전체의 $f_1$ 띠에 대해 캠벨 선도로 운전 회전수를 점검한 뒤, 관절 강성을 포함한 모델을 시험 FRF로 보정하는 순서가 맞다고 본다. 근거는 이 장의 계산에서 고정단 가정이 관절 강성 $\kappa=10$에서 $f_1$을 $12.5\,\%$ 높게 주고(표 ⑨), 가반하중이 $2\to10$ kg으로 변하면 공진 회전수가 $1897\to892$ rpm으로 움직이며(표 ⑧), 반전력 추정의 오차가 감쇠비에 따라 $0.04\to25.6\,\%$로 달라지기 때문이다. 이 순서가 실제 로봇 팔의 시험 결과와 일치하는지는 이 장이 검증하지 않았다.

4.8.6 설계 트레이드오프 — 공진 회피 전략

공진 회피 전략은 다음 중에서 고른다.

• ① 구조 변경(주파수 분리). 강성을 올리고 질량을 줄여 고유 진동수를 가진 대역 밖으로 옮긴다. 효과가 직접적이고 영구적이나 질량·원가·형상 제약과 충돌하며(4.4, 7.4), 비강성이 비슷한 금속 사이의 재료 교체로는 한계가 있다.
• ② 감쇠 추가. 구속 감쇠층·점탄성 재료·마찰 감쇠로 공진 배율을 줄인다. 분리가 어려운 구조에 사용되나 감쇠비는 해석으로 예측하기 어렵고 시험으로 얻어야 하며, 전달률 관점에서는 고주파 절연이 나빠진다(표 ⑦).
• ③ 제어와 운전 회피. 노치 필터·입력 성형·속도 건너뛰기·가속 프로파일로 가진을 구조에 맞춘다(5.9, 12.11). 하드웨어 변경이 없으나 고유 진동수가 가반하중과 자세로 변하는 경우 강건성이 문제가 된다.

견해: 로봇 팔·관절 모듈에는 ①로 1차 고유 진동수를 가진 대역 밖에 두고 ③으로 나머지를 보강하는 것이 합리적이라고 본다. 근거는 이 장의 계산에서 가반하중이 바뀌면 공진 회전수가 $1897\to892$ rpm으로 이동해(표 ⑧) 운전 회피만으로는 범위가 좁고, 감쇠는 해석으로 예측하기 어려워 반전력 추정에도 오차가 있다는 점이다. 세 전략을 같은 제품에 적용해 비용·성능을 비교한 증거는 이 장에 없다.

4.8.7 응용 축 비교

이 표는 이 책의 정리이며 제품별 요구는 다르다. 응용별 수치는 이 장이 확인하지 않았다.

4.8.8 관련 규격·프로토콜

표준 번호 외에 판·조항은 이 장이 단정하지 않는다. 규격 체계 전반은 27.1이다.

4.8.9 관련 규제·인증 요구와의 연결

구조 진동을 직접 규정하는 로봇 규제는 이 장이 확인하지 못했다 (확인 필요). 간접 연결은 다음과 같다. 구조 공진에서 비롯한 제어 불안정이나 위치 오차는 안전 기능이 가정하는 정지 거리와 속도 한계를 흔들 수 있으므로 위험성 평가(26.4)가 고려해야 하는 요인이다. 소음(26.26)과 진동 노출은 사람 접촉 응용의 환경 요건과 닿아 있다. 운송·사용 환경의 진동·충격은 환경 시험(IEC 60068 계열, MIL-STD-810)으로 확인하며 시험 항목 선정은 제품 규격이 정한다(4.13, 27.26). 설계 변경으로 고유 진동수가 이동하면 기존 시험 증빙을 다시 확인해야 한다.

4.8.10 한계·유의점·미해결

• 선형·비례 감쇠. 접합부 마찰·간극·비선형 강성은 반영하지 않았다. 감쇠비는 해석이 아니라 시험으로 얻는다.
• 단순 모델. 균일 단면 보와 3자유도 집중 모델이며 관절 모듈·케이블·감속기 강성을 포함한 실제 팔의 모드는 구하지 않았다(3.9, 4.14).
• 자세 의존. 로봇의 고유 진동수는 자세와 가반하중에 따라 이동한다. 이 장은 가반하중만 변수로 했고 자세 의존은 3.9에 위임했다.
• 과도 응답. 정상 응답의 배율만 다루었다. 공진 통과 때의 최대 응답, 충격 응답 스펙트럼, 랜덤 진동은 계산하지 않았다.
• 회전체. 모터·감속기 회전체의 자이로 효과로 인한 모드 분리(전진·후진 모드)는 다루지 않았다.
• 미해결. 관절 모듈의 실제 비선형 감쇠 모델, 모달 시험과 모델 보정의 정량 절차, 가진 차수의 부품별 목록(모터·감속기)은 이 장이 확인하지 않았다.

4.8.11 체크포인트·연구 질문·출처

체크포인트:

1. 외팔보 $f_1=\frac{3.516}{2\pi}\sqrt{EI/(\rho AL^4)}$에서 $3.516=1.8751^2$임을 보이고, 이 장의 링크(표 입력)에서 $89.84$ Hz를 확인하라.
2. 끝 하중의 정적 처짐 형상 $w=(3x^2L-x^3)/(2L^3)$로 레일리 몫의 계수 $33/140$을 유도하라.
3. 레일리 몫이 항상 $\omega_1^2$ 이상임을 $\psi=\sum c_i\phi_i$ 전개로 증명하라.
4. 허용 배율 $a$에서 $\zeta\to0$일 때 $r_{lo}=\sqrt{1-1/a}$, $r_{hi}=\sqrt{1+1/a}$임을 유도하라.
5. 전달률이 감쇠비와 무관하게 $r=\sqrt2$에서 $1$임을 보이라.
6. 가반하중 $m_p$를 변수로 $f_1\approx\frac1{2\pi}\sqrt{3EI/\big(L^3(m_p+\tfrac{33}{140}m_b)\big)}$이 $3\%$ 안에서 맞는 질량비 범위를 표 ②에서 찾아라.

연구 질문:

• 관절 강성과 링크 강성이 같은 차수인 팔에서 고유 진동수의 자세 의존은 가반하중 의존과 비교해 어느 정도인가(이 장은 계산하지 않았다)?
• 접합부(볼트·감속기)의 마이크로 슬립 감쇠를 가진 진폭의 함수로 어떻게 모델링하고, 이것이 모드 해석의 감쇠비 추정에 어떤 영향을 주는가?
• 서보 대역폭·구조 모드·입력 성형을 함께 고려한 최소 사이클 타임 설계에서 허용 배율 $a$는 어떻게 정해지는가?

1차 출처(조회 2026-11-05).

• COMET-FEniCS, "Cantilever modal analysis" 예제 문서(comet-fenics.readthedocs.io) — 3차원 고체 모델의 고유 진동수와 유효 모달 질량 표를 읽었다(WebFetch 요약). 2차 출처.
• S. Dunkerley, "On the whirling and vibration of shafts," Phil. Trans. R. Soc. A 185, 279–360, 1894, doi:10.1098/rsta.1894.0008 — 서지만 확인.
• D. J. Ewins의 모달 시험 교재·R. D. Blevins의 Formulas for Natural Frequency and Mode Shape·R. W. Clough & J. Penzien의 Dynamics of Structures — 서지만 확인(Crossref에는 서평 형태로만 등재, 책의 판·본문은 읽지 못했다).
• IEC 60068-2-6 — 판매처 목록의 소개문만 읽었다. IEC 60068-2-64·IEC 60068-2-27·ISO 20816-1·ISO 7626·ISO 13373-1·MIL-STD-810 — 번호만, 본문 미확인.
• scikit-fem 12.0.2(ElementLineHermite), SciPy 1.18.1(linalg.eigh·optimize.brentq·integrate.quad·integrate.solve_ivp·linalg.eig), numpy 2.5.4, mpmath 1.4.1.

다음 장. 4.9는 순응 기구(노치·판스프링 유연 힌지·의사강체 모델·정밀 스테이지)를 다룬다. 의도적으로 탄성 변형을 사용하는 기구의 강성·정밀도와 이 장의 고유 진동수는 같은 $k/m$ 관계로 이어진다.