Notch Flexure as Kirigami Cut for Tunable Mechanical Stretchability towards Metamaterial Application

ABSTRACT

Kirigami is an art of paper cutting, which can be used in mechanical metamaterials, actuators, and energy absorption based on its deployable and load-deflection characteristics. Traditional cuts with zero width produce undesirable plastic deformation or even tear fracture due to stress concentration in stretching. This study proposes to enlarge the cut width into a notch flexure, which is applied to an orthogonality-cutted kirigami sheet, which buckles out of plane into a 3D configuration patterns under uniaxial tension. The use of compliant beam as the notch makes the stress distribution around the cuts more uniform in both elastic and elastoplastic regime. The experimental and numerical results show that by tuning the geometric parameters of cuts and material properties of the sheets, the trigger condition of 3D patterns can be adjusted. Potential capability of tunable phononic wave propagation in this kirigami-inspired metamaterial is demonstrated. This design methodology offers a theoretical guide for kirigami-based structures.

Graphical Abstract

KEYWORDS:

• Kirigami design
• compliant mechanism
• flexural hinge
• mechanical metamaterials
• buckling instability

1. Introduction

Kirigami is a paper cutting art that originated in China as early as 1500 years ago and then spread into other Asian countries [1,2]. With these distributed cuts, the stretchability of a paper-like material is greatly improved. Moreover, paper-like materials exhibit a negative Poisson’s ratio character of becoming wider when stretched because of the introduction of these distributed cuts [3–5]. Although there are a variety of fabrication approaches such as 3D printing [6] and 4D printing [7,8] to manufacture mechanical metamaterials, most of them are costly and time-consuming. The emergence of kirigami principle offers a simple and effective approach to design metamaterials toward programmable mechanical properties [9]. Studies demonstrate the capability of the kirigami metamaterials, such as extremely large strain [10,11] and complex 3D structures [12,13], across length scales from macro to nano [14,15] in materials ranging from monocrystalline silicon to plastic [16,17]. Because of their unique behaviors, many innovative designs have been developed, such as structures with programmable patterns and shapes [18], soft actuators and robots [19–24], mechanical logic gates [25,26], wearable devices [27,28], stretchable electronics [29–31], and energy harvesting devices [32–34]. Kirigami also inspires many design methods of mechanical metamaterials with extraordinary response, and one of the most consequential examples is the wave propagation bandgap. The ability to manipulate propagating waves in a periodic medium has developed many applications in wave transmission control [35,36] and vibration attenuation [37].

The cuts in a kirigami works as localized compliant hinges whose deformation in terms of extension, torsion, or twisting provides a set of diverse motions in the developed metamaterials [38–41]. Although it has been found that decreasing thickness leads to instability when the perforation patterns are subjected to tensile loading, most of the previous studies have considered thick sheets and focused on their in-plane deformation [42–44]. Mechanical responses of thin kirigami sheets have been investigated, including the in-plane stiffness and critical buckling strain [45]. Moreover, through combining different hierarchical patterns, mechanical responses of kirigami have been further enriched, which makes encrypting information into kirigami structures possible [46]. However, conventional cuts in a crack-shape were widely employed in previous research. With the narrow width, cuts are easy to produce crack propagation or large plastic deformation in stretching, thus leading to a mechanical failure of the metamaterial. Even though studies indicate that plastic deformation significantly affects the response of the kirigami structures [46], improved design has not been proposed. Here, we design the cuts in the shape of a compliant beam, which works as a flexural hinge. Compliant beam is a typical element in the category of compliant mechanism, which serves as a flexure hinge undergoing elastic deformation relative to adjacent stiffer regions in compliant mechanisms [47]. Instead of a conventional revolute with three components, a compliant beam is capable of mechanical bending motion with a single component only, as well as other advanced features of frictionless, assembling-free, and tolerance clearance-free [48,49]. Normally, the difference of compliance and stiffness is reached by geometric characteristics of deformation regions [50]. Thus, compliant beams have been employed in origami for a new folding joint design that can effectively accommodate the thickness of origami panels [48]. In our study, with a designed cut width in a notch-shaped flexure, each kirigami cut is actually a compliant beam that can bend or twist, which elongates the mechanical strength in the kirigami structures.

Therefore, this study focuses on kirigami sheets comprising of a square array of mutually orthogonal cuts with compliant beam design. This pattern not only reduces stress concentration at the tips of cut but also provides a new way to tune the mechanical response of kirigami. We investigate both the in-plane and out-plane deformation process of this typical kirigami structure with wider cuts. Experimental results validate the analytical and numerical prediction in the tunable mechanical property. Mechanical metamaterials are studied utilizing the kirigami unit cell, illustrating its capability in tunable mechanical property in vibration isolation.

2. Mechanics of complaint beams in kirigami cut

Figure 1(a) shows the difference of a conventional cut and the new design. A conventional cut in the crack-shape will propagate during stretching, leading to a fracture. In the new design, we extend the width of cuts into a notch and model them in the shape of a compliant beam, which works as a flexural hinge and can withstand more mechanical stretching load. Perforation pattern shown in Figure 1(a) introduces a network of square domains of edge $L$ separated by hinges of width $\delta$. Utilizing compliant beam as the notch, we add another geometrical parameter, the width of cuts $2l$, and study its effect on the mechanical behavior. We first study the initial in-plane linear elastic response of the system. In all our calculations, we assume that all deformation is localized at the flexural hinges and the square domains are rigid, which is consistent with the concept of lumped compliant mechanism. Here, we derive an analytical relation for the effective Young’s modulus of the perforated sheets $\bar{E}$, the geometrical parameters $L$, $\delta$, and $2l$, and the Young’s modulus of the sheet $E$.

Figure 1. Kirigami cut design using a compliant beam as a notch flexure. (a) Comparison of conventional cut and compliant beam cut. With the new design, the undeformed and planarly deformed configurations are illustrated. (b) Two types of buckled morphology at different level of applied deformation for perforated sheets with $2l/L=0.1$ and 0.2.

We focus on a unit cell comprising four square domains and stretch it uniaxially along one set of cuts by applying a macroscopic stress:

(1) ${\sigma }_x=\frac{f}{\left(2L+4l\right)t}$(1)

where $f$ is the force applied to the flexural hinges on the vertical boundaries and $\left(2L+4l\right)t$ denotes its cross-sectional area. The applied uniaxial stress generates identical bending moments $M$ at all hinges, which in turn induce the rotation of all square domains by an angle ${\theta }_i$ and the opening of all the cuts by an angle $2{\theta }_i$. Focusing on a single square domain, it is easy to see that

(2) $M=\frac{fL}{2}\times \frac{1}{4}$(2)

Since $M$ must be balanced by the couple induced by the internal loads, we have

(3) $M=\frac{EI}{\rho }$(3)

where $I_i=t{\delta }^3/12$ is the second moment of area of each hinge about the z axis and ${\rho }_i$ denotes the curvature of each bent hinge. When the width of a cut is very narrow, the length of the bent region is assumed as the hinge width $\delta$. However, we noticed that if the cut becomes wider, the effect of rigid squares becomes insignificant. Assuming that when $2l$ increases from 0 to $L/2$, contribution of rigid squares decreases linearly. When cut width reaches $L/2$ and continues to increase and bent region is limited to the hinge area, we obtain

(4) ${\rho }_i\simeq \frac{\left(1-2\delta /L\right)2l+\delta }{2\theta }$(4)

Substitution of Eq (2)(2) $M=\frac{fL}{2}\times \frac{1}{4}$(2) and Eq (3)(3) $M=\frac{EI}{\rho }$(3) into Eq (4)(4) ${\rho }_i\simeq \frac{\left(1-2\delta /L\right)2l+\delta }{2\theta }$(4) yields

(5) ${\sigma }_x=\frac{8EI\theta }{\left(L+2l\right)\left(2l\left(L-2\delta \right)+L\delta \right)t}$(5)

Moreover, since the strain in the loading direction is given by

(6) ${\epsilon }_x=\frac{L\theta \left(cos\theta +sin\theta \right)+2lsin\theta }{\theta \left(L+2l\right)}-1$(6)

In the small deformation regime, we have

(7) ${\epsilon }_x\simeq \frac{L\theta }{L+2l}$(7)

It follows that the effective Young’s modulus of kirigami sheet $\bar{E}$ is

(8) $\bar{E}=\frac{{\sigma }_x}{{\epsilon }_x}=\frac{2E{\delta }^3}{3L\left(2l\left(L-2\delta \right)+L\delta \right)}$(8)

When the out-of-plane deformation of the flexural hinges becomes energetically less costly than their in-plane deformation, buckling of the beams occurs. We define the critical strain ${\epsilon }_c$ at which the out-of-plane deformation is triggered. Then strain energy density can be expressed as $U_i\left({\epsilon }_c\right)+U_o\left({\epsilon }_x\right)$, where $U_i$ and $U_o$ represent the in-plane and out-of-plane strain energy densities, respectively. Energy contribution due to the out-of-plane deformation is not significantly affected by cut width, and similar expression of out-of-plane deformation energy previously obtained for narrow cuts [45] can also be used:

(9) $U_o\left({\epsilon }_x\right)=\frac{E{t}^2{\theta }_o^2}{3{\left(L+2l\right)}^2}$(9)

Where ${\theta }_o$ denotes the opening angle of a square after out-of-plane buckling. Finally, we find that

(10) ${\epsilon }_c=\frac{\left(2l\left(L-2\delta \right)+L\delta \right)L{t}^2}{2{\left(L+2l\right)}^2{\delta }^3}$(10)

3. Experimental validation and FEM study

To validate the prediction of the theoretical model, we characterize the response of our kirigami structures using a mechanical testing machine (WANCE ETM502A). The material considered in this study is polyethylene terephthalate (PET), and its property is obtained by uniaxial tensile test with a 200 N load cell. Samples with a gauge length of 120 mm and width of 10 mm were fully clamped at both ends and stretched with a displacement rate of 6 mm/min up to $\epsilon =0.1$. The stress–strain curve in Figure 2(b) suggests that PET has a Young’s modulus of 4.57 GPa. Moreover, 0.2% offset yield strength is measured as ${\sigma }_y=86.8\mathrm{M}\mathrm{P}\mathrm{a}$, and their plastic strain ${\epsilon }_p$ versus $\sigma$ is plotted up to fully plastic region. Finally, a typical Poisson’s ratio $\nu =0.4$ is assumed for all the plastic sheets. The specimens considered in this study are fabricated by laser cutting an array of $2\times 7$ mutually perpendicular cuts into PET as shown in Figure 2(a). First, we fabricated a line-shaped marker on the samples by laser cutting and measured the width of the cuts using an optical microscope. Multiple measurements show that the thinnest cut width is about 0.2 mm. Then, we take 0.2 mm as the design width of narrow cuts for comparison. The settings of laser cutting can be determined to minimize errors of the cut widths as much as possible. In Figure 2(c), we present results of samples fabricated with different cut widths $2l$ (0.2 mm, 1 mm, and 2 mm). All specimens have the same length ($L=10\mathrm{m}\mathrm{m}$) and thickness ($t=0.075\mathrm{m}\mathrm{m}$). These results are obtained by uniaxial tensile test equipped with a 100 N load cell.

Figure 2. Experimental and FEM study of the kirigami sheets. (a) Experiment setup and view of the perforated sheet. (b) Stress–strain curve of the plastic sheet under uniaxial tension. (c) Snapshots of the sample with $2l/L\in \left[0.02,0.2\right]$ at $\epsilon =0.1$, 0.15 and 0.2. (d) Quantitative comparison between experimental and numerical stress–strain curves. Results for different patterns with $2l/L\in \left[0.02,0.2\right]$, $\delta /L\in \left[0.08,0.12\right]$.

We conduct finite element (FE) simulations using the commercial package ABAQUS/Standard 2017 to further investigate the mechanical stretchability. In all simulations, the models are discretized with 3D shell elements (S4R). Narrow cuts in the sheets are modeled as thin rectangular slits with a thickness of 0.2 mm as observed in the physical samples. Wider cuts are modeled as rectangular slits with design parameters. The material behavior of the sheets is captured using an elastoplastic model with parameters directly extracted from uniaxial tensile tests conducted on the polyester plastic sheets (i.e. Young’s modulus $E=4.57\mathrm{G}\mathrm{P}\mathrm{a}$ and yield stress ${\sigma }_y=86.8\mathrm{M}\mathrm{P}\mathrm{a}$). Specifically, we take the plastic strain versus stress for the plastic sheet according to Figure 2(b). An elasto-plastic model (material models *ELASTIC and *PLASTIC in ABAQUS) is used with the experimentally characterized properties. Moreover, to understand the effect of plastic deformation, we also compare our results to those obtained by using a purely elastic material model with $E=4.57\mathrm{G}\mathrm{P}\mathrm{a}$ and $\nu =0.4$.

To save the computational cost and eliminate the boundary effects, we use representative volume elements (RVEs) to simulate the 2D infinite periodic cut structures. Periodic boundary conditions (PBCs) are applied as displacement boundary conditions on the nodes of the edges while the uniaxial loads are controlled by forces on the nodes of the edges. For a specific RVE, PBCs are applied to the left and right edges, where the displacements along x and y axis on the two vertical edges are kept the same. Uniaxial loads along y-axis are applied to the top and bottom edges, where the displacements along x-axis are kept the same during deformation. On the meshed RVEs, local fine meshes are applied to all the hinge areas since these regions are highly deformed with severe stress concentration. The in-plane linear elastic response of the sheets is simulated by conducting static simulations. For out-of-plane response simulations, we use buckle module in ABAQUS to identify the critical buckling modes. After identifying the critical buckling mode, we then introduce a small imperfection in the form of the critical mode into the mesh to guide the post-buckling analysis through static/Riks algorithm in ABAQUS. We do not need to modify the mesh directly, instead, a result file will be written by buckle module. In the next step, this file containing buckling results of the mesh can be read in post-buckling analysis; meanwhile, we define the value of the small imperfection. In this way, we can obtain the result corresponding to any strain in the post-buckling process. In order to control the strain of the structures more precisely, we also use dynamic/implicit algorithm in ABAQUS, which allows for displacement controlling. Results of the two algorithms are consistent. To validate our unit cell simulations, we compare their results to those of the experimental results. Numerical snapshots for structures using elastoplastic model are reported in Figure 1(b). The critical strain values of the first two buckling modes are almost equal. Either of the two can be chosen to influence the direction of the instability in post-buckling analysis. In the experiment, these two types of buckled patterns occur randomly, which agrees with the simulation results.

As shown in Figure 2(d), we find excellent quantitative agreements between effective Young’s modulus, critical strain, and stress–strain curves, with small discrepancies caused by boundary effects. According to previous observation [45,46], the stress–strain curves can be divided into three stages. At the first stage, all hinges bend in-plane, resulting in a linear response. After then, the hinges generate out-of-plane buckling deformation and cause the stress a sudden departure from linearity to a plateau state, which can be considered as the second stage. Finally, all unit cells are stretched to the maximum extent, and an extra tension to the kirigami structures results in a rapid increase in stress while the strain remains almost unchanged.

We also find excellent agreement between the numerical and theoretical results, as shown in Figure 3(a). Results collected in Figure 3(b) and Figure 3(c) indicate that the stiffness of the perforated sheets is reduced when the width of cuts is increased. As a result, wider cuts lead to an increase in critical strain.

Figure 3. Response of kirigami sheets subjected to uniaxial tension. (a) Comparison between numerical and theoretical effective Young’s modulus (with $\delta /L=0.1$) and critical strain. (b) Effective Young’s modulus of perforated sheet (with $2l/L\in \left[0.02,0.2\right]$,$\delta /L\in \left[0.08,0.12\right]$). (c) Critical strain of perforated sheet (with $2l/L\in \left[0.02,0.2\right]$,$\delta /L\in \left[0.08,0.12\right]$).

4. Plastic deformation in the compliant beam

After verifying the validity of simulations, we examine the effect of plasticity and investigate the response of both elastic and elastoplastic kirigami sheets with cut width $2l/L\in \left[0.02,0.2\right]$. The results reported in Figure 4(a) indicate that with further deformation after the initial linear regime, plastic deformation significantly affects the stress–strain response of the structures at the second deformation stage. Since the large local strain is expected to induce plastic deformation in most materials, it is important to reduce stress concentration by increasing the load bearing area. We then parametrically investigate the effect of geometric parameters (i.e. $2l/L$, $\delta /L$, and $t/L$) on the mechanical response of the system. We simulate the response of kirigami unit cells with $2l/L\in \left[0.02,0.2\right]$, $\delta /L\in \left[0.08,0.12\right]$, and $t/L\in \left[0.0075,0.0188\right]$. The results of this parametric study are displayed in Figure 4(b). For all the considered kirigami sheets, the specific characteristics of the stress–strain curves can be significantly tuned by altering $2l/L$, $\delta /L$, and $t/L$. The results show that our kirigami sheets are capable of constructing a set of buckling-induced 3D morphology through regulating geometric parameters and the applied strain.

Figure 4. Effect of material behavior on the response of kirigami sheets subjected to uniaxial tension. (a) Comparisons between numerical stress–strain curves obtained considering an elastoplastic (solid line) and a purely elastic (dashed line) material model. (b) Effect of the geometric parameters on the stress–strain curves of the elastoplastic kirigami sheets.

Next, we analyze the stress distribution in the structures through FEM simulations using elastoplastic and purely elastic model. Figure 5(b) shows the corresponding simulated von Mises stress contours of kirigami structures at the same stretching strain. By monitoring the distribution of the von Mises stress within the sheets, we find that plastic deformation initiates at the tip of hinges at the onset of buckling and then gradually expand to fully cover the hinges when the sample is fully stretched and the deformation mechanism changes from bending dominated to stretching dominated. Figure 5(a) shows the plastic strain occurs and increases as the sheet is stretched. We find that increasing the width of cuts can delay the generation of plastic strain at the hinge tips while the strain of the sheets is similar, and the value of plastic strain is significantly reduced in the following process. As comparison with the plastic strain model, the maximum von Mises stress at the tip of hinges using the purely elastic material property. It is clear that increasing the width of cuts reduces the stress concentration. Finally, we analyze the strain components at the hinge tips. In the case shown in Figure 5(c), all hinges bend in-plane at zone I, resulting in a linear response. After buckling (i.e. ${\epsilon }_c=0.004$), there is still no plastic deformation at zone II, and the yielding of the hinges begins approximately at $\epsilon =0.04$. This strain is much greater than the critical strain. Then, strain of the sheet changes from elasticity dominated (zone III) to plasticity dominated (zone Ⅳ), which reduces the stretchability of the structure. Our simulations confirm that when the width of cuts is appropriately increased, 3D reconfigurable space of the kirigami sheet is guaranteed.

Figure 5. Increase the width of cuts to improve the strength and stretchability of kirigami structures. (a) Maximum principal plastic strain using an elastoplastic material model and maximum von Mises stress using a purely elastic material model. (b) Simulated von Mises stress contours for elastoplastic sheets with $2l/L\in \left[0.1,0.2\right]$; samples show stress concentration in all the hinge areas. (c) Maximum principal strain components using an elastoplastic material model.

5. Illustration of phononic crystal metamaterials

The ability to tune wave propagation and filter specific bands of frequencies, via Bragg scattering in a periodic medium, has enabled countless applications. Using the kirigami principle is a simple, scalable, and effective approach for creating and controlling periodic structures. The kirigami-based metamaterial can be utilized as a phononic crystal when in-plane strain is induced [11]. For thin plastic sheets, it is very easy to buckle during stretching, which leads to more complex pattern transformation. We focus on the tunable Bragg bandgap characteristics in their buckling-induced 3D morphology. Bloch-wave analysis is conducted using COMSOL Multiphysics 5.6 Solid Mechanics Module to compute the phononic band structure of the metamaterials with different hinges. Bragg bandgap is mainly influenced by material properties and structure parameters. Perforated sheets with $t=0.3\mathrm{m}\mathrm{m}$ are used in order to ensure the stiffness, and the density of the PET is $\rho =1400\mathrm{k}\mathrm{g}/{\mathrm{m}}^3$. Deformed RVE calculated from ABAQUS is imported into COMSOL to conduct eigenfrequency studies, which ensure the external 3D shape of a buckled kirigami sheet. The two buckling forms with different directions of the instability are both taken into consideration. Moreover, the stress distribution is concentrated almost exclusively in hinge regions, and the stress level can be controlled. Compared to the shape of kirigami sheets, stress have little effect on the bandgap characteristics.

Figure 6(a) shows an example of unit cell composed of our rigid squares and connected hinges. It should be pointed out that although the structures we present are 3D after buckling, the periodicity direction is along the y axis. RVEs are periodically arranged with a lattice constant of $a=\left(2L+4l\right)\left(1+\epsilon \right)$, where $\epsilon$ is the strain in the loading direction controlled by boundary conditions. The Irreducible Brillouin Zone of the lattice is $\left[-\pi /a,\pi /a\right]$ and $k$ is the wave number within it. Via COMSOL, floquet boundary condition in the y-axis direction is applied to the deformed RVEs. The phononic bandgap is obtained by sweeping the eigenfrequency analysis. Here, only the wave propagation in the y-axis direction is considered.

Figure 6. The numerical Bragg’s elastic wave propagation bandgaps in the periodic kirigami metamaterials. (a) The periodic kirigami metamaterial with a lattice constant of $a$. (b) Results for the deformed structure with $2l/L\in \left[0.02,0.2\right]$. (c) Results for the deformed structure with $\delta /L\in \left[0.08,0.12\right]$. (d) Results for the deformed structure with different buckling patterns with $2l/L=\delta /L=0.1$. (e) The displacement distribution in the first-order buckling deformed finite structure corresponding to the case in (c) at 500 Hz and 2700 Hz. (f) Result for the in-plane deformed structure with $2l/L=\delta /L=0.1$.

We create geometry from shell elements and rebuild Physics-controlled mesh in COMSOL to simulate 3D structures. We refine the mesh to ensure that the result of bandgap is almost independent of the element size (no more than 0.5%). Figure 6(b) -(d) shows the comparisons of numerically calculated band diagrams and forbidden bandwidths between buckled kirigami metamaterials with different geometric parameters ($L=10\mathrm{m}\mathrm{m}$) of hinges at the same stretching strain ($\epsilon =0.1$). Phononic bandgaps of forbidden frequencies, where waves could not propagate through, are denoted by the gray filled areas. There are two obvious Bragg bandgaps in the results of each structure. Wider cuts induce lower frequencies of bandgaps. Moreover, the bandgaps found in the structure with wider cuts are narrowed. When hinge width is increased, two bandgaps approach to each other. In addition, if we adjust the two buckling patterns through opposite initial perturbation, different bandgaps will be formed, which can be further confirmed by the transmission spectra obtained from the numerical simulation as shown in Figure 6(d). Wider range of bandgap tunability in 3D morphology provides the capability to guide elastic wave. To further demonstrate the wave propagation, we plot the displacement distributions of the finite structure in Figure 6(d) under the first state (frequencies marked with red circles). As shown in Figure 6(e), on the one hand, 500 Hz is on a passband of the guided wave mode, and thus, the elastic waves can propagate. On the other hand, 2700 Hz is in the range of the guided mode bandgap of the deformed structure, and hence, the wave propagation is suppressed at this frequency. At last, we plot the band diagram for in-plane deformed structure as a comparison in Figure 6(f). There is no bandgap in the same frequency range.

6. Conclusions

In conclusion, by taking into account the width of cuts, we design the kirigami cut in terms of a notch flexure with compliant beam instead of a conventional crack-shaped cut. The mechanical responses with different design parameters are studied experimentally, analytically, and numerically. With the introduction of cut width in the compliant beams, the flat kirigami sheets improve the repeatable configuration ability with tunable buckling conditions. Programmable onset of buckling-induced 3D patterns can be predicted and tuned by carefully choosing a series of geometric parameters of the sheets. Moreover, by finite element simulation, we reveal the propagation of plastic deformation from hinge tips, which can be delayed by increasing the width of cuts. Through rational design of the local hinges, we obtain a highly stretchable kirigami metamaterial with enhanced structure strength and controllable deformation as phononic crystal. By the control of hinge geometry and buckling mode in the metamaterials, we are able to filter undesirable frequencies. The insights presented here shall have a broad range of potential applications of mechanical metamaterials, including, but not limited to, noise mitigation, wave transmission and vibration control, acoustic cloaking, non-reciprocity, and even wave-based mechanical computation, with precisely controlled material stiffness and mechanical response.

Acknowledgments

This work was supported by the National Key R&D Program of China (2019YFB1311600), Natural Science Foundation of China (Grant No. 52075411 and U1913213), Shanxi Key Research and Development Program (2020ZDLGY06-11). B. L. thanks XJTU for joint-collaboration project in multi-disciplines (xhj032021014-03).

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was supported by the XJTU joint-collaboration project in multi-disciplines [xhj032021014-03]; National Key Research and Development Program of China [2019YFB1311600]; Shanxi Key Research and Development Program [2020ZDLGY06-11]; and Natural Science Foundation of China [No. 52075411 and U1913213].