Simulation of membrane deployment accounting for the nonlinear crease effect based on absolute nodal coordinate formulation

With the advantages of lightweight and high stowing efficiency, the origami membrane has a broad application prospect in the space structure. As an improvement of the existing deployment simulations, the behavior of crease is considered and both static and dynamic simulation are conducted in this paper. The origami membrane is modeled as a flexible multibody system and its deployment is simulated based on the absolute nodal coordinate formulation. The nonlinear effect of the crease to resist deployment is integrated into the multibody system via the virtual work principle. The facet is modeled by the thin shell element and the crease is constructed as a nonlinear torsional spring with specific position and gradient constraints. The behavior of the crease is parameterized based on the experimental data. The modeling approach is verified by the numerical/experimental reference of the Z-folding structure in the existing literature. For the multi-crease origami membrane model considering the crease effect, a form-finding method to obtain the initial configuration of the flexible origami is proposed based on the principle of minimum potential energy. The deployment of the classical Miura-ori unit membrane structure is analyzed. The magnitude of driving force is estimated based on the force analysis of rigid deployment. Based on the simulation result of flexible deployment, the driving force is planned to achieve a stable deployment with a constant increase speed in the deployment ratio. This work is expected to provide a reference for the design of space membrane structure, the estimation of driving force, and optimizing the deployment strategy of an origami membrane.

and the crease is constructed as a nonlinear torsional spring with specific position and gradient constraints. The behavior of the crease is parameterized based on the experimental data. The modeling approach is verified by the numerical/experimental reference of the Z-folding structure in the existing literature. For the multi-crease origami membrane model considering the crease effect, a form-finding method to obtain the initial configuration of the flexible origami is proposed based on the principle of minimum potential energy. The deployment of the classical Miura-ori unit membrane structure is analyzed. The magnitude of driving force is estimated based on the force analysis of rigid deployment. Based on the simulation result of flexible deployment, the driving force is planned to achieve a stable deployment with a constant increase speed in the deployment ratio. This work is expected to provide a reference for the design of space membrane structure, the estimation of driving force, and optimizing the deployment strategy of an origami membrane.

Introduction
Space tasks such as the deep space exploration and high-resolution Earth imaging propose the requirement of large solar sail and antenna [5,16]. This enlarges the conflict between the volume/weight of spacecraft and the limited payload capability of launch vehicle. With the advantages of lightweight and high stowing efficiency, the deployable space membrane has become the development trend of the large-sized structure of the space equipment. Whether the folded membrane can be successfully deployed on orbit is a key factor that determines the success or failure of the space mission [23]. Researchers have paid much attention to the deployment of the folded/origami membrane structure.
When a membrane is folded along a line, a residual crease mark can be observed. Creases are induced in the folding process of membrane. The residual crease largely alters the shape and behavior of the membrane [2,23]. To date, the analysis method of the membrane deployment can be classified into two types, i.e., the analytical/semi-analytical method and the finite element method. Okuizumi et al. [19] simplified the membrane structure to a mass-spring system. Furuya et al. [10] treated the uniformly loaded Z-folding membrane structure as a one-dimensional structure and analyzed it based on the curved elastic beam theory. Xia et al. [33,34], Satou et al. [24], and Dharmadasa and Jiménez [4] also proposed several beam-theory-based models for the Z-folding structure. Schenk et al. [11,15,25,36] used the so-called bar-hinge method that analyzes the origami structure as an equivalent pin-jointed truss frame. This method models crease as a bar element with a spring hinge. The behavior of the origami is represented by the rotation around the hinge axis. For the dynamic modeling, the inertia is distributed on the vertices. The above methods are classified as the analytical/semi-analytical methods. These methods lead to a significant simplification of the origami structure with a sacrifice in the detailed description of the facet. Specifically, the inertia, bending behavior, and the large deformation of the facet cannot be well represented. With respect to the finite element (FE) method, the facets in the origami are modeled with the beam element [6], the thin shell/membrane element [2,3,20,21], and solid element [7,23,24,32]. Specifically, beam FE models can hardly describe multi-crease origami structure besides the Z-folding structure. In the existing literature, the solid FE models are usually used as an initial reference for the verification or used for the parameter determination of other more computationally efficient methods [32], and the analyzed structure is usually the simple Z-folding structure [28]. In summary, the thin shell/membrane FE model draws the most attention among the adopted FEs to numerically simulate the deployment of an origami membrane. Compared with the beam and solid FE, the thin shell/membrane FE can simulate the deployment of a relatively sophisticated structure with a moderately low computational cost.
The existing thin shell/membrane FE models consist of two groups depending on how to consider the effect of crease. Cai et al. [1] simulated the crease via the effective modulus method. The crease region is implicitly represented from the view of equivalent material behavior, while the other method is inspired by the bar-hinge method and the crease is explicitly modeled as a nonlinear equivalent torsional spring connecting the adjacent facets, which are modeled by thin shell/membrane element. Most of the thin shell/membrane FE models are constructed using the conventional FE and the analyses are performed in the commercial FE software ABAQUS [2,5,28]. The thin shell/membrane FE models are not limited to simple Zfolding structures and extend to complex folding structures including the Miura-ori folding structures [28] and some solar sail [3]. However, most of the above researches [1,2,5,28] did not present dynamic analysis. A relatively dense mesh is used, which brings a demerit in the computational efficiency.
Compared with the conventional FE, the absolute nodal coordinate formulation (ANCF) exhibits an outstanding capability in the analysis of problem with large deformation and large reference motion [26], which makes it a good vehicle for the simulation of membrane deployment. Using the components of the deformation gradient as the nodal coordinates, the ANCF method imposes no restrictions on the magnitude of the deformation or rotation within the FE. With a global-framebased kinematic representation, the ANCF method is straightforward to implement in a multibody system dynamic analysis and the corresponding constant mass matrix results in a concise equation of motion (EoM). Due to its consistency with strain and stress measures from general continuum mechanics, the ANCF method has been widely applied in various nonlinear static/dynamic problems [12,13,30,31]. Luo et al. [17] and Yuan et al. [37,38] analyzed a spinning deployable solar sail, a leaf-in origami membrane, and a leaf-out origami membrane, respectively, using the thin ANCF shell/membrane element. These researches did not consider the nonlinear torque provided by the crease to resist the deployment. This property of crease is called as "anti-deployment" effect in this paper. Later, Tao and Li [29] presented some derivations to integrate the crease effect into the ANCF thin shell model via the Lagrange equation. But the corresponding numerical simulations are not presented. Additionally, the stiffness of the connecting rotational spring is constant and independent of the crease deployment angle. Satou and Furuya [23] pointed out that the crease plays an important role in the deployment of the membrane. The crease exhibits a highly nonlinear behavior associated with the plastic deformation near the crease [2,9,14]. In summary, the ANCF method is a promising approach to simulate the deployment of the folded membrane structure. However, analyses accounting for the nonlinear crease effect based on the ANCF method are still rare.
The objective of this paper is to develop an approach to simulate the deployment of the origami membrane considering the nonlinear anti-deployment effect of crease based on the absolute nodal coordinate formulation. To achieve this, the nonlinear anti-deployment property of the crease is parameterized based on the experiment in [28]. The crease effect is integrated to a multibody system modeled with the ANCF thin shell element via the virtual work principle. The approach is verified by the experimental/numerical reference of the Z-folding structure in [28]. The approach extends to the analysis of a Miura-ori membrane structure and the driving force is planned based on the simulation result. The research is to explore an approach to simulate the membrane deployment with moderate accuracy by exploiting the merit of ANCF method in the nonlinear analysis and by considering the anti-deployment effect of the crease. It is expected to provide suggestions for the design of space membrane structure, the estimation of driving force, and optimizing the deployment strategy of an origami membrane.
The remainder of the paper is organized as follows: In Section 2, the method to integrate the nonlinear crease effect into the origami membrane multibody system model and the calculation of the crease antideployment moment are presented. Section 3 shows how to model the Miura-ori origami unit including the initial configuration finding procedure and the special constraints in the flexible multibody system model. Section 4 presents quasi-static and dynamic analysis of the Z-folding and Miura-ori membrane deployment. Based on the simulation result, the driving force of the membrane deployment is planned to achieve a stable deployment and the corresponding simulation is conducted. Section 5 offers the authors' conclusions.

Multibody system modeling accounting for the anti-deployment effect of crease
This section presents the multibody system modeling method of the origami membrane structure. The nonlinear behavior of crease in deployment is parameterized and integrated into the multibody system model. Then, the specific modeling method of the flexible multibody system based on the ANCF thin shell element is proposed.

Modeling of the nonlinear behavior of crease
The cross section of the membrane in the folding and deployment process is shown in Fig. 1. To save the page space, only half symmetric membrane is shown in the figure. As shown in Step 1, a crease is created by bending the membrane along a line to a tight status so that the distance between the facet and the symmetric plane decreases to a small value d, which is usually comparable with the membrane thickness. During the folding process, permanent deformation appears in the narrow region connecting two facets. We call this region as crease. Once the pressing force is removed, the folded membrane will bounce to an equilibrium state with a unilateral deployment angle ϕ 0 and most regions in the facet keep flat and straight as shown in Step 2. For the sake of compact layout of Fig. 1, the diagram of Step 3 is vertically placed. Under the effect of dragging force, the folded membrane is gradually deployed. The bilateral and unilateral deployment angle is denoted as θ and ϕ, respectively, as shown in Step 3 Fig. 1.
In the deployment process, the crease exhibits a nonlinear behavior with increase in the dragging force. According to [4], the effect of a crease can be characterized as a virtual torsional spring. The nonlinear relationship between the anti-deployment moment M and the unilateral-deployment angle ϕ based on experiment is employed in this paper. In the simulation, the crease is characterized as the virtual torsion spring, whose behavior is parameterized using the experimentbased curve M-ϕ in [28]. To integrate the nonlinear behavior of a crease into the multibody system model, the virtual work principle in the case of the dynamic modeling [22,27] is used as shown in Eq. (1).
where δW i is the virtual work of inertia forces that only appears in the case of dynamic analysis, δW s is the virtual work done by the elastic force associated with flexible deformation that appears in a flexible model and vanishes in a rigid multibody system model, δW c is the virtual work done by the crease to resist the deployment of the membrane, and δW e is the virtual work done by the external force to deploy the membrane. Recall that the bilateral deployment angle θ = 2ϕ, the virtual work δW c is written as where M(ϕ) is obtained from the experiment curve Mϕ in [28]. Based on Eqs. (1) and (2), the equation of motion considering the crease effect can be obtained.

Flexible multibody system model based on ANCF thin shell element
Similar to reference [28], the thin shell element is used in the flexible multibody system model to avoid numerical instability. The thin plate/shell element proposed by Dufva and Shabana [8] is used in this work to model the origami structure. As shown in the left subfigure of As shown in the right subfigure of Fig. 2, a node couple i and j on the crease is enlarged, where the coupled nodes are moved apart to illustrate the corresponding constraint. Normally, the crease region is very narrow and Shen et al. [28] set the distance between the coupled nodes d c as the thickness of the membrane. In this work, d c = 0 is used. The constraints on the node couple i and j are where r k , k = i, j are the position vectors of node i and j, the subscript following a comma denotes partial differentiation with respect to that variable, and α, β = x, y indicate the material coordinate parallel to the crease in the corresponding element in each connected facet. Depending on the element coordinate used in each facet, α is not necessarily same with β in the constraints Eq. (3). The anti-deployment moment of the crease is provided by the virtual torsional spring placed on each node couple. Recalling Eq. (2), the corresponding generalized force Q c is obtained as the derivative of the virtual work W c with respect to the corresponding vector of system nodal coordinate e. With respect to the multibody system, Q c is treated as the internal force and δW c = −Q T c δe.
where the moment associated with the deployment angle θ is obtained from the experiment-based M-ϕ curve in [28]. The deployment angle is related to the angle between the normal surface vectors n i , n j at the two crease couple nodes i and j, i.e., γ =< n i , n j >. By defining the normal surface vectors of involved facets in the same direction before folding, the deployment angle is obtained as θ = π −γ , and further related with the nodal coordinates at node i and j by where the unit normal vectors at crease couple nodes i and j arê where | · | indicates the norm of the corresponding vector. Based on Eqs. (5) and (6), the derivative of θ with respect to the system nodal coordinate e is where the intermediate derivative in Eq. (7) is where the derivatives of the normal vectors and their magnitude on the two couple nodes i and j are where the subscript in the second equation in Eq. (9) indicates the dimension of the matrix, kst = (k − 1) × n cpn + 3, n cpn is the number of coordinate per node, and n f is the total number of system degrees of freedom (DOFs). The column sequence of the intermediate matrix t k corresponds to those of the gradient components of node k. The sequence of the 6 columns in t k in the 3×n f matrix ∂n k /∂e is kst +1 to kst +6. The generalized force Q c can be obtained, based on Eq. (4-9).
In the differential equation solver, the tangential stiffness matrix of Q c needs be obtained and is calculated as The second derivative of θ with respect to e is where the intermediate derivatives are where k, l = i, j, k = l, ∂n k (m)/∂e and ∂ 2nk (m)/∂e 2 are the first and second derivative of mth component of the unit normal vectors with k = i, j and m = 1, 2, 3. The first derivative ∂n k (m)/∂e can be obtained from where the first derivative of mth component of normal vector ∂n k (m)/∂e can be obtained from Eq. (9). The second derivative of mth component of normal vector is denoted as ∂ 2 n k (m)/∂e 2 = A k m . The dimension of A k m is n f × n f and the only nonzero terms in A k m are The stiffness matrix of Q c is obtained based on Eq. (10)(11)(12)(13)(14). Considering the virtual displacement is arbitrary, the equation of motion is obtained based on the virtual work principle Eq. (1) where M is the mass matrix,ë is the second derivative of system nodal coordinate e with respect to time, and Q s is the elastic force. The inertia term Më disappears in the quasi-static analysis. For the ANCF FE flexible model, the mass matrix is where h is the thickness of the thin shell, S is the shape function matrix, and A 0 is the element area in the reference configuration. The elastic force Q s disappears in the case of rigid model. For the ANCF FE flexible model, the elastic force Q s is where ε is the Green-Lagrange strain in the Voigt form, and κ is the curvature vector defined as follows where the components of the strain and curvature vector are where the "0" in the subscript indicates the variable is associated with the reference configuration, the unit normal vector is defined as Eq. (6), r ,αβ , α, β = x, y are the second derivatives of r with respect to α and β. The constitutive matrix in the plane stress case is where E is the Young's modulus and ν is the Poisson's ratio. One can refer to [30,35] for details on the stiffness matrix of the elastic force ∂Q s /∂e.

Modeling of a Miura-ori origami membrane
Miura-ori is a classical origami structure, which is usually used for for small volume folding and large-area deployment of a flat membrane [18]. The four-crease unit of the Miura-ori origami structure in the flat and folding/deployment configuration are shown in Fig. 3.
As shown in the flat configuration of the structure in Fig. 3, set the intersection point as the original point o, and build the axis x i along the creases, where the cyclic index i = 1, 2, 3, 4. The direction of axis z i is defined by the cross product of direction vectors along x i and x i+1 , and direction of axis y i is determined by the right-hand rule. The global coordinate Fig. 3 Miura-ori origami unit frame x 0 y 0 z 0 coincides with the frame x 1 y 1 z 1 in the flat configuration. The angle between x 1 and x 4 , or x 2 is γ 0 . The Miura-ori structure consists one mountain fold and three valley folds represented using solid and dashed lines, respectively, in the flat configuration in Fig. 3. The folding/deployment status of Miura-ori membrane is also shown in Fig. 3. To remove the rigid motion modes in the numerical simulation, the constraints are applied as follows. The position of origin o is fixed, the position component along axis y 0 of node A and C is fixed considering the symmetry, and the position component along axis z 0 of node B is fixed. With three valley folds and one mountain fold, the Miuraori origami unit has only one mechanism DOF. Therefore, the above five constraints and one given actuator can fully determine the configuration of the Miura-ori origami unit. In the following deployment simulation, the deployment angle θ 1 associated with crease one is used as the actuator variable.
Accounting for the anti-deployment moment provided by the crease, the initial equilibrium configuration of the structure needs be calculated beforehand. To this end, the equilibrium configuration in the case of rigid folding should be obtained at the beginning, which is used as the initial configuration to iteratively determine the equilibrium configuration of the corresponding flexible structure by static analysis via Newton-Raphson iteration.

Initial equilibrium configuration in the case of rigid folding
The initial configuration in the case of rigid folding is studied in this section. Denote the angle between axis x 1 and x 0 as γ 1 , and the angle between axes x i , i = 2, 3 and −x 0 as γ i , i = 2, 3 as shown in the fold- According to Eq. (21), the norm of the normal vectors |n 1 | = |n 2 | = sinγ 0 . The angle between two normal vectors is cos <n 1 ,n 2 >=n 1 ·n 2 = cos(2γ 1 ) − cos 2 γ 1 cos 2 γ 2 sin 2 γ 0 Considering the deployment angle associated with crease two θ 2 = π − <n 1 ,n 2 > and substituting Eq. (21) to Eq. (24) Now all the angles have been related with the actuator angle θ 1 . Accounting for the anti-deployment effect of crease, the initial equilibrium configuration can be obtained according to the principle of minimum potential energy. The initial configuration must be at the beginning elastic stage of deployment. According to the experiment curve M-ϕ in [28], the anti-deployment moment is approximately linear to the deployment angle, which indicates the stiffness of the virtual torsion spring is almost constant at the beginning elastic stage. Denote the spring stiffness as k c , the corresponding potential of the whole rigid structure is written as Eq. (26) recalling θ 1 = θ 3 and θ 2 = θ 4 .
The system potential U reaches extremum at the initial equilibrium configuration. Therefore, the derivative of U with respect to θ 1 is zero, which leads to Comparing Eq. (27) and Eq. (28), the actuator angle θ 1 associated with the initial equilibrium configuration is obtained by solving the following equation (θ 1 − θ 0 )cos 2 γ 2 sinθ 2 + (θ 2 − θ 0 )cos 2 γ 1 sinθ 1 = 0 (29) Referring to [28], the angle associated with zero antideployment moment is θ 0 = 9.4 • . This zero-moment angle θ 0 is used as the initial value of θ 1 to start the Newton-Raphson iterative solver of Eq. (29) to obtain the initial equilibrium configuration of rigid folding e r 0 . The actuator angle associated with the initial equilibrium configuration is obtained as θ 1 = 11.04 • , and the corresponding θ 2 = 2.679 • . It is clear that the deployment angle of each crease in the equilibrium configuration deviates slightly from the zero-moment angle.

Flexible multibody system modeling
The flexible model is shown in Fig. 5. Position constraints r i = r j and gradient constraints are applied on the crease node couples. As shown in Fig. 5, constraints r i ,x = r j ,x are applied for the node couple (i, j) on crease x 1 and x 3 and r i ,y = r j ,y for crease x 2 and x 4 . The same constraints used in the rigid deployment analysis that remove the rigid motion modes are adopted. The position of node o is fixed, the position component of node A and C along y 0 direction is fixed, and the position component of node B along z 0 direction is fixed.
Due to the internal elastic force in the thin shell element, there is a unbalanced force in the flexible Miura-ori structure associated with the initial equilibrium configuration obtained from the rigid folding analysis. Therefore, there is a difference between the initial  Fig. 6.
In the simulation of flexible deployment, e r 0 is used as the stress-free reference configuration. Under the anti-deployment effect of crease, there is a pre-stress in the initial configuration e f 0 of flexible multibody system model.

Numerical case
In this section, the deployment of a single-crease Zfolding structure is first simulated by integrating the crease behavior depicted by the M-ϕ curve based on the experiment [28]. The validity of the flexible multi- Fig. 7 Definition of deployment error and ratio body system modeling approach is verified by the comparison on the unilateral deployment angle ϕ. Then, a typical Miura-ori origami unit is analyzed. Based on the simulation result, the driving force is planned to achieve a steady deployment. According to [28], the thickness of the membrane employed in the study is 100 µm, the Young's modulus is 4.8 GPa, and Poisson's ratio is 0.3.
The deployment error d e and ratio ρ are introduced to measure the deployment process of the folding structure. Take the Z-folding structure as an example to show the definition of d e and ρ. As shown in Fig. 7, in the limit case of θ = π , all the facets are deployed to a single plane xoy called 0 , which is defined as the deployment plane. The direction perpendicular to the deployment plane 0 is called the vertical direction. The deployment error d e is defined as the maximum distance between the profile of the origami structure and the deployment plane along the vertical direction. In the special case of Z-folding structure as shown in Fig. 7, d e is the distance between the plane 0 and the plane 1 that is parallel to 0 and contains the mountain fold. The deployment ratio is the ratio between the projection area of the folding structure A f in the deployment plane and the total area of the involved facets A flat in the flat configuration. By means of d e and ρ = A f /A flat , the deployment of the structure is described in the vertical direction and the deployment plane, respectively.

Deployment of a Z-folding structure
A single-crease Z-folding structure consisting of two 25 × 50 mm facets adopted in [28] is analyzed. The Because the membrane is subjected to a larger deformation in the longitudinal direction, a more refined division is used for the longitudinal side, which leads to n y < 2n x for this structure. In this work, 7 meshes as shown in Fig. 9 are used. The deployment of Zfolding structure is simulated in a quasi-static form. The deployment error under different load steps and the final values are shown in Fig. 9. And the unilateral deployment angle ϕ is shown in Fig. 10.
As shown in Fig. 9, models with different meshes show small difference on the deployment error. For the last force step, the difference is even comparable with the membrane thickness. By contrast, as shown in Fig. 10, an evident difference can be observed on the unilateral deployment angle ϕ obtained by models Fig. 10 Unilateral-deployment angle ϕ obtained by models with different meshes with different meshes. Observing the results associated with 10 × 10, 10 × 12, and 10 × 20, it is found that 10 divisions in the lateral direction are enough to obtain a converged solution. With 10 divisions in lateral direction, 8, 10, 12, 16, and 20 divisions in the longitudinal direction are compared. Comparing the result at the last load step, the relative error on ϕ between 16 × 10 and 20 × 10 decreases to 0.6%. The result obtained by 20 × 10 is considered as the converged result and is compared with the experiment result, and the numerical result given by the thin shell element S4R5 model in the commercial FE software ABAQUS used in [28] as shown in Fig. 11.
As shown in Fig. 11, the ANCF model agrees well with the ABAQUS shell element model [28] on the deployment angle. The validity of the modeling approach is verified by the deployment of the Z-folding Fig. 11 Unilateral deployment angle comparison structure. At the final configuration, the difference between the results given by the ANCF thin shell model and the experiment on the deployment angle ϕ is around 7%. The deviation is possibly because of the difference in the boundary condition. In the experiment, the edge sides are clamped to a stiff board. To some extent, the rotation of the edge nodes are constrained. By contrast, in the numerical analysis, the rotation of the edge nodes are not constrained.

Deployment of a Miura-ori origami structure
A unit Miura-ori origami structure with four 50×50 mm diamond facets adopted in [28] is studied in this section. The obtuse angle of the diamond facet is γ 0 = 104 • . Each facet in the flexible Miura-ori folding structure is meshed by n m × n m thin shell elements as shown in Fig. 5, where n m is the division number for each side. To deploy the folding structure, four dragging forces on the corner nodes E, F, G, and H along the diagonal lines of each facet in the ideal final deployment plane are applied as shown in Fig. 5. Due to the symmetric condition, only F 1 and F 4 can be considered as the independent actuator forces. Considering the balance in the force along x 0 direction, F 4 = F 1 tan(γ 0 /2). As a result, F 1 is used as the only independent actuator force. Accordingly, the number of actuator force is equal to the number of mechanism DOF of the fourcrease Miura-ori folding structure. In this section, the force analysis of a rigid deployment is first conducted to estimate the magnitude of actuator force required to fully deploy the structure. Then, the flexible structure is deployed under a slope force F 1 . A quasi-static analysis is adopted to track the configuration at different load steps. At last, the deployment force is planned based on the quasi-static result to achieve a deployment with a constant speed in the increase of deployment ratio and the corresponding dynamic analysis is presented.

Force analysis of rigid deployment
Besides the driving mode, the maximum magnitude of actuator force required to achieve a full deployment is another key issue in the design of the deployable structure. The force analysis of rigid deployment can be a good vehicle to estimate the magnitude of actuator force for the membrane deployment. In the rigid deployment, given a configuration, the actuator force can be obtained via the principle of virtual work. Denote the anti-deployment moment associated with crease i as M i , i = 1, 2, 3, 4. Under a specific configuration, the virtual displacements along directions of four dragging force at corner nodes are δu i , i = 1, 2, 3, 4. Exploiting the symmetry, δu 1 = δu 2 and δu 3 = δu 4 . The principle of virtual work can be written as Substituting F 4 = F 1 tan(γ 0 /2) into Eq. (30), one has Denoting the unit direction vector associated with F i asF i , i = 1, 2, 3, 4, the corresponding virtual displacements δu i , i = 1, 2, 3, 4 are where the position vector of node E, F and the unit direction vectorsF i , i = 1, 4 are Substituting Eq. (33) to Eq. (32), the virtual displacements can be written as Recalling Eq. (21), the relationship between δγ 1 and δγ 2 is Substituting Eq. (35) to Eq. (34), the virtual displacements δu 1 and δu 4 can be written as δu 1 = L (tanγ 2 cosγ 1 − sinγ 2 )sin γ 0 2 + cosγ 2 cos γ 0 2 δγ 2 δu 4 = L (tanγ 2 cosγ 1 + sinγ 2 )cos γ 0 2 + cosγ 2 sin γ 0 2 δγ 2 Recalling Eq. (20), the relationship between δγ 2 and δθ 1 is Substituting Eq. (28), Eqs. (37), and (36) into Eq. (31), the actuator force can be obtained as (38) where the anti-deployment moment M i (θ i ), i = 1, 2 are obtained from the experiment curve in [28]. Considering that the structure is difficult to be deployed to the ideal deployment plane, a close configuration to the ideal plane is adopted to estimate the magnitude of F 1 required for the flexible deployment. Based on the force analysis of rigid deployment, it is found that when γ 1 = 4.14 • , θ 1 = 177.9 • , and θ 2 = 171.5 • , the deployment ratio reaches ρ = 99.7%. Substituting the angels and corresponding moment M i , i = 1, 2 into Eq. (38), the actuator force is obtained as F 1 = 12.73 N, which is close to the value 12.69 N adopted in [28]. For sake of comparison, the magnitude adopted in [28] is used in the following deployment simulation of flexible model.

Quasi-static analysis of flexible deployment
Quasi-static analysis is conducted in this section to track the deployment process under different load steps. The full load is F 1 = 12.69 N. It can be found that the structure deploys fast in the beginning stage with the increase in F 1 . Therefore, small load steps are used at the beginning. The deployment ratio is used to check the convergence of the ANCF model. Different meshes 4×4, 6×6, 8×8, and 10×10 are used for the discretization of each facet. The convergence study is shown in Fig. 12.   Fig. 12 Deployment ratio of flexible Miura-ori structure As shown in Fig. 12, the solution converges at the mesh of 8 × 8. Clearly, the Miura-ori structure under a slope force F 1 exhibits a highly non-uniform deployment. The structure deploys fast at the beginning and slows down at the following period. Specifically, the deployment ratio ρ jumps to 90% with only an approximate 10% increase in F 1 . After the load factor of F 1 increases to around 15%, ρ remains stable. Only a slight rise from 94% to 99.7% in ρ is observed when the load factor of F 1 changes from 0.15 to 1.
It is estimated that the non-uniform deployment will lead to an additional vibration in the dynamic deployment process. Therefore, the deployment force is planned based on the quasi-static analysis solution to achieve a stable deployment. A deployment ratio increase from the initial 0.03 to the final 0.997 with a constant speed is used as the planning goal. The mapping between ρ and F 1 based on the converged solution given by 8 × 8 model is adopted. The planned force is shown in Fig. 13.
As illustrated in Section 3.2, there is a pre-stress in the initial configuration e f 0 as shown in Fig. 14. The maximum von-Mises stress is around 2 MPa and appears in the region near the intersection point o.

Dynamic analysis of flexible deployment
The dynamic deployment analysis under a slope increase of the deployment force is first conducted. The density of the membrane is 1420 Kg/m 3 . The same final magnitude of F 1 = 12.69 N is used and the time period is 10 s. The configuration and the distribution of von-Mises  As shown in Fig. 15, the Miura-ori structure is deployed fast at the beginning stage, which is consistent with the observation of quasi-static analysis.
The stress concentration appears in a very small region near the crease intersection point as shown in Fig. 15a, b or the force dragging corner as shown in Fig. 15c, d. The stress in the stress concentration point reaches around 160 MPa, which is much bigger than that in the remained region. The stress concentration can be alleviated by a better force application means. To better show the stress distribution, the significant stress concentration in the corresponding point is excluded via setting the stress limit to 90 MPa in Fig. 15a, b, and d. A relatively big stress appears in the crease region as shown in Fig. 15. Excluding the stress concentration, the maximum stress is around 70∼80 MPa and appears at some small local region near the crease intersection point. Considering the yield stress of the membrane is 65 MPa [28], the maximum stress will cause plastic deformation in such region, which is consistent with the result in [28]. It should be noticed that the stress in the other most region including the stripes along the crease remain under the yield stress. The deployment ratio obtained by dynamic analysis is compared with that given by the quasi-static analysis in Fig. 16.
As shown in Fig. 16, there are significant vibrations in the dynamic deployment analysis compared with the quasi-static analysis. To alleviate the vibration, the planned force as shown in Fig. 13 is used for the dynamic simulation. The obtained deployment ratio is shown in Fig. 17.  Figure 17 shows that the planned force F 1 can deploy the structure with almost a constant speed of increase in the deployment ratio. Due to the small inertial of the membrane, vibration still occurs during the deployment process. However, the vibration magnitude is much smaller than that in the dynamic analysis under a slope F 1 as shown in Fig. 16.

Conclusion
This paper explores an approach to simulate the membrane deployment accounting for the nonlinear antideployment effect of creases. The origami membrane is analyzed as a flexible multibody system model using the absolute nodal coordinate formulation. The nonlinear behavior of creases are integrated into the multibody system via the principle of virtual work. The facets are meshed with the thin ANCF shell element. The anti-deployment property of a crease is parameterized based on the experimental data in [28]. Creases are explicitly modeled as nonlinear virtual torsional springs with specific constraints to connect the flexible facets. The deployment of a Z-folding membrane is analyzed and the approach is verified by the numerical/experimental reference in [28]. Then, the approach is used to simulate the deployment of a Miura-ori origami unit membrane. The geometric analysis of the corresponding rigid Miura-ori origami structure is conducted, based on which a form-finding method is proposed to find the initial configuration considering the crease effect. Meanwhile, an estimation on the magnitude of the driving force in the flexible deployment is provided. Setting the equilibrium configuration of rigid origami structure as the initial input, the counterpart of the flexible origami is obtained via the Newton-Raphson iterations. The deployment of a unit Miura-ori origami under a slope force is simulated. The quasistatic and dynamic solutions are compared and discussed. Based on the quasi-static result, the driving force is planned to achieve a more stable deployment and the corresponding dynamic analysis is presented. This work shows a possible way to improve the accuracy of membrane deployment simulation by exploiting the merit of ANCF method in nonlinear analysis and by considering the anti-deployment effect of the crease.