Geometry and Kinematics of Cylindrical Waterbomb Tessellation

Origami tessellations are origami obtained by tiling translational copies of a modular crease pattern. Origami tessellations can be folded from flat sheets of paper and transform into various shapes. Even though origami tessellations are corrugated polyhedral surfaces, they sometimes approximate smooth surfaces of various curvatures, including synclastic and anticlastic surfaces [1]. Such macroscopic surfaces exhibit nontrivial behaviors we cannot expect from the periodicity of its crease pattern and recently attracted much attention of scientists and engineers. For example, Schenk and Guest [2] showed that Miura-ori forms anticlastic surfaces while egg-box pattern forms synclastic surfaces through numerical modeling using bar-and-hinge model. Nasser et al. [3,4] analyzed the approximated surface of Miura-ori and egg-box pattern based on the knowledge of differential geometry.

In this paper, we focus on folded surfaces of waterbomb tessellation, specifically, cylindrical folded surfaces called waterbomb tube. Waterbomb tube is known as origami work “Namako” [5] or more commonly, the name “Magic Ball.” Based on the property of waterbomb tube, various engineering applications are attempted, such as origami-stent graft [6], earthworm-like robot [7], soft gripper [8], and robot wheel [9] using the property that the radius of the tube can vary. The dynamics of the stent graft and the wheel inspired by waterbomb tube has been analyzed [10,11]. Also, the kinematics and structural deformation of waterbomb tube has been studied [12–14]. One of the most interesting phenomena reported on waterbomb tube is that it can form wave-like surfaces [15,16] like in Fig. 1. This is the unique phenomenon not yet observed in the folded surfaces of other tessellations such as Miura-ori and egg-box pattern; however, the theoretical reason why wave-like surfaces arise has been unclear.

Here, our objective is to know “why” waterbomb tube produces wave-like surfaces, i.e., to clarify the mathematics behind the behavior. In this paper, we model waterbomb tube as the sequence of modules and focus on its recurrence behavior instead of simultaneously solving the network of 6R spherical linkages as in [12,13,15,17]. We first provide the kinematic model of each module of waterbomb tube and the relation with adjacent modules to obtain the recurrence relation dominating the folded states of modules based on the rigid origami model with symmetry assumption (Sec. 2). Then, through computation and visualization of the folded states of waterbomb tube using the recurrence relation, we observe that the solutions fall into three types: cylinder solution, wave-like solution, and finite solution (Sec. 3). Furthermore, by applying the stability and bifurcation analyses of the dynamical system of waterbomb, we illustrate how the behavior of the system, i.e., whether waterbomb tube can be wave-like surface or not, changes when the crease-pattern parameters are changed (Sec. 4). Finally, we prove the existence of wave-like solution by applying theorems of discrete dynamical systems (Sec. 5).

First, we introduce parameters of the crease pattern of waterbomb tessellation. We can obtain the crease pattern of waterbomb tessellation by tiling translational copies of a unit module. The entire pattern is controlled by four parameters (Fig. 2): a unit module shown in Fig. 2, left, is controlled by two parameters α, β (α ∈ (0, 90°), β ∈ (0, 180° − α)), and the repeating number of modules in column and row directions is given by two integers N and M (⁠$N∈Z>0,M∈Z>0$⁠).

We consider a waterbomb tube as the rigid folded state of the crease pattern of waterbomb tessellation, such that the left and right sides in Fig. 2 are connected to form a cylindrical form. In our model, we assume N-fold symmetry about an axis and mirror-symmetry about N planes passing through the axis as in the existing research [13] (see Fig. 3). However, unlike the previous research [13], we do not assume mirror symmetry with respect to a plane perpendicular to the axis. Here, the behavior of waterbomb tube is governed by three parameters (α, β, N), because M specifies the finite subset interval of infinitely continuing waterbomb tube.

First, we consider the degrees of freedom (DOF) of each module. Because the kinematic DOF of n-valent vertex in general is n − 3 [18,19], the DOF of waterbomb module (without symmetry) is 3. With the assumed mirror symmetry, the DOF drops by one more degree. Therefore, the module has 6 − 3 − 1 = 2 DOF.

Here, we focus on the module with correct mountain and valley (MV) assignment and pop-up pop-down assignment, where we assume that the center vertex of the module is popped-down, i.e., the sum of fold angles (positive for valley and negative for mountain) are positive, and other vertices are popped-down. The definition of pop-up and pop-down follows that of [20], and it is known that rigid origami cannot continuously transform from popped-up and popped-down state without being completely flat.

To represent the folded states of modules, we consider the isosceles trapezoid formed by two pairs of mirror reflected vertices of the module (Figure 4, left). The parallel bottom and top edges are aligned in the hoop direction, while the hypotenuses are along the longitudinal direction. The length of hypotenuses of the trapezoid is fixed at 2cosα, while the half lengths of bottom and top edges, denoted by x and y, change in (0, sinα) by folding the module. For each given pair of x and y in (0, sinα), there exists eight different states of a module (see Fig. 5). However, only one of the folded state has the correct MV and pop-up/pop-down assignment for the following reason, thus we may safely use x and y as the parameters to uniquely represent the folded state of the module.

To show the uniqueness, first notice that the top and bottom modules illustrated in the same column in Fig. 5 are the mirror reflection of each other through the plane containing the isosceles trapezoid, so the top states are popped-down and the bottom states are popped-up. Therefore, the folded state needs to be one of the four top states. Within the top row we compare the first and second left configurations. Vertex P′ is the reflection of vertex P through the plane defined by vertices O, A, and B, so crease OP is mountain and OP′ is valley, thus we choose vertex P (counterclockwise side of cycle OAB) is selected as the correct side. In a similar manner, vertex Q is chosen over Q′. Thus, the top-left configuration with P and Q is chosen. Finally, we can show that the top-left configuration has indeed correct MV-assignment also for other creases; this can be verified by checking that P, D, and C are on the counterclockwise side of OAB and that P is on the symmetry plane between C and D. Therefore, only the top-left configuration has the correct MV-assignment.

Hereafter, we clarify the mathematics behind wave-like surfaces by interpreting the recurrence relation (3) as a discrete dynamical system. Recurrence relation (3) is categorized as a nonlinear discrete dynamical system in two dimension, as functions f, g are nonlinear functions of two variables x_m, y_m. Also, the system (3) is autonomous system because the function f, g is independent of indices m.

To understand the 2-DOF kinematics of waterbomb tube, we computed the sequences (2) under different pairs of initial values (x₀, y₀) and visualized their configurations as the sequence of points in x, y-plane, i.e., the phase space. The sequence ${(xm,ym)=Fm(x0,y0)|m∈Z}$ is called orbit of (x₀, y₀) and the plot of them is called phase diagram. Here, we fix the crease-pattern parameters (α, β, N, M) = (45°, 34.6154°, 24, 100) and show the phase diagram under this parameter in Fig. 7. Under this parameter, we can observe all possible types of solutions we explain below.

It can be observed that the characteristics of the solution change depending on the given initial values: the solutions can be classified into three types; cylinder solution where it forms a constant radius cylinder, wave-like solution where the corresponding waterbomb tube form wave-like surface, and finite solution where the system fails to obtain a solution after some iterations. Note that the behavior of the system (3) can change under different crease-pattern parameters as we discuss in Sec. 4 (i.e., we cannot observe some types of solutions in different parameters). However, we can still classify solutions under different parameters into these three types.

Cylinder solutions shown in black in Fig. 7 are numerically calculated by solving Eq. (6) under crease-pattern parameter (α, β, N, M) = (45°, 34.6154°, 24, 100), which two solutions exist and each have different radius.

The second case corresponds to third to fifth from the top of waterbomb tube shown in Fig. 7, left. The wave-like folded state corresponds to the sequence of points on the phase space moving in a clockwise direction. We call such solution a wave-like solution. In a wave-like solution, the point continues rotating around the same closed curve without divergence or convergence. Note that the points along the closed curve does not coincide, thus the folded state is not periodic. However, the overall folded shape seems to approximate a smooth periodic wave-like surface. These sequences and corresponding cycles are nested around one of two cylinder solutions (smaller one).

In addition, these plots of nested closed curves seemed to be symmetric about the graph of y_m = x_m. The symmetry of the orbits is the result of the reversibility of the system. Generally, the orbit of reversible systems initiated from x₀ = (x₀, y₀), that is defined as the set ${x=Fm(x0)|m∈Z}$⁠, is invariant under the symmetry G when the orbit has a point which is invariant under G [21]. For this reason, the orbits shown in Fig. 7, which initial terms x₀ are invariant under G, are actually symmetric about the graph of y_m = x_m.

The second waterbomb tube from the top in Fig. 7, left, corresponds to the third type, where its points are plotted just a little outside of the above-mentioned concentric plots. The reason plots stop at the middle is that, at some index m, (x_m, y_m) deviate from the region (0, sin α) × (0, sinα), that is, there is no state of modules corresponding to parameter (x_m, y_m). In other words, finite solution appears in the case that the intersection of three spheres become empty at some step, when computing vertices of modules as shown in Fig. 6. Specifically, in the sequence corresponding to second waterbomb tube in Fig. 7, the numerical value of ninth term (x₈, y₈) is (0.689573, 0.724059), which y₈ is greater than sin α = sin 45° ≈ 0.707107. So, these solution terminates at some point, and only a finite portion of the paper can be folded along this solution.

From this visualization, we can observe a single disk region in the phase space as the set of initial values yielding solutions for any m (the gray region in Fig. 7). The region is the union of all wave-like solutions and the smaller cylinder solution. As the region is the configuration space of the mechanism with m → ∞, there exists a 2-DOF rigid folding motion. The motion can be represented by the “amplitude” and “phase” of wave-like surfaces. As the initial configuration gets closer to the cylinder solution, the amplitude of wave shapes gets smaller. We can change the phase by rotating the initial configuration along the closed curve (Supporting Movie).¹ In the example state shown in Fig. 7, each initial value is taken along x₀ = y₀, so the left-most module forms the “valley” of the wave.

From Fig. 7, we found that the system behaves differently around the different cylinder solutions. In Sec. 4, we perform a stability analysis to classify the symmetric fixed points, i.e., cylinder solutions, by the behavior of the system around them. We also investigate how the cylinder solutions emerge and disappear and how the stability changes when the crease-pattern parameters are changed.

We conjecture that quasi-periodic solutions exist around a neutrally stable cylinder solution, which explains the nested closed solutions representing wave-like solutions. In Sec. 5, we give a proof for the conjecture for a particular crease pattern.

In this section, we investigate how the behavior of the system around cylinder solutions changes when the crease-pattern parameters are changed. For this, we numerically analyze the stability of the cylinder solutions and visualize the changes of the stability using the bifurcation diagram. Figure 8 shows the concept of the analysis in this section. From the bifurcation diagram, we found that there are some critical values of crease-pattern parameters where the stability or the behavior of the system changes. Note that the result of this section is obtained numerically (not symbolically) using Mathematica based on the explicit form of system (3).

In Fig. 8, the top-left figure shows the two different cylinder configurations with larger and smaller radius under (α, β, N) = (45°, 40°, 8) represented by the numerically calculated cylinder solutions (w₁, w₁) and (w₂, w₂), respectively. The linear stability of (w₁, w₁) and (w₂, w₂) are unstable saddle and neutrally stable center, respectively.

A fixed point classified as an unstable saddle point in the linearized system (7) is also an unstable saddle point in the original nonlinear system (see Sec. 5). However, a fixed point classified as a neutrally stable fixed point in the linearized system is not necessarily neutrally stable in the original system. Therefore, the linear stability analysis alone does not guarantee the existence of nested closed solutions. Nevertheless, from the phase diagram in Fig. 8, we found that there are nested closed plots around (w₂, w₂). This leads us to conjecture that if the linearized system (7) has a neutrally stable symmetric fixed point, it is also a neutrally stable symmetric fixed point in the original system, around which there are quasi-periodic solutions. We show the conjecture is correct for a particular crease pattern in Sec. 5.

Next, we analyze how the linear stability of cylinder solutions changes when β changes under fixed α and N. For the visualization of the result, we use a bifurcation diagram. The bottom of Fig. 8 shows the bifurcation diagram for α = 45°, N = 8.

To create bifurcation diagram for given α and N, we move the remaining parameter β in the range (0, 180° − α), and compute the cylinder solutions (w, w) of the crease pattern by solving Eq. (6) at each value. Here, we transform parameter w to $Θ≡arcsin(w/(sinα))∈(0,90∘)$⁠, the half of the dihedral angle formed by the two isosceles triangles (Fig. 8), so that the plot range is independent of α. Then, we perform the linear stability analysis and plot the cylinder solutions in β–Θ planes by solid and dotted line if it is neutrally stable and unstable, respectively. In Fig. 8, the vertical line β = 40° intersects with the bifurcation diagram at two different points, (40°, Θ₁) on the dotted line and (40°, Θ₂) on the solid line, shown in red and blue, respectively. The red and the blue points correspond to the unstable cylinder with a larger radius and the neutrally stable cylinder with a smaller radius, respectively.

Now, we describe how the linear stability of the cylinder solutions depends on the parameters using the bifurcation diagram. First, we observe the bifurcation diagram for fixed N = 8. By varying α, we found the critical value α* ≈ 60.4° at which the appearance of the bifurcation diagram changes. So, we use the representative cases α = 45° < α* and α = 65° > α* for further observing the bifurcation diagram for these cases. We found some critical values of β where the number of cylinder solutions or their linear stability change, i.e., whether waterbomb tube can be wave-like surface changes. We also explain the differences between the two cases. Finally, we show that we can apply the observations of the case of N = 8 for the other N.

Here, we fix N = 8. Under fixed N, we can create the 3D bifurcation diagram by varying the value α in range (0, 90°) and arranging the 2D bifurcation diagrams for each α (right in Fig. 9). From the 3D diagram, we found the critical value α* ≈ 60.4° where the behavior of the bifurcation diagram changes. The behaviors of the bifurcation diagram for α < α* are represented by (α, N) = (45°, N = 8) (blue bifurcation diagram on the top-left in Fig. 9), while the bifurcation diagram for α > α* is represented by (α, N) = (65°, N = 8) (red bifurcation diagram on the bottom-left in Fig. 9).

First, we describe the behavior of the representative case of (α, N) = (45°, 8) described in the blue bifurcation diagram on the top-left in Fig. 9. As we change the remaining variable β, there are some critical value of β denoted by β*₁, β*₂, where the linear stability of a cylinder solution changes. As we increase β from 0, at β = β*₁ ≈ 35.8°, the graph of y = f(w, w) and y = w touches and the fixed point that is saddle and center generated (saddle-node bifurcation). The branch of the fixed point that is saddle extends from the bottom of the graph at β = 45°. Then, at β = β*₂ ≈ 46.1°, saddle-node bifurcation occurs again where the center and the saddle coincide and disappear. Thus, the system (3) has no cylinder solutions under β ∈ (0, β*₁), two cylinder solutions (center and saddle) under β ∈ (β*₁, 45°), three cylinder solutions (center and two saddle) under β ∈ (45°, β*₂), and only one cylinder solution (saddle) under β ∈ (β*₂, 135°). Here, Fig. 10 shows the part of the bifurcation diagram and the phase diagrams for β = 35°, 40°, 45.5°, and 50° belonging to (0, β*₁), (β*₁, 45°), (45°, β*₂), and (β*₂, 135°), respectively. Under β = 40° ∈ (β*₁, 45°) and β = 45.5° ∈ (45°, β₂°), there are nested cyclic plots around the center, which suggests the existence of a wave-like configurations. On the other hand, there are no cyclic plots in the other two diagrams, i.e., all solutions are finite solutions; therefore, waterbomb tube under $α=45∘,N=8,β∈(0,β1*)⋃(β2*,135∘)$ cannot be wave-like surfaces. Thus, whether waterbomb tube can be wave-like surface or not changes at the bifurcation values.

Next, for the case of (α, N) = (65°, 8), the red diagram in Fig. 9 shows that there are three critical values for β: β*₁, β*₂, and β*. As we increase β from 0, saddle-node bifurcation occurs at β = β*₁ ≈ 27.8° in the same way as the α = 45°; however, as β is further increased, the linear stability changes at β = β* ≈ 52.3°, which is not observed in the case of α < α*. Furthermore, the branch of the center fixed point extends from the bottom of the graph at β = 65°. Finally, at β = β*₂ ≈ 65.5°, saddle-node bifurcation occurs again where the saddle and the center coincide and disappear. From the above, in the case of α = 65°, the system (3) has zero, two (center and saddle), two (two saddles), three (two saddles and center), and one (saddle) cylinder solutions when β ∈ (0, β*₁), (β*₁, β*), (β*, 65°), (65°, β*₂), and (β*₂, 115°), respectively. Figure 11 shows the part of the bifurcation diagram and the phase diagrams for β = 25°, 40°, 49°, 57°, 65.1°, and 70° belonging to (0, β*₁), (β*₁, β*), (β*₁, β*), (β*, 65°), (65°, β*₂), and (β*₂, 115°), respectively. In the top-right, middle-left, and bottom-left phase diagrams, there are the plots of wave-like solutions around the center; therefore, waterbomb tube has wave-like configurations. However, the wave-like configurations placed in each diagram self-intersect at some parts in the same way as the configuration in the bottom on the left side of Fig. 11, or this is caused by too small radius as for β = 65.1°, which means these wave-like configurations are not realizable. As for the middle-left diagram, we found the unstable nonsymmetric fixed points (x, y) ≈ (0.698488, 0.375731), (x, y) ≈ (0.375731, 0.698488) satisfying f(x, y) = g(x, y) = (x, y), which is not observed in the case of α = 45°, but the configurations of the module and the tube corresponding to the nonsymmetric fixed points are self-intersecting complicatedly as shown in Fig. 11. The other three diagrams, top-left, middle-right, and bottom-right on the right side of Fig. 11, having no wave-like solutions, which suggests that waterbomb tube cannot be wave-like surfaces if $β∈(0,β1*)⋃(β*,65∘)⋃(β2*,115∘)$⁠. Hence, bifurcation values are closely related to the possible forms of waterbomb tubes as in the case of α = 45°.

For cases other than N = 8, the top, middle, and bottom figure in Fig. 12 shows the 3D bifurcation diagrams for N = 6, N = 30, and N = 100, respectively. Comparing 3D diagrams under the different values of N in Figs. 9 and 12, the shape of the diagrams and the critical values α* depends on the values of N. However, the pattern formed by the dotted and solid lines is the same for the four different values of N; therefore, we can apply the features described for N = 8, i.e, the existence of critical values such as α*, β*₁, β₂*, β* and their relation with the behavior of waterbomb tube for N = 6, 30, and 100.

In this section, we claim that there are nested wave-like solution around the neutrally stable cylinder solution and they are quasi-periodic orbits of the system; and we give the proof of existence for fixed crease-pattern parameters (α, β, N) = (45°, 45°, 6). We also show the system behaves differently around the saddle cylinder solution. For the proof of quasi-periodic orbits, we use the Kolmogorov–Arnold–Moser (KAM) theorem that guarantees the equivalence of the behaviors in the nonlinear and linearized systems under some conditions as described later.

On the other hand, at the smaller cylinder solution (x, y) = (w₂, w₂) with a smaller radius, DF(w₂, w₂) has complex conjugate eigenvalues whose magnitudes are exactly 1. Hence, the type of the cylinder solution (w₂, w₂) is center and linear stability of (w₂, w₂) is neutrally stable. Nevertheless, because (w₂, w₂) is an elliptic fixed point, it is not yet guaranteed that the original nonlinear system (3) has nested closed orbits around (w₂, w₂) as we mentioned in Sec. 4.1.

Now, to prove the existence of the nested closed symmetric quasi-periodic orbits around (w₂, w₂) in the original system, we use the one derived form of KAM theorem. The KAM theorem guarantees the existence of such orbits called KAM curve around a fixed point of the reversible dynamical system if the fixed point is elliptic symmetric fixed point and nonresonant, that is, the eigenvalues are not being roots of unity [21,22]. We already know that the system is reversible and (w₂, w₂) is an elliptic symmetric fixed point. Also, (w₂, w₂) is nonresonant because of the following reason. The eigenvalues of DF(w₂, w₂) are the roots of the irreducible polynomial (15). Then, some coefficients of the minimal polynomial, obtained as the monic form of (15), are not integers (e.g., coefficients for seventh-order term is 292/13), so the eigenvalues of (w₂, w₂) are not algebraic integer, which implies that the eigenvalues of (w₂, w₂) are not the roots of unity. Thus, we can apply KAM theorem to the system at (w₂, w₂), and the existence of wave-like solution around the point is proved.

In this paper, we revealed a part of the mathematical structure behind wave-like surfaces of waterbomb tube. First, we have described the kinematic model representing the configuration of each module of waterbomb tube and its entire shape to derive recurrence relation of two variables dominating its kinematics. Then, based on the observation of the plots on the phase space corresponding to the folded states, we classified solutions into three types to identify the wave-like solutions surrounding one of two cylinder solutions. We applied the knowledge of discrete dynamical systems to the recurrence relation, and investigated how the stability of the symmetric fixed point of the system, i.e., the behavior around cylinder solutions, changes when we change the crease-pattern parameters. Then, we observed that around the neutrally stable cylinder solution, there exists nested wave-like solutions, and found critical values where the behavior of the system changes. Finally, we gave the proof of the existence of quasi-periodic solutions, i.e., wave-like solutions, around the neutrally stable cylinder solution.

Note that we gave proof of the existence of quasi-periodic solutions only for fixed crease-pattern parameters; characterizing the existence of such solutions is a future work of this paper. Also, the physical interpretation, i.e., the relation between the mathematical properties of solutions and the mechanical properties are still left to be investigated.

More generally, our proposed approach based on the geometry of a module and applying the knowledge of dynamical system to the recurrence relation can be useful to explain behaviors of various origami tessellations. The analysis of tessellations consisting of modules with higher or lower DOF is of our interest.

This work is supported by JST PRESTO (Grant No. JPMJPR1927).

Here, we derive Eq. (6) imposing cylinder constraint. Assuming that values of crease-pattern parameters are given, we represent some distance of vertices or angle using (α, β, N, M) based on the preservation of edge length.

Here, we derive the function f and g assuming that values of crease-pattern parameters and the state of modules belonging to mth hoop, i.e., x_m and y_m, are given. First, we represent some quantities such as distance between vertices or angles by crease-pattern parameters and x_m and y_m. Next, we derive the function f and g by introducing the coordinate system and representing x_m+1 and y_m+1 by the coordinates of some vertices of waterbomb tube.

First, we represent some quantities using (α, β, N, M) and (x_m, y_m) based on the preservation of edge length. The notation of vertices is defined in Fig. 16, and for simplicity, we omit subscripts of x_m and y_m.

Next, we derive the coordinates of vertex $Bm+1,12$⁠. Vertex $Bm+1,12$ is located at a distance 1, 1, 2 cos α from vertices $Om+1,0,Om+1,1,Cm,02$⁠, respectively. Here, since the vertices O_m+1,0, O_m+1,1 are the mirror-reflection with respect to XZ-plane, the vertex $Bm+1,12$ is contained in the intersection of a circle of radius $1−(Om+1,1y)2$ centered at the vertex Q on the XZ-plane and a circle of radius 2cosα centered at the vertex $Cm,02$⁠. We can obtain the intersection of the circles by Eq. (A4) where $x1=[Cm,02x,Cm,02y,Cm,02z]T,R1=2cosα,x2=[Qx,Qy,Qz]T,R2=$$1−(Om+1,1y)2,n=[0,1,0]T$⁠. In this case, the positive sign is consistent with the prescribed MV-assignment of $Cm,02Bm+1,12$⁠. Therefore, we obtain the coordinates of vertex $Bm+12$⁠. Since $ym+1=[Bm+12x,Bm+12y,Bm+12z]T⋅[0,cosπN,sinπN]T$⁠, the function g is obtained.

There are no conflicts of interest.

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.