999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

Virtual sensing method for monitoring vibration of continuously variable configuration structures using long short-term memory networks

2020-02-22 10:51:22ZhenjingYUELiLIUTengLONGYunchenMA
CHINESE JOURNAL OF AERONAUTICS 2020年1期

Zhenjing YUE, Li LIU,b,*, Teng LONG,b, Yunchen MA

a School of Aerospace Engineering, Beijing Institute of Technology, Beijing 100081, China

b Key Laboratory of Dynamics and Control of Flight Vehicle, Ministry of Education, Beijing 100081, China

KEYWORDS

Abstract Vibration monitoring by virtual sensing methods has been well developed for linear timeinvariant structures with limited sensors. However, few methods are proposed for Time-Varying(TV) structures which are inevitable in aerospace engineering. The core of vibration monitoring for TV structures is to describe the TV structural dynamic characteristics with accuracy and efficiency. This paper propose a new method using the Long Short-Term Memory (LSTM) networks for Continuously Variable Configuration Structures(CVCSs),which is an important subclass of TV structures.The configuration parameters are used to represent the time-varying dynamic characteristics by the‘‘freezing”method.The relationship between TV dynamic characteristics and vibration responses is established by LSTM,and can be generalized to estimate the responses with unknown TV processes benefiting from the time translation invariance of LSTM.A numerical example and a liquid-filled pipe experiment are used to test the performance of the proposed method. The results demonstrate that the proposed method can accurately estimate the unmeasured responses for CVCSs to reveal the actual characteristics in time-domain and modal-domain.Besides,the average one-step estimation time of responses is less than the sampling interval.Thus,the proposed method is promising to on-line estimate the important responses of TV structures.

1. Introduction

Vibration monitoring is a fundamental issue of structural health monitoring, fatigue assessment, flutter suppression and so on.1–5Available techniques for physical measuring vibration responses are limited in practical aerospace engineering structures since vibration responses of certain Degrees of Freedom(DOFs)via physical measurements are costly or even inaccessible. One appealing solution is to use virtual sensing method.This method uses accessible online measurements and mathematical models to estimate a quantity which is difficult or expensive to measure in real time.6–9It is widely used in chemical process fields for monitoring and controlling industrial production,and is a suitable way to estimate those vibration responses of unmeasured DOFs.

As for the mathematical model for the structural dynamics,there are three objects: loads, structural systems, and responses. When one of them is required, we need to know two of them. Response analysis requires the loads and the characteristics of the given structural system; loads can be identified with known structural systems and responses; the identification for the given structural system needs loads and responses information. As for the virtual sensing task, it is practically inaccessible to obtain ambient excitation. If the responses of certain unmeasured DOFs can be estimated by the responses of measured DOFs, the structural system information must be known. How to use the structural system information is the key for virtual sensing methods.

Virtual sensing methods have been applied to monitoring the vibration of linear Time Invariant (TI) structures.Although these works usually claim that they are outputonly, they all use the structural information directly or indirectly.Those methods can be divided into two different classes:namely model-based and data-based. There are two kinds of methods in the model-based virtual sensing methods:one using mode shape information and the other using system matrices via Kalman filter. Generally, the idea of mode shape based methods is that different DOFs in physical domain correspond to the same modal coordinates for each mode. The modal coordinates can be obtained by the responses of measured DOFs and then the responses of unmeasured DOFs in physical domain can be transformed from modal domain via mode shape. Relevant DOFs can be translational,10rotational11or even strain,12,13etc. The structural information can be global or partial.14The main notion of the Kalman filter used in the virtual sensing is the estimation of the states(displacement and velocity) of a partially observable linear system.15,16This method can take full account of the uncertainty about measurements and structural systems. The data-based methods establish the relationship between measured responses and estimated responses based on the off-line dataset that can be built by the full measurements before the service for the structure.17,18This kind of methods do not require the structure system information directly, because the constant structure system information are included in the dataset for linear time invariant structures.

The machine learning has been well developed and widely used in engineering.19–21Further, Deep learning has achieved great success in last ten years.22As an important model in deep learning, the Recurrent Neural Network (RNN) is a powerful tool for the task involving sequential inputs. The sequential inputs transform into hidden layer states by RNN, which preserve the historical information of the whole sequence. However, the challenge of long-term dependencies for RNN due to the gradients vanish has limited the application in practiced long sequences.23Long Short-Term Memory (LSTM) networks introduce self-loops and can consider the influence of long-term input sequences.24LSTM has proved to be effective in sequential problems such as machine translation,25dialogue systems,26video captioning.27

The aforementioned vibration monitoring works by virtual sensing methods can only be applied to the time-invariant structures. However, time-varying characteristics of structures are becoming increasingly inevitable in aerospace engineering.28–32The Continuously Variable Configuration Structures(CVCSs)characterized by configuration parameters is a typical subclass of Time-Varying (TV) structures,33such as morphing aircrafts,30aircrafts in aerial refueling,32launch vehicle with fuel consumption.34Ignoring the time-varying characteristics may lead to reducing accuracy and increasing redundancy.The core issue of vibration monitoring for TV structures is to describe the TV structural information with accuracy and efficiency.

In this paper,we investigate vibration monitoring by virtual sensing method using LSTM for CVCSs.There are three main contributions: (A) the challenges of model-based methods in virtual sensing tasks for the time-varying structures are analyzed; (B) aiming at the characteristics of typical aerospace time-varying structures, a data-driven virtual sensing method is proposed for CVCSs; (C) LSTM is adopted to establish mapping from measured responses and configuration parameters to unmeasured responses, due to its time translation invariance which can handle unknown TV processes under operational conditions. In the next section, the description of CVCSs and the‘‘freezing”method are introduced.In Section 3,the response relationship between measured DOFs and unmeasured DOFs is derived to reveal the challenging of virtual sensing methods for TV structures. In Section 4, the LSTM is briefly introduced and the virtual sensing method for CVCSs is proposed.In Section 5 and Section 6,a numerical example and a liquid-filled pipe experimental example are used to valid the virtual sensing method.Finally,we close the paper with concluding comments and suggestions for future work in Section 7.

2. Continuously variable configuration structures

As an important subclass of TV structures, the description of CVCSs will be briefly introduced, and more details can be found in Ref.33

The governing equation of a continuous-time, lumped parameter, linear TV structure S is

whereu(t)is the displacement vector;F(t)is the external excitation;M(t),C(t)andK(t)are the time-varying mass,damping and stiffness matrices, respectively; [t0,tf] is the considerable time interval.

A sequence TI structures SFTcan be obtained through‘‘freezing” the system matrices at each time instant

where the set SFTis the ‘‘frozen-time” representation of the correspondi ng TV structure; SFT(t′) is the element of the SFT, and stands for a special TI structure at timet′.

For a CVCS ˉS, a configuration parameter vectorp(t) can be defined to determine the state of the structure at each time.The governing equation is

wherePis the set of structural configurations during [t0,tf].

That means TV dynamic characteristics of CVCSs can be described by the ‘‘frozen-configuration” using ‘‘freezing”method.

3. Description of virtual sensing task for CVCSs

The structural dynamic characteristics, such as the modal parameters, can be represented by the structural ‘‘frozen-con figuration” for a continuously variable configuration, linear TV structure ˉS. Therefore, the structural vibration responses can be expressed as

wherep′is the structural ‘‘frozen-configuration” at timet, as defined inEq. (5); φi(p′) andqi(t) are the mode shape vector and modal coordinate of theith mode corresponding at timet; Φ(p′) andq(t) are the mode shape matrix and vector of modal coordinate at timet;Nis the total number of DOFs of the linear system.

By partitioning the responsesu(t) and the mode shape matrix Φ(p′)intothe measured DOFs and unmeasured DOFs,we have

where (·)+is the pseudo-inversion.

The responses of unmeasured DOFsuu(t) can be obtained by the responses of measured DOFsum(t) and the TV mode shape Φ (p′). There are some difficulties in the methods based on mode shape information for the virtual sensing task. On the one hand, there is an urgent need for efficient and robust identification methods or a series of accurate ‘‘frozen-config uration” models to obtain the TV mode shapes; on the other hand the calculation of the pseudo-inversion needs a numerically stable and efficient method. There are similar difficulties in the methods using the system matrices via Kalman filter.To avoid repetition, they are not described in this article. The model-based methods are challenging to achieve the vibration virtual sensing task for CVCSs.

4. Proposed virtual sensing method

In this paper, a data-based method is proposed to achieve the virtual sensing task for CVCSs via LSTM networks,which are an important kind of RNNs.The LSTM can extract the information from the time series data and can represent long term dependencies. Those characteristics are crucial for handling the responses of TV structures. In this section, the recurrent formulation and the basic conception of LSTM are briefly introduced.23,24Then the virtual sensing method using LSTM networks is subsequently proposed.

4.1. Recurrent formulation and time translation invariance

RNNs are a kind of neural networks to deal with sequential problems by the feedback mechanism. They extract the information from the sequence datax1,x2,···,xnand store it in the hidden layer outputh. The recurrent formulation can be expressed as

where thextis the input at timet;θ are the fixed parameters of a given RNN. Then

Therefore,as illustrated in Fig.1,the hidden layer outputhtcomprises the information in the current input and the past sequence. The parameters θ, sharing across the time, decide the nonlinear map from the input sequence to the hidden layer output.This mechanism makes the RNNs have the time translation invariance to process sequences with different lengths.In addition, the sharing parameters is crucial for handling unknown TV processes because the significant information may occur at different time (or different position in the sequence).

4.2. Review of LSTM architecture

Fig.1 Extract information by sharing parameters.

wherextandhtare the current time step input vector and current hidden layer output, respectively;UF,UI,UOandUare respectively the input weights for the forget gate, input gate,output gate and LSTM cell;WF,WI,WOandWare respectively the recurrent weights for the forget gate,input gate,output gate and LSTM cell;bF,bI,bOandbare respectively the biases for the forget gate, input gate, output gate and LSTM cell; σ(·) is the sigmoid function.

The output of the hidden layer is

4.3. Virtual sensing method using LSTM

To solve the core issues of describing TV dynamic characteristics of CVCS, this work adopts the configuration parameters to characterize the structural information, and the LSTMbased network is used to establish the relationship among the structural configuration, responses of measured DOFs and responses of unmeasured DOFs.

The architecture of proposed virtual sensing method for CVCSs is composed of three parts: Input embedding layer,LSTM layers, Full connecting layers, as shown in Fig.2.

The input data contains the response signal data and the configuration data both obtained by measurement at the timet. These data need to be transformed into the unified input space. A linear transformation is used to achieve this goal in the input embedding layer by following equation

In the LSTM layers, the features are extracted in temporal dimension(relationship among the past time steps)and spatial dimension (relationship among the sensor networks). Each LSTM layer can be briefly expressed as Eq. (19), and more details can be found in Section 4.2.

The full connection layers, a full connection feedforward network, further abstract the feature obtained by LSTM,and estimate the responses of unmeasured DOFs.Each nonlinear layer can be defined by following equation

Different from classification problem, the estimation of responses is a typical regression problem. The last layer of the full connection part is a linear layer to obtain the estimate of responses

whereWOandbOare the weight matrix and bias vector of the linear full connection layer; ^ytis the estimate of responses at timet.

The cost function of virtual sensing model is defined by following equation

Fig.2 Architecture of proposed virtual sensing model.

Mean Absolute Error (MAE) and Mean Squared Error(MSE)have been adopted to evaluate the accuracy of the proposed virtual sensing method in time-domain.

whereTNis the numbers of time steps;nis the numbers of unmeasured DOFs;yis the label value; ^yis the estimated value.

5. Numerical example

To valid the vibration virtual sensing method for CVCSs, a numerical example with TV dynamic characteristics under the random vibration is considered in this section.

5.1. Description of numerical example

The numerical example is a 7 DOFs mass-damping-spring TV system,shown in Fig.3,whose mass and damping characteristics are invariant but the stiffness characteristic is time-variant.The stiffnesski(t) of theith spring changes with time as:

wherepi(t)is the configuration parameter in this example.The time constantsTjand phase constants φj(j=1,2,...,16),defined in Eq. (26), are randomly generated by the Latin Hypercube Designs (LHDs) sampling method,37andTj∈[80,200],φj∈[0,2π]. The rest parameters of system,defined in Eq. (25) and Fig.3, are listed in Table 1.

Gaussian white noise excitations are applied to 1 DOF.The excitation frequency lies within the 0–8 Hz,because the highest natural frequency is lower than 4 Hz. The dynamic responses of the system are obtained by the numerical integration method, Runge-Kutta method, with time step Δt=0.015625 s. The considerable time interval is 0–500 s.

In order to roughly illustrate the TV characteristics of the system, the system’s partial configuration parameters and the‘‘frozen” natural frequencies of two cases with differentTjand φjare depicted in Fig.4. The seven curves from the bottom to the top in Figs. 4(a) and 4(c) represent the first to the seventh natural frequencies of system.On one hand,each natural frequency of system varies obviously during 0–500 s, for instance, the first natural frequency within the range 0.3–0.6 Hz and the seventh natural frequency within the range 2.3–3.6 Hz. They have changed more than 50%. On another hand, by adjusting the time constantsTjand phase constants φj, configuration parameters with different changing process(see Figs. 4(b) and 4(d)) can be obtained, which results in different system TV dynamic characteristics (see Figs. 4 (a)and 4(c)) such as the difference between case A and case B.

5.2. Virtual sensing of mass-damping-spring TV system

Fig.3 7-DOF varying stiffness system.

Table 1 Property parameters of the TV mass-damping-spring system.

Fig.4 ‘‘Frozen” characteristics of system.

In this example,we use the displacement responses of 1,3,5,7 DOFs and the configuration parameters to estimate the displacement responses of 2, 4, 6 DOFs. To obtain the accurate virtual sensing model for different TV characteristics of configuration parameters,50 sets of time constantsTjand phase constants φjthat control the changing process of the configuration parameters are generated via LHDs.The LHDs can guarantee samples with good space-filling and projective properties in considered parameter space.In addition,the excitation is independent and identically distributed in each of sets. 50 sets of generated response data are randomly divided into train set(40 sets) and test set(10 sets). The response data of both train set and test set are normalized by the mean and standard deviation of train set responses. In addition, the configuration parameters of both train set and test set are normalized by the maximum of train set configuration parameters.When the virtual sensing method is used to estimate the responses, the normalized responses will be calculated firstly,and then the physical domain responses are obtained via the inverse transformation by the same normalized parameters.

The Input embedding layer has 75 linear units. The LSTM part of virtual sensing model consists of two LSTM layers and each layer has 75 LSTM units. The Full connection part has two full connection nonlinear layers with 50 units in each layer. The generalization factor is 0.01.

The dynamic estimate of responses compares well to the numerical model. Test errors of all ten test cases are given in Table 2.As an example,the 6 DOF estimated result of test case 4 is shown in Fig.5.The main error exists in the initial stage as shown in Fig.5(b).As most recursive algorithms,the principal reason is the initial state of model is not appropriate.However,this error rapidly decays over time. Moreover,the accuracy of the virtual sensing method without using configuration parameters is also given in Table 2. The results demonstrate that the use of configuration parameters can significantly improve the performance for virtual sensing task.That means virtual sensing methods for TV structures should be able to exactly describe the TV dynamic characteristics of TV structures,and the proposed method can achieve this goal succinctly and accurately.

6. Experiment example

To further verify the validity and accuracy of the proposed vibration virtual sensing method for CVCSs, the liquid-filled pipe is excited during the release of liquid in order to collect vibration signals of the real TV structure, and the estimated responses are compared with the known responses. In addition,the estimated results are applied to structural modal analysis to improve the spatial resolution of mode shapes.

6.1. Description of experiment

The experiment structure is a liquid-filled pipe with a 2.4 m height and a 0.05 m radius as shown in Figs.6 and 7.The pipe was made of stainless steel with 1 mm thickness.The pipe was suspended using four polyester ropes. The structural dynamic characteristics are TV with the continuous release of the water in the pipe. The remote controller can control the rudder to determine the start of the water release. The liquid level mea-suring system consists of a liquid surface float and a laser range finder to obtain the liquid level in experiment. In addition, by changing the inner diameter of the outlet,different liquid level changing processes are realized,which obtains the different TV dynamic characteristics of the experimental structure.

Table 2 Estimated errors in test set.

Fig.5 Results comparison and their partial enlarged details of 6 DOFs.

Fig.6 Schematic diagram of experiment system.

Ten piezoelectric accelerometers and an impedance head measure the acceleration responses of the pipe at eleven uniformly distributed positions along the axial direction of the pipe, although only three accelerometers are shown in Fig.6.The impedance head is located at the third position from the bottom to the top, which is the excitation location of the shaker.The shaker and the impedance head use the threaded connection and the impedance head is glued on the pipe.However,the whole structure moves upwards with the decrease of mass in the process of water release,which adversely affects the shaker. In order to reduce the rigid body displacement in vertical direction of the pipe,adding two mass blocks(each 5 kg)to the pipe to reduce the relative variation of mass, and the suspension rope of the pipe adopts the polyester rope with smaller elasticity.A LMS SCADAS III system collects the acceleration response signals, liquid level signals and control exciting force of the shaker.

The liquid-filled pipe is a CVCS due to water release continuously and its dynamic characteristics can be decided by the liquid level. During the experiment, a Gauss white noise excitation within the 0–256 Hz is applied to the pipe.The sampling frequency is 4096 Hz, and the record length is 32 s. The acceleration responses and liquid level are measured.

Fig.7 Liquid-filled pipe experiment.

6.2. Virtual sensing of experimental system

In this experimental example, we use six physical quantities,including the transverse acceleration responses at position 1,3,5,9,11 and the liquid level,to estimate the transverse acceleration responses at position 4,6,7,and 8.The different inner diameters of the outlet can get different water release processes. Fifty experiments been conducted with different inner diameter of the outlet (14.1 mm, 14.4 mm, 14.7 mm, ...,28.5 mm,28.8 mm).Take ten groups of data evenly as the test set (15.3 mm, 16.8 mm, ..., 28.8 mm), the other experiment data are used as train set. The collected data are filtered by a FIR filter with the passband of 256 Hz. Then the normalization setup is same as the Section 4.The transverse acceleration response data of both train set and test set are normalized by the mean and standard deviation of train set responses. In addition, the liquid level data of both train set and test set are normalized by the maximum of the liquid level in train set. When the virtual sensing method is used to estimate the responses, the normalized responses are calculated firstly,and then the physical domain responses are obtained via the inverse transformation by the same normalized parameters.The transverse acceleration response at position 6 and liquid level are shown in Fig.8. In the case with the inner diameter of the outletdo=16.8 mm, start the release of water after starting collecting data about 2 s. Then the liquid level dropped to Position 6 whent≈12 s. Finally, the water was put out at about 31.5 s. However, when the inner diameter of the outlet is increased to the 27.3 mm,the water release process is more rapid.The release began at aboutt≈2 s.Then liquid level dropped to position 6 at aboutt≈6 s, and the end of release occurs at aboutt≈14 s.The pipe is a TI structure in the final 18 s. Fig.8 shows clearly the influence of inner diameter of the outlet on the TV characteristics of the structure.

The Input embedding layer has 75 linear units. The LSTM part of virtual sensing model consists of two LSTM layers and each layer has 75 LSTM units. The Full connection part has two full connection nonlinear layers with 50 units in each layer. The generalization factor is 0.01.

The accuracy of estimated results is favorable. Test errors of all ten test cases are given in Table 3. The measured acceleration response,estimated response and their partial enlarged details of two examples can be seen in Fig.9 fordo=16.8 mm anddo=27.3 mm test cases.They are the maximal MAE and MSE test case,respectively.The maximal error always appears at the start release of water(see Figs.9(c)and 9(d)).The main reason for this error is the abrupt change of the structure when the remote controller controls the rudder to start releasing water.This error decreases with the virtual sensing model handling the data in new stage (see Figs. 9(e) and 9(f)).

Table 3 Estimated errors of test set in experiment example.

The virtual sensing method has promising for on-line estimation because there are only algebraic operation with sufficient computational efficiency in estimation. The sampling interval of the responses and liquid level in experiment is 0.244 ms(4096 Hz),whereas the average estimation time using the virtual sensing method is not more than 0.05 ms with CPU i7-6700 (see Table 3).

6.3. Application in modal analysis

As an application,the proposed virtual sensing method is used to augment data in modal analysis for improving the spatial resolution of mode shapes. The accurate mode shape is beneficial for active vibration control,structural health monitoring and so on.However,it is difficult to arrange enough sensors in actual structures due to plenty of practice limitations. This problem can be settled via estimating the responses at interesting positions using the proposed virtual sensing method. The virtual sensing method can not only obtain accurate timedomain estimated responses, but also the same structural modal characteristics.

To demonstrate the validity of the virtual sensing method,we consider three identification situations with three different response data sets:

(1) All real measurement responses at position 1, 3, 4, 5, 6,7, 8, 9, 11.

(2) All real measurement responses at position 1,3,5,9,11.

(3) Real measurement responses at position 1, 3, 5, 9, 11 and estimated responses at position 4, 6, 7, 8.

A Functional Series Vector Time-Dependent Autoregressive Moving Average (FS-VTARMA) method38is employed for the identification of the liquid-fill pipe. The responses of test casedo=16.8 mm during 3–15 s are used and the responses are re-sampled to the 512 Hz.

Fig.9 Acceleration responses and their partial enlarged details at Position 6 for test case do=16.8 mm and do=27.3 mm.

Fig.10 Natural frequency (the 1st row) and mode shape (the 2nd and 3rd row) at t=6 s of Situation A, B and C.

Fig.10 shows the natural frequency TV process and the mode shapes att=6 s for three identification situations.Comparing the identification results of the three response situations,the first conclusion is that the responses obtained by the virtual sensing method are consistent with those acquired with physical measurement in terms of modal characteristics. The identified frequencies using the virtual sensing method (situation C) are almost identical to those measured with physical sensors (situation A and B). Furthermore, the spatial resolution of mode shape is observably improved via the virtual sensing method.

7. Conclusions

The proposed virtual sensing method for CVCSs employs data-based method for estimating the vibration responses of unmeasured DOFs using LSTM. The key problem, that is how to describe the TV structural dynamic characteristics,has been solved via using configuration parameters and establishing the relationship between TV dynamic characteristics and vibration responses by LSTM.The time translation invariance of LSTM allows the proposed method to estimate responses in unknown TV processes. The proposed method is verified by a 7 DOFs numerical example and a liquid-filled pipe experimental example. This method proves to render accurate enough results in unknown TV processes.The principal error occurs in the initial stage or after abrupt changes,and can rapidly decays over time. In addition, the results demonstrate that this method is reasonably efficient due to the finite algebraic calculation without iteration, and is promising to estimate responses on-line.It should be noted that this method works for those CVCSs whose configuration parameters change within a finite range.For future work,the optimization of physical sensors placement for improving the estimation accuracy under the practical complex constraints should be studied.

主站蜘蛛池模板: 久久久久久久97| 久青草免费在线视频| 青草视频网站在线观看| 日韩欧美在线观看| 国产拍在线| 日韩大乳视频中文字幕| 99偷拍视频精品一区二区| 国产内射一区亚洲| 91网在线| 亚洲视频黄| 全色黄大色大片免费久久老太| 美女免费黄网站| 黄色网在线| 91视频青青草| 日本人又色又爽的视频| 国产成人久久777777| 久久黄色影院| 8090午夜无码专区| 青青久久91| 亚洲色成人www在线观看| 精品乱码久久久久久久| 国产白浆视频| 天天综合天天综合| 精品无码视频在线观看| 国产亚洲精| 网友自拍视频精品区| 国产免费好大好硬视频| 亚洲高清资源| 欧美日韩精品一区二区在线线| 国产色偷丝袜婷婷无码麻豆制服| 9啪在线视频| 奇米精品一区二区三区在线观看| 色婷婷亚洲综合五月| 无码一区中文字幕| 色综合中文综合网| 欧美日韩国产成人在线观看| 国产探花在线视频| 国产欧美日韩一区二区视频在线| 久久亚洲欧美综合| 播五月综合| 国产在线拍偷自揄拍精品 | 在线无码九区| 国产成人免费手机在线观看视频| 国产一区二区人大臿蕉香蕉| 五月天香蕉视频国产亚| 四虎影视国产精品| 亚洲综合经典在线一区二区| 人人91人人澡人人妻人人爽| 国产成人综合日韩精品无码首页| 国产精品美人久久久久久AV| 99热这里只有精品久久免费| 2021天堂在线亚洲精品专区| www.日韩三级| 国产另类乱子伦精品免费女| 国产真实二区一区在线亚洲| 亚洲精品无码专区在线观看| 国内精品手机在线观看视频| www中文字幕在线观看| 欧美日韩第二页| 亚洲男人的天堂久久香蕉网| 有专无码视频| 亚洲看片网| 亚洲美女高潮久久久久久久| 国产爽爽视频| 99国产精品一区二区| 尤物精品国产福利网站| 日韩欧美高清视频| 亚洲成a人在线播放www| 人妻丰满熟妇αv无码| 免费在线看黄网址| 无码综合天天久久综合网| 久久人体视频| 国产91无码福利在线| 欧美一区二区丝袜高跟鞋| 国产精品欧美亚洲韩国日本不卡| 天堂成人在线视频| 99这里只有精品6| 91青青视频| Jizz国产色系免费| 免费jjzz在在线播放国产| 欲色天天综合网| 无码精品福利一区二区三区|