제3부 동역학

3.10 연속체·소프트 로봇 동역학: 코세라 막대·FEM 소프트 모델

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

이 장은 몸 전체가 연속적으로 변형하는 로봇(연속체 로봇·소프트 로봇)의 동역학, 곧 관성력·탄성력·구속력이 운동을 정하는 층을 다룬다. 기구학 장(2.16)은 형상을 곡률 분포의 기하로 적고 끝 하중 하의 평형(엘라스티카)까지 다루었다. 이 장은 같은 코세라 막대를 시간 영역으로 확장하여 운동량 평형식, 반이산화, 시간 적분의 안정 한계를 계산으로 확인하고, 3차원 유한요소(FEM)와 일정 곡률(PCC) 환원 모델을 같은 외팔보 문제로 비교한다. PCC의 매개변수화와 제어는 25.3, SOFA·PyElastica 같은 시뮬레이터와 역모델 QP는 25.4, 연속체 로봇의 설계·구동은 25.6, 접촉은 3.12, 구속 동역학의 일반 이론은 3.11에서 다룬다.

배우는 것:

• 코세라 막대의 동적 평형 — 선운동량·각운동량 평형, 신장·전단·굽힘 구성 방정식, 선형화와 외팔보 진동수.
• 반이산화 — 요소 변형률 에너지와 집중 질량으로 $M\ddot q+\nabla\Pi(q)=Q$를 만들고 기울기·헤시안의 닫힌 식을 확인한다.
• 시간 적분 — 명시(속도 베를레)의 안정 한계 $\Delta t_c=2/\omega_{\max}$와 암시(뉴마크)의 위상·에너지 거동.
• 환원 — 일정 곡률 모델을 레일리-리츠 근사로 위치시키고 진동수의 상계와 수렴을 계산한다.
• 3차원 FEM — 네오-훅 구성식, 접선 강성, 체적 잠금, 구속(라그랑주 승수)에 의한 구동.

선수 지식: 2.16(코세라 운동학·엘라스티카), 3.3(라그랑주 정식화), 1.18(ODE 적분기), 3.9(유연 링크의 가정 모드법), 1.1(선형대수).

특히 중요한 독자: 연속체·소프트 로봇의 동적 모델을 구현하거나 시뮬레이터 결과를 해석하는 개발자와 연구자, 의료용 연속체 로봇의 모델 신뢰성을 검토하는 엔지니어.

출처 방침. 기준 시점은 2026-11-05이다. 직접 읽은 자료: ① Armanini·Boyer·Mathew·Duriez·Renda, "Soft Robots Modeling: A Structured Overview"(arXiv:2112.03645 v2; IEEE T-RO 39(3):1728–1748, 2023)의 3~7절 중 이 장에서 인용한 부분(연속체 역학, 코세라 막대 평형식, FEM, 수치 해법, 일정 곡률 모델). ② Gazzola 외, "Forward and inverse problems in the mechanics of soft filaments"(Royal Society Open Science 5(6):171628, 2018)의 2~3절(지배 방정식, 시간 이산화, 티모셴코 외팔보 검증; 추출 본문에서 수식이 빠져 서술 범위만 읽었다). ③ MuJoCo 문서(modeling 장)의 cable 요소 서술. 이 장의 식(요소 변형률 에너지, 기울기·헤시안, 네오-훅 접선, PCC 질량)은 직접 유도했고 유한 차분과 독립 구현으로 확인했다. 서지만 확인한 것(본문은 읽지 못했다): Rucker·Webster(2011), Till·Aloi·Rucker(2019), Boyer 외(2021), Renda 외(2014, 2018), Della Santina 외(2020), Webster·Jones(2010), Simo·Vu-Quoc(1988), Bergou 외(2008), Duriez(2013), Coevoet 외(2017), Goury·Duriez(2018), Newmark(1959), Hilber·Hughes·Taylor(1977), Timoshenko(1921), Cowper(1966). 읽지 못한 것: 교재 Antman(1995), Bonet·Wood(2008). 소프트웨어 버전(PyElastica 1.0.0, SOFA v26.06.00)은 2026-11-05 조회값이며 (확인 필요).

3.10.2 개념 — 연속체의 동역학은 세 층위의 모델로 나뉜다

결론부터 적으면 다음과 같다. 연속체의 운동은 편미분방정식이고, 계산하거나 제어에 사용하려면 공간을 이산화하거나 환원하여 상미분방정식 $M\ddot q+f_{\text{int}}(q)=f_{\text{ext}}$로 바꾼다. 그림 1은 세 층위를 보인다. ① 3차원 연속체를 절점 변위로 이산화하는 FEM은 일반적이지만 자유도가 많다. ② 단면 내부를 환원한 코세라 막대는 1차원 편미분방정식이며 이산화하면 요소당 자유도가 평면에서 3개이다. ③ 형상 자체를 소수의 계수로 가정하는 환원(PCC·변형률 모드)은 자유도가 가장 적지만, 가정한 형상 공간 밖의 변형(구간 안에서 변하는 하중)을 표현하지 못한다. 시간 영역에서는 가장 빠른 모드가 명시 적분의 걸음 크기를 정한다. 가는 막대에서 이 모드는 굽힘이 아니라 단면의 전단·회전 모드이며, 관심 주파수(첫 굽힘 모드)보다 수백 배 높다.

이 장의 계산 결과를 미리 적으면 다음과 같다. ① 코세라 이산화의 굽힘 진동수가 오일러-베르누이 해석해에 2차로 수렴한다(Kirchhoff 극한에서 $N=20,40,80$의 첫 진동수 오차 $-0.212,\,-0.054,\,-0.014\,\%$). 전단과 회전 관성을 포함한 모델($r/L=0.04$)은 티모셴코 외팔보의 연속체 극한(solve_bvp)과 리처드슨 외삽에서 $2\times10^{-6}$ 안에서 같다. ② 같은 막대의 정적 평형이 2.16의 엘라스티카에 2차로 수렴하고 순수 끝 모멘트에서 일정 곡률 호와 일치한다. ③ 속도 베를레의 안정 한계 $\Delta t_c=2/\omega_{\max}$가 3000걸음 시험에서 그대로 나타난다(0.99배는 안정, 1.01배는 20걸음 만에 발산). ④ 뉴마크($\beta=1/4$, $\gamma=1/2$)는 선형 문제에서 수치 진동수가 닫힌 식과 $10^{-11}$ 안에서 같지만, 큰 변형의 비선형 문제에서는 $\Delta t=0.05$(첫 주기의 1/36)에서 에너지가 67 % 증가한다. ⑤ 일정 곡률 모델은 레일리-리츠 상계이며 $n=1$에서 첫 진동수가 $\sqrt{20}=4.472$(+27.2 %), 중력 해제에서 $n=4$의 끝점 낙하는 코세라 $N=40$보다 1.8 % 작다. ⑥ 평면 변형률 FEM에서 3절점 삼각형(T3)은 체적 잠금으로 $\nu=0.499$의 끝 처짐이 이론의 2~5 %에 그친다.

3.10.3 원리

3.10.3.1 코세라 막대의 운동량 평형

평면 운동으로 제한한다. 막대의 상태는 중심선 $p(s,t)\in\mathbb R^2$와 단면 각 $\theta(s,t)$이다. 단면 좌표계에서 본 신장 $\epsilon$, 전단 $\gamma$, 곡률 $u$는 다음이다.

$$\begin{bmatrix}1+\epsilon\\ \gamma\end{bmatrix}=R(\theta)^{\top}p',\qquad u=\theta'.$$

2.16의 비신장·무전단 막대는 $\epsilon=\gamma=0$인 극한이다. 힘 $n$과 모멘트 $m$에 대한 선운동량·각운동량 평형은 다음이다. 여기서 $p'\times n=p'_xn_y-p'_yn_x$이다.

$$\rho A\,\ddot p=n'+f,\qquad \rho I\,\ddot\theta=m'+p'\times n+l.$$

구성 방정식은 $n=R(\theta)\,\mathrm{diag}(EA,\ k_sGA)\,(\epsilon,\gamma)^{\top}$, $m=EI\,u$이다. 켈빈-보이트 감쇠는 $m$에 $c_v\dot u$를 더하여 넣을 수 있으나 이 장은 구현하지 않았다. 고정단에서 $p=0$, $\theta=0$이고 자유단에서 $n=F_{\text{tip}}$, $m=M_{\text{tip}}$이다. 3차원 형식은 단면 좌표계에서 SE(3)의 변형률 장 $\xi$와 속도 트위스트 $\eta$로 적은 $\mathbb M\dot\eta-\mathrm{ad}_{\eta}^{\top}\mathbb M\eta=\Lambda'-\mathrm{ad}_{\xi}^{\top}\Lambda+\bar F$이다(Armanini 외의 식 (13)을 읽은 대로 옮겼으며 기호는 원문을 따른다). 이 장은 평면 식으로 계산한다.

3.10.3.2 선형화와 외팔보 진동수

직선 상태 $p=(s,0)$, $\theta=0$ 둘레의 선형화는 횡변위 $w$와 단면 회전 $\psi$에 대한 티모셴코 보이다.

$$\rho A\,\ddot w=\bigl(k_sGA\,(w'-\psi)\bigr)',\qquad \rho I\,\ddot\psi=(EI\,\psi')'+k_sGA\,(w'-\psi).$$

전단 강성 $k_sGA\to\infty$이고 회전 관성 $\rho I\to0$이면 오일러-베르누이 보 $\rho A\ddot w+EIw''''=0$이 된다. 외팔보의 고유진동수는 $\omega_k=\Omega_k\sqrt{EI/\rho A}/L^2$, $\Omega_k=\beta_k^2$, $\cos\beta_k\cosh\beta_k=-1$이며 $\Omega_1,\Omega_2,\Omega_3=3.5160,\ 22.0345,\ 61.6972$이다. 이 값을 이 장의 기준으로 삼는다. 전단과 회전 관성은 진동수를 낮추고 고차 모드에서 더 크다. 이 장의 모델($r/L=0.04$ 원형 단면, $\nu=0.5$, 전단 계수 $k_s=0.9$)의 연속체 극한은 $\Omega_1=3.5020$, $\Omega_2=21.442$로 오일러-베르누이보다 각각 0.40 %, 2.69 % 낮다. 전단 계수는 원형 단면의 Cowper(1966) 식 $6(1+\nu)/(7+6\nu)$로 알려져 있다(본문 미열람 (확인 필요)).

3.10.3.3 반이산화 — 요소 변형률 에너지

막대를 길이 $\ell=L/N$인 요소 $N$개로 나눈다. 상태 $q$는 절점 위치 $x_i\in\mathbb R^2$(고정단 제외 $2N$)와 요소 중심의 단면 각 $\varphi_e$($N$개)이다. 요소의 모서리 벡터 $d=x_{e+1}-x_e$를 요소 좌표계로 옮긴 성분 $a=d\cdot(\cos\varphi_e,\sin\varphi_e)$, $b=d\cdot(-\sin\varphi_e,\cos\varphi_e)$에서 신장 $a/\ell-1$과 전단 $b/\ell$이 나온다. 퍼텐셜 에너지는 다음이다($\varphi_{-1}=0$은 고정단 단면).

$$\Pi=\sum_{e=0}^{N-1}\frac{\ell}{2}\Bigl[EA\Bigl(\frac a\ell-1\Bigr)^2+k_sGA\Bigl(\frac b\ell\Bigr)^2\Bigr]+\sum_{e=0}^{N-1}\frac{EI}{2\ell}\bigl(\varphi_e-\varphi_{e-1}\bigr)^2c_e-\Pi_{\text{ext}},\quad c_0=2,\ c_{e>0}=1.$$

질량은 절점당 $\rho A\ell$(양 끝은 절반)과 요소당 회전 관성 $\rho I\ell$이다. 운동방정식은 $M\ddot q+\nabla\Pi(q)=0$이고 기울기와 헤시안은 닫힌 식이다. 기울기는 에너지의 중심 차분과 $5.0\times10^{-10}$, 헤시안은 기울기의 중심 차분과 $2.8\times10^{-9}$ 안에서 같다(상대, 6요소 무작위 형상). 전단 항의 $\varphi$ 기울기 $F_ab-F_ba$가 각운동량 평형의 $p'\times n$ 항에 해당한다. Gazzola 외(2018)는 같은 코세라 막대를 유한 차분과 사다리꼴 구적으로 이산화하고 2차 위치 베를레로 적분한다(읽은 범위). 이 장의 구성은 변형률 에너지에서 출발하는 점이 다르며 두 이산화의 동치 여부는 확인하지 않았다 (확인 필요). 이산 탄성 막대(Bergou 외 2008)와 변형률 모드 환원(Boyer 외 2021)은 서지만 확인했다.

3.10.3.4 시간 적분 — 안정 한계와 에너지 거동

명시 적분(속도 베를레)은 2차이고 심플렉틱이며 선형 안정 조건은 $\omega_{\max}\Delta t\le2$이다. 이 막대에서 $\omega_{\max}$를 지배하는 모드는 굽힘이 아니라 요소 회전 자유도이다(최고 모드 운동에너지의 92~98 %가 $\varphi$ 성분). 이 모드의 강성은 전단 $k_sGA\,\ell$과 이웃 굽힘 $4EI/\ell$, 관성은 $\rho I\ell$이므로 $\omega_{\max}\approx\sqrt{(k_sGA+4EI/\ell^2)/\rho I}$로 추정되며 측정과 3.4 % 이내이다(③의 첫 표). 회전 관성 $\rho I\propto r^2$이 작을수록 $\Delta t_c$가 작아진다. Gazzola 외(2018)는 엄밀한 CFL 조건을 유도하지 않고 경험식 $\Delta t=a\,\delta\ell$($a\sim10^{-2}$ s/m)를 사용했다고 적었다.

암시 적분 뉴마크-$\beta$는 다음이다.

$$u_{n+1}=u_n+\Delta t\,v_n+\Delta t^2\bigl[(\tfrac12-\beta)a_n+\beta a_{n+1}\bigr],\qquad v_{n+1}=v_n+\Delta t\bigl[(1-\gamma)a_n+\gamma a_{n+1}\bigr].$$

매 걸음 $R(u_{n+1})=Ma_{n+1}+f(u_{n+1})=0$을 뉴턴법으로 풀며 유효 접선은 $K_T+M/(\beta\Delta t^2)$이다. 선형 문제에서 $\beta=1/4$, $\gamma=1/2$의 수치 진동수 $\bar\omega$는 $\tan(\bar\omega\Delta t/2)=\omega\Delta t/2$를 만족하고(⑤의 뉴마크 표에서 확인), 상대 주기 오차는 $(\omega\Delta t)^2/12$이다. $\gamma>1/2$는 수치 소산을 도입한다. Hilber·Hughes·Taylor(1977)의 HHT-$\alpha$ 같은 소산 방법은 서지만 확인했다.

3.10.3.5 환원 — 일정 곡률 모델은 레일리-리츠 근사이다

시행함수 $\theta(s)=\sum_{i=1}^nq_i\phi_i(s)$를 둔다. $\phi_i$는 구간 $i$에서 0에서 1까지 선형으로 증가하고 이후 1인 램프 함수이다. 그러면 곡률 $\theta'$가 구간마다 상수 $q_i/\ell$이므로 일정 곡률 모델(2.16)이다. 선형화한 변위 $w=\int\theta\,ds$에서 질량과 강성은 다음이다.

$$M_{ij}=\rho A\int_0^L\psi_i\psi_j\,ds,\quad \psi_i=\int_0^s\phi_i\,d\sigma,\qquad K=\operatorname{diag}(EI/\ell).$$

시행함수가 기하 경계 조건을 만족하는 허용 함수이므로 레일리-리츠 정리에 따라 고유값은 정확한 값의 상계이다. $n=1$이면 $w=s^2/(2L)$이고 레일리 몫이 $\Omega_1^2=20$이다. 비선형 운동방정식은 $\rho A\int J_i^{\top}\ddot p\,ds+\partial V/\partial q_i=0$에서 나오며, $J=\partial p/\partial q$로 $\ddot p=J\ddot q-C(s)$, $C(s)=\int_0^s\dot\theta^2\,t\,d\sigma$($t$는 단위 접선)이므로 다음이다.

$$M(q)\,\ddot q=\rho A\int_0^LJ^{\top}C\,ds-\nabla V(q),\qquad V=\tfrac12\tfrac{EI}{\ell}\textstyle\sum q_i^2+V_g.$$

단일 호($n=1$)의 질량은 닫힌 식 $M(a)=\rho A\ell^3\bigl[\tfrac1{3a^2}-\tfrac{2(\sin a-a\cos a)}{a^5}+\tfrac{2(1-\sin a/a)}{a^4}\bigr]$이고 $a\to0$에서 $\rho A\ell^3/20$이다. Armanini 외는 변형률 장 자체를 모드로 환원하는 방식(기하 변수 변형률 접근, Boyer 외 2021 등)을 서술하며, 이 방식은 시행함수가 경계 조건을 만족하지 않아도 된다고 적었다. 원논문은 읽지 못했고 이 장은 구현하지 않았다. PCC의 매개변수화와 제어는 25.3에서 다룬다.

3.10.3.6 3차원 FEM — 초탄성과 접선 강성

변형 구배 $F=I+\nabla u$, $J=\det F$로 두고 평면 변형률 네오-훅 에너지 밀도와 응력은 다음이다(Lamé 상수 $\mu$, $\lambda$).

$$W=\frac\mu2\bigl(\operatorname{tr}(F^{\top}F)-2-2\ln J\bigr)+\frac\lambda2(\ln J)^2,\qquad P=\frac{\partial W}{\partial F}=\mu\bigl(F-F^{-\top}\bigr)+\lambda\ln J\,F^{-\top}.$$

접선은 $A_{iJkL}=\mu\delta_{ik}\delta_{JL}+(\mu-\lambda\ln J)F^{-1}_{Jk}F^{-1}_{Li}+\lambda F^{-1}_{Ji}F^{-1}_{Lk}$이다. 형상함수 $N_a$로 $u\approx\sum_aN_au_a$를 두면 내부력 $f_a=\int P\,\nabla N_a\,dV$, 접선 강성 $K_{ab}=\int A\,\nabla N_a\nabla N_b\,dV$, 질량 $M_{ab}=\int\rho N_aN_b\,dV$이고 운동방정식은 $M\ddot u+D\dot u+f_{\text{int}}(u)=f_{\text{ext}}$이다. 감쇠는 레일리 감쇠 $D=\alpha M+\beta K$가 흔하며 이 장은 구현하지 않았다. 정적 선형 극한 $F\to I$에서 $P\to\mu(\nabla u+\nabla u^{\top})+\lambda\operatorname{tr}(\nabla u)I$이다. 비압축성($\nu\to0.5$, $\lambda\to\infty$)에서 변위 기반 요소가 휨 변형을 과소 계산하는 체적 잠금(volumetric locking)을 보인다고 알려져 있으며(교재 미열람 (확인 필요)), 이 장의 T3와 T6 실험이 정도를 정량화한다. 혼합(변위-압력) 정식화는 구현하지 않았다.

구동은 구속으로 반영한다. Armanini 외(식 (19))는 선형화한 FEM 평형을 $f_{\text{int}}+K\,dq=f_{\text{ext}}+B^{\top}\lambda$로 적고 $B$를 구동기 구속의 자코비안, $\lambda$를 구동기 힘(케이블 장력·공동 압력)으로 둔다. 자유 응답에서 구속 위반 $\delta_{\text{free}}$를 구한 뒤 $\delta=BK^{-1}B^{\top}\lambda+\delta_{\text{free}}$로 $\lambda$를 정한다. 역모델에서는 $\lambda$가 미지수이고 부등 구속이 있으면 QP가 된다(25.4). 일반 이론은 3.11에서 다룬다.

3.10.3.7 구동과 하중의 모델링

텐던 구동은 두 방식이 있다고 Armanini 외는 서술한다. 하나는 텐던 장력이 막대에 주는 분포 힘과 그 모멘트를 외력으로 넣는 방식(Rucker·Webster 2011의 계열)이고, 다른 하나는 가상 일의 원리로 구성 방정식을 능동 형태로 바꾸는 방식이다. 마찰 없는 텐던의 경우 단위 길이당 힘은 장력과 텐던 접선의 변화율의 곱이고 모멘트는 그 힘의 백본에 대한 모멘트이다. 텐던이 백본과 평행하게 일정 간격 $d$로 지나가면 모든 단면에 모멘트 $\tau d$가 걸리므로 곡률이 $\tau d/EI$로 일정하며, 이 장은 그 핵심(모멘트가 일정하면 곡률이 일정하다)을 순수 끝 모멘트 시험으로 확인한다(②). 공압 구동은 공동 압력을 추종 하중으로 넣거나(코세라) 공동 부피 구속의 승수로 넣는다(FEM). 마찰·접촉과 텐던 경로의 상호작용은 이 장 밖이다(3.12).

3.10.4 알고리즘

코세라 반이산 모델(N 요소, 평면):  q = (φ_e, x_{e+1}, y_{e+1}),  M q̈ = −∇Π(q) + Q
명시(속도 베를레):  v ← v + (Δt/2) a;  q ← q + Δt v;  a ← −M⁻¹∇Π(q);  v ← v + (Δt/2) a   (Δt < 2/ω_max)
암시(뉴마크):       u₁ ← u₀ 에서 반복:  R = M a₁ + f(u₁);  (K_T + M/(βΔt²)) Δu = R;  u₁ ← u₁ − Δu
정적 평형:         ∇Π(q) = Q 를 하중 증분 + 직선 탐색 뉴턴법 (K_T = 헤시안, 밴드 촐레스키)
PCC 레일리-리츠:   M_ij = ρA ∫ ψ_i ψ_j ds;  K = diag(EI/ℓ);  (K − ω² M) x = 0
FEM:               F = I + ∇u;  P, A = 네오-훅(F);  f_int = ∫ P∇N;  K_T = ∫ ∇N A ∇N;  M u'' + f_int(u) = f_ext
구속 구동(선형):    [K −Bᵀ; B 0] [u; λ] = [f; δ*]   ⇔   λ = (B K⁻¹ Bᵀ)⁻¹ (δ* − δ_free)

3.10.5 풀이 예제 — 진동수·평형·시간 적분·환원·FEM

결론부터 적으면 다음과 같다. ① 코세라 이산화의 첫 두 굽힘 진동수가 오일러-베르누이 값에 약 2차로 수렴하고, 최고 고유진동수는 전단·회전 가지가 정한다. ② 정적 평형이 엘라스티카에 2차로 수렴한다. ③ 명시 적분의 안정 한계가 $2/\omega_{\max}$에서 정확히 나타나고, 뉴마크의 걸음 크기는 비선형·고주파 모드 때문에 선형 안정성보다 훨씬 작게 제한된다. ④ 일정 곡률 모델은 상계로 수렴하며 $n=4$에서 낙하 깊이가 1.8 % 작다. ⑤ FEM에서 T3는 잠금이 심하고 T6는 $\nu=0.3$에서 1 % 안이다. ⑥ 행렬·에너지·정적 해는 SciPy·scikit-fem과 $10^{-10}$ 이하로 일치하고, 시간 적분은 DOP853과의 차이를 걸음 크기별로 보였으며, MuJoCo cable은 외삽에서 같은 극한에 도달한다.

구현. make_rod·rod_grad·rod_hess·rod_static·verlet·newmark_step·pcc_ritz·neo_hookean과 Fem 클래스(T3·T6 메시, 내부력·접선·질량). 이 장은 밴드 구조를 사용하므로 1.1의 밀집 LU 대신 밴드 촐레스키(bchol·bsolve)와 야코비 고유값(jacobi_eig)을 스크립트에 구현했다. 아래 코드는 핵심 알고리즘이며 PCC 비선형 동역학, FEM 조립, 부분공간 반복, 그림은 스크립트에 있다.

def gl(a, b):
    """구간 [a, b]의 6점 가우스-르장드르 (점, 가중치)"""
    m, h = 0.5 * (a + b), 0.5 * (b - a)
    return [(m + h * x, h * w) for x, w in _GL]

def rod_grad(q, P, load=None, grav=0.0):
    """이산 코세라 막대의 퍼텐셜 기울기.
    q = [phi_e, x_{e+1}, y_{e+1}] (e = 0..N-1)"""
    N, l, EA, S, EI = (P[k] for k in
                       ('N', 'l', 'EA', 'S', 'EI'))
    g = [0.0] * (3 * N)
    for e in range(N):
        x0, y0 = (q[3*e-2], q[3*e-1]) if e else (0.0, 0.0)
        dx, dy = q[3*e+1] - x0, q[3*e+2] - y0
        c, s = cos(q[3*e]), sin(q[3*e])
        a, b = c*dx + s*dy, -s*dx + c*dy
        Fa, Fb = EA*(a/l - 1.0), S*b/l      # 축력, 전단력
        g[3*e] += Fa*b - Fb*a
        g[3*e+1] += Fa*c - Fb*s
        g[3*e+2] += Fa*s + Fb*c
        if e:
            g[3*e-2] -= Fa*c - Fb*s
            g[3*e-1] -= Fa*s + Fb*c
    g[0] += 2*EI*q[0]/l                     # 고정단 굽힘
    for e in range(1, N):
        d = EI*(q[3*e] - q[3*e-3])/l
        g[3*e] += d
        g[3*e-3] -= d
    for i in range(1, N + 1):
        g[3*i-1] += P['m'][i]*grav          # 중력 -y
    if load:                           # 끝 하중 (Px, Py, M)
        g[3*N-2] -= load[0]
        g[3*N-1] -= load[1]
        g[3*N-3] -= load[2]
    return g

def verlet(q, v, a, dt, acc):
    """속도 베를레 1걸음(2차, 심플렉틱). acc(q) -> 가속도"""
    v = [x + 0.5 * dt * y for x, y in zip(v, a)]
    q = [x + dt * y for x, y in zip(q, v)]
    a = acc(q)
    v = [x + 0.5 * dt * y for x, y in zip(v, a)]
    return q, v, a

def newmark_step(fk, M, hb, u, v, a, dt, be=0.25, ga=0.5):
    """뉴마크-베타 1걸음. R(u1) = M a1 + f(u1) = 0 을
    뉴턴법으로 푼다. fk(u) -> (f, K)"""
    n, c = len(u), 1.0 / (be * dt * dt)
    un = u[:]
    for it in range(40):
        an = [c * (un[i] - u[i] - dt * v[i])
              - (0.5 / be - 1.0) * a[i] for i in range(n)]
        f, K = fk(un)
        Ma = bmul(M, an, hb)
        R = [Ma[i] + f[i] for i in range(n)]
        sc = max(max(map(abs, f)), max(map(abs, Ma)))
        if max(map(abs, R)) <= 1e-10 * sc:
            break
        for i in range(n):          # K + M / (be dt^2)
            for j in range(max(0, i - hb),
                           min(n, i + hb + 1)):
                K[i][j] += c * M[i][j]
        try:
            du = bsolve(bchol(K, hb), R, hb)
        except ValueError:          # 비정부호이면 피벗 소거
            du = dense_solve(K, R)
        un = [x - y for x, y in zip(un, du)]
    else:
        raise RuntimeError('뉴턴법 미수렴')
    vn = [v[i] + dt * ((1 - ga) * a[i] + ga * an[i])
          for i in range(n)]
    return un, vn, an

def neo_hookean(F, mu, lam):
    """평면 변형률 네오-훅: 응력 P와 접선 A[i][J][k][L]"""
    J = F[0][0]*F[1][1] - F[0][1]*F[1][0]
    Fi = [[F[1][1]/J, -F[0][1]/J],
          [-F[1][0]/J, F[0][0]/J]]
    lJ = log(J)
    P = [[mu*(F[i][j] - Fi[j][i]) + lam*lJ*Fi[j][i]
          for j in (0, 1)] for i in (0, 1)]
    A = [[[[mu*(i == k)*(j == l)
            + (mu - lam*lJ)*Fi[j][k]*Fi[l][i]
            + lam*Fi[j][i]*Fi[l][k] for l in (0, 1)]
           for k in (0, 1)] for j in (0, 1)] for i in (0, 1)]
    return P, A

def pcc_ritz(n, EI, rA, L):
    """균등 n구간 일정 곡률(굽힘각 q_i)의
    선형 레일리-리츠 행렬 M, K"""
    l = L / n
    K = [[EI / l * (i == j) for j in range(n)]
         for i in range(n)]
    M = [[0.0] * n for _ in range(n)]
    for k in range(n):              # 구간 k: 4차 다항식 구적
        for s, w in gl(k * l, (k + 1) * l):
            ps = []                 # psi_i(s) = phi_i의 적분
            for i in range(n):
                t = s - i * l
                ps.append(0.0 if t <= 0 else t * t / (2 * l)
                          if t <= l else l / 2 + t - l)
            for i in range(n):
                for j in range(n):
                    M[i][j] += rA * w * ps[i] * ps[j]
    return M, K

① 코세라 이산화의 선형 진동수.

Kirchhoff 극한(왼쪽 네 열)에서 두 진동수 모두 $N$을 두 배로 늘릴 때 오차가 약 1/4로 줄어 2차 수렴이며, 첫 진동수의 수렴 차수는 $N=20\to40$에서 1.98, $40\to80$에서 1.92이다. 전단 포함 모델은 오일러-베르누이보다 첫 진동수가 0.4 %, 둘째가 2.7 % 낮은 값으로 수렴한다. 이 차이는 이산화 오차가 아니라 티모셴코 보의 전단·회전 관성 효과이며, solve_bvp로 푼 티모셴코 외팔보의 연속체 극한 $\Omega_1=3.5020$, $\Omega_2=21.442$와 리처드슨 외삽값이 $2\times10^{-6}$ 안에서 일치한다(⑦).

② 정적 평형과 끝 모멘트.

끝 하중 평형 해의 끝점 오차(2.16 엘라스티카 대비, 단위 $L$)가 $N$을 두 배로 늘릴 때 3.6~4.0배 줄어 2차 수렴한다($N=80$, $\alpha=5$에서 $6.3\times10^{-5}$). 전단 포함 모델($r/L=0.04$)은 전단 변형 때문에 끝점이 $\alpha=2$에서 $1.9\times10^{-3}L$ 더 처진다. 순수 끝 모멘트 $M=EI\kappa$에서는 일정 곡률 호 $(\sin a/\kappa,\ (1-\cos a)/\kappa)$와 끝점이 $\kappa L=\pi$에서 $N=20,40,80$일 때 $6.5\times10^{-4},\,1.6\times10^{-4},\,4.1\times10^{-5}$ 차이(2차 수렴)이다. 이 장의 동적 모델이 정적 극한에서 2.16의 PCC·엘라스티카와 같은 결과를 낸다는 확인이다. 점성 감쇠 $c=2\omega_1$의 동적 이완(속도 베를레, $0.9\,\Delta t_c$)은 $\alpha=2$, $N=20$에서 9000걸음($t=6.46$)에 속도가 $10^{-9}$ 이하가 되었고 뉴턴법 평형 해와 $3.1\times10^{-9}$ 안에서 같다.

③ 시간 적분.

최고 고유진동수의 전단·회전 가지 추정이 측정과 3.4 %($N=20$, 최대)에서 0.6 %($N=80$) 안으로 맞고 축방향 추정 $2c/\ell$($c=\sqrt{EA/\rho A}$)보다 크다. 명시 한계 $\Delta t_c$에서 첫 굽힘 주기는 $N=20$에서 2256걸음, $N=80$에서 7323걸음이다. 직선 상태의 선형 시험에서 $0.9\,\Delta t_c$와 $0.99\,\Delta t_c$는 3000걸음 동안 에너지 $E/E_0$가 1.15 이하로 유지되었고, $1.01\,\Delta t_c$는 20걸음, $1.1\,\Delta t_c$는 7걸음 만에 에너지가 100배를 넘었다. 안정 한계가 $2/\omega_{\max}$에서 정확히 나타난다.

베를레 $0.9\,\Delta t_c$는 에너지 편차 $1.8\times10^{-5}$로 기준(0.1배)과 영점 통과 시각이 $1.3\times10^{-5}$(상대) 안에서 같다. 뉴마크는 $\Delta t=0.01,0.02$에서 에너지 편차 $5.4\times10^{-4},\,3.6\times10^{-3}$이고 시각 오차 $4.4\times10^{-5}$, $1.2\times10^{-3}$(상대)이다. $\Delta t=0.05$(첫 주기의 1/36)에서는 에너지가 67 % 증가하고, $0.1$에서는 뉴턴법이 수렴하지 않는다. 선형 문제의 무조건 안정성은 이 비선형 문제에서 성립하지 않았다. 원인(고주파 모드의 위상 오차, 큰 회전의 비선형)은 분리하지 않았다 (확인 필요). $\gamma=0.6$, $\beta=0.3025$는 안정하게 풀리지만 에너지가 34 %($\Delta t=0.05$)~54 %($0.1$) 소산하고 시각이 1.6~3.9 % 어긋난다. 소산 없이 큰 걸음을 얻는 방법은 이 장에서 시험하지 않았다.

④ 일정 곡률 환원과 중력 해제.

레일리-리츠 진동수는 모든 모드에서 오일러-베르누이 값보다 크다(상계). 첫 진동수 오차는 $n=1,2,4,6$에서 $+27.2,\,+5.1,\,+1.2,\,+0.54\,\%$이고 둘째·셋째 모드는 더 느리게 수렴한다. 최고 진동수는 $n=4$에서 204로 코세라 $N=20$의 $\omega_{\max}$(2509)의 약 1/12이어서 환원은 강성 모드도 함께 제거한다. $n=1$의 질량 $M(a)$ 닫힌 식은 구적과 $a=0.5,1,2,3$에서 $9\times10^{-12}$ 안에서 같고 $a=10^{-3}$에서 $M\to0.05=1/20$이다.

수평 직선에서 정지 출발한 외팔보(중력 $\rho AgL^3/EI=3$)의 끝점 낙하 깊이는 코세라가 $N=10,20,40$에서 $-0.6227,\,-0.6193,\,-0.6182$로 수렴하고, PCC는 $n=1,2,3,4$에서 $-0.4424,\,-0.5721,\,-0.5963,\,-0.6068$로 아래로(깊이가 커지는 쪽으로) 수렴한다. 에너지 표류는 RK4로 $8\times10^{-7}$ 이하이다. 그림 2는 첫 진동수 비를, 그림 3은 끝점 높이의 시간 이력을 보인다.

⑤ 평면 변형률 FEM.

T3는 $\nu=0.3$에서도 휨을 크게 과소 계산하여 20×2 메시의 끝 처짐이 이론의 52 %이고(상수 변형률 요소의 알려진 한계 (확인 필요)) $\nu=0.499$에서 2.3 %(20×2)~5.2 %(40×4)에 그친다. T6는 $\nu=0.3$에서 1 % 안이고 $\nu=0.499$에서 10×1은 85 %, 40×4는 98 %이다. 체적 잠금은 변위 기반 요소를 비압축에 가깝게 사용할 때 메시를 매우 세분해야 완화되며 소프트 로봇의 재료(비압축에 가까움)에서 중요하다. T6 40×4의 1.0041은 미세 메시 기준(skfem P2 80×8)의 1.0048과 0.07 % 안에서 같고, 1에서의 0.5 % 차이에는 전단 변형이 섞여 있다(원인은 분리하지 않았다).

T3의 첫 진동수는 20×2에서 +37.8 %, 40×4에서 +10.5 %로 높고, T6는 10×1에서 +0.33 %, 20×2에서 −0.38 %이다. 둘째 모드의 −3.5 %, −4.9 %는 전단·회전 관성의 영향을 포함한다. 같은 T6 20×2의 값은 미세 메시 기준(skfem P2 80×8)과 $2.2\times10^{-3}$($\Omega_1$)·$2.8\times10^{-3}$($\Omega_2$) 안에서 일치한다(⑦).

선형 문제에서 뉴마크($\beta=1/4$, $\gamma=1/2$) 수치 진동수가 닫힌 식 $\tan(\bar\omega\Delta t/2)=\omega\Delta t/2$와 $3\times10^{-12}$ 안에서 같다. 주기당 20걸음($T_1/20$)에서 진동수가 0.8 % 낮다.

접선 강성은 내부력의 유한 차분과 $6\times10^{-9}$ 이하로 일치하고 내부력은 에너지의 기울기와 $2\times10^{-10}$ 안에서 같다. 큰 변형($\alpha=1$)에서 주기당 20걸음 뉴마크는 에너지 편차가 최대 4.0 %이며 종료 시 +3.6 %이다. 명시 중앙차분의 안정 한계는 선형 문제에서 $0.95\,\Delta t_c$는 유계(변위 비 1.97), $1.05\,\Delta t_c$는 발산($7\times10^{81}$)이다. 구속 구동의 KKT 해와 슈어 보수 해는 $2\times10^{-15}$ 안에서 같다.

⑥ 물리 척도(가정).

실리콘 계열 팔을 가정한 값이며 재료 물성은 이 장이 확인한 데이터가 아니다 (확인 필요). 첫 굽힘 진동수는 약 1 Hz인데 명시 적분 한계는 0.47 ms여서 한 주기에 2256걸음이 필요하다.

⑦ 라이브러리 대조(실행 2026-11-05). SciPy 1.18.1의 eigh(일반화 고유값)·solve_bvp(고유값을 미지 매개변수로)·solve_ivp(DOP853), scikit-fem 12.0.2(P1·P2 평면 변형률, 네오-훅 에너지·잔차·접선 형식, 미세 메시), MuJoCo 3.17.0의 cable 플러그인(mjd_transitionFD)을 같은 입력에서 비교했다. 스크립트는 robotics/fig-src/3.10-libcheck.py, 출력은 robotics/fig-src/data/3.10-libcheck-output.txt에 있다.

코세라 선형 고유값은 eigh와 $10^{-10}$ 안에서, PCC 레일리-리츠는 $10^{-13}$ 안에서 같다. FEM의 변형 에너지·내부력·접선 형식은 scikit-fem과 T3·T6 모두 $3\times10^{-14}$ 안에서 같고 큰 변형 정적 해는 $7\times10^{-15}$(상대), 선형 고유값은 $8\times10^{-11}$ 안에서 같다. solve_bvp의 티모셴코 연속체 극한은 상태 코드 0(수렴)이며 이 장의 $N=80$과 $10^{-4}$~$10^{-3}$, 리처드슨 외삽과 $10^{-7}$~$10^{-6}$ 안에서 같다. 베를레는 DOP853 대비 끝점 차이가 걸음을 반으로 줄일 때마다 $7.7,\,4.1,\,1.4\times10^{-5}$로 감소하고(비 1.9, 2.9; 2차의 4에 접근 중) 영점 통과 시각은 $0.1\,\Delta t_c$에서 $3.6\times10^{-7}$(상대)까지 일치한다.

이상 징후. 뉴마크는 DOP853 대비 끝점 높이 차이가 $\Delta t=0.01$에서도 $2.5\times10^{-3}$이고 $0.02$에서 $6.0\times10^{-3}$인데, 영점 통과 시각의 상대 차이는 $4.4\times10^{-5}$, $1.2\times10^{-3}$이다. 점별 차이가 시각 차이보다 훨씬 크므로, 저주파 위상 오차보다 뉴마크가 실제와 다른 주파수로 진동시키는 고주파 모드 성분이 지배하는 것으로 추정한다. 필터링으로 분리하지 않았다 (확인 필요). 네오-훅 FEM의 뉴마크도 $T_1/20,\,T_1/40,\,T_1/80$에서 $1.5\times10^{-2},\,5.8\times10^{-3},\,2.9\times10^{-3}$으로 줄어, 걸음을 반으로 줄일 때의 감소가 2차의 4배보다 작다. MuJoCo cable의 첫 굽힘 진동수는 구간 수 10, 20, 40, 80에서 $3.905,\,3.701,\,3.606,\,3.560$으로 오일러-베르누이 값 3.516보다 높고, 오차가 구간 수를 두 배로 할 때 약 1/2로 줄어 1차 수렴한다. 이 장의 2차 수렴과 다르며, cable의 첫 요소가 관절 없이 고정되어 고정단에 강체 구간이 생기기 때문으로 추정하나 확인하지 않았다 (확인 필요). 외삽값은 오일러-베르누이와 $3.2\times10^{-4}$($\Omega_1$), $5.3\times10^{-4}$($\Omega_2$) 안에서 같다. MuJoCo 대조는 총 질량으로 정규화한 값이며, 캡슐 끝 반구가 겹쳐 구간 수에 따라 총 질량이 달라진다.

한계. 평면 막대 하나와 평면 변형률 외팔보 하나의 결과이며 3차원 비틀림, 켈빈-보이트·레일리 감쇠, 텐던·공압 구동의 구현, 접촉, 다구간 결합은 구현하지 않았다. 코세라 이산화는 신장·전단을 둔 형식이고 변형률 에너지에서 출발하는 구성이라 PyElastica 등의 구현과 수치적으로 동일하다고 확인하지 않았다 (확인 필요). 시간 적분 비교는 $N=20$, $r/L=0.04$ 한 막대와 한 초기 조건이며 뉴마크가 불안정해지는 걸음 크기의 일반 기준은 구하지 않았다. 혼합 정식화, 비압축 재료의 안정 네오-훅(Smith 외 2018)과 모델 차수 축소(Goury·Duriez 2018)는 구현하지 않았고 서지만 확인했다. 물리 척도 표의 물성은 가정이다.

견해: 연속체·소프트 로봇의 동역학 모델은 용도에 따라 층위를 나누어 고르는 것을 기본으로 둔다고 본다. 실시간 제어와 역기구학에는 환원 모델(PCC 구간 3~4개, 또는 변형률 모드)을, 하중 하 형상과 힘의 검증에는 코세라 이산화를, 복잡한 형상·접촉·초탄성 재료에는 FEM을 사용하고 환원 모델은 FEM이나 코세라 모델과의 대조로 검증한다. 근거는 이 장의 계산에서 ① 같은 문제에서 PCC $n=4$(자유도 4)가 코세라 $N=40$(자유도 120)과 낙하 깊이 1.8 %, $n=1$은 28 % 어긋났고, ② 환원이 최고 진동수를 코세라 $N=20$의 1/12로 낮춰 강성 모드를 제거하며, ③ T3는 $\nu=0.499$에서 20×2~40×4 메시의 처짐이 이론의 2~5 %에 그쳐 요소 선택이 결과를 좌우했다는 점이다. 허용 오차와 계산 예산은 응용마다 다르며 이 장은 기준을 정하지 않았다.

3.10.6 설계 트레이드오프 — 모델 층위의 선택

연속체 로봇의 동적 모델은 정확도, 자유도, 시간 적분의 강성, 구현 비용 사이에서 다음 중에서 고른다.

• ① 코세라 이산화(반이산 막대). 요소당 자유도가 적고 신장·전단·굽힘을 모두 가진다. 하중 하 형상이 정확하나 가는 막대에서 명시 적분 한계가 전단·회전 모드로 정해져 걸음이 작다.
• ② 환원 모델(PCC·변형률 모드). 자유도가 가장 적고 강성 모드가 없어 큰 걸음을 사용할 수 있다. 가정한 형상 밖의 변형을 표현하지 못하며 구간 수·모드 수에 따라 레일리-리츠 상계로 수렴한다.
• ③ 3차원 FEM. 복잡한 단면·재료·접촉을 다루나 자유도가 많고 요소 선택(잠금)과 강성 때문에 암시 적분과 큰 연산 자원이 필요하다.

모델 선택은 층위만이 아니라 문제의 시간 척도로도 갈린다. 관심 주파수가 첫 굽힘 모드 근처이면 환원 모델로 충분하고, 충격·접촉·고주파 응답이 관심이면 코세라나 FEM이 필요하다. 이 장은 접촉 응답을 계산하지 않았다.

3.10.7 응용 축 비교

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

3.10.8 관련 규격·프로토콜

연속체·소프트 로봇의 동역학 모델을 직접 규정하는 국제 표준은 확인하지 못했다 (확인 필요). 표준 색인은 A.3.

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

연속체 로봇 동역학 자체에 적용되는 규제는 없다. 의료용 연속체 로봇(카테터·내시경형, 32.6)이 시뮬레이션 모델을 설계 근거나 시험 대체로 사용할 때에는 모델의 검증·확인·신뢰성 평가가 요구되며, 방법은 41.14과 26.13에서 다룬다. 이 장의 대조 실험(독립 구현과의 일치, 연속체 극한과의 수렴, 보존량 확인)은 모델 검증의 예이지 규제 적합성 입증이 아니다. 법적 판단은 이 책의 범위 밖이다.

3.10.10 한계·유의점·미해결

• 평면 모델이다. 3차원 비틀림·전단의 결합과 텐던 경로의 3차원 효과는 계산하지 않았다.
• 재료는 선형 탄성(막대)과 네오-훅(FEM)이며 점탄성·히스테리시스·소프트 재료의 변형률 속도 의존은 모델에 없다.
• 시간 적분의 일반 안정 기준(비선형 뉴마크)은 구하지 않았다.
• FEM은 변위 기반 요소만 구현했다. 비압축 재료의 혼합 요소와 안정 네오-훅은 서지만 확인했다.
• 구동(텐던·공압)과 접촉의 구현, 다구간 결합, 3차원 FEM(사면체)과 모델 차수 축소는 이 장 밖이다.
• 라이브러리 대조는 선형·정적 극한과 한 비선형 사례에서의 일치이며 실험 데이터와의 비교(물리적 타당성)가 아니다.

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

체크포인트.

1. 코세라 막대의 신장·전단·곡률과 힘·모멘트를 연결하는 구성 방정식을 적고, $k_sGA\to\infty$, $\rho I\to0$ 극한에서 어떤 보가 되는가?
2. 속도 베를레의 안정 조건은 무엇이며, 이 막대에서 $\omega_{\max}$를 정하는 모드는 무엇인가?
3. 일정 곡률 모델의 고유진동수가 상계인 이유는 무엇이며 $n=1$의 첫 진동수 비는 얼마인가?
4. 변위 기반 T3 요소가 $\nu\to0.5$에서 휨을 과소 계산하는 이유는 무엇인가?
5. 구동기를 구속으로 넣을 때 $\lambda$는 어떻게 구하는가?

연구 질문.

• 3차원 코세라 막대(비틀림 포함)에서 전단·회전 가지가 정하는 명시 적분 한계를 해석적으로 구할 수 있는가?
• 비선형 문제에서 뉴마크가 불안정해지는 걸음 크기를 모드 구성으로 예측할 수 있는가?
• 변형률 모드 환원의 모드 수는 응용(제어 대 형상 검증)별로 어떻게 정하는가?
• 혼합 FEM과 안정 네오-훅의 비압축 한계에서 계산 비용은 T3·T6 대비 어떻게 달라지는가?

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

• Armanini, C., Boyer, F., Mathew, A. T., Duriez, C., Renda, F., "Soft Robots Modeling: A Structured Overview," IEEE Transactions on Robotics 39(3):1728–1748, 2023. DOI 10.1109/TRO.2022.3231360. arXiv:2112.03645 v2를 읽었다.
• Gazzola, M., Dudte, L. H., McCormick, A. G., Mahadevan, L., "Forward and inverse problems in the mechanics of soft filaments," Royal Society Open Science 5(6):171628, 2018. DOI 10.1098/rsos.171628.
• Rucker, D. C., Webster III, R. J., "Statics and Dynamics of Continuum Robots With General Tendon Routing and External Loading," IEEE Transactions on Robotics 27(6):1033–1044, 2011. DOI 10.1109/TRO.2011.2160469 (서지만 확인).
• Till, J., Aloi, V., Rucker, C., "Real-time dynamics of soft and continuum robots based on Cosserat rod models," International Journal of Robotics Research 38(6):723–746, 2019. DOI 10.1177/0278364919842269 (서지만 확인).
• Boyer, F., Lebastard, V., Candelier, F., Renda, F., "Dynamics of Continuum and Soft Robots: A Strain Parameterization Based Approach," IEEE Transactions on Robotics 37(3):847–863, 2021. DOI 10.1109/TRO.2020.3036618 (서지만 확인).
• Renda, F., Giorelli, M., Calisti, M., Cianchetti, M., Laschi, C., "Dynamic Model of a Multibending Soft Robot Arm Driven by Cables," IEEE Transactions on Robotics 30(5):1109–1122, 2014. DOI 10.1109/TRO.2014.2325992 (서지만 확인).
• Renda, F., Boyer, F., Dias, J., Seneviratne, L., "Discrete Cosserat Approach for Multisection Soft Manipulator Dynamics," IEEE Transactions on Robotics 34(6):1518–1533, 2018. DOI 10.1109/TRO.2018.2868815 (서지만 확인).
• Della Santina, C., Bicchi, A., Rus, D., "On an Improved State Parametrization for Soft Robots With Piecewise Constant Curvature and Its Use in Model Based Control," IEEE Robotics and Automation Letters 5(2):1001–1008, 2020. DOI 10.1109/LRA.2020.2967269 (서지만 확인).
• Webster III, R. J., Jones, B. A., "Design and Kinematic Modeling of Constant Curvature Continuum Robots: A Review," International Journal of Robotics Research 29(13):1661–1683, 2010. DOI 10.1177/0278364910368147 (서지만 확인).
• Simo, J. C., Vu-Quoc, L., "On the dynamics in space of rods undergoing large motions — A geometrically exact approach," Computer Methods in Applied Mechanics and Engineering 66(2):125–161, 1988. DOI 10.1016/0045-7825(88)90073-4 (서지만 확인).
• Bergou, M., Wardetzky, M., Robinson, S., Audoly, B., Grinspun, E., "Discrete elastic rods," ACM Transactions on Graphics 27(3), 2008. DOI 10.1145/1360612.1360662 (서지만 확인).
• Duriez, C., "Control of elastic soft robots based on real-time finite element method," IEEE ICRA 2013:3982–3987. DOI 10.1109/ICRA.2013.6631138 (서지만 확인).
• Coevoet, E. 외, "Software toolkit for modeling, simulation, and control of soft robots," Advanced Robotics 31(22):1208–1224, 2017. DOI 10.1080/01691864.2017.1395362 (서지만 확인).
• Goury, O., Duriez, C., "Fast, Generic, and Reliable Control and Simulation of Soft Robots Using Model Order Reduction," IEEE Transactions on Robotics 34(6):1565–1576, 2018. DOI 10.1109/TRO.2018.2861900 (서지만 확인).
• Newmark, N. M., "A Method of Computation for Structural Dynamics," Journal of the Engineering Mechanics Division (ASCE) 85(3):67–94, 1959. DOI 10.1061/JMCEA3.0000098 (서지만 확인).
• Hilber, H. M., Hughes, T. J. R., Taylor, R. L., "Improved numerical dissipation for time integration algorithms in structural dynamics," Earthquake Engineering & Structural Dynamics 5(3):283–292, 1977. DOI 10.1002/eqe.4290050306 (서지만 확인).
• Smith, B., De Goes, F., Kim, T., "Stable Neo-Hookean Flesh Simulation," ACM Transactions on Graphics 37(2), 2018. DOI 10.1145/3180491 (서지만 확인).
• Timoshenko, S. P., "On the correction for shear of the differential equation for transverse vibrations of prismatic bars," Philosophical Magazine 41(245):744–746, 1921. DOI 10.1080/14786442108636264 (서지만 확인).
• Cowper, G. R., "The Shear Coefficient in Timoshenko's Beam Theory," Journal of Applied Mechanics 33(2):335–340, 1966. DOI 10.1115/1.3625046 (서지만 확인).
• MuJoCo 문서, modeling 장, cable 요소(https://mujoco.readthedocs.io/en/stable/modeling.html, 조회 2026-11-05).
• Antman, S. S., Nonlinear Problems of Elasticity, Springer, 1995; Bonet, J., Wood, R. D., Nonlinear Continuum Mechanics for Finite Element Analysis, Cambridge University Press, 2008 (읽지 못했다).

다음 장. 3.11은 구속 동역학과 폐연쇄(병렬 기구 동역학·라그랑주 승수)를 다룬다. 이 장의 구속 기반 구동(식 (19)의 $B^{\top}\lambda$)과 슈어 보수가 3.11의 라그랑주 승수와 KKT 계의 특수한 경우이다.