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

    A numerical flat plate friction line and its application*

    2015-02-16 06:50:39WANGZhanzhi王展智XIONGYing熊鷹SHILipan時(shí)立攀LIUZhihua劉志華

    WANG Zhan-zhi (王展智), XIONG Ying (熊鷹), SHI Li-pan (時(shí)立攀), LIU Zhi-hua (劉志華)

    Department of Naval Architecture, Naval University of Engineering, Wuhan 430033, China,

    E-mail: wzz200425@126.com

    A numerical flat plate friction line and its application*

    WANG Zhan-zhi (王展智), XIONG Ying (熊鷹), SHI Li-pan (時(shí)立攀), LIU Zhi-hua (劉志華)

    Department of Naval Architecture, Naval University of Engineering, Wuhan 430033, China,

    E-mail: wzz200425@126.com

    (Received January 6, 2014, Revised May 12, 2014)

    This paper studies the regression of a numerical two-dimensional flat plate friction line using the RANS method with the SST k-wturbulence model. Numerical simulations with different inlet turbulence kinetic energies are first conducted. Comparing with the experimental data, the finest grid and the appropriate inlet turbulence kinetic energy are selected to compute the flat plate friction resistance at 14 Reynolds numbers. Two numerical friction lines are obtained by the least squares root fitting method and one similar to that of the ITTC-1957 line and the cubic polynomials in logarithmic scales, and the results are compared to the friction line proposals available in the open literature. Finally, the full scale viscous resistance predictions of DTMB5415, KVLCC2, SUBOFF are compared between the numerical friction line and the friction line proposals available in the open literature based on the form factor approach. It is shown that the form factor keeps relatively constant via the numerical friction line for the bare hull, but the form factor concept in extrapolating the model test results is not appropriate for appended hulls. It is suggested that for computing a form factor numerically, it is best to use a numerical friction line.

    friction line, form factor, flat plate, RANS, uncertainty analysis

    Introduction

    The full scale ship resistance prediction is a fundamental issue in the ship design. At present, the model testing is an effective way to determine the resistance of a ship. The forces measured in the model scale are “extrapolated” to the full scale ship by a special procedure. The so-called form factor approach suggested by the International Towing Tank Conference (ITTC) is a popular procedure used by many modern basins. In this method, the drag of a ship is split into two independent components, the viscous resistance and the wave resistance. The viscous resistance is a function of the Reynolds number,Re , while the wave resistance is affected by the Froude number, Fr.

    The viscous resistance coefficient Cvis assumed to be proportional to the frictional resistance coefficient given by a ship-model correlation line (based on equivalent smooth flat plate flow experiments) where(1+k )is the form factor and CFis the equivalent smooth flat plate frictional resistance coefficient.

    In practice, the form factor(1+k )is assumed to be scale-independent, and the equivalent smooth flat plate resistance is computed from the ITTC-1957 friction line, then the full scale ship resistance is determined. Therefore, the validity of this approach heavily relies on the independence of the form factor on the Reynolds numberRe and the accuracy of the equivalent flat plate resistance. Unfortunately, the form factor does depend on the Reynolds number if the ITTC-1957 correlation line is chosen. Grigson[1]reanalyzed the resistance tests of the Lucy Ashton and Victory and found that the form factor is obviously scale-dependent when the ITTC-1957 friction line is adopted. Garcia[2]analyzed the experimental data of four geosim models and concluded that one sees a significant scale effect by using the ITTC-1957 frictionline, and proposed an empirical formula to correlate the form factor from the model scale to the full scale. Kouh et al.[3]investigated the scale effect of the form factor via a numerical method and found that the form factor has a near linear and increasing dependence on the Reynolds number. During the last decades, the ITTC-1957 correlation line was put into doubt and new ones were proposed[4,5]. Through the two-dimensional boundary-layer equations in zero pressure gradients using the boundary-layer velocity distribution based on Coles, the plate friction resistance could be obtained.

    With the rapid development of the CFD, it is possible to derive a plate friction line by the RANS method ranging from the model scale Reynolds number to the full scale Reynolds number. Eca and Hoekstra[6,7]presented a numerical study of the twodimensional flat plate friction resistance as a function of the Reynolds number via seven eddy-viscosity models. They found that the differences between the numerical friction lines are smaller than the difference between the four lines proposed in the open literature. However because of the low level inlet turbulence kinetic energy, the calculations show that in a region at the leading edge of the plate, one sees a transition from laminar to turbulent flows, so the predicted friction resistance is lower than that in the open literature, and the model test should be carried out in a turbulent flow condition.

    The present paper studies the regression of a numerical two-dimensional flat plate friction line using the RANS method with the SST k-wturbulence model. Numerical simulations with different inlet turbulence kinetic energies and numerical uncertainty studies are carried out. Comparing with the experimental data, the finest grid and the appropriate inlet turbulence kinetic energy are selected to compute the flat plate friction resistance at 14 Reynolds numbers. Two numerical friction lines are obtained using the least squares root fitting method: from one similar to that of the ITTC-1957 friction line and cubic polynomials in logarithmic scales. The results are compared to the friction line proposals available in the open literature. The full scale viscous resistance predictions of DTMB5415, KVLCC2, SUBOFF are compared between the numerical friction line and the friction line proposals available in the open literature based on the form factor approach, and the direct CFD result. It is shown that the form factor keeps relatively constant with respect to the numerical friction line for a bare hull.

    1. Numerical model

    1.1 Governing Equations field over a two-dimensional flat plate are the instantaneous conservation of mass (continuity equation) and momentum (Reynolds averaged Navier-Stokes equation, RANS). These equations can be expressed as follows:

    where, all variables are time-averaged,t,ui,ρ,p,are the time, the velocity, the fluid density, the static pressure, the dynamic viscosity, the body force per unit volume and the Reynolds stress, respectively.

    1.2 Turbulence model

    An additional equation is needed in order to obtain the unknown Reynolds stress. The turbulence model is the closure equation which combines the fluctuating quantities and the time averaged ones. Here, the SST k-w turbulence model is selected. The SST k-wturbulence model was developed by Menter[8]to effectively blend the robust and accurate formulation of the standardk-wmodel in the near-wall region with the free-stream independence of the standardkεmodel in the far field. The SST k-wturbulence model refines the standardk-wmodel in the following manner: firstly, in the SSTk-wturbulence model, a damped cross-diffusion derivative term is incorporated in thewequation, secondly, the definition of the turbulent viscosity is modified to account for the transport of the turbulent shear stress. These features make the SST k-wturbulence model more accurate and reliable for a wider class of flows than the standardk-wmodel.

    The turbulence kinetic energykand the specific dissipation rateware obtained from the following transport equations:

    The governing equations for the turbulent flow

    and

    where

    1.3 Computational domain and boundary conditions

    The plate length is Lp, the undisturbed inlet flow velocity isU∞, and the Reynolds number is defined as Re=U∞Lp/ν, andνis the kinematic viscosity.

    The computational domain is a rectangle of 4Lpin length with the inlet boundary located 0.5Lpupstream of the leading edge of the plate, the outlet boundary2.5Lpdownstream of the trailing edge and the top boundary0.5Lpaway from the plate. The structured grids are generated to discretize the flow domain. The grid node distribution in the longitudinal(x)direction follows an uneven distribution along the plate where the leading edge of the plate is refined and a one-sided geometric law is applied upstream and downstream of the plate. In the normal(y)direction, a one-sided geometric law is applied. On the plate, the number of grid nodes in the longitudinal direction, Nx=2800, and that in the normal direction,Ny= 200, 89% of that in the longitudinal direction.

    The medium is water with the density of 998.2 kg·m-3and the fluid viscosity of 0.001003 kg·m-1·s-1. The fluid control and turbulence equations are solved numerically using the finite volume method, which is realized by the Ansys Fluent code. The inlet boundary is set as the velocity inlet, the velocity varies according to the Reynolds number, the outlet boundary is set as the pressure outlet, the top field boundary is set as the symmetry plane. The plate is a wall assumed to be hydrodynamically smooth. The governing equations and the turbulence model are discretized by the finite volume method with a second order upwind scheme and the pressure-velocity coupling is realized by the SIMPLEC method.

    2. Results and discussions

    2.1 Influence of inlet turbulence kinetic energy

    In practical turbulent flows, the transition from laminar to turbulent flows may be promoted by the turbulence kinetic energy of the undisturbed flow,k∞. Experimental data are available for the flat plate flow field with different levels of inlet turbulence kinetic energies according to the ERCOFTAC Classic Database.

    In the transition region from laminar to turbulent flows, there is a bouncing of the local skin friction coefficient,Cf, along the plate. The local skin friction coefficient,Cf, is defined as

    Figure 1 presents Cfalong the flat plate for the calculations atlg(Re)=6.25. The plot includes the experimental results, the numerical predictions with different inlet turbulence kinetic energies (k∞=and the Blasius formula for laminar flow,

    Fig.1 Local friction coefficient,Cf, along the flat plate at lg(Re)=6.25

    From Fig.1, it is found that:

    (1) both the experimental and CFD results agree well with the Blasisus formula in the laminar region,

    (2) the numerical simulations for all different inlet turbulence kinetic energies produce a smallerRexthan where such departure appears in the experiment with k∞=at which transition occurs,Rex= U∞x/ν, andxis the distance to the leading edge of the plate,

    (3) when the inlet turbulence kinetic energy increases, the transition occurs for lowerRex,

    (4) with the increase of the inlet turbulence kinetic energy, the amplitudes and the range of the local skin friction in the transition region are reduced,

    (5) the value ofCfis nearly independent of the inlet turbulence kinetic energy for high Rex.

    It must be noticed that the model test is carried out in a relatively fully developed turbulent condition via the the turbulence stimulator, and it is appropriate to set a higher inlet turbulence kinetic energy for an earlier transition at lower Reynolds numbers to obtain the flat plate friction resistance.

    Fig.2 Local friction coefficient,Cf, along the flat plate at lg(Re)=9

    2.2 Influence of y+at high Reynolds number

    Under the high Reynolds number flow condition, the boundary layer becomes much thinner compared to the low Reynolds number flow, thus it is not appropriate to place certain cells in the viscous sublayer because too large cell aspect ratios to the wall may lead to numerical problems, and the wall function approach is a good choice to resolve the problem. The wall functions are a set of semi-empirical formulae and functions that in effect “bridge” or “l(fā)ink” the solution variables at the near-wall cells and the corresponding quantities on the wall. The nondimensional spacing, y+, is defined by

    where?yis the distance of the first grid node to the wall,uτis the friction velocity, and uτ=τwis the shear stress at the wall.

    The wall function works well when the first grid node is located in the logarithmic boundary layer, 11.225≤y+≤300.

    In order to study the influence of the nondimensional spacing,y+, on the flat plate friction resistance, 4 sets of grids are designed to approach the target of y+=50,y+=100,y+=200,y+=300. Figure 2 is the comparison of the local friction coefficient along the plate between theexperimentaldata[9]and the results of Karman-Schoenherr equation,

    where

    From Fig.2, it is found that:

    (1) the numerical results agree well with the results of the Karman-Schoenherr equation and the experimental data, indicating that the wall function approach is capable of simulating high Reynolds number flows equivalent to the full-scale ship Reynolds number,

    (2) the local friction coefficient along the flat plate is nearly independent of y+when the first grid node is located in the logarithmic boundary layer, in which the wall function works well, indicating that the wall function approach is not sensitive to the grid distribution in calculating the friction resistance.

    2.3 Numerical uncertainty analysis

    In the numerical simulation of the two-dimensional flat plate flow field using the RANS method, the uncertainty analysis must be carried out for the solution and the computational grid. The RANS method is a steady simulation method, thus the simulation uncertaintyUSNis composed of the grid uncertainty UGand the iterative uncertainty

    The assessment of UIis made by observing the oscillatory amplitude of the calculated value with time Zhang et al.[10], the iterative uncertainty is smaller than the grid uncertainty by 2 orders of magnitude, than the grid uncertainty by 2 orders of magnitude, and can be neglected. Thus, the numerical simulation uncertainty,USN, is equal to the grid uncertainty,UG.

    The uncertainty analysis is carried out for the flat plate friction resistance coefficient. Many studies made CFD uncertainty analyses[11-14]. Based on the recommended Procedure 7.5-03-01-01 from the International Towing Tank Conference, three sets of computational grids with different grid scales (i.e., fine, medium, and coarse grids) could be built for the numerical simulation of the flat plate flow field. The refinement ratio,rG, is2in each direction of the coordinate. The cell numbers of the grids are 560 000, 281 160 and 140 000.

    Table1 The flat plate friction coefficient for three grids

    Table1 shows the calculated flat plate friction coefficients for the three grids from lg(Re) =6.25to lg(Re)=9.5, and the uncertainty analysis procedure is illustrated below.

    The values calculated by the coarse, medium, fine grids are denoted as SGC,SGMand SGF, respectively. Grid changes(ε)in the coarse-to-medium and medium-to-fine grids at every Reynolds number are defined by

    The grid convergence ratio,PG, is defined by

    The order of accuracy,PG, is defined by

    The correction factor,CG, is defined by

    The grid uncertainty at every Reynolds number, UG, is defined by

    The parameters used for the uncertainty verification of the flat plate friction coefficient at every Reynolds number are shown in Table 2. The grid convergence ratio,RG, for every Reynolds number is less than 1, indicating that the computational grids ensure a monotonic convergence.

    Table2 The flat plate friction coefficient for three grids

    2.4 Curve fit of the numerical data

    In order to obtain a numerical two-dimensional flat plate line, a curve fit must be performed to the numerical data for the 14 Reynolds numbers. Here, the cubic polynomial in the logarithmic scales and a formula like that of the ITTC-1957 are tested.

    The parameters,a0,a1,a2,a3,b1,b2,b3, in the two equations are estimated using a weighed least squares method, and they are listed in Table 3.

    Table3 The parameters in Eqs.(14) and (15)

    There are many analytical equations proposed for the friction resistance coefficient of a flat plate as a function of the Reynolds number. In the present study, 5 alternative proposals from the open literature are selected for comparison of the numerical two-dimensional flat plate friction lines.

    Fig.3 Comparison of the numerical friction lines and proposals from the open literature

    The first is the ITTC-1957 correction line, the second is the line derived by Katsui et al.[5], the third is the line proposed by Grigson[4], the fourth is the Schoenherr line, which is very popular in the ATTC, and the last is recently proposed by Eca et al.[8], which is also a numerical friction line obtained via the RANS method coupled with the SST k-wturbulence model.

    Figure 3 shows the comparison of the numerical two-dimensional flat plate friction lines and proposals from the open literature, part (b) is the zoom-in in the case when lg(Re)≤7.75.

    From Fig.3, we found that:

    (1) at high Reynolds numbers, there is an excellent agreement between the numerical line and the proposals of ITTC-1957, Katsui, Schoenherr and Eca, while the line derived by Grigson is slightly higher than the others.

    (2) at low Reynolds numbers, the differences between these lines become much larger. Obviously, the ITTC-1957 line has the largest slope and makes the values of CFlarger than those of the others because it is not a merely turbulent flat plate friction line, and it is modified for the purpose of the resistance extrapolation from the model scale to the full scale. Eca’s line has the least slope and gives the lowest value of CFbecause there is a laminar flow region at the edge of the plate in the numerical simulation. The numerical result is very close to the Grigson line when lg(Re)<7.

    (3) the line derived by Grigson has a discontinuity of slope atlg(Re)=7, while the others are continuous.

    3. Application of the numerical friction line

    In the towing-tank test, the form factor is determined by the Prohaska’s method, measuring the total resistance at low speed, where the wave resistance contribution vanishes or can be estimated. This method might be quite inaccurate. Compared with the experimental procedure, the RANS method is more straightforward for determining the form factor: the deformation of the water surface is simply not taken into account, and a “double-body flow” computation is done. The scale effect on the form factor using different friction lines can be simply determined if making such computations in both the model scale and the full scale.

    3.1 DTMB5415

    The DTMB 5415 is first tested for determining whether the numerical friction line is suitable for the model-ship extrapolation method. It was conceived as a preliminary design for a navy surface combatant. The hull geometry includes both a sonar dome and a transom stern. The propulsion is provided by twinopen-water propellers driven by shafts supported by struts. Since 1996, the DTMB 5415 model has been adopted by the ITTC as a recommended benchmark for the CFD validation for resistance and propulsion[15,16]. The hull’s geometry is shown in Fig.4.

    Fig.4 Geometry of DTMB 5415

    The model speed is 2.097 m/s, and the Reynolds number based on the model length and the model speed is 1.19× 107and the Froude number is Fr= 0.28. The full scale ship speed is 10.288 m/s, and the Reynolds number is 1.45× 109. This is one of test cases in the Gothenburg 2010 Workshop on Numerical Ship Hydrodynamics.

    In the “double-body flow” computation at the model scale, the near-wall spacing is set to ensure the nondimensional spacing being at a target of y+=100 at the middle of the hull. The computational domain has 300×60×90nodes (1.62× 106cells) in the longitudinal, wall-normal and girthwise directions. In the full scale, the near-wall spacing is set to ensure the nondimensional spacing being at a target ofy+=300 at the middle of the hull. The computational domain has 750×200× 140nodes (2.1× 107cells) in the longitudinal, wall-normal and girthwise directions. All computational grids are constructed of the HO topology, refined longitudinally towards the bow and the stern to resolve the large velocity gradients. The computational grids are shown in Fig.5.

    Fig.5 Computational grids of DTMB 5415

    The inlet boundary is set as the velocity inlet, the outlet boundary is set as the pressure outlet, and the top field boundary is set as a symmetry plane. The ship‘s wall is assumed to be hydrodynamically smooth. The governing equations and the turbulence model are discretized by the finite volume method with a second order upwind scheme and the pressure-velocity coupling is achieved by the SIMPLEC method. The full scale computation is carried out in parallel processing in 32 cores (Intel Xeon E5-2670, 2.6 GHz), realized in the Dawning TC4600 high performance computer.

    The computed resistance components in the model scale and full scale ships are summarized in Table 4, where,CVis the viscous resistance coefficient,CFis the frictional resistance coefficient and CVPis the viscous pressure resistance coefficient.

    In the form factor approach, the viscous resistance coefficient of a full scale ship is assumed to be proportional to the friction coefficient of a flat plat. In the so-called ITTC-1978 method, the ITTC-1957 model-ship correlation line is chosen for the plate friction line. However, it must be noticed that the ITTC-1957 model-ship correlation line is not a real friction line, it is a modified one for better prediction of the ship friction resistance if no form factor method is used.

    Evidently, the form factors depend on the plate friction coefficients assumed. In the last section, Fig.4 shows the numerical two-dimensional flat plate friction lines and other proposals from the open literature. Although the difference in the Re-dependence is limited but it is sufficient for a very different performance predictions. The form factors for model scale and full scale ships based on these lines are listed in Table 5, where, the ratio is the full scale ship form factor to the model scale ship form factor.

    From Table 5, the following observations can be made.

    (1) The form form factor(1+k )increases substantially from the model scale to the full scale ships if the ITTC-1957 line is chosen. For the DTMB 5415, an extrapolation using a fixed form factor would underestimate the full scale viscous resistance by 4.86%.

    (2) The form factor tends to be constant if the numerical line is chosen.

    Table5 Form factors of DTMB5415 for model scale and full scale ships based on different friction lines

    (3)The proportionality between the ship viscous resistance and the flat plate resistance is well reproduced if one of the modern friction lines is used (Katsui line and Grigson line).

    3.2 KVLCC 2 tanker

    The KVLCC2 is the second test case and is the second variant of the MOERI tanker with more U-shaped stern frame-lines. The local flow characteristics around the VLCC hull were extensively studied experimentally and the results were documented in detail and used as the validation resource for numerical solutions[17]. The hull’s geometry is shown in Fig.6.

    Fig.6 Geometry of KVLCC2

    The model speed is 1.047 m/s, the Reynolds number based on the model length and the model speed is 4.6× 106and the Froude number is Fr= 0.142. The full scale KVLCC2 speed is 7.973744 m/s, and the Reynolds number is 2.54× 109. This is also one of test cases in the Gothenburg 2010 Workshop on Numerical Ship Hydrodynamics.

    Fig.7 Computational grid of KVLCC2

    In the model scale, the computational domain is meshed with 400×60× 150nodes (3.6× 106cells) in the longitudinal, wall-normal and girthwise directions. While in the full scale, it is meshed with800×200× 300 nodes (4.8× 107cells) in the longitudinal, wallnormal and girthwise directions. All computational grids are constructed in the HO topology, refined longitudinally towards the bow and the stern to resolve the large velocity gradients. The computational grid is shown in Fig.7.

    The boundary condition is the same as the DTMB5415.The full scale computation is carried out in parallel processing in 64 cores (Intel Xeon E5-2670, 2.6 GHz) in the Dawning TC4600 high performance computer.

    Table6 Computed resistance components of KVLCC2 for model scale and full scale ships

    The computed resistance components of the KVLCC2 for the model scale and full scale ships are summarized in Table 6, which shows that the viscous resistance coefficient decreases by a factor of 0.4276 from the model scale to full scale ships.

    The form factors of the KVLCC2 for the model scale and full scale ships based on different friction lines are listed in Table 7.

    Table7 shows that the Schoenherr line gives the best prediction of the full scale ship viscous resistance if a constant form factor is used, the numerical line and the Katsui line give a fairly constant form factor(-1.2%, -2.9%), although not as precise as that of the DTMB5415, while the ITTC-1957 line and the Grigson line yield a significant increase (3.5%, 4.5%) of the from factor1+k from the model scale to the full scale. This phenomenon is encouraging.

    3.3 Darpa suboff

    The Suboff model was built by the David Taylor Research Center as a recommended benchmark for the CFD validation for the flow field and the resistance[18]. It has an overall length of 4.356 m. The hull is composed of a fore-body of 1.016 m in length, a parallel mid-body of 2.229 m in length, an aft-body of 1.111 m in length, and an end cap of 0.095 m in length.The scale ratio is 24.

    Table7 Form factors of KVLCC2 in model scale and full scale based on different friction lines

    The model speed is 3.0452 m/s, and the Reynolds number based on the model length and the model speed is 1.32× 107. The full scale Suboff speed is 14.9186 m/s, and the Reynolds number is1.55× 109.

    For the bare hull in the model scale, the computational domain is meshed with350×60×50nodes (1.05× 106cells) in the longitudinal, wall-normal and girthwise directions. While for the bare hull in the full scale, it is meshed with700×200× 100nodes (1.4× 107cells) in the longitudinal, wall-normal and girthwise directions. For the hull with full appendages in the model scale, the computational domain is meshed with480×70× 100nodes (3.6× 106cells) in the longitudinal, wall-normal and girthwise directions. While for the hull with full appendages in the full scale, it is meshed with850×200×200nodes (3.2× 107cells) in the longitudinal, wall-normal and girthwise directions. All the computational grids are constructed of the HO topology, refined longitudinally towards the bow, the stern and the appendages to resolve the large velocity gradients. The computational grids are shown in Figs.8 and 9.

    Fig.8 Computational grids of Suboff bare hull

    Fig.9 Computational grids of Suboff with full appendages

    The boundary condition is the same as the DTMB5415. The full scale computation is carried out in parallel processing in 32 cores (Intel Xeon E5-2670, 2.6 GHz) in the Dawning TC4600 high performance computer.

    The computed resistance components of the Suboff bare hull and those with full appendages for model scale and full scale ships are summarized in Table 8.

    Table8 shows the numerical results of viscous resistance coefficients, which agree quite well with the experimental ones for both the bare hull and those with full appendages. For the bare hull, the viscous resistance coefficient decreases by a factor of 0.5276 from the model scale to the full scale ships; while the factor is 0.5591 for the full appendage Suboff.

    The form factors of the Suboff bare hull and those with full appendages for model scale and full scale ships based on different friction lines are listed in Table 9.

    Table9 shows that for the bare hull, the numerical line gives almost constant form factor (-0.54%),the Katsui and Grigson lines give a fairly constant form factor (-1.12% and -1.17%, respectively). Again, the ITTC-1957 line will yield a significant increase (4.04%) of the form factor(1+k )from the model scale to the full scale. While for the Suboff with full appendages, the form factor increases substantially from the model scale to the full scale ships for all friction lines, the ITTC-1957 line suffers a great scale effect. The assumption of a constant form factor seems not reasonable for the hull with appendages, which may be due to the different local Reynolds number of the appendages in the full scale. Anyway, the extrapolation method from the model scale to the full scale ships with appendages needs to be studied furthers.

    4. Conclusions

    The present work is devoted to a numerical plate friction line via the RANS method with the SST k-w turbulence model. The inlet turbulence kinetic energy and numerical uncertainty analyses are conducted. Comparing with the experimental data, the finest grid and the appropriate inlet turbulence kinet ic energy a re selected to compute the flat plate frictionresistanceat14 Reynolds numbers. Two numerical friction lines are obtained using the least squares root fitting method, from one similar to that of the ITTC-1957 line and cubic polynomials in logarithmic scales, and are compared to the friction line proposals available in the open literature.

    Table8 Computed resistance components of Suboff for model scale and full scale ships

    Table9 Form factors of Suboff for model scale and full scale ships based on different friction lines

    The form factor of DTMB5415, KVLCC2, SUBOFF are compared between the model scale and the full scale ships using different friction lines. It is shown that the numerical line sees little scale effect from the model scale to the full scale ships, while the assumption of a constant form factor using the ITTC-1957 line is not reasonable. It is suggested that when computing a form factor numerically, it is best to use a numerical friction line. Several more ships need to be studied to validate the numerical friction line for a fundamental support to the form factor concept.

    [1] GRIGSON C. W. A reanalysis of Lucy Ashton and victory experiments[J]. Transactions of RINA, 1996, 138(1): 117-130.

    [2] GARCIA G. A. On the form factor scale effect[J]. Ocean Engineering, 2000, 26(1): 97-109.

    [3] KOUH J. S., CHEN Y. J. and CHAU S. W. Numerical study on scale effect of form factor[J]. Ocean Engineering, 2009, 36(1): 403-413.

    [4] GRIGSON C. W. A planar friction algorithm and its use in analyzing hull resistance[J]. Transactions of RINA, 1999, 141: 76-115.

    [5] KATSUI T., ASAI H. and HIMENO Y. The proposal of a new friction line[C]. Fifth Osaka Colloquium on Advanced CFD Application to Ship Flow and Hull Form Design. Osaka, Japan, 2005.

    [6] ECA L., HOEKSTRA M. On the accuracy of the numerical prediction of scale effects on ship viscous resistance-computational method in marine engineering[C]. International Conference on Computational Methods in Marine Engineering. Oslo, Norway, 2005.

    [7] ECA L., HOEKSTRA M. The numerical friction line[J]. Journal of Marine Science and Technology, 2008, 13(4): 328-345.

    [8] MENTER F. R. Two-equation eddy-viscosity turbulence models for engineering applications[J]. AIAA Journal, 1994, 32(8): 1598-1605.

    [9] WATSON R. D., HALL R. M. and ANDERS J. B. Review of skin friction measurements including recent high-Reynolds number results from NASA Langley NTF[C]. Proceedings of AIAA Paper 2000-2392. Denver, USA, 2000.

    [10] ZHANG Nan, SHEN Hong-cui and YAO Hui-zhi. Uncertainty analysis in CFD for resistance and flow field[J]. Journal of Ship Mechanics, 2008, 12(2): 211-224(in Chinese).

    [11] STERN F., COLEMAN H. W. and PATERSON E. G. Comprehensive approach to verification and validation of CFD simulations-Part 1: Methodology and procedures[J]. Journal of Fluids Engineering, 2001, 123(4): 793-802.

    [12] STERN F., WILSON R. and SHAO J. Quantitative V&V of CFD simulations and certification of CFD codes[J]. International Journal for Numerical Methods in Fluid, 2006, 50(11): 1335-1355.

    [13] ECA L., HOEKSTRA M. and HAY A. Verification of RANS solves with manufactured solutions[J]. Engineering with Computers, 2007, 23(4): 253-270.

    [14] YAO Zhen-qiu, SHEN Hong-cui and GAO Hui. A new methodology for the CFD uncertainty analysis[J]. Journal of Hydrodynamics, 2013, 25(1): 131-147.

    [15] OLIVIERI A., PISTANI F. and AVANAINI A. et al. Towing tank experiments of resistance, sinkage and trim, boundary layer, wake, and free surface flow around a naval combatant INSEAN 2340 model[R]. Iowa, USA: The University of Iowa, 2001.

    [16] LARSSON L., STERN F. and BERTRAM V. Benchmarking of computational fluid dynamics for ship flows: The Gothenburg 2000 workshop[J]. Journal of Ship Research, 2003, 47(1): 63-81.

    [17] KIM W. J., VAN S. H. and KIM D. H. Measurement of flows around modern commercial ship models[J]. Experiment in Fluids, 2001, 31(5): 567-578.

    [18] BULL P. The validation of CFD predictions of nominal wake for the SUBOFF fully appended submarine geometry[C]. Proceedings of 21st Symposium on Naval Hydrodynamics. Trondheim, Norway, 1996.

    * Project supported by the National Natural Science Foundation of China (Grant No. 51179198, 51479207), the High Technology Marine Scientific Research Project of Ministry of Industry and Information Technology of China (Grant No. [2012]534).

    Biography: WANG Zhan-zhi (1986- ), Male, Ph. D., Lecturer

    Correspondindg author: XIONG Ying,

    E-mail: ying_xiong28@126.com

    成人一区二区视频在线观看| 成人毛片a级毛片在线播放| 国产一区二区三区视频了| 亚洲 国产 在线| 91av网一区二区| 精品久久国产蜜桃| 禁无遮挡网站| 亚洲男人的天堂狠狠| 成人精品一区二区免费| 久久久久国产精品人妻aⅴ院| 欧美激情国产日韩精品一区| 成人永久免费在线观看视频| 午夜日韩欧美国产| 色5月婷婷丁香| 91狼人影院| 搡老熟女国产l中国老女人| 国产探花在线观看一区二区| 免费观看精品视频网站| 欧美高清成人免费视频www| 18+在线观看网站| 91狼人影院| 亚洲精品一卡2卡三卡4卡5卡| 久久精品夜夜夜夜夜久久蜜豆| 香蕉av资源在线| 精品午夜福利在线看| 欧美一区二区国产精品久久精品| 黄色视频,在线免费观看| 69av精品久久久久久| 久久久久性生活片| 麻豆av噜噜一区二区三区| 少妇人妻精品综合一区二区 | 精品久久久久久久末码| 欧美一区二区亚洲| 男女之事视频高清在线观看| 嫩草影院入口| 最后的刺客免费高清国语| 国产成人影院久久av| 国产成人av教育| 亚洲欧美日韩高清在线视频| 国产伦一二天堂av在线观看| 久久精品91蜜桃| 熟女电影av网| 午夜福利在线在线| 大型黄色视频在线免费观看| 免费看a级黄色片| 日本免费a在线| 99热这里只有是精品在线观看| 乱人视频在线观看| 亚洲七黄色美女视频| 欧美黑人巨大hd| 一个人观看的视频www高清免费观看| 亚洲人成网站高清观看| 三级毛片av免费| 桃红色精品国产亚洲av| 一个人观看的视频www高清免费观看| 亚洲精华国产精华精| 国模一区二区三区四区视频| 最近最新免费中文字幕在线| 毛片一级片免费看久久久久 | 久久亚洲真实| 午夜福利在线观看吧| 春色校园在线视频观看| 春色校园在线视频观看| 狂野欧美白嫩少妇大欣赏| 热99在线观看视频| 午夜福利在线在线| 亚洲国产精品sss在线观看| 一级a爱片免费观看的视频| 精品乱码久久久久久99久播| 女生性感内裤真人,穿戴方法视频| 搡女人真爽免费视频火全软件 | 女人被狂操c到高潮| 深夜精品福利| 日日撸夜夜添| 亚洲精品日韩av片在线观看| 日韩中字成人| 日韩人妻高清精品专区| 日韩精品青青久久久久久| 自拍偷自拍亚洲精品老妇| 久久这里只有精品中国| 女生性感内裤真人,穿戴方法视频| 亚洲av第一区精品v没综合| 亚洲欧美日韩高清专用| 性欧美人与动物交配| 中文字幕人妻熟人妻熟丝袜美| 亚洲午夜理论影院| 99久久九九国产精品国产免费| 又粗又爽又猛毛片免费看| 国产一区二区亚洲精品在线观看| 97超级碰碰碰精品色视频在线观看| 国产精品人妻久久久久久| 人人妻人人澡欧美一区二区| 69av精品久久久久久| 久久精品国产99精品国产亚洲性色| 成人美女网站在线观看视频| 午夜日韩欧美国产| 中文字幕av在线有码专区| 欧美+日韩+精品| 国产大屁股一区二区在线视频| 亚洲欧美日韩高清在线视频| 精品久久久噜噜| 国产在线精品亚洲第一网站| 熟妇人妻久久中文字幕3abv| 在线天堂最新版资源| 国产一区二区三区av在线 | 欧美精品国产亚洲| 熟女人妻精品中文字幕| 国产毛片a区久久久久| 久久亚洲精品不卡| 日韩亚洲欧美综合| 亚洲欧美日韩高清在线视频| 国产乱人视频| 久久久久久大精品| 特大巨黑吊av在线直播| 欧美三级亚洲精品| 免费高清视频大片| 女人被狂操c到高潮| 色综合色国产| 国产精品免费一区二区三区在线| 三级毛片av免费| 欧美激情国产日韩精品一区| 亚洲av免费在线观看| 最近在线观看免费完整版| 国产激情偷乱视频一区二区| 黄片wwwwww| 午夜福利在线观看吧| 精品久久久噜噜| 日韩一本色道免费dvd| 亚洲最大成人av| 人妻久久中文字幕网| 变态另类成人亚洲欧美熟女| 国产精品爽爽va在线观看网站| 亚洲成a人片在线一区二区| a级一级毛片免费在线观看| 亚洲电影在线观看av| 国产免费一级a男人的天堂| 一级黄色大片毛片| 国产精品久久久久久久电影| 国产男靠女视频免费网站| 色综合站精品国产| 深爱激情五月婷婷| 久久国内精品自在自线图片| 一级av片app| 免费一级毛片在线播放高清视频| 亚洲电影在线观看av| 午夜a级毛片| 国产一区二区亚洲精品在线观看| 一区二区三区四区激情视频 | 中文字幕av在线有码专区| 真实男女啪啪啪动态图| 久久久久免费精品人妻一区二区| 日韩中文字幕欧美一区二区| 日本色播在线视频| 国产一区二区三区在线臀色熟女| 欧美高清性xxxxhd video| 成人特级黄色片久久久久久久| 国产精品日韩av在线免费观看| 夜夜夜夜夜久久久久| 少妇高潮的动态图| 国产精品人妻久久久影院| aaaaa片日本免费| 国内精品一区二区在线观看| 长腿黑丝高跟| 欧美极品一区二区三区四区| ponron亚洲| 亚洲av一区综合| 亚洲久久久久久中文字幕| 国产精品一区二区性色av| 国产精品爽爽va在线观看网站| 午夜福利成人在线免费观看| 国产女主播在线喷水免费视频网站 | av天堂中文字幕网| 极品教师在线免费播放| 日本三级黄在线观看| 我的老师免费观看完整版| 少妇裸体淫交视频免费看高清| 别揉我奶头 嗯啊视频| 亚洲中文字幕一区二区三区有码在线看| or卡值多少钱| 99热这里只有精品一区| 国产高潮美女av| 一夜夜www| 亚洲天堂国产精品一区在线| 亚洲成人精品中文字幕电影| 午夜精品在线福利| 国产精品99久久久久久久久| 久久亚洲精品不卡| 最近最新免费中文字幕在线| 欧美+日韩+精品| 国产私拍福利视频在线观看| 欧美3d第一页| 久久香蕉精品热| 久久久久久大精品| 色尼玛亚洲综合影院| 永久网站在线| 女人十人毛片免费观看3o分钟| 最近最新中文字幕大全电影3| 欧美+亚洲+日韩+国产| 一个人观看的视频www高清免费观看| 亚洲欧美日韩高清在线视频| 黄片wwwwww| av女优亚洲男人天堂| 成人高潮视频无遮挡免费网站| 99热只有精品国产| 校园春色视频在线观看| 欧美一区二区精品小视频在线| 欧美性猛交╳xxx乱大交人| 日本在线视频免费播放| 日本与韩国留学比较| 最近最新免费中文字幕在线| 内射极品少妇av片p| 欧美黑人欧美精品刺激| 中文字幕免费在线视频6| 免费在线观看影片大全网站| 嫩草影院入口| 伦精品一区二区三区| 88av欧美| 在线播放国产精品三级| 亚洲成人中文字幕在线播放| 久久九九热精品免费| 国产伦在线观看视频一区| 美女高潮喷水抽搐中文字幕| 欧美另类亚洲清纯唯美| 日本五十路高清| 男插女下体视频免费在线播放| 亚洲,欧美,日韩| 免费不卡的大黄色大毛片视频在线观看 | 国产精品久久久久久久久免| 国产乱人视频| 国产精品嫩草影院av在线观看 | 最近视频中文字幕2019在线8| 中文在线观看免费www的网站| av在线观看视频网站免费| 国产成人影院久久av| 亚洲成av人片在线播放无| 美女黄网站色视频| 色尼玛亚洲综合影院| 又爽又黄无遮挡网站| 久久精品91蜜桃| 亚洲欧美日韩卡通动漫| 欧美+日韩+精品| 天堂√8在线中文| 成熟少妇高潮喷水视频| 国产午夜精品久久久久久一区二区三区 | 在线观看美女被高潮喷水网站| 真实男女啪啪啪动态图| 久久99热这里只有精品18| 国产中年淑女户外野战色| 在线观看66精品国产| 在线观看美女被高潮喷水网站| 日韩大尺度精品在线看网址| 久久久成人免费电影| 免费在线观看成人毛片| 22中文网久久字幕| 给我免费播放毛片高清在线观看| 亚洲经典国产精华液单| 欧美激情在线99| 午夜日韩欧美国产| 国产精品三级大全| 久久国产精品人妻蜜桃| 成人国产麻豆网| 非洲黑人性xxxx精品又粗又长| 欧美日韩综合久久久久久 | 人人妻人人看人人澡| 国产精品国产三级国产av玫瑰| 成人鲁丝片一二三区免费| 国产一级毛片七仙女欲春2| 日本a在线网址| 午夜视频国产福利| 香蕉av资源在线| 99久国产av精品| 国内少妇人妻偷人精品xxx网站| 在线观看午夜福利视频| 两性午夜刺激爽爽歪歪视频在线观看| 成年女人永久免费观看视频| 精品国内亚洲2022精品成人| 亚洲国产欧洲综合997久久,| 日本一本二区三区精品| 久久99热这里只有精品18| 欧美中文日本在线观看视频| 22中文网久久字幕| 成人av在线播放网站| 免费av不卡在线播放| 日韩欧美免费精品| 国产精品人妻久久久久久| 亚洲国产欧洲综合997久久,| 国产在线精品亚洲第一网站| 他把我摸到了高潮在线观看| 亚洲男人的天堂狠狠| 色在线成人网| 成人性生交大片免费视频hd| 看黄色毛片网站| 亚洲专区中文字幕在线| 亚洲无线在线观看| 91久久精品国产一区二区成人| 久久久久精品国产欧美久久久| 国产高清视频在线播放一区| 最近视频中文字幕2019在线8| 91麻豆av在线| 一级黄色大片毛片| 变态另类丝袜制服| 亚洲欧美激情综合另类| 伦理电影大哥的女人| 美女黄网站色视频| 中文字幕精品亚洲无线码一区| 特级一级黄色大片| 高清毛片免费观看视频网站| 国产色爽女视频免费观看| 久久亚洲精品不卡| 日韩欧美在线二视频| 久久人人精品亚洲av| 欧美xxxx性猛交bbbb| 我的老师免费观看完整版| 高清毛片免费观看视频网站| 国产亚洲精品久久久久久毛片| 免费搜索国产男女视频| 黄色配什么色好看| 久久久久久久精品吃奶| 看免费成人av毛片| 日韩av在线大香蕉| 看十八女毛片水多多多| 最近最新免费中文字幕在线| 日韩精品中文字幕看吧| 狠狠狠狠99中文字幕| 亚洲成av人片在线播放无| 国产一区二区激情短视频| 国产色婷婷99| 久久中文看片网| 欧美高清性xxxxhd video| 日韩中字成人| 精品久久久久久成人av| 国产精品av视频在线免费观看| 特大巨黑吊av在线直播| 一a级毛片在线观看| 91麻豆av在线| 日韩精品有码人妻一区| 亚洲av电影不卡..在线观看| 最新中文字幕久久久久| 别揉我奶头 嗯啊视频| avwww免费| 国产精品久久久久久久久免| 欧美bdsm另类| 国产精品永久免费网站| 国产私拍福利视频在线观看| 国产在线精品亚洲第一网站| 国产精品一及| 亚洲在线自拍视频| 床上黄色一级片| 淫妇啪啪啪对白视频| 91麻豆精品激情在线观看国产| 别揉我奶头 嗯啊视频| 国内揄拍国产精品人妻在线| 精品久久久久久成人av| 丰满的人妻完整版| 最近最新中文字幕大全电影3| 久久久久久久久大av| 国产精品久久电影中文字幕| 婷婷六月久久综合丁香| 国产aⅴ精品一区二区三区波| 久久这里只有精品中国| 亚洲不卡免费看| 日韩中字成人| 色尼玛亚洲综合影院| 久久婷婷人人爽人人干人人爱| 欧美日韩亚洲国产一区二区在线观看| 国产精品无大码| 麻豆精品久久久久久蜜桃| 欧美+亚洲+日韩+国产| 无人区码免费观看不卡| 国产三级在线视频| 亚洲欧美精品综合久久99| 亚洲精品成人久久久久久| 女人被狂操c到高潮| 亚洲国产精品sss在线观看| 夜夜夜夜夜久久久久| 午夜福利18| 色视频www国产| 亚洲自拍偷在线| 两人在一起打扑克的视频| 神马国产精品三级电影在线观看| 成人鲁丝片一二三区免费| 精品久久久久久,| 成人二区视频| 精品国产三级普通话版| 国产又黄又爽又无遮挡在线| 国产成年人精品一区二区| 又黄又爽又刺激的免费视频.| 免费电影在线观看免费观看| 亚洲av成人精品一区久久| 波多野结衣高清作品| 精品国内亚洲2022精品成人| 日韩强制内射视频| 欧美日韩国产亚洲二区| 嫩草影院新地址| 日韩精品中文字幕看吧| 日韩欧美精品免费久久| 国产免费一级a男人的天堂| 午夜爱爱视频在线播放| 日本免费一区二区三区高清不卡| 久久精品国产亚洲网站| 国产男靠女视频免费网站| 久久久国产成人精品二区| 国产精品无大码| 老司机福利观看| 草草在线视频免费看| 日韩欧美国产在线观看| 久久精品影院6| 午夜视频国产福利| 国产伦人伦偷精品视频| 波野结衣二区三区在线| 在现免费观看毛片| 国产高潮美女av| 中文字幕精品亚洲无线码一区| 最好的美女福利视频网| 亚洲avbb在线观看| 亚州av有码| 亚洲人与动物交配视频| 日韩欧美精品v在线| 成熟少妇高潮喷水视频| 国产91精品成人一区二区三区| 欧洲精品卡2卡3卡4卡5卡区| 国产黄色小视频在线观看| 久久久精品大字幕| 久久久久久国产a免费观看| 91在线观看av| 日韩av在线大香蕉| 久久精品91蜜桃| 人妻少妇偷人精品九色| 亚洲成人精品中文字幕电影| 久久精品国产亚洲av天美| 色综合亚洲欧美另类图片| 久久久久久久亚洲中文字幕| 成人特级av手机在线观看| 亚洲成av人片在线播放无| 嫩草影院新地址| av福利片在线观看| 日本欧美国产在线视频| 婷婷精品国产亚洲av在线| 91精品国产九色| 搡老岳熟女国产| 如何舔出高潮| 精品乱码久久久久久99久播| 小蜜桃在线观看免费完整版高清| 国产黄色小视频在线观看| 亚洲经典国产精华液单| 欧美色视频一区免费| 成人特级黄色片久久久久久久| 免费在线观看成人毛片| 亚洲av中文字字幕乱码综合| 亚洲自偷自拍三级| 亚洲va在线va天堂va国产| 在线观看66精品国产| 国国产精品蜜臀av免费| 精品国内亚洲2022精品成人| 亚洲欧美日韩无卡精品| 国产欧美日韩一区二区精品| 欧美最新免费一区二区三区| 少妇裸体淫交视频免费看高清| 日韩一区二区视频免费看| 婷婷六月久久综合丁香| 精品无人区乱码1区二区| 国产极品精品免费视频能看的| a级毛片a级免费在线| 黄色欧美视频在线观看| 国产精品久久久久久精品电影| 女人被狂操c到高潮| 国产男靠女视频免费网站| 一进一出抽搐动态| 亚洲成人久久性| 性色avwww在线观看| 亚洲最大成人中文| 中文字幕熟女人妻在线| 小蜜桃在线观看免费完整版高清| 国产精品三级大全| 欧美最黄视频在线播放免费| 美女高潮喷水抽搐中文字幕| 国产 一区 欧美 日韩| 午夜a级毛片| 中文字幕免费在线视频6| 日韩欧美精品v在线| 性色avwww在线观看| 91狼人影院| a级毛片a级免费在线| 亚洲男人的天堂狠狠| 两个人的视频大全免费| 国产精华一区二区三区| 波野结衣二区三区在线| 毛片一级片免费看久久久久 | 99久久精品热视频| 美女被艹到高潮喷水动态| 搡老妇女老女人老熟妇| 熟女电影av网| 97超级碰碰碰精品色视频在线观看| 国产精品av视频在线免费观看| 免费在线观看日本一区| 亚洲七黄色美女视频| 亚洲第一区二区三区不卡| 日本撒尿小便嘘嘘汇集6| 草草在线视频免费看| 国产欧美日韩精品一区二区| 女人十人毛片免费观看3o分钟| 午夜福利在线观看免费完整高清在 | 国产av不卡久久| 久久精品国产亚洲av涩爱 | 91久久精品国产一区二区成人| 国产一区二区三区视频了| 干丝袜人妻中文字幕| 在线免费观看的www视频| 国产精品不卡视频一区二区| 亚洲精品影视一区二区三区av| 亚洲成人久久性| 如何舔出高潮| 哪里可以看免费的av片| 两性午夜刺激爽爽歪歪视频在线观看| 国产亚洲av嫩草精品影院| 午夜福利在线在线| 国产在视频线在精品| 亚洲最大成人av| 变态另类成人亚洲欧美熟女| 99热6这里只有精品| 国产不卡一卡二| 欧美一级a爱片免费观看看| 午夜亚洲福利在线播放| 不卡一级毛片| 久久精品国产99精品国产亚洲性色| 久久久久久伊人网av| 亚洲午夜理论影院| a在线观看视频网站| 亚洲av中文av极速乱 | 国产精品日韩av在线免费观看| 婷婷亚洲欧美| 亚洲国产色片| 一区福利在线观看| 少妇的逼水好多| 精品久久久久久久末码| 人人妻,人人澡人人爽秒播| 精品久久久久久久久久久久久| 久久久久久国产a免费观看| 国产白丝娇喘喷水9色精品| 在线免费十八禁| 最新中文字幕久久久久| 赤兔流量卡办理| 欧美xxxx黑人xx丫x性爽| 亚洲精品在线观看二区| 观看免费一级毛片| 偷拍熟女少妇极品色| 日本黄色视频三级网站网址| 老司机午夜福利在线观看视频| 亚州av有码| 日本熟妇午夜| 国产成人av教育| 国产一区二区三区av在线 | 国产av一区在线观看免费| 舔av片在线| 欧美区成人在线视频| 成年免费大片在线观看| 精品午夜福利在线看| 看免费成人av毛片| 一进一出抽搐动态| 久久久久久伊人网av| 91狼人影院| 久久久色成人| 一区福利在线观看| 亚洲午夜理论影院| 国产亚洲av嫩草精品影院| 亚洲成人中文字幕在线播放| av专区在线播放| 午夜福利在线观看吧| 欧美日本亚洲视频在线播放| 精品久久久噜噜| 丝袜美腿在线中文| 波多野结衣高清无吗| 最近最新中文字幕大全电影3| 欧美日韩乱码在线| 一区二区三区高清视频在线| 日本在线视频免费播放| 国产精品女同一区二区软件 | 夜夜夜夜夜久久久久| 深夜精品福利| 99久久九九国产精品国产免费| 午夜激情福利司机影院| 国产精品国产三级国产av玫瑰| 在线播放无遮挡| 日本一二三区视频观看| 日韩中字成人| 亚洲在线观看片| 国产欧美日韩精品亚洲av| 韩国av一区二区三区四区| 狠狠狠狠99中文字幕| 一级黄色大片毛片| 熟女人妻精品中文字幕| 久久99热这里只有精品18| 国产伦一二天堂av在线观看| av女优亚洲男人天堂| www.色视频.com| 国产免费男女视频| 国产美女午夜福利| 丰满人妻一区二区三区视频av| 黄色配什么色好看| 亚洲人成伊人成综合网2020| 久久中文看片网| 又黄又爽又免费观看的视频| 18禁黄网站禁片午夜丰满| 美女黄网站色视频| 色综合色国产| 久久99热6这里只有精品| 天堂动漫精品| 亚洲精品一卡2卡三卡4卡5卡| 欧美3d第一页| 久久久久久国产a免费观看| 欧美高清成人免费视频www| 成人无遮挡网站| 精华霜和精华液先用哪个| 麻豆一二三区av精品| 99久久精品国产国产毛片|