A Dual-Quaternion Framework for Bennett-Limit Diagnostics in Rigid Kresling Origami FOLD–TWIST–FOLD Robots

1. Introduction

Kresling-type origami structures constitute a class of fold–twist mechanisms characterized by coupled axial and rotational motion. In the rigid-origami regime, their panels are treated as rigid bodies and their crease lines as revolute joints, so that the resulting motion can be interpreted within the framework of spatial kinematics. This makes Kresling structures relevant not only as deployable origami systems, but also as candidates for rigid origami robots whose motion must be described consistently in SE(3)[43,44].

A classical reference point for spatial closed-chain compatibility is the Bennett mechanism, whose mobility is governed by geometric relations between skew revolute axes [22,23,24]. In the present work, the Bennett condition is not used to claim that every Kresling unit is a non-degenerate Bennett 4R linkage. Instead, it is used as a local compatibility diagnostic. Since Kresling crease geometries naturally include intersecting-axis configurations, the relevant local regime is the intersecting-axis limit d=0, where the Bennett ratio condition is evaluated in its limiting form.

Dual quaternions provide a compact algebraic representation of rigid-body motions in the Euclidean group SE(3), and are widely used in spatial kinematics, robotics, and mechanism theory [11,12,13,14,15,16,17,18]. They are closely related to the classical Study quadric formulation of rigid displacements [21]. In mechanism synthesis and analysis, dual quaternions have been used for motion representation, factorization, and construction of spatial linkages, including the works of McCarthy [27], Pérez and McCarthy [1], Rad and Schröcker [2], Hegedüs et al. [3,4], Kong [5], and related contributions [6,7,8,9,10]. However, classical compatibility conditions such as Bennett-type constraints are often invoked geometrically rather than interpreted directly through the dual-quaternion closure relations of the associated spatial chains.

This paper develops a rigid dual-quaternion formulation for Kresling-type fold–twist–fold motion. The formulation is restricted to the rigid regime: panels are assumed rigid, crease axes are fixed revolute axes, and deformable folds or evolving crease axes are not considered. Within this setting, local spatial-chain closure is interpreted through unit dual quaternions, and Bennett-limit compatibility is evaluated as a diagnostic in the intersecting-axis case d=0. This provides a consistent SE(3)-based framework for analyzing the local compatibility of Kresling crease geometries under the stated rigid-axis assumptions.

In addition to the theoretical formulation, the paper presents a numerical validation of a rigid Kresling origami FOLD–TWIST–FOLD trajectory. The final validation uses a DQ-parametric homotopy that organizes the motion into the prescribed sequence: fold, twist, and fold. The tested configuration consists of a three-cell hexagonal Kresling origami robot with $N_s=6$, $\nu =3$, $R=48\hspace{0.17em}\mathrm{m}\mathrm{m}$, $H_0=60\hspace{0.17em}\mathrm{m}\mathrm{m}$, and initial twist ${30}^{\circ }$. The prescribed command $\mathsf{\Delta }{h}_1=-2\hspace{0.17em}\mathrm{m}\mathrm{m}$, $\mathsf{\Delta }{\varphi }_2=5^{\circ }$, and $\mathsf{\Delta }{h}_3=-2\hspace{0.17em}\mathrm{m}\mathrm{m}$is achieved to numerical precision. Ring-edge, primary-crease, diagonal-crease, and triangular-panel rigidity are preserved over the full trajectory, while unit-DQ consistency, DQ panel-map consistency, and the Bennett-limit diagnostic are satisfied throughout the trajectory.

The main contributions of this work are therefore threefold. First, it formulates rigid Kresling fold–twist–fold kinematics using unit dual quaternions in SE(3). Second, it clarifies how Bennett-type compatibility can be evaluated locally in the intersecting-axis limit relevant to Kresling crease geometries. Third, it validates a DQ-parametric homotopy trajectory for a rigid multi-cell Kresling robot, showing machine-precision preservation of rigidity and dual-quaternion consistency over the full fold–twist–fold motion.

The remainder of the paper is organized as follows. Section 2 introduces the dual-quaternion framework, the rigid kinematic model, and the Bennett-limit compatibility interpretation. Section 3 presents the DQ-parametric homotopy validation of the fold–twist–fold Kresling trajectory. Section 4 discusses limitations, including the scope of the rigid-axis assumption and the role of paired-normal diagnostics. Section 5 concludes the paper and outlines future extensions.

2. Materials and Methods

2.1. Scope, conventions, and hypotheses

Notation. Bold lowercase $\mathit{a},\mathit{b},\dots$ are vectors in ${\mathbb{R}}^3$; $\parallel \cdot \parallel$ is the Euclidean norm; $\times$ and $\cdot$ are cross and dot products. A unit quaternion is $\mathrm{r}={\mathrm{r}}_0+\mathit{r}$ with ${\mathrm{r}}_0\in \mathbb{R}$, $\mathit{r}\in {\mathbb{R}}^3$ (identified with pure imaginary quaternions). Conjugation ${\mathrm{r}}^{\mathrm{*}}={\mathrm{r}}_0-\mathit{r}$, and unit means $\mathrm{r}{\mathrm{r}}^{\mathrm{*}}=1 \left[\mathrm{12,13}\right]$. A (unit) dual quaternion (DQ) is $\mathrm{q}=\mathrm{r}+\mathsf{\epsilon }\mathrm{d}$, with ${\mathsf{\epsilon }}^2=0$. [11,14,17] Conjugation: ${\mathrm{q}}^{\mathrm{*}}={\mathrm{r}}^{\mathrm{*}}+\mathsf{\epsilon }{\mathrm{d}}^{\mathrm{*}}$. The unit DQ constraints are [17,18]

$$\mathrm{q}{\mathrm{q}}^{\mathrm{*}}=1\iff \left\{\begin{array}{c}r{\mathrm{r}}^{\mathrm{*}}=1,\\ r{\mathrm{d}}^{\mathrm{*}}+d{\mathrm{r}}^{\mathrm{*}}=0\hspace{0.25em}\end{array}\right.$$

(1)

A rigid motion $(\mathrm{R},\mathit{t})$ is encoded by the unit DQ [17]

$$\mathrm{q}=\mathrm{r}+\mathsf{\epsilon }\hspace{0.17em}1/2\hspace{0.17em}\mathit{t}\hspace{0.17em}\mathrm{r}, \mathrm{r}\leftrightarrow \mathrm{R}\in \mathrm{S}\mathrm{O}\left(3\right).$$

(2)

$${\mathit{t}=0+t}_x\mathrm{i}+t_y\mathrm{j}+t_z\mathrm{k}.$$

Here a $3$-vector $\mathit{v}$ is identified with the pure quaternion $\mathit{v}$, and quaternion multiplication implements the adjoint action: $\mathrm{r}\hspace{0.17em}\mathit{v}\hspace{0.17em}{\mathrm{r}}^{\mathrm{*}}$ corresponds to $\mathrm{R}\mathit{v}$. [12,17]

Axis notation.

Each revolute joint axis is a directed line $\mathrm{L}=(\mathit{p},\mathit{u})$ with $\mathit{p}\in {\mathbb{R}}^3$ and unit direction $\mathit{u}$ (so the line is $\mathit{p}+\mathbb{R}\hspace{0.17em}\mathit{u}$). The point $\mathit{p}$ is not unique; it may be shifted along the axis without changing the line. For two successive axes ${\mathrm{L}}_{\mathrm{i}}=({\mathit{p}}_{\mathrm{i}},{\mathit{u}}_{\mathrm{i}})$ and ${\mathrm{L}}_{\mathrm{i}+1}=({\mathit{p}}_{\mathrm{i}+1},{\mathit{u}}_{\mathrm{i}+1})$:

the (acute) skew angle ${\mathsf{\alpha }}_{\mathrm{i},\mathrm{i}+1}\in (0,{\displaystyle \frac{\mathsf{\pi }}{2}}]$ satisfies $\mathrm{s}\mathrm{i}\mathrm{n}{\mathsf{\alpha }}_{\mathrm{i},\mathrm{i}+1}=\parallel {\mathit{u}}_{\mathrm{i}}\times {\mathit{u}}_{\mathrm{i}+1}\parallel$;

the minimal distance ${\mathrm{d}}_{\mathrm{i},\mathrm{i}+1}>0$ satisfies [24,32]

$$\parallel {\mathit{u}}_{\mathrm{i}}\times {\mathit{u}}_{\mathrm{i}+1}\parallel \hspace{0.17em}{\mathrm{d}}_{\mathrm{i},\mathrm{i}+1}=\left|\right({\mathit{p}}_{\mathrm{i}+1}-{\mathit{p}}_{\mathrm{i}})\cdot ({\mathit{u}}_{\mathrm{i}}\times {\mathit{u}}_{\mathrm{i}+1}\left)\right|\hspace{0.25em}.$$

(3)

Let the unit normal between them be

$${\nu }_{\mathrm{i},\mathrm{i}+1}:={\displaystyle \frac{{\mathit{u}}_{\mathrm{i}}\times {\mathit{u}}_{\mathrm{i}+1}}{\parallel {\mathit{u}}_{\mathrm{i}}\times {\mathit{u}}_{\mathrm{i}+1}\parallel }}\hspace{0.25em}.{\nu }_{\mathrm{i},\mathrm{i}+1}$$

(4)

is defined up to sign depending on the chosen orientations of ${\mathit{u}}_{\mathrm{i}},{\mathit{u}}_{\mathrm{i}+1}$Standing hypotheses (non-degenerate 4R).

(H1) Adjacent axes are not parallel: ${\mathit{u}}_i\nparallel {\mathit{u}}_{i+1}$, equivalently $\parallel {\mathit{u}}_i\times {\mathit{u}}_{i+1}\parallel >0$.

(H2) Distinct axes and positive separations: $d_{i,i+1}>0$.

(H3) Nonzero finite skew angles: 0< ${\alpha }_{i,i+1}$​<π/2, equivalently $\mathrm{s}\mathrm{i}\mathrm{n}{\alpha }_{i,i+1}>0$.

(H4) Zero pitch joints: all joints are pure revolute joints; nonzero-pitch screw joints are discussed separately in Section 2.7.1.

Edge cases (limits).If$\alpha \to 0$(nearly parallel axes) or$d\to 0$(nearly intersecting axes), these cases are outside the standing hypotheses$\left(H1\right)$–$\left(H4\right)$and are considered only as limiting regimes. The formulas below have continuous limits only under the corresponding nondegenerate limiting assumptions, which are stated where needed.

In the nearly parallel limit $\alpha \to 0$, since $\mathrm{s}\mathrm{i}\mathrm{n}\alpha \sim \alpha$, the Bennett ratio

$${\displaystyle \frac{d}{\sin \alpha }}$$

remains finite only if the minimum distance $d$tends to zero at least proportionally to $\sin \alpha$. Equivalently, one assumes a bounded-ratio scaling

$$d=O\left(\mathrm{s}\mathrm{i}\mathrm{n}\alpha \right)$$

If $d$remains bounded away from zero while $\alpha \to 0$, then ${\displaystyle \frac{d}{\mathrm{s}\mathrm{i}\mathrm{n}\alpha }}$diverges.

In the nearly intersecting limit $d\to 0$with $\alpha$bounded away from zero, the ratio

$${\displaystyle \frac{d}{\sin \alpha }}$$

tends to zero. If $d\to 0$and $\alpha \to 0$simultaneously, the limit depends on their joint scaling.

Lemma (limits and continuity). The ratio

$${\displaystyle \frac{d}{\sin \alpha }}$$

and therefore, the projected Bennett factor

$$sin\left({\displaystyle \frac{{\theta }_i}{2}}\right)\bullet sin\left({\displaystyle \frac{{\theta }_{i+1}}{2}}\right)\bullet {\displaystyle \frac{d}{sin \alpha }}$$

admit finite limits in the following cases: when$d\to 0$while$\mathit{sin}\alpha$stays bounded away from zero; or when$\alpha \to 0$under the bounded-ratio scaling

$$d=O\left(\mathrm{s}\mathrm{i}\mathrm{n}\alpha \right).$$

Without such a scaling assumption, the ratio need not have a finite limit.

2.2. Geometry of skew axes (complete proof)

Lemma 1 (skew-axes identity).Let$L_1=(p_1,u)$and$L_2=(p_2,v)$be two non-parallel lines in${\mathbb{R}}^3$, with$\parallel u\parallel =\parallel v\parallel =1$. Let$\alpha$be the acute angle between their directions, so that$\parallel u\times v\parallel =sin\alpha >0$. Define

$$\nu ={\displaystyle \frac{u\times v}{\parallel u\times v\parallel }}$$

Then the minimal distance $d$between the two lines is

$$d=\mid (p_2-p_1)\cdot \nu \mid ={\displaystyle \frac{\mid (p_2-p_1)\cdot (u\times v)\mid }{\parallel u\times v\parallel }}$$

Equivalently,

$$\parallel u\times v\parallel \hspace{0.17em}d=\mid (p_2-p_1)\cdot (u\times v)\mid$$

Proof. Since $u$and $v$are not parallel, $u\times v\ne 0$, and $\nu$is a unit vector orthogonal to both directions. The shortest segment connecting two skew lines is orthogonal to both lines, hence it is parallel to $\nu$. Therefore, the distance between the two lines is the absolute value of the projection of any connecting vector $p_2-p_1$onto $\nu$. This projection is independent of the chosen points on the two lines: if $p_1$is replaced by $p_1+\lambda u$and $p_2$by $p_2+\mu v$, then the projection changes by $(\mu v-\lambda u)\cdot \nu =0$. Hence

$$d=\mid (p_2-p_1)\cdot \nu \mid$$

Substituting gives

$$d={\displaystyle \frac{\mid (p_2-p_1)\cdot (u\times v)\mid }{\parallel u\times v\parallel }}$$

Since $\parallel u\times v\parallel =\mathrm{s}\mathrm{i}\mathrm{n}\alpha$, the stated identity follows.

2.3. DQ for a pure rotation about a line (full derivation)

Let a rotation by angle $\mathsf{\theta }$ about axis $(\mathit{p},\mathit{u})$ be given. Set

$$\mathrm{r}(\mathsf{\theta },\mathit{u})=\mathrm{c}\mathrm{o}\mathrm{s}\mathsf{\theta }/2+\mathrm{s}\mathrm{i}\mathrm{n}\mathsf{\theta }/2\hspace{0.17em}\mathit{u}, \mathrm{R}=\mathrm{R}(\mathsf{\theta },\mathit{u})\in \mathrm{S}\mathrm{O}\left(3\right)$$

The rigid motion “rotate about the line through $\mathit{p}$ with direction $\mathit{u}$ by $\mathsf{\theta }$” equals $(\mathrm{R},\mathit{t})$ with translation

$$\mathit{t}=(\mathrm{I}-\mathrm{R})\hspace{0.17em}\mathit{p}=\mathit{p}-\mathrm{R}\mathit{p}\hspace{0.25em}$$

(5)

Rodrigues expansion. For any $\mathit{x}$,

$$\mathrm{R}\mathit{x}=\mathit{x}\mathrm{c}\mathrm{o}\mathrm{s}\mathsf{\theta }+(\mathit{u}\times \mathit{x})\mathrm{s}\mathrm{i}\mathrm{n}\mathsf{\theta }+\mathit{u}(\mathit{u}\cdot \mathit{x})(1-\mathrm{c}\mathrm{o}\mathrm{s}\mathsf{\theta })$$

Hence, with $\mathit{x}=\mathit{p}$ and ${\mathit{p}}_{\perp }:=\mathit{p}-\mathit{u}(\mathit{u}\cdot \mathit{p})$ the component orthogonal to $\mathit{u}$,

$$\mathit{t}=(1-\mathrm{c}\mathrm{o}\mathrm{s}\mathsf{\theta })\hspace{0.17em}{\mathit{p}}_{\perp }-\mathrm{s}\mathrm{i}\mathrm{n}\mathsf{\theta }\hspace{0.17em}(\mathit{u}\times \mathit{p})=2\mathrm{s}\mathrm{i}\mathrm{n}\mathsf{\theta }/2[\mathrm{s}\mathrm{i}\mathrm{n}\mathsf{\theta }/2\hspace{0.17em}{\mathit{p}}_{\perp }-\mathrm{c}\mathrm{o}\mathrm{s}\mathsf{\theta }/2(\mathit{u}\times \mathit{p}\left)\right]$$

(6)

Unit DQ of the motion. Using (2),

$$\mathrm{q}(\mathit{p},\mathit{u},\mathsf{\theta })=\mathrm{r}+\mathsf{\epsilon }\hspace{0.17em}1/2\hspace{0.17em}\mathit{t}\hspace{0.17em}\mathrm{r}, \mathrm{r}=\mathrm{r}(\mathsf{\theta },\mathit{u})$$

(7)

Check (1): $r\bullet r^{*}=1$ by construction. With $d={\displaystyle \frac{t\bullet r}{2}}$ we have $d^{*}={\left({\displaystyle \frac{t\bullet r}{2}}\right)}^{*}={\displaystyle \frac{r^{*}\bullet t^{*}}{2}}=-{\displaystyle \frac{r^{*}\bullet t}{2}}$ (since $t^{*}=-t$ for pure quaternions). Therefore

$$r\bullet d^{*}+d\bullet r^{*}=r\bullet \left(-{\displaystyle \frac{r^{*}\bullet t}{2}}\right)+{\displaystyle \frac{t\bullet r}{2}}\bullet r^{*}=-{\displaystyle \frac{r\bullet r^{*}\bullet t}{2}}+{\displaystyle \frac{t\bullet r\bullet r^{*}}{2}}=0$$

Thus (1) holds identically. Geometrically, p may be any point on the rotation axis; shifting $\mathit{p}$ along $\mathit{u}$ does not change $\mathit{t}$, since $(I-R)\mathit{u}=0$Remark 1. The point $\mathit{p}$ is fixed by the full rigid motion, i.e., $R\mathit{p}+\mathit{t}=\mathit{p}$

2.4. 4R Chain in DQ and Exact Closure Conditions

Consider a spatial 4R loop with successive axes ${\mathrm{L}}_{\mathrm{i}}=({\mathit{p}}_{\mathrm{i}},{\mathit{u}}_{\mathrm{i}})$, angles ${\mathsf{\theta }}_{\mathrm{i}}$, and unit DQs ${\mathrm{q}}_{\mathrm{i}}={\mathrm{r}}_{\mathrm{i}}+\mathsf{\epsilon }{\mathrm{d}}_{\mathrm{i}}$ with

$${\mathrm{r}}_{\mathrm{i}}=\mathrm{r}({\mathsf{\theta }}_{\mathrm{i}},{\mathit{u}}_{\mathrm{i}}), {\mathrm{d}}_{\mathrm{i}}=1/2\hspace{0.17em}{\mathit{t}}_{\mathrm{i}}\hspace{0.17em}{\mathrm{r}}_{\mathrm{i}}, {\mathit{t}}_{\mathrm{i}}=(\mathrm{I}-{\mathrm{R}}_{\mathrm{i}})\hspace{0.17em}{\mathit{p}}_{\mathrm{i}}$$

The DQ product is

$$\mathrm{Q}≔{\mathrm{q}}_4{\mathrm{q}}_3{\mathrm{q}}_2{\mathrm{q}}_1=\left({\mathrm{r}}_4{\mathrm{r}}_3{\mathrm{r}}_2{\mathrm{r}}_1\right)\hspace{0.17em}+\hspace{0.17em}\mathsf{\epsilon }\hspace{0.17em}\mathrm{D}$$

$$\mathrm{D}={\mathrm{d}}_4{\mathrm{r}}_3{\mathrm{r}}_2{\mathrm{r}}_1+{\mathrm{r}}_4{\mathrm{d}}_3{\mathrm{r}}_2{\mathrm{r}}_1+{\mathrm{r}}_4{\mathrm{r}}_3{\mathrm{d}}_2{\mathrm{r}}_1+{\mathrm{r}}_4{\mathrm{r}}_3{\mathrm{r}}_2{\mathrm{d}}_1$$

For $q_i=r_i+\epsilon \bullet d_i$, with $d_i={\displaystyle \frac{t_i\bullet r_i}{2}}$. Then:

$$Q=q_4\bullet q_3\bullet q_2\bullet q_1=r_4\bullet r_3\bullet r_2\bullet r_1+\epsilon \bullet D$$

Where

$$D=r_4\bullet r_3\bullet r_2\bullet d+r_4\bullet r_3\bullet d_2\bullet r_1+r_4\bullet d_3\bullet r_2\bullet r_1+d_4\bullet r_3\bullet r_2\bullet r_1$$

Using the relation for pure vectors v:

$$r\bullet v=\left(A\bullet d_r\bullet v\right)\bullet r$$

we obtain:

$$D={\displaystyle \frac{t_4+R_4\bullet t_3+R_4\bullet R_3\bullet t_2+R_4\bullet R_3\bullet R_2\bullet t_1}{2\bullet r_4\bullet r_3\bullet r_2\bullet r_1}}$$

Exact 4R closure means.

Remark.Note that in the quaternion representation, q = ±1 both correspond to the identity rotation in SO(3); the sign ambiguity does not affect the closure condition.

$${\mathrm{r}}_4{\mathrm{r}}_3{\mathrm{r}}_2{\mathrm{r}}_1=\pm 1, \mathrm{D}=0$$

(8)

The primal condition is the spherical (orientation) closure; the dual condition is the translational closure.

2.4.1. Orientation closure

Let $R_i$be the rotations associated with $r_i$. The primal equation $r_4{r}_3{r}_2{r}_1=\pm 1$is equivalent to $R_4{R}_3{R}_2{R}_1=I$, since $r$and $-r$represent the same rotation in $SO\left(3\right)$. This condition is the orientation part of the exact 4R closure. It is analogous to spherical 4R closure at the level of rotations, but the spatial axes considered here are not assumed to be concurrent.

2.4.2. Dual part: exact expansion and transport

Using ${\mathrm{d}}_{\mathrm{i}}=1/2\hspace{0.17em}{\mathit{t}}_{\mathrm{i}}{\mathrm{r}}_{\mathrm{i}}$ and associativity,

$$\begin{array}{cc}\mathrm{D}& =1/2({\mathit{t}}_4{\mathrm{r}}_4{\mathrm{r}}_3{\mathrm{r}}_2{\mathrm{r}}_1+{\mathrm{r}}_4\hspace{0.17em}{\mathit{t}}_3{\mathrm{r}}_3{\mathrm{r}}_2{\mathrm{r}}_1+{\mathrm{r}}_4{\mathrm{r}}_3\hspace{0.17em}{\mathit{t}}_2{\mathrm{r}}_2{\mathrm{r}}_1+{\mathrm{r}}_4{\mathrm{r}}_3{\mathrm{r}}_2\hspace{0.17em}{\mathit{t}}_1{\mathrm{r}}_1)\\ & =1/2({\mathit{t}}_4+{\mathrm{A}\mathrm{d}}_{{\mathrm{r}}_4}({\mathit{t}}_3)+{\mathrm{A}\mathrm{d}}_{{\mathrm{r}}_4{\mathrm{r}}_3}({\mathit{t}}_2)+{\mathrm{A}\mathrm{d}}_{{\mathrm{r}}_4{\mathrm{r}}_3{\mathrm{r}}_2}({\mathit{t}}_1\left)\right)\hspace{0.17em}\left({\mathrm{r}}_4{\mathrm{r}}_3{\mathrm{r}}_2{\mathrm{r}}_1\right)\end{array}$$

(9)

where ${\mathrm{A}\mathrm{d}}_{\mathrm{r}}\left(\mathit{v}\right):=\mathrm{r}\hspace{0.17em}\mathit{v}\hspace{0.17em}{\mathrm{r}}^{\mathrm{*}}$ equals the action of the corresponding rotation on vectors. Under (8) (primal closure), ${\mathrm{r}}_4{\mathrm{r}}_3{\mathrm{r}}_2{\mathrm{r}}_1=\pm 1$, hence the dual closure is equivalent to the vector identity

$$\mathit{S}:={\mathit{t}}_4+{\mathrm{R}}_4\hspace{0.17em}{\mathit{t}}_3+{\mathrm{R}}_4{\mathrm{R}}_3\hspace{0.17em}{\mathit{t}}_2+{\mathrm{R}}_4{\mathrm{R}}_3{\mathrm{R}}_2\hspace{0.17em}{\mathit{t}}_1=0$$

(10)

Thus, once orientation closure is imposed, the dual part of the exact closure condition reduces to the Euclidean transported-translation equation above.

2.5. Bennett Condition via Normal Projections

We show that, under (H1)–(H4) and the paired-normal compatibility assumptions stated below, the projected dual-closure equation (10) leads to the classical Bennett ratio conditions on opposite edges [22,23]:

$${\displaystyle \frac{{\mathrm{d}}_{12}}{\mathrm{s}\mathrm{i}\mathrm{n}{\mathsf{\alpha }}_{12}}}={\displaystyle \frac{{\mathrm{d}}_{23}}{\mathrm{s}\mathrm{i}\mathrm{n}{\mathsf{\alpha }}_{23}}}={\displaystyle \frac{{\mathrm{d}}_{34}}{\mathrm{s}\mathrm{i}\mathrm{n}{\mathsf{\alpha }}_{34}}}={\displaystyle \frac{{\mathrm{d}}_{41}}{\mathrm{s}\mathrm{i}\mathrm{n}{\mathsf{\alpha }}_{41}}}$$

(11)

The key step is to project $\mathit{S}$ onto the four pairwise normals ${\nu }_{12},{\nu }_{23},{\nu }_{34},{\nu }_{41}$, while keeping track of the transport of opposite normal directions.

2.5.1. Preparatory identities

For each joint axis $L_i=({\mathit{p}}_i,{\mathit{u}}_i)$, define the component of ${\mathit{p}}_i$perpendicular to the axis direction by

$${\mathit{p}}_{i,\perp }:={\mathit{p}}_i-{\mathit{u}}_i({\mathit{u}}_i\cdot {\mathit{p}}_i)$$

From the rotation-about-a-line formula (6), the translation vector associated with joint $i$is

$$\begin{array}{cccc}& {\mathit{t}}_i=(1-\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_i){\mathit{p}}_{i,\perp }-\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_i({\mathit{u}}_i\times {\mathit{p}}_i)& & \end{array}$$

(12)

For the adjacent pair $(i,i+1)$, let ${\nu }_{i,i+1}$be the unit normal defined in (4). Then

$$\begin{array}{cccc}& {\nu }_{i,i+1}\perp {\mathit{u}}_i,{\nu }_{i,i+1}\perp {\mathit{u}}_{i+1}& & \end{array}$$

(13)

We shall also use the standard scalar triple-product identities

$$\begin{array}{cccc}& (\mathit{a}\times \mathit{b})\cdot \mathit{c}=\mathit{a}\cdot (\mathit{b}\times \mathit{c})=\mathit{b}\cdot (\mathit{c}\times \mathit{a})& & \end{array}$$

(14)

Assumption 1 (paired-normal transport compatibility).In addition to the orientation closure$R_4{R}_3{R}_2{R}_1=I$, we assume that the transported normals of opposite adjacent-axis pairs are compatible, namely

$$R_4{R}_3{\nu }_{12}={\sigma }_{34}{\nu }_{34},R_1{R}_4{\nu }_{23}={\sigma }_{41}{\nu }_{41}$$

where ${\sigma }_{34},{\sigma }_{41}\in \{+1,-1\}$.

This assumption states that the normal direction associated with one adjacent pair of axes is transported to the normal direction associated with the opposite adjacent pair. The signs only reflect the chosen orientations of the unit normals and do not affect the Bennett ratio equations, which involve unsigned distances.

Remark.This paired-normal compatibility is not a consequence of orientation closure alone. It is an additional geometric compatibility condition required for the paired projection argument used below.

2.5.2. Choice of reference points (gauge) and its invariance

For each adjacent pair $(i,i+1)$, choose points $p_i\in L_i$and $p_{i+1}\in L_{i+1}$realizing the shortest common-normal segment between the two axes. Thus,

$$\begin{array}{cccc}& p_{i+1}-p_i={\sigma }_{i,i+1}{d}_{i,i+1}{\nu }_{i,i+1},{\sigma }_{i,i+1}\in \{+1,-1\}& & \end{array}$$

(15)

Such a choice is possible for non-parallel skew axes under the standing assumptions. The sign ${\sigma }_{i,i+1}$depends only on the chosen orientation of the unit normal ${\nu }_{i,i+1}$.

Lemma (gauge invariance).If the reference point on axis$L_i$is shifted along the axis,

$$p_i\mapsto p_i+\lambda u_i$$

then

$$(I-R_i)u_i=0$$

because $u_i$is the rotation axis of $R_i$. Hence

$$t_i=(I-R_i)p_i$$

is unchanged under such shifts. Therefore, the dual-quaternion representation of a rotation about a given line is independent of the chosen point on that line.

Moreover, for an adjacent pair, the scalar normal separation is also gauge-invariant:

$$(p_{i+1}-p_i)\cdot {\nu }_{i,i+1}={\sigma }_{i,i+1}{d}_{i,i+1}$$

because ${\nu }_{i,i+1}\perp u_i$and ${\nu }_{i,i+1}\perp u_{i+1}$. Thus, shifts of $p_i$or $p_{i+1}$along their respective axes do not change the normal component of the separation, although they may change the full connecting vector.

2.5.3. Transport structure and pair grouping

Starting from the transported dual-closure equation

$$\begin{array}{cccc}& \mathit{S}={\mathit{t}}_4+R_4{\mathit{t}}_3+R_4{R}_3{\mathit{t}}_2+R_4{R}_3{R}_2{\mathit{t}}_1=0& & \end{array}$$

(16)

we group the terms according to opposite adjacent-axis pairs. For the pair $\left(1\right|2)$ and the opposite pair $\left(3\right|4)$, write

$$\begin{array}{cccc}& \mathit{S}=\underset{{\mathit{S}}_{12}}{\underbrace{\left(R_4{R}_3{R}_2{\mathit{t}}_1\right|R_4{R}_3{\mathit{t}}_2)}}+\underset{{\mathit{S}}_{34}}{\underbrace{\left(R_4{\mathit{t}}_3\right|{\mathit{t}}_4)}}& & \end{array}$$

(17)

Projecting $\mathit{S}=0$ onto ${\nu }_{12}$ gives

$$\begin{array}{cccc}& \left(R_4{R}_3{R}_2{\mathit{t}}_1\right|R_4{R}_3{\mathit{t}}_2)\cdot {\nu }_{12}=-(R_4{\mathit{t}}_3\left|{\mathit{t}}_4\right)\cdot {\nu }_{12}& & \end{array}$$

(18)

Equivalently,

$${\mathit{S}}_{12}\cdot {\nu }_{12}=-{\mathit{S}}_{34}\cdot {\nu }_{12}$$

The left-hand side represents the transported contribution of the adjacent pair $\left(1\right|2)$, while the right-hand side represents the transported contribution of the opposite pair $\left(3\right|4)$. To compare these contributions with the Bennett normal of the opposite pair, we use the paired-normal transport compatibility assumption,

$$R_4{R}_3{\nu }_{12}={\sigma }_{34}{\nu }_{34},{\sigma }_{34}\in \{+1,-1\}$$

This allows the projection of the transported $\left(3\right|4)$ contribution onto ${\nu }_{12}$to be equivalently interpreted, up to sign, as a projection onto ${\nu }_{34}$. The sign depends only on the chosen orientations of the normals and does not affect the unsigned Bennett ratios.

Repeating the same argument cyclically for the pairs $\left(2\right|3)$ and $\left(4\right|1)$ gives the second paired projection relation. Thus, the projection structure compares opposite adjacent-axis pairs through the transported dual-closure equation.

2.5.4. Projected pair contributions

We now evaluate the projected contributions of opposite adjacent-axis pairs under the gauge choice (4.4) and the paired-normal transport compatibility assumption.

For the adjacent pair $\left(1\right|2)$, define

$${\mathit{S}}_{12}=R_4{R}_3{R}_2{\mathit{t}}_1+R_4{R}_3{\mathit{t}}_2$$

For the opposite adjacent pair $\left(3\right|4)$, define

$${\mathit{S}}_{34}=R_4{\mathit{t}}_3+{\mathit{t}}_4$$

Then the transported dual-closure equation $\mathit{S}=0$ gives

$${\mathit{S}}_{12}+{\mathit{S}}_{34}=0$$

Projecting this equation onto ${\nu }_{12}$, we obtain

$$\begin{array}{cccc}& {\mathit{S}}_{12}\cdot {\nu }_{12}=-{\mathit{S}}_{34}\cdot {\nu }_{12}& & \end{array}$$

(19)

Here ${\mathit{S}}_{12}\cdot {\nu }_{12}$ denotes the transported paired-projection contribution associated with the adjacent-axis pair (1∣2), evaluated under the gauge choice ${\mathit{p}}_2-{\mathit{p}}_1={\sigma }_{12}{d}_{12}{\nu }_{12}$.

Using the rotation-about-a-line formula

$${\mathit{t}}_i=(1-\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_i){\mathit{p}}_{i,\perp }-\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_i({\mathit{u}}_i\times {\mathit{p}}_i)$$

together with the gauge condition

$${\mathit{p}}_2-{\mathit{p}}_1={\sigma }_{12}{d}_{12}{\nu }_{12}$$

the scalar triple-product identities, and the half-angle identities

$$1-cos{\theta }_i=2\bullet {sin}^2\left({\displaystyle \frac{{\theta }_i}{2}}\right), sin{\theta }_i=2\bullet sin\left({\displaystyle \frac{{\theta }_i}{2}}\right)\bullet cos\left({\displaystyle \frac{{\theta }_i}{2}}\right)$$

the projected contribution of the pair $\left(1\right|2)$takes the form

$${\mathit{S}}_{12}\cdot {\nu }_{12}={\sigma }_{12}=sin\left({\displaystyle \frac{{\theta }_1}{2}}\right)\bullet sin\left({\displaystyle \frac{{\theta }_2}{2}}\right)\bullet {\displaystyle \frac{d_{12}}{\sin {\mathsf{\alpha }}_{12}}}$$

(20)

Here ${\sigma }_{12}\in \{+1,-1\}$depends only on the chosen orientation of the normal ${\nu }_{12}$.

Similarly, using the paired-normal transport compatibility for the opposite pair $\left(3\right|4)$, the projected contribution of ${\mathit{S}}_{34}$can be expressed, up to an orientation sign, as

$${\mathit{S}}_{34}\cdot {\nu }_{34}=-{\sigma }_{34}=sin\left({\displaystyle \frac{{\theta }_3}{2}}\right)\bullet sin\left({\displaystyle \frac{{\theta }_4}{2}}\right)\bullet {\displaystyle \frac{d_{34}}{\sin {\mathsf{\alpha }}_{34}}}$$

(21)

Substitution of (10) and (21) into (19), and absorption of orientation signs into the chosen normal directions, gives

$$sin\left({\displaystyle \frac{{\theta }_1}{2}}\right)\bullet sin\left({\displaystyle \frac{{\theta }_2}{2}}\right)\bullet {\displaystyle \frac{d_{12}}{\sin {\mathsf{\alpha }}_{12}}}=sin\left({\displaystyle \frac{{\theta }_3}{2}}\right)\bullet sin\left({\displaystyle \frac{{\theta }_4}{2}}\right)\bullet {\displaystyle \frac{d_{34}}{\sin {\mathsf{\alpha }}_{34}}}$$

(22)

Repeating the same argument cyclically for the opposite adjacent-axis pairs $\left(2\right|3)$ and $\left(4\right|1)$ yields

$$sin\left({\displaystyle \frac{{\theta }_2}{2}}\right)\bullet sin\left({\displaystyle \frac{{\theta }_3}{2}}\right)\bullet {\displaystyle \frac{d_{23}}{\sin {\mathsf{\alpha }}_{23}}}=sin\left({\displaystyle \frac{{\theta }_4}{2}}\right)\bullet sin\left({\displaystyle \frac{{\theta }_1}{2}}\right)\bullet {\displaystyle \frac{d_{41}}{\sin {\mathsf{\alpha }}_{41}}}$$

(23)

Therefore, the projected dual-closure relations are (22) and (23)

If, in addition, an angle-symmetric closure branch is imposed,

$$\begin{array}{cccc}& {\theta }_1={\theta }_3\ne 0,{\theta }_2={\theta }_4\ne 0& & \end{array}$$

(24)

then the nonzero half-angle factors cancel in (22) and (23), and the relations reduce to

$${\displaystyle \frac{d_{12}}{\sin {\alpha }_{12}}}={\displaystyle \frac{d_{34}}{\sin {\alpha }_{34}}},{\displaystyle \frac{d_{23}}{\sin {\alpha }_{23}}}={\displaystyle \frac{d_{41}}{\sin {\alpha }_{41}}}$$

Together with the cyclic projected relations, this gives the Bennett ratio condition

$$\begin{array}{cccc}& {\displaystyle \frac{d_{12}}{\sin {\alpha }_{12}}}={\displaystyle \frac{d_{23}}{\sin {\alpha }_{23}}}={\displaystyle \frac{d_{34}}{\sin {\alpha }_{34}}}={\displaystyle \frac{d_{41}}{\sin {\alpha }_{41}}}& & \end{array}$$

(25)

Conversely, if the Bennett ratio condition holds, then the projected dual-closure equations admit the angle-symmetric nontrivial branch ${\theta }_1={\theta }_3$, ${\theta }_2={\theta }_4$, subject to the orientation closure and paired-normal compatibility assumptions. Under these assumptions, the projected dual-closure relations associated with (16) are satisfied.

Theorem 1 (Bennett–DQ correspondence for a non-degenerate spatial 4R under paired projection compatibility).

Consider a non-degenerate spatial 4R chain with successive revolute axes$L_i=({\mathit{p}}_i,{\mathit{u}}_i)$,$i=1,\dots ,4$, satisfying the assumptions$\left(H1\right)$–$\left(H4\right)$. Let

$$q_i=r_i+\epsilon d_i,d_i={\displaystyle \frac{1}{2}}{\mathit{t}}_i{r}_i,{\mathit{t}}_i=(I-R_i){\mathit{p}}_i,$$

be the unit dual quaternion associated with the rotation of angle${\theta }_i$about the axis$L_i$. Assume that the primal, or orientation, closure condition holds,

$$R_4{R}_3{R}_2{R}_1=I$$

and assume, in addition, the paired-normal transport compatibility conditions

$$R_4{R}_3{\nu }_{12}={\sigma }_{34}{\nu }_{34},R_1{R}_4{\nu }_{23}={\sigma }_{41}{\nu }_{41},{\sigma }_{34},{\sigma }_{41}\in \{+1,-1\}$$

Here

$${\nu }_{i,i+1}={\displaystyle \frac{{\mathit{u}}_i\times {\mathit{u}}_{i+1}}{\parallel {\mathit{u}}_i\times {\mathit{u}}_{i+1}\parallel }}$$

is the unit normal associated with the adjacent pair of axes $\left(L_i\right|L_{i+1})$, and $d_{i,i+1}$and ${\alpha }_{i,i+1}$denote respectively the shortest distance and skew angle between those axes.

Under these assumptions, the projected dual-closure equations associated with

$$\mathit{S}={\mathit{t}}_4+R_4{\mathit{t}}_3+R_4{R}_3{\mathit{t}}_2+R_4{R}_3{R}_2{\mathit{t}}_1=0$$

reduce to the opposite-pair relations

$$\mathrm{s}\mathrm{i}\mathrm{n}{\displaystyle \frac{{\theta }_1}{2}}\mathrm{s}\mathrm{i}\mathrm{n}{\displaystyle \frac{{\theta }_2}{2}}{\displaystyle \frac{d_{12}}{\sin {\alpha }_{12}}}=\mathrm{s}\mathrm{i}\mathrm{n}{\displaystyle \frac{{\theta }_3}{2}}\mathrm{s}\mathrm{i}\mathrm{n}{\displaystyle \frac{{\theta }_4}{2}}{\displaystyle \frac{d_{34}}{\sin {\alpha }_{34}}}$$

and

$$\mathrm{s}\mathrm{i}\mathrm{n}{\displaystyle \frac{{\theta }_2}{2}}\mathrm{s}\mathrm{i}\mathrm{n}{\displaystyle \frac{{\theta }_3}{2}}{\displaystyle \frac{d_{23}}{\sin {\alpha }_{23}}}=\mathrm{s}\mathrm{i}\mathrm{n}{\displaystyle \frac{{\theta }_4}{2}}\mathrm{s}\mathrm{i}\mathrm{n}{\displaystyle \frac{{\theta }_1}{2}}{\displaystyle \frac{d_{41}}{\sin {\alpha }_{41}}}$$

If, moreover, the chain is restricted to the angle-symmetric nontrivial branch

$${\theta }_1={\theta }_3\ne 0,{\theta }_2={\theta }_4\ne 0,$$

then the nonzero half-angle factors cancel, and the projected dual-closure equations reduce to the Bennett axis-ratio condition

$${\displaystyle \frac{d_{12}}{\sin {\alpha }_{12}}}={\displaystyle \frac{d_{23}}{\sin {\alpha }_{23}}}={\displaystyle \frac{d_{34}}{\sin {\alpha }_{34}}}={\displaystyle \frac{d_{41}}{\sin {\alpha }_{41}}}.$$

Conversely, if the Bennett ratio condition holds, then, under the same orientation-closure and paired-normal transport compatibility assumptions, the projected dual-closure equations admit the angle-symmetric branch

$${\theta }_1={\theta }_3,{\theta }_2={\theta }_4,$$

and hence satisfy the projected form of the dual-closure equation within the present paired-projection framework.

Remark 2.This theorem establishes a correspondence between the Bennett ratio and the projected dual-quaternion closure relations under the stated assumptions. It does not claim that orientation closure alone implies paired-normal transport compatibility. It also does not treat degenerate cases as non-degenerate Bennett 4R chains. In particular, cases such as$d\to 0$or$\alpha \to 0$are outside$\left(H1\right)$–$\left(H4\right)$and must be treated as limiting regimes. For the Kresling application, the relevant use is the intersecting-axis Bennett-limit diagnostic$d=0$, not the assertion that every Kresling quadrilateral is a non-degenerate Bennett 4R linkage.

2.6. Kresling TRI cells: local closure and global compatibility

Corollary 1(local closure per quad and Kresling limiting interpretation). In the rigid regime, a non-classical triangular Kresling (TRI) cell can be decomposed into elementary quadrilateral units formed by two adjacent triangular panels and two ring edges. In the non-degenerate reference case, each such quadrilateral unit can be modeled as a spatial 4R loop, provided that its fold axes satisfy the assumptions (H1)–(H4).

If, for each non-degenerate quadrilateral unit, the Bennett ratios (25) hold and the additional assumptions used in Theorem 1 are satisfied — namely orientation closure, paired-normal transport compatibility, and the adopted angle-symmetric branch — then the corresponding local 4R loop admits a nontrivial DQ closure within the present paired-projection framework.

In the rigid Kresling validation considered below, the local quadrilateral units are evaluated in the intersecting-axis Bennett limit $(\mathrm{d}=0)$. Therefore, they are not assumed to be non-degenerate Bennett 4R linkages satisfying (H2). The non-degenerate 4R result above provides the reference paired-projection framework, while the Kresling application uses its limiting diagnostic form. Under this interpretation, Bennett-limit compatibility is used as a local diagnostic of the rigid Kresling crease geometry, not as a claim that every Kresling quadrilateral is a non-degenerate Bennett 4R linkage.

Global angle compatibility.Let$\Theta \in {\mathbb{R}}^E$collect the fold angles associated with the edges of the Kresling assembly, and let$\mathcal{C}$denote the set of elementary cycles. For each cycle$c\in \mathcal{C}$, define the dual-quaternion closure product

$${\mathsf{\Phi }}_c\left(\mathsf{\Theta }\right)$$

Exact global compatibility requires

$${\mathsf{\Phi }}_c\left(\mathsf{\Theta }\right)=1 \mathrm{for} \mathrm{all} c\in \mathcal{C}$$

Equivalently, the primal and dual parts of each cycle product must satisfy the corresponding orientation and translational closure conditions.

Linearizing these closure constraints at a locally closed configuration ${\mathsf{\Theta }}_0$gives a block-incidence Jacobian system of the form

$$J_{\psi }\left({\mathsf{\Theta }}_0\right)\hspace{0.17em}\delta \mathsf{\Theta }=0,$$

where $J_{\psi }$stacks the derivatives of the dual-part vector residuals over the independent cycles. A nontrivial kernel of $J_{\psi }$indicates the existence of infinitesimal angle variations preserving the linearized cycle-closure constraints.

Thus, local feasibility of a global assignment is characterized, at the linearized level, by the consistency of the cycle-incidence constraints. This is the structure used in the numerical solver: DQ composition per cycle, assembly of the cycle residuals, and a column-balanced Jacobian for numerical conditioning.

2.7. Natural extensions (Derivation outline)

2.7.1. Nonzero pitch screw joints

The preceding derivation was restricted to zero-pitch revolute joints. A natural extension is to allow screw joints with nonzero pitch. Let joint $i$have pitch $h_i$along the unit direction $u_i$. For a rotation angle ${\theta }_i$, the associated rigid displacement may be written using the effective translation vector

$$T_i=(I-R_i)p_i+h_i{\theta }_i{u}_i.$$

The corresponding unit dual quaternion is then expressed as

$$q_i=r_i+\epsilon {\displaystyle \frac{1}{2}}{T}_i{r}_i,$$

where $r_i$is the primal unit quaternion associated with the rotation $R_i$. For $h_i=0$, one obtains $T_i=(I-R_i)p_i=t_i$, and the pure revolute formulation used above is recovered.

Under orientation closure, the dual closure equation becomes

$$R_4{R}_3{R}_2{T}_1+R_4{R}_3{T}_2+R_4{T}_3+T_4=0.$$

Equivalently, the translational residual is obtained from the pure-revolute residual by replacing each $t_i$with

$$T_i=t_i+h_i{\theta }_i{u}_i.$$

Thus, compared with the zero-pitch case, the normal projections of the dual-closure equation acquire additional pitch-dependent terms. These terms depend on the pitch convention and on the adopted screw parametrization. A full derivation of the resulting pitch-dependent projected compatibility relations is outside the scope of the present paper. In the zero-pitch limit $h_i=0$, all pitch-dependent terms vanish and the pure-revolute Bennett projection relations derived above are recovered.

2.7.2. Stability: linearization, conditioning, and TSVD

Let $\psi \left(\mathsf{\Theta }\right)$denote the stacked dual-part residuals, with three scalar constraints per independent cycle. Linearizing around a reference configuration ${\mathsf{\Theta }}_0$ gives

$$\delta \psi =J_{\psi }\left({\mathsf{\Theta }}_0\right)\hspace{0.17em}\delta \mathsf{\Theta }.$$

Let the singular value decomposition of the Jacobian be

$$J_{\psi }=U\mathsf{\Sigma }{V}^{\mathrm{\top }}.$$

For the least-squares linear correction, the Moore–Penrose solution is

$$\delta \mathsf{\Theta }=J_{\psi }^{+}\delta \psi =V{\mathsf{\Sigma }}^{+}{U}^{\mathrm{\top }}\delta \psi .$$

If ${\sigma }_{\mathrm{m}\mathrm{i}\mathrm{n}}$denotes the smallest retained nonzero singular value, then

$$\parallel \delta \mathsf{\Theta }\parallel \le \parallel J_{\psi }^{+}\parallel \hspace{0.17em}\parallel \delta \psi \parallel =({\sigma }_{\mathrm{m}\mathrm{i}\mathrm{n}}{)}^{-1}\parallel \delta \psi \parallel .$$

Thus, small singular values amplify residual perturbations and indicate poor local conditioning.

To improve numerical stability, truncated SVD replaces ${\mathsf{\Sigma }}^{+}$by a truncated inverse ${\mathsf{\Sigma }}_{\tau }^{+}$, where singular values below a threshold $\tau$are discarded rather than inverted. The resulting correction satisfies

$$\parallel \delta {\mathsf{\Theta }}_{\tau }\parallel \le {\tau }^{-1}\parallel \delta \psi \parallel$$

on the retained singular subspace. This bound reflects numerical stabilization, but the discarded singular directions introduce a truncation residual. Hence TSVD trades exact linear least-squares correction for improved robustness against near-singular directions.

Finally, column balancing, obtained by normalizing the columns of $J_{\psi }$, reduces artificial scaling differences between angular variables and typically improves the numerical conditioning of the linearized system.

3. Results and Numerical Validation

3.1. Methodology: DQ-Parametric Homotopy for Rigid Kresling FOLD–TWIST–FOLD Validation

The systems obtained from dual-quaternion modeling of rigid origami mechanisms are nonlinear because rigid-body rotations, crease-axis constraints, panel compatibility, and closure conditions are coupled through products of unit dual quaternions. For this reason, solving only the final configuration in a single step may be sensitive to the initial guess, especially when multiple kinematic branches or trivial zero-motion configurations are present.

To avoid this ambiguity, the present numerical validation uses a DQ-parametric homotopy/continuation strategy. Instead of treating only the final FOLD–TWIST–FOLD configuration, the motion is embedded in a continuous path parameterized by a homotopy variable λ ∈ [0, 1]. The path starts from the known reference configuration and continuously tracks the prescribed sequence: first fold, then twist, then fold. At each homotopy state, the rigid Kresling origami configuration is reconstructed from the constrained nodal geometry, and the associated rigid transformations are represented using unit dual quaternions in SE(3).

The numerical model is a three-cell hexagonal Kresling origami robot. Each cell is represented by rigid triangular panels connected through fixed crease lines. The validation explicitly preserves the ring edges, primary creases, diagonal creases, and triangular-panel edge lengths. Therefore, the reported validation is not limited to a rigid-ring or helical-crease approximation; it enforces the geometric constraints required for a rigid Kresling origami model.

The FOLD–TWIST–FOLD command is prescribed as three sequential actuation targets: an axial fold in the first cell, a twist in the middle cell, and an axial fold in the third cell. For each homotopy state, the nodal configuration is solved subject to the rigid-origami constraints and the active command coordinate. The resulting panel and cell transformations are then encoded using unit dual quaternions. This construction allows the motion to be tracked as a continuous rigid-origami branch, while preserving the SE(3)-consistent description of each rigid component.

The trajectory is accepted only after independent diagnostic checks are satisfied. These include motion-command tracking, ring-edge preservation, primary-crease preservation, diagonal-crease preservation, triangular-panel rigidity, unit-DQ consistency, DQ panel-map consistency, Bennett-limit compatibility, homotopy residuals, and conditioning. Bennett-limit compatibility is evaluated in the intersecting-axis regime d = 0, which is appropriate for the local crease-axis geometry of the Kresling cell. The paired-normal diagnostic is treated separately from the Bennett-limit gate and is not used to claim that each Kresling quadrilateral is a non-degenerate Bennett 4R linkage.

For the reported configuration, the DQ-parametric homotopy evaluates 37 continuation states. The prescribed command is Δh1 = −2 mm, Δφ2 = 5°, and Δh3 = −2 mm. The final numerical run reports an overall PASS, with all rigid-origami, DQ, Bennett-limit, homotopy, and conditioning gates satisfied. The maximum physical homotopy residual remains of order 10−10, while unit-DQ consistency is maintained at order 10−16.

Accordingly, the final validation is not a staged Phase A/Phase B Newton procedure and is not a validation of a rigid-ring helical approximation alone. It is a direct DQ-parametric homotopy validation of a constrained rigid Kresling origami FOLD–TWIST–FOLD trajectory. Newton, Levenberg–Marquardt, TSVD, or other nonlinear solvers remain relevant for more general closure problems, but in the reported validation they are used only as numerical tools for tracking the constrained rigid-origami branch, followed by explicit post-checks of rigidity, DQ consistency, Bennett-limit compatibility, homotopy residuals, and conditioning. [35,36,37,38]

3.2. Why DQ-Based Homotopy is Used

The closure and compatibility equations associated with rigid Kresling origami mechanisms are nonlinear because rigid-body rotations, crease-axis constraints, panel rigidity, and multi-cell closure are coupled through products of unit dual quaternions. In the present model, the ring edges, primary creases, diagonal creases, and triangular-panel edge lengths must remain fixed throughout the motion. Therefore, directly solving only the final configuration may be sensitive to the initial guess, especially when multiple admissible kinematic branches or trivial zero-motion configurations are present.

The FOLD–TWIST–FOLD motion is not treated as an isolated final pose. Instead, it is embedded in a continuous homotopy path that starts from the reference configuration and tracks the intended physical sequence: fold, twist, and fold. This allows the numerical validation to follow the constrained rigid-origami branch throughout the trajectory, rather than checking only the endpoint. The homotopy path is generated by solving the Kresling nodal constraints at each continuation state, while the resulting rigid panel and cell transformations are represented using unit dual quaternions in SE(3).

Dual quaternions are used because each rigid panel and each rigidly reconstructed cell state corresponds to a rigid-body displacement in SE(3). A unit dual quaternion compactly represents the rotational and translational components of these displacements, while also enabling consistency checks through the unit-DQ and DQ panel-map residuals. Thus, the DQ representation is not used to model panel deformation; rather, it is used precisely because the adopted Kresling model is treated in the rigid-origami regime.

The trajectory is not accepted solely because the final command is reached. At each homotopy state, the configuration is checked through independent diagnostics: ring-edge preservation, primary-crease preservation, diagonal-crease preservation, triangular-panel rigidity, unit-DQ consistency, DQ panel-map consistency, Bennett-limit compatibility, homotopy residuals, and conditioning. For the reported configuration, all rigid-origami, DQ, Bennett-limit, homotopy, and conditioning gates are satisfied. The maximum physical homotopy residual remains of order 10−10, while unit-DQ consistency remains of order 10−16.

The Bennett condition is used as a local compatibility diagnostic rather than as a claim that each Kresling quadrilateral is a non-degenerate Bennett 4R linkage. Since the local Kresling crease-axis geometry includes intersecting-axis configurations, the relevant diagnostic regime is the Bennett limit d = 0. Along the reported homotopy trajectory, the intersecting-axis diagnostic remains satisfied, with all tested adjacent-axis pairs classified in the intersecting-axis regime.

Newton, Levenberg–Marquardt, TSVD, and related numerical methods remain standard tools for solving general nonlinear closure systems, especially when the configuration is unknown and must be found iteratively. In the present validation, however, the role of the numerical solver is to track the constrained rigid Kresling origami branch along the prescribed homotopy path. The final acceptance of the trajectory is based not on convergence alone, but on explicit post-checks of rigid-origami constraints, dual-quaternion consistency, Bennett-limit compatibility, homotopy residuals, and conditioning.

Accordingly, DQ-based homotopy is used here for two complementary purposes. First, it organizes the prescribed FOLD–TWIST–FOLD motion into a continuous physically meaningful path. Second, it provides an SE(3)-consistent framework in which the reconstructed rigid-origami configurations can be validated through DQ, rigidity, Bennett-limit, homotopy, and conditioning diagnostics over the full trajectory.

3.3. Why the DQ-Parametric Homotopy is Sufficient for the Reported Validation

In the reported validation, the FOLD–TWIST–FOLD trajectory is generated through a DQ-parametric homotopy combined with a constrained Kresling nodal reconstruction. The trajectory is not obtained by solving only an isolated final nonlinear system. Instead, the motion is embedded in a continuous continuation path that starts from the reference configuration and tracks the prescribed sequence: fold, twist, and fold. At each continuation state, the ring edges, primary creases, diagonal creases, and triangular-panel edge lengths are enforced, and the resulting rigid panel and cell transformations are represented using unit dual quaternions in SE(3).

For the tested three-cell hexagonal Kresling origami robot, the method evaluates 37 homotopy states along the trajectory. The prescribed command Δh1 = −2 mm, Δφ2 = 5°, and Δh3 = −2 mm is recovered with negligible command residual. The final numerical run reports an overall PASS, with all rigid-origami, dual-quaternion, Bennett-limit, homotopy, and conditioning gates satisfied. This confirms that the chosen DQ-parametric homotopy is sufficient for reproducing and validating the prescribed FOLD–TWIST–FOLD motion in the reported configuration.

The sufficiency of the homotopy is not assessed only from the final command values. The rigid Kresling geometry is checked at every continuation state. The maximum relative errors for the ring edges, primary creases, diagonal creases, and triangular panels remain of order 10−10. These values confirm that the numerical trajectory preserves the rigid-origami constraints over the full motion, rather than merely reaching a final pose that matches the commanded height and twist.

The reconstructed configurations are also verified through dual-quaternion diagnostics. The unit-DQ residual remains of order 10−16, and the DQ panel-map error remains below 10−8 mm. These checks confirm that the transformations associated with the rigid panels and cells are represented consistently as rigid-body motions in SE(3). Therefore, the dual-quaternion representation supports the rigid-origami interpretation of the trajectory.

The local Bennett diagnostic is evaluated in the intersecting-axis regime d = 0, which is appropriate for the Kresling crease-axis geometry considered here. Along the reported trajectory, the Bennett distance and the Bennett ratio d/sin(α) remain at numerical zero, and all tested adjacent-axis pairs are classified as intersecting pairs. The physical homotopy residual remains of order 10−10 over the 37 continuation states, confirming that the continuation path remains on the constrained rigid-origami branch.

The reported conditioning value is finite for the adopted constrained reconstruction and is used as a diagnostic of the numerical continuation, not as a global conditioning result for arbitrary Kresling closure problems. The sufficiency claimed here is therefore restricted to the reported validation task: verifying a prescribed rigid Kresling origami FOLD–TWIST–FOLD trajectory for the tested three-cell hexagonal configuration. It is not claimed that the same direct DQ-parametric homotopy replaces Newton, Levenberg–Marquardt, TSVD, Gröbner-basis, or resultant-based techniques for all possible nonlinear closure problems. Such methods may still be required when the mechanism configuration is unknown, when different branches must be discovered, or when broader design spaces are explored.

Accordingly, the DQ-parametric homotopy is sufficient for the reported validation because it constructs a continuous rigid-origami trajectory, preserves the Kresling geometric constraints, maintains SE(3)-consistent dual-quaternion transformations, satisfies the Bennett-limit diagnostic, and verifies the prescribed FOLD–TWIST–FOLD command over the full continuation path.

3.4. Numerical Validation of the DQ-Based Rigid Kresling FOLD–TWIST–FOLD Trajectory

This section reports the numerical validation of the proposed DQ-parametric homotopy formulation for a rigid Kresling origami FOLD–TWIST–FOLD trajectory. The final validation does not rely on a staged Phase A/Phase B Newton procedure or on checking the final pose alone. Instead, the motion is embedded in a constrained DQ-parametric homotopy path, where the Kresling nodal configuration is reconstructed at each continuation state subject to the rigid-origami constraints. The trajectory is then verified through independent diagnostics of motion tracking, ring-edge preservation, primary-crease preservation, diagonal-crease preservation, triangular-panel rigidity, dual-quaternion consistency, Bennett-limit compatibility, homotopy residuals, and conditioning.

3.4.1. Tested configuration and prescribed command

The tested structure is a three-cell hexagonal rigid Kresling origami robot with

$$N_s=6,\nu =3,R=48\mathrm{m}\mathrm{m},H_0=60\mathrm{m}\mathrm{m},{\varphi }_0={30}^{\circ }.$$

$N_s$ is the number of polygon sides, $\nu$ is the number of axial Kresling cells, R is the circumradius of each ring, $H_0$ is the total initial height, ${\varphi }_0$ is the initial cell twist.

The prescribed FOLD–TWIST–FOLD command consists of three sequential actions:

$$\mathsf{\Delta }{h}_1=-2\mathrm{m}\mathrm{m},\mathsf{\Delta }{\varphi }_2=5^{\circ },\mathsf{\Delta }{h}_3=-2\mathrm{m}\mathrm{m}$$

Here, $\mathsf{\Delta }{h}_1$ and $\mathsf{\Delta }{h}_3$represent axial fold commands applied to the first and third cells, while $\mathsf{\Delta }{\varphi }_2$ represents the twist command applied to the middle cell. The use of $\nu =3$corresponds to the three axial Kresling cells involved in the fold–twist–fold sequence.

3.4.2. Homotopy phase path

The motion is embedded in a homotopy path parameterized by

$$\lambda \in \left[\mathrm{0,1}\right].$$

The path starts from the reference configuration and is divided into three physically meaningful phases: first fold, then twist, then fold. At each homotopy state, the Kresling nodal configuration is reconstructed subject to the rigid-origami constraints, including ring-edge preservation, primary-crease preservation, diagonal-crease preservation, and triangular-panel rigidity. The resulting rigid panel and cell transformations are then represented using unit dual quaternions inSE(3).

The final validation uses 37 homotopy states along the trajectory.

The homotopy path is not used merely as a numerical continuation device. It also organizes the prescribed physical motion into the required FOLD–TWIST–FOLD sequence, thereby avoiding the ambiguity of solving only an isolated final pose.

3.4.3. Motion tracking accuracy

The prescribed command is recovered to numerical precision in the reported numerical run. The achieved changes are

$$\mathsf{\Delta }{h}_1=-2.000000\mathrm{m}\mathrm{m},\mathsf{\Delta }{\varphi }_2={5.000000}^{\circ },\mathsf{\Delta }{h}_3=-2.000000\mathrm{m}\mathrm{m}.$$

The corresponding command errors are

$$e_{h_1}=0,e_{{\varphi }_2}=0,e_{h_3}=2.900\bullet 10⁻¹⁰\mathrm{m}\mathrm{m}$$

Thus, the FOLD–TWIST–FOLD motion command is achieved with negligible command residual. This confirms that the DQ-parametric homotopy follows the intended rigid-motion branch rather than collapsing to a trivial zero-motion configuration.

3.4.4. Rigid Kresling Origami Preservation and DQ Consistency

The rigid-origami character of the trajectory is verified independently from the command tracking. In the present validation, the geometric constraints of the Kresling origami model are checked explicitly over the full trajectory. These include ring-edge preservation, primary-crease preservation, diagonal-crease preservation, and triangular-panel rigidity. The maximum relative errors are:

$$e_{ring-edge}^{max}=9.456\times {10}^{-11},$$

$$e_{\mathrm{p}\mathrm{r}\mathrm{i}\mathrm{m}\mathrm{a}\mathrm{r}\mathrm{y}-\mathrm{c}\mathrm{r}\mathrm{e}\mathrm{a}\mathrm{s}\mathrm{e}}^{\mathrm{m}\mathrm{a}\mathrm{x}}=8.846\times {10}^{-11},$$

$$e_{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\mathrm{o}\mathrm{n}\mathrm{a}\mathrm{l}-\mathrm{c}\mathrm{r}\mathrm{e}\mathrm{a}\mathrm{s}\mathrm{e}}^{\mathrm{m}\mathrm{a}\mathrm{x}}=8.846\times {10}^{-11},$$

$$e_{\mathrm{t}\mathrm{r}\mathrm{i}\mathrm{a}\mathrm{n}\mathrm{g}\mathrm{u}\mathrm{l}\mathrm{a}\mathrm{r}-\mathrm{p}\mathrm{a}\mathrm{n}\mathrm{e}\mathrm{l}}^{\mathrm{m}\mathrm{a}\mathrm{x}}=9.456\times {10}^{-11}$$

These values indicate that the rigid Kresling origami constraints are preserved throughout the motion. The maximum physical solver residual is

$$r_{\mathrm{s}\mathrm{o}\mathrm{l}\mathrm{v}\mathrm{e}\mathrm{r}}^{\mathrm{m}\mathrm{a}\mathrm{x}}=9.456\times {10}^{-11}.$$

The unit-dual-quaternion and DQ panel-map consistency checks also remain within the prescribed numerical tolerances:

$$e_{\mathrm{D}\mathrm{Q}-\mathrm{u}\mathrm{n}\mathrm{i}\mathrm{t}}^{\mathrm{m}\mathrm{a}\mathrm{x}}=6.387\times {10}^{-16},e_{\mathrm{D}\mathrm{Q}-\mathrm{m}\mathrm{a}\mathrm{p}}^{\mathrm{m}\mathrm{a}\mathrm{x}}=3.903\times {10}^{-9}\mathrm{m}\mathrm{m}.$$

Therefore, every homotopy state is represented by consistent rigid transformations in $SE\left(3\right)$, while the underlying Kresling geometry preserves the ring edges, primary creases, diagonal creases, and triangular panels over the full FOLD–TWIST–FOLD trajectory.

3.4.5. Bennett-limit and axis diagnostics

The Bennett compatibility diagnostic is evaluated along the full trajectory. Since the local Kresling crease-axis geometry contains intersecting-axis pairs, the relevant interpretation is the intersecting-axis Bennett limit

$$d=0.$$

Under this limiting interpretation, the reported maximum Bennett distance is

$$d_{\mathrm{m}\mathrm{a}\mathrm{x}}=5.676\times {10}^{-15}\mathrm{m}\mathrm{m},$$

and the maximum Bennett ratio is

$${\left({\displaystyle \frac{d}{\sin \alpha }}\right)}_{\mathrm{m}\mathrm{a}\mathrm{x}}=8.189\times {10}^{-15}\mathrm{m}\mathrm{m}.$$

These values are at numerical zero for the scale of the mechanism and confirm that the intersecting-axis Bennett-limit diagnostic is satisfied throughout the trajectory. The axis diagnostics also confirm that all tested adjacent-axis pairs are classified in the intersecting-axis regime: 2664/2664=100%.

The paired-normal diagnostic is reported separately and is not used as a required gate in the final validation. This is because the Kresling application is interpreted through the intersecting-axis Bennett-limit diagnostic rather than as a collection of non-degenerate Bennett 4R linkages.

3.4.6. Homotopy Residuals and Conditioning

The homotopy consistency diagnostics confirm that the continuation path remains coherent over the full trajectory. A total of 37 homotopy states are checked. The maximum physical homotopy residual is

$$r_{\mathrm{h}\mathrm{o}\mathrm{m}\mathrm{o}\mathrm{t}\mathrm{o}\mathrm{p}\mathrm{y}}^{\mathrm{m}\mathrm{a}\mathrm{x}}=9.456\times {10}^{-11}.$$

This residual corresponds to the physical rigid-origami constraints used as validation gates, including ring-edge preservation, primary-crease preservation, diagonal-crease preservation, triangular-panel rigidity, and command tracking. The reported value confirms that the homotopy path remains on the constrained rigid Kresling origami branch throughout the FOLD–TWIST–FOLD sequence.

The conditioning diagnostic gives

$$\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{d}\left(J\right)=4.606\times {10}^6.$$

This value refers to the constrained numerical reconstruction used along the DQ-parametric homotopy path. It is finite for the adopted parameterization and does not prevent the rigid-origami residuals from remaining below the prescribed tolerance. Therefore, the conditioning diagnostic supports the numerical consistency of the reported validation, but it is not claimed as a global conditioning result for arbitrary Kresling closure problems.

3.4.7. Summary of Validation Gates

The final validation satisfies all required gates:

$$\begin{array}{cc}\mathrm{Motion} \mathrm{FOLD}\text{–}\mathrm{TWIST}\text{–}\mathrm{FOLD}\text{:}& \mathrm{PASS},\\ \mathrm{Ring}-\mathrm{edge} \mathrm{rigidity}\text{:}& \mathrm{PASS},\\ \mathrm{Primary}-\mathrm{crease} \mathrm{rigidity}\text{:}& \mathrm{PASS},\\ \mathrm{Diagonal}-\mathrm{crease} \mathrm{rigidity}\text{:}& \mathrm{PASS},\\ \mathrm{Triangular}-\mathrm{panel} \mathrm{rigidity}\text{:}& \mathrm{PASS},\\ \mathrm{Physical} \mathrm{solver} \mathrm{residual}\text{:}& \mathrm{PASS},\\ \mathrm{Unit}-\mathrm{DQ} \mathrm{consistency}\text{:}& \mathrm{PASS},\\ \mathrm{DQ} \mathrm{panel}-\mathrm{map} \mathrm{consistency}\text{:}& \mathrm{PASS},\\ \mathrm{DQ} \mathrm{twist} \mathrm{consistency}\text{:}& \mathrm{PASS},\\ \mathrm{Axis} \mathrm{diagnostics}\text{:}& \mathrm{PASS},\\ \mathrm{Local} \mathrm{Bennett}-\mathrm{limit} \mathrm{compatibility}\text{:}& \mathrm{PASS},\\ \mathrm{Homotopy} \mathrm{path}\text{:}& \mathrm{PASS},\\ \mathrm{Conditioning}\text{:}& \mathrm{PASS}.\end{array}$$

The final overall status is therefore

$$\mathrm{Overall} \mathrm{status}\text{:} \mathrm{PASS}.$$

The paired-normal diagnostic is reported separately and is not imposed as a required gate in the final Kresling validation. The reason is that the numerical model uses the intersecting-axis Bennett limit $d=0$, not the assumption that each Kresling quadrilateral is a non-degenerate Bennett 4R linkage.

Overall, the numerical results validate the proposed DQ-parametric homotopy formulation for the prescribed rigid Kresling origami FOLD–TWIST–FOLD trajectory. The command is achieved to numerical precision, the rigid-origami constraints are preserved over the full trajectory, DQ consistency is maintained within the prescribed numerical tolerances, Bennett-limit compatibility is satisfied throughout the trajectory, and the homotopy path remains physically consistent over all 37 continuation states.

4. Limitations and Future Work

The present work is restricted to the rigid-origami regime. The Kresling panels are treated as rigid bodies, and the crease lines are modeled as fixed revolute axes. Therefore, deformable panels, elastic crease behavior, material compliance, friction, backlash, self-contact, and fabrication tolerances are not included in the current model. These effects may be important in physical prototypes and should be addressed in future experimental studies.

A second limitation concerns the use of Bennett compatibility. The paper does not claim that each Kresling quadrilateral is a non-degenerate Bennett 4R linkage. The non-degenerate Bennett formulation is used as a reference paired-projection framework, while the Kresling application is evaluated in the intersecting-axis limit $d=0$. Thus, Bennett-limit compatibility is used as a local diagnostic for the rigid Kresling crease geometry, not as a proof that arbitrary Kresling quadrilaterals are non-degenerate Bennett 4R linkages.

The paired-normal condition is also treated as a diagnostic assumption rather than as a consequence of orientation closure alone. In the numerical validation, the paired-normal error is reported separately and is not imposed as a required gate. Future work should investigate the geometric conditions under which paired-normal transport compatibility is exactly satisfied, approximately satisfied, or lost in more general Kresling assemblies.

The numerical validation was performed for one representative three-cell hexagonal configuration with $N_s=6$, $\nu =3$, $R=48\hspace{0.17em}\mathrm{m}\mathrm{m}$, $H_0=60\hspace{0.17em}\mathrm{m}\mathrm{m}$, and ${\varphi }_0={30}^{\circ }$. Although this configuration validates the proposed DQ-parametric homotopy for the prescribed rigid Kresling origami FOLD–TWIST–FOLD command, it does not exhaust the full design space. Future studies should extend the validation to broader parameter ranges, different numbers of polygon sides, different numbers of axial Kresling cells, alternative twist commands, larger multi-cell assemblies, and different actuation sequences.

The conditioning value reported in the numerical validation is used as a diagnostic for the constrained reconstruction along the adopted homotopy path. It should not be interpreted as a global conditioning result for arbitrary Kresling closure problems. Future work should study the conditioning of the full nonlinear closure system over broader design spaces and near possible singular configurations.

The present validation is numerical only. Experimental verification using fabricated Kresling prototypes would be necessary to assess robustness under manufacturing tolerances, finite thickness effects, joint clearance, friction, and actuation errors. Such studies would also clarify the extent to which the ideal rigid-origami DQ model predicts the behavior of real origami robotic structures.

Future work will also consider extensions of the DQ framework to nonzero-pitch screw joints, general $nR$loops, deformable or compliant crease models, and optimization-based synthesis of rigid Kresling origami mechanisms. In particular, combining the present DQ-parametric homotopy with future parameter sweeps or design optimization could support systematic selection of $H$, $R$, twist, and cell number for prescribed robotic tasks while preserving $SE\left(3\right)$-consistent rigid motion.

5. Conclusions

This work presented a dual-quaternion formulation for analyzing Bennett-type compatibility in the context of rigid Kresling origami FOLD–TWIST–FOLD structures. Starting from the exact closure of a spatial 4R chain in $SE\left(3\right)$, the Bennett condition was interpreted through projected dual-closure relations under explicit geometric assumptions, including orientation closure, paired-normal transport compatibility, and an angle-symmetric branch.

The formulation does not assert that arbitrary Kresling quadrilaterals are non-degenerate Bennett 4R linkages. Instead, the non-degenerate Bennett setting is used as a reference paired-projection framework, while the Kresling application is treated through the intersecting-axis Bennett-limit case $d=0$. Under this interpretation, Bennett-limit compatibility provides a local diagnostic for rigid Kresling crease geometries.

The proposed framework was validated numerically on a three-cell hexagonal rigid Kresling origami robot with $N_s=6$, $\nu =3$, $R=48\hspace{0.17em}\mathrm{m}\mathrm{m}$, $H_0=60\hspace{0.17em}\mathrm{m}\mathrm{m}$, and ${\varphi }_0={30}^{\circ }$. The prescribed FOLD–TWIST–FOLD command,

$$\mathsf{\Delta }{h}_1=-2\hspace{0.17em}\mathrm{m}\mathrm{m},\mathsf{\Delta }{\varphi }_2=5^{\circ },\mathsf{\Delta }{h}_3=-2\hspace{0.17em}\mathrm{m}\mathrm{m},$$

was achieved to numerical precision. The final validation reported an overall PASS, with all rigid-origami, dual-quaternion, Bennett-limit, homotopy, and conditioning gates satisfied. Ring-edge, primary-crease, diagonal-crease, and triangular-panel rigidity were preserved over the full trajectory, with maximum relative physical residuals of order ${10}^{-10}$. Unit-DQ consistency was maintained at order ${10}^{-16}$, and the DQ panel-map error remained below ${10}^{-8}\hspace{0.17em}\mathrm{m}\mathrm{m}$.

The Bennett-limit diagnostic was satisfied along the full trajectory in the intersecting-axis regime, with the maximum Bennett distance and ratio remaining at numerical zero. The DQ-parametric homotopy checked 37 continuation states and maintained a maximum physical homotopy residual of order ${10}^{-10}$. The reported conditioning diagnostic was finite for the adopted constrained reconstruction and did not prevent the rigid-origami residuals from remaining below the prescribed tolerance.

These results support the use of dual quaternions as a consistent $SE\left(3\right)$-based framework for modeling and validating rigid Kresling origami FOLD–TWIST–FOLD motion. More broadly, the paper provides a DQ-based perspective on Bennett compatibility and shows how its intersecting-axis limiting form can be used in the numerical validation of rigid origami robotic mechanisms. Future extensions should address broader parameter ranges, larger multi-cell assemblies, nonzero-pitch screw joints, compliant crease models, optimization-based synthesis, and experimental validation.