Three-Dimensional-Printed Multistable Mechanical Metamaterials With a Deterministic Deformation Sequence

Mechanical metamaterials are materials whose effective properties arise from the underlying architecture, rather than from the bulk behavior of their constituents [1]. For example, acoustic metamaterials can exhibit unusual acoustic behavior, such as the ability to guide or stop elastic wave propagation along a desired path [2–4]. Other interesting and unconventional properties of metamaterials include negative bulk modulus [5–8], negative Poisson's ratio [9,10], high specific energy absorption [11,12], or negative (dynamic) mass [5,8,13,14].

However, the functionality of a mechanical metamaterial is limited if its properties cannot be tuned after fabrication of the metamaterial. Recent studies have shown great interest in developing tunable metamaterials by exploiting elastic instabilities that can cause configurational changes in the microarchitecture of these metamaterials [15]. For example, reversible changes in the microarchitecture geometry due to elastic instabilities have been used to tune the bandgaps in acoustic [16–20] and photonic [21,22] metamaterials.

A snap-through instability is a kind of elastic instability in which a structure instantaneously jumps from one configuration to another configuration when an applied stimulus reaches a critical level [23]. The design of mechanical metamaterials with snap-through instabilities has been the focus of active research in recent years [24–28]. For example, Correa et al. [25] and Restrepo et al. [26] developed metamaterials that are made up of multiple unit cells that each includes an initially curved beam. Through the elastic snap-through of high number of unit cells arranged in series, this kind of metamaterial can exhibit large hysteresis loops in response to cyclic loads. Because of these hysteresis loops, these metamaterials dissipate a large amount of mechanical energy. In contrast to traditional honeycombs that rely on plastic energy dissipation, the deformation can be fully recovered in these metamaterials. While some of the proposed metamaterials with snap-through instabilities are monostable, i.e., they return to their undeformed configuration after removal of the load (for example, the negative stiffness honeycombs designed by Correa et al. [25]), other designs are multistable, i.e., they have multiple stable configurations. For example, Shan et al. [24] developed multistable metamaterials using bistable tilted beams. These metamaterials can achieve large deformations and offer significant energy absorption capacity due to the ability of the material to trap elastic energy in a stable deformed configuration in response to an impact. Haghpanah and Salari-Sharif [28] investigated another design for multistable metamaterials, based on bistable structural elements designed to be considerably stronger than the previous designs. While these previous studies focus on multistable metamaterials in response to compressive loads, Rafsanjani et al. [27] showed that the design of architectures that exhibit multistability in response to tensile loads is possible.

In multistable metamaterials, a sequence of stable configurations is obtained as the stimulus amplitude is increased. However, this sequence (which we will call the deformation sequence) is unpredictable if, as in most previous works, the metamaterials consist of identical unit cells [24,26]. Imperfections in the geometry, material properties [29–31] or in the boundary conditions determine the deformation sequence in the case of a multistable metamaterial with identical unit cells. Although previous papers about multistable metamaterials give important insight into the design and mechanics of these metamaterials, these previous works have not investigated strategies that would make it possible to control their deformation sequence.

The objective of this paper is to develop and analyze approaches to tune the deformation sequence of one kind of multistable metamaterial. To achieve this objective, a rigorous theoretical model, based on an analytical model for a single unit cell and on the minimization of the total strain energy, is developed to analyze the nonlinear mechanics and deformation sequence of multistable metamaterials. Experiments with 3D-printed samples are used to test the theory. One possible application of metamaterials with a deterministic deformation sequence could be the control of elastic waves. In a multistable metamaterial with a deterministic deformation sequence, a preload could be applied to cause a predictable switch to another stable configuration; this switch in the geometry of the microarchitecture would affect wave propagation due to a predictable change in the frequency bandgaps of the metamaterials. In contrast to previous tunable metamaterials which requires to maintain the stimulus [16,19], these metamaterials would remain in a deformed configuration after removal of the preload due to multistability. Thus, although tuning the deformation sequence only requires snap-through buckling, multistable metamaterials are the focus of this paper because multistability could have advantages in applications such as the development of tunable metamaterials with switchable properties.

The organization of this article is as follows. First, the theoretical model and the experimental methods are presented. Then two methods, both based on small variations in the unit cell geometry, are proposed to obtain a deterministic deformation sequence. These methods are validated using theoretical analysis, finite-element simulations, and experiments performed on 3D-printed samples.

The mechanics of a specific kind of multistable metamaterial, shown in Fig. 1(a), is investigated using a theoretical model. The metamaterial is made up of multiple bistable unit cells (that are similar to [26]), whose geometrical parameters are shown in Fig. 1(b). The unit cell consists of thick horizontal and vertical elements and a thin curved part. As in the previous studies [24,26], because the loading is intended to be only in the y-direction, these metamaterials are called one-dimensional (1D) metamaterials.

where U is the strain energy, and an overbar denotes a normalized quantity. The procedure and equations needed to obtain the relationship between the normalized force, normalized strain energy, and normalized deformation are given in Appendix  A.

The theoretical model was validated using nonlinear finite-element analysis (FEA) in abaqus/standard using the static/general nonlinear procedure. The FEA model is built with four-node bilinear quadrilateral plane-stress element with reduced integration (CPS4R in abaqus) and an elastic, isotropic constitutive model. The model is stimulated using displacement control, such that the vertical (y) displacement on the top edge of the unit cell is prescribed. The other boundary conditions are the following: the bottom edge of the unit cell is fixed, and the top edge of the unit cell has a fixed x-direction displacement.

Results of Figs. 1(d)–1(g) show that there is excellent agreement between the theoretical model and the FEA simulation, both in the case a₃ = 0 which corresponds to the model developed by Qiu et al. [32] (Figs. 1(d) and 1(e)), and in the case a₃ ≠ 0 which corresponds to the model extension of this paper (Figs. 1(f) and 1(g)). It can be seen in the normalized force versus normalized deformation curve that there is one critical point whose normalized force value is $f¯cr$⁠, from which point the beam begins to exhibit a negative stiffness. In a force control simulation, this point would correspond to the point at which the unit cell would snap through. In the normalized strain energy versus normalized deformation curve, there is a global minimum (point I in Figs. 1(e) and 1(g)) when $d¯=0$⁠, and a local minimum (point III in Figs. 1(e) and 1(g)) when $d¯≈2$⁠. These are the two stable equilibrium positions that correspond to points I and III in Figs. 1(d) and 1(f). Moreover, there is one local maximum (point II in Figs. 1(e) and 1(g)) in the strain energy versus deformation curve, which is an unstable equilibrium position corresponding to point II in Figs. 1(d) and 1(f).

The mechanics of a multistable metamaterial that is comprised of a series of bistable unit cells can be derived from the model for a single unit cell. Consider the metamaterial with n unit cells in series (Fig. 2(a)) and the corresponding theoretical model shown in Fig. 2(b). In the theoretical model, the beams are connected by a rigid link. The bottom beam (beam 1) has clamped–clamped boundary conditions; the left and right sections of the other beams (beams 2, 3,…, n) are only allowed vertical (y-direction) displacement. The thicknesses and mode imperfection sizes of the curved part of each unit cell are t_i and a_3,i, respectively. For all unit cells, the initial height of the curved parts is h. The displacements of the midpoint of the curved part of unit cell i is u_i (see Fig. 2(b)); the deformations of the unit cells can be expressed as d₁ = u₁, $di=ui−ui−1$ for i = 2,…, n, and the total deformation is $dtot=∑i=1ndi$⁠. The total deformation, d_tot, the total strain energy, U_tot, and the loading force, f, are normalized with respect to the parameters of unit cell 1, i.e., the following definitions are used for the normalized force, $f¯$⁠, the normalized deformation, $d¯tot$⁠, and the normalized strain energy, $U¯tot$⁠:

where $f¯i, d¯i$⁠, and $U¯i$ are the normalized force, deformation, and strain energy for unit cell i, respectively.

This optimization problem is solved using a local optimization algorithm (the interior-point algorithm using the fmincon function in matlab). Using a local minimization algorithm is critical; during loading the metamaterial does not switch to a state that corresponds to the global minimum but remains in the neighborhood of the current configuration.

Once the values of $d¯i$ are obtained at all steps for the loading (or unloading), the values of $d¯i$ can be plotted as a function of $d¯tot$ to analyze the sequence of deformed configurations the metamaterial goes through when the loading parameter is increased. In the case of a metamaterial with 2 unit cells in series, the curve of $d¯1$ as a function of $d¯tot$ completely characterizes how the material deforms during loading. This curve will be referred to as the deformation path of the metamaterial.

In addition to the theoretical model, a corresponding FEA model with elastic material is built using the boundary conditions shown in Fig. 2(a), i.e., the bottom edge of unit cell 1 is fixed, the left and right sides of these unit cells have a fixed x-direction displacement, and the top edge of unit cell n also has a fixed x-direction displacement. The deformation path and the force versus deformation curves predicted by a static nonlinear FEA simulation with displacement control are used to validate the theoretical model.

Samples were fabricated by a multimaterial 3D printer (Objet Connex 260, Stratasys, Edina, MN). The printing material is DM9895, which is a digital material derived by mixing two base materials. One of the two base materials is TangoblackPlus, which is a rubbery material at room temperature; the other one is Verowhite, which is a rigid plastic at room temperature. The in-plane dimensions of the printed samples are about 10 cm × 10 cm; the out-of-plane thickness is 1 cm.

In order to determine the deformation sequences of printed samples, compression tests were conducted on a universal Material Testing System (MTS, Model Insight 10, Eden Prairie, MN) in a displacement control manner with a 10 kN load cell. The samples were compressed (at room temperature) using a customized compression fixture (Fig. 6(b)) at a testing velocity of 10 mm/min until all unit cells collapse. A digital camera was used to record the whole testing process.

In this section, the theoretical model for elastic multistable metamaterials is first applied to a metamaterial that consists of two unit cells in series to demonstrate the validity of the model. The same algorithm is then applied to design metamaterials with a larger number of unit cells.

Using the theoretical method from Sec. 2.1, the total strain energy contours of the metamaterial with two unit cells in series can be obtained as a function of $d¯tot$ and $d¯1$⁠, as shown in Fig. 3. The strain energy contours show that there are four local minima, A, B, C, and D (that correspond to stable equilibria) and one local maximum of the total strain energy (that corresponds to an unstable equilibrium). The presence of multiple minima indicates the multistability of this metamaterial. The global minimum A is the initial undeformed configuration. The local minima B and C correspond to cases when only unit cell 2 has collapsed and only unit cell 1 has collapsed, respectively. The local minimum D is when both unit cell 1 and unit cell 2 have collapsed. Because the two unit cells are identical, the local minima B and C have the same strain energy value, which means the deformation path (the curve $d¯1(d¯tot)$⁠) during loading can follow the sequence A–B–D or A–C–D. In an experiment, the presence of imperfections, for example, caused by the manufacturing process or in the experimental setup (for instance in the application of the boundary conditions), would determine which of the deformation paths is preferred.

Based on the analysis of the metamaterial with identical unit cells, it can be seen that the sequence of stable configurations the metamaterial switches in response to a given stimulus is unpredictable if it consists of identical unit cells. In order to get a deterministic deformation sequence, small variations in the unit cell geometry are utilized, specifically varying the thicknesses or initial shapes of the curved parts of the unit cells. These variations will give different values of the critical forces (as defined in Sec. 2.1.1) for each unit cell. When the unit cells are assembled together, the unit cell with smallest critical force will deform first, followed by the unit cell with next smaller critical force. These two methods are analyzed next.

One method to obtain a deterministic deformation sequence is to use thickness variation, i.e., varying the thickness value of the curved part of the unit cell t_i from row to row; for example, t₂ = 1.5t₁ for a metamaterial architecture with two unit cells (Fig. 2(a)). The results for single unit cells show that the critical force value of unit cell 2, $f¯cr,2$⁠, is higher than the critical force value of unit cell 1, $f¯cr,1$⁠, which means unit cell 2 requires a larger loading force to collapse. Moreover, the strain energy of unit cell 2 is also greater than that of unit cell 1.

If the unit cells are combined in series, the theoretical results of the total strain energy contour, deformation path and normalized force versus deformation curves can be determined using the algorithm of Sec. 2.1.2 as the solid lines in Figs. 4(c)–4(e), respectively. As shown in the contour plot of this metamaterial's total strain energy (Fig. 4(c)), the four local minima of the strain energy are different, such that a unique deformation path is obtained. In order to validate the theoretical results, the simulation results for the corresponding FEA model are given in Figs. 4(d)–4(e), as the blue dashed lines. Note the excellent agreement for the deformation path between the theoretical model and the FEA simulation (Fig. 4(d)). Both the theoretical model and the FEA simulation show that the deformation sequence is A–C–D for loading and D–B–A for unloading. As shown in Figs. 4(c) and 4(d), when the normalized total deformation $d¯tot$ is about 3.4, there will be one snap-back of unit cell 1 (i.e., the value of $d¯1$ drops instantaneously) during loading when unit cell 2 snaps through. In the normalized force versus normalized deformation curve 4e, the first peak has a smaller amplitude than the second peak due to the lower stiffness of unit cell 1; futhermore, the force goes below the horizontal axis (which means a tensile force is observed) twice during loading.

Besides the method of varying the thickness, the deformation sequence can also be tuned by varying the mode shape imperfection size a₃ (where the imperfection in the initial shape of the beam, δW(x), is given by Eq. (4)). In order to prove that the deformation sequence can be tuned by mode shape imperfection, the force-deformation (Fig. 5(b)) and strain energy-deformation (Fig. 5(c)) curves are first compared for a unit cell without a mode shape imperfection (i.e., a₃ = 0) and a unit cell with a 10% third-mode imperfection (a₃ = 0.1). The results of both the theoretical model and FEA simulations show that the critical loading force and the strain energy of a single unit cell decrease when the model has a third-mode imperfection.

After obtaining the results for a single unit cell, the results for the metamaterial with two cells with different mode imperfection sizes can be determined using the same procedure as was used in Sec. 3.2. Consider the metamaterial with two unit cells in series, where unit cell 2 is perfect (a_3,2 = 0) and unit cell 1 has a 10% third-mode imperfection (a_3,1 = 0.1). In contrast to the case of identical unit cells, the local minima of the strain energy are different, which makes the deformation path unique (Fig. 5(d)). The deformation sequence is A–C–D for loading, and D–B–A for unloading. As in the case of the metamaterial with two cells using thickness variation, there will be one snap-back of the unit cell 1 during the initial snap-through of unit cell 2 (Fig. 5(d)), when the normalized total deformation $d¯tot$ is about 3.3. As in the case of thickness variation, the first peak in the normalized force versus normalized deformation curve has a smaller amplitude due to the smaller critical force in unit cell 1 (Fig. 5(f)); furthermore, a tensile force is observed twice during the loading. However, in contrast to the case of the thickness variation, the first peak has a smoother shape than the second peak because adding an imperfection significantly affects the shape of the normalized deflection versus normalized deformation curve for a single unit cell (see Fig. 5(a)).

Several multistable architectures were fabricated using a 3D printing technique. These samples were tested using a customized compression fixture (shown in Fig. 6(b)) to experimentally validate that the methods of thickness variation and mode shape imperfection can be used to tune the deformation sequence. Moreover, considering the viscoelasticity of the samples, FEA models with a viscoelastic constitutive model were built for comparisons with the experiments (see Appendix  B). In these models, the bottom edge is fixed and the top edge has a fixed x-direction displacement, while its y-displacement is prescribed to move at a constant velocity; the left and right edges are left free. The response to stimulation of the metamaterial at a constant strain rate was simulated using nonlinear quasi-static FEA simulations (VISCO procedure in abaqus) with displacement control.

One of the printed samples is the multistable metamaterial shown in Fig. 6(a) (metamaterial A), which consists of five rows, each with 5 unit cells (in the rest of the paper, we will call a metamaterial with N rows, each with M unit cells as an N × M metamaterial). For this 5 × 5 multistable metamaterial, the thickness is constant within each row, while the thickness varies from row to row (t₁ < t₅ < t₂ < t₃ < t₄, where the unit cells are ordered from bottom to top). Figures 6(b)–6(c) shows snapshots of the response of the multistable architecture at different effective strains ϵ for both experiment and finite-element simulation. These snapshots demonstrate that the deformation sequence matches what is expected from the previous analysis (i.e., row 1 collapses first, then row 5, row 2, row 3, and row 4).

Another printed sample is the 5 × 5 multistable metamaterial shown in Fig. 7(a) (metamaterial B). The unit cells of this 5 × 5 metamaterial have a uniform mode imperfection size within each row, but the imperfection size varies from row to row. The thickness of each unit cell's curved part is 1 mm. Snapshots of the response at different strains ϵ for both experiment and finite-element simulation (Figs. 7(b)–7(c)) demonstrate that the deformation sequence matches the expected sequence (i.e., row 1 collapses first, then row 5, row 2, row 3, and row 4).

The experiments were also analyzed more quantitatively by plotting the value of the normalized deformation of each unit cell, $d¯i$⁠, as a function of the total deformation, $d¯$⁠, in Figs. 8(a) and 8(b). The normalized deformation was computing using frames of the movies recorded using the digital camera during the compression test. For the experiments, FEA simulations (with the same boundary conditions as in Sec. 3.4.1) and theoretical model, the same clear deformation sequence can be seen by inspecting the results (i.e., row 1 collapses first, then row 5, row 2, row 3, and row 4 in both Figs. 8(a) and 8(b)). While both the FEA simulations and the theoretical model are able to qualitatively capture the deformation sequence, only the FEA simulations are able to quantitatively match the experiments. In particular, the unit cells in the theoretical model deform earlier than in the FEA models and experiments. Moreover, there is an instantaneous snap-back in theoretical model; however, a more gradual and limited snap-back is observed in the FEA model and in experiments (for example, see ϵ = 18% for $d¯1$ in Fig. 8(a) and ϵ = 15% for $d¯1$ in Fig. 8(b)). The inaccuracies of the theoretical model are due to two sources: (1) the viscoelasticity of the printed materials (see Appendix  B for the relaxation modulus of the material) and (2) the assumption regarding the kinematics in the theoretical model (a small rotation assumption is used and parts of the unit cell are assumed to be rigid). The smoothing of the curves seen in the experiments and FEA simulations (with a viscoelastic model) are primarily due to the viscoelasticity of the printed samples. Simulations with finite element with an elastic constitutive model (not shown in the manuscript) do show a discontinuous instability jump similar to what is seen in the theoretical model. In spite of the quantitative differences between the theory and experiments, the theoretical model does capture the important characteristics of the deformation paths of this multistable metamaterial and can be used to inform the design of metamaterials with a deterministic deformation sequence.

The effective stress versus strain curve was also obtained during the compression test and compared with the FEA model with viscoelastic material as shown in Figs. 9(a) and 9(c). The effective stress is calculated from the compression force divided by the bottom area of the model, and the effective strain is calculated from the total deformation divided by the initial height of the model. For both metamaterials, there is a good quantitative agreement between the experiments and FEA simulation results. This excellent quantitative agreement is interesting because of the difference in the boundary conditions applied in the FEA and in the experiments: the top edge of the sample was prescribed to move at a constant velocity in the FEA simulations; in the experiments, the top edge of the compression fixture was prescribed to move at a constant velocity and the top edge of the sample was not attached to the compression fixture. In the case of an elastic sample, a separation of the top edge of the sample from the compression fixture would have been expected when one the row of the samples snap-through. However, the presence of viscoelasticity significantly affects the kinetics of the snap-through process, such that the sample remains in contact with the compression fixture. The experiments and FEA simulations with a viscoelastic constitutive model do not show that the force crosses the horizontal axis, which could be wrongly interpreted as the absence of multistability. However, the multistability was experimentally verified by letting the metamaterial relax to any of the expected stable configurations. Furthermore, the first peak has a higher amplitude than the second peak, in contrast to what has been observed in Figs. 4 and 5.

Because of the viscoelasticity of the 3D-printed sample, the theoretical model cannot be directly compared to the experiments or the FEA simulations with a viscoelastic material; hence, the theoretical model was compared to FEA simulations with an elastic constitutive model (Figs. 9(b) and 9(d)). For both metamaterials, the theoretical result has excellent agreement with the FEA with elastic material and captures the major characteristics of this metamaterial, showing its multistability (i.e., crossings of the horizontal axis) and growing stress peaks. This implies that the absence of a tensile force and the fact that the second stress peak has a lower amplitude than the first peak are due to the viscoelasticity of the 3D-printed samples (note, however, that obtaining a tensile force would also require to use an adhesive between the compression fixture and the sample, as done, for example, in Ref. [24]). Because of the assumption of small rotation, some parts being rigid in the theoretical model, the values of the strain corresponding to maxima and minima are somewhat shifted in the theoretical model compared to FEA with an elastic constitutive model.

In this paper, we propose to use the methods of small geometric variations to tune the deformation sequence of one kind of metamaterial that consists of a quasi-periodic microarchitecture with multiple stable states. This work overcomes the limit shown in the previous studies that when the unit cells used in a multistable metamaterial are identical, the stable configuration that the metamaterial switches to is unpredictable due to the effects of imperfections in the manufacturing process and in the boundary conditions. Two different methods of thickness variation and mode shape imperfection are analyzed and used to obtain a deterministic deformation sequence. In order to obtain a deterministic deformation sequence, a rigorous theoretical model is developed to analyze the mechanics and help to design this metamaterial. The theoretical model is validated by comparison to finite-element simulations of a model that consists of multiple bistable unit cells in series. Both the loading and unloading response of the theoretical model are in excellent agreement with finite-element simulations.

Moreover, the results of finite-element simulations and experiments on 3D-printed samples demonstrate that the deformation sequence of multiple unit cells can be tuned using either varied thickness or higher mode shape imperfection of the unit cell curved parts. Despite viscoelasticity, which is not taken into account in the theoretical model, the theory does predict the experimentally observed deformation sequence. Excellent quantitative match between experiments and FEA simulations is observed when viscoelasticity is taken into account. Despite the presence of manufacturing imperfections, a deterministic deformation sequence is obtained for thickness variations or mode shape imperfections of moderate sizes. The validity of the proposed model was demonstrated for systems with a small number of unit cells and with simple variations in the model parameters. However, the proposed model is general and could be applied to more complex metamaterials with a large number of unit cells with spatial variations in multiple geometrical and material parameters and might be particularly useful for the analysis of multistable materials of larger sizes, which could be useful for the design multifunctional metamaterials with programmable and switchable properties.

This research was supported in parts by start-up funds from Georgia Tech, by NSF Award No. CMMI-1462894 to H. Jerry Qi and by the China Scholarship Council (Chao Yuan).

For the theoretical model of the initially curved beam, the second mode is constrained, i.e., only symmetric solutions were considered (because it is observed that the deformation mode of the unit cell is almost symmetric in FEA simulations), such that B₂ = 0.

In order to obtain the normalized force versus normalized deformation curve, a vector is first formed for the normalized compression force, $p¯$⁠. Equation (A10) is solved for each entry of this vector. If the discriminant of Eq. (A10), Δ, is negative, it means that the value of the normalized compression force is impossible. If Δ is positive, then two roots are obtained for the normalized force, $f¯$⁠. For each of these roots, the mode coefficients B_i can be obtained from Eq. (A9), and the normalized deformation can be derived from Eq. (A5). Thus, the value of normalized force can be plotted as a function of the normalized deformation for each root. Connecting the points obtained for each entry of the vector for $p¯$ that corresponds to Δ > 0, a curve $f¯$ versus $d¯$ is obtained for each of the two families of roots; these two curves converge to the point that corresponds to Δ = 0. A single normalized force versus normalized deformation curve is obtained by connecting the two curves.

If there is no third-mode imperfection (a₃ = 0), then, the solutions are comprised of two forms of solutions.

• Solution 1: if $p¯<N32$⁠, a solution of the form as in the case a₃ ≠ 0 is obtained. The equations obtained in the general case can be applied by setting a₃ = 0.

• Solution 2: If $p¯=N32$⁠, then another form of solution is obtained. In that case, the value of B₃ cannot be found directly using Eq. (A8), while the values of all other B_i are the same as in the case a₃ ≠ 0. Using Eq. (A5) and the equation $p¯=N32$⁠, a linear equation that relates the normalized force, $f¯$⁠, and the normalized deformation, $d¯$⁠, is obtained

  $f¯=−−N12N32−N12+d¯−1∑i=1,5,9,13…1Ni2(N32−Ni2)$

  (A12)

In the model without imperfections, the normalized force versus normalized deformation curve is obtained by connecting the two forms of solutions, as described in Ref. [32].

where G_∞ is the long-term shear modulus, G_i is the shear modulus of branch i, and τ_i is the relaxation time constant. These parameters were fit to relaxation data of one printed sample (3 mm × 1 mm × 10 mm) using the printed material (DM9895) at room temperature, as shown in Fig. 10. The Poisson's ratio was assumed to be equal to 0.495 and to be time-independent. The value of G_∞ is set to 0.0167 MPa, and the values of G_i and τ_i are shown in Table 1.