• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    Compressor geometric uncertainty quantification under conditions from near choke to near stall

    2023-04-22 02:05:12JunyingWANGBotongWANGHeliYANGZhenzhongSUNKiZHOUXinqinZHENG
    CHINESE JOURNAL OF AERONAUTICS 2023年3期

    Junying WANG, Botong WANG, Heli YANG, Zhenzhong SUN,Ki ZHOU, Xinqin ZHENG,,*

    a State Key Laboratory of Automotive Safety and Energy, School of Vehicle and Mobility, Tsinghua University, Beijing 100084, China

    b Department of Aerodynamics and Thermodynamics, Institute for Aero Engine, Tsinghua University, Beijing 100084, China

    KEYWORDS Compressor;Geometric uncertainty quantification;Interpretable machine learning;Multiple conditions;Neural network

    Abstract Geometric and working condition uncertainties are inevitable in a compressor,deviating the compressor performance from the design value.It’s necessary to explore the influence of geometric uncertainty on performance deviation under different working conditions.In this paper,the geometric uncertainty influences at near stall,peak efficiency, and near choke conditions under design speed and low speed are investigated.Firstly,manufacturing geometric uncertainties are analyzed.Next,correlation models between geometry and performance under different working conditions are constructed based on a neural network.Then the Shapley additive explanations (SHAP)method is introduced to explain the output of the neural network.Results show that under real manufacturing uncertainty,the efficiency deviation range is small under the near stall and peak efficiency conditions.However, under the near choke conditions, efficiency is highly sensitive to flow capacity changes caused by geometric uncertainty, leading to a significant increase in the efficiency deviation amplitude, up to a magnitude of -3.6%.Moreover, the tip leading-edge radius and tip thickness are two main factors affecting efficiency deviation.Therefore, to reduce efficiency uncertainty,a compressor should be avoided working near the choke condition,and the tolerances of the tip leading-edge radius and tip thickness should be strictly controlled.

    1.Introduction

    As the core component of an aero-engine or gas turbine, the aerodynamic performance of a compressor has an important influence on the engine performance.However, due to manufacturing errors or in-service degradation,1there are a lot of unavoidable geometric variable deviations in a compressor,which have strong uncertainty characteristics.The geometric uncertainties lead to the performance deviating from the design value and presenting a probability distribution.Moreover, with continuous improvement of compressor performance indices, blade loading and therefore the sensitivity of performance to blade geometry are increasing.Therefore, it is becoming more and more important to investigate the influence of geometric uncertainty on performance and find out critical geometric variables2.

    The geometric uncertainty quantification problem in compressors has characteristics of high randomness, high dimensionality, and strong nonlinearity.In relative research, it is a core difficulty to evaluate the quantitative relationship between geometry and performance efficiently and accurately.To solve this problem, in recent years, a series of Uncertainty Quantification(UQ)methods,including the direct sampling method,3–5the sensitivity-based method,6,7and the surrogate model-based method,8–11has achieved rapid development and application.Among all the methods, the surrogate model-based method,which uses a fast correlation model as a low computational cost alternative to a high-fidelity model such as Computational Fluid Dynamics (CFD), has the dual advantages of saving computing resources and retaining high accuracy of model results.In terms of the selection of a surrogate model, there exists a trade-off between the complexity and interpretability of the model.12As a highly nonlinear model,a neural network can predict performance accurately.However,due to the opacity of the internal structures of neural networks, usually they can only be used as‘‘black boxes”.To solve this problem,Shapley additive explanations (SHAP), an interpretable machine learning method based on game theory,13,14was firstly introduced to turbomachinery research by Ref.15to explore the multi-geometrical variable coupling effect in a transonic compressor.In this study, the SHAP method was adopted to extract the influence rule and clarify sensitive uncertainty variables from a highly nonlinear neural network.On one hand, results from SHAP analysis can provide guidance for compressor tolerance design and control;on the other hand,they can help explore the influence mechanism of critical geometric uncertainty variables on performance.

    There has been a series of studies on compressor geometric uncertainty quantification based on different UQ methods.Based on the Latin Hypercube Sampling (LHS) method,Lange et al.16explored the influences of 18 manufacturing geometric uncertainty variables of a 1.5-stage high-pressure compressor blade on performance.Results showed that the leading-edge thickness was the most important geometric parameter affecting efficiency.Based on the same method,similar studies were carried out on asymmetric compressors17and multistage axial compressors.18Based on the adjoint method,Zhang et al.19evaluated the influences of geometric uncertainties on the compressor performance.Moreover, Ju and Zhang20used the Support Vector Regression (SVR) model as a surrogate model to rapidly predict the cascade performance under pollution conditions.Most of the existing studies are quantitative analysis under a single working condition.However,due to different matching criteria or operating environments, a compressor may be used under different working conditions at different rotational speeds.Under different working conditions,compressor mass flow rates,loading characteristics, and internal flow field structures are different, so the influence of geometry on performance may also have significant differences.Schnell et al.21studied the distribution of the outflow angle β2caused by geometric uncertainties under different incident angles β1based on a 2-dimensional blade profile.Results showed that, with the same geometric uncertainty, the deviation range of β2showed an obvious growth trend with β1.At present, research on geometric uncertainty quantification under different working conditions is still very limited, and it is not very clear how geometry and different flow field structures interact to influence performance.

    The motivation of this paper is to systematically explore the differences in the influences of real manufacturing uncertainties of geometric variables at different speeds and under different working conditions.Based on the results, guidance for compressor design and tolerance design and control can be provided.The rest of the paper is organized as follows.Firstly,the distribution characteristics of geometric variables with manufacturing uncertainty are extracted from measured data.Then, the uncertainty quantification method adopted in this paper is introduced, including CFD, fully connected neural network, and the SHAP method.Afterward, the influences of geometric uncertainties are analyzed and discussed, while performance distribution and parameter sensitivity are compared under different working conditions.Finally, the flow mechanism behind the influence of the critical uncertainty variable is discussed based on CFD simulation results.

    2.Case overview

    2.1.Compressor geometry and parameterization

    The research object of this study is the first-stage rotor of a three-stage compressor, and the illustration of the computation domain is presented in Fig.1.

    To study the influences of different geometric uncertainty variables on performance, the rotor blade geometry should be parameterized first, as shown in Fig.2.Firstly, the number of blades and the meridional profile of hub and shroud are parameterized.Secondly, five spanwise sections with 0, 0.19,0.51,0.83,and 1.00 spanwise are parameterized,and the stacking lines between different sections are defined.Finally, each spanwise section of the blade is defined by parameters including stagger angle G, leading-edge radius Rle, trailing-edge radius Rte, camber curve shape, and blade thickness distribution T1-T5.

    2.2.Compressor geometric manufacturing uncertainty

    Fig.1 Illustration of compressor rotor computation domain.

    Fig.2 Schematic of geometrical uncertainty variables of compressor rotor.

    For the compressor rotor,the three-coordinate measuring data of more than 100 new blades manufactured are statistically analyzed.For each blade, the profiles in 0.19, 0.51, and 0.83 spanwise sections are measured and represented by the coordinates of blade profile points.An automatic parameterization program is written to extract the geometric variables (Rle,Rte, T, and G) from coordinate points.The thickness T,leading-edge radius Rle, and trailing-edge radius Rteare the radii of the corresponding inner circles of the blade profile.To calculate the stagger angle G, firstly, the center points of the inner circles are connected as the camber curve.Then,two intersection points of the camber curve and the blade profile are defined as the leading/trailing edge points.Finally, the angle between the line formed by the leading/trailing points and the axial direction is defined as the stagger angle G.

    The absolute deviations of the measured data from the design values are marked as ΔRle,ΔRte,ΔT,and ΔG(note that the range of ΔT is the union set of the ranges of ΔT1-ΔT5).Table 1 lists the variation ranges of the deviation values at 0.19, 0.51, and 0.83 spanwise sections (marked as L, M, and H sections, respectively).It can be seen that: the relative deviations of Rleand Rteare very obvious and on the order of ±50% relative to the design values; the relative deviations of T are smaller and on the order of ±6% relative to the design values; the absolute deviation of the stagger angle is on the order of ±0.1°.

    In the following study, the variation ranges of the 0.19,0.51, and 0.83 spanwise sections are given directly based onmeasured results, while in other sections,geometric deviations are given by B-spline interpolation.This interpolation can be considered reasonable referring to Refs.,16,22which present smooth distributions of geometric deviations along blade heights.

    Table 1 Geometric variable deviation ranges of rotor due to manufacturing uncertainty.

    3.Uncertainty quantification method

    3.1.Method overview

    Based on the statistical analysis results of compressor geometric uncertainty, the influence of uncertainty is quantified and analyzed through the following three main modules:

    I Design of Experiment (DoE).3In this study, it is assumed that geometric uncertainty variables conform to uniform distribution within the variation range, and the LHS3method is adopted to generate 365 blade samples with different geometric variable values.Note that for each sample, ΔT1-ΔT5are assumed the same and marked as ΔT.Therefore, there are totally 12 independent geometric uncertainty variables in this study.The total number of samples is more than 30 times the total number of geometric uncertainties, which ensures the validity of sampling.For each sample, its performance is obtained through CFD simulation, and the detailed settings of the CFD method will be described in Section 3.2.

    II Surrogate model construction.Based on the database generated in the DoE process, the fully connected neural network23is trained to map the correlation between geometric variables and compressor performance, and is used as the surrogate model of 3-dimensional CFD simulation to rapidly predict compressor performance.The detailed setting of the neural network will be described in Section 3.3.

    III Sensitivity analysis.Based on the neural network surrogate model, the interpretable machine learning method SHAP can be used to extract the influence rule and clarify sensitive uncertainty variables from the highly nonlinear neural network.The principles of the SHAP method will be described in Section 3.4.

    Based on Modules I and II, the distribution of compressor performance due to uncertainty can be obtained.This part of study aims at providing guidance for compressor design.Based on Module III, critical geometric variables can be extracted from a large number of variables.This part of study aims at providing guidance for compressor tolerance design and control.

    3.2.CFD method and validation

    FINE/Turbo (Version 15.1) software is used for CFD simulation.The computation is set as steady and compressible.The multi-grid procedure is adopted to improve the convergence rate of simulation,and the Spalart-Allmaras turbulence model is adopted.24The computational domain is a single rotor passage with periodic boundary conditions.At the inlet, the total pressure (101325 kPa), total temperature (288.15 K), and velocity direction (given according to the mean outlet flow angle of the inlet guide vane of the 3-stage compressor) are defined as the boundary condition.At the outlet, under complete choking conditions, the static pressure is defined, while under other conditions, the mass flow rate is defined as the boundary condition.Moreover, all the solid walls are set as non-slip wall boundaries.The blade and hub surfaces are set as rotating, while the shroud surface is set as stationary.

    The software AutoGrid is used to discretize the computational domain.To quantify the discrete error, mesh independence analysis is conducted.The number of grid nodes in all dimensions increases to about 1.2–1.5 times uniformly when the mesh is refined.For all the meshes,the thickness of the first layer grid is fixed as 0.003 mm, corresponding to an average y+value below 3 and a maximum y+value below 10.The details of the five meshes are summarized in Table 2.

    CFD simulations under 100 % design speed are conducted using the five sets of meshes,and Fig.3 shows how the normalized choke mass flow rate mc,Nand normalized peak efficiency ηp,Nvary with the total number of grid points ng.Based on Fig.3, the Fine grid is selected in this study considering the trade-off between computational accuracy and cost.Compared with the UltraFine results,the relative error of the choke mass flow rate and the absolute error of the peak efficiency are -0.01% and -0.11%, respectively.

    CFD simulation validation is conducted based on the experimental results of the 3-stage compressor.It needs to be noted that, in the simulation of the 3-stage compressor, the above simulation settings are adopted, and the mesh grid settings of each blade row are basically consistent with that of the above Fine grid.Table 3 presents the compressor performance errors of CFD simulation compared with experiment results, which shows that the simulation can underestimate the choke mass flow rate (mc) and peak efficiency (ηp) of the 3-stage compressor.The relative error of the choke mass flow is -0.62% and the absolute error of the peak efficiency is -0.41%, which meet the numerical simulation accuracy requirement in this study.

    Fig.3 Compressor performance simulated with different mesh sizes.

    Table 3 Comparison of 3-stage compressor performance between experimental and CFD results.

    It is worth noting that 0.41%is the CFD error of the whole 3-stage compressor instead of the single rotor alone.However,in the subsequent UQ simulation and analysis, the research object is only a single rotor blade (R1).Generally speaking,the fact that CFD errors for multistage compressors are significantly higher than that of a single blade row25,26is mainly resulted from the following two reasons: (A) errors induced by the mixing plane setting at the rotor/stator interface27; (B)stage-by-stage amplification from upstream to downstream blade rows.27,28In addition,CFD should be used on a comparative basis, especially in the situation of UQ or optimization study.27Refs.29,30have demonstrated the ability of CFD to capture the relative deviation of performance change, bycomparing the CFD and experimental results of compressors with geometric deviations.

    Finally, the performance map of the compressor rotor,including the pressure ratio-mass flow rate and efficiencymass flow rate performance curves at 100% and 70% design rotational speeds, is obtained through CFD simulation, as shown in Fig.4.In the figure,the mass flow rate(m),pressure ratio(π),and efficiency (η)are normalized.The normalization methods are as follows:

    In the following part of this study, on one hand, at 100%design speed with supersonic/transonic flow fields,three working conditions are selected:near stall point(A),peak efficiency point (B), and near choke point (C).On the other hand, at 70 % design speed with subsonic flow fields, three working conditions are selected: near stall point (D), highest efficiency point(E),and near choke point(F).For each of the 6 working conditions, the influence of geometric uncertainty on performance will be quantitatively analyzed.

    3.3.Fully connected neural network

    The fully connected neural network is a widely used model to extract nonlinear mapping from a large amount of data.Generally,a neural network consists of an input layer,several hidden layers, and an output layer.Each layer has several neurons, which are connected to the neurons of the upstream and downstream layers in a linear combination.The nonlinearity of the model is introduced by the activation function of each intermediate neuron.The structure of the neural network model mimics that of the human brain.

    In this study,the neural network establishment and training process is based on the PyTorch framework.The neural network consists of one input layer, two hidden layers, and one output layer, as shown in Fig.5.The number of nodes in the input layer is 12 (i.e., nin= 12), where each node corresponds to a geometric uncertainty variable.There are 6 nodes in the output layer (i.e., nout= 6), where each node corresponds to the efficiency of a working condition.The hyperbolic tangent function is chosen as the activation function.The root mean square error is used as the loss function.For all 365 samples,300 samples are used as the training set and 65 samples are used as the validation set.The training process of the neural network is driven by the adaptive moment estimation optimizer, and the number of iterative steps is set to 20000, which is sufficient to find the optimal parameters of the neural network to make the model with the highest accuracy.More details on the structure and algorithms of neural networks can be found in Ref.31.

    The node numbers of two hidden layers (i.e., nh1and nh2)are the hyperparameters, which need to be specified before neural network training.To determine the optimal hyperparameters,a grid search is performed for nh1and nh2.Moreover,the K-fold cross-validation(K=6 to find the balance between bias and variance) and determination coefficients (R2) are adopted to evaluate the performance of neural networks.The concrete steps are:

    Step 1.The original training set is partitioned into K subsets(folds) of size 365/K.

    Step 2.For each of the K fold, a neural network is trained on the union of the other K - 1 folds (training set).Then the error of its output is estimated using the fold (validation set).In this study, the determination coefficient (R2) between the ηNvalues predicted by CFD and the neural network is chosen as the evaluation index of the network error.

    Step 3.After Step 2, K neural networks are trained and R2on the validation set is calculated K times.Then the average R2of all K folds, marked as mean(R2), is used to measure the accuracy of the neural network.

    Fig.6 presents the contour of the mean(R2) as nh1and nh2change.According to Fig.6, the point with nh1= 8 and nh2=10 is chosen as the optimal hyperparameter point,which is marked as a yellow star in the figure.Using the hidden layer nodes above, the mean(R2) of the neural network on the validation set is higher than 0.98, indicating that the neural network model has established an effective mapping between geometric uncertainty variables and the compressor efficiency.

    Fig.4 Rotor performances at 100% and 70% design rotational speeds.

    Fig.5 Schematic of fully connected neural network structure.

    In addition, when nh1= 8 and nh2= 10, the neural network of Fold 3 has the highest R2, therefore it is selected for the subsequent result analysis.To provide more information about the precision of the neural network, the Root Mean Squared Error (RMSE) on the validation set is calculated,and the RMSE values are 0.01%, 0.01%, 0.08%, 0.01%,0.01%, and 0.05% at Conditions A - F, respectively.This result further shows that the prediction accuracy of the neural network can meet the needs of this study.

    3.4.Shapley additive explanations

    To extract the critical geometric variables from a large number of variables and provide guidance for compressor tolerance design and control, sensitivity analysis is needed, as the geometric uncertainty quantification problem in compressors has the characteristics of high dimensionality and strong nonlinearity.In general, the sensitivity of a geometric variable is not uniform at different points in the uncertainty space.Moreover,the impact of one geometric variable on performance can be significantly affected by other geometric variables15.

    To provide accurate sensitivity analysis results and effective guidance for tolerance control, both local and global sensitivity indices are needed.Additionally, for each sensitivity index,both univariate and multivariable effects should be taken into consideration accurately.Facing such demands, traditional sensitivity analysis methods, such as using the coefficients of linear/polynomial regression as a sensitivity index,cannot provide comprehensive results which satisfy the above requirements.

    Therefore, in this study, an interpretable machine learning method named SHAP is introduced.The principle and mathematical derivation of the SHAP method are as follows:

    (1) Local interpretable model-agnostic explanations.

    There is a tradeoff between model complexity and model interpretability.On one hand,some interpretable models,such as linear regression models,decision trees,or rule lists,are simple in structure, but have insufficient prediction accuracy for complex problems.On the other hand, a well-trained neural network can provide sufficient predictions, but it’s difficult to interpret prediction results due to the structural complexity of the model.

    Fig.6 Contours of mean(R2) scores with different numbers of hidden layer nodes under Conditions A - F.

    That is, the neural network is usually used as a black box.For example, given any combination of geometric uncertainty variables,the neural network obtained in the above section can predict the corresponding performance.However, the contribution of each geometric uncertainty variable to performance uncertainty cannot be obtained directly.Facing such contradiction, a Local Interpretable Modelagnostic Explanations (LIME) method is proposed, and its main idea is using simple local surrogate models to explain the predictions of single instances of complex black box models.32A local surrogate model can be trained to approximate the predictions of an underlying black box model in the neighborhood of an instance to be explained.Mathematically, a local surrogate model with explanatory constraints can be expressed as follows:

    (2) SHAP method: local explanation for a single instance.

    SHAP13,15is a method used to explain the output of machine learning models based on the game theory.The goal of SHAP is to explain the prediction of a certain instance by calculating the contribution of each feature to the prediction.

    The SHAP method combines ideas from the game theory and the LIME method.Based on the game theory,the Shapley value of each input feature can be calculated and used as a measure of its contribution to the output for an instance.Shapley values tell us how to evenly distribute the prediction between features.An important feature of this method is that Shapley values are additive.Therefore, a linear model can be constructed between the output value and the Shapley values of each input features.Based on the idea of the LIME method,this linear model can be used as the local surrogate model g to explain the output of the black box model.

    In the SHAP method, the explanation of X is specified as

    where EYk(S ) is the prediction for feature values in the subset that are marginalized over features that are not included in S,i.e.,

    As can be seen from Eq.(7),the difference between the output prediction for a subset S with and without feature {xi} is used as a measurement of its contribution to the output.Therefore,for all subsets S,the corresponding δi(S )can be calculated.Then Eq.(6), the Shapley value calculation formula,can be perfectly derived by weighting and summing up all possible combinations of feature values.Finally, the obtained Shapley value can be used to measure the contributions of input features to the output.

    Since φiis obtained through the calculation of the marginal contribution of variable {xi}, φihas both univariate and multivariable effects.More details about how to extract the multivariable coupling effect from φican be found in Ref.15.

    To sum up,by adopting the SHAP method,both local and global sensitivity indices can be calculated accurately.Moreover,for each sensitivity index, both univariate and multivariable effects can be taken into consideration.Therefore, it is chosen in this study.

    4.Results and analysis

    4.1.Performance distribution analysis

    To obtain the performance distribution due to manufacturing geometry uncertainty, all geometric variables are assumed to be gaussian distributions, which are close to the distributions obtained by the real measured data.For each geometric variable,its mean value is the midpoint of the corresponding measured variation range in Table 1, and it is assumed that the measured variation range in Table 1 corresponds to ± 3σ(standard deviation).

    In addition, for any two geometric variables, the coefficients R2are calculated based on the measured data of more than 100 samples.Results show that for most groups,the coefficients R2are smaller than 0.5,indicating that the correlations between different geometric variables are not significant.Therefore, for these groups, the variables are assumed to be independent on each other.However, the coefficient R2between ΔGMand ΔGHis 0.52,and the coefficient R2between ΔTMand ΔTHis 0.51.Therefore,for each of these two groups,a multivariate Gaussian distribution is adopted as the joint probability distribution based on the real covariance matrix.Based on the above distribution, 20000 samples are generated based on the LHS method and gaussian distribution in the geometric uncertain space.

    For each of the samples, the corresponding performance is predicted by the neural network obtained in Section 3.3.The distributions of normalized efficiencies ηNof all samples are shown as the red distribution histograms in Fig.7.

    It can be indicated from Fig.7 that:

    (1) Under most of the conditions, the efficiencies of most samples are lower than the design value, and the efficiency decrease amplitude is greater than or equal to the efficiency increase amplitude, indicating that uncertainty will lead to the decline of average performance.

    (2) At 100% design rotational speed, the efficiency deviation range is small under the near stall and peak efficiency conditions caused by manufacturing geometric uncertainty,with a maximum efficiency decline amplitude of -0.3% and -0.4%,respectively, while under the near choke condition, the efficiency is highly sensitive to flow capacity changes caused by geometric uncertainty, leading to a significant increase in the efficiency deviation range, with a maximum decrease of -3.9%.Therefore, to reduce the efficiency deviation range caused by geometric uncertainty, the compressor should be avoided working near the choke condition.

    (3)At 70%design rotational speed,the efficiency deviation range is small under the near stall and peak efficiency conditions caused by manufacturing geometric uncertainty, with a maximum efficiency decline amplitude of -0.2%and -0.2%, respectively, while under the near choke condition, the efficiency is highly sensitive to flow capacity changes caused by geometric uncertainty, leading to a significant increase in the efficiency deviation range, with a maximum decrease of -2.5%.Therefore, to reduce the efficiency deviation range caused by geometric uncertainty, the compressor should be avoided working near the choke condition.

    (4) A comparison is conducted between the efficiency deviations at 100% design rotational speed with mainly supersonic/transonic flow field and at 70% design rotational speed with mainly subsonic flow fields in the compressor.Results show that despite there exist more complex flow structures such as shock wave in the supersonic/transonic flow field,which do not exist in the subsonic flow fields,there is no significant difference between the efficiency deviation ranges under high and low speeds.

    4.2.Sensitivity analysis

    Under all the 6 working conditions, for each of the 365 samples,the Shapley value of each geometric variable is calculated and used as a measurement of its contribution to the efficiency deviation.Fig.8 presents the Shapley value information of all samples under each working condition.The bar chart on the left is the SHAP feature importance figure, and the horizontal axis represents the mean|Shapley| value, which is the absolute mean of Shapley values for all samples and can be used to measure the average contribution of the geometrical variables on the total efficiency deviation.In addition, the scatter plot on the right is the Shapley value summary figure for 365 samples.Each point in the scatter plot represents a sample, with the color of the point representing the geometric feature value and the axial position of the point representing the corresponding Shapley value.

    As can be seen from the SHAP feature importance plots on the left of the 6 working conditions in Fig.8, under near stall Conditions A and D and peak efficiency Conditions B and E,the blade tip leading-edge radius Rle,Hplays a dominant role in efficiency deviation.In working Conditions A and D,the average contribution rate of Rle,Hto total efficiency deviation can reach 34%and 41%,respectively,while in working Conditions B and E, the average contribution rate of Rle,Hto total efficiency deviation can reach 49% and 54%, respectively.More importantly, in the near choke Conditions C and F, the sensitivity of the blade tip thickness THincreases significantly,which together with the leading-edge radius Rle,Hbecomes the two main factors affecting efficiency deviation.In working Conditions C and F,the sum of the average contribution rates of THand Rle,Hto total efficiency deviation can reach 36%and 25%, respectively.In addition, the influences of the stagger angle and thickness are also important under near stall Conditions A and D and near choke Conditions C and F, but not obvious under peak efficiency Conditions B and E.However,the trailing-edge radius of all spanwise sections has almost zero influence on efficiency under all working conditions.Therefore, among a variety of geometric variables, the tip leadingedge radius and tip thickness should be strictly controlled in terms of manufacturing tolerances.

    Fig.7 Distribution of compressor efficiency deviation considering geometric uncertainty under Conditions A - F.

    Fig.8 SHAP feature importance (left) and Shapley value summary (right) under Conditions A - F.

    According to the Shapley value summary plots on the right of the 6 working conditions in Fig.8, the Shapley value range of geometric features with greater significance is wider than that of insensitive geometric features.In addition, by investigating the relationship between feature values and Shapley values,most geometric features show an obvious monotone effect on efficiency.For example, for the leading-edge radii Rle,H,Rle,M, and Rle,L, the Shapley value decreases with an increase of the radius value in almost all working conditions, which means that when the leading-edge radius increases, the efficiency will be reduced.What’s more, in the choke conditions,efficiency also presents a significant trend of decline with an increase of the stagger angle (GH, GM, GL) or the thickness(TH, TM, TL).For a specific compressor case with unqualified performance,the local Shapley value can be used to guide how to repair and make it qualify.

    4.3.Flow mechanism analysis

    4.3.1.Comparison of influence mechanisms under different working conditions

    To explore the physical mechanism behind the influence of the critical geometric uncertainty variable on efficiency, two samples with a maximum and a minimum value of Rle,Hare added,with other geometric variables kept as the design values.Simulation results of the two samples under different working conditions at 100%and 70%design rotational speeds are shown in Fig.9.It can be indicated that at both 100 % and 70 %design rotational speeds, an increase of Rle,Hwill lead to a decrease of the choke flow rate and efficiencies, and the efficiency reduction amplitude in the near choke condition is significantly higher than those in the peak efficiency and near stall conditions.This result corresponds to the Shapley value of Rle,H under different working conditions in Fig.8.

    Based on the CFD results, a flow field analysis is conducted.To explore the mechanisms behind the efficiency decline, the entropy loss coefficient ζsis adopted33, which is defined as

    Fig.9 Influence of blade tip leading-edge radius on compressor efficiency with 100% and 70% rotational speeds.

    Fig.10 shows the distribution of the efficiency deviation coefficient ΔηMmax{Rle,H}-min {Rle,H} along the streamwise direction.The streamwise locations LEhuband TEhubare the leading-edge and trailing-edge locations at the hub section of the rotor blade, respectively.The black, red, and blue curves in the figure correspond to working Conditions A, B, and C,respectively.Through Fig.10, the streamwise position which is obviously affected by a change of Rle,Hcan be located.It can be seen that in Conditions A and B, the efficiency reduction mainly occurs around the leading edge of the blade, and the position with the largest efficiency deviation is denoted as Section I, while in working Condition C, on one hand,the efficiency decrease occurs around the leading edge of the blade; on the other hand, it occurs near the rear part of the blade.The endpoint of the region with the largest increase in the efficiency deviation coefficient is marked as Section II.

    To further investigate the mechanism of efficiency deviation under different working conditions, the detailed flow field at the 0.80 spanwise section is analyzed.Fig.12 presents the contour of the relative Mach number (Ma) at the 0.80 spanwise section.It’s indicated that when Rle,Hincreases, the strength of the detached shock wave near the leading edge is enhanced significantly, resulting in a higher shock wave loss, which corresponds to the increase of the efficiency deviation near the leading edge in Fig.10.

    Fig.13 presents the distribution of the static pressure rise coefficient Cpalong the blade surface at the 0.80 spanwise section, and the definition of Cpis

    where p01and p1are the average total and static pressure at the compressor inlet, respectively.

    In Fig.13, M’is the normalized coordinate of the streamwise position from the blade inlet to outlet at the 0.80 spanwise section.It can be indicated that under Conditions A,B,and C,when Rle,Hchanges, the pressure distribution near the leading edge changes significantly, corresponding to loss changes caused by acceleration and deceleration process changes near the pressure surface.In addition, under Conditions B and C,when Rle,Hincreases,the intensity of the spike near the suction surface increases, leading to an increase of the boundary layer loss.34More importantly, it can be seen from Fig.12 and Fig.13 that for Condition C, when Rle,Hincreases, the shock wave moves downstream, and the Ma after the passage shock wave increases obviously.From the general point of view of transonic rotor design, when the leading-edge radius changes,the area distribution along the streamwise direction changes,and then the flow acceleration and deceleration process in the passage changes, resulting in a change of the shock wave strength and a flow loss.Meanwhile, the wake region near the blade trailing edge increases due to a change of Ma after the shock wave,leading to an increase of the blade profile loss.

    Fig.11 Efficiency deviation distributions along spanwise direction at Sections I and II under Conditions A, B, and C.

    Fig.12 Relative Mach number contours at 0.80 spanwise section under Conditions A, B, and C.

    Fig.13 Cp distribution along blade surface at 0.80 spanwise section under Conditions A, B, and C.

    Fig.14 Efficiency deviation distribution along streamwise direction under Conditions B and E.

    Fig.15 Efficiency deviation distribution along spanwise direction at Section III under Conditions B and E.

    4.3.2.Comparison of influence mechanisms under different rotational speeds tion.The black and red curves in the figure correspond to working Conditions B and E respectively.Through Fig.14, the streamwise position which is obviously affected by a change of Rle,Hcan be located.It can be seen that for Condition B at 100% rotational speed, the efficiency reduction starts from a long distance upstream from the blade leading edge and lasts to a short distance downstream from the leading edge.For Condition E at 70% rotational speed,the efficiency reduction starts relatively downstream and only occurs in a short distance upstream and downstream from the leading edge.The position where the efficiency deviation is almost the largest in both conditions is marked as Section III.

    Figs.16(a) and 16(b) present the distributions of the static pressure rise coefficient Cpalong the blade surface at the 0.80 spanwise section under Conditions B and E, respectively.In the figure, M’is the normalized coordinate of the streamwise position from the blade inlet to outlet at the 0.80 spanwise section.It can be indicated that under Conditions B and E, when Rle,Hchanges, the pressure distribution near the leading edge changes significantly, corresponding to loss changes caused by acceleration and deceleration process changes near the pressure surface.Fig.16(c) presents the contour of the relative Mach number (Ma) at the 0.80 spanwise section under Condition E.Results show that different from working Condition B at 100% rotational speed, there is no shock wave near the leading edge under Condition E at 70% rotational speed.It can be seen from Fig.12 that in Condition B, there is a detached shock wave near the leading edge, and the shock wave strength varies significantly with Rle,H, which could explain the difference between the starting points of efficiency reduction under Conditions B and E in Fig.14.

    Fig.16 Cp distribution along blade surface and Ma contours at 0.80 spanwise section under Conditions B and E.

    5.Conclusions

    In this paper,the geometric uncertainty influences at near stall,peak efficiency, and near choke conditions under design speed and low speed are investigated, based on the real manufacturing geometric uncertainty of a compressor rotor blade.Conclusions are shown as follows:

    (1) At design speed and low speed, the efficiency deviation ranges caused by manufacturing geometric uncertainties are small under the near stall and peak efficiency conditions.For the case studied in this paper, the maximum efficiency decline amplitude is on the order of -0.4%.Meanwhile, under the near choke condition, the efficiency is highly sensitive to flow capacity changes caused by geometric uncertainty, leading to a significant increase in the efficiency deviation range.For the case studied in this paper, the maximum efficiency decline amplitude is on the order of -3.9%.Therefore, to reduce the efficiency deviation range caused by geometric uncertainty,the compressor should be avoided working near the choke condition.

    (2) Under the near stall and peak efficiency conditions, the blade tip leading-edge radius plays a dominant role in efficiency deviation.For the case studied in this paper,the average contribution rate of the tip leading-edge radius to total efficiency deviation can reach 34%-54%.Meanwhile,in the near choke conditions, the sensitivity of the blade tip thickness increases significantly, which together with the tip leadingedge radius becomes the two main factors affecting efficiency deviation.For the case studied in this paper, the sum of the average contribution rates of the tip leading-edge radius and blade tip thickness to total efficiency deviation can reach 25%-36%.Therefore,among a variety of geometric variables,the tip leading-edge radius and tip thickness should be strictly controlled in terms of manufacturing tolerances.

    (3) A comparison is conducted between the efficiency deviations at 100% design rotational speed with mainly supersonic/transonic flow field and at 70% design rotational speed with mainly subsonic flow fields in the compressor.Results show that despite there exist more complex flow structures such as shock wave in the supersonic/transonic flow field,which do not exist in the subsonic flow fields,there is no significant difference between the efficiency deviation ranges under high and low speeds.

    Declaration of Competing Interest

    The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

    Acknowledgement

    This research was supported by the National Science and Technology Major Project, China (No.2017-II-0004-0016).

    99热这里只有精品一区| 亚洲精品国产色婷婷电影| 91成人精品电影| 免费不卡的大黄色大毛片视频在线观看| 国产在线视频一区二区| 国产精品蜜桃在线观看| 中文字幕精品免费在线观看视频 | 国产一区有黄有色的免费视频| 在现免费观看毛片| √禁漫天堂资源中文www| 老熟女久久久| 午夜福利影视在线免费观看| 男女高潮啪啪啪动态图| 国产精品国产三级专区第一集| 汤姆久久久久久久影院中文字幕| 人妻夜夜爽99麻豆av| 精品人妻一区二区三区麻豆| 久久青草综合色| 国产免费现黄频在线看| 在线观看www视频免费| 91久久精品国产一区二区成人| 亚洲内射少妇av| 亚洲内射少妇av| 狂野欧美激情性xxxx在线观看| 大码成人一级视频| 亚洲精品色激情综合| 国产精品一国产av| 高清黄色对白视频在线免费看| 老司机影院毛片| 性高湖久久久久久久久免费观看| 久久久国产精品麻豆| 国产午夜精品一二区理论片| 在线免费观看不下载黄p国产| 欧美日韩国产mv在线观看视频| 视频区图区小说| 久久精品人人爽人人爽视色| 伦精品一区二区三区| 中文字幕免费在线视频6| 热99国产精品久久久久久7| 熟女人妻精品中文字幕| 中文字幕人妻熟人妻熟丝袜美| 精品亚洲乱码少妇综合久久| videos熟女内射| 成人二区视频| 一级毛片我不卡| 国产精品人妻久久久影院| 午夜福利网站1000一区二区三区| 免费大片黄手机在线观看| 中文字幕人妻熟人妻熟丝袜美| 日本黄色片子视频| 亚洲欧美一区二区三区国产| 18禁在线播放成人免费| 26uuu在线亚洲综合色| 永久网站在线| 91精品一卡2卡3卡4卡| 亚洲图色成人| 午夜91福利影院| 午夜久久久在线观看| 成人黄色视频免费在线看| 国产av码专区亚洲av| 99久久精品一区二区三区| 99热6这里只有精品| 制服人妻中文乱码| 少妇人妻 视频| 成人毛片a级毛片在线播放| 午夜免费男女啪啪视频观看| 大香蕉97超碰在线| 人人妻人人澡人人看| 国产成人精品福利久久| 考比视频在线观看| 久久人人爽人人片av| 高清不卡的av网站| 大香蕉久久网| av黄色大香蕉| a级毛片黄视频| 毛片一级片免费看久久久久| 欧美另类一区| 高清av免费在线| 一级片'在线观看视频| 免费不卡的大黄色大毛片视频在线观看| a级毛片在线看网站| 边亲边吃奶的免费视频| 亚洲成色77777| 午夜免费男女啪啪视频观看| 亚洲精品乱码久久久v下载方式| 欧美 亚洲 国产 日韩一| 国产视频首页在线观看| 3wmmmm亚洲av在线观看| 久久精品人人爽人人爽视色| 丁香六月天网| 一边摸一边做爽爽视频免费| 国产精品久久久久久av不卡| 美女内射精品一级片tv| 日本av手机在线免费观看| 日韩免费高清中文字幕av| 一本大道久久a久久精品| 午夜免费鲁丝| 九色成人免费人妻av| 色婷婷av一区二区三区视频| 少妇人妻久久综合中文| 两个人的视频大全免费| 你懂的网址亚洲精品在线观看| 伦理电影大哥的女人| 亚洲精品乱码久久久v下载方式| 国产精品国产三级国产专区5o| 亚洲高清免费不卡视频| 狂野欧美激情性xxxx在线观看| 久久久久视频综合| 多毛熟女@视频| 大片免费播放器 马上看| 午夜激情福利司机影院| 99久久综合免费| 久热这里只有精品99| 一区二区日韩欧美中文字幕 | 亚洲国产欧美在线一区| 亚洲综合精品二区| 久久久国产欧美日韩av| 亚洲av二区三区四区| 日本av免费视频播放| 丰满乱子伦码专区| 美女国产高潮福利片在线看| 久久这里有精品视频免费| 亚洲欧美精品自产自拍| 2018国产大陆天天弄谢| 精品少妇内射三级| 晚上一个人看的免费电影| 色5月婷婷丁香| 日韩av在线免费看完整版不卡| 日韩大片免费观看网站| 欧美少妇被猛烈插入视频| 日韩不卡一区二区三区视频在线| 午夜免费观看性视频| 亚洲三级黄色毛片| 精品人妻一区二区三区麻豆| 中文字幕久久专区| 赤兔流量卡办理| 高清视频免费观看一区二区| 欧美老熟妇乱子伦牲交| av线在线观看网站| 韩国高清视频一区二区三区| 99九九在线精品视频| 亚洲精品国产av成人精品| 热99国产精品久久久久久7| 插阴视频在线观看视频| 日韩成人av中文字幕在线观看| 亚洲,一卡二卡三卡| 成人国语在线视频| 天美传媒精品一区二区| 我要看黄色一级片免费的| 亚洲精品中文字幕在线视频| 日韩一区二区视频免费看| 亚洲丝袜综合中文字幕| 久久国产精品男人的天堂亚洲 | 黑丝袜美女国产一区| 久久人妻熟女aⅴ| 亚洲伊人久久精品综合| 精品国产一区二区三区久久久樱花| 看免费成人av毛片| 久久精品国产亚洲av涩爱| 人人妻人人爽人人添夜夜欢视频| 久久亚洲国产成人精品v| 夜夜爽夜夜爽视频| 成人综合一区亚洲| 亚洲精品一区蜜桃| 国内精品宾馆在线| 又粗又硬又长又爽又黄的视频| 免费大片18禁| videossex国产| 男人爽女人下面视频在线观看| 97精品久久久久久久久久精品| 美女cb高潮喷水在线观看| 大片免费播放器 马上看| 精品久久久精品久久久| 精品少妇久久久久久888优播| 欧美精品国产亚洲| 校园人妻丝袜中文字幕| 国产免费视频播放在线视频| 欧美 亚洲 国产 日韩一| 精品人妻熟女av久视频| 99热全是精品| 日韩欧美精品免费久久| av在线播放精品| 国产精品一区www在线观看| 日韩欧美一区视频在线观看| 亚洲怡红院男人天堂| 女性被躁到高潮视频| 人妻一区二区av| 免费播放大片免费观看视频在线观看| 中文字幕最新亚洲高清| 夜夜骑夜夜射夜夜干| 熟女人妻精品中文字幕| 久久久a久久爽久久v久久| 人人妻人人澡人人爽人人夜夜| 日韩精品有码人妻一区| av在线播放精品| 久久人人爽人人片av| 在线播放无遮挡| 日韩av在线免费看完整版不卡| 一本—道久久a久久精品蜜桃钙片| 一二三四中文在线观看免费高清| 在线 av 中文字幕| 最近的中文字幕免费完整| 纯流量卡能插随身wifi吗| 韩国av在线不卡| 香蕉精品网在线| 国产精品蜜桃在线观看| 三上悠亚av全集在线观看| 亚洲av综合色区一区| 亚洲精品成人av观看孕妇| 亚洲国产色片| 老司机影院成人| 一本大道久久a久久精品| 能在线免费看毛片的网站| 精品国产国语对白av| 美女xxoo啪啪120秒动态图| xxxhd国产人妻xxx| 精品亚洲成国产av| 波野结衣二区三区在线| 80岁老熟妇乱子伦牲交| 在线播放无遮挡| 国产日韩欧美在线精品| 久久这里有精品视频免费| 亚洲人成77777在线视频| 黄片无遮挡物在线观看| 汤姆久久久久久久影院中文字幕| 亚洲成色77777| 欧美老熟妇乱子伦牲交| 80岁老熟妇乱子伦牲交| kizo精华| av专区在线播放| 国产精品人妻久久久久久| 久热这里只有精品99| 亚洲精品美女久久av网站| 免费看光身美女| 人妻一区二区av| 中文字幕人妻熟人妻熟丝袜美| 日日啪夜夜爽| 成人综合一区亚洲| 高清视频免费观看一区二区| 在线天堂最新版资源| 少妇高潮的动态图| 久久精品国产亚洲av涩爱| 人妻制服诱惑在线中文字幕| 亚洲av欧美aⅴ国产| 日产精品乱码卡一卡2卡三| 免费人妻精品一区二区三区视频| 成人毛片60女人毛片免费| 人体艺术视频欧美日本| 美女大奶头黄色视频| 国产永久视频网站| 我的女老师完整版在线观看| 久久久国产一区二区| 9色porny在线观看| 卡戴珊不雅视频在线播放| 久久午夜福利片| 国内精品宾馆在线| 美女国产高潮福利片在线看| 久久精品国产亚洲网站| 在线播放无遮挡| 五月伊人婷婷丁香| 国产精品免费大片| 少妇猛男粗大的猛烈进出视频| 中文字幕最新亚洲高清| 午夜91福利影院| 国产成人aa在线观看| 中文字幕亚洲精品专区| 一级爰片在线观看| 久久鲁丝午夜福利片| 最新的欧美精品一区二区| 国产欧美日韩综合在线一区二区| 亚洲av男天堂| av免费在线看不卡| 视频中文字幕在线观看| 婷婷色综合www| 亚州av有码| av在线播放精品| 免费久久久久久久精品成人欧美视频 | 少妇的逼水好多| 三上悠亚av全集在线观看| 黄色欧美视频在线观看| 插逼视频在线观看| 亚洲精品乱码久久久久久按摩| 亚洲精品乱码久久久久久按摩| 少妇的逼好多水| 国产精品国产三级国产专区5o| 97在线人人人人妻| 国产av码专区亚洲av| 最近最新中文字幕免费大全7| 热re99久久国产66热| 国产欧美日韩一区二区三区在线 | 久久99蜜桃精品久久| 欧美国产精品一级二级三级| 麻豆精品久久久久久蜜桃| 麻豆成人av视频| 亚洲美女黄色视频免费看| 国产亚洲最大av| 国产69精品久久久久777片| 亚洲精华国产精华液的使用体验| 日产精品乱码卡一卡2卡三| 久久人人爽av亚洲精品天堂| 亚洲成人一二三区av| 免费人成在线观看视频色| 九九在线视频观看精品| 色吧在线观看| 曰老女人黄片| 国产精品人妻久久久久久| 日本色播在线视频| 丰满少妇做爰视频| 亚洲经典国产精华液单| 国产精品久久久久久av不卡| 成人亚洲精品一区在线观看| 欧美性感艳星| 国产乱人偷精品视频| 亚洲一区二区三区欧美精品| 久久99热这里只频精品6学生| 亚洲国产精品国产精品| 国产av精品麻豆| 国产极品粉嫩免费观看在线 | 欧美另类一区| 日韩一区二区三区影片| 少妇被粗大的猛进出69影院 | 少妇高潮的动态图| 国产一区二区三区av在线| 一级毛片黄色毛片免费观看视频| 国产精品国产三级专区第一集| 久久 成人 亚洲| 精品人妻在线不人妻| 国产一区二区在线观看av| 国产 精品1| 免费观看av网站的网址| 美女cb高潮喷水在线观看| 老司机影院毛片| 亚洲成人手机| 国产毛片在线视频| 搡老乐熟女国产| 亚洲人成网站在线观看播放| 国产高清国产精品国产三级| 日韩成人伦理影院| 国产精品女同一区二区软件| 亚洲性久久影院| 亚洲国产精品一区三区| 在线精品无人区一区二区三| 91在线精品国自产拍蜜月| 亚洲久久久国产精品| 黑人欧美特级aaaaaa片| 日韩欧美一区视频在线观看| 亚洲国产精品专区欧美| 久久精品夜色国产| .国产精品久久| 伦理电影免费视频| 国产在线一区二区三区精| 午夜日本视频在线| 老司机影院成人| 国产成人精品久久久久久| 久久国内精品自在自线图片| 久久99一区二区三区| 老司机影院成人| 三级国产精品欧美在线观看| 简卡轻食公司| 亚洲精品自拍成人| 国产亚洲精品久久久com| 日韩熟女老妇一区二区性免费视频| 91久久精品国产一区二区三区| 亚洲精品乱码久久久v下载方式| 亚洲成人av在线免费| 亚洲综合色惰| 三级国产精品片| 国产av码专区亚洲av| 黄色怎么调成土黄色| 一个人看视频在线观看www免费| 一边亲一边摸免费视频| 久久这里有精品视频免费| 国产精品欧美亚洲77777| 一区二区日韩欧美中文字幕 | 国产男人的电影天堂91| 日韩欧美一区视频在线观看| 搡老乐熟女国产| 两个人免费观看高清视频| 亚洲av成人精品一二三区| 久久久久国产精品人妻一区二区| 国产成人精品在线电影| 多毛熟女@视频| 涩涩av久久男人的天堂| 免费播放大片免费观看视频在线观看| 日日啪夜夜爽| 777米奇影视久久| 成人综合一区亚洲| a级毛片免费高清观看在线播放| av在线老鸭窝| 国产精品秋霞免费鲁丝片| 久久久精品区二区三区| 成人毛片a级毛片在线播放| 色视频在线一区二区三区| 亚洲丝袜综合中文字幕| 最近最新中文字幕免费大全7| 国产精品熟女久久久久浪| 国产成人精品福利久久| 大码成人一级视频| 中文字幕制服av| 午夜福利网站1000一区二区三区| 热re99久久国产66热| 日本黄色日本黄色录像| 国产免费一区二区三区四区乱码| 久久人妻熟女aⅴ| a级毛片黄视频| 日本黄大片高清| 欧美日韩av久久| 婷婷色麻豆天堂久久| 在线观看www视频免费| 欧美人与性动交α欧美精品济南到 | 亚洲色图 男人天堂 中文字幕 | 精品午夜福利在线看| 国产成人精品一,二区| 国产亚洲欧美精品永久| 久久久精品免费免费高清| 久久久国产一区二区| 少妇被粗大的猛进出69影院 | 亚洲成人手机| 永久网站在线| 麻豆乱淫一区二区| 日韩伦理黄色片| 老熟女久久久| 久久久国产一区二区| 另类亚洲欧美激情| 啦啦啦中文免费视频观看日本| 国产白丝娇喘喷水9色精品| 99九九在线精品视频| 欧美xxⅹ黑人| 麻豆乱淫一区二区| 99热国产这里只有精品6| 精品亚洲乱码少妇综合久久| xxx大片免费视频| 精品人妻在线不人妻| 亚洲情色 制服丝袜| 久久精品国产亚洲av涩爱| 久久久久人妻精品一区果冻| 中文字幕久久专区| 成人国产麻豆网| 麻豆乱淫一区二区| 精品熟女少妇av免费看| 一区二区三区四区激情视频| 午夜免费男女啪啪视频观看| 午夜久久久在线观看| 在线观看一区二区三区激情| 看非洲黑人一级黄片| 夜夜骑夜夜射夜夜干| 欧美激情极品国产一区二区三区 | 亚洲,一卡二卡三卡| 久久久久国产网址| 国产精品久久久久久av不卡| 免费观看在线日韩| 久久精品夜色国产| 少妇被粗大猛烈的视频| 国产乱人偷精品视频| 大片电影免费在线观看免费| 国产精品不卡视频一区二区| 夫妻午夜视频| 日韩一区二区视频免费看| 制服人妻中文乱码| 色94色欧美一区二区| 亚洲国产最新在线播放| 久久久久久久久久久久大奶| 男人爽女人下面视频在线观看| 在线天堂最新版资源| 精品久久久精品久久久| 中文字幕久久专区| 亚洲伊人久久精品综合| 国产精品偷伦视频观看了| 特大巨黑吊av在线直播| 欧美日韩在线观看h| a级毛色黄片| 制服诱惑二区| 精品亚洲乱码少妇综合久久| 色婷婷av一区二区三区视频| 乱人伦中国视频| 亚洲av男天堂| kizo精华| 精品国产国语对白av| 在线看a的网站| 天堂俺去俺来也www色官网| 免费播放大片免费观看视频在线观看| 久久国产精品男人的天堂亚洲 | 性色avwww在线观看| 日韩av不卡免费在线播放| 18禁动态无遮挡网站| 桃花免费在线播放| 另类亚洲欧美激情| 亚洲成人手机| 久久久亚洲精品成人影院| 嘟嘟电影网在线观看| 国产成人91sexporn| 久热这里只有精品99| 国精品久久久久久国模美| 制服丝袜香蕉在线| 啦啦啦啦在线视频资源| 成人国产av品久久久| a级毛片免费高清观看在线播放| 欧美激情国产日韩精品一区| 美女cb高潮喷水在线观看| 汤姆久久久久久久影院中文字幕| 乱人伦中国视频| 欧美另类一区| 久久女婷五月综合色啪小说| 狂野欧美激情性bbbbbb| 亚洲美女视频黄频| 80岁老熟妇乱子伦牲交| av卡一久久| 寂寞人妻少妇视频99o| 啦啦啦啦在线视频资源| 爱豆传媒免费全集在线观看| 欧美+日韩+精品| 最近手机中文字幕大全| 免费高清在线观看视频在线观看| 精品人妻一区二区三区麻豆| 青春草亚洲视频在线观看| 久久国产精品大桥未久av| 97精品久久久久久久久久精品| 成人国产麻豆网| 99久久精品一区二区三区| 久久97久久精品| 三级国产精品欧美在线观看| 如何舔出高潮| 日韩av免费高清视频| 久久久久久久久久成人| 一本一本综合久久| 69精品国产乱码久久久| 日韩电影二区| 卡戴珊不雅视频在线播放| 美女大奶头黄色视频| 久久久久久人妻| 亚洲综合色网址| 亚洲精品一二三| 99九九在线精品视频| 最新的欧美精品一区二区| 日日啪夜夜爽| 9色porny在线观看| 日本-黄色视频高清免费观看| 超碰97精品在线观看| 97超视频在线观看视频| 精品人妻熟女毛片av久久网站| 亚洲精品国产av成人精品| 中文欧美无线码| 精品亚洲成国产av| 人人妻人人添人人爽欧美一区卜| 欧美xxⅹ黑人| 亚洲av不卡在线观看| 久久久久久久久大av| 午夜影院在线不卡| 黑人猛操日本美女一级片| 日本欧美国产在线视频| 欧美激情 高清一区二区三区| 一级毛片电影观看| 另类亚洲欧美激情| 中文字幕最新亚洲高清| 久久av网站| 久久97久久精品| 黑丝袜美女国产一区| 99久久精品国产国产毛片| 人妻夜夜爽99麻豆av| 国产精品久久久久久精品电影小说| 亚洲精品乱码久久久久久按摩| 97精品久久久久久久久久精品| 免费高清在线观看视频在线观看| 热re99久久精品国产66热6| av播播在线观看一区| 秋霞在线观看毛片| 亚洲精品,欧美精品| 亚洲精品456在线播放app| 观看av在线不卡| 极品少妇高潮喷水抽搐| 成人漫画全彩无遮挡| 国产成人精品久久久久久| 中文字幕亚洲精品专区| 男女国产视频网站| 高清午夜精品一区二区三区| 另类精品久久| 国产亚洲欧美精品永久| 人人妻人人澡人人爽人人夜夜| 国产精品国产三级国产专区5o| 亚州av有码| 精品酒店卫生间| 亚洲精品成人av观看孕妇| 一本—道久久a久久精品蜜桃钙片| 久久精品国产亚洲av涩爱| 久久久久精品性色| 国产伦精品一区二区三区视频9| 91精品一卡2卡3卡4卡| 亚洲av免费高清在线观看| 国产极品粉嫩免费观看在线 | 少妇被粗大的猛进出69影院 | 麻豆乱淫一区二区| 各种免费的搞黄视频| 蜜桃在线观看..| 视频中文字幕在线观看| 久久女婷五月综合色啪小说| 夜夜骑夜夜射夜夜干| 国产女主播在线喷水免费视频网站| 51国产日韩欧美| 日产精品乱码卡一卡2卡三| 亚洲人成网站在线播| 亚洲国产av新网站| 亚洲欧美成人精品一区二区| 国产精品久久久久成人av| 亚洲av中文av极速乱| 亚洲国产成人一精品久久久| 国产成人免费观看mmmm| 国产色婷婷99| 黄色怎么调成土黄色| 91午夜精品亚洲一区二区三区| videos熟女内射| 高清不卡的av网站| 蜜臀久久99精品久久宅男| 男人添女人高潮全过程视频| 国模一区二区三区四区视频| 欧美日韩在线观看h|