Abstract

In this paper, the effect of a fractional constitutive model on the rheological properties of fluids and its application in numerical simulation are investigated, which is important to characterize the rheological properties of fluids and physical characteristics of materials more accurately. Based on this consideration, numerical simulation and analytical study of unsteady fractional Oldroyd-B viscoelastic flow are carried out. In order to improve the degree of accuracy, the mixed partial derivative including the fractional derivative in the differential equation is converted effectively by integrating by parts instead of by direct discretization. Then, the stability, convergence, and unique solvability of the difference scheme are verified. The validity of the finite difference method is tested by making comparisons with analytical solutions for two kinds of fractional Oldroyd-B viscoelastic flow. Numerical results obtained using the finite difference method are in good agreement with analytical solutions obtained via the variable separation method. Viscoelastic characteristics of the unsteady Poiseuille flow are similar to the second-order fluid or integer-order Oldroyd-B fluid when the parameter is close to 0 or to 1. Oscillation characteristics of fractional viscoelastic oscillatory flow are similar to those of the classical viscoelastic fluid under the same condition. Compared with the previous research, the present work studies the rheological properties of fluids with the finite difference method, and the application of fractional constitutive models in describing the rheological properties of fluids is developed. Meanwhile, more cases reflecting the fractional-order characteristics are given. The present method can accurately capture the flow characteristics of unsteady fractional Oldroyd-B viscoelastic fluid and is applicable for the generalized fractional fluid.

1. Introduction

Recently, non-Newtonian fluid models including fractional derivatives have aroused great interest because of their extensive application in practical engineering problems. Dassios and Baleanu [1] studied singular linear systems of finite difference equations by using the Caputo fractional derivative and its two alternative forms. Alsharif and Elmaboud [2] investigated the physical characteristics of electroosmotic flow by solving the fractional Cattaneo model. Khan et al. [3, 4] obtained some analytical solutions for the fractional Oldroyd-B fluid with different initial and boundary conditions. Based on the Fourier sine transform, Nadeem [5] presented exact solutions for the fractional Oldroyd-B fluid in periodic motion. Qi and Xu [6] obtained the analytical solution for the Oldroyd-B viscoelastic flow by using the discrete Laplace transform. However, it is very difficult to carry them out in the field of engineering applications. On the one hand, analytical solutions always appear in complex series form. On the other hand, the Mittag-Leffler function or the Fox-H function is often included in the exact solution. To overcome the above problems, researchers studied the Oldroyd-B viscoelastic flow with other methods. In combination with the numerical inversion of Laplace transforms, Huang et al. [7] analyzed the motion characteristics of the generalized second fluid in a double barrel rheometer. Yao and Zhang [8] set up a kind of implicit difference scheme for the viscoelastic fluid following a generalized Maxwell fractional differential equation. Wang et al. [9] studied the unsteady Poiseuille flow of the fractional Oldroyd-B viscoelastic fluid between two parallel plates and used the Stehfest algorithm to carry out numerical inversion of the Laplace transform. Liu and Zheng [10] used the Laplace transform coupling with the Hankel transform to study unsteady helical flows and obtained exact solutions of a generalized Oldroyd-B fluid between two infinite concentric cylinders. Yu [11] adopted a high-order compact finite difference method to solve the generalized fractional Oldroyd-B fluid model and mainly investigated the corresponding stability and convergence for the finite difference method. However, it should be pointed out that the related work on the numerical analysis is much less compared with that on viscoelastic rheological behaviors and analytical solutions. With the help of a discretization of the Caputo time-fractional derivative and the finite difference method, Zhang et al. [12] obtained an approximate implicit difference method to solve a two-dimensional multiterm time-fractional Oldroyd-B equation. Al-Maskari and Karaa [13] considered the numerical approximation of a generalized fractional Oldroyd-B fluid problem with a semidiscrete scheme based on the piecewise linear Galerkin finite element method. In other engineering fields, Raslan et al. [14] present a computational scheme using the bi-finite difference method (Bi-FDM) to solve the hyperbolic telegraph equation in two dimensions. The fluid flow model governed by partial differential equations is reduced to a system of nonlinear ordinary differential equations for solving the calculation in reference [15]. However, no relevant literature has been seen in which the method is applied to fractional-order models. It shows that the finite difference method has great advantage in solving the fractional Oldroyd-B fluid problem.

Large quantities of experiments and engineering practices have shown that the fractional derivative model can describe the constitutive relation of viscoelastic materials and their rheological properties more accurately than other models. Moreover, traditional constitutive models can be used as special cases of fractional derivative constitutive models. Therefore, the fractional derivative is widely used in the expanding research of the constitutive model. However, the research of applying the fractional-order derivative to the Oldroyd-B model is still relatively rare, especially when the analytical and numerical solutions are given in specific arithmetic cases and compared. The fractional-derivative-type viscoelasticity principal structure relationship established by the theory on fractional derivatives was a newly developed principal structure model nearly two decades ago. Compared to the Riemann-Liouville [16] fractional derivative, Caputo fractional derivatives [14] are more suitable and convenient for numerical analysis. Khan and Rasheed [17] have adopted the Galerkin finite element method to determine the approximate solution that is uniform with the finite difference approximation of the Caputo fractional time derivative. In another article by them, the Caputo fractional derivative is used to simulate the physical properties of the fractional unsteady magnetohydrodynamic flow [18]. Zaid and Dumitru [19] studied and proposed a new system of fractional differential equations in conjunction with the generalized Caputo fractional derivative. Erturk et al. [20] developed a mathematical analysis for a corneal-shaped model in combination with the Caputo fractional derivative. These show that the Caputo fractional derivative is widely used in mathematical and physical problems Moreover, the initial and boundary conditions defined by the Caputo fractional derivatives are specific in physical meaning, which can provide great convenience for computation and application in the field of physics and engineering. Therefore, this paper adopts Caputo fractional derivatives directly.

In summary, it has been known that the fractional derivative is widely used in the expanding research of the constitutive model to solve mathematical and physical problems, but the research of applying the fractional-order derivative to the Oldroyd-B model is still relatively rare, especially when the analytical and numerical solutions are given in specific arithmetic cases and compared. This research idea is executed in this paper, and four sections are mainly included here. The finite method is given in Section 2, including different schemes, the proof of the stability, and convergence, as well as the analysis of unique solvability. Taking the unsteady Poiseuille flow and oscillatory flow in a circular pipe as example, numerical simulations are developed in Section 3. Then, the final conclusion is given in Section 4. Moreover, the proof of the numerical simulation method and analytical solution is given in the appendix.

2. Finite Difference Method

2.1. Governing Equations with Caputo Derivative

The fundamental equations describing the unsteady motion of an incompressible fluid are given [21]: in which is the velocity vector and , , and that are, respectively, along the , , and directions in the Cartesian coordinate system. is the uniform density, and denotes the material derivative with respect to time. is the Cauchy stress tensor, in which represents the pressure, is the unit tensor, and is the extra stress tensor.

The constitutive equation of the fractional Oldroyd-B fluid [3] is where is the viscosity and and are the relaxation time and retardation time, respectively. and are the order of the fractional derivatives with the condition . is the first Rivlin-Ericksen tensor given by

where and are defined as

Here, the fractional differential operator based on Caputo’s definition [8] is expressed as

Here, is the Gamma function.

Besides, the Laplace transform of the time-fractional Caputo derivative is given as and the time-fractional Caputo derivative is defined as

The definition is stricter than that of the Riemann-Liouville fractional-order derivative, in which is an integrable function of order .

For the unidirectional flow, the velocity and stress are considered as follows:

where is the unit vector along the direction. Equations (5) and (9) satisfy the initial condition ; thus, is established. Corresponding to Equation (1), the following expressions are obtained:

In the absence of body forces, the momentum Equation (2) for the unidirectional flow of the fractional Oldroyd-B fluid is written as

According to Equations (12) and (14), the motion equation of the fractional Oldroyd-B fluid is expressed as

where .

2.2. Difference Scheme

In the space region , let us set , , is the space step, and is the time step. Here, represents the approximate value of , and it is assumed that is a -dimensional vector.

For convenience, the following marks are introduced:

The Crank-Nicolson (C-N) scheme is used to discretize partial derivative terms in Equation (11). At node , it is satisfied as . The time partial and space derivative with two orders are, respectively, expressed as

In combination with Equation (9) and the above C-N scheme, the implicit difference scheme of Equation (13) at node is expressed as

The detailed derivation process is given in Appendix A.

In addition, the initial and boundary conditions are also discretized according to the above steps.

2.3. Stability

Lemma 1. (see [22]).
The function is defined in ; if is satisfied, then is established.

Theorem 2. The implicit difference scheme of Equation (18) is unconditionally stable.

Proof. Supposing is the exact solution of Equation (18). Multiplying Equation (18) by , and summing up for from 1 to ; then, for from 1 to , it can be obtained as The following emphasis is focused on analyzing Equation (19). In combination with Equations (14) and (15), it can be concluded that The detailed proof of Equation (20) is given in Appendix B. According to Lemma 1 [22], we can obtain Therefore, the implicit difference scheme is unconditionally stable.

2.4. Convergence

Theorem 3. Equation (18) is convergent when it is satisfied by .

Proof. Suppose that is the exact solution of Equation (13) at node ; is the numerical solution, and the error between two solutions is defined as . Here, assume that there exists no error for the initial data, according to Equations (17) and (18); then, Here, . The initial and boundary values of the error are It follows from Equation (22) that We have Because is established, it can be further obtained that where is a constant which has nothing to do with and . Therefore, according to , when it is satisfied that and , then is established, and the convergence is proven.

2.5. Unique Solvability

Theorem 4. Equation (18) is uniquely solvable.

Proof. Equation (16) can be written as Here, it is satisfied that .
Equation (27) contains equations, and the part of the coefficient can be written as a matrix form: is the velocity vector and is the right-hand side (R.H.S.) of Equation (27), where Because the coefficient matrix is strictly diagonally dominant, the difference scheme Equation (18) is uniquely solvable.

3. Numerical Study of Poiseuille Flow with Finite Difference Method

Figure 1 shows the two-dimensional Poiseuille flow between two parallel plates. In order to verify the finite difference method, the analytical solution of the Poiseuille flow is firstly given: where

The above analytical solution satisfies the following initial and boundary conditions:

Because the form of the analytical solution is very complex, we only consider the approximate analytical solution under the condition that is a small number. In order to execute the numerical simulation, it is supposed that and . In Figure 2, the velocity on the centerline of the two parallel plates is considered, and numerical comparison between the finite difference method and the approximate analytical solution is shown. Here, three groups of parameters are chosen:

In the first two group parameters, the variation induced by and is compared. In the first and third group parameters, the variation induced by and is considered. As displayed in Figure 2, numerical results illustrate that the finite difference method is in good agreement with the approximate analytical solution under three group parameters. It can be preliminarily determined that the velocity amplitude shows a decreasing trend with the above four parameters. The velocity based on the parameter of group 1 is larger than other two group parameters. Meanwhile, Table 1 displays the relative error of the velocity on the centerline with time for three group coefficients. For example, Error 1 corresponds to group 1. It demonstrates that the relative error gradually turns to be within the scope of 2% as time increases, especially for group 2 and group 3. It powerfully states that the finite difference method is stable and reliable. The relative error of the velocity on the centerline with time for the above three group coefficients is given in Table 1, and the small relative error further reflects the accuracy of the numerical algorithm.

Next, the influence of four parameters on the velocity and fluid physical properties is considered with the finite difference method. Figure 3 displays the influence of the parameter on the velocity at the center line with the other three fixed parameters; here, it is supposed that and . When the parameter or , the velocity represents the moving velocity of the second-order fluid or integer-order Oldroyd-B fluid. It can be observed that the velocity distribution with and is in good agreement with that based on integral-order constitutive Equation (18). When it satisfies , the velocity of the fractional Oldroyd-B fluid is between that of the second-order fluid and that of the integer-order Oldroyd-B fluid. It also shows in Figure 3(a) that the velocity gradually tends to be stable with time, and maximum values of the velocity increase with . This phenomenon is called velocity overshoot, which appears as the increasing elasticity and a fluctuation in velocity. As displayed in Figure 3(b), the distribution of flow is filed symmetric along the direction with different . Meanwhile, the closer the parameter is to 1, the more obvious is the velocity overshoot. The physical characteristics of the fluid is manifested as an increase in fluid elasticity and appearances of velocity fluctuations, and the flow field variation is also closer to that of the integer-order Oldroyd-B fluid. In summary, the fractional Oldroyd-B fluid displays the viscoelastic characteristics similar to the second-order fluid or the integer-order Oldroyd-B fluid when the parameter is close to 0 or to 1.

Figure 4 mainly studies the influence of the order parameter on the velocity at the center line with fixed parameters; here, it is satisfied as and . As shown in Figure 4(a), the velocity increases gradually at the steady state, but the variation of the amplitude reduces with . It illustrates that the fluid viscosity reduces with ; this caused the amplitude of the velocity to reduce gradually. As shown in Figure 4(b), when is large, the variation of the velocity is basically the same. The parameter almost has no influence on the amplitude of velocity and the time to steady state. There appears a velocity overshoot phenomenon which is identical with Figure 3. It can be concluded that variations of parameters and can lead to velocity overshoot; this is because an increase in causes an increase in the elasticity of the fluid, while an increase in causes a decrease in the viscosity of the fluid.

4. Numerical Simulations of the Oscillatory Flow in a Circular Pipe

The oscillatory flow in a circular pipe exists extensively in the blood, petroleum, and polymer solutions. This kind of flow not only has the viscosity of the fluid but also shows the elasticity of the solid. Therefore, it is difficult to use one model to represent all the characteristics of viscoelastic fluids, and many viscoelastic models and constitutive equations have emerged. Fractional Oldroyd-B has made a breakthrough in the research of material constitutive theory. In general, only a simple analysis of the velocity and stress is performed in the literature, but the variation of periodic velocity and stress amplitude is rarely studied. In this paper, the velocity and shear stress of the oscillatory flow in a circular pipe is studied in detail.

The laminar flow in a circular pipe with an infinite radius is considered, and it is supposed that (1) the flow is incompressible, (2) the velocity along the radius direction is equal to zero, (3) the flow is axisymmetric about the pipe, and (4) the axial velocity is only related with the radius.

In order to verify the correctness of the finite difference method, the analytical solution of the oscillatory flow is firstly derived in detail, and then, numerical simulations of the velocity and shear stress are executed and compared between two kind of methods. The cylindrical coordinate system with is established here, in which the axis is along the central axis of the circular pipe. Here, the processing of the difference format for shear stress is similar to that for the velocity.

Starting from Equations (10) and (11), the variable separation method is used to derive the analytical solution. Here, the dimensionless analytical amplitudes of the velocity and shear stress on the pipe wall are given:

The detailed process is contained in Appendix C.

According to the expressions of and , the amplitudes of velocity and stress are related with the pipe radius , oscillation frequency , and specific value , as well as the order of the fractional Oldroyd-B constitutive model. In the numerical simulation, the influence of these parameters is studied.

4.1. The Influence of Parameters on the Velocity

Next, numerical simulations of the oscillatory flow in a circular pipe are executed. Firstly, the finite difference method is verified by comparing numerical results and the analytical solution. When the order parameters are chosen as , the fractional Oldroyd-B constitutive model degenerates to the classical Oldroyd-B constitutive model. Figure 5 shows the variation of the velocity amplitude with the oscillation frequency. It demonstrates that numerical results are in good agreement with the analytical solution for different pipe radii. The velocity amplitude basically shows a decreasing trend with the pipe radius. In Figure 5(a), the resonance phenomenon appears at , because the peak value of the velocity amplitude is much larger than at other frequencies. In Figures 5(b)5(d), the number of the formant is gradually increasing with the pipe radius. Meanwhile, the peak value substantially reduces especially for the velocity amplitude at . Compared with the case of , all formants concentrate at . The curve of the velocity amplitude gradually approximates to a monotonic decreasing curve with , which is very close to the curve of the Newtonian fluid.

Secondly, a different order parameter is chosen to study the variation of the velocity amplitude. Figure 6 shows the velocity amplitude with and , and Figure 7 chooses the order parameter as and . It can be observed that the formant also appears in the velocity amplitude curve of the fractional Oldroyd-B fluid similar to the classical Oldroyd-B fluid. On the one hand, the first peak appears in Figure 7(a) corresponding to , and the number of the formant increases, but the peak value substantially reduces with . On the other hand, the location of the peak value moves left along the horizontal axis. Otherwise, when the parameter is chosen as and in Figure 7, the velocity amplitude also decreases quickly comparing with Figures 6(a) and 6(b). In brief, the velocity amplitude reduces gradually with and the parameters and . Comparison with Figure 5 shows that the physical variation laws of fractional-order Oldroyd-B fluids are similar to those of classical Oldroyd-B fluids. It further shows that the oscillation peak and oscillation frequency can be reduced by reasonably controlling the above parameters.

In addition, we discuss the influence of and on the velocity amplitude. In Figure 8, the influence of with fixed and the influence of with fixed are all shown. It can be found that there are multiple formants in the curve with , but only one formant exists with and in Figure 8(a). However, there only exists one formant with different in Figure 8(b). The same characteristic is that the first peak with different and corresponds to the same frequency. Meanwhile, the peaks are gradually decreasing with and . On the whole, the velocity amplitude increases with the increase of and decreases with at low frequencies.

4.2. The Influence of Parameters on the Shear Stress

Zhu et al. [23] studied physical characteristics of the Maxwell flow oscillating in a pipe. Their research shows that the resonance phenomenon occurred on the amplitude of the wall shear stress at certain frequencies, and this discovery is helpful to improve the core oil displacement efficiency in the practical engineering field. Therefore, the variation of the shear stress with different parameters and frequency is further studied in this part.

At first, the parameters are chosen as and ; the influence of the frequency on the shear stress is given in Figure 9. It can be observed that numerical results are in good agreement with analytical solutions. Similar to the velocity image, the resonance phenomenon also exists in the amplitude of the wall shear stress, especially for as displayed in Figure 9(b). For different , the variation tendency of is identical with , such as and . The first peak value reaches the maximum value, and other peak value shows a decreasing trend on the whole. Meanwhile, the number of peaks gradually increases with . The first peak value decreases and moves to the left along the horizontal axis. It can be inferred that when is greater than a certain value, the cure will approximate to a monotone decreasing cure. Finally, after is greater than a certain value, the stress change curve will approach a monotonically decreasing curve, close to that of a Newtonian fluid.

Secondly, the influence of parameters on the wall shear stress is considered. For Figures 10 and 11, , and , are chosen, respectively. Figure 10 shows the same variation trend as Figure 9, and the resonance peak increases with . The difference is that the amplitude of is smaller compared with that in Figure 9. It can be concluded that the amplitude of varies with , and . Here, has the greatest impact compared to the other two parameters. The larger is, the faster decreases with little and . The curve corresponding to is basically monotonically decreasing with frequency for different .

As Figure 11 shows, when both and are small, the curves obtained for different are monotonically decreasing. And the larger the , the faster the magnitude decreases. This physical phenomenon fully illustrates that parameters and influence the appearance of the resonance phenomenon. The classical viscoelastic model can be considered as a special case of the fractional-order viscoelastic model, and this is theoretically consistent with Zhu’s findings.

Next, the parameters are chosen as and , and the influence of different order parameters on the wall shear stress is studied in detail. As shown in Figure 12(a), the shear stress with results in a resonance phenomenon at , and the peak value increases dramatically with at this frequency. Multiple peak values appear especially for , and the amplitude decreases with when the frequency satisfies . The similar phenomenon also exists in Figure 12(b), where the amplitude decreases with for . However, there is only one peak value at , and this indicates that the parameter almost has no influence on the resonance phenomenon. However, as decreases, the amplitude of the wall shear stress decreases very slightly at low frequencies, but there is a significant increase in the amplitude at high frequencies.

5. Conclusions

In the present paper, the difference scheme based on the differential equation with fractional derivatives is obtained. Then, the stability, convergence, and unique solvability are proven by the energy method. Finally, the finite difference method is verified with the Poiseuille flow and oscillatory flow in a circular pipe. The following conclusions are obtained: (1)The fractional mixed partial derivative in the differential equation is solved reasonably by integration, successfully avoiding the problem of low-order accuracy caused by discretizing the mixed derivative directly(2)The consistency of the numerical and exact solutions indicates that the method is reliable, and it is effective and concise to deal with applicable engineering problems. Moreover, according to the results in this paper, the finite difference method can be applied to differential equations with fractional mixed partial derivatives similar to Equation (13)(3)The finite difference method accurately captures the velocity overshoot phenomenon in the Poiseuille flow and resonance phenomena that appeared in the oscillatory flow. Based on these studies, it can be inferred that the same phenomena that occur in the classical viscoelastic fluid flow also occur in the fractional viscoelastic fluid flow. Otherwise, the above phenomena have a significant dependence on the order of fractional derivatives. Meanwhile, the classical viscoelastic model can be regarded as a special case of the fractional viscoelastic model. For example, the simulation of the oscillatory flow can well describe the physical deformation of materials, which is very helpful for the study of the properties of several materials in terms of damping

Appendix

A The Proof of Finite Difference Scheme

For Equation (13), the time partial derivative with the order is given:

Then, the following formula is obtained: where

Based on the following formulas,

Equation (A.2) is further simplified as

Supposing , we have

The first item in the square brackets on the right-hand side of Equation (A.6) is

Similar to the discrete method of the preceding item , the second term is discretized. Here,

In combination with Equations (A.2), (A.5), (A.7), and (A.8), the partial derivative of with order is given

Inserting Equations (15), (16), (A.2), and (A.8) into Equation (11), we obtain the expression at node .

where and is the local truncation error.

The expression (A.10) is further arranged as

B The Proof of Equation (20)

For the second item on the LHS of Equation (A.11), it is satisfied:

And the RHS of the above inequality is further simplified as

For the formula (B.2), the following inequality is established:

Therefore, formula (B.1) can be equally expressed as

Meanwhile, the equality is found for the RHS of Equation (B.4):

Then, the following result can be easily obtained:

Because is a negative, we can obtain

Similar to Equation (B.1), the third item on the RHS of Equation (A.10) can be dealt with

Taking into account Equations (B.7) and (B.8), it can be concluded that

Namely,

C The Derivation of Analytical Solution of the Oscillatory Flow in a Circular Pipe

Based on the hypothesis mentioned in Section 4, the model equation of the fractional Oldroyd-B fluid is equally expressed as

where is the component of shear stress and is the velocity component along z direction. The influence of the volume force is ignored, and the momentum equations with the nonslip boundary condition are obtained

In combination with Equations (C.2) and (C.3), the velocity equation along the direction is obtained:

In order to verify the numerical method, the analytical solution of the velocity and stress are given first. The expression of the pressure gradient is given:

The flow considered here has the characteristics of a simple harmonic oscillation, so all physical quantities change periodically like the frequency. Therefore, the velocity and shear stress are expressed, respectively, as

Substituting Equations (C.5) and (C.6) into Equation (C.4), the velocity equation is simplified as

where

Here, the new variables are introduced:

Then, Equation (C.8) is equally expressed as with the nonslip boundary condition

Obviously, Equation (C.11) is the first kind of zero-order Bessel equation, so that the exact solution of can be given as The coefficient is obtained with the nonslip boundary condition:

Furthermore, in combination with Equations (C.13) and (C.15), is solved and we can directly give the exact solution of the velocity:

With the same method, the exact solution of the shear stress for the oscillating flow is obtained:

In the following paper, the numerical simulation is executed with dimensionless analysis. Here, the dimensionless quantities are introduced:

Then, we can obtain

where is the steady velocity corresponding to the pressure gradient .

The dimensionless velocity is expressed as where .

Based on Equation (C.19), the amplitude of dimensionless velocity is obtained:

With the same method, the dimensionless shear stress and its amplitude on the pipe wall are also obtained:

Nomenclature

:The velocity vector, and , , and are, respectively, along the , , and directions
:The uniform density
:The pressure
:The Cauchy stress tensor
:The extra stress tensor
:The shear stress component along the -axis in the plane corresponding to the -axis
:The material derivative of time
:The unit tensor
:The unit vector along the direction
:The first Rivlin-Ericksen tensor
:Interpolation coefficient
:The dimensionless analytical amplitude of the velocity
:Pipe radius
:The fractional differential operator
:The order of the fractional derivatives
:The relaxation time
:The retardation time
:The viscosity
:Gamma function
:The discrete point of the mesh
, :The space step and the time step
:Grid node
:The gradient of the velocity vector
:The transpose form of
:The error between numerical result and exact solution at the -th time step
:The interpolation operator
:The width of two parallel plates
:The dimensionless analytical amplitude of the stress tensor
:Oscillation frequency.

Data Availability

All the data supporting the finding of the present study are included within the article.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This research was supported by the Natural Science Foundation of Ningxia (no. 2021AAC03209) and the research project of the National Key Laboratory of Science and Technology on Aerodynamic Design and Research (no. 614220120030115).