Abstract

To clarify the influence of the unbalanced coefficient and the change in lay-up angle on the failure characteristics of the laminate in the static aeroelastic problem of the aircraft, numerical simulations were performed based on the classical laminate theory and the Tsai-Wu failure criterion, as well as the fluid-structure coupling calculation method. The structure’s stress-strain and failure curves are found to decrease as the unbalanced coefficient increases. The stress curve’s slope is relatively stable, whereas the strain and failure curves’ slopes change three times, indicating that strain may be the primary cause of structural failure. Unbalanced coefficient laminates are classified into three types based on their mechanical properties low unbalanced coefficient laminates (Unbalanced coefficient 0.2 to 0.3), quasi-balanced coefficient laminates (Unbalanced coefficient 0.4 to 0.6, Balanced laminate when the Unbalanced coefficient is 0.5), and high unbalanced coefficient laminates (Unbalanced coefficient 0.7 to 0.8). Within their respective spanning intervals, the mechanical properties of the three types of laminates remain relatively stable. An increase in the ply angle reduces both the elastic deformation of the structure and the failure factor. The variation patterns of structural strain and failure at 45° and 60° ply angles decrease as unbalanced coefficients increase, whereas the opposite is true for 30° ply angles. Finally, a two-level optimization method based on “equalized plies” and “equal-angle plies” was developed, resulting in a 23.93% reduction in elastic deformation and a 37.04% reduction in the laminated structure's failure coefficient when compared to the preoptimization results.

1. Introduction

Composite materials as an emerging technology in the last few decades have developed rapidly in the engineering field. Composite materials offer many advantages over traditional metal materials, including high specific strength, high specific stiffness, good corrosion resistance, and directional desirability. Take the aerospace industry as an example, the composite structure of large aircraft in the civil aviation market accounts for about 50% of the fuselage structure. In military aircraft, the number of composite materials in the world’s most advanced military aircraft accounts for about 20% to 50% of the total structural weight of the aircraft [13]. In the aircraft design phase, to ensure the flight performance and structural performance of the aircraft, researchers have proposed aircraft structural shear design, which provides various designs and optimization of composite structure or lay-up parameters. Aircraft design engineers can achieve a more robust and lightweight aircraft structure by modifying the composite layer design based on the actual situation. It may be argued that the application of composite science and technology has fostered the rapid development of modern product structure design. At the same time, with the massive use of composite materials in aircraft structures, in addition to the conventional constraints such as structural strength and stiffness of products, new problems such as structural failure of composite materials should also be considered in the design phase of aircraft structures.

To further develop the theory of failure and structural optimization of composite materials, much of the fundamental research has been done by relevant researchers and experts. In the exploration of the failure law of composite materials, Huang and Zhu [4] took the aircraft composite long truss as the research object and analyzed the out-of-plane failure of the long truss based on the Hoffman material failure criterion and using the finite element method. Wang and Kehua [5] investigated the gradual failure of an open curved composite plate and obtained the process of extending the damage near the opening. Qingquan et al. [6] from Northwestern Polytechnic University focused on the analysis of progressive damage and failure of 2.5 D woven composites and accurately predicted the stress-strain response of material singlet under compressive loading by experimental and finite element simulations and were able to simulate the damage initiation and evolution processes in woven structures. Based on the continuous damage mechanics and bonded cell model, Li et al. [7] conducted a low-velocity impact and postimpact compression simulation analysis of laminates using ABAQUS software to analyze the impact damage and compression failure behavior of laminates. Luo et al. [8] established a fine mechanic's representative volume element model based on the macrofine mechanics cross-scale analysis method and calculated the failure envelope of the composite. Cadieu et al. [9] investigated the effect of temperature variables on the damage of composite laminates. Caminero et al. [10] analyzed the effect of layer thickness and layer sequence on fracture and damage characteristics.

Others have focused their research on the study of the mechanical and structural properties of laminated plates. For example, Belbachir et al. [11] used the four-variable optimized plate theory and proposed an improved theory of laminates to study the buckling of antisymmetric cross laminates under nonlinear thermal and force loads, and the original variables were reduced in order to derive and validate the analytical solutions for simply supported plates. Abualnour et al. [12] analyzed the effect of structural parameters on the performance of composite laminates and Sandwich panels using a four-variable triangular shear deformation model. Draiche et al. [13] developed an analytical model for the static force analysis of composite laminates under sinusoidal uniform load based on an improved first-order shear deformation theory and compared the results of other plate theories to prove the accuracy and efficiency of the theory. These theories may bring new ideas for the structural calculation of composite laminate plates.

Many methods and theories have also emerged in the structural optimization of composite laminates. Ding and Xu [14] proposed a new discrete-continuous optimization model for parallel optimization of the macrostructure topology and local fiber angles of composite structures. Deveci and Artem [15] investigated the problem of estimation and optimization of multi-axial fatigue life prediction models. The fatigue life problem was predicted for graphite/epoxy, glass/epoxy, carbon/epoxy, and carbon/PEEK composite laminates. Subsequently, a hybrid algorithm combining the particle swarm algorithm and generalized pattern search algorithm was utilized to optimize the lay-up sequence of the laminate materials. Kaveh et al. [16] used several hybrid algorithms to solve the computation with the maximization of flexural load as the optimization objective, and finally, by comparing the optimization of various algorithms, it was determined that the optimization capability of the JAYA algorithm, in this case, Cui and Li [17] calculated the general analytical conditions for laminates under multiple coupling conditions. Then four laminates with bending and torsional coupling structures were designed and optimized. Xu et al. [18] used the relative density as a design variable to describe the distribution of different materials, and the sum of the stiffnesses of the intermediate structures at a specific degradation time step as the optimization objective to optimize the topology of the composite microstructure configuration.

In summary, among the three types of research directions in composite laminates, the first type is based on some classical laminate failure theories (e.g., maximum stress criterion, Tsai-Hill criterion, Tsai-Wu criterion, Hoffman criterion, and Hashin criterion) using experimental or numerical simulation methods to analyze the failure of composite structures under multiple working conditions and various engineering contexts. The second category is the development and analysis of computational methods for laminates under various loads (thermal, nonlinear, fluid loads, etc.), which will provide some theoretical basis for optimal laminate design. The last category is the use of existing mature optimization algorithms, based on single or hybrid algorithms for the optimization of the objected objective in question. In this paper, we first discuss and simulate the failure effects of different lay-up angles and different lay-up unbalanced coefficients on the wing structure, and then develop a two-stage optimization method to optimize the laminate structure parameters based on this background. The results show that this optimization idea can effectively improve the maximum deformation and failure coefficient of the unbalanced composite structure wing. The first part of the article introduces the theoretical solution method; the second part describes the establishment of the finite element model of the unbalanced laminate of the wing and the detailed design of the lay-up scheme; the third part analyzes the specific effects of each variable on the structural failure of the laminate for different unbalanced coefficients and different lay-up angles. The two-level optimization method is also used to optimize the parameters of the laminate structure.

2. Computational Theory

For the static aeroelastic problem, the structural deformation is considered to be slow, so the inertial and aerodynamic loads caused by the structure and other loads and aerodynamic forces acting on the aircraft are a tiny amount compared to the negligible law. Then the static aeroelasticity calculation does not consider time as an independent variable, but directly uses alternating iterations of structural deformation and constant aerodynamic forces until convergence. The calculation process is shown in Figure 1. In the calculation, the initial conditions are first set to solve the aerodynamic force. When the aerodynamic force converges, the aerodynamic load is transmitted to the grid node to solve the structural deformation. Finally, the structural deformation is fed back to the flow field to recalculate the aerodynamic force. The above process is iterated repeatedly until the wing structural deformation converges.

2.1. Aerodynamic Calculation

The Navier–Stoke equation is used as the control equation to solve the flow field, and the integral form in the Cartesian coordinate system can be expressed as follows [19, 20]:where is the density of the grid center; is Reynolds number; are the velocities in directions; is pressure; is the total energy; is the stress tensor; is the heat flux; as the enthalpy; is the flow field conservation variable matrix; is the convection flux matrix; is the dissipated flux matrix; is the gas microcluster control body; is the boundary of the control volume unit; is the control volume element; and is the external normal area vector of the area element.

2.2. Structural Calculation

The structural equilibrium equation for the static aeroelastic calculations can be expressed as follows [21, 22]:

For linear deformation problems, the external load and stiffness matrices are determined before calculation, and the external load and structural deformation are linearly related by solving the problem only once to derive the structural deformation. When a geometrically nonlinear problem involving large structural deformation is involved, the stiffness matrix can be expressed as follows:where is the force vector; is the main tangential stiffness matrix; is the large displacement stiffness matrix; and is the initial load matrix. For the set of equations described above, the principal tangential stiffness matrix is related to the initial configuration of the structure; the large displacement stiffness matrix reflects the changing characteristics of the structure geometry; and the initial load matrix represents the changing load direction. When the equations are solved, the structural stiffness matrix is constantly iteratively corrected. Therefore, a stepwise loading approach is required to solve the equations. In this paper, the Newton–Raphson method is used to achieve the solution of the nonlinear system of equations.

For the composite wing lay-up problem, it is required that according to the classical laminate theory, the internal force and internal bending moment acting on the surface of the plate for the whole laminate can be expressed as follows:

According to the forces shown in Figure 2 for the laminated plate, where the internal forces and moments are of the following specific form:

After integration, the matrix form can be expressed as follows:

Among them,where is the tensile stiffness matrix, is the coupling stiffness matrix, and is the bending stiffness matrix. We use for out-of-plane; the displacement is called the bending curvature of the laminated plate.

2.3. Failure Calculation

The Tsai-Wu failure criterion, as a relatively mature failure theory for composites, is more predictive of the failure of composite laminates. Its specific form can be expressed as follows [23]:

Among them,where are the tensile and compressive strength of the fiber in the longitudinal direction; are the tensile and compressive strength in the transverse direction; is the in-plane shear strength.

In practical engineering applications, the value of can be simplified accordingly. In combination with references [24, 25], the value of is made to be 0. The final laminate failure can be judged by the sum of the left-end terms of the equation and whether the value is within the interval [−1,1] by solving the equation.

2.4. Optimization Calculation

In the optimization design, the maximum force deformation of the laminate is taken as the optimization objective, the constraint is the failure strength constraint, and the optimization variables are the lay-up angle and the lay-up sequence. With reference to the actual process, the lay-up angles are discrete as ±30°, ±45°, ±60°, 0°, and 90°. Thus, the optimization model of unbalanced laminate can be abstracted as the following mathematical models:where equation (10) is the objective function, which is the displacement of the flat plate after the force in this paper; equation (11) is used to define the inequality constraint, here is the maximum failure factor; equation (12) is used to specify the upper and lower boundaries of the angle of the discrete pavement.

Usually, multiple layers of lay-up angles are set in laminates, and there are many variables of lay-up angles and lay-up sequences, which will make the optimization work very difficult. Therefore, this paper develops a two-level optimization method based on the Screening method to optimize the lay-up angle and sequence of the laminate. In the first-level optimization, all the laminates with the same angle [] are defined as one “balanced layer,” so that each laminate can be divided into several “balanced layers,” and each “balanced layer The angular variation of each “equalization layer” is ±30°, ±45°, and ±60°. The remaining single laminates that are not defined as “equalization layers” are divided separately and their lay-up angles can be defined as ±30°, ±45°, ±60°, 0°, and 90°. Therefore, in the first-level optimization, if the “balanced layer” appears in pairs and the angles of the remaining plies are only 0° and 90°, then the whole laminate is a balanced laminate; if there are other plies besides the “balanced layer” with the lay-up angles of ±30°, ±45° or ±60°, then the whole laminate is a balanced laminate.

The lay-up angles of composite laminates can be obtained according to the first-level optimization. On the basis of the first-level optimization, the layers with identical lay-up angles are then integrated into one “equal-angle layer,” which can be divided into eight types of “equal-angle layers” depending on the lay-up angle. They are “0° equal-angle layer,” “30° equal-angle layer,” “30° equal-angle layer,” “45° equal-angle layer,” “−45° equal-angle layer,” “60° equal -angle layer,” “−60° equal-angle layer,” and “90° equal-angle layer.” Next, each “equal-angle layer” is sieved, iterated, and selected as the best lay-up angle. For the optimization of the ply order, the final ply angles already screened are directly sorted and optimized based on the secondary optimization. In summary, an optimal deformation solution can be obtained after two-level of laminate design optimization.

3. Computational Model

The flat wing model used in the example has a half-spread of 600 mm, and the airfoil is NACA 4412. The structure of the wing is the classic double-beam multi-rib, upper and lower skin structure. The incoming flow characteristics are as follows: the flight speed at sea level is 0.5 Mach number, and the angle of attack of the wing is 1.5°. Among them, the wing beam and wing ribs are made of aluminum alloy with a yield limit of 280 MPa, safety factor of 2, Young’s modulus of elasticity of 71 GPa, and Poisson’s ratio of 0.33. It is shown in Figure 3.

On the premise of not changing the overall stiffness and affecting the calculation accuracy, other structures on the wing skin are simplified accordingly, and the solid shell is used for geometric modeling. Among them, T700 carbon fiber material is selected for the upper and lower layers of the skin, and A550 composite material is selected for the central core. The specific material parameters are shown in Tables 1 and 2.

In the initial structural design, in order to investigate the failure effect of the laminate with different lay-up angles and different unbalanced coefficients. The lay-up scheme is defined to be the same for both upper and lower panel lay-up areas, and the unbalanced coefficient in the whole lay-up scheme can be adjusted by subsequently changing the ratio of −α to α in the lay-up sequence. Here the values of α are 30°, 45°, and 60°, respectively. The 0° fiber reference direction is set to the Z-axis direction for the pavement design. The layer thickness of T700 carbon fiber material is 0.1 mm per layer, and the thickness of A550 core material is 1 mm. The specific initial lay-up scheme is shown in Table 3. The final laminate lay-up geometry and lay-up direction are shown in Figure 4.

In the flow field solution, the grid type is the structural grid, and the total number of grids is 1,057,176. The Spalart–Allmaras one-equation turbulence model [26] was used in the flow field solution, and the control equations were discretized using the finite volume method. In the spatial discretization term, the viscous flux vector is chosen in the central format and the convective flux vector is chosen in the Roe-FDS format. In the time discretization term, the LU-SGS implicit time discretization method is used to advance the solution [27, 28]. The boundary conditions at the wing root are set as symmetric boundary conditions, the wing surface is set as no-slip boundary conditions, and the rest are set as pressure far-field boundary conditions. The flow field calculation domain is shown in Figure 5.

For structural calculations, the laminate structure is divided using a hexahedral mesh with a total mesh count of 162,885. The laminate structure is divided using a hexahedral mesh. The fixed constraint boundary condition is set at the wing root. Subsequently, the failure of the composite can be deduced with known structural deformation. The structural limit element model is shown in Figure 6.

4. Result Analysis

4.1. Effect of Unbalanced Coefficients on Structural Failure

Taking the normal ±45° lay-up angle as an example, Figure 7 shows the stress, strain, deformation, and the relative failure distribution of the laminate when the unbalanced coefficient increases from 0.2 to 0.8. From the figure, it can be seen that the maximum stress and strain of the wing structure are mainly concentrated at the wing root, and the stress and strain show a decreasing trend as the unbalanced coefficient increases. The difference is that the slope of the stress curve is relatively stable, while the slope of the strain curve changes three times. The inflection point of the first change occurs near the unbalanced coefficient of 0.3, the slope change of the second curve occurs at the unbalanced coefficient of about 0.45, and the slope change of the third curve occurs at the unbalanced coefficient of 0.6. The fastest decreasing interval is mainly concentrated in the stage of the unbalanced coefficient of 0.4 to 0.6. This indicates that when the laminate changes from unbalanced to balanced laminate, it will have a major effect on the strain of the structure, while the effect on the stress is relatively small. On the other hand, when the unbalanced coefficient increases, the total deformation magnitude of the structure will show a basically symmetric distribution with an unbalanced coefficient of 0.5. When the unbalanced coefficient is 0.5 (balanced lay-up), the structural deformation is minimal, and subsequently either decreasing or increasing the unbalanced coefficient will lead to an increase in the deformation of the laminate. This is in agreement with the description of the paper [28].

The relationship between the failure variation of the structure and the unbalanced coefficient is shown in Figure 7(d). As the unbalanced coefficient increases, the failure coefficient increases and then decreases. The maximum failure coefficient occurs at the unbalanced coefficient of 0.5. This indicates that the laminates with a balanced structure are more susceptible to structural failure than the laminates with an unbalanced structure. It is noteworthy that the structural strain case after lay-up shows a strong agreement with the failure case. The abrupt change points of both are concentrated at the unbalanced coefficients of 0.3 and 0.6, as well as the extreme values of both are also found at the unbalanced coefficient of 0.5. This also implies that strain may be the main cause of structural failure rather than structural stress under this lay-up.

Figure 8 gives a cloud of the failure variation after the wing skin lay-up. Under aerodynamic loading, the failure areas of both upper and lower airfoils are concentrated at the wing root, but the failure of the lower airfoil is more serious than that of the upper airfoil. Comparing the strain clouds of the wing structure, it can be seen that the strain distribution area of the lower airfoil is larger than that of the upper airfoil, and the distribution of the failure area is basically the same.

4.2. Effect of Lay-Up Angle on Structural Failure

The effects of different lay-up angles on the deformation, stress, strain, and failure of the laminate were compared, as shown in Figure 9. For the deformation of the laminate, the deformation size distribution of the laminate composed of three different lay-up angles is symmetrically distributed along the unbalanced coefficient of 0.5, and the overall elastic deformation decreases from the unbalanced coefficient of 0.2 to the unbalanced coefficient of 0.5 and then increases in the opposite direction, and the increase of the lay-up angle has a significant change on the elastic deformation of the laminate. The elastic deformation caused by the increase of ply angle will decrease. From Figure 9(b), the maximum stresses of the laminates with the three lay-up angles decreased with the increase of the unbalanced coefficient, and the basic changes showed a linear trend. It is noteworthy that when comparing the structural strains and failures, it is found that the 45° and 60° lay-up angles show consistent variation patterns, while the 30° lay-up angle shows the opposite variation to the above two lay-up angles. This can be interpreted as a change in the stiffness matrix in the composite calculations due to the inversion of the trigonometric relationship. However, the strains and failures of the three angles of the laminate remain consistent, i.e., when the strain of the structure increases, the failure factor of the structure decreases.

Interestingly, the inflection point of the change in all three lay-up angles occurs at an unbalanced coefficient of 0.4, and the progressive interval of the slope change is from 0.4 to 0.6 for the unbalanced coefficient. This also shows that for the three angles mentioned above, a large change in material properties occurs when the unbalanced plate is converted to the balanced plate and then to the unbalanced plate. Here, the laminates below the unbalanced coefficient of 0.4 are recorded as low unbalanced coefficient laminates, the laminates with unbalanced coefficients of 0.4-0.5 and 0.5-0.6 are recorded as quasi-balanced laminates, and the laminates above the unbalanced coefficient of 0.6 are recorded as high balanced coefficient laminates. According to the numerical simulation results, the mechanical properties of the three types of laminates will also remain relatively stable within their respective spanning intervals.

4.3. Structural Optimization Analysis of Unbalanced Laminates

In order to obtain the minimum elastic deformation of the wing structure under the laminate failure constraint, the optimization analysis is carried out by using the two-level ply optimization concept described in Section 2, based on the clear influence of the ply structure parameters on the structure and failure of the laminate. The initial lay-up scheme is selected as a 45° laminate with an unbalanced coefficient of 0.5, and the upper and lower lay-up schemes are the same. The total number of initial samples is defined as 1000, and to optimally filter this data. The optimization results were obtained after two-level of optimization screening as shown in Table 4.

The unbalanced laminate with mixed angles is obtained after the first optimization, where the unbalanced coefficients are 0.5 for 30° ply angle, 0.6 for 45° ply angle, and 0.8 for 60° ply angle. Compared with the initial pavement solution, the deformation of the first optimization is reduced by 12.27%, and the failure coefficient is increased by 37.06%. The deformation of the second level optimization is reduced by 13.29% and the failure coefficients are increased by 00.39% compared with the first optimization, while the deformation of the initial scheme is reduced by 23.93% and the failure coefficients are reduced by 37.04% compared with the initial scheme. Since the static aeroelasticity will lead to negative torsional of the wing, the trailing edge of the wing will produce more obvious deformation at this time. Therefore, the lift direction of the trailing edge of the wing (i.e., the Y-axis direction) was taken as the deformation output here. A comparison of the deformation cases obtained from the final optimization is shown in Figure 10.

The optimized failure case is shown in Figure 11. Here, the failure of the upper airfoil is denoted as F-U, the failure of the lower airfoil is denoted as F-L, the strain case of the upper airfoil is denoted as S-U, and the strain case of the lower airfoil is denoted as S-L. It can be seen from the clouds that the failure area at the upper and lower airfoils of the wing structure gradually decreases as the optimization scheme changes. The main change area of strain also corresponds to the main change area of failure. This is consistent with the case in the above subsection.

It should be noted that for the wing structure other than the skin, the comparison of the results without optimization and twice optimization reveals no significant change in the location of the strain and stress at the wing root structure. Only the magnitude of the stresses and strains changed due to the change in the structural stresses. Both optimizations resulted in increased stresses and strains in the metal wing ribs and beams. The stress-strain clouds of the wing beam and wing rib are shown in Figure 12.

Finally, a schematic diagram comparing the performance of the laminate before and after optimization in polar coordinates is given as shown in Figure 13. According to the figure below, the laminate has a larger but smaller range of tensile strength in the vertical and horizontal directions during the initial scheme phase. The first optimization increased the orthogonal tensile strength in quadrants one and three, and quadrants two and four, as well as the interlaminar shear strength. The second optimization further expands the range of tensile and shear strengths in all directions of the laminate based on the first optimization. Although this results in a reduction of the strength magnitude, the resistance of the laminate in all directions is enhanced.

After several iterations of optimization, the original single-direction laminate is gradually converted into an unbalanced and quasi-isotropic laminate due to the change of lay-up method and the variety of lay-up angles. This type of laminate is also widely used in the aeroelastic shear design of aircraft. These laminates have relatively good resistance in the face of bending, torsion, and shear deformation. It is also assumed that the change in a laminate also results in a change in the strain and failure regions and in the magnitude of the values taken. Compared with traditional optimization methods, the optimization method based on the idea of integrating variables discussed in this paper performs fast optimization of structural parameters (balanced/unbalanced, equal/nonequal angle) of the laminate, reducing the inconvenient optimization problems caused by a large number of parameters and improving the optimization efficiency of the wing skin lay-up to a certain extent.

5. Conclusion

The influence of the unbalanced coefficient and the change of the laminate’s lay-up angle on structural failure in the context of the static aeroelastic application of the vehicle is investigated in this paper using the computational method of fluid-structure coupling, and the structural parameters of the laminate are optimized at the end using the two-level optimization idea.

In the first part of the article, the main study of the relationship between the variation of the unbalanced coefficient of the laminate and the mechanical properties is established. Both the stress and strain curves are decreasing. The slope of the stress curve was relatively stable, while the slope of the strain curve changed three times. The failure variation is the same as the stress variation pattern, implying that strain may be the main cause of structural failure rather than structural stress.

The second part of the study focuses on the lay-up angle of the unbalanced laminate. It is found that the elastic deformation, stress-strain, and failure coefficients decrease with increasing ply angle. 45° and 60° ply angles show the same variation in structural strain and failure, while the 30° ply angle shows the opposite variation. For the three different lay-up angles, the laminates can be classified as low unbalanced coefficient laminates, quasi-equilibrium coefficient laminates and high unbalanced coefficient laminates according to the magnitude of the unbalanced coefficient, and the mechanical properties of the three types of laminates will remain relatively stable within their respective span.

In the final part, Based on the Screening method, the two-level optimization method of “balanced layer” and “equal-angle layer” is used to optimize the structural parameters of the laminate, which results in a 23.93% reduction in the elastic deformation of the optimized laminate compared with the initial one, and a 37.04% reduction in the failure coefficient. Compared with the traditional optimization method, this method reduces the inconvenience of optimization caused by a large number of parameters and improves the optimization efficiency of the wing skin laminate to a certain extent.

Data Availability

Some of the data and models used to support the results of this study can be obtained from the authors upon reasonable request. All simulation data were used only during the study and do not involve any commercial profit situation, which includes the Excel worksheet used for documentation. The experimental data in this paper include calculated stress-strain, deformation and failure data, design solutions for the lay-ups, and the wing model parameters used. Information on the computational grid and initial computational conditions are also available through the authors.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

The authors gratefully acknowledge support from the Shaanxi Provincial Department of Science and Technology and the Shaanxi Provincial Key Laboratory of Industrial Automation. This study was supported by the General Project of Shaanxi Provincial Department of Science and Technology (No.2020JM-600); Shaanxi Province technology Innovation Guidance Special Project (No.2021QFY05-03); and Shaanxi University of Technology Graduate Innovation Fund(No. SLGYCX2232).