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

    Review: Recent Development of High-Order-Spectral Method Combined with Computational Fluid Dynamics Method for Wave-Structure Interactions

    2020-07-15 06:29:54YuanZhuangandDechengWan

    Yuan Zhuang and Decheng Wan

    (Computational Marine Hydrodynamics Lab (CMHL), State Key Laboratory of Ocean Engineering, School of Naval Architecture,Ocean and Civil Engineering, Shanghai Jiao Tong University, Shanghai 200240, China)

    Abstract: The present paper reviews the recent developments of a high-order-spectral method (HOS) and the combination with computational fluid dynamics (CFD) method for wave-structure interactions. As the numerical simulations of wave-structure interaction require efficiency and accuracy, as well as the ability in calculating in open sea states, the HOS method has its strength in both generating extreme waves in open seas and fast convergence in simulations, while computational fluid dynamics (CFD) method has its advantages in simulating violent wave-structure interactions. This paper provides the new thoughts for fast and accurate simulations, as well as the future work on innovations in fine fluid field of numerical simulations.

    Keywords: potential-viscous flow; high-order-spectral (HOS) method; computational fluid dynamics (CFD) method

    1 Introduction

    In recent years, much research attention has been devoted to wave-structure interactions in open seas and nonlinear waves. The numerical research on wave-structure interaction take interest in large-scale model and real-sea state simulations. This gives a big challenge on modeling the nonlinear wave propagation in large domain and the application of computational fluid dynamics (CFD). With the development of computer science, the CFD method becomes a popular approach due to its ability to resolve the complicated wave-structure interaction problems such as impact pressure on structures and wave patterns around high-speed ship. The ability of the CFD method to deal with large-scale model and refinement fluid simulations makes it become the main method in the future. However, the cost of propagating open-state waves or freak waves in the CFD method cannot be ignored. As is known, the viscous effect in incident wave is small, therefore employing the potential theory to do the fast wave generation is a good choice in numerical simulations. The traditional potential theories are mature in generating nonlinear waves, but a spectral method called high-order-spectral (HOS) method has salient ability on dealing with long-time wave evolution and large domain propagation in an efficient way. Yet its shortcomings are obvious: it has difficulties in simulating large amplitude motion of structures and complicated wave phenomena such as wave breaking. Therefore, combining the HOS method with the CFD method seems to be the best way to solve complex wave-structure interactions efficiently and accurately. This paper reviews the development of the HOS method and includes the coupling method in recent years. Besides, based on the strength as well as shortcomings of the combined method, this paper discusses the future research work of this combined method.

    2 HOS Method

    The HOS method is a robust numerical method for solving nonlinear wave equations. With the application of pseudo-spectral method and Fast Fourier Transform, the HOS method embodies the feature of efficiency and fast convergence. This method has another name called Dirichlet-Neumann Operator[1]. Its accelerated version[2]is equivalent to the HOS method[3]. Ducrozet et al.[4]gave a comparison and conclusion that HOS models are more efficient than advanced models with finite-difference discretization. Dommermuth and Yue[5]compared spectral method with fully nonlinear mixed-Eulerian-Lagrangian (MEL)[6], and this spectral method accomplishes a more rapid convergence in calculation as well as total energy conservation. With the computational effort increasing only linearly with modes and orders, the spectral method costs less computational time.

    2.1 Theory of the HOS Method

    Dommermuth and Yue[5]gave the numerical description of the HOS method. They extended Zakharov equation[7-8]/model coupling idea[9-11]to steeper waves, and modified this equation and the approach to a higher order for gravity waves. It can simulate the wave nonlinearity into a specified order and use a large number of free Fourier modes in nonlinear simulations. The computational effort increases only proportional to the order and modes. The fluid is considered to be irrotational, homogeneous, incompressible, and inviscid. They considered a velocity potential which satisfies Laplace’s equation. With the pseudo-spectral method, the equation only accounts for the free surface, thus they defined the surface potential asφs(x,t)=φ(x,η(x,t),t). The kinematic and dynamic boundary conditions on the free surface ofφsare[5]

    ηt+xφs·xη-(1+xη·xη)φz(x,η,t)=0

    (1)

    (2)

    They assumed a measure of the wave steepnessεas a small parameter.φandηareO(ε) quantities. For a given orderMinε,φcan be written in a perturbation series based onε[5]. Thusφ(m)is the same quantity withO(εm). They expanded eachφ(m)onz=ηin a Taylor series, then[5]

    (3)

    When expanding Eq. (3)[12-13], they got a sequence of equations ofφ(m):

    Collectingηin the same order, Eq. (3) is a Dirichlet boundary condition. The surface potentialφs(x,0) andη(x,0) are known in

    (4)

    wherem=1,φs=φ(m)(x,0,t). Eq. (4) is the final formulation of the HOS method. While in a mode-coupling approach, except Dirichlet free-surface condition (4), Eq. (3) can be solved by using eigenfunction[5].

    The eigenfunctionφnhas analytical solutions for simple geometries. For example, the 2π-periodic boundary conditions can be represented as[5]

    For constant finite deep water, with depthh[5],

    The study of Dommermuth and Yue[5]is essential for developing the the HOS method. With its high efficiency and accuracy, it opens up a path in simulating nonlinear gravity waves for succeeding researchers. With the existence of the exponential profile, the HOS method can provide rapid convergence and consume less time. The boundary-value problems are solved rapidly through Fast Fourier Transform, which connects spatial derivatives in spectral representation and nonlinear products in physical space.

    West et al.[14]presented using the HOS method to study the free and bound waves on nonlinear ocean surface. The descriptions of the HOS method is similar to Dommermuth and Yue[5]. However, there is a little difference in dealing with the vertical velocity in West[14].They first assumed a power series inε:

    (5)

    (6)

    In this way, all the terms under boundary condition which contain the vertical velocity are consistent with the orderingε, making the comparison with analytical data easier.

    Although the advantages of the HOS method are obvious, some errors exist due to its features. Because the velocity potential is single value in the HOS method, it cannot simulate overturn phenomenon in waves. Besides, the HOS method has a limit in wave steepness, and it will fail when wave steepness is very large.

    The first error is the truncation in the number of modesNand orderM. For the singularity in the analytic continuation of velocity potential in the HOS method, the convergence breaks when the wave steepness is large. The second is amplification of round-off errors. The error at any order in high wavenumber modes will be amplified. The unstable growth of this error in many nonlinear free-surface simulations will occur when high wave number is used[15-16]. Dommermuth and Yue[5]used Longuet-Higgins and Cokelet’s methods to make a low-pass filter to eliminate such high-wavenumber instabilities.

    The aliasing errors[17-18]happen because the nonlinear products in free-surface boundary conditions are performed in physical space instead of in spectral space. West et al.[14]removed this error by changing the numbers of mesh pointsNx,Nyin the physical plane to satisfy the condition

    Nx>(M+1)kmax,Ny>(M+1)lmax

    De-aliased computation can be performed by extending spectra in zero padding[19-20]. The zero-padding applies half-rule and are defined on a number of points:

    2.2 Development and Application of the HOS Method

    Liu et al.[21]extended the HOS method to study nonlinear wave-body interactions. They considered the submerged body as a dipole distributed on the free surface and a source affected on the body surface. They showed the ability and accuracy of the spectral method in simulating submerged body. However, they also pointed out the method is not available when the body submergence is small and in wave breaking case.

    Unlike many numerical models[8,10,22]which have difficulties in simulating deep-water dispersive waves due to the complexity of the nonlinear terms, the HOS method is more suitable for strong nonlinearity terms. However, it has difficulty in simulating the progress of nonlinear wave which may give spurious high-frequency waves. It can be eliminated by either giving sufficient time when initialized with linear wave or changing initialization procedures. Dommermuth[23]developed an adjustment for simulating the turbulent free-surface flows. In case of turbulent free surface flows, Dommermuth[24]also developed an adjustment that allows nonlinear free-surface simulations to be initialized with linear solutions. By introducing the vertical component velocityφz|z=η=Ws(x,t), and letting the first order vertical velocity denote the leading-order, the free-surface Eqs. (1)-(2) can be expressed as[24]

    ηt-W(1)=-W(1)-xφs·xη+
    (1+xη·xη)Ws

    Then the nonlinear terms were all collected on the right-hand side. They applied a termTawhich is the duration of the transition period. They replaced right-hand terms into terms withTa(namedFandG), which smoothed the linear initial conditions to nonlinear computations. As the partial derivatives of the nonlinear termsFandGare equal to zero at any order, the procedure reduces the imbalance under the initial conditions. Dommermuth[23,25]also used atmosphere pressure terms to balance the dynamic free-surface boundary conditions. It needs to be mentioned that the work of Refs. [23] and [25] requires the adjustment of dynamic boundary condition, while that of Ref. [24] requires adjustment of both dynamic and kinematic boundary condition.

    The solution of the initial condition and less computational time of the HOS method make extreme waves with long time generation in open seas more available. Brandini and Grilli[26]used the HOS method to do the first stage of growth of wave modulations and then applied the results as the initial conditions for extreme waves generation. It is worth noting that Brandini and Grilli presented the 2D HOS method as initial condition in transversal and longitudinal direction and it showed typical 3D wave modulations. Tanaka[27]employed the HOS method as the basis numerical scheme, and in order to study the energy spectrum or time evolution of skewness and kurtosis of wave elevation, they introduced a complex amplitude function[8]. The complex amplitude functionb(k,t)is defined as[8]

    φsandηare described withb(k,t):

    The complex amplitudeb(k) is given by the value ofkin discretized space, and for eachk,

    The energy densityEcan be given according tobkon the wave field:

    whereEcan also be given in terms of directional spectrumΦ(ω,θ) as

    (7)

    ΔSk=Δkx×Δky

    (8)

    Eqs. (7)-(8) gives the relation between |bk| andΦ(ω,θ):

    (9)

    Through Eq. (9), the directional spectrumΦ(ω,θ) can be determined by calculating the norm |bk(t=0)|, and the phase ofbk(t=0) is determined by a random number in [0,2π].Φ(ω,θ) can be expressed with JONSWAP spectrum.

    Ducrozet et al.[28]followed the initializing condition of Dommermuth[23]and simulated a forced focusing event wave based on the method of Tanaka[27], and then a natural freak waves during a long time simulation were considered as well, as shown in Fig.1. Their study gives a definite description on generating freak wave in an open-sea state, discusses the selection of parameters in directional spreading, and gives a proper solution to it. This encourages the exploration on the freak wave mechanics and the developments of the HOS method. Li et al.[29]computed the sloshing wave in 3D tank which shows the efficiency of the HOS method in simulating focusing wave.

    Fig.1 3-D surface elevation at t=20.3Tp[28](Reprinted with the permission from Ref.[28] copyright ?Author(s) 2007. This work is Licensed under a Great Commons License)

    Zhao et al.[30-32]applied an improved Adams-Bashforth-Moulton (ABM) numerical integration instead of Runge-Kutta forth (RK4) scheme to reduce simulation time and studied focusing wave simulation in experimental scale and real scale. The improved ABM scheme can reach a truncation error ofO(Δt6), while RK4 scheme can only reachO(Δt5). However, as mentioned in Dommermuth and Yue’s work[5], they employed RK4 instead of ABM because RK4 can provide a lower global truncation error and a larger stability region. The improved ABM only increases the precise order, and they did not mention the global truncation error and numerical stability in their work. Therefore, it needs to be discussed in the future which numerical integration method has wider application scope.

    Sergeeva and Slunyaev[33]studied the characteristics of rogue waves by the HOS method. They produced space-time wave data sheets using the HOS method, and observed the detailed picture of individual rogue waves. They simulated 10 km field with 20 min sheets of data with surface elevation, fluid velocity, and so on. Fig.2 illustrates rogue waves and rogue events, where markers indicate rogue waves and rectangles assemble them into rogue events. Rogue waves are difficult to observe for its long-time duration, short time appearance, and high crests with shallow troughs. Their study applies the HOS method to give the probability of appearance and lifetime of rogue waves, showing the possibility to study the rogue waves and wave-structure interactions in rogue waves. The application of the HOS method makes the study of rogue waves possible.

    Fig. 2 Spatio-temporal sheet of the surface elevation[33](Reprinted with the permission from Ref.[33] ?Author(s) 2013.CC Attribution 3.0 License)

    Xiao et al.[34]studied the occurrence and dynamics of rogue waves employing the HOS method in three dimensions. They compared the HOS results to nonlinear Schrodinger equations, giving a quasi-stationary evolution in long-time prediction, which cannot be predicted by the nonlinear Schrodinger equations. They also introduced the definition of rogue wave occurrence which related to the occurrence probability, and discussed different situations of narrow spread seas. Along with the study mentioned above, the research on probability of rogue wave occurrence provides valuable reference for studying the structures which were emerged and changed in the rogue waves.

    Those studies focused on unbounded domains of open-sea states, with periodic boundary conditions applied. However, unbounded domain means that the wave propagates without wave generation, absorption, and beach reflections, which makes experimental validations more difficult. Bonnefoy et al.[35]developed a non-periodic HOS technique and presented a linear wave maker motion in the numerical wave tank. They considered an additional potential[36]added on the original potential to make sure the potential satisfied the wave maker condition. By numerically solving the additional potential also in time domain, they extended the idea of Agnon and Bingham[36], and adjusted a linear piston wavemaker to a generic wavemaker, as shown in Fig.3. They presented 2D cases of focusing wave packet with extending second order initial decomposition and target phase signal iterative corrections.

    Fig.3 Extended basin[35](Reprinted with the permission from Ref.[35] Copyright ?2014 ISOPE)

    Li[37]employed the high-order-spectral tank (HOST) derived from Bonnefoy et al.[35]to simulate 2-D focusing wave and made comparison with the experimental results. Li and Liu[38]also applied the HOST method to generate directional focusing waves. Fig.4 shows the comparison between experimental results and 2-D focusing wave generated by HOST. The surface elevation of 3-D focusing wave generation is displayed in Fig.5. They studied the properties of the focusing wave and gave conclusions that shift of focusing point and maximum crest decreases when spreading parameter decreases in 3-D case. In addition, Li et al.[39]used the HOST method to study the influence of current on the focusing wave. The potential velocity are expressed as the uniform current velocity added with wave potential velocity in their study. The numerical model revealed strong nonlinear interaction in the wave-current interaction.

    Fig.4 Comparison between experimental data and numerical results[38](Reprinted with the permission from Ref.[38] ?2015 Chinese Ocean Engineering Society and Springer-Verlag Berlin Heidelberg)

    Based on the previous works on HOST[35]and its successful results compared with experiments, Ducrozet et al.[40]developed this method into second-order wave maker motion. They combined the model developed by Bonnefoy et al.[13,15]with the HOS method wave tank, and then extended HOST to second order. Fig. 6 shows the calculation process of this new HOST. A 2D focusing wave packet embedded in irregular wave fields was considered and showed apparent improvement in results compared with those in linear HOST.

    Fig.5 Surface elevation of 3-D focusing waves at selected time[38](Reprinted with the permission from Ref.[38] ?2015 Chinese Ocean Engineering Society and Springer-Verlag Berlin Heidelberg)

    Fig. 6 The 2nd HOST calculation setup (in Runge Kutta)[40](Reprinted with the permission from Ref.[40] Copyright ?by the International Society of offshore and Polar Engineers)

    Bonnefoy et al.[41]applied both second-order wave maker motion and fully nonlinear numerical wave tanks and gave the validation results to experimental complex sea-states in 3-D.

    Those improvements and developments of the HOS method lays the foundation for more application of the HOS method. In order to make the method more conveniently applied in many fields, Ecole Centrale Nantes, LHEEA Lab developed an open-source HOS model[42]named HOS-Ocean. The software is available on the Github (https://github.com/LHEEA/ HOS-ocean/wiki). The HOS method is realized in a numerical simulation way, and provides possibility to couple the CFD software.

    Assuming the simulation based on a rectangular fluid domain with periodic boundary condition, the velocity potential of surface boundary condition can be written as Eqs. (5)-(6), wherekn=n(2π)/Lxis the wave numbers. Fig. 7 presents the brief description of code working in HOS-Ocean according to Ref.[42]. For the wave spectra used in HOS-Ocean is only JONSWAP spectrum, Song et al.[43]made an attempt to add two new wave spectrum P-M spectrum and ITTC 2-parameter spectrum in it.

    Meanwhile, Ducrozet et al.[44]developed the open-source numerical wave tank based on HOS (HOS-NWT). The HOS-NWT is also available on Github (https://github.com/LHEEA/HOS-NWT/ wiki). They extended wave maker modeling into three-order generation.

    Fig.7 Procedure of time-stepping in HOS-Ocean according to Ref.[42]

    Seiffert and Ducrozet[45]included the wave breaking mechanism into HOS-NWT. As the HOS method assumes free-surface as single-valued, they approximated the broken surface as a single value. They suggested a threshold of water particle velocity to crest velocity between a broken and an unbroken wave in HOS solvers. Seiffert et al.[46]implemented a wave breaking onset criteria in HOS-NWT, which predicts good wave-breaking onset time and location compared with the experiments. Seiffert and Ducrozet[47]validated this wave-breaking mechanism through energy dissipation by adding an eddy viscosity parameter under boundary conditions. This kind of wave breaking onset can be incorporated with other nonlinear potential solver to reduce the risk of numerical instabilities, thus increasing the scope of application.

    The development of the HOS method shows a perspective in ocean and ship engineering. The ability of generating large-domain wave field with fast convergence and accuracy gives researchers the chance to study rogue waves. The difficulties in studying rogue waves are due to its transient time appearance and long-time with large space simulation, which requires huge investment in resources. With the help of the HOS method, the parameters of rogue waves and the position of its occurrence can be found out. This lays the foundation for simulating severe conditions in ocean and sail engineering.

    The applications and extensions of the HOS method are not limited to wave simulations. Many physical mechanisms are included by employing and validating the the HOS method, which include nonlinear energy transfers[48], bi-modal seas[49], modulation instabilities[50-51], and so on (cf.[42]). There are also some developments based on the HOS method, including wind forcing effects[52], jet current effects[53], variable bathymetry[54-58], and so on.

    3 Combined Method of the HOS and CFD Method

    In order to simulate the wave-structure interaction in a more realistic and flexible way, the CFD method is carried out. This approach functions extremely well in resolving complex surface problem, large amplitude motion of structures, and so on. However, even though the development of computational science is rapid, the time duration and cost of the CFD method is still nonnegligible. In fact, when wave-structure interaction is considered, wave loads extending on the structures can be focused even more. Besides, the viscous effects are not obvious in wave simulation, and long-time wave simulation in the CFD method may cost a long time and cause numerical dissipation. Therefore, combining the potential and viscous flow together can improve the simulation in an efficient and realistic way.

    Li[59]gave a view of two groups which divides the combined potential and viscous flow: Domain Decomposition and Function Decomposition, the schematic diagram is shown in Fig.8. He mentioned Domain Decomposition sets viscous domain and potential domain apart, while Function Decomposition sets viscous part and potential part in viscous domain.

    3.1 Function Decomposition

    Many studies have been conducted[60-64]to solve a wave-structure interaction employing function decomposition approaches. These studies applied BEM as potential solver and RANS or Lattice Boltzmann Method (LBM) as viscous solver. Later on, a new approach was carried out focused on solving diffracted flow only, which is named Spectral Wave Explicit Navier-Stokes Equations (SWENSE). It is derived from a frame of potential theory[65-66]and splits all unknown variables into incident terms and diffracted terms[67]. SWENSE modifies the variables in RANS equations as follows[68]:

    whereα∈{1,2}, the unknown variablesu,p, andhstands for Cartesian components of velocity, pressure, and free-surface elevation, respectively. The subscripts I and D represent incident and diffracted variables, respectively.

    Fig.8 View of two groups of combined potential and viscous flow according to Ref.[59]

    According to Navier-Stokes equations[59],

    uα=0

    and Euler equations[59],

    The transport equations are obtained by substituting the Euler equations from the NS equations[59],

    With those equations, the incident variables were explicitly computed, thus named the solver SWENSE Equations. The incident nonlinear regular wave was obtained by a stream function model[69]while the nonlinear wave trains were derived from a spectral formulation combined with BEM[70]first, and then from the HOS method. The SWENSE method using the spectral formulation combined with BEM was validated both in 2D cases[71-72]and 3D cases[73].

    Luquet et al.[74]employed the HOS method to play the role as incident wave model. Compared with the previous potential theories, the HOS method makes the simulation of arbitrary complex sea-state or freak waves more convenient and efficient. They used this method to simulate a Tension Leg Platform (TLP) in regular and irregular waves.

    Other studies such as ship maneuvering with 6DOF, CALM buoy, and ship motion in 2DOF have been done in SWENSE model[75-78]. A validation of ship motion in SWENSE in irregular head wave was carried out (Fig.9).

    However, those original SWENSE methods are based on one-phase CFD solver, which solve the free surface through the boundary condition. Li et al.[79-80]extended the SWENSE model into two-phase SWENSE model. They modified Euler momentum equation by introducing an incident pressure to balance the governing equations both in water and air phases.

    The SWENSE model divides the total field into two parts: incident field and diffraction field. In this way, the incident field was solved by potential theory, thus the computational grids can be a coarse mesh while only the grids near diffraction field need to be refined. It not only improves efficiency when using the HOS method as incident model, but also decreases the cost when simulating wave structure interactions due to less grids in far field.

    Gatin et al.[81-82]combined HOS with CFD RANS equation to solve 6DOF ships in extreme waves,as shown in Fig.10 and Fig.11. The combination was performed based on open source CFD software OpenFOAM[83]. They obtained free surface elevation and velocity from the HOS method, thus the velocity field is presented as[81-82]

    wherez'=zd/(d+η)+d(d/(d+η)-1).

    Fig.9 Comparison of ship motion with experimental results[77](Reprinted with the permission from Ref.[77] Copyright ? 2009 ISOPE)

    They coupled the HOS method with the the CFD method using a two-phase SWENSE method[84-86]in level set method[87]. This two-phase SWENSE method is different from that in Li[79-80]. It deals with incident wave as a source term in CFD meshes. Therefore, the incident domain in CFD mesh still needs to be refined to keep accurate.

    Fig.10 Full-scale KCS in extreme waves[81](Reprinted with permission from Ref.[81])

    Fig.11 Extreme wave encountering the barge[82](Reprinted with permission from Ref.[82])

    3.2 Domain Decomposition

    The domain decomposition sets computational domains apart into viscous flow sub-domain and potential flow sub-domain. In general, with advantages of solving complex free surface and viscous effect around structure, the CFD method takes the inner domain and potential flow method takes the outer domain. The two methods solve separately in each sub-domain and communicate through the common boundary. Therefore, the domain decomposition method can be two-way coupling or one-way coupling.

    Before the emergence of the HOS method, either two-way coupling[88-91]or one-way coupling[92-94]employed BEM method in potential flow sub-domain. Kim et al.[95]applied an overlay zone to do the communication between CFD and Euler wave, which is named Euler-Overlay Method (EOM). They combined EOM with mooring and riser model to do the simulation of Vortex-Induced Motion (VIM). Paulsen et al.[96]combined Oceanwave3D[97]with the CFD method and gave validations under different sea bed conditions.

    Choi et al.[98]combined the HOS method with OpenFOAM (foamStar) through a HOS wrapper program named Grid2Grid[99]which can be obtained from GitHub (https://github.com/LHEEA/ Grid2Grid). Grid2Grid reconstructs HOS wave results and interpolates at demanded time and position in viscous flow field. The frame of program of Grid2Grid is shown in Fig.12. Choi et al.[98]tested the new solver in empty wave tank and validated the stable and accuracy of this solver.

    Zhuang et al.[100]developed the HOS method into an in-house CFD solver naoe-FOAM-SJTU[101-108]. They applied relaxation zone[109]as the common boundary to communicate information, as shown in Fig.13. They tested the ability and accuracy of the new combined solver in empty wave tank, and gave good results. Song et al.[110]also applied this new solver to simulate freak waves in the empty tank. The CFD zone can obtain wave elevation and velocity from the HOS method in arbitrary time spot and positions, thus the new combined solver reduces the simulation cost in the CFD method.

    Fig.12 Frame of program in Grid2Grid according to Ref.[99]

    Fig.13 Sketch of the HOS method combined with naoe-Foam-SJTU[100](Reprinted with permission from Ref.[100])

    The CFD solver naoe-FOAM-SJTU is based on OpenFOAM and develops many modules and functions in recent years, such as 6DOF module[111], overset grids[112-114], internal flow coupled with external flow[115-117], and so on. Combined with the HOS method, the new solver has advantages on solving structures in open-sea states in an efficient way. A number of simulations have been carried out[118]by the new combined solver, presenting the ability to simulate the cylindrical structure in multi-directional sea-state (Fig.14) as well as focusing wave (Fig.15).

    Fig.14 Wave pattern of a cylinder in directional wave in HOS-CFD domain (left) and in CFD domain (right)[118](Reprinted with permission from Ref.[118])

    Fig.15 Focusing wave passing through a cylinder in CFD domain

    In general, a ship sailing in the seas is assumed to be sailing in the encountered added current. But a ship sailing in a large domain constitutes realistic ship maneuvering. With the help of overset grids and the HOS method, ship can sail in open seas without large computational domain requirement, as shown in Fig.16. KCS ship model in model scale is chosen as the ship model and the Lpp of the ship is 6.070 2 m. The ship speed is 2.017 m/s, and the wave is first order Stokes wave inλ/L=1.15. Through the results of heave motion and pitch motion as shown in Ref. [118], the coupled method is closer to experimental results than pure CFD method.

    Quinn[119]combined HOS-NWT with OpenFOAM to study the breaking waves. The HOS-NWT replicates a 2D waves conducted in experiments with a 3D long-crest irregular wave and the results were analyzed with different breaking onset criteria. However, the breaking process was not visible in the coupled method. Quinn pointed out that this may due to the diminishing of crest value in HOS-NWT, which provided flows in OpenFOAM. Meanwhile, Quinn also simulated freak and rouge waves in HOS-NWT, and gave a conclusion that with the increase of the wave steepness, the discrepancy between HOS results and experimental results becomes large.

    Fig.16 A ship maneuvering in the regular wave ((a) and (b) show the ship sailing in HOS domain; (c)-(f) show the details of the ship in a period in CFD domain)[118](Reprinted with permission from Ref.[118])

    The review of application with the coupled method gives a huge possibility in wave-structure interactions. Nowadays, with the help of overset grids, the large amplitude of structures is achieved in the CFD method. The coupled method solved a big problem which the CFD method cannot solve at this moment. The HOS method gives the input information of large-scale or large domain of wave or sea states, the CFD method can simulate any situations in ocean and naval engineering. The given results of the coupled method are not enough though. Many research such as large-scale simulations, fine wave fields, and ship sailing in freak waves need to be done.

    4 Conclusions

    The HOS method has been developed in the 1980s. With the feature of efficient wave generation and fast numerical convergence, an increasing number of applications using the HOS method have been carried out. The HOS method is popular in recent years for its ability in simulating large-domain wave generation and long-time evolution. From the very beginning, the HOS method was applied for studying nonlinear gravity wave and wave-wave interaction. After that, the attempt in studying wave-body interaction in the HOS method was carried out. The extreme wave and long-time simulation was completed from 2D to 3D, along with the observation of occurrence of freak wave in open-sea state. Then the application of the HOS method extended to bounded boundary condition which is wave tank simulation. With the development of computer science, the open source software of using the HOS method is well developed. Although many studies have been conducted to eliminate numerical error or employed assumptions to solve the shortfalls, development of the HOS method should not be stopped.

    The combined method of the HOS and CFD method is a new theory in numerical simulations. The idea of combining the HOS with CFD method was put forward in recent decades. The development and application of the integrated method is all about fixed structure to test the accuracy of the combination. Some studies have been done to simulate ships in freak waves, showing the ability of the combined method in dealing with large motion of structure in freak waves.

    Although the combined method of the HOS and CFD method only rises in recent years, its achievement is a progress in numerical simulation. The ability and advantages of the combined method are listed as follows:

    1)Because the HOS method can give incident wave information in arbitrary space and time spot, the CFD method can choose required position and time spot in wave propagation. For example, the HOS method can be employed to do a real-sea state generation. After observing the appearance of freak wave, the CFD method can be used to receive the freak wave signal to simulate the freak wave-structure interactions. It is well known that the capture of freak wave needs thousands of seconds of simulation, which is impossible for the CFD method, while the combined method provides the possibility to do the freak wave structure interactions.

    2)The large-scale model simulation in the HOS method coupled with the CFD method can also be achieved. The HOS method has the advantage in generating wave in large field efficiently, thus it can be an economic way to simulate large-scale model, especially in two-way coupling. The domain of the CFD method can be minimized around the large-scale model, and the wave field outside can be solved by the HOS method. Therefore, it will save lots of time and resources compared with pure CFD method.

    3)The combined method has the ability and advantages in simulating fine wave field. When the simulations are based on focusing wave or fully developed nonlinear wave, the CFD method needs to spend lots of time to make sure the nonlinear wave or focusing wave are fully developed. The viscous effects in the CFD method may diminish the wave amplitude during the wave propagation. In this way, it costs splendid time and resources when the fine wave field is considered. With the help of HOS wave field which is a fully developed wave field, the time during wave propagation can be ignored in the CFD method. Therefore, the meshes in CFD can be as fine as possible to achieve the accuracy as required.

    As stated above, the combined method has great ability in simulating some kinds of wave-structure interactions, which include ships or platforms in freak waves that may encounter impact pressure, breaking wave, and air gap phenomenon; large-scale model simulations; realistic wave-structure interactions such as real-sea state simulation, and so on.

    5 Future Works

    Some research perspectives and corresponding solutions are given below.

    1)There are two kinds of coupled method of the HOS and CFD method: one is function decomposition of two-phase SWENSE[59,79-80], the other is SWENSE in Naval hydropack[81-82,84-86]and domain decomposition[98,100]. The first two-phase SWENSE model is limited to fixed structure up till now, therefore a 6DOF module and dynamic mesh needs to be applied. These modules can be implemented through OpenFOAM. The second SWENSE and the domain decomposition have the 6DOF module and dynamic/overset grids, but they have to refine the mesh in incident wave domain.

    2)For coupled method of the HOS and CFD method in domain decomposition, it is one-way coupling at present. It does not need extra iterations between potential and viscous flow. However, it requires larger domain to realize complete radiation flow of wave-interaction influence than that in two-way coupling. Considering the ship maneuvering in the coupled method, the one-way coupling will end with misplaced wave patterns in the view of whole domain. Nevertheless, the two-way coupling of the HOS and CFD method needs to change the formulation of the HOS method and improvement of the common boundary condition, thus it is still under consideration.

    3)The simulations based on the HOS method[28,42,44]are conducted on single-processor, while simulations in the CFD method often need multiple processors. This leads to separate calculation in the HOS method and the CFD method. In domain decomposition, the CFD calculation in multiple processors will duplicate the results of the HOS method. When the internal memory of HOS results is large, the duplication in parallel running may cause a huge work in CFD simulation. The HOS method has no apparent solutions on parallel running at present, therefore the CFD method can solve this problem by sparing an extra processor for the HOS method. This will prevent the duplication of HOS results files.

    4)The studies of breaking waves are often included in wave-structure interactions. With the assumption of single value of breaking waves[45-47], the HOS method has the ability to predict breaking wave onset criteria. However, the assumption is not realistic and cannot be obviously visible in CFD domain. Therefore, the study of breaking phenomena needs to be done in viscous flow in the future. The CFD method can capture arbitrary time spot of HOS results, thus the thought is to map the breaking time spot of HOS to CFD domain, thus the breaking in waves happens under viscous flow solution. It requires further tests, because the incident information of the HOS method is absent, and the breaking wave in CFD domain may collapse or propagate in a wrong way.

    5)There are some numerical instabilities with the combined HOS and CFD method. As all the variables are solved on free surface in the HOS method, the interaction with the CFD method will cause instability in sharp free surface. For function decomposition of two-phase SWENSE[59,79-80], it applies Ghost Flow Method (GFM)[120-122]to keep discontinuity in the CFD method. For domain decomposition of the combined method, the accumulated velocity on sharp free surface can be solved by an adding inversed velocity.

    9热在线视频观看99| 亚洲国产成人一精品久久久| 亚洲精品456在线播放app| 色婷婷久久久亚洲欧美| 久久 成人 亚洲| 免费观看性生交大片5| 亚洲内射少妇av| 一级毛片电影观看| 欧美精品av麻豆av| 精品一区在线观看国产| 色哟哟·www| 亚洲成国产人片在线观看| 久久久久精品久久久久真实原创| 国产成人精品无人区| 女性被躁到高潮视频| 亚洲欧美一区二区三区国产| 国产 精品1| 天天躁夜夜躁狠狠躁躁| 999精品在线视频| 激情五月婷婷亚洲| 日本色播在线视频| 国产熟女欧美一区二区| 欧美 日韩 精品 国产| 成人漫画全彩无遮挡| 各种免费的搞黄视频| 国产精品免费大片| 黑人猛操日本美女一级片| 黑人高潮一二区| 在线观看免费日韩欧美大片| 国产女主播在线喷水免费视频网站| 国产精品.久久久| 考比视频在线观看| 久久午夜福利片| 日本av免费视频播放| 亚洲久久久国产精品| 亚洲精品久久成人aⅴ小说| 国产xxxxx性猛交| 夜夜爽夜夜爽视频| h视频一区二区三区| 久久午夜福利片| 久久精品久久精品一区二区三区| 久久午夜综合久久蜜桃| 国产福利在线免费观看视频| 久久狼人影院| 老熟女久久久| 丝袜在线中文字幕| 曰老女人黄片| 国产精品麻豆人妻色哟哟久久| 欧美日韩国产mv在线观看视频| 国产一区亚洲一区在线观看| 亚洲一级一片aⅴ在线观看| 亚洲精品aⅴ在线观看| 国产精品国产三级专区第一集| 这个男人来自地球电影免费观看 | 女人精品久久久久毛片| 插逼视频在线观看| 国产男女超爽视频在线观看| 夫妻性生交免费视频一级片| 成年人免费黄色播放视频| 午夜福利,免费看| 久久国产精品大桥未久av| 老司机影院毛片| av网站免费在线观看视频| 亚洲第一区二区三区不卡| 男女免费视频国产| 不卡视频在线观看欧美| 日韩,欧美,国产一区二区三区| 久久久久久久久久成人| 久久毛片免费看一区二区三区| 乱人伦中国视频| 欧美成人精品欧美一级黄| 欧美日韩国产mv在线观看视频| 水蜜桃什么品种好| 王馨瑶露胸无遮挡在线观看| 最黄视频免费看| 国产成人精品无人区| 日本免费在线观看一区| 国产免费福利视频在线观看| 在线观看免费视频网站a站| 亚洲欧美日韩另类电影网站| 国产成人精品婷婷| 最近最新中文字幕免费大全7| 国产精品久久久久久av不卡| 丝袜脚勾引网站| 一个人免费看片子| 免费看不卡的av| 国产精品人妻久久久久久| 伦理电影大哥的女人| 欧美日韩综合久久久久久| 男女下面插进去视频免费观看 | 在现免费观看毛片| 少妇人妻精品综合一区二区| 久久久久久久大尺度免费视频| 精品亚洲成a人片在线观看| 午夜视频国产福利| 久久久久久久久久久免费av| 国产免费一级a男人的天堂| 午夜日本视频在线| 一本色道久久久久久精品综合| 国产成人免费无遮挡视频| 免费黄网站久久成人精品| 免费av不卡在线播放| 一区二区av电影网| 男男h啪啪无遮挡| 自线自在国产av| 人体艺术视频欧美日本| 啦啦啦在线观看免费高清www| 曰老女人黄片| 亚洲精品美女久久久久99蜜臀 | 激情视频va一区二区三区| 99国产精品免费福利视频| 人人妻人人澡人人看| 丰满少妇做爰视频| 国产精品人妻久久久久久| 香蕉国产在线看| 久久 成人 亚洲| videos熟女内射| 日本黄大片高清| 免费av中文字幕在线| 激情视频va一区二区三区| 免费大片黄手机在线观看| 亚洲欧美成人综合另类久久久| 丰满少妇做爰视频| 宅男免费午夜| 视频区图区小说| av电影中文网址| 国产xxxxx性猛交| 视频在线观看一区二区三区| 午夜免费鲁丝| 午夜影院在线不卡| 蜜臀久久99精品久久宅男| xxxhd国产人妻xxx| 亚洲第一av免费看| 久久99蜜桃精品久久| 高清不卡的av网站| 国产色婷婷99| 免费高清在线观看视频在线观看| 18禁在线无遮挡免费观看视频| 国产日韩欧美在线精品| 91国产中文字幕| 久久国产亚洲av麻豆专区| 精品人妻一区二区三区麻豆| 欧美精品人与动牲交sv欧美| 一二三四中文在线观看免费高清| 国产xxxxx性猛交| 丝袜人妻中文字幕| 久久韩国三级中文字幕| 国产av一区二区精品久久| 91aial.com中文字幕在线观看| 久久韩国三级中文字幕| 男女无遮挡免费网站观看| 国产精品成人在线| 大片免费播放器 马上看| av在线老鸭窝| 国产精品无大码| 丝袜美足系列| 国产av一区二区精品久久| 国产精品一国产av| 日韩欧美精品免费久久| 日韩成人伦理影院| 少妇的丰满在线观看| 51国产日韩欧美| 捣出白浆h1v1| 欧美 日韩 精品 国产| 日本黄色日本黄色录像| 国产精品国产三级国产av玫瑰| 九色亚洲精品在线播放| www.av在线官网国产| 一本色道久久久久久精品综合| 欧美xxxx性猛交bbbb| 在线天堂中文资源库| 女人精品久久久久毛片| 欧美xxxx性猛交bbbb| 黄色配什么色好看| 国产精品人妻久久久影院| 亚洲五月色婷婷综合| 国产欧美亚洲国产| 一级片'在线观看视频| 蜜桃国产av成人99| 人成视频在线观看免费观看| 91精品国产国语对白视频| 欧美精品高潮呻吟av久久| 91午夜精品亚洲一区二区三区| 国产av国产精品国产| av线在线观看网站| 久久久亚洲精品成人影院| 精品少妇久久久久久888优播| 一个人免费看片子| 美女主播在线视频| 侵犯人妻中文字幕一二三四区| 久久鲁丝午夜福利片| 成人国产av品久久久| 日本猛色少妇xxxxx猛交久久| 天堂俺去俺来也www色官网| 亚洲成色77777| 夫妻午夜视频| 久久久久久久大尺度免费视频| 亚洲欧美精品自产自拍| 一二三四中文在线观看免费高清| 久久久精品免费免费高清| 国产欧美日韩综合在线一区二区| 侵犯人妻中文字幕一二三四区| 久久精品人人爽人人爽视色| 久久综合国产亚洲精品| 久久国产精品男人的天堂亚洲 | 久久久久久久国产电影| 亚洲精品视频女| 热re99久久国产66热| 久久婷婷青草| 最近最新中文字幕免费大全7| 亚洲欧美日韩另类电影网站| 欧美日韩综合久久久久久| av天堂久久9| 国产日韩欧美亚洲二区| 免费高清在线观看视频在线观看| 亚洲综合精品二区| 九草在线视频观看| av一本久久久久| 国产爽快片一区二区三区| 日本黄大片高清| 久久精品久久精品一区二区三区| 久久精品aⅴ一区二区三区四区 | 国产精品免费大片| 久久人人97超碰香蕉20202| av网站免费在线观看视频| 91精品伊人久久大香线蕉| 秋霞伦理黄片| 国产一区亚洲一区在线观看| 国产在线一区二区三区精| 丝袜人妻中文字幕| 亚洲国产色片| 女性生殖器流出的白浆| 国产片内射在线| 亚洲精品一区蜜桃| 街头女战士在线观看网站| 春色校园在线视频观看| 国产一区二区三区综合在线观看 | 日韩成人伦理影院| 色哟哟·www| 乱人伦中国视频| 国产精品国产三级专区第一集| 夜夜骑夜夜射夜夜干| 在线精品无人区一区二区三| 男女国产视频网站| 日韩视频在线欧美| www.色视频.com| a级毛色黄片| 香蕉丝袜av| 天美传媒精品一区二区| 久久久欧美国产精品| 亚洲伊人久久精品综合| 国产成人欧美| 欧美日韩成人在线一区二区| 熟女人妻精品中文字幕| 在线天堂最新版资源| 观看美女的网站| 久久精品熟女亚洲av麻豆精品| 日韩不卡一区二区三区视频在线| 麻豆乱淫一区二区| 亚洲精品aⅴ在线观看| av天堂久久9| 久久99热6这里只有精品| 九九在线视频观看精品| 夫妻午夜视频| 中文字幕制服av| 一级,二级,三级黄色视频| av电影中文网址| 看免费成人av毛片| 欧美性感艳星| 久久婷婷青草| 午夜激情久久久久久久| tube8黄色片| 婷婷色综合大香蕉| 亚洲av综合色区一区| 欧美精品人与动牲交sv欧美| 在线 av 中文字幕| 香蕉丝袜av| 黄色毛片三级朝国网站| 日韩,欧美,国产一区二区三区| av一本久久久久| 在线观看一区二区三区激情| 亚洲综合色惰| 亚洲av在线观看美女高潮| 国产 一区精品| 黄色视频在线播放观看不卡| 寂寞人妻少妇视频99o| 美女福利国产在线| 免费大片黄手机在线观看| 制服人妻中文乱码| 国产精品人妻久久久久久| 亚洲av福利一区| av线在线观看网站| 少妇人妻 视频| 丰满饥渴人妻一区二区三| 亚洲精品aⅴ在线观看| 美国免费a级毛片| 九草在线视频观看| 国产综合精华液| 亚洲人成网站在线观看播放| 久久人人爽人人爽人人片va| 成人综合一区亚洲| 亚洲国产欧美在线一区| 国产精品久久久av美女十八| 日日撸夜夜添| 97在线人人人人妻| 亚洲伊人久久精品综合| 精品亚洲乱码少妇综合久久| 精品国产一区二区久久| 18禁动态无遮挡网站| 18在线观看网站| 精品一区二区三区视频在线| 在线看a的网站| 日韩欧美精品免费久久| 精品国产露脸久久av麻豆| 啦啦啦视频在线资源免费观看| 精品一区二区三区视频在线| 国产探花极品一区二区| 丰满迷人的少妇在线观看| 一本一本久久a久久精品综合妖精 国产伦在线观看视频一区 | 成人国产av品久久久| 又黄又粗又硬又大视频| 最近中文字幕2019免费版| 久久韩国三级中文字幕| 成人亚洲精品一区在线观看| 黄色怎么调成土黄色| 亚洲精品一二三| 有码 亚洲区| 免费少妇av软件| 日日撸夜夜添| 免费黄频网站在线观看国产| 一本色道久久久久久精品综合| 乱人伦中国视频| 人人妻人人澡人人看| 国产精品蜜桃在线观看| 中文欧美无线码| 少妇 在线观看| 天堂中文最新版在线下载| 欧美日韩国产mv在线观看视频| 五月开心婷婷网| 最近中文字幕高清免费大全6| 免费在线观看完整版高清| 男女国产视频网站| 日韩在线高清观看一区二区三区| 精品第一国产精品| 99re6热这里在线精品视频| 国产精品一区二区在线不卡| 色5月婷婷丁香| 中文字幕最新亚洲高清| 久久亚洲国产成人精品v| 精品一区在线观看国产| 久久久精品免费免费高清| 国产激情久久老熟女| 亚洲欧美色中文字幕在线| 日韩大片免费观看网站| 亚洲第一av免费看| 国产激情久久老熟女| 男人添女人高潮全过程视频| 精品久久国产蜜桃| 亚洲人成网站在线观看播放| 成人漫画全彩无遮挡| 亚洲欧美成人精品一区二区| 国产成人av激情在线播放| 欧美激情 高清一区二区三区| 伦理电影免费视频| 欧美激情 高清一区二区三区| 搡老乐熟女国产| 亚洲精华国产精华液的使用体验| 国产精品欧美亚洲77777| 日韩欧美一区视频在线观看| av天堂久久9| 国产精品不卡视频一区二区| 青青草视频在线视频观看| 热99久久久久精品小说推荐| 国产综合精华液| 91精品三级在线观看| 天美传媒精品一区二区| 国产精品.久久久| 久久国产亚洲av麻豆专区| 国产在线免费精品| 久久精品国产综合久久久 | 狂野欧美激情性bbbbbb| 丝袜人妻中文字幕| 亚洲精品中文字幕在线视频| 精品国产乱码久久久久久小说| 国产欧美日韩一区二区三区在线| 搡女人真爽免费视频火全软件| 乱人伦中国视频| 亚洲av.av天堂| videos熟女内射| 日韩大片免费观看网站| 精品视频人人做人人爽| 国产高清国产精品国产三级| 99国产精品免费福利视频| 啦啦啦中文免费视频观看日本| 捣出白浆h1v1| 9191精品国产免费久久| 成年美女黄网站色视频大全免费| 9热在线视频观看99| 日韩一本色道免费dvd| 亚洲欧美一区二区三区国产| 亚洲一区二区三区欧美精品| 妹子高潮喷水视频| 精品国产乱码久久久久久小说| 免费黄网站久久成人精品| 亚洲国产av影院在线观看| 乱码一卡2卡4卡精品| av又黄又爽大尺度在线免费看| 巨乳人妻的诱惑在线观看| 日本猛色少妇xxxxx猛交久久| 欧美变态另类bdsm刘玥| 午夜免费鲁丝| 欧美最新免费一区二区三区| 久久国内精品自在自线图片| 欧美精品一区二区大全| 国产亚洲一区二区精品| 婷婷色综合大香蕉| 国产日韩欧美亚洲二区| 中文字幕亚洲精品专区| 亚洲欧洲国产日韩| 欧美精品高潮呻吟av久久| 久久久国产一区二区| 赤兔流量卡办理| 青春草亚洲视频在线观看| 91精品三级在线观看| 国产免费视频播放在线视频| 亚洲第一区二区三区不卡| 9热在线视频观看99| 国产黄色免费在线视频| 如何舔出高潮| 我的女老师完整版在线观看| 国产淫语在线视频| 国产日韩一区二区三区精品不卡| 欧美人与善性xxx| 国产1区2区3区精品| 久久女婷五月综合色啪小说| 欧美国产精品va在线观看不卡| 制服丝袜香蕉在线| 亚洲五月色婷婷综合| 午夜av观看不卡| 18禁观看日本| 亚洲成国产人片在线观看| 久久精品熟女亚洲av麻豆精品| 91久久精品国产一区二区三区| 少妇的逼水好多| 午夜免费鲁丝| 精品国产一区二区三区久久久樱花| 色婷婷久久久亚洲欧美| 只有这里有精品99| 韩国精品一区二区三区 | 亚洲,欧美精品.| 亚洲色图综合在线观看| 国产av码专区亚洲av| 在线 av 中文字幕| 久久99蜜桃精品久久| 国产爽快片一区二区三区| 一本色道久久久久久精品综合| 在线观看免费视频网站a站| 久久热在线av| 亚洲经典国产精华液单| 男人爽女人下面视频在线观看| 丝袜人妻中文字幕| 九九爱精品视频在线观看| 亚洲欧美中文字幕日韩二区| 日本免费在线观看一区| 久久人人爽人人片av| 日本色播在线视频| 老女人水多毛片| 2018国产大陆天天弄谢| 99re6热这里在线精品视频| 99久久综合免费| 国产福利在线免费观看视频| 99热6这里只有精品| 99香蕉大伊视频| 亚洲少妇的诱惑av| 国产精品三级大全| 不卡视频在线观看欧美| 春色校园在线视频观看| 国产极品天堂在线| 欧美丝袜亚洲另类| 七月丁香在线播放| 赤兔流量卡办理| 久久97久久精品| 香蕉丝袜av| 亚洲欧美日韩另类电影网站| 国产国语露脸激情在线看| 久久久久久久亚洲中文字幕| 国产视频首页在线观看| 久久久久久久久久久免费av| a级毛片黄视频| av播播在线观看一区| 亚洲国产av新网站| 日韩一区二区视频免费看| 欧美日韩视频高清一区二区三区二| 免费高清在线观看日韩| 最近的中文字幕免费完整| 久久精品国产自在天天线| 国产色爽女视频免费观看| 午夜激情久久久久久久| 精品人妻熟女毛片av久久网站| 9191精品国产免费久久| 男女边摸边吃奶| 狠狠精品人妻久久久久久综合| 日韩大片免费观看网站| 欧美激情 高清一区二区三区| 毛片一级片免费看久久久久| 久久久久久久久久久久大奶| 欧美成人精品欧美一级黄| 免费不卡的大黄色大毛片视频在线观看| 国产av码专区亚洲av| 丰满少妇做爰视频| 丝袜人妻中文字幕| 99久久人妻综合| 夫妻午夜视频| videosex国产| 最近的中文字幕免费完整| 免费播放大片免费观看视频在线观看| 成人综合一区亚洲| 亚洲精品久久午夜乱码| 午夜福利在线观看免费完整高清在| 黄网站色视频无遮挡免费观看| 亚洲国产精品专区欧美| 少妇被粗大的猛进出69影院 | 午夜福利,免费看| videossex国产| 午夜日本视频在线| 久久精品国产鲁丝片午夜精品| 黑人猛操日本美女一级片| 日韩熟女老妇一区二区性免费视频| 下体分泌物呈黄色| 国产国语露脸激情在线看| 午夜福利乱码中文字幕| 国产成人a∨麻豆精品| 美女大奶头黄色视频| 三上悠亚av全集在线观看| 免费久久久久久久精品成人欧美视频 | 18禁国产床啪视频网站| 亚洲婷婷狠狠爱综合网| freevideosex欧美| videos熟女内射| 9色porny在线观看| 国产有黄有色有爽视频| 韩国av在线不卡| 色婷婷av一区二区三区视频| 18+在线观看网站| 亚洲国产最新在线播放| 超色免费av| 精品亚洲成a人片在线观看| 成人18禁高潮啪啪吃奶动态图| 国产精品一区二区在线观看99| 51国产日韩欧美| 精品一区二区免费观看| 久久久欧美国产精品| 丰满少妇做爰视频| 中文欧美无线码| 国产又爽黄色视频| 啦啦啦中文免费视频观看日本| 边亲边吃奶的免费视频| 韩国精品一区二区三区 | 色94色欧美一区二区| 少妇人妻精品综合一区二区| 色婷婷av一区二区三区视频| 午夜福利,免费看| av不卡在线播放| 99香蕉大伊视频| 久热这里只有精品99| 99久久综合免费| 欧美亚洲 丝袜 人妻 在线| 少妇精品久久久久久久| 天美传媒精品一区二区| 亚洲av欧美aⅴ国产| 丝瓜视频免费看黄片| av天堂久久9| 亚洲第一区二区三区不卡| 欧美成人精品欧美一级黄| 精品久久国产蜜桃| 免费高清在线观看日韩| 欧美丝袜亚洲另类| 久久久国产一区二区| 国产免费现黄频在线看| 99热这里只有是精品在线观看| 日韩三级伦理在线观看| www.av在线官网国产| 51国产日韩欧美| 亚洲av.av天堂| 亚洲精品美女久久久久99蜜臀 | 波多野结衣一区麻豆| 久久婷婷青草| 午夜91福利影院| 黄片无遮挡物在线观看| 成年女人在线观看亚洲视频| 王馨瑶露胸无遮挡在线观看| 黄片无遮挡物在线观看| 精品一品国产午夜福利视频| 成人无遮挡网站| 亚洲美女黄色视频免费看| 久久精品熟女亚洲av麻豆精品| 丝袜美足系列| 欧美精品一区二区免费开放| 精品一区二区三卡| 亚洲成人手机| 只有这里有精品99| 亚洲三级黄色毛片| 好男人视频免费观看在线| 国产精品无大码| 精品人妻在线不人妻| 丰满少妇做爰视频| 哪个播放器可以免费观看大片| 午夜免费观看性视频| 亚洲精品国产av成人精品| 69精品国产乱码久久久| 国产av一区二区精品久久| 一本久久精品| 最黄视频免费看| 亚洲欧美成人精品一区二区| 亚洲精品日本国产第一区|