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

    NUMERICAL SIMULATION OF SOLITARY WAVE RUN-UP AND OVERTOPPING USING BOUSSINESQ-TYPE MODEL*

    2012-08-22 08:32:14TSUNGWenShuo

    TSUNG Wen-Shuo

    Department of Hydraulic and Ocean Engineering, National Cheng Kung University, Tainan 701, E-mail: slavik_1942@hotmail.com

    HSIAO Shih-Chun

    Department of Hydraulic and Ocean Engineering, National Cheng Kung University, Tainan 701 LIN Ting-Chieh

    Tainan Hydraulics Laboratory, National Cheng Kung University, Tainan 701

    (Received May 17, 2012, Revised August 30, 2012)

    NUMERICAL SIMULATION OF SOLITARY WAVE RUN-UP AND OVERTOPPING USING BOUSSINESQ-TYPE MODEL*

    TSUNG Wen-Shuo

    Department of Hydraulic and Ocean Engineering, National Cheng Kung University, Tainan 701, E-mail: slavik_1942@hotmail.com

    HSIAO Shih-Chun

    Department of Hydraulic and Ocean Engineering, National Cheng Kung University, Tainan 701 LIN Ting-Chieh

    Tainan Hydraulics Laboratory, National Cheng Kung University, Tainan 701

    (Received May 17, 2012, Revised August 30, 2012)

    In this article, the use of a high-order Boussinesq-type model and sets of laboratory experiments in a large scale flume of breaking solitary waves climbing up slopes with two inclinations are presented to study the shoreline behavior of breaking and non-breaking solitary waves on plane slopes. The scale effect on run-up height is briefly discussed. The model simulation capability is well validated against the available laboratory data and present experiments. Then, serial numerical tests are conducted to study the shoreline motion correlated with the effects of beach slope and wave nonlinearity for breaking and non-breaking waves. The empirical formula proposed by Hsiao et al. for predicting the maximum run-up height of a breaking solitary wave on plane slopes with a wide range of slope inclinations is confirmed to be cautious. Furthermore, solitary waves impacting and overtopping an impermeable sloping seawall at various water depths are investigated. Laboratory data of run-up height, shoreline motion, free surface elevation and overtopping discharge are presented. Comparisons of run-up, run-down, shoreline trajectory and wave overtopping discharge are made. A fairly good agreement is seen between numerical results and experimental data. It elucidates that the present depth-integrated model can be used as an efficient tool for predicting a wide spectrum of coastal problems.

    Boussinesq equations, solitary wave, run-up, shoreline, scale effect, overtopping

    Introduction

    The terrific catastrophe of Great Indian Ocean tsunami in 2004[1-3]and the 2011 Tohoku earthquake tsunami[4]raised significant concern all over the world. Based on field evidence, the water inundation and debris flow that accompany a tsunami wave affect water motion in coastal areas. The run-up/run-down heights are thus extensively studied in tsunami wave research.

    Due to their simplicity, tsunami-like solitary waves and shallow water wave equations (SWEs) are frequently utilized to investigate the run-up and rundown heights of tsunami waves[1,5,6]. The informationobtained in these studies has led to significant progress in warning systems established for tsunami disaster mitigation. However, owing to the physical limitations of SWEs (i.e., hydrostatic pressure and uniform horizontal velocity profile assumptions), SWE-based models are unsuitable for problems in which both frequency dispersion and wave nonlinearity are of great concern[7]. In addition, large-scale experiments on solitary waves propagating over a mild slope have rarely been conducted due to the limits of laboratory facilities. Hence, the reliability of numerical model analyses is questionable[1,6].

    Researchers have developed highly robust and accurate solvers for run-up and overtopping with the SWEs[8]. Nevertheless, few Boussinesq-based models have been used to study overtopping, one such model was proposed by Stansby and Feng[9]. One potential reason is that to simulate dispersive wind waves, the Boussinesq model usually requires a complicatednumerical scheme for accuracy which does not lend itself to capturing the complex flow patterns that are common with overtopping. Navier-Stokes-based approaches[10,11]have been shown to be very accurate in predicting overtopping rates on a small scale. However, solving the Navier-Stokes equations has led to an extremely high computational cost, and their using range is generally restricted to a low number of precise waves and structure conditions.

    The objectives of this article are thus twofold: to investigate the shoreline properties caused by breaking and non-breaking solitary waves on plane beaches with various inclinations and to understand whether the Boussinesq-type equations can well simulate the wave overtopping the coastal structures. Numerical simulations based on a high-order Boussinesq model (COULWAVE, Cornell University long and intermediate wave modeling) by Kim et al.[12]and sets of laboratory experiments are conducted in the super flume of breaking solitary waves on 1:20 and 1:60 plane slopes are reported. We not only present new laboratory data of maximum run-up heights on two comparably gentle slopes but also measure the corresponding shoreline trajectory of swash motions for the numerical model calibration. The properties of scale effect on run-up heights of two comparable beaches between small flume[5]and large tank (1:20 and 1:60) are briefly discussed. In particular, present study intends to extend the results of Li and Raichlen[6]based on the nonlinear shallow water wave equations for an extensive concern of maximum run-up heights due to breaking and non-breaking solitary waves on plane beaches. By using the present numerical model with the available experiments and the given laboratory works, their results will be confirmed and extended.

    The rest of this study is organized as follows. The Boussinesq-type model and numerical scheme are briefly described in Section 1. The present laboratory methods and scale effects are given in Section 2. In Section 3, the model simulation capability is examined with the available laboratory data and present experiments. Both breaking and non-breaking waves are considered. For the non-breaking cases, the energy dissipation algorithms built in this Boussinesq model are calibrated against the laboratory data by Synolakis[5], Li and Raichlen[6]and the given experiments, respectively. Section 4 then outlines a series of numerical experiments to widely discuss the shoreline kinematics of breaking and non-breaking solitary waves on sloping beaches. The properties of run-up, run-down and shoreline trajectory are also reported. The applicability of the empirical formula proposed by Hsiao et al.[1]for generally predicting the maximum run-up height for breaking solitary waves on plane slopes is also confirmed. Furthermore, solitary waves impacting and overtopping an impermeable seawall at various water depths are investigated. Comparisons are made with experimental data. Finally, Section 5 concludes this paper.

    1. Numerical model

    The fundamental Boussinesq modeling approach is introduced in this section. The governing equations with viscous effect, and the numerical scheme are presented.

    1.1 Governing equations

    The fully nonlinear, weakly dispersive, weakly rotational depth-integrated Boussinesq-type equations in conservation form are expressed in one horizontal dimension as given by Kim et al.[12]

    where H=ζ+h is the total water depth,andrepresent the second-order terms in the depthintegrated x and y horizontal momentum equations respectively, and Hccontains the second-order terms in the continuity equation. The terms are given by

    where

    where (ψx,ψy)=ψThe second-order terms in the continuity equation are given by

    The physical meanings of the high-order terms can be found in Ref.[12].

    In the present study, a quadratic friction equation was used to approximate the bottom stress

    For simulating subgrid-scale dissipation, the horizontal eddy viscosity is modeled using the Smagorinsky model, which is given by

    The vertical eddy viscosity is given by

    1.2 Numerical scheme

    In this subsection, only a brief conceptual overview of the numerical solution is given. For more details, refer to Ref.[12].

    Considering the extended Boussinesq-type equations, in which there are first- to third-order spatial derivatives, the time integration should be fourth-order accurate to prevent numerical truncation errors from the same order as included derivatives. A third-order Adams-Bashforth predictor and the fourth-order Adams-Moulton corrector scheme are used for the time integration.

    The predictor step is

    The corrector step is

    where Ch=κ/6 is used, u*' is the friction velocity, κ is the Von Kármán constant, a value of 0.4 being adopted in this study.

    After each predictor and corrector step, a matrix solve r is used to solve P and Q. All the computed physical values are cell-averagedvalues since a cellaveraged finite-volume method is implemented to represent the spatial numerical grid.

    To calculate the leading-order (shallow water) terms in the governing equations, a fourth-order compact Monotone Upstream-centered Scheme for Conservation Laws-Total Variation Diminishing (MUSCL TVD) scheme is used to establish the interfacial physical values. Further, the Courant number determines the calculation time step used in the model, a value of 0.5 generally leads to stability and convergence, but for simulations with higher nonlinear waves, a value of as low as 0.1 may be required for stability. A value of 0.25 is utilized as default in this study.

    1.3 Initial and boundary conditions

    In the COULWAVE model, for a solitary wave, thecorresponding free surface elevations and velocities are described using conventional Boussinesq theory, in which the characteristic wave length is determined based on 95% solitary wave volume. In addition, the wave crest is located at a user-specified position. In the limited computational domain of a numerical model, the appropriate boundary conditions are needed. In this model, the sponge layer of the absorbing boundary used here absorbs both mass and wave energy, and has been shown to be an excellent absorber for all types of waves, with negligible reflection. Furthermore, in a numerical wave flume, the target structure must be placed far enough from the wavemaker to prevent the reflected wave from touching the left boundary (i.e., incident side). For this reason, the absorbing boundary condition is utilized for the left boundary to avoid re-reflection wave. In the numerical model, the users can switch between the Finite Difference Method (FDM) and the Finite Volume Method (FVM) solvers. The FDM is traditionally used, and is the scheme used by the original numerical code. TheFDM provides considerably high accuracy with relative fast computation. The downside of the FDM is that it is very sensitive to incisive fronts and shocks, making it prone to crashing. The FVM uses a highorder, shock-capturing, approximate Riemann solver for the leading-order flux terms, and is extremely stable and accurate. This scheme adds numerical dissipation for situations with barely resolved shocks, where the wave is not resolved enough for breaking to occur, the shock-capturing properties of the solver preserve the steep front and dissipate energy numerically. Note that the FVM approach requires 50% to 100% or more computational time compared with the FDM solution. Considering its stability, the FVM approach is used here in all overtopping cases.

    2. Laboratory setup and procedure

    2.1 Solitary waves over sloping beaches

    The following variables and the given notations areutilized throughout this paper. x and z represent the abscissa and ordinate in the Cartesian coordinate system, respectively, H0is the initial wave height, η is the free surfaceelevation, h0is the still water depth, h is the local water depth,β is the slope of plane beach, Ruand Rdare the vertically maximum deviation of initial shoreline for the run-up and run-down motions, respectively. Also, we note that in the following investigationsg the gravitational acceleration, R the time history of water shoreline trajectory, ε=H/h0the wave nonlinearity, the superscript “*” the normalized physical quantity. Moreover, the dimensionless variables are defined as

    Table 1(a) Experimental setup and wave condition[1]

    The solitary wave run-up laboratory experiments were carried out in a super wave flume with dimensions of 300 m long, 5.0 m width and 5.2 m deep at the Tainan Hydraulics Laboratory (THL) in National Cheng Kung University. The target solitary wave was generated at one end of the wave flume by a programmable wavemaker. Two impermeable concrete slopes with inclinations of tan-1β=20 and tan-1β=60 were utilized. The correspondingwave conditions and experimental setup are summarized in Table 1(a). Totally, 15 and 80-92 wave gauges were deployed along the wave flume to record the local free surface elevation for the tan-1β=20 and tan-1β= 60, respectively. The reference time is defined as the time at which the incident wave crest passed the reference gauge. The trajectory of shoreline movements were accumulatively recorded by the six five-meter segment wires which were placed along the slope bottom. High coincidence of maximum run-up height is given compared the gauge data with the visual observations determined based on the instantaneous location of the wet-dry intersection between the water layer and the beach slope[1]. It should be noted that all waves employed in the present experiments broke during run-up and run-down courses based on laboratory observations. More measured data and analyses followed by the present numerical studies will be performed in the next section. Further laboratory information and discussions on data can be ascertained in the study by Hsiao et al.[1]

    Fig.1Definition sketc[2h]of a solitary wave interacting with a sloping seawall

    Table 1(b) Experime ntal con[2d]itions of a solitar y wave over a sloping seawall

    2.2 Solitary waves impinging and overtopping an impermeable seawall

    Figure 1 presents the layout of the experimental wave flume and the symbols for the corresponding physical variables used in analyses. This experiment was conducted by Hsiao and Lin[2]. Additionally,cR is the freeboard in Fig.1, which is defined as the vertical distance from the still water level to the seawall crown. The corresponding wave condition and experimental setup are summarized in Table 1(b). The overtopping experiments were conducted in a two-dimensional flume with a scale of 22 m (length), 0.5 m(width), 0.75 m (depth) at THL. The experimental topography was composed of two sections. One is a uniform and impermeable aluminum 1:20 slope starting at 10 m from the wavemaker (i.e., x=10m), the slope surface of which is smooth plexiglas, which significantly reduces friction. The other section is a trapezoidal object made of plexiglas with seaward and landward slopes of 1:4 and 1:1.8, respectively. The seawall model was rigidly mounted on the slope starting at a horizontal distance of 3.6 m from the beach toe (i.e., x=13.6m). Silicone was used to fill the gaps between the seawall, slope, and side glasswall to prevent infiltration. The free surface elevation along the flume during experiments was measured using 9 wave gauges. A reference gauge was placed at 1.1 m in front of the beach slope (i.e., x=8.9m).

    Table 2 Typical comparison of scale effect on maximum run-up height

    Fig.2 Comparison of run-up heights between small and large wave flumes. Present modeling (Δx=0.2m, Δt= 0.02 s, f=0.0025)

    2.3 Scale effect on run-up height

    The property of scale effect is of great interest for various coastal hydrodynamics. This section utilizes our data measured in a large scale tank, the available data obtained in a small flume and present simulated results to discuss the scale effect on run-up height. Table 2 compares run-up data correlated with similar normalized wave nonlinearities of breaking solitary waves climbing uppresent slopes(1:20and 1:60) and the well-known measured data bySynolakis[5]which concern breaking solitary waves travelling upon a 1:19.85 uniform slope in a small wave flume. Moreover, our numerical results for both of corresponding small- and large-scale problems are also incorporated for discussion. Note that the comparisons for gentle and relatively steep slopes are named A and B for the following descriptions, respectively.

    Clearly, favorably good agreements in B are given (also see Fig.2 for the comparis ons with the formulaproposed by Synolakis[5]and our numerical results). The maximum deviation of run-up data is approximately 4.3%. The error may be partly due to the different bottom frictions in the two laboratory environments. According to the comparisons between the experimental data and present simulations, these results imply that the scale effects on the run-up problem could be possibly ignored at least for these two particular setups. However, it is also found that the data listed in Table 2 by Synolakis[5]are somewhat scattering for the issue of scale effect. In particular, the increase of water depth with the same wave nonlinearity the corresponding run-up height increase (i.e., =ε 0.188, ~3.8%) and furthermore, with the increase of wave nonlinearity the corresponding run-up height does not significantly increase (i.e., ε=0.188 and 0.193). Although the discrepancy is not prominent, these phenomena seem to respond to theexperimental observations of Jensen et al.[13]that they attributed the discrepancies of run-up heights due to the nearly the same wave nonlinearities with different water depths to the scale effects correlated with the viscosity and

    Fig3 Normalized wave evolution (a1and b1) and shoreline trajectory (a2and b2) of a non-breaking solitary wave climbing up a 1: 2.75 plane slope

    su rface tension (see Jensen et al.[13], p182).

    Interestingly, for a comparison of Aand B we could not observe the same physical phenomena in A as was noted above (i.e., ε=0.07). It can be seen that the discrepancies of the run-up data correlated with similar normalized wave-nonlinearity in A are almost negligible. It reveals that for the same wave flume in our cases the scale effect could still be reasonably ignored. However, it must be also noted that since no available data performed in a small scale wave tank with a slope of 1:60, direct comparisons with ours are not feasible.

    Overall, the analysesin this section suggest that the scale effect on run-up height in our large scale flume is reasonably negligible compared to the data from Synolakis[5]. In addition, the similar tendency is also observed in the present numerical simulations. However, for a small wave flume the scale effect on run-up height seems more prominent and causes the data somewhat scattered[5,13]. Whether or not, that this discrepancy is related to scale effect needs further investigation.

    3. Numerical model validation

    A series of[7,1n2u]merical tests using the COULWAVE model are presented to interpret the applicability of COULWAVE under one horizontaldime nsion conditions in this section. Both published laboratory works and the present experiments are used to examine the accuracy and stability of the model. Note that only one-layer model, which is identical to the extended Boussinesq equations, is applied in this paper. The numerical setup was the same as that employed in the laboratory experiments. Details of the experimental setups can be found in the original papers. The numerical parameters chosen in the model calculations are revealed in the corresponding figure captions. Besides, we emphasized that our numerical results are not grid dependent.

    Fig.4 Normalized maximum run-up height of anon-breaking solitary wave on a 1: 2.08 sloping beach

    Fig.5 Normalized wave profile snapshot of a breaking solitary wave on a sloping beach

    3.1 Non-breaking solitary waves

    Figure 3 shows the comparisonof present model results and the laboratory data of Zelt[14]. The numerical results of Zelt[14]based on the conventional Boussinesq equations are also introduced in the figures for a discussion. As was observed by Zelt[14]for the two chosen wave conditions, both of which did not break during the run-up process but one of which strongly broke during the run-down course (i.e., =ε 0.2). Apparently, excellent agreements are found between our numerical results and the laboratory data for the time histories of water surface elevations measured respectively in front of the slope (see Figs.3(a1) and 3(b1)) and for the corresponding shoreline trajectories (see Figs.3(a2) and 3(b2)). The discrepancies betw[1e4e]n COUL WAVE and the numerical results of Zelt are clearly observed especially for the case of relatively strong wave nonlinearity. The results indicate that both models can favorably simulate the wave evolution until maximum run-up heights are achieved while COULWAVE can much better describe the wave deformation during the run-down course when the wave breaks (i.e., oscillating tails). In fact, due to the imbalance of wave nonlinearity and frequency dispersivity on the sloping bottom, oscillating tails with dimensionless wavelength khO(2) at the reference gauge are generated. With such dimensionless wavelength, it is not surprising that the phase of oscillating tails based on the conventional Boussinesq equations does show large discrepancy with measureddata.

    Furthermore, it is also worth noting that even our model satisfactorily captures the maximum run-up heights but somewhat underestimates the maximum run-do wn heights compared to the laboratory data for both the cases of small and strong nonlinearities (see Figs.3(a2) and 3(b2)). One reason for causing this is that for a non-breaking wave on a steep slope and the large curvature of the shoreline would lead to larger error[7]. Nevertheless, the overall results are quite favorable compared with the laboratory data.

    Figure 4 further gives a comparison of maximum run-up heights between the present model results and the laboratory data from Li[15]. The waves employed in the experiments did not break for the complete evolution courses as was addressed by Li[15]. Three different bottom friction coefficients are introduced to discuss their influence on maximum run-up heights. Clearly, the maximum run-up heights are insensitive to the chosen bottom frictions for a steep slope and additionally the results indicate that COULWAVE can satisfactorily capture the maximum run-up heights.

    In short, the model results have proved that the moving shoreline scheme and the energy dissipation algorithm employed in COULWAVE are validated to be satisfactory for a solitary wave climbing up a steep slope. Both wave nonlinearity and frequency dispersion can be well captured during the swash cycle by the present numerical model.

    3.2 Breaking solitary waves

    The shoreline behaviors of breaking and nonbreaking solitary waves are quite different. For a breaking wave event, the propagat ing bore after wave breaking will continue to climb up beach slope and the wave reflection phenomenon is relatively insignificant compared to a non-breaking wave[5]. Additionally, as waves propagate over a mild slope the effects of bottom friction play an important role on the swash motion owing to the thin layer of leading bore fronts dominates the run-up tongue and the long course of swash motion[7].

    Figure 5 depicts the comparisons of wave profile snapshots between our numerical simulations and the laboratory data of breaking solitary waves climbing up two different sloping bottoms. The waves employed in Fig.5all broke during the run-up and run-down stages according to the laboratory observations[5,1]. The numerical model of Li and Raichlen[6]based on the nonlinear shallow water-wave equations (in this paper referred to WENO, see Li and Raichlen[6]for details) are also included for a discussion. Note that in the two experimental simulations a bottom friction coefficient f=0.0025 is chosen in the corresponding model calculations. Clearly, Fig.5 reveals that COULWAVE with the employed numerical parameters and the energy dissipation algorithms are validated to be favorable for both gentle slopes in comparison with the laboratory data. The wave life cycle on sloping beaches including propagation, breaking, bore formation, run-up and run-down can be reasonably recreated (see Hsiao et al.[1]for plentiful descriptions of wave evolution properties and physical deformation mechanisms). It can also be seen that for the case of tan-1β= 19.85 the COULWAVE shows a superior capability for the wave deformation modeling compared to the WENO model of Li and Raichlen[6]. The reason is similar to those given in Section 3.1, that is, the WENO model could not well capture the balance and kinematics of wave nonlinearity and phase dispersion effects of a strong nonlinear wave travelling upon a slope owing to the physical limitations of SWE[6].

    Fig.6 Normalized maximum run-up height of a breakingsolitary wave on a sloping beach

    Figure 6 represents the comparison of maximum run-up heights of the present model results with the laboratory data of breaking solitary waves moving on threesloping bottoms. Note th at the bottom friction coefficient for the three simulation cases is 0.0025. Obviously, the results indicate that COULWAVE can excellently simulate the maximum run-up heights for the slopes of 1:20 and 1:60 but slightly overestimate the maximum run-up heights for the slope of 1:15. The discrepancy is partially due to the same coefficients of the bottom friction effects used in two different laboratory environments. Nevertheless, fairly good agreements between the model results and measured data in the Fig.6 imply that the chosen bottom friction coefficient is reasonable for the available experimental environments.

    Fig .7 Normalized shoreline trajectory time series of a breaking solitary wave on a plane slope

    Moreover, Fig.7 compares the computed results with the laboratory data of the time series of shoreline motions for the cases of tan-1β=15, 20 and 60, res pectively. It is noted that the dataof shoreline trajectory(i.e., Fig.7(a), ε=0.3) were not measured by Synolakis[5]but were done by Li and Raichlen[6]with the same experimental conditionsusing a high-speed camera photography. Overall, favorably good agreements are given from the comparisons. We note that from Fig.7 the trajectories of the swash motions on the gentle and relatively gentle slopes are notably different. The maximum shoreline deviation during the run-down processes for the case of tan-1β=60 is nearly close to the original shoreline while for the cases of tan-1β=19.85 and 20 the hydraulic jumps are formed (see Figs.5(a6), 5(b6) and 7). It is because that as the leading bore front after wavebreaking continues to travel upon a gentle beach until the maximum run-up height and the saturation of wave kinetics and potential will be increasingly balanced owing to the incessant effects of gravitational force and the bottom frictions. Therefore as the retreated wave is dragged down the slope and the moving momentum is relatively insufficient compared to the run-up stage. Interestingly, we also observe that for the same slope the arrival time of the maximum run-up are almost the same for different incident wave heights (see Figs.7(b) and 7(c)). More discussions on shoreline movement will be given below.

    4. Results and discussion

    So far, the aforementioned section elucidates that COULWAVE can simulate propagating solitary wave kinematics on mild and steeper sloping beaches with high stability and accuracy. A seriesof numerical investigations was conducted to widespread study the shoreline properties for non-breaking and breaking solitary waves on various beach slopes.

    A numerical slope was constructed with a range of 1≤tan-1β≤60 starting at x=20m in the computational domain. The offshore water depth is fixed at h=1.0m and three wave nonlinearities (ε=0.1, 0.2 and 0.3) were used to determine the characteristics of run-upand run- down heights. Areference gauge was dynamically placed at a distance from the beach toebased on the definition given by Synolakis[5](i.e., Eq.(21)) for mitigating wave reflection. The bottom friction effect on the shoreline movement of a solitary wave is less important for steeper bottoms where wave breaking does not occur. Furthermore, Hsiao et al.[1]experimentally confirmed that the breaking criterion by Grilli et al.[16]is appropriate for a fairly mild slope. Thus, the inclusion of bottom friction effects depends on whether wave breaking occurs during the swashing cycle. The formula given by Grilli et al.[16]is adopted in the numerical calculations

    Fig.8 Normalized (a) run-up and (b) run-down of breaking and non-breaking solitary waves on sloping beaches

    4.1 Run-up/run-down and shoreline properties

    Figure 8 summarizes the calculated results of run-up and run-down of the present numerical experiments. The breaking limitation by Li[15], the run-up law by Synolakis[5], the modified run-up law byLi[15], the co rresponding laboratory data (see the caption of Fig.8) and the empirical formula proposed by Hsiao et al.[1]for generally predicting the maximum run-up heights of breaking solitary waves on sloping beaches (Eq.(22)) are also incorporated into the figures for a comparison.

    Overall, the results in Fig.8 reveal that the maximum run-up/run-down height increases and decreases with the decreasing slope angle for non-breaking solitary waves and breaking waves, respectively. The same conclusion was drawn by Li and Raichlen[6], but only for the run-up. Good agreement can be found between our numerical modeling, the available experimental data, results from analytic approaches, and results calculated using Eq.(22) shown in Fig.8(a). Even though the results are quite promising, there are no available data to validate the run-down results for such a wide range of wave conditions and beach slopes. Nevertheless, it is instructive to observe how the run-down changes with the nonlinearity and slope. Clearly, the run-down height is insensitive to the initial wave conditions for beach slopes below 1:25. This is not surprising because the long course of propagation on a mild slope leads to energy dissipation. Therefore, near saturation on the shoreline is achieved. The run-down height is greatly influenced by wave nonlinearity for the slope of tan-1β<20. In particular, the run-down height increases and decreases with decreasing beach slope until the maximum is reached for non-breaking and breaking events, respectively. This tendency is similar to that for run-up.

    Fig.9 Normalized time series of shoreline motion

    Figure 9 describes the wave nonlinearityeffects on the shoreline trajectories for breaking an d nonbreaking waves on sloping beaches. The difference of run-up/run-down heights between these two slopes increases with increasing wave nonlinearity. In particular, while the rundown for non-breaking wave is approximately proportional to wave nonlinearity, the wave nonlinearity appears to have no influences on run-down for the case of breaking wave. Furthermore, there is a clear time delay of wave arrival time with decreasing wave nonlinearity for the case of breaking wave. In contrast, the wave arrival time of run-up phenomenon is insensitive to wave nonlinearity which is identical to the findings in Fig.7.

    4.2 Impingement and overtopping on a seawall due to solitary waves

    In this subsection, numerical simulations were conducted to evaluate the ability of COULWAVE to simulate fluid motion under complex circumstances, such as a wave impacting and overtopping a seawall or breakwater. To our knowledge, only very few studies that used SWEs to simulate solitary wave overtopping are documented[17]. Therefore, more extensiveand comprehensive c[2]mparisons with experimental data of Hsiao and Lin were made in this subsection. Acceptable agreement with experimental data would validate the use of COULWAVE to predict overtopping volume in coastal areas.

    Fig.10 Time series of free surface elevation due to wave overtopping the impermeable seawall

    Three typical solitary wave cases with different breaking locations were used tosimulate approaching tsunami waves. In the first case, a violent bore rushes landward and subsequently contacts and overtops the seawall. In the second case, a wave directly collapses on the seawall with overtopping flow immediately generated. In the third case, a wave straightforwardly overtops the seawall crown and breaks behind the seawall. These cases were chosen because they exhibit the characteristics of tsunami waves illustrated in tsunami disaster reports.

    T o acqu ire suit abl e compar iso ns, t he numerical flumesetupshouldbeidenticaltotheexperimentalsetup. For computational efficiency, an adequate numerical flume scale must be established. In wave condition of type 1, 150 grids were used to describe one wave length (i.e., Δx=0.015m ) and the slope starts at x=15.0m. 100 and 250 grids were used for simulating one wave length (i.e., Δx=0.027m , Δx=0.02m) in wave types 2 and 3, respectively. Therefore, the numerical reference gaugewas settled at x =13.9m. For describing the bottom friction effects in the laboratory, a value of roughness height Ks=10-4m was adopted throughout all simulations. Additionally, the total physical simulation time was 30 s. The corresponding time step was automatically adjusted during the calculations to satisfy the stability constrain by the advection process, in which the maximum time step was chosen to match the sampling rate of the experimental measurements. The aforementioned computational conditions were utilized in this portion.

    Table 3(a) Total overtopping discharge for type 1 wave condition obtained from simulation and experiments[2]

    Table 3(b)Total overtopping discharge for type 2 wave condition obtained from simulation and experiments

    Table 3(c)Total overtopping discharge for type 3 wave condition obtained from simulation and experiments

    Figure 10 shows the agreement betweenlaboratory data and numerical results of a time series in local free surface elevation along the flume. A sequence of oscillatory waves due to reflection by the seawall at about time t=11.5s can be clearly observed (see Figs.10(a2), 10(b2) and 10(c2)). The reflected wave amplitudes decrease as the freeboard decreases. A reasonable explanation is that in the bore formation and impinging case (i.e.,type 1), the wave breaks earlier compares to the other two cases, the wave breaks in front of the seawall and losses considerable wave energy. Most of the wave energy is thus blocked before the seawall and continuously interacts with the seawall front. In the meantime, a reflected tail wave is produced by the retreating fluid that does not overtop the crown completely. Furthermore, the time series data of free surface elevation shown in Figs.10(c4) and 10(c6) obviously respond the second wave breaking behind the seawall, in that distinct surface variations take place after the peak in the laboratory data. Even ifthis Boussinesq model spreads out an excellent ability to simulate wave breaking and overtopping at most conditions. There still were oscillating shocks in the free surface elevation of Fig.10(b6). This appearance of the oscillating shocks is due to the highly nonlinear processes caused by strong wave collapsing and backwashing. Nevertheless, the crucial and complex circumstances such as wave breaking and overtopping are always difficult and challenging problems for the numerical model with depth-integrated equations to resolve.

    Fig.11 Comparison of calculated overtopping discharges with experimental data for three solitary wave conditions

    In the experiment, the overtopping discharge is the water volume which gets across the seawall crest and reaches the leeward side of the seawall. Hence, after the wave overtops, the water body accumulates between the slope and the seawall. Thus, the experimen tal measurement of the water level is settling a wave gauge at the toe at the leeward side of the seawall. Once the water level is obtained, the overtopping discharge per unit width can be estimated using a basic geometrical relation. In the simulation, the setup must be identical to the experimental conditions for a valid comparison. The simulation results obtained using the Boussinesq equations and the experimental data are compared in Tables 3(a), 3(b), 3(c) and Fig.11. A fairly good match can be observed. The Boussinesq model is considerably accurate in predicting the propagation of a solitary wave overtopping and its computational cost is much lower than those of other numerical models. As a result, the current Boussinesq model can be adequately used to predict the process of solitary waves overtopping the seawall with highly efficiency and acceptable accuracy.

    5. Concluding remarks

    The shoreline motions of non-breaking and breaking solitary waves on plane beaches, and the interactions between tsunami-like solitary waves and an impermeable seawall on a 1:20 sloping beach have been investigated using a modified COULWAVE of Kim et al.[12]. Sets of laboratory data of run-up, shoreline motion and surface elevation of breaking solitary waves climbing up two plane slopes (tan-1β=20 and 60) are reported. In addition, the experimental data of solitary waves overtopping an impermeable seawall is also presented. The scale effect on run-up height of two comparably gentle slopes is discussed. Our analysis suggests that the scale effectcould be reasonably ignored based on limited comparison of experimental data and present numerical results, however, the data from Jensen et al.[13]show inconsistently with our analysis. No decisive conclusions can be drawn and further experimental investigations are needed to verify this issue. In the numerical part for both breaking and non-breaking waves, the model simulation capability is well-validated to be satisfactory by recreating the available experiments and the given laboratory works. The moving boundary scheme and energy dissipation algorithm built in the model are reasonably validated. Numerical simulations have been used to widely study the run-up and run-down heights of non-breaking and breaking solitary waves on various slopes. The maximum run-up/run-down height increases and decreases with increasing slope angle for non-breaking and breaking waves, respectively. The empirical formula proposed by Hsiao et al.[1]was confirmed to be suitable for applications to different nonlinearity tsunami waves on beach slopes with a wide range. The shoreline properties between non-breaking and breaking solitary waves are distinct owing to the wave saturation effects correlated to the energy dissipation processes. The effect of wave non-linearity and freeboard caused by interactions between waves and the seawall was identified. The overtopping discharge monotonically increased with the increasing wave nonlinearity for a given water depth. The results elucidate that the present depth-integrated model can be used as an efficient tool for predicting a wide spectrum of coastal problems.

    Acknowledgements

    This work was financially supported by the National Science Council (Grant NSC 101-2628-E-015-MY3). The authors acknowledge Dr. Hwang K. S. for his laboratory assistance and discussion. The authors gratefully thank the research team at the Tainan Hydraulics Laboratory of National Cheng Kung University for their support in providing experimental data.

    [1] HSIAO Shih-Chun, HSU Tai-Wen and LIN Ting-Chieh et al. Onthe evolution and run-up of breaking solitary waves on a mild sloping beach[J]. Coastal Engineering, 2008, 55(12): 975-988.

    [2] HSIAOShih-Chun, LIN Ting-Chieh. Tsunami-like solitary waves impinging and overtopping an impermeable seawall: Experiment and RANS modeling[J]. Coastal Engineering, 2010, 57(1): 1-18.

    [3] WU Yun-Tai, HSIAO Shih-Chun and HUANG Zhi-Cheng et al. Propagation of solitary waves over a bottom-mounted barrier[J]. Coastal Engineering, 2012, 62(4): 31-47.

    [4]MORI N., TAKAHASHI T. and THE 2011 TOHOKU EARTHQUAKE TSUNAMI JOINT SURVEY GROUP. Nationwide post event survey and analysis of the 2011 Tohoku earthquake tsunami[J]. Coastal Engineering Journal, 2012, 54(1): 125001.

    [5] SYNOLAKISC. E. The run-up of solitary waves[J]. Journal of Fluid Mechanics, 1987, 185: 523-545.

    [6] LI Y., RAICHLEN F. Non-breaking and breaking solitary waves run-up[J]. Journal of Fluid Mechanics, 2002, 456: 259-318.

    [7] LYNETT P. J., WU T.-R. and LIU P. L.-F. Modeling wave runup with depth-integrated equations[J]. Coastal Engineering, 2002, 46(2): 89-107.

    [8] HU K., MINGHAM C. G. and CAUSON D. M. Numerical simulation of wave overtopping of coastal structure using non-linear shallow water equation[J]. Coastal Engineering, 2000, 41(4): 433-465.

    [9]STANSBY P., FENG T. Surf zone wave overtopping a trapezoidal structure: 1-D modelingand PIV comparison[J]. Coastal Engineering, 2004, 51(5-6): 483-500.

    [10]SHAO S., JI C. and GRAPAM D. I. et al. Simulation of wave overtopping by an incompressible SPH model[J]. Coastal Engineering, 2006, 53(9): 723-735.

    [11] LIN Ting-Chieh, HWANG Kao-Shu and HSIAO Shih-Chun et al. An experimental observation of a solitary wave impingement, run-up and overtopping on a seawall[J]. Journal of Hydrodynamics, 2012, 24(1): 76-85.

    [12] KIM D.-H., LYNETT P. J. and SOCOLOFSKY S. A. A depth-integrated model for weakly dispersive, turbulent, and rotational fluid flows[J]. Ocean Modelling, 2009, 27(3-4): 198-214.

    [13]JENSEN A., PEDERSEN G. K. and WOOD D. J. An experimental study of wave run-up at a steep beach[J]. Journal of Fluid Mechanics, 2003, 486: 161-188.

    [14]ZELT J. A. The run-up of non-breaking and breaking solitary waves[J]. Coastal Engineering, 1991, 15(3): 205-246.

    [15] LI Y. Tsunamis: non-breaking and breaking solitary wave run-up[D]. Doctoral Thesis, Pasadena, USA: California Institute of Technology, 2000.

    [16] GRILLI S. T., SVENDSEN I. A. and SUBRAMANYA R. Breaking criteria and characteristics for solitary waves onslopes[J]. Journal of Waterway, Port, Coastal, and Ocean Engineering, 1997, 123(3): 102-112.

    [17]ORSZAGHOVA J., BORTHWICK A. G. L. and TAYLOR P. H. From the paddle tothe beach–A Boussinesq shallow water numerical wave tank based on Madsen and S?rensen’s equations[J]. Journal of Computational Physics, 2012, 231(2): 328-344.

    10.1016/S1001-6058(11)60318-1

    * Biography: TSUNG Wen-Shuo (1987-), Male, Ph. D. Candidate

    HSIAO Shih-Chun,

    E-mail: schsiao@mail.ncku.edu.tw

    精品熟女少妇av免费看| 亚洲欧洲国产日韩| 在线观看美女被高潮喷水网站| 久久国内精品自在自线图片| 午夜久久久久精精品| 亚洲最大成人中文| 亚洲乱码一区二区免费版| 麻豆国产97在线/欧美| 一本久久精品| 亚洲美女搞黄在线观看| 国产精品国产高清国产av| 神马国产精品三级电影在线观看| 久久久久久久久久黄片| 好男人在线观看高清免费视频| 国产精品久久电影中文字幕| 婷婷亚洲欧美| 99久久精品热视频| 精品久久久久久久末码| 内地一区二区视频在线| 久久久久久久久久久免费av| 日本五十路高清| 欧美一区二区国产精品久久精品| 国产中年淑女户外野战色| 日韩大尺度精品在线看网址| 1000部很黄的大片| 一级毛片电影观看 | 久久精品国产亚洲网站| 91狼人影院| 欧美又色又爽又黄视频| 一级二级三级毛片免费看| 狂野欧美白嫩少妇大欣赏| 99精品在免费线老司机午夜| 男人舔奶头视频| 国产一区二区激情短视频| 国产老妇伦熟女老妇高清| 久久久a久久爽久久v久久| 亚洲成人av在线免费| 一区福利在线观看| 性色avwww在线观看| 别揉我奶头 嗯啊视频| 九草在线视频观看| 丰满的人妻完整版| 禁无遮挡网站| 国产激情偷乱视频一区二区| 午夜福利在线在线| 午夜亚洲福利在线播放| 蜜桃久久精品国产亚洲av| 岛国毛片在线播放| avwww免费| 日韩av在线大香蕉| 国产色爽女视频免费观看| 欧美日韩一区二区视频在线观看视频在线 | 亚洲欧美日韩无卡精品| 免费观看人在逋| 夜夜夜夜夜久久久久| 麻豆一二三区av精品| 国内少妇人妻偷人精品xxx网站| 又粗又硬又长又爽又黄的视频 | 99久久九九国产精品国产免费| 日韩中字成人| 少妇猛男粗大的猛烈进出视频 | 99久久精品一区二区三区| 国产精品一区二区性色av| 黑人高潮一二区| 男插女下体视频免费在线播放| 国产又黄又爽又无遮挡在线| 国产高潮美女av| 色哟哟·www| 欧美色视频一区免费| 日韩大尺度精品在线看网址| 人人妻人人澡欧美一区二区| 欧美成人a在线观看| 国产精品久久电影中文字幕| 日本黄大片高清| 亚洲国产欧美在线一区| 麻豆乱淫一区二区| 亚洲最大成人中文| 一本—道久久a久久精品蜜桃钙片 精品乱码久久久久久99久播 | 久久精品国产亚洲av涩爱 | 午夜久久久久精精品| 亚洲中文字幕日韩| 国产精品日韩av在线免费观看| 国产人妻一区二区三区在| 亚洲丝袜综合中文字幕| 美女大奶头视频| av在线播放精品| 一卡2卡三卡四卡精品乱码亚洲| www.av在线官网国产| 国产不卡一卡二| 久久韩国三级中文字幕| 久久精品国产鲁丝片午夜精品| 免费av不卡在线播放| 久久中文看片网| 久久国内精品自在自线图片| 亚洲精华国产精华液的使用体验 | 欧美日本亚洲视频在线播放| avwww免费| 亚洲七黄色美女视频| 国产精品蜜桃在线观看 | 欧美一区二区亚洲| 亚洲久久久久久中文字幕| .国产精品久久| 国产精品一及| 熟女人妻精品中文字幕| 国产伦在线观看视频一区| 久久久久久九九精品二区国产| 日韩中字成人| 免费电影在线观看免费观看| 日韩大尺度精品在线看网址| a级毛色黄片| av黄色大香蕉| 国产黄片美女视频| 少妇被粗大猛烈的视频| 久久综合国产亚洲精品| 狂野欧美白嫩少妇大欣赏| 国内精品久久久久精免费| 永久网站在线| 久久中文看片网| 中国美女看黄片| 亚洲图色成人| 国产成人freesex在线| 22中文网久久字幕| 国产精品av视频在线免费观看| 99在线人妻在线中文字幕| 久久综合国产亚洲精品| 色播亚洲综合网| 亚洲欧美日韩高清专用| 国产精品伦人一区二区| 亚洲在线观看片| 两性午夜刺激爽爽歪歪视频在线观看| 在线天堂最新版资源| 国产精品无大码| 精品久久国产蜜桃| 中国美女看黄片| 欧美成人免费av一区二区三区| 黄色配什么色好看| 精华霜和精华液先用哪个| 偷拍熟女少妇极品色| 特大巨黑吊av在线直播| av.在线天堂| 国产精品三级大全| 欧美成人免费av一区二区三区| 一级黄色大片毛片| 国产一区二区激情短视频| 国产69精品久久久久777片| 白带黄色成豆腐渣| 麻豆国产97在线/欧美| 国产视频内射| 国产精品久久久久久久电影| 噜噜噜噜噜久久久久久91| 晚上一个人看的免费电影| 黄色日韩在线| 国产一区二区三区在线臀色熟女| 午夜久久久久精精品| 22中文网久久字幕| 精品久久久久久久末码| 国产欧美日韩精品一区二区| 最新中文字幕久久久久| 久久午夜亚洲精品久久| 久久人妻av系列| 免费av不卡在线播放| 欧美日本亚洲视频在线播放| 亚洲国产欧洲综合997久久,| 最近视频中文字幕2019在线8| 18禁在线播放成人免费| 亚洲va在线va天堂va国产| 国产成人影院久久av| 禁无遮挡网站| 国产人妻一区二区三区在| 热99re8久久精品国产| 啦啦啦啦在线视频资源| 国产伦在线观看视频一区| 午夜福利在线观看免费完整高清在 | 成人毛片a级毛片在线播放| www.色视频.com| 观看美女的网站| 国产精品一区www在线观看| 中文字幕制服av| 国产精品不卡视频一区二区| 亚洲第一区二区三区不卡| 男女做爰动态图高潮gif福利片| 国产黄色小视频在线观看| 99热全是精品| 国产在线男女| 神马国产精品三级电影在线观看| 欧美在线一区亚洲| 国产精品人妻久久久久久| 少妇丰满av| 亚洲最大成人手机在线| 日韩成人伦理影院| 国产成年人精品一区二区| 色视频www国产| 日本成人三级电影网站| 欧美一区二区亚洲| 亚洲丝袜综合中文字幕| 日日撸夜夜添| 一级av片app| 亚洲国产高清在线一区二区三| 亚洲无线在线观看| 久久久a久久爽久久v久久| 免费看a级黄色片| 国产精品久久视频播放| 国产成人一区二区在线| 亚洲欧洲日产国产| 午夜久久久久精精品| av天堂在线播放| 久久午夜福利片| a级一级毛片免费在线观看| eeuss影院久久| 精品欧美国产一区二区三| 麻豆成人av视频| 日韩欧美三级三区| 91精品国产九色| 身体一侧抽搐| 欧美xxxx黑人xx丫x性爽| 中文字幕av在线有码专区| 亚洲欧美成人综合另类久久久 | 免费人成在线观看视频色| 91av网一区二区| 麻豆成人午夜福利视频| 女同久久另类99精品国产91| 狠狠狠狠99中文字幕| 大又大粗又爽又黄少妇毛片口| 美女大奶头视频| 麻豆成人午夜福利视频| 99视频精品全部免费 在线| 偷拍熟女少妇极品色| 99热这里只有是精品50| 菩萨蛮人人尽说江南好唐韦庄 | 看黄色毛片网站| 99久久精品热视频| 特大巨黑吊av在线直播| 国产av在哪里看| 美女cb高潮喷水在线观看| 亚洲性久久影院| 黑人高潮一二区| 国产精品蜜桃在线观看 | 亚洲精品国产av成人精品| 欧美人与善性xxx| 又粗又硬又长又爽又黄的视频 | 国产三级在线视频| 国产精品美女特级片免费视频播放器| 全区人妻精品视频| av卡一久久| 人人妻人人看人人澡| 精品一区二区免费观看| 一卡2卡三卡四卡精品乱码亚洲| 中文字幕av成人在线电影| 丝袜美腿在线中文| 欧美人与善性xxx| 变态另类丝袜制服| 高清在线视频一区二区三区 | 久久99热这里只有精品18| 久久久久性生活片| 日本黄大片高清| 卡戴珊不雅视频在线播放| 国产在线男女| 国产黄a三级三级三级人| 亚洲不卡免费看| 亚洲真实伦在线观看| 日本欧美国产在线视频| 九九在线视频观看精品| 久久韩国三级中文字幕| 日本熟妇午夜| 偷拍熟女少妇极品色| 啦啦啦观看免费观看视频高清| 国产不卡一卡二| 精品无人区乱码1区二区| 欧美日韩一区二区视频在线观看视频在线 | 欧美色欧美亚洲另类二区| 成人特级黄色片久久久久久久| 一夜夜www| 国产在线男女| 亚洲熟妇中文字幕五十中出| 十八禁国产超污无遮挡网站| 久久精品国产清高在天天线| 国产精品久久电影中文字幕| 黄片wwwwww| 在线免费观看不下载黄p国产| 欧美区成人在线视频| 一夜夜www| 国产午夜精品一二区理论片| 久久久久久国产a免费观看| 丝袜喷水一区| 99久久中文字幕三级久久日本| 极品教师在线视频| 亚洲国产精品sss在线观看| 欧美xxxx性猛交bbbb| 免费观看在线日韩| 一级毛片我不卡| 国产精品乱码一区二三区的特点| or卡值多少钱| 人人妻人人看人人澡| 久久精品夜夜夜夜夜久久蜜豆| videossex国产| 久久久欧美国产精品| 国产爱豆传媒在线观看| 中文亚洲av片在线观看爽| 伦精品一区二区三区| 黄片无遮挡物在线观看| 听说在线观看完整版免费高清| 久久午夜亚洲精品久久| 国产激情偷乱视频一区二区| 只有这里有精品99| 亚洲av第一区精品v没综合| 欧美成人免费av一区二区三区| 国产毛片a区久久久久| 日韩制服骚丝袜av| 久久国内精品自在自线图片| 国产成人午夜福利电影在线观看| 精华霜和精华液先用哪个| 国产精品久久视频播放| 青春草视频在线免费观看| 日产精品乱码卡一卡2卡三| 日韩av在线大香蕉| 午夜久久久久精精品| 22中文网久久字幕| 日本av手机在线免费观看| 精品久久久久久久久久久久久| 精品人妻偷拍中文字幕| 美女内射精品一级片tv| 国产精品麻豆人妻色哟哟久久 | 在线a可以看的网站| 欧美激情在线99| 亚洲无线观看免费| 听说在线观看完整版免费高清| 国产伦精品一区二区三区四那| 亚洲国产欧洲综合997久久,| 青春草视频在线免费观看| 日韩三级伦理在线观看| 天堂网av新在线| 国产精品免费一区二区三区在线| 十八禁国产超污无遮挡网站| 亚洲国产欧美在线一区| 啦啦啦韩国在线观看视频| 亚洲18禁久久av| 男女那种视频在线观看| 国产精品一区二区三区四区久久| 精品久久久噜噜| 国产精品国产三级国产av玫瑰| 国产av一区在线观看免费| 精品国产三级普通话版| 久久精品国产鲁丝片午夜精品| 热99在线观看视频| 精品人妻一区二区三区麻豆| 久久久久久伊人网av| 国内久久婷婷六月综合欲色啪| 精品国产三级普通话版| 又爽又黄a免费视频| 婷婷六月久久综合丁香| 成人亚洲精品av一区二区| 精华霜和精华液先用哪个| 午夜激情福利司机影院| 熟女电影av网| 国产成人91sexporn| 久久精品久久久久久噜噜老黄 | 男女边吃奶边做爰视频| 免费在线观看成人毛片| 美女脱内裤让男人舔精品视频 | 亚洲欧美精品综合久久99| 99久久久亚洲精品蜜臀av| 国产黄a三级三级三级人| 亚洲精品乱码久久久v下载方式| 黄色一级大片看看| 3wmmmm亚洲av在线观看| 亚洲综合色惰| 亚洲精品影视一区二区三区av| 1024手机看黄色片| 免费观看精品视频网站| 国产一区二区三区av在线 | 欧美最新免费一区二区三区| 亚洲欧洲日产国产| 内地一区二区视频在线| 麻豆久久精品国产亚洲av| 91在线精品国自产拍蜜月| 99热精品在线国产| 午夜a级毛片| 美女 人体艺术 gogo| 国产高清不卡午夜福利| 99九九线精品视频在线观看视频| 久久人妻av系列| 国产乱人偷精品视频| 日本与韩国留学比较| 久久午夜亚洲精品久久| 人妻制服诱惑在线中文字幕| 国产国拍精品亚洲av在线观看| 真实男女啪啪啪动态图| 国产亚洲91精品色在线| 黄色欧美视频在线观看| 久久久久久久午夜电影| 国产成人一区二区在线| 女人十人毛片免费观看3o分钟| 婷婷亚洲欧美| 性色avwww在线观看| 国产精品久久久久久久电影| 在线免费十八禁| 22中文网久久字幕| 美女 人体艺术 gogo| 亚洲精品久久久久久婷婷小说 | 九九爱精品视频在线观看| 国产伦精品一区二区三区视频9| 精品人妻一区二区三区麻豆| 亚洲国产欧洲综合997久久,| 小蜜桃在线观看免费完整版高清| 久久综合国产亚洲精品| 免费无遮挡裸体视频| 国产精品久久视频播放| 日本黄大片高清| 成人av在线播放网站| 成人毛片a级毛片在线播放| 国产麻豆成人av免费视频| av天堂在线播放| 性欧美人与动物交配| 久久99蜜桃精品久久| 一级av片app| 久久精品久久久久久噜噜老黄 | 国产熟女欧美一区二区| 亚洲国产精品成人久久小说 | 一卡2卡三卡四卡精品乱码亚洲| 男女做爰动态图高潮gif福利片| 国产在视频线在精品| 亚洲国产欧美在线一区| 18+在线观看网站| 18禁在线播放成人免费| 看十八女毛片水多多多| 99久国产av精品国产电影| 一区二区三区高清视频在线| 韩国av在线不卡| 成人三级黄色视频| 直男gayav资源| 日韩三级伦理在线观看| 人人妻人人看人人澡| 亚洲av一区综合| 黄色配什么色好看| 午夜免费激情av| 日韩欧美 国产精品| 韩国av在线不卡| 哪个播放器可以免费观看大片| 亚洲av电影不卡..在线观看| 爱豆传媒免费全集在线观看| 99久久精品一区二区三区| 美女xxoo啪啪120秒动态图| www.av在线官网国产| 最后的刺客免费高清国语| 寂寞人妻少妇视频99o| 人妻少妇偷人精品九色| 亚洲国产精品sss在线观看| 蜜桃久久精品国产亚洲av| 国产精品国产三级国产av玫瑰| 国产欧美日韩精品一区二区| 天美传媒精品一区二区| 九九在线视频观看精品| 狂野欧美白嫩少妇大欣赏| 精品人妻视频免费看| 内地一区二区视频在线| 长腿黑丝高跟| 91狼人影院| 男女视频在线观看网站免费| av在线蜜桃| 久久精品国产亚洲网站| 国产69精品久久久久777片| 男女啪啪激烈高潮av片| 中出人妻视频一区二区| 我的女老师完整版在线观看| 国产精品野战在线观看| 国产黄色小视频在线观看| 欧美最新免费一区二区三区| 成人av在线播放网站| 亚洲欧美精品专区久久| 高清毛片免费看| 国产精品一区二区三区四区久久| 国产亚洲精品av在线| 免费人成视频x8x8入口观看| 欧美三级亚洲精品| 欧美色欧美亚洲另类二区| 免费观看a级毛片全部| 桃色一区二区三区在线观看| av在线观看视频网站免费| 最近2019中文字幕mv第一页| 边亲边吃奶的免费视频| 熟女人妻精品中文字幕| av在线蜜桃| 美女内射精品一级片tv| 亚洲,欧美,日韩| 亚洲欧洲日产国产| 国产又黄又爽又无遮挡在线| 成人国产麻豆网| 欧美成人免费av一区二区三区| 国产精品无大码| 久久精品久久久久久久性| 亚洲真实伦在线观看| 秋霞在线观看毛片| 亚洲精品日韩在线中文字幕 | 12—13女人毛片做爰片一| 亚洲四区av| 日本黄色视频三级网站网址| 久久久久九九精品影院| 熟女电影av网| 99riav亚洲国产免费| 波野结衣二区三区在线| 亚洲高清免费不卡视频| 国产精品国产三级国产av玫瑰| 高清在线视频一区二区三区 | 国产三级中文精品| 三级毛片av免费| 亚洲真实伦在线观看| 一个人观看的视频www高清免费观看| 性欧美人与动物交配| www日本黄色视频网| 国产av一区在线观看免费| 一进一出抽搐gif免费好疼| 亚洲欧美中文字幕日韩二区| 成人毛片60女人毛片免费| 可以在线观看的亚洲视频| 人人妻人人澡人人爽人人夜夜 | av在线老鸭窝| 国产国拍精品亚洲av在线观看| 国产在线男女| 亚洲av一区综合| 网址你懂的国产日韩在线| 色播亚洲综合网| 中文亚洲av片在线观看爽| 成人av在线播放网站| 久久99精品国语久久久| 成熟少妇高潮喷水视频| 久久99热这里只有精品18| 丝袜喷水一区| 一本—道久久a久久精品蜜桃钙片 精品乱码久久久久久99久播 | 免费看av在线观看网站| 婷婷色综合大香蕉| 欧美极品一区二区三区四区| 18禁黄网站禁片免费观看直播| 亚洲av电影不卡..在线观看| 国产不卡一卡二| 九九爱精品视频在线观看| 最新中文字幕久久久久| 人妻久久中文字幕网| 久久婷婷人人爽人人干人人爱| 晚上一个人看的免费电影| 亚洲成人中文字幕在线播放| 91精品国产九色| 亚洲精品日韩av片在线观看| 国产一级毛片在线| 97人妻精品一区二区三区麻豆| 国产精品一二三区在线看| 亚洲精品成人久久久久久| 国产v大片淫在线免费观看| 91av网一区二区| 国产一级毛片在线| 成人三级黄色视频| 久久草成人影院| 18禁黄网站禁片免费观看直播| 如何舔出高潮| eeuss影院久久| 啦啦啦观看免费观看视频高清| 国产成年人精品一区二区| 午夜精品一区二区三区免费看| 国产精品人妻久久久影院| 久久久久久久亚洲中文字幕| 久久中文看片网| 日本一本二区三区精品| 色吧在线观看| 国产精品麻豆人妻色哟哟久久 | 精品久久久久久久人妻蜜臀av| 国产黄a三级三级三级人| 欧美潮喷喷水| 色综合站精品国产| av在线亚洲专区| 国产精品久久视频播放| 国内精品美女久久久久久| 久久人人爽人人片av| 少妇丰满av| 亚洲在线观看片| 久久久久久久午夜电影| 欧美高清成人免费视频www| 欧美性感艳星| 婷婷精品国产亚洲av| 日韩高清综合在线| 欧美精品一区二区大全| 99热全是精品| 日韩中字成人| 老师上课跳d突然被开到最大视频| 色哟哟·www| 欧美+日韩+精品| av在线播放精品| 岛国毛片在线播放| 精品免费久久久久久久清纯| 在线观看免费视频日本深夜| 99在线视频只有这里精品首页| av在线观看视频网站免费| 亚洲美女搞黄在线观看| 国产亚洲91精品色在线| 在线播放无遮挡| av卡一久久| 国产午夜福利久久久久久| 久久久久网色| 日本免费一区二区三区高清不卡| 亚洲成a人片在线一区二区| 男人舔女人下体高潮全视频| 国产又黄又爽又无遮挡在线| a级毛片免费高清观看在线播放| 国产精品人妻久久久久久| 日本爱情动作片www.在线观看| h日本视频在线播放| 精品欧美国产一区二区三| 日日啪夜夜撸| 欧美另类亚洲清纯唯美| 91麻豆精品激情在线观看国产| av专区在线播放| 免费观看人在逋| 99久久九九国产精品国产免费| 春色校园在线视频观看| 熟女电影av网| 熟女人妻精品中文字幕| 国产真实伦视频高清在线观看|