CFRP Origami Metamaterial with Tunable Buckling Loads: A Numerical Study

Origami has played an increasingly central role in designing a broad range of novel structures due to its simple concept and its lightweight and extraordinary mechanical properties. Nonetheless, most of the research focuses on mechanical responses by using homogeneous materials and limited studies involving buckling loads. In this study, we have designed a carbon fiber reinforced plastic (CFRP) origami metamaterial based on the classical Miura sheet and composite material. The finite element (FE) modelling process’s accuracy is first proved by utilizing a CFRP plate that has an analytical solution of the buckling load. Based on the validated FE modelling process, we then thoroughly study the buckling resistance ability of the proposed CFRP origami metamaterial numerically by varying the folding angle, layer order, and material properties, finding that the buckling loads can be tuned to as large as approximately 2.5 times for mode 5 by altering the folding angle from 10° to 130°. With the identical rate of increase, the shear modulus has a more significant influence on the buckling load than Young’s modulus. Outcomes reported reveal that tunable buckling loads can be achieved in two ways, i.e., origami technique and the CFRP material with fruitful design freedoms. This study provides an easy way of merely adjusting and controlling the buckling load of lightweight structures for practical engineering.


Introduction
Buckling is a failure mode that often occurs when the structure is under compressive stress. The typical characteristic of buckling is the sudden lateral displacement of the structural members, causing structural instability and leading to the critical importance of designing structures with buckling resistance ability, particularly for thin rods or thinwalled structures [1].
Generally, there are three main methods to compute the buckling loads or enhance one structure's buckling resistance ability, i.e., analytical methods, numerical methods, and test methods. The premise of utilizing analytical methods is that the structure's geometry is relatively simple, for example, rods, beams, and shells [2]. Test methods generally involve expensive costs, so relatively few experimental tests are conducted, which are often employed to validate analytical methods and numerical methods [3]. Various numerical methods have been proposed to deal with the buckling problems with complex geometry topology [4][5][6][7][8][9]. Although topology optimization can achieve novel designs that possess excellent buckling resistance ability, the buckling load is determined and hard to change once designed. However, a tunable mechanical property is in higher demand, which is generally realized by structural design with multiple materials [10,11] or the origami technique [12,13]. The latter is considered in this study.
Origami, originated from traditional culture, has become a powerful tool for designing novel engineering structures with various intriguing mechanical properties. An origami structure can be divided into two categories, according to whether it has an open or arranging the unit cells in the x and y directions, the Miura sheet's geometric topology is further determined by the folding angle, θ. These four parameters can be seen in Figure 1B, and have the following rigorous mathematical relations [29]: cos γ = sin 2 β cos θ + cos 2 β (1A) cos δ = sin 2 β cos 2 (θ/2) − cos 2 β sin 2 β cos 2 (θ/2) + cos 2 β (1B) w = 2b sin(γ/2) (1C) h = a cos(δ/2) (1D) It should be noted that, in addition to these four parameters, other parameters in Equations (1A-F) can refer to reference [27]. To be brief without losing generality, we choose the Miura sheets with one unit cell in the y direction and n unit cells in the x direction to construct the origami metamaterial. Therefore, the designed origami tube was finally identified by five independent parameters, i.e., the length a, width b, acute angle β, the folding angle θ, and the number of the unit cells in the x direction, n. For example, Figure 1B shows a typical case when n is equal to 4. In this study, we employed carbon fiber reinforced plastic (CFRP) rather than homogeneous material, e.g., metal, to constitute the origami metamaterial to further improve the design freedom, thereby emerging the CFRP origami metamaterial as depicted in the lower part of Figure 1A. It should be noted that, in addition to these four parameters, other parameters in Equations (1A-F) can refer to reference [27]. To be brief without losing generality, we choose the Miura sheets with one unit cell in the y direction and n unit cells in the x direction to construct the origami metamaterial. Therefore, the designed origami tube was finally identified by five independent parameters, i.e., the length a, width b, acute angle β, the folding angle θ, and the number of the unit cells in the x direction, n. For example, Figure 1B shows a typical case when n is equal to 4.
In this study, we employed carbon fiber reinforced plastic (CFRP) rather than homogeneous material, e.g., metal, to constitute the origami metamaterial to further improve the design freedom, thereby emerging the CFRP origami metamaterial as depicted in the lower part of Figure 1A.

Numerical Simulations
We conducted numerical simulations to identify the buckling loads of the CFRP origami metamaterial by the commercial finite element code Abaqus/Linear perturbation [38]. Two identical Miura sheets were rigidly connected at the junction surface, as described in Section 2.1. Regarding the CFRP origami metamaterial boundary conditions, all degrees of freedom on the left side were constrained. Four axial compression forces of 100 N were applied on the right end nodes; the loading was directly applied in Abaqus in the buckle step, and the type of loading was concentrated force. The whole CFRP origami metamaterial was meshed by 12,800 four-node shell elements with reduced integration (S4R) and an approximate size of 0.5 mm, as shown in Figure 2. The shell element was a four node thick shell, with reduced integration, hourglass control, and finite membrane strains. There were a total of eight layers of CFRP, symmetrically laid; the sequence of composite lamination was [45/−45/90/0], with a thickness of each layer of 0.125 mm. The material properties are listed in Table 1. E 1 , E 2 , and E 3 represent the material's elastic modulus in the main direction of elasticity 1, 2, and 3, respectively; υ 12 , υ 13 , and υ 23 are the negative values of the ratio of strain in the 1 direction to strain in the 2 direction, strain in the 1 direction to strain in the 3 direction, and strain in the 2 direction to strain in the 3 direction, respectively. G 12 , G 13 , and G 23 are the shear modulus in the 1-2 plane, 1-3 plane, and 2-3 planes, respectively. Material properties were set in the material editing interface, i.e., [mechanical]/[Elasticity]/[Elastic]/[Lamina]. The geometric parameters to determine the topology of the CFRP origami metamaterial are: a = 10 mm, b = 10 mm, β = 55 • , θ = 130 • , and n = 4. It should be underlined that the buckling load was obtained by multiplying the axially applied load by the Eigenvalue, and the Eigenvalues can be extracted from the numerical simulations, e.g., the buckling load of mode 1 can be calculated by multiplying the axially applied load by the Eigenvalue corresponding to mode 1; other situations are similar.
(S4R) and an approximate size of 0.5 mm, as shown in Figure 2. The shell element was a four node thick shell, with reduced integration, hourglass control, and finite membrane strains. There were a total of eight layers of CFRP, symmetrically laid; the sequence of composite lamination was [45/-45/90/0], with a thickness of each layer of 0.125 mm. The material properties are listed in Table 1. E1, E2, and E3 represent the material's elastic modulus in the main direction of elasticity 1, 2, and 3, respectively; υ12, υ13, and υ23 are the negative values of the ratio of strain in the 1 direction to strain in the 2 direction, strain in the 1 direction to strain in the 3 direction, and strain in the 2 direction to strain in the 3 direction, respectively. G12, G13, and G23 are the shear modulus in the 1-2 plane, 1-3 plane, and 2-3 planes, respectively. Material properties were set in the material editing interface, i.e., [mechanical]/[Elasticity]/[Elastic]/[Lamina]. The geometric parameters to determine the topology of the CFRP origami metamaterial are: a = 10 mm, b = 10 mm, β = 55°, θ = 130°, and n = 4. It should be underlined that the buckling load was obtained by multiplying the axially applied load by the Eigenvalue, and the Eigenvalues can be extracted from the numerical simulations, e.g., the buckling load of mode 1 can be calculated by multiplying the axially applied load by the Eigenvalue corresponding to mode 1; other situations are similar.

Theoretical Prediction of the Axial Buckling Loads for a Classical CFRP Plate
Considering the complex geometry and the complicated constitutive relationship of the proposed CFRP origami metamaterial, we employed a classical CFRP plate, which has a theoretical solution, to verify the FE modelling process in terms of the buckling load.
When the unidirectional composite is in the principal axis direction, the elastic constant is calculated as follows. The relationship between the stress and the strain is given as [39]:

Theoretical Prediction of the Axial Buckling Loads for a Classical CFRP Plate
Considering the complex geometry and the complicated constitutive relationship of the proposed CFRP origami metamaterial, we employed a classical CFRP plate, which has a theoretical solution, to verify the FE modelling process in terms of the buckling load.
When the unidirectional composite is in the principal axis direction, the elastic constant is calculated as follows. The relationship between the stress and the strain is given as [39]: When the unidirectional composite is in the principal off-axis direction, the elastic constant is calculated as follows. The stress and strain have the following relationship [39]: Defining the angle between the coordinate axis and the fiber direction as θ, one gets: Elastic properties of laminated plates can be calculated as: in-plane strain. K represents the change in curvature and torsion of the middle plane before and after the deformation. A is the tensile stiffness matrix, B is the tensile and bending coupling stiffness matrix, Axial buckling load can be finally computed as: where N x is the axial buckling load per unit length, γ is the half wave number of buckling in the direction of the plate, and l and w are the length and width of the plate, respectively.

Validation of the FE Modelling Process for Calculating the Buckling Load
The CFRP plate was 300 mm in length and 200 mm in width, subject to a compressive load of 100 N/mm at the right side. The plate structure's left end was fixed, and the right end was fixed with other degrees of freedom, except for the axial direction. Based on the theory presented in Section 2.3, one can get: Note that the units for the matrix A, B, and D are MPa·mm, MPa·mm 2 , and MPa·mm 3 , respectively. The number of half-waves can be determined from the buckling mode, as shown in Figure 3. The first five buckling modes are depicted in Figure 4, with the corresponding Eigenvalues as 5.77633 × 10 −2 , 6.57919 × 10 −2 , 7.32183 × 10 −2 , 9.82133 × 10 −2 , and 0.13141, respectively. The first five buckling loads extracted from theoretical solutions and numerical simulations are compared in Figure 4. It can be seen that the axial buckling loads obtained from the theoretical computation were 1248.4 N, 1399.6 N, 1570.8 N,  2097.4 N, and 2795.2 N, respectively, while the axial buckling loads acquired from the numerical analysis were 1155.3 N, 1315.8 N, 1464.4 N, 1964.3 N, and 2628.2 N, respectively. The numerical simulation's relative errors compared with the theoretical computation were then calculated as −7.46%, −5.99%, −6.77%, −6.35%, and −5.97%, respectively, indicating that the FE modelling has relatively high accuracy. Therefore, we will employ the FE modelling process to predict the axial buckling loads of the proposed CFRP origami metamaterial in the succeeding paragraphs.     Figure 5 shows the influences of the folding angle on buckling loads of the CFRP origami metamaterial by varying the folding angle θ from 10 • to 130 • at an interval of 40 • and with other parameters fixed. The geometric topologies of these CFRP origami metamaterials are also depicted on the right side of Figure 5. Figure 5 26 kN), respectively, as θ enlarges from 10 • to 130 • . It can be seen that the buckling loads can be widely tuned by merely varying the folding angle, especially for high-order modes. Specifically, the larger the folding angle, the larger the buckling load for modes 3-6, while an unusual law is discovered for the first two modes, i.e., the buckling load decreased first (θ from 10 • to 50 • ) and then increased (θ from 50 • to 130 • ), and the buckling load changed more drastically for mode 2 than for mode 1 when θ exceeded 90 • . It is interesting to find that, for mode 1 and mode 3, the buckling loads were approximately equal when θ = 10 • and θ = 50 • , respectively, enriching the design freedom of the CFRP origami metamaterial in terms of the ability to resist buckling. Moreover, it can be surprisingly found that the buckling loads can be tuned by as large as approximately 2.5 times for mode 5. Thus, the widely tuned buckling loads can be realized for the proposed origami metamaterial by altering the folding angle with the base material properties unchanged.   Figure 5 shows the influences of the folding angle on buckling loads of the CFRP origami metamaterial by varying the folding angle θ from 10° to 130° at an interval of 40° and with other parameters fixed. The geometric topologies of these CFRP origami metamaterials are also depicted on the right side of Figure 5. Figure 5 26 kN), respectively, as θ enlarges from 10° to 130°. It can be seen that the buckling loads can be widely tuned by merely varying the folding angle, especially for high-order modes. Specifically, the larger the folding angle, the larger the buckling load for modes 3- 6, while an unusual law is discovered for the first two modes, i.e., the buckling load decreased first (θ from 10° to 50°) and then increased (θ from 50° to 130°), and the buckling load changed more drastically for mode 2 than for mode 1 when θ exceeded 90°. It is interesting to find that, for mode 1 and mode 3, the buckling loads were approximately equal when θ = 10° and θ = 50°, respectively, enriching the design freedom of the CFRP origami metamaterial in terms of the ability to resist buckling. Moreover, it can be surprisingly found that the buckling loads can be tuned by as large as approximately 2.5 times for mode 5. Thus, the widely tuned buckling loads can be realized for the proposed origami metamaterial by altering the folding angle with the base material properties unchanged.  Figure 6. The layers were symmetrically laid, with a total of eight layers of CFRP. Each layer had a uniform thickness of 0.125 mm. Figure 7A-D shows the influences of the layer order on the buckling load when the folding angles are equal to 10°, 50°, 90°, and 130°, respectively.  When θ = 10 • , obvious influences could be found, particularly for high-order modes. The first-order buckling loads were 2.45 kN, 2.16 kN, 2.12 kN, and 2.07 kN, respectively, for these four cases, identifying that relatively slight changes were found for the last three cases. The buckling load corresponding to modes 2-6 had an identical change in law, i.e., it became bigger first, then decreased, and then became larger; CFRP origami metamaterial with the layer order [0/90/−45/45] and layer order [45/−45/0/90] had the largest and the least buckling load, respectively. It can also be seen that, by altering the layer order, the buckling load can be tuned by as large as roughly 50% (corresponds to mode 5).

Influences of the Folding Angle
Compared with the CFRP origami metamaterial with θ = 10 • , the buckling loads of the other three cases when the folding angles equaled 50 • , 90 • , and 130 • were found to be relatively less sensitive to the alteration of the layer order. For example, when θ = 50 • , all of the first-order and second-order buckling loads were about 1.1 kN and 2.2 kN, respectively. The differences in the buckling loads were almost within 5% for modes 3-6 as one turned the layer order. However, the situation changed for CFRP origami metamaterials with θ = 90 • and 130 • . For instance, when θ = 90 • , the least and largest buckling loads were (5.4 kN, 6.3 kN), (5.5 kN, 6.5 kN), (5.7 kN, 6.6 kN), and (5.9 kN, 7.0 kN) for modes 3-6, realizing 17%, 18%, 16%, and 19% tunable properties, respectively. CFRP origami metamaterial with the folding angle of 130 • had an apparent tunable buckling load performance only for high-order modes, i.e., mode 5 and mode 6. Specifically, the buckling load can be tuned from approximately 7.7 kN to 9.2 kN and from 7.9 kN to 9.9 kN, respectively, for mode 5 and mode 6 by adjusting the layer order. It was also interesting to find a significant jump in the buckling load from mode 2 to mode 3 in Figure 7C,D. This phenomenon may be caused by a large change in the buckling half-wave number in this adjacent mode. However, considering the complex geometry of the CERP origami metamaterial, this reason needs further exploration. Overall, the buckling load of the proposed CFRP origami metamaterial can be largely tuned via simply modifying the CFRP laminate layer order. When θ = 10°, obvious influences could be found, particularly for high-order modes. The first-order buckling loads were 2.45 kN, 2.16 kN, 2.12 kN, and 2.07 kN, respectively, for these four cases, identifying that relatively slight changes were found for the last three cases. The buckling load corresponding to modes 2-6 had an identical change in law, i.e., it became bigger first, then decreased, and then became larger; CFRP origami metamaterial with the layer order [0/90/-45/45] and layer order [45/-45/0/90] had the largest and the least buckling load, respectively. It can also be seen that, by altering the layer order, the buckling load can be tuned by as large as roughly 50% (corresponds to mode 5).
Compared with the CFRP origami metamaterial with θ = 10°, the buckling loads of the other three cases when the folding angles equaled 50°, 90°, and 130° were found to be relatively less sensitive to the alteration of the layer order. For example, when θ = 50°, all of the first-order and second-order buckling loads were about 1.1 kN and 2.2 kN, respectively. The differences in the buckling loads were almost within 5% for modes 3-6 as one turned the layer order. However, the situation changed for CFRP origami metamaterials with θ = 90° and 130°. For instance, when θ = 90°, the least and largest buckling loads were

The Material Properties
We finally explored how the buckling load of the proposed CFRP origami metamaterial changed if the base material parameters, i.e., Young's modulus and shear modulus, were altered. It should be underlined that the material properties (Young's modulus and shear modulus) selection followed the following guidelines: the second group was the reference group (parameters used in the previous section); the values of the first group were half of the reference group; and the values of the third group were the first group plus the second group. Figure 8A shows the influences of Young's modulus on the buckling loads with other material properties unchanged. Three representative cases were considered, i.e., CFRP origami metamaterial with Young's modulus of (E 1 = 72. 35 GPa, E 2 = 4.825 GPa), (E 1 = 144.7 GPa, E 2 = 9.65 GPa), and (E 1 = 217.05 GPa, E 2 = 14.475 GPa), respectively. It can be found that increasing the Young's modulus values leads to the increase of the buckling load for all of the modes, revealing that Young's modulus has a positive influence on the ability to resist buckling instability. For example, the first-order buckling load can be enlarged by roughly 112% (from 1.7 kN to 3.6 kN) as Young's modulus tends to be larger. It is also interesting to find that there are two special cases in which the buckling loads were almost kept constant, although Young's modulus was enlarged by 1.5 times (increase in Young's modulus from E 1 = 144.7 GPa and E 2 = 9.65 GPa to E 1 = 217.05 GPa and E 2 = 14.475 GPa for mode 2 and mode 4, respectively). This unusual phenomenon can allow for the designing of the CFRP origami metamaterial with specific buckling loads. Figure 8B depicts the influences of the shear modulus on the buckling load, in which a similar change in the law can be found. The buckling load became larger with the shear modulus increase, except for mode 1, which showed a relatively insignificant change. Unlike Young's modulus, no approximately identical buckling load could be discovered for the CFRP origami metamaterial when the shear modulus varied. Moreover, it can be found that the six-order buckling load can be improved by almost 174% when the shear modulus is increased from (G 12 = 2.6 GPa, G 13 = 2.6 GPa, G 23 = 1.7 GPa) to (G 12 = 7.8 GPa, G 13 = 7.8 GPa, G 23 = 5.1 GPa). By comparing Figure 8A,B, it can be found that, with the same increase rate, the shear modulus has a more significant impact on the buckling load than Young's modulus. In summary, the buckling load of the proposed CFRP origami metamaterial can be easily tuned by simply changing Young's modulus and the shear modulus of the base material.

Conclusions
This study has proposed a novel origami metamaterial with a wide-range tunable buckling load based on the Miura sheet and carbon fiber reinforced plastic (CFRP). The tunable buckling load property of the proposed CFRP origami metamaterial is thoroughly investigated using numerical simulations, whose finite element modelling process is verified by employing a CFRP plate with theoretical solutions. Results reveal that the ability to resist buckling instability for the CFRP origami metamaterial can be widely tuned through two ways, i.e., altering the folding angle by utilizing the origami merit and changing the layer order and base material properties by the fruitful design freedoms of the CFRP. Moreover, it seems that, with the same rate of increase, the shear modulus has a more significant influence on the buckling load than Young's modulus. This research provides a solution for the design of a lightweight structure with the demand for tunable buckling resistance in actual engineering.

Conclusions
This study has proposed a novel origami metamaterial with a wide-range tunable buckling load based on the Miura sheet and carbon fiber reinforced plastic (CFRP). The tunable buckling load property of the proposed CFRP origami metamaterial is thoroughly investigated using numerical simulations, whose finite element modelling process is verified by employing a CFRP plate with theoretical solutions. Results reveal that the ability to resist buckling instability for the CFRP origami metamaterial can be widely tuned through two ways, i.e., altering the folding angle by utilizing the origami merit and changing the layer order and base material properties by the fruitful design freedoms of the CFRP. Moreover, it seems that, with the same rate of increase, the shear modulus has a more significant influence on the buckling load than Young's modulus. This research provides a solution for the design of a lightweight structure with the demand for tunable buckling resistance in actual engineering.