An experimental study and model determination of the mechanical stiffness of paper folds

Over the past few decades folding paper has extended beyond the origami deployable applications to reach the engineering field. Nevertheless, mechanical information about paper behavior is still lacking, especially during folding/unfolding. This article proposes an approach to characterize the paper fold behavior in order to extract the material data that will be needed for the simulation of folding and to go a step further the single kinematics of origami mechanisms. The model developed herein from simple experiments for the fold behavior relies on a macroscopic local hinge with a non-linear torsional spring. Though validated with only straight folds, the model is still applicable in the case of curved folds thanks to the locality principle of the mechanical behavior. The influence of both the folding angle and the fold length are extracted automatically from a set of experimental values exhibiting a deterministic behavior and a variability due to the folding process. The goal is also to propose a methodology that may extend the simple case of the paper crease, or even the case of thin material sheets, and may be adapted to other identification problems. This is a postprint of an article that was published in its final form in Journal of Mechanical Design, volume 138, issue 4, ASME Digital Collection URL: https://asmedigitalcollection.asme.org/ mechanicaldesign/article-abstract/138/4/ 041401/474813, ASME © 2016. This manuscript version is made available under the CC-BY distribution license https://creativecommons.org/licenses/by/4.0/ doi: 10.1115/1.4032629


Nomenclature α
Folding angle Residual variability component σ Singular value L Fold length n α Number of discretization points for α n L Number of different possible lengths t Folding torque V Variance

Introduction
Research on mechanisms made from assemblies of thin plates have been going on for a few decades now. The mechanical studies were mainly focused on deployable panels for spatial devices [5,35] or adjustable architecture [4,14], but also in other fields such as medical prosthesis (stents in [34,24]) and biomimicry studies (confinement in buds [33,23]). Many of the studied systems provide a geometric pattern of creases simply repeated in the plane, in a way that could be pursued to infinity in every direction of this plane, and called tessellation. They were discovered through the origami art and experimentally studied with paper in the 1960s and 1970s [27,15]. Their kinematic and dynamic abilities were explored mainly under the pseudo-rigidity assumption [32] considering the rigid and creases as the concentration zone of strain and movement. It simplifies modeling but does not fully describe the degrees of freedom of tessellation structures. It is now known that the movement of those depends not only on the crease pattern but also on the flexibility of faces [8,25,31]. Concerning engineering design, paper models are also efficient for prototyping foldable structures made with other materials. The difficulty is to proceed to a change in scale and material. Once a paper model has been tested experimentally, similitude is a useful tool to predict the real-scale and real material engineering structure design. Nevertheless the behavior of both materials (hence the paper fold behavior) have to be established. If simulation is used, the numerical models also need to be fed by a fold behavior model. This is the case when designing paper-made structures, for instance for packaging with lockings and pop-ups that take advantage of the flexibility of the paper to open the structure (spring joints and their elastic energy are used as distributed unilateral motors). The same situation occurs for deployable light-weighted spatial structures when a single opening operation is required: though they are usually not made with paper, some paper prototypes are useful in early design stages.
Paper is mainly designed in order to resist to industrial printing processes [22], endure moisture [20] or cycling folding that cause tearing [21]. . . These criteria lack the description of paper as a thin plate and the fold as a mechanical joint. This article proposes a methodology to characterize paper, and more specifically the folds.
First, paper is considered as a classic orthotropic material in a thin fiber composite plate. Second, the crease behavior with a simple macroscopic model will be considered. Nevertheless, this article is only based on an application to a special case: the paper used for these experiments is a 120 g/m 2 ECF woodfree pulp uncoated paper: the Arjowiggins Conqueror CX22™.

Experiments and results
Two different behaviors are studied: (i) the tensile test of paper faces (in-plane traction) which is a classical mechanical test for industrial paper and (ii) the behavior of folds at a macroscopic scale, i.e. the torque-angle relationship for a monotonic loading. Smooth bending of paper that can be modeled with a thin plate behavior and identified with dedicated tests, is not taken into account herein. Nevertheless some classical in-plane behaviors are of interest since they can influence the fold behavior. They are therefore tested and reported in the following, while focus remains on the fold behavior. The creasing process itself is not of main concern, though it may influence the subsequent fold stiffness.

Preliminary data: grammage and thickness
Classical parameters for industrial paper are the thickness and the surfacic density. Mean thickness e = 129 µm has been obtained with sampling points and a precision ∆e ≤ 5 µm on 15 specimens whose mean area a = 100 cm 2 is measured with a precision ∆a ≤ 2 cm 2 . Mass m = 11.62 g is measured with a precision ∆m ≤ 0.01 g, at an ambient temperature of 25°C and a relative humidity (RH) of 40 %. This leads to relative accuracies of ∆e/e ≤ 3.9 % for the thickness and ∆ρ s /ρ s ≤ 2.1 % for the surfacic density ρ s = 116.2 g/m 2 . This is in accordance with the statistical variations in the measurements (∆e/e = 4.8 % and ∆ρ s /ρ s = 1.7 %) which are therefore not significant.

In-plane tensile tests
Paper faces are modeled as thin plates. The behavior characterization assumes that these faces remain plane and that only tension solicitations are applied (compression would lead to buckling). Due to the presence of a microstructure (roughly, a weakly entangled and in-plane layered fiber arrangement component, within a softer matrix) and a specific manufacturing process, the behavior is expected to be anisotropic. Indeed, the fiber orientation is not uniformly distributed in the paper plane [30]. Two paper orientations are tested: (i) the so-called machine direction (MD) where the stress is oriented along the main fiber direction, and (ii) the cross-machine direction (CD), perpendicular to the main fiber direction.
Experimental device. A Lloyd LF Plus machine equipped with a 1 kN load cell is used. Tests are performed with a 1 mm/min velocity. The samples are all 30 mm wide and 125 mm long (between jaws). The force is measured with the load cell, while a camera is used to determine the local strain field via digital image correlation (using the Vic2D analysis software) [28,2,9]. To do so, a black paint speckle pattern is first applied on one face of the specimen. A comparison between virgin and speckled sample behaviors allows to assess that the sprayed paint does not influence the mechanical behavior of the paper. For each paper direction, 8 samples are tested, with a RH between 40 % and 45 %, and a temperature between 22°C and 23.5°C.
Experimental results. Figure 1 reports the obtained traction curves, with the average behavior and bars corresponding to the minimum and maximum obtained values, using 8 specimens for each traction direction. The results confirm the anisotropy of the paper. Data dispersion is quite small (the ratio of maximum shift in amplitude to the maximum stress is less than 2.6 %) allowing to conclude to a good reproducibility of the experiment.

Crease behavior
The crease behavior is studied at a macroscopic point of view: it is modeled as a hinge joint with a torsional spring (possibly non linear). The movement is expected to be sufficiently slow to discard viscous effect and the experiment is based on static equilibrium. The fold process is indeed complex at the microscopic scale [13], but is not under the scope of this article. The aim is to observe the rigid movement of two faces linked by a simple straight fold. If the identification of a local fold behavior is performed (i.e. independently of its length and orientation), the model would be applicable to curved folds for which the bending of the faces should also be taken into account [19,6,7]. The criteria for designing the testing device are the following: a possible large range for the folding angle, a sufficient control precision for small loads (to test somehow low grammage papers), and a uniform loading of the fold along its length. A load control with calibrated weights, with a simple device, and not the standard creasing test rigs, is used. Protocol. A sample with a single crease is prepared, splitting the paper in two identical parts and boundary conditions are: a known load prescribed by increasing weights and a low stiffness hinge (in order to neglect its stiffness with respect to the fold stiffness, so to consider it as a perfect hinge), see figure 2. The opening angle of the sample for each weight is measured with a protractor allowing a precision ∆α ≤ 1°.
As a quantity of interest, the torque applied on the crease is estimated with the assumption that the boundary hinge is perfect, and considering the results presented in [25], a face length assumed to limit the influence of the face bending is selected: samples are 100 mm long, folded in their middle, while their widths may vary, in order to test the influence of the length of the crease. The load is applied at L p = 50 mm of the crease, on the free face and on its symmetry axis in order to impact as homogeneously as possible the whole length of the crease. When neglecting the face bending, figure 2, the torque is given with t = mgL p cos(α/2) where m is the loading weight, gravity is g, α is the opening angle. A dispersion analysis provides the precision estimate for the torque: ∆t/t ≤ ∆m/m + ∆L p /L p + 1 2 (tan α)∆α; depending on the angle value, this relative precision is between 2.2 % and 6 %. For the present measurement campaign, three samples are tested for each of the nine categories: the fold orientation could be longitudinal (parallel to MD), transverse (perpendicular to MD) or at 45°, and the fold length L could be 30 mm, 40 mm or 50 mm.
Concerning the creasing process, the experimental protocol is the following: each specimen if completely folded in half, with the thumbnail, from the fold center towards each fold extremity, completed by a sweep along the whole length from one extremity to the other. Then, the paper is completely unfolded, using a 5 kg roller, and folded again using the same roller. Finally, the fold is progressively opened within the experimental device, using prescribed marked weights and measuring the opening angle α. The experiments where performed with an ambient temperature between 20.5°C and 21.5°C, and a RH between 38 % and 47 %.
Experimental results. Figure 3 reports the obtained torques as functions of opening angle, fold orientation and length. It exhibits a large dispersion, much larger that the measurement precision that can be neglected (see e.g. the convex hull span of figure 3a where several specimens of same length and orientation are considered). This dispersion is considered herein as an intrinsic part of the model, and should be identified together with the deterministic model itself. Figure 3b depicts the influence of the fold length: similar dispersions and a trend to increase the torque with the length that can be guessed. Finally, figure 3c depicts the influence of the fold orientation, that is clearly not as important as for the traction curves.

Advanced model identification
A modeling approach may consist in an a priori choice of the form of the dependence of the results on the inputs. Then a data fitting provides the value of the model parameters and the fitting error. Another approach relies on the latent model discovery from the raw data, inspired by data mining. It is used herein, with a small amount of a somehow structured data, which simplifies the approach: the principal component analysis (PCA) or the singular value decomposition (SVD) [17], can provide a numerically identified model.
For the present study, the folding torque t is assumed to depend both on the folding angle α and orientation, and on the fold length L. Additionally, for a given fold orientation, we assume that it can be modeled as a separated function of the folding angle and the fold length. Therefore it can be expressed as t = k(α)f (L) without a priori knowledge on the two functions k and f . Since variability is intended to be non negligible, it is also a parameter of the model, and ε will denote its associated random variable, so that t = k(α)f (L) + ε. Finally, due to the small discrepancies on measurements when compared to the variability, the measurement errors are neglected in all of the following.
Since each specimen i (i ∈ {1, . . . , n}) takes its fold length L i in a single list of n L values, the data associated to this length is much structured. On the other hand, since the tests are controlled by prescribing weights, each specimen i has its own list of n i angles and torques which are therefore non-structured data. As a consequence, for a given orientation of the straight fold, the data is composed for each specimen i of: (i) the length of the fold L i , (ii) a list of opening angles (α ij ) j=1...ni , and (iii) a list of associated torques (t ij ) j=1...ni . If the variability is assumed not to be present for a specific single specimen (but only between different specimens), and if the torque is related to the angle with a single smooth dependence, an interpolation of this angle-torque experimental curve from sampled points is meaningful for each specimen independently.

Restructuring data: interpolating and averaging
To get a single list of discretized n α values of angles (a priori independent of the numbers n i ), the angle range is merely discretized in n α − 1 elements as for a 1D regular finite element (FE) mesh.
To get the corresponding values for the torque, for each set independently, the following framework is used: (i) interpolate the available data (angle-torque curve) specimen per specimen in the current set (for instance by linear interpolation); (ii) average the obtained torque values for a same angle, for the specimens of the current set. It therefore ends up with a single list of angles and as many lists of torques as there are lengths. These data can then be collected into a rectangular snapshot matrix C whose size is n α × n L . The averaging procedure filters out a priori the variability and depends on the modeling parameter n α , i.e. the refinement in the discretization for the function k(α), provided by the user.
All-in-all, the snapshot matrix is a sampling of the functional model t(α, L), with a single list of angles α and lengths L. A standard SVD may therefore be used for assessing an emerging deterministic model: A thin-SVD provides n p = min(n α , n L ) modes in two orthogonal matrices U of size n α × n p and V of size n L × n p storing the tabulated functions u p and v p in their columns and a diagonal matrix Σ storing the positive singular values σ p , p = 1, ..., n p ordered by decreasing amplitudes. In practice, the experimental data may not cover all the range of discretized angles α within a single set, so that some entries in matrix C may lack. This is due to the load control of the experiments with a fixed set of weights leading to a range of values for the angle which is specimen dependent, see Figure 3(a), while the range for the discretized function k(α) is fixed and unique. To address this issue, the optimality theorem of Eckart-Young [10] is considered to recast the SVD analysis in a minimization procedure, including a weight for each entry of C in a boolean matrix W (a mask): a missing value in C corresponds to a null weight in W . The proposed approach is very similar to the so-called iterative gappy POD [11,26]. A rank-one approximation of matrix C is of the form uv T , with column vectors u of size n α and v of size n L obtained with the minimization problem (well posed as soon as W does not exhibit any null full row nor column): The underdeterminacy in the product uv T is reduced by prescribing an arbitrary normalization, for instance v T v = 1. The corresponding generalized singular value is obtained with σ = √ u T u. This procedure can then be applied recursively, on the matrix C − uv T , to find a second rank correction, and so on, to provide the successive modes of this generalized SVD with weighting. With this approach, the snapshot matrix contains only deterministic values (the raw data has been filtered with interpolation and averaging) and the rank-one corrections are also all deterministic functions of angle and length.
For the transverse orientation of the fold, n = 18 specimens are tested, 6 for each of the n L = 3 values of the fold length. There were n i = 13 measured points per specimen. n α is a user parameter, related to the precision of the discretization of the function of the angle. The obtained singular values are reported in table 1. The small ratio σ 2 /σ 1 (less than 2.6%) validates the assumption on the separation of variables α and L in the model, so that a rankone approximation is sensible. The various modes, scaled by the associated singular value are plotted in figure 4 for n α = 5 and for n α = 31 (actually, √ n α u/ √ σ and √ n L v √ σ are depicted). If n α = 2, a linear evolution is provided and when n α is too large, some oscillations may appear.
The assumption that the torque is proportional to the length L can be checked by looking at the mode 1 for v 1 (L) that is close to a linear function of L. For instance, with n α = 5, the RMS error e RMS for the mode 1, once a linear regression is performed, is less than 3.4%. A non-linear dependence of the torque on α leads to model a non-constant stiffness with respect to the opening angle.
The first modes u 1 (α) obtained for different values of n α are compared on figure 5(a), where the functions v 1 (L) have been normalized to best fit the identity function.

Direct analysis of raw data
To avoid the a priori interpolation and averaging pre-processings, the identification of both the variability and the deterministic parts in a single shot procedure is now proposed. A numerical determination of the deterministic model is performed by discretizing the corresponding modeling function, as in the previous SVD analysis. A specific norm is used (the residual least square as for a classical regression) for the random part, while a Frobenius norm is still used for the other variables. The angles and torques are now collected in n L sets of specimens with the same length, leading to a wider list of angles and torques, containing the variability component: (α ji , C ji ), i = 1 . . . m j for the set of specimens of the same length L j (j = 1 . . . n L ). If a separated variable approximation for the torque is used, t(α, L) ≈ k(α)f (L), the residuals are ε ji = C ji − N (α ji )κf j , N (α) being the FE-like shape function vector (once the space of α is discretized), and κ the column vector of nodal values of the discretized function k. The least square error is, for the discrete set of values, e 2 j = mj i=1 ε 2 ji and e 2 = n L j=1 e 2 j . The problem of the automatic determination of the discretized functions k and f is therefore min κ,fj e 2 . This leads to a small coupled system of equations: Note also that an underdetermination in this system still requires an additional arbitrary normalization equation, for instance n L j=1 f 2 j = 1. k(α) and f (L) are deterministic functions, since the variability could be defined as lying in the residual ε = t(α, L) − k(α)f (L). This assumption could be checked if, once the previous minimization problem is solved (for instance with a fixed point, or an alternating projection method), the procedure is iterated for determining the next order correction (as provides the SVD): t(α, L) − k(α)f (L) ≈k(α)f (L), etc. The validation lies in the fact that the second order deterministic corrections should be small when compared to the first order ones. The SVD usually compares the singular value amplitudes of the successive corrections; these singular values correspond here to σ = Algorithm 1 Fixed point algorithm to find the best rank one approximation Inputs: Data set (α ji , C ji ), L j and discretization of α {Compute constant elementary contributions} Return: κ and f j Euler equations corresponding to the minimization for a second couple of correction,κ andf j , lead to and, with the minimization for the previous correction: If the various α ji correspond exactly to discretization points of α, then the matrices M j are identity, and the previous orthogonality conditions read: [ j f jfj ](κ Tκ ) = (κ Tκ )[ jf jfj ] = 0 leading naturally to the classical SVD orthogonality conditions: j f jfj = 0 or κ Tκ = 0, and two possible normalizations for the two unknowns: jf 2 j = 1 orκ Tκ = 1. The present case is a generalization, with a cross-unknowns orthogonality (1).
Hierarchical modeling. When using the rank-one approximation, once the parameters κ and f j are identified, the residual ε ji = C ji − N (α ji )κf j contains the random component. Therefore, the partition between a deterministic part k(α)f (L) and a random part ε, is a modeling choice driven by the choice of the  Table 2: Evolution of the first singular values for the longitudinal orientation, vs the number of elements chosen to discretize the angle function. e RMS is the error in the linearity for the function f (L); V is the variance with respect to the deterministic function. Fold orientation is transverse. Dimensional quantities are in Nmm.
discretization of α. For instance, using only two nodes for the discretization leads to a classical linear regression. Therefore the number of nodes can be increased, in order to increase the size of the deterministic description space. Nevertheless, n α should not be too large, first to avoid singular matrices M j (otherwise, a Tikhonov-like regularization should be added [16]), and second, to avoid oscillations in function k(α) denoting that the variability part is not correctly filtered, figure 5(b). One can note the averaging property of this approach: the successive solutions k(α) present each an additional correction to the previous ones, while preserving the same generalized averages. In the following, the representation of the function k(α) with n α = 5 is selected.
Variability issue. For the random part, the averageε of the random variable ε is zero due to the presence of a constant function in the minimization test function space. The variance can then be estimated with: j=1 m j being the total number of torque values and e is the function to be minimized. The deterministic model is therefore identified thanks to the minimization of the variance. For a rankone approximation: leading to the maximization of the predominance of deterministic mode number one. This identification is also expected to be effective for behaviors with smaller variability. The comparison between the raw data and the model is depicted in figure 6. The variability is assessed by the cumulative density function (cdf) that can be drawn from the data contained in ε ji . Once this experimental cdf is available, figure 7(a), the limits of the model in terms of level of confidence can be computed as well: with a level of 90%, bounds on the model are the enveloping surfaces on figure 6.
The other fold orientations can be dealt with in the same way. Figure 8 reports the evolution of the functions k(α) for n α = 5, and for the three tested orientations. Though the traction tests   exhibited differences depending on the orientation, this is much reduced when considering the fold stiffness. One can therefore model all the orientations in a single shot procedure. Gathering all the data leads to a set of 43 specimens and m = 611 experimental points. Analyzing the whole set leads to a best fit linear error for f (L) of 2% and a ratio of singular values σ 2 /σ 1 = 3.1%, a standard deviation of √ V = 0.289 Nmm, and the cdf of Figure 7(b). Indeed, with a single identification of the whole experimental set, there is a significant increase in the variability.

Discussion
Anisotropy. The in-plane tensile tests show that the machine direction resists up to a 50 % higher stress, but the cross-machine direction can elongate twice as much, assessing the level of anisotropy for the paper behavior. The fold behavior does not exhibit such a dependency on its orientation. The fiber orientation is not uniformly distributed in the paper plane contributing to the stiffness anisotropy observed in the in-plane tensile tests [30]. The case of the fold is different due to its complex but localized behavior [13,3]: the involved physical phenomena (large strains, damage...) when creasing occurs at a small scale (the paper thickness) are depending on the local paper constitution close to the fold. The fold behavior is driven by the cross-section along the fold line, with a complex process at this scale that may allow local fiber rotations. Consequently the fraction of the crossed fibers is the first order influent parameter, leading to a similar behavior independently of the fold orientation. Indeed, image analysis of transverse cuts along different directions shows almost identical covariances [29]. This weak dependence on the fold orientation can be further noticed in several literature results, such as the figure 31 in [18] for experiments on paperboard. The previous justification is nevertheless to take with care; indeed, the fold-affected zone is not so thin (presently it is 2.3 times the thickness, and 10 times the fiber section diameter), and further insight is needed to identify the true mechanism of orientation independence.
Variability. The reproducibility of the tensile test illustrates the fact that the paper making process is controlled up to a high precision. The fold behavior nevertheless exhibits a large dispersion which is a consequence of the material and of the creasing process which is difficult to control. With a given folding protocol, the variability is therefore a part of the model itself and has to be identified as well.
Model identification. The fact that torque is proportional to length allows to confirm that the test rig produced a quite uniform repartition of torque along the fold. The identification is therefore suited to obtain a local stiffness of the fold, allowing the same model to be used for a curved crease. The one-shot proposed identification leads to the automatic determination of characteristic functions (of the opening angle and of the fold length). The user has to specify the number of degrees of freedom needed to represent these evolutions; a hierarchical approach is herein suggested: beginning with a linear interpolation, the discretization could be successively refined up to the point where the function evolution is stable, but before exhibiting oscillations. This cut-off frequency is the transition between the deterministic part and the random part of the model. Experimental devices. Black speckle pattern on white paper gives optimal conditions to realize digital image correlation (contrast, light. . . ). Nevertheless, tests need to be performed with great precision: load repartition, sample installation and clamping. . . The same technique could also be helpful to observe the deformation of an entire paper model structure. The care devoted to boundary conditions has to be particularly important for the fold characterization in order to produce the fold intrinsic behavior. Moreover, with the different movements of an origami tessellation, paper folds are used on a wider range than the tested one, typically from 0°°(completely closed) to approximately 160°. Finally, the folding/unfolding process may be repeated in the tessellation use, leading to an evolution of the folding torque that can be due to an evolution of the microstructure (damage...). In such a case, a model with some internal variables to track the history would be useful.

Conclusions
The proposed methodology was first intended to develop a phenomenological macroscopic model of the paper fold, in order to use it on an origami tessellation structure. It could be useful for other thin material characterizations for structural applications, e.g. polymer films [1], or other thin films [12]. Nevertheless, it is herein developed for a monotonic unfolding path, on a limited angle range, without consideration of material fatigue nor damage. The model identification relies on a deterministic model and a variability component, both to be identified. A principal component analysis taking into account both of these aspects has been proposed and studied. Without many assumptions on the parameter influence on the model, it allows to determine the evolution of the output of interest on the whole range of parameter variations (here: the linear dependency with the length, and the non-linear dependency with the angle). Fold orientation appears to have only a second order influence, that could have been determined by an experimental design as well.