Emilio Turco

Emilio Barchiesi

Francesco Dell'isola

Nonlinear dynamics of origami metamaterials: energetic discrete approach accounting for bending and in-plane deformation of facets

Keywords: Mathematics Subject Classification. 74A05, 74A60, 74B20, 74G65, 74H45, 74K15 Nonlinear dynamics, Origami metamaterials, Time-stepping analysis

come

Nonlinear dynamics of origami metamaterials: energetic discrete approach accounting for bending and in-plane deformation of facets

Introduction

In this paper, structural elements having an origami-type microstructure and micro-kinematics are considered. The model used in order to get quantitative predictions about their nonlinear dynamical behaviour is finite dimensional: the origami microstructure configuration is characterized by means of the vertices positions of the constituent triangular discrete elements. All introduced triangular elements are suitably interconnected via holonomic constraints to constitute the considered metamaterial. A Lagrangian energy function is introduced to establish the evolution equations of the considered origami structures and the used time integration scheme allows us to explore the promising dynamical behaviour of an origami beam (for a survey of existing models for exploring the dynamical behaviour of metamaterial see [START_REF] Del Vescovo | Dynamic problems for metamaterials: review of existing models and ideas for further research[END_REF]).

An ancient technique revived for modern technological applications

The Japanese word origami originally stands for the art of folding sheets of paper in order to produce three-dimensional shapes. The idea behind this art has been the source of several scientific proposals in the field of metamaterials [START_REF] Barchiesi | A review of some selected examples of mechanical and acoustic metamaterials[END_REF][START_REF] Barchiesi | Mechanical metamaterials: a state of the art[END_REF]. Specifically, during the 80's, it has been the source of a proposal [START_REF] Miura | Method of Packaging and Deployment of Large Membranes in Space[END_REF] on the design of periodic metamaterials whose unit cell can be folded, making use of typical origami techniques. Such a proposal has been then further developed in the papers [START_REF] Miura | Map fold a la Miura style, its physical characteristics and application to the space science[END_REF][START_REF] Miura | A note on intrinsic geometry of origami[END_REF][START_REF] Kawasaki | On the relation between mountain-creases and valley-creases of a flat origami[END_REF].

Initially, the most promising application for such a kind of metamaterials was its use-in conjunction with solar panels-for space missions. More recently, several applications have been thought of for origami metamaterials in the biomedical field. Indeed, as an instance, they could be employed for building medical stents, see [START_REF] Silverberg | Using origami design principles to fold reprogrammable mechanical metamaterials[END_REF]. To this end, cylinders possessing the so-called Kresling pattern seem a very attractive solution [START_REF] Kresling | Natural twist buckling in shells: from the hawkmoth's bellows to the deployable Kresling-pattern and cylindrical Miura-ori[END_REF], as they exhibit spontaneous buckling under torsional load. More generally, origami metamaterials can be used for building devices with programmable shape [START_REF] Dudte | Programming curvature using origami tessellations[END_REF].

The possibilities opened by modern 3D-printing technology

Albeit the first papers dealing with origami metamaterials date back to the 80's, owing to the lack of appropriate manufacturing technologies, the realization of these metamaterials has been long considered infeasible. Thanks to new developments in additive manufacturing like 3D-printing, which can reach extremely high precision and allows for creating complex assemblies without the need of performing an actual a posteriori assembling, origami-based materials are currently manufactured and studied. A wide scientific literature has been attempting at studying them from the purely geometric viewpoint. Indeed, folded sheets can only approximate surfaces, and moreover, passing from a flat to a curved surface is not always possible. It has to be remarked [START_REF] Turco | Enhanced Piola-Hencky discrete models for pantographic sheets with pivots without deformation energy: numerics and experiments[END_REF][START_REF] Golaszewski | Metamaterials with relative displacements in their microstructure: technological challenges in 3D printing, experiments and numerical predictions[END_REF][START_REF] Spagnuolo | Phenomenological aspects of quasi-perfect pivots in metallic pantographic structures[END_REF] that it has been possible to print 3D microstructures in which relative rotations and displacements of their different parts are made possible. This technological capacity clearly opened the path towards the construction of specimens of metamaterials, which, before, could be only conceived as a logical possibility.

Models for origami metamaterials

An increasingly wide literature is interested in describing the mechanical behaviour of origami metamaterials, aiming at characterizing their mechanical properties, mostly making use of (semi-)discrete and/or generalized continuum modelling [START_REF] Dell'isola | At the origins and in the vanguard of peridynamics, non-local and highergradient continuum mechanics: an underestimated and still topical contribution of Gabrio Piola[END_REF][START_REF] Abali | Theory and computation of higher gradient elasticity theories based on action principles[END_REF][START_REF] Abali | Strain gradient elasticity with geometric nonlinearities and its computational evaluation[END_REF], this last being required by the occurrence of (quasi-)zero-energy deformation modes [START_REF] Rivlin | Plane strain of a net formed by inextensible cords[END_REF], as well as local instabilities and size effects [START_REF] Wang | Energy-minimizing deformations of elastic sheets with bending stiffness[END_REF][START_REF] Steigmann | Equilibrium of elastic nets[END_REF].

It is worth to note that, addressing the so-called synthesis problem [START_REF] Dell'isola | Synthesis of fibrous complex structures: designing microstructure to deliver targeted macroscale response[END_REF][START_REF] Milton | Which elasticity tensors are realizable?[END_REF], having the goal of developing a deployable metamaterial-namely a three-dimensional material that, exploiting mechanisms, can be compressed into surfaces or a two-dimensional material that can be compressed into curves-many solutions have been proposed in the literature, like pantographic metamaterials [START_REF] Dell'isola | Pantographic metamaterials: an example of mathematically driven design and of its technological challenges[END_REF][START_REF] Cuomo | First versus second gradient energies for planar sheets with two families of inextensible fibres: investigation on deformation boundary layers, discontinuities and geometrical instabilities[END_REF][START_REF] Dell'isola | Advances in pantographic structures: design, manufacturing, models, experiments and image analyses[END_REF][START_REF] Barchiesi | Wave dispersion in non-linear pantographic beams[END_REF][START_REF] Barchiesi | Pantographic beam: a complete second gradient 1D-continuum in plane[END_REF], which exploit a scissor-like mechanism.

It has to be also underlined that while the notion of convergence can be still used [START_REF] Allaire | Homogenization and two-scale convergence[END_REF], the standard homogenization techniques used in the literature for finding continuum models cannot be extended easily to origami-based metamaterials. Actually, it seems to us that the class of generalized continua to be introduced as target models in the searched homogenization process must have a rather richer kinematics than those used systematically, up to now, in the literature. Therefore, we decided to focus our analysis on discrete Lagrangian models for origami beams, in order to describe their mechanical behaviour in a very thorough and detailed way, to get a logical basis for conjecturing the structure of a future effective continuum model to be used for their overall mechanical description.

A brief description of Miura-ori and Kresling origami microstructures

A plethora of geometries for origami metamaterials-essentially borrowed from classical origami tessellationshave been proposed in the literature and consist of the periodic repetition in a two-dimensional array of a foldable unit cell. The most studied geometries, as long as mechanical behaviours are concerned, are the ones based on the so-called Miura-ori cell and the egg box tessellation. The Miura-ori tessellation owes its name to the fact that it was first introduced in the scientific literature by Miura [START_REF] Miura | Method of Packaging and Deployment of Large Membranes in Space[END_REF]. Indeed, as recalled above, Miura originally proposed this tessellation as an efficient geometry for packing and deploying solar panel sheets for space missions. It is however worth to mention that this kind of tessellation can be recognized in many biological systems, like leaves and embryonic intestine, which show geometric patterns that resemble the Miura-ori tessellation. Interestingly, such a tessellation also emerges in deformation patterns of surfaces subjected to biaxial compression, which somewhat points to its optimal properties.

To our knowledge, the study of the optimality of Miura-ori tessellation has not been fully explored yet: however, one can rather straightforwardly note that, if folding were allowed (only) in the most deformed regions, then the global stiffness would decrease, thus yielding an extremely compliant material that, once deformed, stores low deformation energy. It is in this regard relevant to mention that a pattern emerges when torsion tests are performed on thin-walled cylinders in paper or polypropylene: a thin sheet is wrapped around two coaxial mandrels, leaving a gap among the two of them and, when the mandrels are rotated about their axis, a folding pattern appears in the gap. Such a folding pattern is made up of inclined and stretched parallelograms, which are the mountain-folds, divided on their long diagonal by a valley-fold. This tessellation is known as the Kresling origami pattern [START_REF] Kresling | Natural twist buckling in shells: from the hawkmoth's bellows to the deployable Kresling-pattern and cylindrical Miura-ori[END_REF] and it is seen that, for any given gap size between the mandrels, a left-hand or right-hand twist-buckling pattern appears. Still regarding folding patterns emerging from mechanical tests on thin sheets, it is also worth to mention that Kresling has demonstrated how a cone, when smashed, can fold itself into a structure that seems to mimic the spiral geometries of pinecones and mollusk shells [START_REF] Kresling | The growing turbinate shell -model for a deployable technical shell[END_REF].

Related geometric and synthesis problems in the theory of origami metamaterials

In Miura's pioneering work, the original engineering problem, i.e. that of packing in a compact manner solar panels, was indeed overlying a fundamental mathematical-geometrical, specifically-problem which is at the very core of origami.

Such a problem reads as follows: given a surface, what are, if any, the transformations that are able to map rigidly, i.e. without stretch or compression, an intrinsically flat surface into a curved surface, i.e. a surface with intrinsic curvature? To better understand the previous problem statement it is worth to mention that by intrinsically flat surface we mean a surface that is completely flat and not a curved surface oscillating around a flat mid-plane.

This problem is strictly connected to the one of seeking to develop transformations mapping intrinsically flat surfaces into folded structures with apparent intrinsic curvature. This last problem can be regarded as analogous to the metamaterials' synthesis one, although also mechanics would in that case be involved. Such a problem can be also regarded as an inverse one that, generally speaking, is very difficult to solve. Notwithstanding these difficulties, amplified by limitations due to inherent developability constraints, the rational design of the crease and cut patterns proved several times successful in enabling the shape-shifting of (nearly) inextensible sheets into geometries with apparent intrinsic curvature, i.e. the folded sheets. Note that by apparent intrinsic curvature we mean the curvature possessed by the mid-surface of the folded sheets.

Mechanics of the Miura-ori tessellation

In the case of the Miura-ori tessellation, a doubly corrugated surface, obtained via the repetition in a two-dimensional array of a unit cell as represented in Fig. 1, is considered [START_REF] Schenk | Geometry of miura-folded metamaterials[END_REF][START_REF] Zhou | Origami mechanical metamaterials based on the miura-derivative fold patterns[END_REF]. The fundamental cell is made up of four congruent parallelograms, each one characterized by the lengths of its two sides sharing one common vertex and an angle. The four parallelograms are arranged in such a way that the sum of the four angles which meet at one vertex is equal to 2π, which allows to have a developable surface [START_REF] Kawasaki | On the relation between mountain-creases and valley-creases of a flat origami[END_REF]. As it can be seen in Fig. 1, the edges of the parallelograms form a re-entrant lattice. The geometry of the tessellation pattern can be described in terms of four parameters: the two lengths and the angle, specifying the parallelograms, and the angle between the flat plane and one of the parallelograms in the folded configuration. For Miura-ori tessellations, we can distinguish between in-plane and out-of-plane deformations, where the words "in-plane" and "out-of-plane" are referred to the plane passing through the valleys and mountains, i.e. resulting from the average of the mountain and valley array. As far as in-plane deformations are concerned, the behaviour of the Miura-ori tessellation is described by its inplane Poisson's ratio: one can immediately notice that the in-plane Poisson's ratio is always negative [START_REF] Greaves | Poisson's ratio and modern materials[END_REF], which can be easily deduced from the remark above: the edges of the parallelograms form a re-entrant lattice [START_REF] Yasuda | Reentrant origami-based metamaterials with negative Poisson's ratio and bistability[END_REF].

When out-of-plane deformations are considered, it is extremely difficult to predict the behaviour of the Miura-ori tessellation without making use of simulations based on a mechanical model. Among the several models that have been developed in the course of the past years, it is here worth to mention the simple-yet extremely successful and powerful-mechanical model introduced by Schenk and Guest in [START_REF] Schenk | Geometry of miura-folded metamaterials[END_REF][START_REF] Schenk | Folded Shell Structures[END_REF], which replaces the sheet with a system of bars interconnected via spherical pin-joints, and whose spirit is very close to the approach adopted in the present paper, as well as in relevant recent literature. Clearly, different descriptions are possible and, as an instance, the whole sheet could be modelled as a single thin shell or facets could be considered to move rigidly.

The presented description of the mechanical properties of the Miura-ori tessellation makes clear why it is so difficult to conceive a continuous model suitable to describe origami metamaterials based on the Miura-ori folding pattern at a macroscopic level: as an instance, macro-elongations are related mainly to relative micro-rotations and the macro-description of the system requires the possibility of accounting for these rotations. Moreover, the geometrical complexity of the Miura-ori tessellation makes very intricate the relationship between macro-and micro-kinematical descriptors, which mades the use of the micromacro-kinematical identification rather difficult in homogenizing such a microstructure [START_REF] Pideri | A second gradient material resulting from the homogenization of an heterogeneous linear elastic medium[END_REF][START_REF] Abali | Additive manufacturing introduced substructure and computational determination of metamaterials parameters by means of the asymptotic homogenization[END_REF].

Other origami microstructures can be useful in metamaterials design

As mentioned above, a plethora of origami-based metamaterials can be obtained by exploiting different origami tessellations and, besides the Miura-ori pattern, the so-called egg box tessellation pattern (see Fig. 2) is one of the most studied ones. This pattern is obtained by alternating along a two-dimensional array pyramids with a quadrilateral base having opposite orientations. Unlike the Miura-ori, this pattern cannot be obtained by suitably folding a flat sheet, since the sum of the angles at each vertex is less than 2π. From the point of view of its mechanical properties, differently from the Miura-ori tessellation, this pattern presents a positive in-plane Poisson's ratio. However, the previously mentioned models, and in particular the approach introduced by Schenk and Guest [START_REF] Schenk | Origami folding: a structural engineering approach[END_REF], can be applied also to the egg box pattern [START_REF] Schenk | Folded Shell Structures[END_REF]. It is worth to remark that an approach has been recently developed in [START_REF] Zhou | Origami mechanical metamaterials based on the miura-derivative fold patterns[END_REF], which allows to find the surfaces that can be tessellated via an egg box pattern, in the limit for the characteristic length of every pyramid going to zero.

Deformation energy for Miura-ori microstructures with non-rigid facets

The present paper attempts at studying computationally the nonlinear mechanical behaviour of origami metamaterials with non-rigid facets-facets are regions which the sheet is partitioned into, according to the crease pattern-making use of numerical tools, as well as of discrete spring model formulations, developed by the same authors in previous contributions to the literature on mechanical metamaterials [START_REF] Turco | Hencky-type discrete model for pantographic structures: numerical comparison with second gradient continuum models[END_REF][START_REF] Turco | A Lagrangian Hencky-type non-linear model suitable for metamaterials design of shearable and extensible slender deformable bodies alternative to Timoshenko theory[END_REF][START_REF] Turco | Forecasting nonlinear vibrations of patches of granular materials by elastic interactions between spheres[END_REF]. The sheet is assumed to have zero thickness, which is equivalent to state that we assume no deformation to be undergone by the sheet in the thickness direction. Notwithstanding the generality of the approach, the formulation is specialized to the case of Miura-ori tessellations, where facets, which are parallelograms, are modelled as two adjacent triangular sheets having as common edge a diagonal of the parallelogram (see Fig. 1). In the present approach, triangles are allowed to undergo (only) affine deformations in the space. Therefore, the position of their vertices in the current configuration makes up a minimal set of global Lagrangian coordinates for the mechanical system. From a mechanical viewpoint, the current lengths of these triangles are sufficient to write the deformation energy of the system due to the in-plane stretch and shear experienced by the triangles. In particular, for a homogeneous, isotropic sheet, the elastic energy of a triangle (with zero thickness) undergoing (in-plane) affine deformations can be equivalently thought of as if its sides were extensional Hooke's springs [START_REF] Steigmann | Two-dimensional models for the combined bending and stretching of plates and shells based on three-dimensional linear elasticity[END_REF], which will be accordingly called in the sequel stretch springs. Torsional springs obeying Hooke's law are placed in-between adjacent triangles, altogether called a panel, belonging to the same facet or to adjacent facets. The torsional springs placed in-between adjacent triangles belonging to the same facet are meant to model the bending behaviour of the facet, while those placed in-between adjacent triangles belonging to adjacent facets are meant to model the folding behaviour. In the sequel, according to such a distinction among torsional springs, we will refer to bending and folding springs, respectively.

For the sake of simplicity, the present contribution deals only with the nonlinear dynamical behaviour of origami metamaterials obtained as the periodic repetition of a Miura-ori unit cell along a single direction, i.e. with a one-dimensional array of Miura-ori unit cells. The present approach shares some features with the one presented recently by Liu and Paulino [START_REF] Liu | Nonlinear mechanics of non-rigid origami: an efficient computational approach[END_REF], even if, unlike [START_REF] Liu | Nonlinear mechanics of non-rigid origami: an efficient computational approach[END_REF], the present contribution introduces also inertial forces and, accordingly, uses a specific time-integration strategy to recover the evolution of the system.

Outline of the present paper

Having in mind the aforementioned goals, after this section, Sect. 2 introduces the kinematics as a global maximal set of Lagrangian coordinates of the system-node, i.e. triangles' vertices, positions-as well as the expression for the strain energy, consisting of contributions due to affine deformation of triangles in space, bending of facets, and folding of the sheet about creases, accounted for by stretch, bending, and folding springs, respectively. The expression of the kinetic energy is also introduced in Sect. 2. Expressions for strain and kinetic energies are discussed and motivated. Successively, in Sect. [START_REF] Barchiesi | Mechanical metamaterials: a state of the art[END_REF], what has been discussed in the previous section is used as a building block to formulate the time-stepping scheme, which is later used to perform numerical simulations presented and thoroughly discussed in Sect. 4. Concluding remarks and future challenges follow in Sect. 5. We finally note that the results obtained in the present paper suggest that the considered formulation can be fruitfully applied and developed further, aimed at exploring the class of behaviours that can be exhibited by origami metamaterials, possibly having tessellations different from the Miura-ori one and being obtained as the periodic repetition in two directions of unit cells that can be folded.

Specifying the model by means of its kinematics and associated energies

Aimed at exploring the mechanical behaviour of a generic origami metamaterial, a discrete model can be employed, having as fundamental brick a triangular element. Hence, strain and kinetic energies of the whole origami metamaterial can be computed by summing up over all the triangles, in which the origami geometry has been partitioned, these quantities.

In the discrete modelling that we adopt in this contribution, the reference configuration shall be completely defined by the vector X collecting the nodal coordinates, i.e. the coordinates of the vertices of all the triangles which form the origami, in the reference configuration. Since the nodes of an origami are, generally speaking, not lying all on the same plane, we cannot choose coordinate axes such that, a priori, the vector X has less than 3N n entries, where N n is the number of nodes. In the sequel, we will denote with N t the number of triangles that the origami has been decomposed into.

The current configuration shall be completely defined by the vector x collecting the nodal coordinates in the current configuration. Reference and current configurations allow to compute the nodal displacement vector u = x -X, having 3N n entries, which are the Lagrangian coordinates used to describe the motion in our approach.

In such a framework, each triangle, in which the origami metamaterial geometry has been decomposed into, is mapped into a (different) triangle in the current configuration, i.e. each triangle deforms affinely. In order to completely define the mechanical model, besides the kinematics, we need to define for given Lagrangian coordinates and velocities, i.e. the rate-in time-of Lagrangian coordinates, the strain and kinetic energies of the assembly of triangles forming the origami. Concerning the strain energy, we can distinguish:

1. The energy stored in the triangles-regarded as elastic bodies-because of the change in their shapes; 2. The energy stored in torsional springs, which can be in turn subdivided into bending and folding contributions, storing energy for relative rotations of adjacent triangles of the same facet about their common side and for relative rotations of adjacent triangles of two different facets about their common side, respectively; mimicking Hencky's idea for approximating through discrete modelling Euler's beam, see [START_REF] Hencky | Über die Angenäherte Lösung von Stabilitätsproblemen im Raum mittels der elastischen Gelenkkette[END_REF], triangles can only deform in plane and bending of facets, as well as folding around creases, is modelled in a lumped-parameter fashion (for a distributed-parameter approach see [START_REF] Hilgers | Energy-minimizing deformations of elastic sheets with bending stiffness[END_REF]). The two contributions to the energy mentioned above are analogous to-and are, actually, the natural extension to origami metamaterials of-those used to approximate the stretching and bending energies of an Euler-Bernoulli beam in two dimensions, see, e.g. [START_REF] Turco | Hencky-type discrete model for pantographic structures: numerical comparison with second gradient continuum models[END_REF]. In the present case, since the mechanical model describes two-dimensional structures moving in the three-dimensional space, it is necessary, even more than for the case of two-dimensional beams, to define the kinematics and the strain energy of the model as simply as possible.

Energy associated to in-plane deformation of triangular elements

Let us consider the e-th triangle of the origami and denote its vertices with the indices i, j, and k, all belonging to the set [1; N n ]. In the reference configuration, the coordinates of the vertices i, j, and k expressed in the global reference coordinate system are collected in the vectors X i , X j , and X k , respectively, see Fig. 3.

The hypothesis that triangular elements deform affinely means that strain is uniform along each triangle. In order to define the strain energy of each triangle in the simplest way, an ingenious idea introduced by Argyris and co-workers, see [START_REF] Argyris | TRIC: a simple but sophisticated 3-node triangular element based on 6 rigidbody and 12 straining modes for fast computational simulations of arbitrary isotropic and laminated composite shells[END_REF], where the deformation state is characterized by imagining to place three strain gauges on the triangle sides, i.e. a rosette of strain gauges, shall be exploited.

As an instance, let us consider what will be henceforth called the k-side, i.e. the side opposite to node k. The direction and length of this side are encoded within the vector dX k = X j -X i in the reference configuration and within the vector dx k = x jx i in the current one. The Green's strain measure g k of the k-side is defined as

g k = 2 k -L 2 k 2L 2 k , ( 1
)
where

2 k = dx k • dx k and L 2 k = dX k • dX k
are the squares of the current and the reference length of the k-side, respectively. Current positions of nodes i and j are x i = X i + u i and x j = X j + u j , respectively.

We define the stretching energy associated to the k-side as

S k = A k L k 1 2 Eg 2 k , (2)
where we have introduced the Young's modulus E of the material the triangle is made of and the crosssectional area A k of a fictitious rod placed along the k-side. Essentially, Argyris' idea consists of replacing the considered affinely-deforming triangle with an equivalent determinate truss whose elements are the sides of the triangle. Clearly, the total stretching energy associated to the three sides of the considered triangle is given by the sum of the stretching energies associated to each side

S e = S i + S j + S k , (3)
A practicable way to identify the parameters A i , A j , and A k consists of equating the strain energy E t of the considered triangle with that associated to the three sides. Remind that, for the sake of simplicity, the triangles' thickness s does not vary within each triangle and is the same for all triangles which the origami geometry has been partitioned into.

Concerning the Green's strain measure associated to the considered triangle's side, say the i-side, we can write

2 i -L 2 i = 2dX i • GdX i , ( 4
) where G denotes the Green-Saint-Venant strain tensor associated to the considered triangle, which is obtained by assuming that no deformation occurs along the thickness of the triangles.

Let us now consider a local Cartesian reference system with origin in the barycentre of the considered triangle, and two axes lying in the triangle plane. Then, the previous relationship can be written as

2 i -L 2 i = 2dX i • GdX i , ( 5
)
where the superimposed bars denote vectors and tensors in the local reference system. In such a reference system, only the strain measures G 11 , G 22 , and G 12 = G 21 , associated to the two axes lying in the triangle plane, are non-vanishing and, hence, needed to compute component-wise the Green-Saint-Venant strain tensor associated to the considered triangle. Therefore, we can write the following system of equations

2 i -L 2 i = 2dX i • GdX i , 2 j -L 2 j = 2dX j • GdX j , 2 k -L 2 k = 2dX k • GdX k , (6)
to find G 11 , G 22 , and G 12 = G 21 starting from the knowledge of the i-, j-, and k-side stretch. By utilizing the vector g = [G 11 G 22 G 12 ] T , the energy associated to in-plane deformation of the considered affinely-deforming triangular element can be written as

E e = 1 2 sA e g • Cg , ( 7
)
where A e and s denote the area1 and the thickness of the considered triangle, respectively. The matrix C represents in Voigt notation-in the local basis consisting of the two axes lying in the triangle plane-the elasticity tensor associated to the material the triangle is made up of, which is assumed to obey Hooke's law, to be homogeneous and isotropic, and be the same for all triangles. Such a matrix encodes the constitutive law, i.e. the material behaviour, for the plane strain case and reads as

C = E (1 + ν)(1 -2ν) ⎡ ⎢ ⎣ 1 -ν ν 0 ν 1 -ν 0 0 0 1 -2ν 2 ⎤ ⎥ ⎦ , ( 8
)
where the Poisson's ratio ν of the material has been introduced. By equating the total stretching energy associated to the sides of the considered triangle, given by Eq. ( 3), and the energy associated to its in-plane deformation, see Eq. ( 7), we can identify the crosssectional area of the fictitious rods making up the determinate truss system associated to the considered triangle so as to reach an exact energy equivalence between the two mechanical descriptions.

Bending and folding energies

The strain energy associated to bending of a facet and that associated to folding of a crease can be chosen in an approximate, yet simple, way following the key idea proposed by Hencky in his celebrated doctoral thesis, see [START_REF] Hencky | Über die Angenäherte Lösung von Stabilitätsproblemen im Raum mittels der elastischen Gelenkkette[END_REF], aimed at estimating the buckling load of an Euler's beam. A revival of Hencky's respectively, then the triangle area can be computed as

Ae = p(p -L i )(p -L j )(p -L k ), where p = 1 2 (L i + L j + L k ),
i.e. it is equal to half the triangle perimeter. idea can be found in a series of articles aiming at modelling pantographic structures [START_REF] Turco | Hencky-type discrete model for pantographic structures: numerical comparison with second gradient continuum models[END_REF], beams [START_REF] Turco | Discrete is it enough? The revival of Piola-Hencky keynotes to analyze three-dimensional Elastica[END_REF][START_REF] Turco | In-plane dynamic buckling of duoskelion beam-like structures: discrete modeling and numerical results[END_REF][START_REF] Baroudi | Nonlinear dynamics of uniformly loaded elastica: experimental and numerical evidence of motion around curled stable equilibrium configurations[END_REF], and granular materials [START_REF] Turco | Forecasting nonlinear vibrations of patches of granular materials by elastic interactions between spheres[END_REF][START_REF] Turco | In-plane shear loading of granular membranes modeled as a Lagrangian assembly of rotating elastic particles[END_REF][START_REF] Turco | A nonlinear Lagrangian particle model for grains assemblies including grain relative rotations[END_REF]. Roughly speaking, Hencky's idea consists of stating the equivalence between the bending-and, in our case, folding-energy stored in a mechanical element by subdividing such an element in two adjacent elements hinged through cylindrical hinges in such a way that only relative rotations around the (considered) bending (or folding) axis is allowed. Resistance to such a relative rotation is added by placing a-Hooke's in the classical approach, which we adopt also in the present contribution-torsional spring in-between the two elements, thus yielding a lumped-parameter model. In this paper, we adopt this very same concept to evaluate the energy stored in a panel, i.e. two triangles sharing a side, because of facet bending or origami folding about creases. Figure 4 reports a generic panel of the considered origami in the current configurations, where m and n stand for the indices associated to the triangles in the panel.

The vertices of the triangles in the considered panel, also called nodes, are labelled with the indices i, j, k, and h. Reference and current positions of each node are denoted by X and x, respectively, and labelled with the associated index, i.e. i, j, k, or h. On each triangle we can consider three mutually orthogonal unit vectors, say {d All the sets of unit vectors can be evaluated from node positions by using standard geometry formulas. As an instance, the three mutually orthogonal unit vectors related to the m-th triangle in the current configuration can be computed as

d 1m = x j -x i x j -x i , d 2m = d 3 × d 1 , d 3m = (x k -x j ) × (x j -x i ) (x k -x j ) × (x j -x i ) . ( 9
)
Both in the reference and current configurations of the m-th triangle, the unit vector in Eq. ( 9) having 1 as subscript is parallel to the k-side of the triangle, that with subscript 3 is orthogonal to the plane containing the considered triangle and, finally, that with subscript 2 is such that

d 2m = d 3m × d 1m .
We now aim at defining a strain measure associated to bending of facets or folding of the origami around creases. The underlying idea is to measure, for each panel, the relative rotation of the two adjacent triangles making up a single panel. It is in this regard worth to be mentioned that, since the local bases {d 1m , d 2m , d 3m } and {d 1n , d 2n , d 3n } are built making use of one-remark, not the shared one-side of the stretched triangles m and n, respectively, the rotation leading from {d 1m , d 2m , d 3m } to {d 1n , d 2n , d 3n } is not occurring around the side shared by the triangles m and n, i.e. its eigenvector does not identify the direction of the shared side, thus not representing an ideal strain measure associated to bending of facets or folding of the origami around creases in the Hencky-type approach mentioned at the beginning of this subsection. To get a measure of the relative rotation of the two adjacent triangles making up a single panel that is fulfilling this desired property, we seek to write the rotation tensors associated to passing from {D 1m , D 2m , D 3m } and {D 1n , D 2n , D 3n } to {d 1m , d 2m , d 3m } and {d 1n , d 2n , d 3n }, respectively, by filtering out the (spurious) rotations induced by the in-plane affine deformation of the triangles m and n, respectively. In order to do so, we make use of the right polar decomposition of the deformation gradient tensor F associated to the deformation in three space dimensions of the triangles m and n.

Referring to Fig. 4 the tensor F m for the m-th triangle can be computed by making use of the following equations:

F m (X j -X i ) = x j -x i F m (X k -X j ) = x k -x j F m D 3m = d 3m (10)
that represent a compatible linear system of nine scalar linearly independent algebraic equations in the nine independent components of F m , which are unknown. Such a system admits an elementary closedform solution. The same holds for F n . At this point, it is only necessary to recall the right polar decomposition formula

F = RU , ( 11
)
where R is a proper orthogonal tensor representing the rotational part of the deformation gradient F and U is a positive definite symmetric tensor representing the stretch part of the deformation gradient F.

Since F T F = U T R T RU = U T U = U 2 (12)
then U = √ F T F. Reminding that the square root of a tensor can be computed by means of its spectral decomposition, that is

√ F T F = V √ ΛV T , ( 13
)
being V and Λ the orthogonal matrix collecting the eigenvectors of the symmetric tensor F T F and the diagonal matrix whose elements are the eigenvalues of F T F arranged in the same order as that of the corresponding eigenvectors in the matrix V, respectively. Finally, when U is known, we can compute the pure rotation matrix R = FU -1 .

Having hence defined the rotation matrices R m and R n for the m-th and n-th triangles in the p-th panel, the sought finite relative rotation tensor ΔR can be defined as

ΔR = R T m R n . ( 14
)
Rodrigues' formula allows to rewrite in a more insightful way a rotation tensor R as R = cos ϕI + (1cos ϕ)e ⊗ e + sin ϕE [START_REF] Abali | Theory and computation of higher gradient elasticity theories based on action principles[END_REF] being E and e the tensor and its axial vector fulfilling the relationship Ev = e × v for any vector v. From Eq. ( 15), simple computations give

2 cos ϕ = tr(R) -1 , 2 sin ϕ e = R × , ( 16
)
where tr(•) is the trace operator and R × is the axial vector of the skew symmetric part of R.

For computational purposes, it is convenient to introduce the rotation vector

ϑ ϑ ϑ = 2 tan ϕ 2 e , ( 17
)
which allows to rewrite the rotation tensor R as

R = 1 4 + ϑ 2 (4 -ϑ 2 ) I + 2ϑ ϑ ϑ ⊗ ϑ ϑ ϑ + 4E ϑ ϑ ϑ , ( 18
)
being

ϑ 2 = ϑ ϑ ϑ •ϑ ϑ ϑ and E ϑ ϑ ϑ v = ϑ ϑ ϑ × v.
The version of the Rodrigues' formula in Eq. ( 18) is equivalent to [START_REF] Abali | Theory and computation of higher gradient elasticity theories based on action principles[END_REF], but allows to avoid singularities that would otherwise arise when algebraic manipulations are performed, necessary to define quantities such as the strain energy and its first and second derivatives. In addition, the use of R defined by equation ( 18) makes the rotation vector vary in the same range, from -∞ to +∞, of the displacements vector and therefore, in principle, is beneficial for the well-conditioning of the equilibrium equations to be solved. Since the relative finite relative rotation tensor ΔR is orthogonal, by using ( 16) and the half-angle trigonometric formula, we can write the relative finite rotation vector

Δϑ ϑ ϑ = 2(ΔR) × 1 + tr(ΔR) (19)
and, finally, the folding strain energy for the p-th panel

B p = 1 2 b X k -X j Δϑ ϑ ϑ • Δϑ ϑ ϑ , (20)
being b the stiffness per unit length which controls the relative rotation in-between the two triangles in the considered panel. In principle, the quantity b can assume different values according to whether it is associated to bending or folding.

From the bending, folding, and stretching strain energies, we can compute the structural reaction vector s and the tangent stiffness matrix K of the system. Since these strain energies are additive quantities, we can define the total strain energy just for a single generic panel p, i.e. for two triangles, say m and n-see Fig. 4-sharing a side.

Let us denote with the symbol B the sum over the whole origami of the bending and folding strain energies associated to each panel and let us denote with the symbol S the sum over the whole origami of the sides' stretching energy. We define the total strain energy associated to the origami as the sum B + S. Collecting the displacements of the nodes i, j, k, and h in u p , we can write the structural reaction vector of the p-th panel as

s p = d(B + S) du p , ( 21
)
while its tangent stiffness matrix reads accordingly as

K p = ds p du p . ( 22
)
Definitions ( 21) and ( 22) provide us with the tools to evaluate automatically the structural reaction vector and the tangent stiffness matrix associated to the total strain energy of the p-th panel. A simple assembly operation, following rules identical to those used in finite element codes, is necessary to build the structural reaction vector and the tangent stiffness matrix of the whole origami.

Kinetic energy

The kinetic energy associated to the e-th triangle of the origami can be written as

K e = 1 2 Ve ρ u • u dV , ( 23
)
where ρ, V e , and u are the mass density, the three-dimensional domain the triangle is associated tonamely the extruded triangle-and the velocity of the elementary volume dV , respectively. We suppose the mass density and, as mentioned before, the thickness s of the triangle as uniform along the triangle.

Let us choose a linear interpolant to approximate the velocity of the elementary volume dV . Then the kinetic energy of the considered triangular element can be approximated as

K e ≈ 1 2 Ae ρ s (B e ue ) • (B e ue ) dA = 1 2 ρs ue • ⎛ ⎝ Ae B T e B e dA ⎞ ⎠ ue , ( 24
)
where A e indicates the considered triangle. The shape function matrix B e and the nodal velocity vector ue allow to compute the velocity vector in each point of the triangle, i.e. u = B e ue . The shape function matrix B e for a generic triangle can be defined in an indirect way by defining it just for a specific triangle, as it going to be done in Eq. ( 26), while the velocity vector reads as

ue = ⎡ ⎣ ui uj uk ⎤ ⎦ . ( 25
)
The computation of the mass matrix M e of the considered triangular element, considering that the factor ρs is uniform in space, requires the evaluation of the integral over the triangular element of the term B T e B e . Such an evaluation becomes straightforward if the integral is performed on the so-called reference triangle, i.e. the triangle having, in the Cartesian plane, as vertices the points with coordinates (0, 0), (1, 0) and (0, 1). For such a special triangle, the shape function matrix reads in formulas as

B e = ⎡ ⎣ (1 -ξ -η) 0 0 ξ 0 0 η 0 0 0 ( 1-ξ -η) 0 0 ξ 0 0 η 0 0 0 ( 1 -ξ -η) 0 0 ξ 0 0 η ⎤ ⎦ ( 26
)
where 0 ≤ ξ ≤ 1 and 0 ≤ η ≤ 1.

Having introduced these tools, we are now able to write

) 27
where Âe denotes the special triangle specified above and, with an abuse of notation, A e denotes the area of the e-th triangle. More explicitly, the triangular element mass matrix reads as

M e = ρsA e 12 ⎡ ⎣ 2 I 3×3 I 3×3 I 3×3 I 3×3 2 I 3×3 I 3×3 I 3×3 I 3×3 2 I 3×3 ⎤ ⎦ (28)
where, aimed at keeping as short as possible the formula, the 3 × 3 identity matrix, denoted with I 3×3 , has been utilized.

Stepwise approximate solution of the nonlinear equilibrium problem

By referring to the original work of Casciaro [START_REF] Casciaro | Time evolutional analysis of nonlinear structures[END_REF] or to the successive paper by Turco [START_REF] Turco | Stepwise analysis of pantographic beams subjected to impulsive loads[END_REF], where formulas are derived in a detailed way, commented, and light is shed on some otherwise unaddressed issues, we briefly recall here the main elements of the time integration scheme that we have utilized to solve in a stepwise fashion the nonlinear equilibrium problem deriving from the modelling presented in the previous sections and obtain the results that will be presented in the next one. First of all, we shall specify that the integration scheme is based on subdividing in N time steps of equal length Δt the desired time horizon of the simulation. Starting from the knowledge of the initial conditions, i.e. the initial displacement vector u 0 = u(t = 0) and the initial velocity vector u0 = u(t = 0), we want to reconstruct the solution at times t j+1 = t j + Δt (j = 0, 1, . . . , N). Since we are considering a recursive time integration scheme, and hence we can assume as known the solution for t j , i.e. we know the vectors u j = u(t j ) and uj = u(t j ), the scheme can be characterized by describing how to compute the solution u j+1 and uj+1 .

For each time step, since we are considering as unknowns the vectors u j+1 and uj+1 , two sets of equations are needed so as to equate the number of independent scalar equations and that of the scalar unknowns. Referring to Casciaro's scheme, the following sets of equations are considered:

u j+1 = u j + 1 2 -β uj + 1 2 + β uj+1 Δt , ( 29
)
and

M( uj+1 -uj ) + 1 2 -α (s j -f j ) + 1 2 + α (s j+1 -f j+1 ) Δt = 0 . ( 30
)
Equation ( 29) reflects the assumption that the unknown displacement vector uj+1 depends, besides the displacement vector u j , upon the velocity vectors uj and uj+1 weighted by two dimensionless factors that, in turn, are depending upon the parameter β, which has to be chosen suitably. Readers who are familiar with B-splines can recognize that Eq. ( 29) is resembling that kind of interpolation. Indeed, in that equation, the displacement vector is extrapolated quadratically using its time derivative, i.e. the velocity, evaluated both at the beginning and at the end of the considered time step. 2Equation ( 30) prescribes the equilibrium by means of the so-called momentum-impulse relationship. Indeed, the first term of Eq. ( 30) can be interpreted as a finite-difference approximation of the rate of momentum, while it can be recognized that the second term is the average net impulse. This last is weighted by two dimensionless factors depending upon the parameter α, which has to be selected suitably.

To summarize, Casciaro's scheme makes use of a quadratic B-spline interpolation for the displacements and of the momentum-impulse relationship, introducing two dimensionless parameters α and β that can be finely tuned to obtain the best performances from the integration scheme.

These two dimensionless parameters, following Casciaro's spirit, have to be tuned according to the guess for the first, T 1 , and for the last, T N , natural periods of the considered system (with essential boundary conditions) and the chosen time-step length Δt that, generally speaking, must be in turn chosen on the basis of the considered problem (the first natural period, usually). Therefore, operationally speaking, an eigenvalue analysis limited to the first and the last natural periods computed in the reference configuration has to be performed. This preliminary analysis involves the mass matrix M and the stiffness matrix K computed in the reference configuration, i.e. K(u = 0). In Casciaro's paper [START_REF] Casciaro | Time evolutional analysis of nonlinear structures[END_REF], a rule for choosing in an optimal way the dimensionless parameters α and β is presented. Such a rule is based on the analysis of a one-degree-of-freedom system in the linear regime.

However, the work [START_REF] Casciaro | Time evolutional analysis of nonlinear structures[END_REF] illustrates only a few details concerning the derivation of this rule, while more extensive discussions are carried out in a series of reports published some years before, which cannot be easily found. This, along with the fact that the original formulas provided by Casciaro contained some typos, caused, most likely, that Casciaro's idea was eventually forgotten. 3The optimal values of the dimensionless parameters α and β are

β = -α = - 1 4 + 1 γ 2 - 1 + 1 + tan 2 γ 2 tan 2 γ , 0 < γ < π 2 , β = -α = - 1 4 + 1 γ 2 - 1 -1 + tan 2 γ 2 tan 2 γ , π 2 < γ < π , (31)
having introduced γ = 2πΔt/T n . We remark that Eq. ( 31), when Δt → 0, leads to β = -α = 1/ √ 6. Instead, when Δt → T n /2, β = -α = 1/π is obtained.

Conversely, when Δt > T n /2, the optimal values assume the form

α = - T n 2πΔt + c 3 1 + 2c 3 , β = T n 2πΔt + c 3 1 + 2c 3 , c = 2Δt -T n T 1 -T n . ( 32
)
Formulas ( 32) contain both the sum of two terms: the first one is a hyperbola branch, which links the optimal values of α and β at the end of the interval Δt < T n /2; the second one is the ratio between two cubic polynomials, which tends-when Δt → ∞-to the value 1/2. Hence, the values of α and β corresponding to this case are able to reproduce a quasi-static solutions. 4 At this point, we shall sketch the solution strategy, which we have utilized, encoded in an algorithm, to solve-for each time step-the nonlinear equations prescribing the equilibrium, thus overcoming the difficulties due to the geometric nonlinearities included in structural reaction. The key ingredient of the solution strategy is the first order Taylor expansion of the structural reaction s j+1 = s(u j+1 ). In formulas, we have:

s j+1 ≈ s j + K j (u j+1 -u j ) , (33)
where the stiffness matrix K j is computed at the beginning of the considered time interval, i.e.:

K j = ds du uj . ( 34
)
Remark that the velocity variables collected in the vector uj+1 can be eliminated by using Eq. (29) from the system formed by Eqs. ( 29) and [START_REF] Zhou | Origami mechanical metamaterials based on the miura-derivative fold patterns[END_REF]. Therefore, Equations ( 29)- [START_REF] Zhou | Origami mechanical metamaterials based on the miura-derivative fold patterns[END_REF], along with the Taylor expansion [START_REF] Schenk | Folded Shell Structures[END_REF], give the necessary tools to build a Newton-like algorithm for the considered problem. Starting from the formula expressing the reminder calculated at the i-th iteration

r i = M( uj+1,i -uj ) + 1 2 -α (s j -f j ) + 1 2 + α (s j+1,i -f j+1 ) , (35)
an update estimate of the solution, using a Newton-like strategy, is given by the recursive relation

u j+1,i+1 = u j+1,i -H -1 j r j , (36)
where the iteration matrix H j is defined as

H j = ∇r| j = 1 1 2 + β M + 1 2 + α K j,i . (37)
For addressing the solution of problems involving origami metamaterials, modelled as explained in the previous sections, we used a slight modification of Casciaro's scheme. Essentially, we used as primary variable the displacement vector instead of the velocity vector. This does not change significantly the properties of the scheme but, in practice, in view of addressing complex problems such as those involving origami metamaterials, it allows to get an algorithm that is significantly more stable.

At the end of this section, it is worth remarking that, for dynamic problems-even complex ones as those analysed in the present contribution-using the Newton method-as we do here to predict the behaviour of origami metamaterials-is not a major limitation. This could seem counter-intuitive since, as it is well known, when static problems are considered, Newton's method is not able to get around limit points, i.e. configurations that make singular the stiffness matrix K, because for such a method the iteration matrix coincides with the tangent stiffness matrix. Indeed, for this reason, its use for the static analysis of systems exhibiting limit points is not recommended. When dynamic analyses are concerned, since the iteration matrix also includes an additive contribution depending on the mass matrix (which does not vanish when the stiffness matrix is singular), see [START_REF] Turco | Hencky-type discrete model for pantographic structures: numerical comparison with second gradient continuum models[END_REF], the problem is instead regularized, thus making the Newton scheme usable to address complex systems exhibiting limit points, which is the case we are going to analyse in the next section.

Numerical simulation of origami beams

Using the model that has been sketched in the previous sections, we critically present in this section the results of some numerical simulations concerning straight and curved origami beams, in short o-beams, i.e. beams obtained as the repetition of a foldable unit cell in a mono-dimensional array, highlighting the peculiarities of such a kind of beams. Specifically, concerning the mechanical problems that we are going to study, we shall address the problem of a straight beam subjected to (1) a transverse shearing force, (2) an axial compression force, and (3) a torque at the free end, which are varying in time with an assigned law. Additionally, we analyse the case of a simply supported circular arch subjected to a time-varying shearing force applied in its mid-point.

Remark that, in this section, all the data and results, when not in non-dimensional form, are expressed in the MKS system of units.

Cantilever o-beam subjected to a shearing force at its free end

From now on, we study the straight origami beam sketched in Fig. 5 (where X-, Y -, and Z-view means view from the X, Y , and Z axis, respectively), obtained repeating 25 times along a straight unidimensional array the Miura-ori unit cell. In this subsection, let us focus on the so-called cantilever problem. The kinematic boundary conditions applied at the nodes of the clamped end prescribe vanishing nodal displacements. At the free end of the o-beam, which consists of the three nodes having vanishing X coordinate (see Fig. 5d), a shearing force acting in the plane having Y as normal, varying in time according to the law reported in Fig. 6, is applied on the mid-node.

As described in Sect. 3, the utilization of Casciaro's scheme requires the choice of two parameters. In order to choose optimally these parameters, we evaluate the first and the last natural periods of the considered o-beam for the studied problem. For the problem studied in this subsection, with geometrical and material parameters chosen as those reported in Table 1, we obtain that the first natural period is T 1 = 1.76, while the last one is T n = 1.90 × 10 -4 . In Fig. 7, we report the first nine natural modes.

It is worth to remark that, among the first nine modes in Fig. 7, there are modes associated to bending of the o-beam in its plane of symmetry, having as normal vector the Y -axis, specifically the 1st, 3rd, 5th, and 7th mode. The 2nd, 4th, and 6th mode are instead associated to bending in the plane having as normal vector the Z-axis. The 8th and 9th mode are instead associated to extension and torsion, respectively. The time evolution of the system, for the shearing load p max = 5000, t i = 1, and t d = 0.25, computed using a time integration step Δt = 0.001 is summarized in Fig. 8. The quantities u, v, and w in Fig. 8 are the displacement components along the axes X, Y , and Z, respectively, of the mid-node of the free end of the o-beam.

Henceforth, colours will be used in the stroboscopic pictures to indicate the time evolution. Specifically, ten different shots, equally spaced in time, will be taken to build the stroboscopic pictures. For each shot, the employed percentage of blue colour indicates the progression of time: grey corresponds to the reference configuration, light violet corresponds to the last shot.

It is worth to be remarked that the applied load causes an unexpected oscillation in the Y direction, as it can be clearly seen from the yellow line in Fig. 8a,d. This oscillation is triggered by the unavoidable presence of numerical errors. While the numerical scheme is unable to deal with this phenomenon, as no numerical viscous dissipation has been introduced to damp undesired oscillations, the observed oscillations in the Y direction hint to the high sensitivity of the o-beam to imperfections.

A videoclip showing the cantilever o-beam oscillations induced by the time-varying shearing load applied at the mid-point of its free end, considered in this subsection, can be found in the supplementary material, see the file c25bending.avi.

Cantilever o-beam subjected to compression forces at its free end

Using the same data of the previous test, we now study the problem of an o-beam with cantilever kinematic boundary conditions, subjected to a compression force along the X direction acting at its free end. Specifically, such a force is applied to all the nodes at the free end.

The maximum amplitude reached by the load, which is varying in time according to the law shown in Fig. 6, is equal to 8000. Using a time integration step Δt = 0.001, we obtain the plots shown in Fig. 9.

As for the previous test problem, also in this case the applied load produces an unexpected oscillation in the Y direction, due to numerical errors, as it can be clearly seen by looking at the yellow line in Fig. 9(a,d).

A videoclip showing the cantilever o-beam oscillations induced by the time-varying compression loads applied at the nodes of its free end, considered in this subsection, can be found in the supplementary material, see the file c25buckling.avi. As it can be seen looking at the red line in Fig. 10a, the torque applied at the free end produces oscillations also in the X direction, which hints at the nonlinear coupling between torsional and extensional deformation modes in the o-beam.

A videoclip showing the o-beam oscillations observed from the X-axis for the case analysed in the present subsection can be found in the supplementary material, see the file c25twisting.avi.

Simply-supported origami arch subjected to a shearing force at its mid-point

We now study the origami arch obtained by repeating along a curved line 25 unit cells, see Fig. 12. Specifically, the bottom-most nodes of the origami arch are placed along a circumference of radius 100. Material parameters considered in this subsection for the origami arch are the same reported in Table 1. Displacement of the mid-nodes of the two ends is completely prevented. For the remaining four nodes of the two ends, we prevent only displacement in the vertical direction, i.e. along the Z-axis. Since the geometry and kinematic boundary conditions have changed with respect to those analysed in the previous tests, we compute the first and the last natural periods to find the optimal parameters of the time-integration scheme. The eigenvalue analysis gives for the first and the last natural periods the values T 1 = 0.221 and T n = 7.97 × 10 -5 , respectively.

Figure 12 reports the first nine natural modes and periods. We remark that the first nine modes are all flexural modes. The odd modes correspond to deformations out of the plane which the origami is lying onto in the initial configuration, i.e. the X -Y plane, while the even modes lie in the initial origami plane.

The arch is loaded at its mid-point with a vertical force, directed towards the interior of the circumference, that is increasing in time with constant, unit velocity. Looking at Fig. 13, we observe an instability which resembles that experienced by the Von Mises arch when loaded quasi-statically. In the present case, the instability is much more complex, since occurs in dynamics and both in the plane normal to the Y -axis and in the plane normal to the Z-axis. It is worth to remark that the employed Newton scheme, which is used for a given time step to solve the nonlinear equilibrium equations, has proven successful in dealing with limit points since the iteration matrix also includes an additive contribution depending on the mass matrix.

Also for the case addressed in the present subsection a videoclip showing the deformation in time of the origami arch is included in the supplementary material, see the file arch25.avi.

Concluding remarks and future challenges

In this paper, we have presented a viable way to perform numerical simulations of origami metamaterials taking into account inertial effects and nonlinearities entering into play when large displacements are concerned. Employing a time integration algorithm inspired by a contribution by Casciaro, almost unrecognized in the literature, we were able to reconstruct the entire nonlinear equilibrium path. The key ingredients of the approach presented in this contribution are the definition of strain and kinetic energies of the fundamental brick of the origami metamaterial, i.e. what we have called the panel. For each triangle which the facets have been partitioned into, only affine transformations are considered. Therefore, the strain energy related to the in-plane deformation of triangles can be computed by making use of information concerning the stretch of their edges only, while that related to bending and folding, associated to a panel, i.e. two triangles sharing an edge, can be computed by knowing only the relative rotation of the two triangles around their common edge.

The general formulas reported in [START_REF] Liu | Nonlinear mechanics of non-rigid origami: an efficient computational approach[END_REF][START_REF] Van Schaik | A structure refinement method based on molecular dynamics in four spatial dimensions[END_REF][START_REF] Bekker | Molecular Dynamics Simulation Methods Revised[END_REF], which are making use of the change in the dihedral angle between two triangles in a panel, are simple and computationally efficient; these characteristics are essential for a code that should be able to deal with problems involving a large number of degrees of freedom. However, the approach presented in [START_REF] Liu | Nonlinear mechanics of non-rigid origami: an efficient computational approach[END_REF][START_REF] Van Schaik | A structure refinement method based on molecular dynamics in four spatial dimensions[END_REF][START_REF] Bekker | Molecular Dynamics Simulation Methods Revised[END_REF] is affected by some singularities that pose difficulties in computing the structural reaction and the stiffness matrix for some configurations using automatic symbolic computation. In this regard, it is worth to mention in this concluding section that, similarly to what has been done in [START_REF] Turco | Modeling of three-dimensional beam nonlinear vibrations generalizing Hencky's ideas[END_REF], a tensorial definition of such a relative rotation has been introduced in this paper, which is not biased by in-plane affine deformation of triangles and can get around the singularities affecting the approach presented in [START_REF] Liu | Nonlinear mechanics of non-rigid origami: an efficient computational approach[END_REF][START_REF] Van Schaik | A structure refinement method based on molecular dynamics in four spatial dimensions[END_REF][START_REF] Bekker | Molecular Dynamics Simulation Methods Revised[END_REF].

In other words, we have generalized some ideas developed to model beams and beam lattices in the so-called Hencky's approach, presented in works devoted to formulate simple, yet effective and computationally efficient, mechanical models for metamaterials [START_REF] Turco | Hencky-type discrete model for pantographic structures: numerical comparison with second gradient continuum models[END_REF]. We have shown that numerical simulations performed for origami beams folded according to the Miura-ori pattern reveal interesting features. The more relevant one is the sensitivity to imperfections, which causes unexpected vibrations, as it is possible to see by looking at the results presented in Sect. [START_REF] Miura | Method of Packaging and Deployment of Large Membranes in Space[END_REF].

Concerning the outlooks of the present work, it is worth to mention the following future challenges that could be tackled:

1. This work has presented, in addition to the stepwise solution strategy, a quite general tool to compute the kinetic and elastic deformation energies of a triangle moving in space, as well as the bending and folding energies of a panel; summing up these contributions, it is straightforward to obtain the total energy associated to a unit cell and, consequently, that associated to the whole origami metamaterial; such a tool could be used in the future to analyse the vibrations of origami plates and shells (see Fig. 14) with an arbitrary folding pattern, as an instance the so-called Kresling pattern [START_REF] Kresling | Origami-structures in nature: lessons in designing "smart" materials[END_REF] (see Fig. 15, where we show a depiction of the so-called Kresling tower ), which appears interesting for its wide application field. Regarding Kresling folding patterns, it is worth to mention that they show some analogies with the so-called duoskelion beam, introduced in previous works by the authors of the present contribution [START_REF] Barchiesi | Equilibria determination of elastic articulated duoskelion beams in 2D via a Riks-type algorithm[END_REF][START_REF] Turco | A numerical investigation on impulse-induced nonlinear longitudinal waves in pantographic beams[END_REF]-exploiting chirality of the microstructure [START_REF] Lorato | The transverse elastic properties of chiral honeycombs[END_REF]-in that they have a very sensitive buckling behaviour, see [START_REF] Eremeyev | Enriched buckling for beam-lattice metamaterials[END_REF]. In this regard, future works elaborating on top of the present contribution could deal with the exploration of the role of local Euler buckling in Kresling origami structures during pattern formation, induced by external loading [START_REF] Kresling | The fifth fold: complex symmetries in Kresling-origami patterns[END_REF]. 2. The extremely simple form of bending and folding energies used in this wok does not ensure that relative rotations of two triangles in a panel making the triangles overlap are avoided, nor kinematic variables are constrained so to prevent such an occurrence; the way proposed by [START_REF] Liu | Nonlinear mechanics of non-rigid origami: an efficient computational approach[END_REF] to address this issue, by making use of an energy barrier, is interesting, but alternative routes could be explored in the future.

3. Despite the fact that the presented modelling approach has been here developed to perform dynamical analyses of origami metamaterials, the same approach, which is intrinsically discrete, could be exploited in the future to analyse numerically standard shells, taking advantage of its computational efficiency, see [START_REF] Turco | Elasto-plastic analysis of Kirchhoff plates by high simplicity finite elements[END_REF][START_REF] Oñate | Rotation-free triangular plate and shell elements[END_REF]. We expect this to be particularly profitable when more complex material behaviours, as an instance damage and plasticity, are considered [START_REF] Turco | Elasto-plastic analysis of Kirchhoff plates by high simplicity finite elements[END_REF][START_REF] Timofeev | Hemivariational continuum approach for granular solids with damageinduced anisotropy evolution[END_REF][START_REF] Placidi | Micromechanics-based elasto-plastic-damage energy formulation for strain gradient solids with granular microstructure[END_REF]. 4. Finally, the authors believe that results reported in the recent paper [START_REF] Eremeyev | Minimal surfaces and conservation laws for bidimensional structures[END_REF] concerning the computation of some conservation laws along edges might improve the computational efficiency of the code built to analyse the nonlinear behaviour of origami-like structures.

Fig. 1 .

 1 Fig. 1. Standard Miura unit cell: before (left) and after (right) the folding

Fig. 2 .

 2 Fig. 2. Pyramidal egg box unit cell (left) and pattern (right)

Fig. 3 .

 3 Fig. 3. In-plane deformation of a triangular element

Fig. 4 .

 4 Fig. 4. Origami panel in the current configuration: nodes and normal unit vectors

  1 , d 2 , d 3 } m for the m-th triangle and {d 1 , d 2 , d 3 } n for the n-th triangle in the current configuration. The dihedral angle in-between the planes containing the triangles is denoted with the greek letter δ. An analogous figure can be traced for the reference configuration. In this case, we would use capital letter to indicate the mutually orthogonal unit vectors, i.e. {D 1 , D 2 , D 3 } m for the m-th triangle and {D 1 , D 2 , D 3 } n for the n-th triangle. The dihedral angle in the reference configuration is denoted with δ 0 .

Fig. 5 .Fig. 6 .

 56 Fig. 5. Views of the o-beam unit cell-X-(a), Y -(b), Z-(c), and 3D-view (d)-and whole o-beam made up of 25 unit cells (e)

Fig. 7 .

 7 Fig. 7. First nine natural periods and modes for the o-beam subjected to cantilever boundary conditions. Gray and red colors indicate the reference configuration and the modes, respectively

Fig. 8 .

 8 Fig. 8. Origami cantilever beam subjected to a shearing force on the free end mid-point. Evolution in time of displacements (a) and velocities (b) along the X, Y , and Z directions of the mid-point at the free end of the o-beam. Stroboscopic pictures of the o-beam deforming in time observed from -Y (c) and -X (d)

Fig. 9 .

 9 Fig. 9. Origami cantilever beam subjected to compression forces on the free end nodes. Evolution in time of displacements (a) and velocities (b) along the X, Y , and Z directions of the mid-point at the free end of the o-beam. Stroboscopic pictures of the o-beam deforming in time observed from -Y (c) and -X (d)

Fig. 11 .

 11 Fig. 11. Origami arch with 25 unit cells

Fig. 12 .Fig. 13 .

 1213 Fig. 12. First nine natural periods and modes for the simply-supported origami arch. Gray and red colors indicate the reference configuration and the modes, respectively

Fig. 14 .

 14 Fig. 14. Plate (a), vault (b), shell (c), and half-sphere (d) with standard Miura-ori folding pattern

Fig. 15 .

 15 Fig. 15. Depiction of the so-called Kresling tower

Table 1 .

 1 Utilized dataset of geometrical and material parameters in the MKS system of units

	Cells	dX	dY	dZ	s	ρ	E	ν	b 1	b 2
	25	5	5	5	0.1	10 -2	10 7	0.3	10 6	10 6

In order to compute the area of the e-th triangle, Heron's formula might be used, which makes use of the lengths of the triangle's sides. If L i , L j and L k denote the lengths of the sides that are opposite to the triangle's vertices i, j and k,

The idea of representing a quadratic polynomial, or a polynomial of higher order, using both the values of the polynomial and of its derivatives is widespread in the field of B-splines, see[START_REF] Greco | B-Spline interpolation of Kirchhoff-Love space rods[END_REF], and NURBS, see[START_REF] Cazzani | Isogeometric analysis of plane curved beams[END_REF]. It is also worth to mention that this concept was used by Aristodemo to propose a highly efficient two-dimensional finite element for plane stress and plane strain problems, see[START_REF] Aristodemo | A high-continuity finite element model for two-dimensional elastic problems[END_REF].

For the interested scholars, a more complete sketch of Casciaro's scheme along with correct formulas for optimal values of the parameters α and β is reported in[START_REF] Turco | Stepwise analysis of pantographic beams subjected to impulsive loads[END_REF].

Acknowledgements

Emilio Turco gratefully acknowledges the support of the University of Sassari (Fondo di Ateneo per la ricerca 2020).

Cantilever o-beam subjected to a torque at its free end

In this subsection, we study the case of a cantilever origami beam subjected to a time-varying torque produced by two opposite and equal-in-magnitude time-varying forces acting in the Z direction, in the plane normal to the X-axis, acting on the opposite nodes of the free end. Such forces producing a torque are varying in time according to the law shown in Fig. 6 and reach, at their peak, a magnitude equal 10 5 .

The time evolution of the system for the case analysed in this subsection, computed using a time integration step Δt = 0.001, is summarized in Fig. 10.

In Fig. 10a, the time evolution of the displacement components u, v, and w along the axes X, Y , and Z, respectively, are shown for the first node of the free end, i.e. that having coordinates X = Y = Z = 0 (see Fig. 5). Figure 10b shows the evolution of the rate of the displacement components plotted in Fig. 10a, namely the evolution of the corresponding velocities, while Fig. 10c shows a stroboscopic picture of the o-beam deforming in time.

Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Springer Nature or its licensor (e.g. a society or other partner) holds exclusive rights to this article under a publishing agreement with the author(s) or other rightsholder(s); author self-archiving of the accepted manuscript version of this article is solely governed by the terms of such publishing agreement and applicable law.


Gilles Dowek
email: gilles.dowek@polytechnique.fr

Skolemization in Simple Type Theory: the Logical and the Theoretical Points of View

 a new way to prove analogs of the Miller theorem for different, but equivalent, formulations of Simple Type Theory.

In this paper, that does not contain new technical results, I try to show that the history of the skolemization problem and of its various solutions is an illustration of a tension between two points of view on Simple Type Theory: the logical and the theoretical points of view.

Skolemization 1.The Skolem theorem

Let T be a theory in first-order predicate logic containing an axiom of the form ∀x 1 ...∀x n ∃y A and T ′ be the theory obtained by replacing this axiom by

where f is a function symbol not used in T . Then, the Skolem theorem asserts that the theory T ′ is a conservative extension of T and, in particular, that one theory is contradictory if and only if the other is.

Extending Skolem theorem to Simple Type Theory

The Skolem theorem plays a key role in automated theorem proving because it permits to eliminate quantifier alternation in the proposition to be proved, or refuted, and this alternation is often delicate to manage.

in 1971, seeking for a generalization of the Resolution method to Simple Type Theory [START_REF] Peter | Resolution in Type Theory[END_REF], Peter Andrews has proposed the problem of finding an analog of the Skolem theorem for Simple Type Theory.

Following the Skolem theorem for first-order predicate logic, we can try to replace an axiom of the form ∀x 1 ...∀x n ∃y A where x 1 , ..., x n and y are variables of type T 1 , ..., T n and U , by

∀x 1 ...∀x n (((f x 1 ... x n )/y)A)
where f is now a new constant of type T 1 → ... → T n → U . Unfortunately, the theory T ′ obtained this way is not always a conservative extension of the theory T . For instance, it is not possible to prove the proposition ∃g ∀x (P x (g x)) from the axiom ∀x ∃y (P x y), because, as proved again by Peter Andrews in 1972 [START_REF] Peter | General Models, Descriptions, and Choice in Type Theory[END_REF], the axiom of choice (∀x ∃y (P x y)) ⇒ (∃g ∀x (P x (g x)))

is not provable in Simple Type Theory. But this proposition is obviously provable from the skolemized form of this axiom: ∀x (P x (f x)). Thus, although this naive skolemization can be used in Simple Type Theory extended with the axiom of choice, it cannot be used in the usual formulation of Simple Type Theory, without the axiom of choice.

A more restricted form of skolemization has been proposed in 1983 by Dale Miller. In the Miller theorem, the Skolem symbol f is not only given the type T 1 → ... → T n → U but also the arity ⟨T 1 , ..., T n , U ⟩ and, unlike the usual symbols of a functional type, the Skolem symbols are not terms per se. To form a term with a Skolem symbol f of arity ⟨T 1 , ..., T n , U ⟩, it is necessary to apply it to terms t 1 , ..., t n of type T 1 , ..., T n , called the necessary arguments of the symbol f . This reflects the intuition that, in a model of the negation of the axiom of choice, we may have for each n-uple t 1 , ..., t n the object (f t 1 ... t n ), without having the function f itself as an object.

As it is not a term, the symbol f cannot be used as a witness to prove the proposition ∃g ∀x (P x (g x)). However, this restriction is not sufficient, because, although the symbol f cannot be used as a witness, the term λz (f z), can, yielding, after normalization, the same result: ∀x (P x (f x)). This explains why Miller has introduced a second restriction: that the variables free in the necessary arguments of a Skolem symbol cannot be bound by a λ-abstraction, higher in a term substituted for a variable in a quantifier rule.

With these two restrictions, Miller has been able to give a syntactic proof of an analogous of the Skolem theorem: any proof using the axiom

∀x 1 ...∀x n (((f x 1 ... x n )/y)A)
and whose conclusion does not use the symbol f can be transformed into a proof using the axiom ∀x 1 ...∀x n ∃y A

The logical and the theoretical point of view

Simple Type Theory is otherwise known as higher-order logic. This duality of names reveals that this formalism is sometimes seen as a theory and sometimes as a logic. As a theory, it should be compared to other theories, such as arithmetic or set theory. As a logic, to other logics, and in particular to first-order predicate logic, that, historically, is a restriction of it.

From the logical point of view, it is natural to try to extend the theorems and algorithms known for first-order predicate logic, such as the Gödel completeness theorem, the Skolem theorem, the Resolution method, ... to higher-order predicate logic. But, from the theoretical point of view, it is more natural to try to express Simple Type Theory as an axiomatic theory in first-order predicate logic and apply these theorems and algorithms to this particular theory.

This tension between the logical and the theoretical point of view is illustrated in a 1968 discussion between J. Alan Robinson and Martin Davis [START_REF] Robinson | New Directions In Mechanical Theorem Proving[END_REF]5]. Robinson calls first-order predicate logic: "the restricted predicate calculus" and higher-order logic: "the full predicate calculus", while Davis, replying that it is a simple matter to express Simple Type Theory as a theory in first-order predicate logic, calls such a theory expressed in first-order predicate logic a "theory with standard formulation".

The logical point of view has dominated the history of Simple Type Theory: Henkin models, Higher-order resolution, ... have been designed specifically for Simple Type Theory and not for a class of theories formulated in first-order predicate logic, among which Simple Type Theory is one instance. It is only recently that presentations of Simple Type Theory in first-order predicate logic have shed a new light on Henkin models, as anticipated by Martin Davis [5], on higher-order unification, on proof search algorithms, on cut elimination theorems, on functional interpretations of constructive proofs, ... and on the Miller theorem.

Indeed, as shown in [START_REF] Dowek | HOL-lambda-sigma: an Intentional First-order Expression of Higher-order Logic[END_REF], expressing Simple Type Theory as a theory in firstorder predicate logic allows to apply the Skolem theorem and this way to reconstruct the Miller theorem.

3 Simple Type Theory as a theory in first-order predicate logic

We shall express Simple Type Theory as a theory in many-sorted first-order predicate logic with equality. If a single-sorted formulation were needed, the usual relativization method could be applied.

Sorts

As, in Simple Type Theory, functions and predicates are objects, the sorts of the theory are not only the base types ι and o of Simple Type Theory but all its types:

ι, o, ι → ι, ι → o, o → o, ...
As usual in many-sorted first-order predicate logic, terms are assigned a sort and function and predicate symbols are assigned a tuple of sort called an arity or a rank. If f is a function symbol of arity ⟨T 1 , ..., T n , U ⟩ and t 1 , ..., t n are terms of sort T 1 , ..., T n , then f (t 1 , ..., t n ) is a term of sort U and if P is a predicate symbol of arity ⟨T 1 , ..., T n ⟩ and t 1 , ..., t n are terms of sort T 1 , ..., T n , then P (t 1 , ..., t n ) is an atomic proposition.

Symbols

It is well-known that making a predicate P an object requires to introduce a copula ∈ and write a ∈ P what was previously written P (a). In the same way making a function f an object requires to introduce an application symbol α and write α(f, a) what was previously written f (a).

In Simple Type Theory, we need a function symbol α T,U of arity ⟨T → U, T, U ⟩ for each pair of sorts T, U and a single predicate symbol ε of arity ⟨o⟩ to promote a term t of sort o to a proposition ε(t).

For instance, if P is a term of sort T → o and t a term of sort T , then the proposition usually written P (t) is not written t ∈ P , like in set theory, but ε(α T,o (P, t)) where P is first applied to t using the function symbol α T,o to build a term of sort o, that is then promoted to a proposition using the predicate symbol ε. In the same way, the proposition usually written

∀P (P ⇒ P ) is written ∀p (ε(p) ⇒ ε(p)).
Then, we need symbols to construct terms expressing functions and predicates. Introducing the binding symbol λ is not possible in first-order predicate logic and we have to use a first-order encoding of λ-calculus. A simple solution is to use the combinators S and K. Thus, for each triple of sorts T, U, V , we introduce a constant S T,U,V of sort (T → U → V ) → (T → U ) → T → V and for each pair of sorts T, U , we introduce a constant K T,U of sort T → U → T . Finally, we need similar combinators to build terms of sort o. Thus, we introduce constants =T of sort T → T → o, ⊤ and ⊥ of sorts o, ¬ of sort o → o, ∧, ∨ and ⇒ of sort o → o → o, ∀T and ∃T of sort (T → o) → o. Of course, some of these symbols are redundant and could be defined from others using de Morgan's law. It is also possible to define all these symbols from equality =T , following the idea of Leon Henkin and Peter Andrews [START_REF] Henkin | A theory of propositional types[END_REF][START_REF] Peter | A Reduction of the Axioms of the Theory of Propositional Types[END_REF].

Indices may be omitted when they can be reconstructed from the context.

Axioms

Finally, we need axioms expressing the meaning of these symbols. Besides the axioms of equality, we take the axioms

∀x∀y∀z (α(α(α(S, x), y), z) = α(α(x, z), α(y, z))) ∀x∀y (α(α(K, x), y) = x) ∀x∀y (ε(α(α( =, x), y)) ⇔ (x = y)) ε( ⊤) ⇔ ⊤ ε( ⊥) ⇔ ⊥ ∀x (ε(α( ¬, x)) ⇔ ¬ε(x)) ∀x∀y (ε(α(α( ∧, x), y)) ⇔ (ε(x) ∧ ε(y))) ∀x∀y (ε(α(α( ∨, x), y)) ⇔ (ε(x) ∨ ε(y))) ∀x∀y (ε(α(α( ⇒, x), y)) ⇔ (ε(x) ⇒ ε(y))) ∀x (ε(α( ∀, x)) ⇔ ∀y ε(α(x, y))) ∀x (ε(α( ∃, x)) ⇔ ∃y ε(α(x, y)))
To these axioms, we may add, as usual, the extensionality axioms, the axiom of infinity, the description axiom and the axiom of choice.

To recall the choice we have made to express terms with combinators, we call this theory HOL-SK.

Properties

An easy induction on the structure of t permits to prove that for each term t there exists a term u such that the proposition α(u, x) = t is provable. The term u is often written λx t. In the same way, an easy induction on the structure of P permits to prove that for each proposition P there exists a term u such that the proposition

ε(u) ⇔ P is provable.
This way, we can translate all the terms of the usual formulation of Simple Type Theory with λ-calculus to terms of HOL-SK and prove that if t is a term of sort o and t ′ its translation, then t is provable in the usual formulation of Simple Type theory if and only if ε(t ′ ) is provable in HOL-SK. However, this theorem requires that the extensionality axioms are added to both theories, because the combinators S and K do not simulate λ-calculus exactly, but only up to extensionality.

Skolemization

As a consequence of the Skolem theorem, we get that if T is a theory containing an axiom of the form ∀x 1 ...∀x n ∃y A where x 1 , ..., x n and y are variables of sorts T 1 , ..., T n and U , and T ′ is the theory obtained by replacing this axiom by

∀x 1 ...∀x n ((f (x 1 , ..., x n )/y)A)
then HOL-SK ∪ T ′ is a conservative extension of HOL-SK ∪ T . The symbol f is a function symbol of arity ⟨T 1 , ..., T n , U ⟩ and not a constant of sort T 1 → ... → T n → U , first because the Skolem theorem for first-order predicate logic introduces a function symbol and not a constant and then because it ignores the internal structure of sorts. Thus, the symbol f alone is not a term, and the term obtained by applying the symbol f to the terms t 1 , ..., t n is the term f (t 1 , ..., t n ) and not the ill-formed term α(...α(f, t 1 ), ..., t n ). In contrast, when the sort U has the form V → W , the term f (t 1 , ..., t n ) can be further applied to a term t n+1 using the application symbol: α(f (t 1 , ..., t n ), t n+1 ).

In short, the Skolem symbols are not at the level of the symbols S or K, but at the level of the symbols α T,U .

We prove this way an analogous of the Miller theorem for HOL-SK. For this formulation Miller's second condition vanishes as there is no binder λ. All that remains is the first condition: the fact that Skolem symbols must be applied. Notice however that the term λx t cannot be defined for all terms t containing the symbol f , but only those that do not have an occurrence of the variable x in an argument of the symbol f . An advantage of expressing Simple Type Theory as a theory in first-order predicate logic is that the proof of the Miller theorem is simplified as it is then proved as a consequence of the Skolem theorem. Moreover this shows that arities are not a feature of the Skolem symbols only, but that all function symbols have arities, in particular the symbols α T,U and the Skolem symbols.

Miller's second condition

We may wonder if it is possible to go one step further and reconstruct Miller's second condition as a consequence of the Skolem theorem for first-order predicate logic. As we shall see, this is possible but this requires to use a more precise first-order encoding of the λ-calculus: the λ-calculus with nameless dummies introduced by Nicolaas de Bruijn [START_REF] Nicolaas | Lambda Calculus Notation with Nameless Dummies, a Tool for Automatic Formula Manipulation, with Application to the Church-Rosser Theorem[END_REF].

De Bruijn indices

In a λ-term, we may add, to each occurrence of a bound variable, a natural number expressing the number of abstractions separating the occurrence from its binder. For instance, adding indices to the term λxλy (x λz (y x z))

yields λxλy (x 2 λz (y 2 x 3 z 1 ))
Indeed, the index of the occurrence of the variable z is 1 because the binder λz is just one level above in the term, while the index of the second occurrence of the variable x is 3 as the the binder λx is three levels above. Once the indices are added this way, the names of the bound variables are immaterial and can be dropped. We thus get the term λλ ( 2 λ ( 2 3 1 ))

Notice that the term λxλy (x λz (y x z)) is closed but that its subterm λy (x λz (y x z)) contains a free variable x. In the same way, in the term λλ ( 2 λ ( 2 3 1 )) all the de Bruijn indices refer to a binder, but its subterm λ ( 2 λ ( 2 3 1 )) contain two "free indices", that exceeded the number of binders above them, and correspond to the former occurrences of the variable x.

The λ-calculus with de Bruijn indices is a first-order language as, once variables names have been dropped, the symbol λ is not a binder anymore. This language is formed with an infinite number of constants 1 , 2 , 3 , ..., a binary function symbol α for application and a unary function symbol λ. Closed terms, such as λ ( 2 λ ( 2 3 1 )) or 1 , may contain free indices, exceeding the number of binders above them.

Sorts

When type-checking an open λ-term à la Church, i.e. with explicitly typed bound variables, for instance λy ((ι→ι)→ι)→ι→ι (x λz ι (y x z))

it is necessary to have a context defining the type of the free variables such as x. Indeed, if x is assigned the type (ι → ι) → ι, then the term is well-typed, but not if it is assigned, for instance, the type ι.

In the same way, when typing a term that may contains free de Bruijn indices, even if this term is closed (i.e. does not contain named variables), we need a context defining the types of the de Bruijn indices exceeding the number of binders above them. This context is a finite list of types: the type of the indices exceeding by 1 the number of binders above them, that of the indices exceeding by 2 the number of binders above them, ... For instance the term

λ ((ι→ι)→ι)→ι→ι ( 2 λ ι ( 2 3 1 )) has type (((ι → ι) → ι) → ι → ι) → ι in the context [(ι → ι) → ι].
Thus, even if it is closed, a term t can be assigned a type T , only relatively to a context Γ. In other words, it can be absolutely assigned an ordered pair formed with a context Γ and a type T . We write such a pair Γ ⊢ T .

When we express Simple Type Theory as a first-order theory, using the λcalculus with de Bruijn indices as a first-order encoding of the λ-calculus, the sort are not just the simple types, like in HOL-SK, but such pairs Γ ⊢ T formed with a list Γ of simple types and a simple type T .

In this formulation, the quantified variables, in contrast to the λ-abstracted ones, are not replaced by de Bruijn indices but, as in all theories expressed in first-order predicate logic, they are kept as standard named variables and they are assigned a sort of the form ⊢ T , where the context is the empty list.

The last step of the construction of this theory would be to introduce explicit substitutions, hence its name HOL-λσ. I do not want to go into these details here and the interested reader can refer to [START_REF] Dowek | HOL-lambda-sigma: an Intentional First-order Expression of Higher-order Logic[END_REF]. But, I want to insist on two points. First, the fact that this theory is intentionally equivalent to the usual presentation of Simple Type Theory with λ-calculus, i.e. even if the extensionality axioms are not assumed. Second, that, as explained above, the sorts are pairs Γ ⊢ T , i.e. that they contain scoping information. In particular terms of a sort ⊢ T , with an empty context, are de Bruijn-closed, i.e. they do not contain indices that may be bound higher in the term by a λ-abstraction.

Skolemization

As a consequence of the Skolem theorem, we get that if T is a theory containing an axiom of the form ∀x 1 ...∀x n ∃y A where x 1 , ..., x n and y are variables of sort ⊢ T 1 , ..., ⊢ T n and ⊢ U , and T ′ is the theory obtained by replacing this axiom by

∀x 1 ...∀x n ((f (x 1 , ..., x n )/y)A)
then HOL-λσ ∪ T ′ is a conservative extension of HOL-λσ ∪ T . The symbol f is a function symbol of arity ⟨⊢ T 1 , ..., ⊢ T n , ⊢ U ⟩. Thus, not only f is a function symbol, but a function symbol whose arguments must not contain indices that may be bound by a λ-abstraction higher in a term. We get this way exactly Miller's conditions. The first is rephrased as the fact that f is a function symbol and the second as the fact that the arguments of these function symbol must have a sort with an empty context. We obtain this way an alternative proof of the Miller theorem for HOL-λσ, as consequence of the Skolem theorem.

A difficult point in the Miller theorem is the discrepancy between the general formation rules for the terms and the propositions and the more restricted ones for the terms substituted for variables in quantifier rules. As we have seen, the variables free in the necessary arguments of the Skolem symbols cannot be bound in a term substituted for a variable. But this restriction does not apply to the terms and propositions in general, because, if it did, the skolemized axiom itself would not be well-formed, as the arguments of the Skolem symbol in this skolemized axiom are universally bound variables.

In HOL-λσ the situation is slightly different. The formation rules for the terms substituted for variables in quantifier rules are the same as those of the other terms. But the λ-bound variables and the quantified variables are treated differently. As we have seen, the λ-bound variables are replaced by de Bruijn indices, but the quantified variables are kept as named variables, as in all theories expressed in first-order predicate logic. Thus the λ-bound variables cannot appear in the arguments of a Skolem symbol but the quantified variables can and, this way, the skolemized axiom is well-formed.

From the theoretical point of view

With this example of skolemization, we have seen two advantages of the theoretical point of view on Simple Type Theory. First the proofs of the theorems are simplified, because we can take advantage of having already proved similar theorems for first-order predicate logic. Then some difficult points of these theorems are explained.

We can mention several other theorems and algorithms that have been simplified and explained this way. The Henkin completeness theorem, already mentioned in 1968 by Martin Davis, independence results, in particular of the extensionality axioms, cut elimination theorems, both model-based ones and reduction-based ones, proof search algorithms and the functional interpretation of constructive proofs.

This succession of points of view on Simple Type Theory is itself an example of a common back and forth movement in the development of science: first, new objects and new results are discovered, breaking with the old framework, they are supposed not to fit in. Then, after a possible evolution of the general framework, these new objects are integrated back. We may try to understand what kind of evolution of first-order predicate logic, the integration of Simple Type Theory requires or suggests.

The first evolution is the shift from the single-sorted first-order predicate logic to the many-sorted one. Although it is always possible to relativize a many-sorted theory to a single-sorted one, the natural framework to express Simple Type Theory is many-sorted first-order predicate logic. It is interesting to see that early papers on many-sorted first-order predicate logic, e.g. [START_REF] Hao | Logic of Many-sorted Theories[END_REF], already motivate the introduction of many-sorted first-order predicate logic by the will to express Simple Type Theory in it.

A second evolution is to take into account the possibility to consider terms and propositions up to reduction. The axioms of HOL-SK (and those of HOLλσ) can easily be transformed into rewrite rules: α(α(α(S, x), y), z) -→ α(α(x, z), α(y, z)) α(α(K, x), y) -→ x ε(α(α( =, x), y)) -→ x = y ε( ⊤) -→ ⊤ ε( ⊥) -→ ⊥ ε(α( ¬, x)) -→ ¬ε(x) etc. and identifying equivalent propositions is more natural than keeping these axioms as such, exactly like in the usual formulation of Simple Type Theory, β-normalizing terms and propositions is more natural than keeping β-conversion as an axiom. This idea has lead to the development of Deduction modulo, that again was initially motivated by the will to express Simple Type Theory in it.

A third evolution is the introduction of binders in first-order predicate logic. Although this has been a hot topic recently, it is fair to say that we have not yet a completely satisfactory extension of first-order predicate logic with such binders.

To conclude, I want to mention the influence of Peter Andrews' work, course and book [START_REF] Peter | An Introduction to Mathematical Logic and Type Theory: to Truth through Proofs[END_REF], on the emergence of the theoretical point of view on Simple Type Theory (although I do not attribute any point of view to anyone except myself). First, it is striking that in his course, Peter Andrews compared Simple Type Theory not only to first-order predicate logic but also to set theory in particular that he insisted on the fact that both theories are restriction of the naive, inconsistent, set theory with full comprehension. More technically, the alternative characterization of Henkin models, given by Peter Andrews [START_REF] Peter | General Models, Descriptions, and Choice in Type Theory[END_REF] showed the way out of Henkin mysterious conditions that all λ-terms must have a denotation to a standard condition that the axioms S and K must be valid and thus to the idea that Henkin models were just models of some theory expressed in first-order predicate logic.

Acknowledgements

The author wants to thank Claude Kirchner and Dale Miller for helpful comments on a previous version of this paper.