## Abstract

Nonuniform slot air bearing (NSAB) systems have two major advantages, the external air supply and slot restrictor design, and their inherent multidirectional supporting forces and stiffness that provide excellent rotational stability. However, NSAB systems are prone to vibration from nonperiodic or chaotic motion caused by nonlinear pressure distribution within the gas film, gas supply imbalance, or simply inappropriate design. It is necessary to determine under which conditions these nonperiodic motions arise, and to design a NSAB system that will resist these vibrations. The dynamic behavior of a rotor supported by an NSAB system was studied using spectral response, bifurcation, Poincaré map, and the maximum Lyapunov exponent. The numerical results showed that chaos in an NSAB system occurred within specific ranges of rotor mass and bearing number. For example, the chaotic regions where the maximum Lyapunov exponents were greater than zero occurred in the intervals of rotor mass 20.84 ≦ m_{f} < 24.1 kg with a bearing number of *Λ* = 3.45. In addition, the coupling effect of rotor mass and bearing number was also investigated. To predict chaotic behavior, ensemble regression, and the back propagation neural network were used to forecast the occurrence of chaos. It was found that ensemble regression using dataset of 26 × 121 gave the best results and most accurate prediction for this NSAB system. The results may make a valuable contribution to the design of NSAB systems for use in a wide variety of industrial applications.

## Introduction

The mathematical theory of thin fluid lubrication was first derived by Reynolds [1] in 1886, in his famous partial differential equation related to pressure, density, relative motion, and velocity. The Reynolds equation established the basis of fluid lubrication theory. In 1961, Malanoski and Loeb [2] pioneered passive throttling when he used capillarity and small orifices to control pressure distribution in hydrostatic bearings and to improve the problem of insufficient oil film stiffness. Mayer and Shaw [3] used a variable throttling device to bring better rigidity to the bearing system and achieved the same, or even better results, than other throttling methods such as those using membranes. Later, related scholars [4] discussed the advantages and applications of this type of throttling method, until in 2004, Park and Kim [5] developed a slot-restricted gas journal bearing with nonuniform slot clearance in the circumferential direction that proved to be much more stable than a conventional slot-restricted bearing. The experimental instability threshold speed was compared to the theoretical results and confirmed the validity of numerical analysis. Slit-type air-supply bearings have now become quite common in precision equipment. In 2013, Garg and Kumar [6] found that a non-Newtonian fluid in the bearing system undergoes a viscosity change that is caused by a rise in temperature. This has a significant impact on the performance of an orifice throttling hybrid journal bearing system. Kushare and Sharma [7] studied the effect of a non-Newtonian fluid on the stability of a double-lobed hybrid bearing with symmetrical orifice throttling. In 2015, Garg [8] did a theoretical study of the stability parameters and rheological effects of lubricating fluid on various groove-type hybrid journal bearing systems. He used the Reynolds equation and the finite element method to calculate the bearing fluid clearance space and the flow relationship with the throttle. This was especially so for non-Newtonian fluids, where the stiffness coefficient, damping coefficient and rotation frequency of the bearing system are calculated in the mixed journal bearings of different configurations of throttle, and compared with the behavior of a Newtonian fluid. From the above references, it can be seen that the use of appropriate throttles in these precision bearing applications can definitely improve performance. Active throttles are more effective than passive ones. However, most of the previous studies focus on oil film bearing systems, and research into the behavior of air film bearings is relatively rare. Almost all the studies of vibration, dynamic performance, and prediction of bearing rotor behavior have been made with oil film bearings.

Experiments into rotor dynamic behavior carried out by Bently [9] in 1974 showed that oil film bearing systems exhibited second and third subharmonic vibration. Childs [10] used analytical methods to demonstrate the existence of subharmonic vibrations in rotor bearing systems. Wang and Lin [11] and Wang et al. [12] investigated the air film pressure and dynamic performance of several different air pressure bearing systems, and found that periodic, quasi-periodic, or subharmonic motion, and even chaos, arose in the rotor and journal center under different operating conditions. The effect of chaos on an air film bearing system is unpredictable and it must be avoided because it is especially dangerous to precision bearing systems. Therefore, in this study, the main objective has been the analysis of the dynamic behavior of nonuniform slot air bearing (NSAB) systems under different operating conditions. The maximum Lyapunov exponents (MLE) [13] characteristics of the NSAB system can be used to determine the presence of chaotic phenomena and accurately predict the dynamic trajectory of the system. There are many different machine learning (ML) models that can be used to solve the MLE prediction problem: multilayer perceptron (MLP) [14], ensemble regression [15], the back propagation neural network (BPNN) [16], the support vector machine (SVM) [17], the decision tree (DT) [18], and random forest (RF) [19]. All these provide excellent performance in many other research areas as well [20–22]. However, they have rarely been used to predict MLE values.

The major contributions made by this paper are: (1) The derivations of the MLE calculation process; (2) the MLE prediction model to reduce computation time; (3) validation of the feasibility of the proposed framework; and (4) comparison of the performance of the machine models.

The design and analysis of the NSAB system are introduced in The Design and Analysis of a Nonuniform Slot Air Bearing System section. The system architecture and derivations of the MLE calculation process are detailed in the Rotor Dynamic and Numerical Simulation section. Results and Discussion section presents the experimental results and suggestions and the Conclusion section presents the conclusions.

## The Design and Analysis of a Nonuniform Slot Air Bearing System

Pneumatic bearings often have the problem of insufficient support force and the need for high-speed stability and for these reasons the application of slit air-supply throttle bearings has received considerable attention. In this study, a new type of nonuniform slit air supply throttling shaft system has been proposed which is easy to manufacture. The new shaft system has a hybrid bearing that supports both axial and radial loads. Coupled with the adjustable slit width air supply method, this type of bearing has higher air film stiffness and stability. The structure is shown in Fig. 1.

The new type of nonuniform slot air bearing has slots of two different widths in the circumferential direction of the bearing, see Fig. 1. As shown in the figure, the width of the second zone (*a*_{2}) is larger than that of the first zone (*a*_{1}). The range of the first zone depends on the nonuniform angle ($\alpha a$). When $\alpha a$ is 0 deg or 180 deg, this new type of bearing is the equivalent of a radial shaft bearing system. The main function of the nonuniform slot is adjustment of the pressure of the external air supply. In addition to the avoidance of air hammer, the second area also acts as an air supply cavity that stabilizes air flow.

The relationship between the pressure of the internal fluid and the thickness of the gas film when the shaft system is supporting the rotor, must be based on the Reynolds equation, and the governing equations in the bearing gap and the slot area are different, as shown in Eqs. (1) and (2) below:

where $P\u0303n$, $h\u0303n$, $Bs$, *Λ*, and *τ* are dimensionless air pressure, dimensionless air film thickness, slot width ratio, bearing number, and dimensionless time. $\theta ,\u2009\u2009\zeta $ are the angular direction and the dimensionless axial coordinate system *R _{n}*

*=*

*r/r*is the dimensionless radial coordinate system.

_{o}*Q*is the flowrate in the slot,

_{s}*B*is the bearing parameter,

_{p}*R*

_{s}*=*

*r*is the air supply pressure, and

_{s}/r_{o}, P_{s}*Q*must satisfy the Eq. (7)

_{s}where $\kappa $ is the isotropic coefficient.

The pressure distribution function, bearing capacity, mass flow, stiffness, and other performance coefficients of the bearing system can be obtained by calculating the Reynolds equation through Eqs. (1)–(7). At the same time, pressure changes inside the bearing will be used to solve the dynamic performance and rotor motion behavior.

## Rotor Dynamic and Numerical Simulation

### Rotor Dynamics Analysis.

The analysis presented in this study considers a rotor-bearing system comprised of a perfectly balanced, flexible rotor of mass *m _{f}* supported symmetrically by two identical NSABs which are mounted on rigid pedestals. Given these idealizations, the current analysis can be confined to a single bearing supporting a rotating rotor of mass mf with 2 deg of translatory oscillation in the transverse plane.

The iterative computation procedure can be summarized and shown in Fig. 2.

### Numerical Simulation.

*τ*, and hence Eqs. (1) and (2) become

and “$\u2297$” denotes the convolution operation in the *K* domain. The finite difference method was then used to discretize Eqs. (11) and (12) with respect to the *θ* and $\zeta $ directions. Note that Eqs. (11) and (12) are discretized using the second-order-accurate central-difference scheme for both the first and the second derivatives.

From Eq. (16), $P\u0303n(l)i,j,k$ was obtained for each time interval, where *i*, *j*, and *k* indicate the node position and *l* indicates the *l*th term. Note that in this study, the time-series data of the first 1000 revolutions were excluded from the dynamic behavior investigation to ensure that the analyzed data corresponded to the steady-state conditions. The data included the orbital paths and velocity of the rotor center *O*_{2} and the orbits of the journal center *O*_{3}.

## Results and Discussion

### Verification Analysis for Different Numerical Methods.

Two different numerical methods, perturbation and the hybrid method, combined with finite difference and differential transformation, were used to solve the rotor displacement in the new nonuniform slot air bearing system. The results showed that the simulation values obtained by the hybrid and perturbation methods were consistent. This proved that the calculation of the Reynolds equation using the hybrid method was correct and this verified the applicability of the hybrid method at the same time.

The accuracy of these two numerical methods were further compared, and the results showed that displacements of the rotor center in the horizontal and vertical directions, calculated by the perturbation method, were prone to cause unstable variation in numerical value. The accuracy of the perturbation method can reach three decimal digits, as shown in Table 1. It was also observed that the precision of displacement value calculated by the hybrid method, regardless of the calculation time interval, can reach five decimal digits. Therefore, the hybrid method was more accurate than the perturbation method for the NSAB system.

X_{2}(nT) | Y_{2}(nT) | ||||
---|---|---|---|---|---|

Numerical method | Operational parameter | $H\u0303$ = 0.001 | $H\u0303$ = 0.01 | $H\u0303$ = 0.001 | $H\u0303$ = 0.01 |

Perturbation method | m = 6.45 kg_{f} | 0.071055 | 0.071626 | 0.086724 | 0.086736 |

Hybrid method | Λ = 3.1 | 0.071164 | 0.071165 | 0.087715 | 0.087718 |

Perturbation method | m = 7.98 kg_{f} | 0.082046 | 0.082351 | 0.099710 | 0.099870 |

Hybrid method | Λ = 3.1 | 0.088910 | 0.889119 | 0.098734 | 0.098733 |

Perturbation method | m = 9.8 kg_{f} | 0.099159 | 0.099393 | 0.178628 | 0.178891 |

Hybrid method | Λ = 5.47 | 0.099840 | 0.099840 | 0.178655 | 0.178652 |

Perturbation method | m = 9.8 kg_{f} | 0.099200 | 0.099239 | 0.199745 | 0.199420 |

Hybrid method | Λ = 5.5 | 0.099604 | 0.099603 | 0.199746 | 0.199745 |

X_{2}(nT) | Y_{2}(nT) | ||||
---|---|---|---|---|---|

Numerical method | Operational parameter | $H\u0303$ = 0.001 | $H\u0303$ = 0.01 | $H\u0303$ = 0.001 | $H\u0303$ = 0.01 |

Perturbation method | m = 6.45 kg_{f} | 0.071055 | 0.071626 | 0.086724 | 0.086736 |

Hybrid method | Λ = 3.1 | 0.071164 | 0.071165 | 0.087715 | 0.087718 |

Perturbation method | m = 7.98 kg_{f} | 0.082046 | 0.082351 | 0.099710 | 0.099870 |

Hybrid method | Λ = 3.1 | 0.088910 | 0.889119 | 0.098734 | 0.098733 |

Perturbation method | m = 9.8 kg_{f} | 0.099159 | 0.099393 | 0.178628 | 0.178891 |

Hybrid method | Λ = 5.47 | 0.099840 | 0.099840 | 0.178655 | 0.178652 |

Perturbation method | m = 9.8 kg_{f} | 0.099200 | 0.099239 | 0.199745 | 0.199420 |

Hybrid method | Λ = 5.5 | 0.099604 | 0.099603 | 0.199746 | 0.199745 |

The results show that the hybrid method used in this study had excellent convergence and applicability for the NSAB system. An analysis of the stability of the numerical value at time intervals of 10^{−2} and 10^{−3} was also carried out. Changes in displacement under different operating conditions were compared with the numerical value of the Poincaré cross section diagram (PCSD). The test results in Table 1 show that the numerical displacements for different operational parameters were sufficiently precise to at least 10^{−2} or 10^{−3}. This meant a time interval value of 10^{−2} would be accurate enough. Therefore, 0.01 was used as the calculation time interval to ensure efficiency of the subsequent calculation of the bifurcation characteristics and the dynamics of the NSAB system.

### Analysis of the Dynamic Behavior of the NSAB System Using Rotor Mass as the Bifurcation Parameter.

The circumferential slot in the NSAB has two different widths with ratios *B _{s}*

_{1}=

*a*

_{1}

*/C*and

_{r}*B*

_{s}_{2}=

*a*

_{2}

*/C*which were set at 0.4 mm and 0.8 mm, and the nonuniform angle ($\alpha a$) was set to 105 deg, see Fig. 1. The mass of the rotor affects the strength of air floatation and the overall stability of the system. In this section, the influence of the rotor mass on the NSAB system has been analyzed. The bearing number of the system was assumed to be

_{r}*Λ*= 3.45, and the rotor mass was used as the bifurcation parameter:

#### Rotor Dynamic Behavior.

Figures 3.1(a) and 3.1(b)–3.3(a) and 3.3(b) show the phase plane of rotor center (*X*_{2}, *Y*_{2}) and the orbit of journal center (*X*_{3}, *Y*_{3}). Figures 3.1(c) and 3.1(d)–3.3(c) and 3.3(d) reveal the dynamic response of the rotor center in the horizontal and vertical directions. From the results, the rotor exhibited a regular and periodic pattern at low rotor mass (*m _{f}* = 8.4 kg). However, when the mass increased to 15.8, nonperiodic phenomena replaced regular motion. When the mass increased further, to 18.05 kg, 2T subharmonic motion appeared.

#### Chaotic Analysis and Verification.

The rotor mass range was set to between 0.1 and 25.0 kg for actual operating conditions, as shown in the bifurcation diagrams in Fig. 4. It can be seen that when *m _{f}* < 15.8 kg, the center of the system rotor exhibited T periodic motion in both the horizontal and vertical directions. This can also be proven in Fig. 5(a) of PCSD. However, when the mass increased to 15.8 kg, rotor center motion became quasi-periodic. It can be found that this motion was maintained over the range of 15.8 ≦

*m*< 18.05 kg. As the mass rose to 18.05 kg, the system became 2T-subharmonic motion and remained so within the range of 18.05 ≦

_{f}*m*< 18.38 kg. When the mass reached 18.38 kg, the 2T-subharmonic motion bifurcated into quasi-periodic motion again. This motion persisted when 18.38 ≦

_{f}*m*<19.08 kg, but became 2T-subharmonic motion when the mass increased to 19.08 kg. This motion continued for 19.08 ≦

_{f}*m*< 20.84 kg. When the rotor mass increased to 20.84 kg, chaotic motion appeared and this motion persisted to 20.84 ≦

_{f}*m*<24.1 kg. When the rotor mass rose to 24.1 kg, chaotic motion was replaced by stable periodic motion and this motion persisted to 24.1 ≦

_{f}*m*≦ 25.0 kg.

_{f}The maximum-Lyapunov-exponents (MLE) were used to prove the occurrence of chaotic motion, and for interpretation and further verification of the study findings. The MLE tended to zero or less than zero, the system did not show chaotic behavior as shown in Fig. 5(b).

It is clear from these findings that movement of the rotor presents quite a complicated pattern. Table 2 shows a general summary of the movement patterns of the system as they varied with rotor mass.

### Analysis of the Dynamic Behavior of the NSAB System Using Bearing Number as the Bifurcation Parameter.

In an NSAB system, the rotor bearing number (speed) directly affects pressure distribution in the bearing and also influences the relative stability and performance of the entire system. Therefore, in this section, the bearing number (*Λ*) has been used as the bifurcation parameter, and the rotor mass was set to *m _{f}* = 5.6 kg. The relevant dynamic behavior of the bearing system was analyzed and discussed:

#### Chaotic Analysis and Verification.

The bearing number *Λ* was used as the main analysis parameter to explore the influence of different bearing numbers on the NSAB system and it was set to a range between 1.0 and 6.0 for actual operating conditions. When the bearing number rose to 4.48 ≦ *Λ* < 5.47, the system generated 2T and chaotic motion alternately. However, as it increased to 5.47 ≦ *Λ* ≦6.0, the chaotic state began to fade and the system tended to exhibit stable T periodic motion. All motion behavior generated in the interval of 1.0 ≦ *Λ* ≦ 6.0 was summarized in Table 3.

Λ | [1.0, 3.58) | [3.58, 4.47) | [4.47, 4.48) | [4.48, 4.49) | [4.49, 4.5) |

Behavior | T | 2T | Chaos | 2T | Chaos |

Λ | [4.5, 4.58) | [4.58, 4.59) | [4.59, 4.61) | [4.61, 4.62) | [4.62, 4.67) |

Behavior | 2T | Chaos | 2T | Chaos | 2T |

Λ | [4.67, 4.68) | [4.68, 5.21) | [5.21, 5.47) | [5.47, 6.0] | |

Behavior | Chaos | 2T | Chaos | T |

Λ | [1.0, 3.58) | [3.58, 4.47) | [4.47, 4.48) | [4.48, 4.49) | [4.49, 4.5) |

Behavior | T | 2T | Chaos | 2T | Chaos |

Λ | [4.5, 4.58) | [4.58, 4.59) | [4.59, 4.61) | [4.61, 4.62) | [4.62, 4.67) |

Behavior | 2T | Chaos | 2T | Chaos | 2T |

Λ | [4.67, 4.68) | [4.68, 5.21) | [5.21, 5.47) | [5.47, 6.0] | |

Behavior | Chaos | 2T | Chaos | T |

To confirm chaotic behavior, the MLE was calculated with the bearing number as shown in Fig. 6. The results showed that chaos and periodic motion alternated (occurred in the interval of 4.67 ≦ Λ < 4.68), and this was consistent with the findings in Table 3. The area of chaotic motion, where the MLE was greater than zero (red curve), occurred in the interval of 5.21 ≦ *Λ* < 5.47, the remaining blue curves are in the stable area and are fully consistent with the results shown in Table 3.

### Verification and Prediction of the Occurrence of Chaotic Behavior in an NSAB System.

#### Establishment of Stable and Unstable Areas in an NSAB System.

To gain an understanding of the behavior of the bearing system under different parameters, a range of different rotor masses and bearing numbers were used in an analysis of their interactive influence on the bearing system. The MLE was used to determine regions of stability and instability. In Fig. 7, yellow to deep red in the figure mark areas of instability and chaos. The blue to deep blue areas are regions of stability. When the bearing number was high, system motion was nonperiodic and sometimes chaotic with different rotor masses. This kind of three-dimensional graphical representation is highly effective for distinguishing stable and unstable regions of operation for these air film bearing systems and can be used as an important reference for bearing system design.

#### Prediction of Chaotic Motion.

Ensemble regression and the back propagation neural network (BPNN) were used to analyze the MLE of the bearing system to predict the occurrence of chaos and prevent the NSAB system from becoming unstable.

Figure 8 shows the ensemble regression analysis results using three different training datasets to predict the MLEs of the NSAB system. A comparison with Fig. 7, showed that chaotic behavior in the range of 20.0 ≦ *m _{f}* <25.0 kg and 5.0 ≦ $\Lambda \u2009$< 6.0 was not obvious, but the chaos that occurred at 4.0 ≦

*m*<5.0 kg and 1.0 ≦ $\Lambda \u2009$< 2.5 could be accurately predicted.

_{f}Figure 9 shows the results of BPNN analysis using the three different training datasets to predict the MLEs of the NSAB system. In Fig. 9(a), it can be seen that some part of the chaotic behavior could not be forecast with precision using the 26 × 61 dataset when 4.0 ≦*m _{f}* < 5.0 kg and 1.0 ≦ $\u2009\Lambda \u2009$< 2.5. In Figs. 9(b) and 9(c), for chaos in the ranges of 20.0 ≦

*m*< 25.0 kg and 5.0 ≦ $\Lambda \u2009$< 6.0, the 26 × 121 dataset predicted a larger chaotic area than the 26 × 81 dataset.

_{f}A comparison was made between the ensemble regression and BPNN with the three different training datasets and the values of root-mean-square error (RMSE), mean-square error (MSE), mean absolute error (MAE), and coefficient of determination (*R*^{2}) were applied to evaluate the performance of prediction for chaotic motions shown in Table 4. *R*^{2} value is between 0 and 1, and the closer *R*^{2} is to 1, the closer the predicted value is to the true value.

Method | Data | RMSE | R^{2} | MSE | MAE |
---|---|---|---|---|---|

Ensemble regression | 26 × 61 | 0.9267 | 0.7279 | 0.8587 | 0.5947 |

26 × 81 | 1.0093 | 0.6772 | 1.0187 | 0.6223 | |

26 × 121 | 0.8722 | 0.7589 | 0.7607 | 0.5454 | |

BPNN | 26 × 61 | 1.1465 | 0.5835 | 1.3145 | 0.7273 |

26 × 81 | 1.0034 | 0.6810 | 1.0068 | 0.6374 | |

26 × 121 | 1.4390 | 0.7538 | 2.0708 | 0.8487 |

Method | Data | RMSE | R^{2} | MSE | MAE |
---|---|---|---|---|---|

Ensemble regression | 26 × 61 | 0.9267 | 0.7279 | 0.8587 | 0.5947 |

26 × 81 | 1.0093 | 0.6772 | 1.0187 | 0.6223 | |

26 × 121 | 0.8722 | 0.7589 | 0.7607 | 0.5454 | |

BPNN | 26 × 61 | 1.1465 | 0.5835 | 1.3145 | 0.7273 |

26 × 81 | 1.0034 | 0.6810 | 1.0068 | 0.6374 | |

26 × 121 | 1.4390 | 0.7538 | 2.0708 | 0.8487 |

For using ensemble regression, the prediction result by 26 × 81 dataset is the worst compared with 26 × 61 and 26 × 121 datasets by evaluating RMSE, MSE, and MAE. The values of RMSE, MSE, and MAE for 26 × 81 dataset are all larger than those of the other two datasets shown in Table 4 and it means that the accuracy of prediction by 26 × 81 dataset is the worst. Meanwhile, the value of *R*^{2} by 26 × 81 dataset is lower than the other two datasets and it also means the forecast by 26 × 81 dataset is poor.

It was also found that the ensemble regression prediction made using the 26 × 61 dataset was better than that using BPNN. Meanwhile, the value of *R*^{2} by ensemble regression is closer to 1 than BPNN and verifies that the ensemble regression is more suitable than BPNN for the forecast of chaos in an NSAB system.

In Table 4, it is worth noting that the 26 × 121 dataset of ensemble regression presented the smallest errors including RMSE, MSE, and MAE than any other dataset by two machine learning methods. The value of *R*^{2} by 26 × 121 dataset of ensemble regression was 0.7589 and it was the highest value in all results. Accordingly, with respect to RMSE, MSE, MAE, and *R*^{2} detections, ensemble regression has the advantage over BPNN for the prediction of motion conditions. It is suggested that ensemble regression can be applied with dataset of 26 × 121 for the forecast of chaotic motion and obtained the best results and most accurate predictions for an NSAB system.

## Conclusion

The objective of this study was an analysis and prediction of the dynamic behavior of a flexible rotor supported by a nonuniform slot air bearing (NSAB) system. The finite difference, perturbation, and a hybrid method were used to solve the pressure distribution at the highest nonlinearity in the system. After this the dynamic equations of the flexible rotor center were used to obtain the orbits of the rotor center. Analysis was then conducted on the orbit data to generate spectrum diagrams, Poincaré maps, bifurcation diagrams, and maximum Lyapunov exponents. These showed that the orbit of the rotor center changed along with the rotor mass and the bearing number and complicated motion was synchronously generated in the horizontal and vertical directions. These included normal and subharmonic vibration, as well as periodic, quasi-periodic, and chaotic motion. The hybrid method applied in this paper yielded better accuracy and precision for the verification of values than the finite difference and perturbation methods. The chaotic regions where maximum Lyapunov exponents are greater than zero occur in the intervals of rotor mass 20.84 ≦ *m _{f}* < 24.1 kg with the bearing number

*Λ*= 3.45 and in the six intervals of bearing number with a rotor mass of

*m*= 5.6 kg. Furthermore, ensemble regression and BPNN were used to forecast the occurrence of chaos. It was found that using ensemble regression with the suggested 26 × 121 dataset gave the best results and most accurate prediction for an NSAB system. The results obtained in this study can be used as a basis for future NSAB system design and the prevention of instability.

_{f}## Funding Data

Ministry of Science and Technology (MOST), Taiwan (Award Nos. MOST-110-2221-E-167-019 and MOST-110-2622-E-167-006: Funder ID: 10.13039/501100004663).