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

    NUMERICAL STUDIES OF INTERNAL SOLITARY WAVE GENERATION AND EVOLUTION BY GRAVITY COLLAPSE*

    2012-08-22 08:31:49LINZhenhuaSONGJinbao

    LIN Zhen-hua, SONG Jin-bao

    Institute of Oceanology, Chinese Academy of Sciences, Qingdao 266071, China

    Key Laboratory of Ocean Circulation and Waves, Chinese Academy of Sciences, Qingdao 266071, China, E-mail: linzhenhuaouc@gmail.com

    (Received December 22, 2011, Revised April 21, 2012)

    NUMERICAL STUDIES OF INTERNAL SOLITARY WAVE GENERATION AND EVOLUTION BY GRAVITY COLLAPSE*

    LIN Zhen-hua, SONG Jin-bao

    Institute of Oceanology, Chinese Academy of Sciences, Qingdao 266071, China

    Key Laboratory of Ocean Circulation and Waves, Chinese Academy of Sciences, Qingdao 266071, China, E-mail: linzhenhuaouc@gmail.com

    (Received December 22, 2011, Revised April 21, 2012)

    In this study, an analysis on the internal wave generation via the gravity collapse mechanism is carried out based on the theoretical formulation and the numerical simulation. With the linear theoretical model, a rectangle shape wave is generated and propagates back and forth in the domain, while a two-dimensional non-hydrostatic numerical model could reproduce all the observed phenomena in the laboratory experiments conducted by Chen et al. (2007), and the related process realistically. The model results further provide more quantitative information in the whole domain, thus allowing an in depth understanding of the corresponding internal solitary wave generation and propagation. It is shown that the initial type of the internal wave is determined by the relative height between the perturbation and the environmental density interface, while the final wave type is related to the relative height of the upper and lower layers of the environmental fluid. The shape of the internal wave generated is consistent with that predicted by the KdV and EKdV theories if its amplitude is small, as the amplitude becomes larger, the performance of the EKdV becomes better after the wave adjusts itself to the ambient stratification and reaches an equilibrium state between the nonlinear and dispersion effects. The evolution of the mechanical energy is also analyzed.

    gravity collapse, internal wave, energy analysis

    Introduction

    The internal waves are the gravity waves that occur within, rather than on the surface of, a stratified fluid. The density differences within the fluid are usually much smaller than the density of the fluid itself, thus the internal waves typically have much lower frequencies and larger amplitudes than their surface counterparts. They are ubiquitous in the ocean and the atmosphere, and may induce oscillations of the density interface, to alter the magnitude and the direction of acoustic transmissions; may also bring nutrients up from the deeper ocean into the euphotic layer and contributes to the primary production. The internal wave breaking is a major driver in the local mixing process.

    The stratified fluid and the perturbation are nece-ssary to excite the Internal Solitary Waves (ISW). One common generation mechanism is the interaction of the tidal flow and the irregular bottom topography of the ocean floor[1-4]. Other generation mechanisms were also explored, for example, through river plumes[5], and the nonlinear steepening from an initial tilted density interface[6]. The ISW is of a depression type if the thickness of the lower layer is larger than that of the upper layer, after it passes the turning point where the thickness of the deeper layer is equal to that of the upper layer, the internal wave of depression gradually turns to one of an elevation type, finally, the ISW breaking occurs at coastal boundaries and contributes to the local mixing.

    Various methods, e.g., the theoretical formulation[7], the in situ observation[8], the satellite image[9], the numerical model[10]and the laboratory experiment[11-15], were applied to study the ISW generation, propagation, interaction and breaking processes. Since it is often difficult to thoroughly investigate the characteristics of ISWs in field observations, and certain assumptions such as weak nonlinearity are commonly used in the theoretical formulation, some alternativeapproaches combined with laboratory tests to study the ISWs were developed[16,17]. Moreover, the laboratory tests are superior to other methods because of low cost, and good repeatability and reliability of results. However, unlike in oceanic scenarios, where the ISWs are mostly generated by stratified flows over an irregular topography, the gravity collapse is one of the widely used methods in the laboratory experiments to excite the ISWs. This method is simpler than others, such as those by the oscillating boundary or the movable body, and, therefore, is commonly used in the laboratory. Although the ISWs propagation and transformation properties were much studied via this generation mechanism[16], the ISWs generation and propagation by gravity collapse remains an interesting topics of research.

    In the laboratory experiments conducted by Chen et al.[16], the ISWs are generated by gravity collapse in a two-layer stratified fluid system in a flume. The laboratory settings can be found in detail in their paper. The main conclusions are: (1) if the upper layer is thinner than the lower layer (i.e., the environmental condition of H1<H2), the ISWs of depression type may be generated, on the other hand, the ISWs of elevation type are generated when H1>H2, (2) the ISW in the wave flume maintains its soliton feature if the environment condition permits, otherwise, a train of wave oscillations occurs following the leading wave.

    The conclusions of Chen et al.[16]are very instructive for later laboratory experiments on the ISW generation via gravity collapse mechanism, however, there are two apparent shortcomings inherent in their studies. Firstly, the density interface displacement at the generation site is determined via video snapshots using a patch of dye transfused at the interface. As the motion starts, it would not be easy to identify the interface because the dye is very diffusive and discontinuous. Secondly, there is no velocity field data measured in the experiment, which makes it impossible to discern the energy evolution characteristics of this process. The results presented in Chen et al.[16]are valuable from a qualitative point of view, but further studies are necessary to understand this ISW generation and propagation mechanism quantitatively. In this study, a high-resolution two-dimensional non-hydrostatic numerical model is used to fill this gap, which is based on the Reynolds Averaged Navier-Stokes equations. It should be stressed that the numerical method used in this study is different from that based on the theoretical simplified model[17-19], which generally contains certain assumptions, such as the inviscid fluid, weakly nonlinear and weakly dispersive etc., thus would provide more robust results.

    The goal of the present paper is to better understand the ISW generation, propagation and reflection by the gravity collapse mechanism. This paper is arranged as follows. In Section 1 the mathematical problem is formulated. The model used is introduced in Section 2 and the numerical results are analyzed in Section 3. Finally conclusions are drawn in Section 4. The current study may provide some food of thought for future laboratory experiments on how to arrange the fluid with different densities to get the desired ISW type.

    Fig.1 Schematic representation of the initial condition of the gravity collapse scenario with the main region in the left and the perturbation region in the right. The dashed line marks the initial interface distribution, and the dashdotted line marks the final equilibrium interface position

    1. Mathematical formulation

    In the present study, the effects of the earth’s rotation is neglected and the system is subject to only a single forcing event. The schematic representation of the system is shown in Fig.1. The flow is driven by the baroclinic pressure gradient that results from a density step. The interface displacement for a twolayer stratified fluid is governed by the linear non-dispersive wave equation under the rigid lid approximation and the Boussinesq approximation[6],

    where the horizontal coordinate axis x is fixed along the bottom boundary (positive rightward), the vertical axis z is positive upward, t is the time, (,)xtη is the vertical displacement of the density interface,is the linear long-wave speed,is the reduced gravity, and1H and2H are the depths of the upper and the lower layers, respectively.

    For the two superposed incompressible fluids with different densities, the no-flux boundary conditions apply on two ends, thus it is required that the horizontal velocity vanishes and the fluid interface remains perpendicular to the boundary[6],

    The initial conditions consist of a piecewise distribution of the interface displacement and no initial motion,

    where1f and2f are constants representing the interface heights at different positions,0x is the position of the gate. The general solution for the interface displacement is found in the following form,

    wherenk are the coefficients in the Fourier cosine series,

    The general analytical result shown above involves the infinite summation of trigonometrical functions, and certain procedures are required to simplify it to facilitate a direct understanding of what the solution represents. For simplicity, in the following we assume L/2<x0<L, then Eq.(4) can be simplified as

    The detail of the mathematical transformation is given in appendix. It should be noted that the above solution is only valid for 0<t<(L-x0)/c0. Similar tactics can be used to derive the alternative forms of Eq.(4) within other parts of the period 0<t<2L/c0. The solutions for one period are summarized as follows

    Now the solution (4) is reformulated into a group of piecewise spatial functions in different time ranges, allowing a direct and accurate identification of what kinds of motions they represent. Without loss of generality, we further assume that f1>f2, Figs.2(a)-2(h) shows the temporal evolution of the density interface displacement for a whole period. From the figure it can be seen that from the initial condition 2, two opposite directional interface motions with the same amplitude are generated as shown in Fig.2(b), in Fig.2(c) it is shown that, after the right limb reaches the right boundary, it is reflected and then moves leftward, a rectangle type wave is formed with a constant propagation speed0c, as inferred from Eqs.(7), the amplitude of the wave is equal to a half of the initial density interface depth difference across the step position0x. When the left limb reaches the left boundary, it is reflected and moves rightward as shown in Fig.2(d). Finally, the two branches meet at L-x0and then the inversion process continues. It is worth noting that from the initial state, as shown in Fig.2(a), to the steady wave propagation state, as shown in Fig.2(c), a half of the available potential energy is lost, which means that for this theoretical model, the kinetic energy of the system is equal to the Available Potential Energy (APE) when the wave is formed.

    Fig.2 Temporal evolution of the density interface obtained by the analytical solution

    The above solution shows the case for L/2<x0<L, a similar method can be used to obtain the solution for 0<x0<L/2. It can also readily be inferred that for x0=L/2, a standing wave solution is obtained, however, in the current work this scenario is not further analyzed, since the propagation waves are desired in almost all laboratory experiments.

    2. The numerical model

    In this section, a two-dimensional non-hydrostatic model is constructed to simulate the ISWs generated by the gravity collapse. Our main objective here is to use the laboratory observations made by Chen et al.[16]to validate our numerical results, and the model results in turn could provide more reliable information,such as the temporal and spatial distributions of the interface displacements, the horizontal and vertical velocity fields across the whole domain, for a better understanding of the ISW generation, propagation and reflection processes.

    The model amounts to solving the two-dimensional Reynolds Averaged Navier-Stokes (RANS) equations for an incompressible Newtonian fluid with the Boussinesq approximation,

    the tracer transport equation and the equation of state,

    where (,)xz are the horizontal and vertical coordinates and (,)uw are the associated velocity vector, p and ρ are the pressure and density fields, respectively,0ρ is the reference density,tν is the turbulent viscosity coefficient, S is the salinity andsσ is the turbulent Schmidt number, g is the gravitational acceleration in the vertical direction. The standard kε- turbulence model is used to describe the sub-grid process.

    Table 1 Initial layer depths for different numerical experiments

    Fig.3 The temporal evolution of the interface displacement for Case 1

    Fig.4 The temporal evolution of the interface displacement in Case 2

    Fig.5 The temporal evolution of the interface displacement for Case 3

    Fig.6 ISW generation processes for Case 1 from 0 s to 7 s, the time interval is 1 s

    The model settings in our numerical simulations are shown in Fig.1, similar to those in laboratory experiments performed by Chen et al.[16], but with some differences based on numerical considerations. The steel-framed wave flume used by Chen et al.[16]is 12 m long, with a cross section of 0.7 m high by 0.5 m wide. The flume has a slope to facilitate the study of its interaction with the internal waves, which is beyond the scope of the present study, and the total water depth is 0.4 m for all their laboratory experiments, thus our numerical model domain is set to be 6 m long and 0.4 m high. The salinity of the upper layer is 0 psu (fresh water) and 39 psu (brine water) for the lower layer, the corresponding densities are about 999 kg/m3and 1 029 kg/m3for the upper and lower layers, respectively, producing a stratified fluid environment of the Boussinesq parameter (the ratio of density difference and reference density) of 0.03, a typical value in the ocean environment[16]. For all our model studies,0x is set to be 5.7 m, the width of the right minor water compartment (perturbation region) is 0.3 m, the same value as in the laboratory experiment configuration. Three numerical experiments are performed, as shown in Table 1, initialized with different interface displacements. The differential equations for the momentum and salinity are discretized and implemented with the finite volume method within Open Field Operation and Manipulation (OpenFOAM) framework, the initial velocity field is zero in the domain, no flux and no slip conditions are applied along all the boundaries, thus the initial static state is driven by the initial baroclinic pressure gradient force caused by the inhomogeneous horizontal salinity distribution. All numerical experiments are integrated for 50 s, though our concern here is mainly on the internal wave generation, propagation and transformation processes, some reflection features are also briefly analyzed, as shown in the next section.

    3. The analysis of results

    3.1 General model results

    Fig.7 ISW generation processes for Case 2 from 0 s to 7 s, the time interval is 1 s

    The temporal evolution of the density interfaces (defined as the isopycnals of the averaged density of the upper and lower layers) in the whole domain are shown in Fig.3 (Case 1), Fig.4 (Case 2) and Fig.5 (Case 3), respectively. From Figs.3 and 4 it can be seen that the general internal wave generation, the propagation and reflection characteristics of Case 1 and Case 2 are similar. At first, the ISWs of depression are generated at the right side of the domain, with a leading ISW and small scale trailing interface oscillations, the ISW propagates to the left and reflects at the left boundary. It can be seen that the ISW propagation speed is well predicted by the linear theoretical model in the current scenarios. However, a distinct result shows in Fig.5 for Case 3, where it can be seen that firstly a wave of elevation is generated, as it propagates to the left, the amplitude of the leading elevation wave diminishes gradually and a train of depression ISWs emerge behind, the wave propagation speed is reduced as compared with that predicted by the linear model. The above numerical observations conform to the findings in the laboratory experiments of Chen et al.[16]. It may be more understandable if their conclusions are rephrased as follows: the initial type of the ISW is determined by the relative height of the perturbation density interface and the environmental density interface (a wave of depression is generated if the position of the perturbation density interface is lower than that of the environmental density interface, otherwise, a wave of elevation is generated.), while the final type of wave is determined by the position of the environmental density interface, namely, the relative thickness of the upper and lower layers (if the thickness of the upper layer is smaller than that of the lower layer, the wave of depression emerges finally, otherwise, the wave of elevation type is generated).

    3.2 Internal wave generation

    Fig.8 ISW generation processes for Case 3 from 0 s to 7 s, the time interval is 1 s

    Chen et al.[16]show snapshots of video images to exemplify the internal wave generation via the gravity collapse, the density interface is indicated by the dye injected beforehand. This method is commonly used in laboratory experiments to visualize the density interface evolution, however, as can be discerned from Chen et al.’s demonstration, two issues of deficiencies are evident, firstly, the dye injection is not continuous along the density interface, e.g., at the step location, which renders it impossible to know the temporal evolution of the interface; and secondly, the dye is extremely diffusive, thus it is hard for us to discern the interface after the initialization. These two issues can be easily addressed by the high-resolution numerical results, as shown in Fig.6 (Case 1), Fig.7 (Case 2) and Fig.8 (Case 3). In these three figures, the temporal evolutions of the density interface during the first seven seconds are presented to see the difference of the ISW generation processes, and the temporal interface evolutions predicted by the theoretical model are also presented for comparison. From these figures it can be seen that the ISW generation mechanisms are similar among cases, which resemble the lock-exchange phenomenon, while some distinctions are also evident. For example, the vortex structures shown in Figs.6(d) and 6(f) are obvious, while these phenomena are less evident at the corresponding time in Fig.7 and even less in Fig.8. These vortex structures might be caused by the Kelvin-Helmholtz instability at the stratified fluid interface. For all our numerical simulations, the stratifications are identical, the differences among these cases lie in the relative height between the environmental density interface and the perturbation density interface. The larger the relative height (as Case 1), the larger the shear across the density interface, the more prone to form the vortex structures it will be. It should be stressed that the accurate description of the ISW generation process via this mechanism is beyond the capability of any theoretical two-layer stratified fluid model, because the temporal evolution of the density interface, is not a function of the horizontal coordinate, as shown in Figs.6(d), 6(f), 7(c) and 7(d).

    It can be seen from Figs.6 and 7 that, although the theoretical model fails to predict the vortex structures during the initialization processes because of its inherent linear feature, the final wave position and the polarity feature are generally in good agreement with numerical results, as shown in Figs.6(h) and 7(h). This is, however, not true for Case 3, as shown in Fig.8(h), where the wave generated is of elevation type, while the environmental stratification favors the wave type of depression, which slows the corresponding propagation speed.

    3.3 Temporal evolution of interface

    The temporal evolution of the density interface at the middle of the computational domain is shown in Fig.9, where one sees the perfect ISW of depression and the accompanied trailing oscillations for Cases 1 and 2, and an elevation-type internal wave followed by the trailing internal wave of depression type for Case 3, in keeping with the recorded wave profiles in Chen et al.’s laboratory experiments qualitatively, which further validates our numerical results.

    3.4 The shape of ISW

    The propagation of ISWs can be described to the first order of magnitude in the amplitude via the weakly nonlinear Korteweg-de Vries (KdV) equation[20],

    A particular solution for Eqs.(13) is the solitary wave of the following form,

    where the phase speed c and the length scale are given by

    The classic KdV equation is derived under a strong restriction of weak nonlinearity and is likely not to give right results in certain circumstances when the solitons are extremely nonlinear[21]. The Extended Korteweg-de Vries (EKdV) equation accounts for a higher degree of nonlinearity,

    A solitary wave solution of Eqs.(16) is of the following form,

    where ν is a nonlinearity parameter that falls in the range from 0 to 1.

    The numerical results provide an opportunity for us to evaluate the performance of the KdV and EKdV equations in predicting the shape of the ISWs generated. Here the wave profiles at instants 15 s and 30 s are chosen for the first two cases. The reason to choose these two specific instants is that they are not affected by the ISW generation and reflection processes, and could be representative as the “young” and“mature” ISWs for their respective case scenario. The comparisons are shown in Fig.10.

    From the figure, one can see some common features among the two cases, e.g., the shape of the wave obtained by the KdV theory is steeper than that obtained by the EKdV theory for all comparisons. Distinct differences are also evident. We will make the comparisons for Case 2 first, as shown in Figs.10(b) and 10(d). Though the ISW amplitude decreases a little with the time, the differences from the model results, obtained by the KdV and EKdV theories, are very small. This can be explained by the fact that the initial forcing for Case 2 is weak, thus produces small amplitude ISWs. In that case, the effect of nonlinearity does not dominate, thus the wave shapes inferred from both KdV and EKdV equations conform to that produced from the numerical output. However, this does not hold for Case 1, as shown in Figs.10(a) and 10(c). In Fig.10(a), none of the KdV and EKdV theories predicts the ISW shape well, while for a later instant, the wave shape obtained by the EKdV theory agrees well with that obtained by the model, as shown in Fig.10(c). The reason for the distinct difference of the ISW shape obtained by the KdV and EKdV theories is that the amplitude is larger in Case 1 than that in Case 2, the nonlinearity effect is more important in such circumstance. Since both the classic KdV and EKdV theories are valid for steady state ISWs, when the effect of nonlinearity is equal to the effect of dispersion. Theinconsistency between the model output and the theoretical predictions in Fig.10(a) might be due to the fact that the initially generated ISW still adjusts itself to the environmental stratification, finally the wave shape predicted by the EKdV theory agrees better with the numerical result in Case 1 as the amplitude is larger. As a short conclusion, the applicability of the EKdV theory is wider than the KdV theory from the results shown in this study, as it incorporates a higher order nonlinear term.

    Fig.10 Comparison of model predicted wave shape and those predicted by the KdV and EKdV models

    3.5 Energy analysis

    One undisputable merit of the numerical simulation is that it could provide quantitative results for all quantities, such as the velocity, the temperature, the salinity, the turbulent kinetic energy and its dissipate rate in the whole computational domain, which are inaccessible for laboratory measurements as a whole. The availability of these quantitative data helps us better understand the related phenomena in depth, e.g., the energy evolution process, which will be discussed in this subsection.

    Fig.11 The temporal evolutions of total kinetic, potential energy and total mechanical energy (E) in Case 1, normalized by its initial available potential energy

    Fig.12 The temporal evolutions of the kinetic energy in Cases 1, 2 and 3, normalized by their initial available potential energy, respectively

    The temporal evolutions of the kinetic energy and the potential energy as well as the total mechanical energy of the system in Case 1 are shown in Fig.11. The overall tendency shown in this figure indicates that this is mainly a laminar process, where the turbulent mixing is not important and the total mechanical energy is approximately conserved. From the initial static state, the potential energy drops sharply and is converted to the kinetic energy. As the ISW forms and propagates leftward from the right boundary, the amplitude of the ISW decreases slightly, accompanied with a conversion from the kinetic energy to the potential energy. When the internal wave reaches the left boundary, the kinetic energy diminishes swiftly and is converted into the potential energy, an extrema of theinterface displacement is reached. The kinetic energy of the system increases again when the reflection process initiates.

    Fig.13 The spatial distribution of kinetic energy in Case 1 (upper panel), Case 2 (middle panel) and Case 3 (lower panel). The bin size along the x axis is 0.01 m

    The temporal evolution of the kinetic energy in different cases are shown in Fig.12, with the general trends of development being similar, as analyzed above, however, three distinct features are evident. Firstly, the kinetic energy booms up to more than 50% of its initial APE during the generation phase, and then drops a little to around 50%, which could probably attribute to the inertial effect during the gravity collapse phase and the later self adjustment. Secondly, the system’s kinetic energy is approximately one half of its initial available potential energy during the ISW propagation phase, accompanied with a slight conversion from the kinetic energy to the potential energy as the ISW adjusts to the environmental stratification. Thirdly, it is worth noting that the occurrence of the maximum total kinetic energy is not at the ISW generation phase, but at the ISW reflection stage.

    The spatial distributions of the kinetic energy in the three cases at instants 10 s, 20 s and 30 s are presented in Fig.13, from which it can be seen that in Cases 1 and 2 the kinetic energy packets are approximately conserved and propagate with the ISW speed, though their spatial distributions flatten slightly with time. This indicates that the kinetic energy of the system comes mainly from the ISW induced currents. The kinetic energies contained in the ISWs in Case 1 (Case 2), defined as the integration of the kinetic energy from the left end to the extrema location just behind leading wave, are approximately 95.09%, 90.49%, 89.13% (94.98%, 90.31%, 88.31%) of their total kinetic energies at instants 10s, 20s, and 30s, respectively. Though the temporal evolution of the total kinetic energy in case 3 resembles that in Cases 1 and 2, as shown in Fig.12, its spatial distribution diverges dramatically as shown in the lower panel of Fig.13. The kinetic energy distribution collapses spatially as the time elapses.

    4. Conclusions

    In this paper, systematic theoretical and numerical studies are carried out for the ISW generation and propagation characteristics by the gravity collapse mechanism. The numerical results could reproduce all the phenomena observed in laboratory experiments performed by Chen et al.[16]and the numerical model is used to overcome the deficiencies inherent in their laboratory experiments, permitting a quantitative understanding of this ISW generation mechanism and evolution features.

    The theoretical model indicates that a rectangletype wave is generated, which propagates back and forth in the domain. The numerical results produce more realistic outputs, from which it can be seen that the initial ISW type is determined by the relative height of the perturbation density interface and the environmental density interface, while the final ISW type depends on the relative height of the upper layer and the lower layer of the environmental fluid. One notable phenomenon is that the propagation speed of the ISW inferred from the theoretical model conforms to that of the numerical results considerably well if no ISW polarity conversion process occurs. The difference of the wave shapes predicted by the KdV and EKdV theories is obvious when the amplitude of the ISW is relatively large, the prediction based on the EKdV theory is more appropriate under this circumstance as it incorporates a higher-order term of nonlinearity. The result shows that the total mechanical energy is generally conserved, a rapid conversion from the potential energy to the kinetic energy occurs during the ISW generation phase. As the ISW forms and propagates away from its source region, the total kinetic energy is approximately 50% of its initial APE, though a slight conversion to the potential energy is accompanied. The spatial distributions of the kinetic energy tend to flatten slowly as the ISW adjusts itself to the environmental stratification if no polarity conversion occurs, on the other hand, the spatial distribution of the kinetic energy collapses very quickly and covers a larger area for the polarity conversion case.

    The present study demonstrates that the laboratory scale numerical model is an indispensable tool that can further extend our understanding of the interested process.

    Acknowledgement

    This work was supported by the State Key Laboratory of Tropical Oceanography, South China Sea Institute of Oceanology, Chinese Academy of Sciences (Grant No. LTO1104).

    [1] LI Qun, XU Zhao-ting. A nonhydrostatic numerical model for density stratified flow and its applications[J]. Journal of Hydrodynamics, 2008, 20(6): 706-712.

    [2] CAI S., XIE J. A propagation model for the internal solitary waves in the northern South China Sea[J]. Journal of Geophysical Research, 2010, 115(12): C12074.

    [3] CAI Shuqun, HE Jianling and XIE Jieshuo. Recent decadal progress of the study on internal solitons in the South China Sea[J]. Advances in Earth Science, 2011, 26(7): 703-710(in Chinese).

    [4] LI Juan, GU Xing-fa and YU Tao et al. Simulation investigation on the nonlinear interaction of internal waves in Luzon Strait via the variable-coefficient model[J]. Chinese Journal of Hydrodynamics, 2011, 26(2): 157-166(in Chinese).

    [5] NASH J., MOUM J. N. River plumes as a source of large-amplitude internal waves in the coastal ocean[J]. Nature, 2005, 437(7057): 400-403.

    [6] BOEGMAN L., IVEY G. N. and IMBERGER J. The energetics of large-scale internal wave degeneration in lakes[J]. Journal of Fluid Mechanics, 2005, 531: 159-180.

    [7] CHEN C. Y., HSU J. R. C. Interaction between internal waves and a permeable seabed[J]. Ocean Engineering, 2005, 32(5-6): 587-621.

    [8] YANG Y. J., TANG T. Y. and CHANG M. H. et al. Solitons northeast of Tung-Sha Island during the ASIAEX pilot studies[J]. IEEE Journal of Oceanic Engineering, 2004, 29(4): 1182-1199.

    [9] ZHAO Z., ALFORD M. H. and GIRTON J. et al. Internal tides around the Hawaiian Ridge estimated from multisatellite altimetry[J]. Journal of Geophysical Research, 2011, 116(12): C12039.

    [10] BUIJSMAN M. C., KANARSKA Y. and MCWILLIAMS J. C. On the generation and evolution of nonlinear internal waves in the South China Sea[J]. Journal of Geophysical Research, 2010, 115(2): C02012.

    [11] CHEN C. Y., HSU J. R. C. and CHEN H. H. et al. Laboratory observations on internal solitary wave evolution on steep and inverse uniform slopes[J]. Ocean Engineering, 2007, 34(1): 157-170.

    [12] CHEN C. Y., HSU J. R. C. and CHENG M. H. et al. An investigation on internal solitary waves in a two-layer fluid: Propagation and reflection from steep slopes[J]. Ocean Engineering, 2007, 34(1): 171-184.

    [13] CHEN C. Y., Hsu J. R. C. and CHENG M. H. et al. Experiments on mixing and dissipation in internal solitary waves over two triangular obstacles[J]. Environmental Fluid Mechanics, 2008, 8(3): 199-214.

    [14] CHEN C. Y. Amplitude decay and energy dissipation due to the interaction of internal solitary waves with a triangular obstacle in a two-layer fluid system: The blockage parameter[J]. Journal of Marine Science and Technology, 2009, 14(4): 499-512.

    [15] CHEN C. Y. Statistical and dynamical analyses of propagation mechanisms of solitary internal waves in a two-layer stratification[J]. Journal of Marine Science and Technology, 2011, 16(1): 100-114.

    [16] CHEN C. Y., HSU J. R. C. and CHEN C. W. et al. Generation of internal solitary wave by gravity collapses[J]. Journal of Marine Science and Technology, 2007, 15(1): 1-7.

    [17] CHENG M. H., HSU J. R. C. and CHEN C. Y. Laboratory experiments on waveform inversion of an internal solitary wave over a slope-shelf[J]. Environmental Fluid Mechanics, 2011, 11(4): 353-384.

    [18] LYNETT P. J., LIU P. L.- F. A two-dimensional, depthintegrated model for internal wave propagation over variable bathymetry[J]. Wave Motion, 2002, 36(3): 221-240.

    [19] CHEN C. Y., HSU J. R. C. and CHEN C. W. et al. Numerical model of internal solitary wave evolution on impermeable variable seabed in a stratified two-layer fluid system[J]. China Ocean Engineering, 2006, 20(2): 303-313.

    [20] HORN D. A., IMBERGER J. and IVEY G. N. et al. A weakly nonlinear model of long internal waves in closed basins[J]. Journal of Fluid Mechanics, 2002, 467: 269-287.

    [21] SMALL R. J., HORNBY R. P. A comparison of weakly and fully non-linear models of the shoaling of a solitary internal wave[J]. Ocean Modelling, 2005, 8(4): 395-416.

    Appendix

    Transformation of (,)xtη

    where

    10.1016/S1001-6058(11)60276-X

    * Project supported by the National Natural Science Foundation of China (Grant Nos. 61072145, 41176016), the Fund for Creative Research Groups by National Natural Science Foundation of China (Grant No. 41121064).

    Biography: LIN Zhen-hua (1981-), Male, Ph. D.

    av在线蜜桃| 国产精品99久久99久久久不卡| 亚洲男人的天堂狠狠| 亚洲国产欧美一区二区综合| or卡值多少钱| 一个人免费在线观看电影 | 动漫黄色视频在线观看| 好看av亚洲va欧美ⅴa在| 99久国产av精品| 很黄的视频免费| 午夜亚洲福利在线播放| 19禁男女啪啪无遮挡网站| 日韩欧美在线乱码| 老司机午夜福利在线观看视频| 99riav亚洲国产免费| 亚洲五月婷婷丁香| 19禁男女啪啪无遮挡网站| 99久久精品热视频| 午夜福利18| 亚洲欧美精品综合一区二区三区| 精品国产乱子伦一区二区三区| bbb黄色大片| 国产精品亚洲av一区麻豆| 最新在线观看一区二区三区| 桃色一区二区三区在线观看| 巨乳人妻的诱惑在线观看| 91av网一区二区| 日韩欧美精品v在线| 18禁观看日本| 可以在线观看毛片的网站| 级片在线观看| 又爽又黄无遮挡网站| 亚洲人与动物交配视频| 给我免费播放毛片高清在线观看| a级毛片在线看网站| 两个人视频免费观看高清| 在线观看一区二区三区| 搡老妇女老女人老熟妇| 亚洲精品粉嫩美女一区| 我要搜黄色片| 一个人免费在线观看电影 | 90打野战视频偷拍视频| 色综合亚洲欧美另类图片| 天天躁狠狠躁夜夜躁狠狠躁| 日韩成人在线观看一区二区三区| 亚洲国产精品999在线| 99久久久亚洲精品蜜臀av| 国产高清视频在线播放一区| 国产一区二区三区视频了| 大型黄色视频在线免费观看| av欧美777| 变态另类丝袜制服| 91久久精品国产一区二区成人 | 听说在线观看完整版免费高清| 叶爱在线成人免费视频播放| 午夜精品在线福利| 亚洲第一欧美日韩一区二区三区| 脱女人内裤的视频| 亚洲国产精品合色在线| 午夜视频精品福利| 五月伊人婷婷丁香| 91麻豆精品激情在线观看国产| 麻豆国产97在线/欧美| 亚洲最大成人中文| 久久久久性生活片| 一进一出抽搐动态| 精品久久蜜臀av无| 99久久99久久久精品蜜桃| 搡老熟女国产l中国老女人| 久久久久国内视频| 国产真实乱freesex| 成人性生交大片免费视频hd| 国产精品自产拍在线观看55亚洲| 不卡一级毛片| 亚洲av成人一区二区三| 国产真实乱freesex| 欧美日本亚洲视频在线播放| 国产精品98久久久久久宅男小说| 亚洲中文字幕日韩| 国产乱人视频| 久久精品91蜜桃| 欧美日本亚洲视频在线播放| 在线看三级毛片| 国产高潮美女av| 色综合亚洲欧美另类图片| 亚洲成人久久性| 小说图片视频综合网站| 日韩欧美在线乱码| 国产精品 欧美亚洲| 日本黄色片子视频| 99久国产av精品| 国产亚洲欧美98| 免费观看的影片在线观看| 久久久久久久久久黄片| 长腿黑丝高跟| 日韩欧美一区二区三区在线观看| 久久精品国产清高在天天线| 狠狠狠狠99中文字幕| 性色avwww在线观看| 成人亚洲精品av一区二区| 国产熟女xx| 日韩免费av在线播放| 久久精品国产清高在天天线| 亚洲黑人精品在线| 老司机午夜福利在线观看视频| 小蜜桃在线观看免费完整版高清| 99久久久亚洲精品蜜臀av| avwww免费| 国产黄色小视频在线观看| 精品欧美国产一区二区三| 亚洲avbb在线观看| 一个人看的www免费观看视频| 亚洲国产精品999在线| 欧美日韩中文字幕国产精品一区二区三区| 久99久视频精品免费| 久久人妻av系列| 亚洲av成人一区二区三| 99久久久亚洲精品蜜臀av| 亚洲 国产 在线| av天堂在线播放| 午夜免费成人在线视频| 欧美乱色亚洲激情| 人妻丰满熟妇av一区二区三区| 国产精品亚洲美女久久久| 99国产精品99久久久久| 国产激情久久老熟女| 老司机午夜福利在线观看视频| 老司机福利观看| 国产成人aa在线观看| x7x7x7水蜜桃| 久久久久性生活片| 99热6这里只有精品| 精品午夜福利视频在线观看一区| 国产高清videossex| 国产在线精品亚洲第一网站| 久久午夜综合久久蜜桃| 99久久精品一区二区三区| 国产成人精品久久二区二区免费| a级毛片在线看网站| 九九热线精品视视频播放| 国产亚洲欧美在线一区二区| 男人舔女人下体高潮全视频| 国产又色又爽无遮挡免费看| 日韩欧美精品v在线| 日本精品一区二区三区蜜桃| 亚洲熟妇中文字幕五十中出| 亚洲无线在线观看| а√天堂www在线а√下载| 午夜免费激情av| 国产视频内射| 欧美日韩精品网址| 午夜福利免费观看在线| av黄色大香蕉| 日韩成人在线观看一区二区三区| 后天国语完整版免费观看| 在线观看舔阴道视频| 国产一区在线观看成人免费| 精品一区二区三区视频在线 | 亚洲精华国产精华精| 国产高清有码在线观看视频| 亚洲黑人精品在线| 久久天躁狠狠躁夜夜2o2o| 欧美成人性av电影在线观看| 亚洲av熟女| 日韩有码中文字幕| 日韩有码中文字幕| 欧美3d第一页| 精品国内亚洲2022精品成人| 成人特级黄色片久久久久久久| 美女高潮喷水抽搐中文字幕| 欧美激情在线99| 97超级碰碰碰精品色视频在线观看| 欧美日本视频| 村上凉子中文字幕在线| 51午夜福利影视在线观看| 老司机福利观看| 国语自产精品视频在线第100页| 欧美午夜高清在线| 欧美成狂野欧美在线观看| 久久人妻av系列| 国产爱豆传媒在线观看| a级毛片a级免费在线| 国产精品一区二区免费欧美| 真实男女啪啪啪动态图| 黑人操中国人逼视频| 一卡2卡三卡四卡精品乱码亚洲| 真人做人爱边吃奶动态| 少妇熟女aⅴ在线视频| 好男人电影高清在线观看| www日本黄色视频网| 天堂动漫精品| 成人av在线播放网站| 真实男女啪啪啪动态图| 超碰成人久久| 久久人人精品亚洲av| 久久精品综合一区二区三区| 又大又爽又粗| 色视频www国产| 国内久久婷婷六月综合欲色啪| 欧美日韩精品网址| 久久天躁狠狠躁夜夜2o2o| 亚洲av中文字字幕乱码综合| 一a级毛片在线观看| 真人一进一出gif抽搐免费| 成人永久免费在线观看视频| 国产一区二区在线av高清观看| 人妻久久中文字幕网| 精品国产超薄肉色丝袜足j| 最近视频中文字幕2019在线8| 级片在线观看| 97碰自拍视频| 男插女下体视频免费在线播放| 亚洲成av人片在线播放无| 国产主播在线观看一区二区| 在线免费观看不下载黄p国产 | 日韩欧美免费精品| 波多野结衣巨乳人妻| 久久99热这里只有精品18| 精品国产超薄肉色丝袜足j| 在线观看美女被高潮喷水网站 | 午夜精品在线福利| 一二三四在线观看免费中文在| 亚洲熟妇熟女久久| 亚洲第一电影网av| 国产真人三级小视频在线观看| 成人av在线播放网站| 婷婷亚洲欧美| 在线观看美女被高潮喷水网站 | 免费看十八禁软件| 国产1区2区3区精品| 国产一区在线观看成人免费| 成年免费大片在线观看| 国产精品99久久久久久久久| 国产精品香港三级国产av潘金莲| 欧美日韩中文字幕国产精品一区二区三区| 亚洲 欧美 日韩 在线 免费| 18禁国产床啪视频网站| 欧美激情在线99| 午夜福利成人在线免费观看| 一二三四在线观看免费中文在| 国产又色又爽无遮挡免费看| 男女之事视频高清在线观看| 成人三级做爰电影| 在线观看免费视频日本深夜| 五月玫瑰六月丁香| 别揉我奶头~嗯~啊~动态视频| 最新美女视频免费是黄的| 国产一区二区三区视频了| 视频区欧美日本亚洲| 老司机午夜福利在线观看视频| 在线观看66精品国产| www.999成人在线观看| 18禁黄网站禁片午夜丰满| 亚洲自偷自拍图片 自拍| 99久久精品热视频| 无人区码免费观看不卡| 在线观看舔阴道视频| 婷婷六月久久综合丁香| 免费观看的影片在线观看| 9191精品国产免费久久| 人人妻,人人澡人人爽秒播| 亚洲av熟女| 亚洲av成人精品一区久久| 久久精品人妻少妇| 亚洲av成人av| 亚洲人成电影免费在线| 在线国产一区二区在线| 1000部很黄的大片| www.999成人在线观看| 免费在线观看亚洲国产| 欧美不卡视频在线免费观看| 91av网一区二区| 这个男人来自地球电影免费观看| 国产欧美日韩精品一区二区| 亚洲成人中文字幕在线播放| 欧美成狂野欧美在线观看| 久久久久免费精品人妻一区二区| 草草在线视频免费看| 精品久久久久久久末码| 看片在线看免费视频| 精品久久久久久,| 少妇裸体淫交视频免费看高清| 国产成人一区二区三区免费视频网站| 亚洲色图av天堂| av欧美777| 国产欧美日韩一区二区精品| 国产精品亚洲美女久久久| 国产精品久久久人人做人人爽| 久9热在线精品视频| 搞女人的毛片| 亚洲18禁久久av| 久久久国产成人精品二区| av福利片在线观看| 久久久久久久久久黄片| 国产精品美女特级片免费视频播放器 | 三级毛片av免费| 无限看片的www在线观看| 男女做爰动态图高潮gif福利片| 午夜福利在线观看吧| 制服人妻中文乱码| 国产爱豆传媒在线观看| 老汉色av国产亚洲站长工具| 国产单亲对白刺激| 亚洲电影在线观看av| 99国产综合亚洲精品| 亚洲av电影不卡..在线观看| 一区二区三区激情视频| 哪里可以看免费的av片| 精品久久蜜臀av无| 精品国产超薄肉色丝袜足j| 亚洲精品美女久久av网站| 成人性生交大片免费视频hd| svipshipincom国产片| 一个人看的www免费观看视频| 国产免费av片在线观看野外av| 日本 欧美在线| 99久国产av精品| 不卡av一区二区三区| 精华霜和精华液先用哪个| 欧美黑人欧美精品刺激| 久久九九热精品免费| 久久国产精品人妻蜜桃| 亚洲无线观看免费| 午夜福利18| 18禁美女被吸乳视频| 麻豆国产av国片精品| 在线国产一区二区在线| 无遮挡黄片免费观看| 亚洲av第一区精品v没综合| 免费看日本二区| 午夜福利在线观看免费完整高清在 | 国产精品久久电影中文字幕| 99久久精品一区二区三区| 国产高清视频在线播放一区| 欧美黄色淫秽网站| 久久欧美精品欧美久久欧美| 精品一区二区三区视频在线 | 深夜精品福利| 天堂影院成人在线观看| 婷婷六月久久综合丁香| 久久久色成人| 亚洲国产精品合色在线| 久久中文字幕一级| 亚洲一区二区三区不卡视频| 黄色视频,在线免费观看| 天堂av国产一区二区熟女人妻| 美女黄网站色视频| 午夜成年电影在线免费观看| 免费av毛片视频| 法律面前人人平等表现在哪些方面| 成人三级黄色视频| 成人鲁丝片一二三区免费| 欧美又色又爽又黄视频| 亚洲国产色片| 国产激情欧美一区二区| 女生性感内裤真人,穿戴方法视频| 亚洲人成伊人成综合网2020| 成人国产一区最新在线观看| 国产真实乱freesex| 欧美日韩黄片免| 18禁裸乳无遮挡免费网站照片| 手机成人av网站| а√天堂www在线а√下载| 99国产精品一区二区三区| 成在线人永久免费视频| 亚洲无线观看免费| 两个人视频免费观看高清| 亚洲专区字幕在线| 一个人免费在线观看电影 | 99国产综合亚洲精品| 熟女人妻精品中文字幕| 中文亚洲av片在线观看爽| 色吧在线观看| 男插女下体视频免费在线播放| 午夜精品久久久久久毛片777| 国模一区二区三区四区视频 | 亚洲人成电影免费在线| 亚洲男人的天堂狠狠| 亚洲国产高清在线一区二区三| 精品99又大又爽又粗少妇毛片 | 精品国产超薄肉色丝袜足j| 国产精品av久久久久免费| 久久人妻av系列| 搡老岳熟女国产| 在线永久观看黄色视频| 久9热在线精品视频| 日韩国内少妇激情av| 九色成人免费人妻av| 欧美大码av| 久久天堂一区二区三区四区| 国产人伦9x9x在线观看| 精品乱码久久久久久99久播| 国产 一区 欧美 日韩| 88av欧美| 欧美一级a爱片免费观看看| 色吧在线观看| 国产aⅴ精品一区二区三区波| 久久精品国产清高在天天线| 久久久久久久久免费视频了| 国产蜜桃级精品一区二区三区| 色av中文字幕| 国产91精品成人一区二区三区| 亚洲av第一区精品v没综合| 欧美黄色片欧美黄色片| 亚洲第一电影网av| 99热只有精品国产| 国产精品野战在线观看| 露出奶头的视频| 一级毛片女人18水好多| 琪琪午夜伦伦电影理论片6080| 99国产精品99久久久久| 午夜亚洲福利在线播放| 曰老女人黄片| 免费观看的影片在线观看| 午夜福利在线观看免费完整高清在 | 日韩人妻高清精品专区| x7x7x7水蜜桃| 欧美在线黄色| 国产av一区在线观看免费| 欧美黄色片欧美黄色片| 可以在线观看毛片的网站| 欧美一区二区国产精品久久精品| 亚洲国产欧美网| 国产伦人伦偷精品视频| 精品福利观看| 男人和女人高潮做爰伦理| av在线蜜桃| 美女 人体艺术 gogo| 亚洲精品在线观看二区| 熟女少妇亚洲综合色aaa.| 一a级毛片在线观看| 久久久久精品国产欧美久久久| 国产91精品成人一区二区三区| 人妻丰满熟妇av一区二区三区| 19禁男女啪啪无遮挡网站| 91在线观看av| 久久久久久久精品吃奶| 国产亚洲精品久久久久久毛片| 久久九九热精品免费| 国产高潮美女av| 两人在一起打扑克的视频| 少妇人妻一区二区三区视频| 亚洲国产精品合色在线| 精品国产乱码久久久久久男人| 搞女人的毛片| 最近最新免费中文字幕在线| 国产午夜福利久久久久久| 午夜精品在线福利| 亚洲精品一区av在线观看| 成年免费大片在线观看| 久久99热这里只有精品18| 国产精品一区二区免费欧美| 中文字幕熟女人妻在线| 亚洲欧美日韩东京热| 在线观看午夜福利视频| 一区二区三区高清视频在线| 首页视频小说图片口味搜索| 国产一区二区在线观看日韩 | 丁香六月欧美| 免费大片18禁| 国产黄片美女视频| 给我免费播放毛片高清在线观看| 久久久久亚洲av毛片大全| 亚洲精品456在线播放app | 天堂av国产一区二区熟女人妻| 亚洲性夜色夜夜综合| 老司机深夜福利视频在线观看| 母亲3免费完整高清在线观看| 好看av亚洲va欧美ⅴa在| 啪啪无遮挡十八禁网站| 12—13女人毛片做爰片一| 久久这里只有精品19| 国产男靠女视频免费网站| 午夜久久久久精精品| 日本与韩国留学比较| 亚洲欧美日韩高清专用| 亚洲欧美精品综合一区二区三区| 香蕉丝袜av| 人妻丰满熟妇av一区二区三区| 国产午夜精品论理片| 99热只有精品国产| 精品日产1卡2卡| 人人妻人人澡欧美一区二区| 亚洲精品久久国产高清桃花| 亚洲午夜精品一区,二区,三区| 俄罗斯特黄特色一大片| 真人做人爱边吃奶动态| 精品乱码久久久久久99久播| 日本三级黄在线观看| 亚洲美女黄片视频| 在线免费观看不下载黄p国产 | 91久久精品国产一区二区成人 | 悠悠久久av| 精品一区二区三区四区五区乱码| 国产成人系列免费观看| svipshipincom国产片| 亚洲精品粉嫩美女一区| 日日摸夜夜添夜夜添小说| or卡值多少钱| 国产毛片a区久久久久| 老鸭窝网址在线观看| 亚洲成人免费电影在线观看| 很黄的视频免费| 少妇人妻一区二区三区视频| 亚洲天堂国产精品一区在线| 久久久国产精品麻豆| 国产精品国产高清国产av| 国产精品久久视频播放| 国内少妇人妻偷人精品xxx网站 | 五月玫瑰六月丁香| 亚洲无线观看免费| 嫁个100分男人电影在线观看| www.自偷自拍.com| 美女扒开内裤让男人捅视频| 五月玫瑰六月丁香| 成人一区二区视频在线观看| 51午夜福利影视在线观看| 日本黄色视频三级网站网址| 三级毛片av免费| 久久国产乱子伦精品免费另类| 色av中文字幕| www.999成人在线观看| 国产伦精品一区二区三区四那| 午夜免费观看网址| 看片在线看免费视频| 99久久无色码亚洲精品果冻| 一区二区三区高清视频在线| 国产激情久久老熟女| 免费高清视频大片| 欧美乱码精品一区二区三区| 亚洲国产欧美一区二区综合| www.自偷自拍.com| 亚洲精品乱码久久久v下载方式 | 18禁黄网站禁片午夜丰满| 99re在线观看精品视频| 伦理电影免费视频| 黄片小视频在线播放| 成人一区二区视频在线观看| 91在线精品国自产拍蜜月 | 久久久久亚洲av毛片大全| 欧美成人免费av一区二区三区| 男女床上黄色一级片免费看| 亚洲欧美日韩东京热| 久久九九热精品免费| 成人三级黄色视频| 精品国内亚洲2022精品成人| 日日干狠狠操夜夜爽| 国产午夜福利久久久久久| 欧美黄色淫秽网站| www.999成人在线观看| 搞女人的毛片| 午夜精品久久久久久毛片777| 18禁黄网站禁片免费观看直播| 亚洲av成人精品一区久久| 我要搜黄色片| 色哟哟哟哟哟哟| 欧美黄色淫秽网站| 免费高清视频大片| 不卡av一区二区三区| 午夜精品久久久久久毛片777| 午夜a级毛片| 中文字幕久久专区| 亚洲七黄色美女视频| 国产精品久久久久久久电影 | 欧美最黄视频在线播放免费| 国产成人精品久久二区二区91| 亚洲中文字幕日韩| 欧美3d第一页| 久久久水蜜桃国产精品网| 一个人看视频在线观看www免费 | 成人欧美大片| 最近最新中文字幕大全免费视频| 国产精品综合久久久久久久免费| 特大巨黑吊av在线直播| 色综合站精品国产| 日本成人三级电影网站| 国产单亲对白刺激| 麻豆国产97在线/欧美| 老司机在亚洲福利影院| 国内精品久久久久久久电影| 国产精品,欧美在线| 亚洲最大成人中文| 欧美日韩福利视频一区二区| 丰满的人妻完整版| 成人精品一区二区免费| 亚洲精品久久国产高清桃花| 国产三级在线视频| 人妻丰满熟妇av一区二区三区| 久久亚洲真实| 亚洲专区中文字幕在线| 国产综合懂色| 精品免费久久久久久久清纯| 亚洲 国产 在线| av黄色大香蕉| 91在线精品国自产拍蜜月 | 成在线人永久免费视频| 中文亚洲av片在线观看爽| 婷婷丁香在线五月| АⅤ资源中文在线天堂| 久久国产精品影院| 欧美黄色淫秽网站| 久久香蕉国产精品| 亚洲成人免费电影在线观看| www国产在线视频色| 欧美一区二区国产精品久久精品| 小说图片视频综合网站| 中出人妻视频一区二区| 国产三级黄色录像| 亚洲国产高清在线一区二区三| 亚洲av电影不卡..在线观看| 午夜福利免费观看在线| 国产99白浆流出| 亚洲国产中文字幕在线视频| 叶爱在线成人免费视频播放| 人妻夜夜爽99麻豆av|