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

    Hydrodynamic interactions of water waves with a group of independently oscillating truncated circular cylinders

    2016-11-04 08:53:18XiaohuiZengMinShiShanlinHuang
    Acta Mechanica Sinica 2016年5期

    Xiaohui Zeng·Min Shi·Shanlin Huang

    ?

    RESEARCH PAPER

    Hydrodynamic interactions of water waves with a group of independently oscillating truncated circular cylinders

    Xiaohui Zeng1·Min Shi1·Shanlin Huang1

    In thisstudy,we examine the waterwave radiation by arrays of truncated circular cylinders.Each cylinder can oscillate independently in any rigid oscillation mode with a prescribed amplitude,including translational and rotational modes such as surge,sway,heave,pitch,roll,and their combinations.Based on the eigenfunction expansion and Graf’s addition theorem for Bessel functions,we developed an analytical method that includes the effects of evanescent modes in order to analyze such arrays of cylinders.To investigate the effects of several influential factors on convergence,our objective is to dramatically reduce the number of tests required and determine the influencing relationshipsbetween truncation number and convergence behavior for different factor combinations.We use the orthogonal test method to fulfill the objective.Lastly,we present our results regarding the effects of evanescent modes on hydrodynamic coefficients.

    Radiation·Water wave·Arrays of truncated circular cylinders oscillating independently·Translational oscillation·Rotational oscillation

    1 Introduction

    In recent years,various offshore floating structures designed for the exploitation of space and ocean resources have drawn increasing attention.These structures mainly comprise a sin-gle circular cylinder or an array of circular cylinders,such as deep-sea platforms,floating bases for wind power devices,or wave power extraction devices with oscillating cylinders[1-6].

    ? Xiaohui Zeng

    zxh@imech.ac.cn

    1Key Laboratory for Mechanics in Fluid Solid Coupling Systems,Institute of Mechanics,Chinese Academy of Sciences,Beijing 100190,China

    Many studies have focused on the hydrodynamic interactions between cylinders in waves,using numerical and analyticalmethods.Forexample,Wang etal.[7]investigated radiation by a group ofcylindersusing a finite-element-based numerical method.Zhou et al.[8]investigated the nonlinear radiated waves generated by a floating structure in forced motion by using a higher-order boundary element method(HOBEM).Babarit[9]investigated the energy absorption problem in both regular and irregular waves,and considered wave interactions in arrays of two surging or heaving cylinders,whose hydrodynamic coefficientswere calculated using the boundary element method(BEM).In addition,Williams et al.[10,11]used a modified plane-wave approach to investigate the diffraction and radiation problems associated with an array of truncated circular cylinders.Furthermore,Kagemoto and Yue[12]presented an accurate algebraic method(including the interactions of evanescent waves)for multiple bodies by combining and extending the matrix method and the multiple-scattering concept,based on the diffraction characteristics of the isolated member.Also,Linton and Evans[13]developed an exact theory for wave diffraction by vertical circular cylinders and achieved a major simplification.Additionally,Yilmaz[14]obtained analytical solutions for the wave scattering and radiation of truncated cylinders using the accurate algebraic method developed by Kagemoto and Yue[12].Moreover,Siddorn and Eatock Taylor[15]also used an exact algebraic method to study radiation and diffraction in an array of cylinders oscillating independently.Then,Chatjigeorgiou[16]investigated wave diffraction by arrays of truncated elliptical cylin-ders subjected to regular waves and presented an analytical solution.and Guillot[22]examined the feasibility ofan intelligentsensor fusion to estimate on-line surface finish and dimensional deviations with the orthogonal arrays used for experiments planning.In addition,Green et al.[23]applied orthogonal arrays to conjoint analysis in marketing.Also,Lee et al.[24]used an iterative optimization algorithmwith orthogonal arrays in discrete space for structure design.The feasibility of a method that combines orthogonal experiment theory with the finite element method(FEM)was discussed by Wang et al.[25],who used this method to compute and evaluate the three dimensional(3-D)natural stress field.In view of the high efficiency of orthogonal design in reducing the number of required tests or computations,we use this method in this study to investigate the effect of different influentialfactors(e.g.,oscillation mode and waterdepth)on convergence behaviors and to identify primary and secondary factors.Using the orthogonal test method,with a relatively small amount of calculation,we can identify how the numberofitemsretained in the expansion affects the convergence behavior in circumstances with different parameter combinations.Accordingly,we can then develop more reasonable calculation procedures for other cases.

    In this study,we classify arrays of floating cylinders into two categories.The first is an“array of cylinders with no relative motion”(hereafter referred to as the“C1 array”)because the connection between the cylinders is so strong that the cylinders can be regarded as being rigidly connected as a whole.The second array category is an“array of cylinders with relative motions”(hereafter referred to as the“C2 array”)because there isno rigid connection between the component members and each component member can oscillate independently.
    Studies of the hydrodynamic interactions between cylinders have mainly focused on the“C1 array,”and there are relatively few studies of the“C2 array”(especially using analytical methods).Published results for the“C2 array”are mainly concerned with concrete examples in which each cylinder oscillates in translational modes(i.e.,surge,sway,heave)rather than there being an array of cylinders oscillating independently in rotational modes(i.e.,roll or pitch),let alone in multiple modes that include both translation and rotation.Asmore and more“C2 array”characteristicsemerge in new offshore structures,further investigation is merited regarding the hydrodynamic interactions of“C2 array”oscillations with multiple modes,including both translation and rotation.Therefore,in this study,we investigate water wave radiation and diffraction by the“C2 array”in the context of linear wave theory.
    With respect to computation,the infinite series in equations must be truncated,so it is important to determine the appropriate number of terms retained in the expansion after truncation with lessthan about1%error.An isolated cylinder does not require much computation to obtain the appropriate truncation number by trial and error.However,it would be onerous and time consuming to do so for an array of cylinders,because the hydrodynamic problems are much more complex and there are more factors influencing convergence behavior.For the same reasons,it is very hard to obtain explicit mathematical expressions to analyze convergence behavior.If convergence were to be examined by trial and error for cases involving different factors,the costs would be prohibitive.Therefore,another alternative must be sought to determine the relationships between different influential factors and convergence behaviors.In this study,we adopted the orthogonal test method to determine these relationships.
    The orthogonal test(experimental design)method is an efficient testing strategy for identifying the relationships between influential factors and the test index,as well as for differentiating the primary from the secondary factors in fewer trials.Since the 1940s,this method has been adopted by industry and academia[17-21],and has found worldwide applications.There are many examplesofthe successful application of the orthogonal design.For example,Azouzi

    2 Solutions for velocity potentials

    Due to the axial symmetry of a circular cylinder,the yaw motion introduces no hydrodynamic forces in a perfect fluid. Therefore,only five modes-surge,sway,heave,roll,and pitch-are considered in this study.

    The radius and draught of an isolated truncated cylinder are a and h,respectively,and the water depth is d.The fluid is divided into two regions(core and exterior),as depicted in Fig.1.

    Figure 2 shows the six oscillation modes of a single truncated circular cylinder.

    The array(Fig.3)consists of N verticaltruncated circular cylinders.We employ a right-handed Cartesian space fixed coordinate system XOY with the plane xoy coinciding with the still-water level(SWL),and the positive z-axis pointing vertically upwards.The coordinate of the origin Oiof each cylinderin the overallcoordinate system is(xi,yi).The radius of each cylinder is ai,i=1,2,...,N.The local cylin-drical coordinate system Oiriθiz is centered on the origin of each cylinder.The coordinate of the origin Ojof cylinder j in the cylindrical coordinate system centered on the origin Oiof cylinder i is(Lij,αij,z),i,j=1,2,...,N.

    Fig.1 Sketch of a single truncated circular cylinder

    Fig.2 Array of truncated cylinders and six oscillation modes

    Fig.3 Local cylindrical coordinate system for an array of truncated circular cylinders

    Assuming that the fluid is inviscid and incompressible and that the flow is irrotational,we can employ the potential theory.We can consider the linear steady-state oscillations since the wave steepness and amplitudes of motion are small enough.Therefore,physical quantities such as the velocity potentials and displacements are periodic functions of time with angular frequency ω0.Then,the total velocity potential can be written as the product of the space and time factors,as follows

    The potential ?RDfor an array of cylinders with independent prescribed motions is the sum of the radiation and scattering potentials,which can be written as

    The velocity potential in a linear water wave problem can be decomposed into the incident potential(due to environmental waves or other sources,such as waves emanating from other cylinders due to wave diffraction or radiation),diffraction potential(fixed cylinder but an incident wave),and radiation potential(oscillating cylinder without incident wave).For the problem considered in this study,is the diffraction potentialofcylinder j;isthe radiation potential of cylinder j;with-the sum ofthe potentialofwaves emanating from othercylinders(i=1,2,...,N,except cylinder j)due to wave diffraction and radiation-acting as the incident potential on cylinder j.

    The instantaneous displacement of cylinder i in smallamplitude periodic oscillation with angular frequency ω0is

    where ζ(i)srepresents the oscillation amplitude of mode s for cylinder i.The corresponding velocity is

    When cylinder i oscillates independently with amplitude ζ(i)sin mode s,the corresponding radiation potential is

    According to the impermeable condition on the body surface,the normal velocity of the cylinder is equal to that ofthe adjacent fluid,and the equation is written as

    Suppose the coordinate of the centroid of the cylinder isAs shown in Fig.1,the fluid domain is divided into two regions:the core region and exteriorregion.The radiation potentialsfor the exterior region and thefor the core region satisfy the impermeability conditions on the body surface,as follows

    Apart from the body surface conditions,the velocity potentials for the exterior and core regions must satisfy Laplace’s equation,the free-surface condition,the seabed condition,and the radiation condition atinfinity,respectively,as follows

    where Hm=Jm+iYmis an m-th-order Hankel function of the first kind.Jmand Ymare m-th-order Bessel functions of the first and the second kinds,respectively.Kmis a modified Bessel function of the second kind.

    The wave number k0and the angular frequency ω0satisfy the usual dispersion relation

    The eigenvalues kn(n=1,2,…)are the positive real roots of the following equations

    In Eq.(10),Z0(z)and Zn(z)are expressed asin Eq.(13)can be solved using Eq.(A1)in the appendix

    Equation(10)can be expressed in matrix form as

    where“D-E”is short for“diffraction-exterior region”andrepresents the column vector of the partial radiation waves in the exterior region of cylinder i(except for the amplitude,the function form ofradiation and diffraction partial waves are the same in the exterior region,so we adopt the notationfor both radiation and diffraction partial waves in this paper).

    Using Graf’s addition theorem for the Bessel functions(valid for rj<Rij)and replacing l by-l,Eq.(15)is rewritten as

    where Ilis the l-th-ordermodified Besselfunction ofthe first kind and Eq.(16)can be expressed in matrix form as

    ψIjis an n0(2l0+1)column vector representing incident partial-wave functions written as

    The diffraction potentialofthe exteriorregion in the vicinity of cylinder i,which is generated from radiation waves emanating from other cylinders,is

    and its matrix form is

    where ARiis an unknown column vector of n0(2m0+1) elements,and each element is,representing the amplitudes of the partial diffraction waves in the exterior region caused by cylinder oscillations.

    The total incident potential in the vicinity of cylinder j is the sum ofthe radiation potentials caused by the independent oscillation ofother N-1 cylinders and the diffraction potentials due to the hydrodynamic interactionsbetween cylinders,as follows

    For cylinder j,the diffraction potential(Eq.(21))and the totalincidentpotential?Ij(Eq.(22))can be connected by the diffraction transfer matrixin the exterior region as

    Given the prescribed amplitudeand the knownand Risvalues,the unknown vector ARj(n,m)can be solved from Eq.(23).

    Thus,the total exterior potential in the vicinity of cylinder j as each cylinder oscillates in different modes(s= 1,2,3,4,5)is obtained as

    The total core potential of cylinder j is

    C0mand Cnmin Eq.(26)are obtained from Eq.(A1)in the appendix.in Eq.(25)is the column vector of the partial-wave functions in the core regionin Eq. (25) is the diffraction transfer matrix for cylinder j in its core region, representing the amplitudes of the partial diffraction waves caused by unit-amplitude partial incident waves,which can also be found in the appendix.

    Therefore,we obtain the velocity potential adjacent to cylinder j as each cylinder simultaneously oscillates with independent amplitudeIn Eq.(24),is the velocity potential for the exterior region in the vicinity of cylinder j.In Eq.(25),is the velocity potential for the core region of cylinder j.

    3 Hydrodynamic forces and moments

    Then,we can obtain the horizontal hydrodynamic forceson cylinder j as each cylinder oscillates in different modes

    The verticalhydrodynamic forceon cylinder j is

    where Sbis the side surface of a cylinder,with unit normal vector written as=(-cosθj,-sinθj,0);Sbhis the bottom surface,with unit normal vector=(0,0,1).

    4 Comparisons and verifications

    Fig.4 1-4 heave-heave and heave-surge.a Real parts.b Imaginary parts

    Fig.5 Comparisons of hydrodynamic coefficients in pitch mode between present study and that of Ali and Khalil[26].a Added mass.b Damping coefficients

    Based on the method developed in the preceding sections,we wrote a computer program and then verified the program priorto performing furthercalculationsforthe hydrodynamic interactions between truncated circular cylinders in an array with some oscillating cylinders.

    Note that Siddorn and Eatock Taylor[15]provided an example of the hydrodynamic interactions of an array of cylinders caused by one cylinder oscillating in translational mode with a prescribed amplitude.We then compared the results of our study with those of Siddorn and Eatock Taylor(Fig.14 in Ref.[15]),as shown in Fig.4a and b,in which the circularand square points denote the results of Siddorn and Eatock Taylor[15],and the solid and dashed-dotted lines denote the results of our study.Using the nomenclature adopted by Siddorn and Eatock Taylor[15],i-j designates interactions between cylinders i and j.Forexample,the label 1-4 heave-surge refers to the heave force on cylinder 1 due to unit surge velocity of cylinder 4.We observe good agreement,as the results of our study and Ref.[15]are almost identical,with maximum differences of 0.2%.

    We also compared the results from this study with those from previous studies using the numerical method(3-D source-sink method)from Ali and Khalil[26].For the cases shown in Fig.5a and b,the square and triangular points denote the results of Ali and Khalil[26]and the solid lines denote the results of this study.Comparisons of the results of this study with those of Ali and Khalil(Fig.14 in Ref.[26])are shown in Fig.5a and b,respectively.For the cases shown in Fig.5a and b,there are three cylinders,andrefers to the pitch moment on cylinder 2 due to the unit pitch velocity of cylinder 1.We can see that the variation trends of the present and reference studies are identical,with the majority of differences in the results ranging from 0.2%to 9.9%.The numerical 3-D source-sink method was adopted in Ref.[26]and the accuracy of their results is dependent on many factors(such as the mesh discretization).The results of our study agree well with those of Ref.[26].

    5 Investigations of truncation number using orthogonal test method

    For an array of truncated cylinders,there are many factors that can affect convergence behavior,such as column number,column spacing,oscillation mode,water depth,draught,and wave number,and the appropriate truncated number(corresponding to convergent results)for different factorcombinations also differs.Itis very difficultto obtain explicit mathematical expressions suitable for an analysis of convergence.The trial-and-error examination of convergence for cases with different factor combinations would involve prohibitive costs.Therefore,we introduce the orthogonal test method in this study,which investigates(with a minimum number of tests)the relationships of influence between different influential factors and convergence behaviors,and identifies the primary and secondary factors.In this section,we present concrete examples for an array consisting of four cylinders placed at the vertices of a square.

    For this study,we use the mixed-level orthogonal array L50(510×101),where L denotes a Latin square,50 is the numberoflevelcombinations,510indicates thatthere are ten factors and five levels per factor,and 101indicates that there is one factor that has ten levels.

    After carrying out the tests according to the orthogonal array,we analyzed the data using a direct analytical method known as range analysis.The range(generally denoted by R)is the maximum average numerical difference of each factorateach level.Because the factors have differentlevels,a conversion coefficient is needed(generally denoted by d). Then the formula of the revised range(denoted by R′A)is,where rAisthe frequency ofoccurrence of factor A in one column according to each level.Thus the primary and secondary factors can be analyzed by comparing the sizes of R′A.Table 1 shows the conversion coefficient d for different levels.

    The influence of differentfactors(oscillation mode,water depth,draught,column spacing,and wave number)and theircombinationson convergence mustbe comprehensively investigated.In the array,the wave numberhas ten levels and other factors have five,respectively.Table 2 shows the factors,which are dimensionless,as well as their levels.

    Once the factorsand theirlevels in the corresponding position are filled in the orthogonal array,then the calculation is carried outaccording to the orthogonalarray(here,we selectcolumns 1,2,4,7,and 11 in the orthogonal array).We note that to investigate the convergence of the results,we compare the results solved from smallertruncation items(m=6,n=25,p=25)with those from large enough truncation items(m=6,n=60,p=60).Before introducing the orthogonal test method to determine the appropriate truncated number,we performed many calculations by trial and error,and the resultsshow thatthere is almostno difference if the value of m is greaterthan 6.The maximum errorbetween m=6 and m=7 is about 0.0002%in the region of interest,so 6 is a sufficient value for m.The orthogonal array and calculation results are listed in Table 3.With the same truncation items,cases with large percentage errors require moretruncation items to achieve the same convergence accuracy,which means a slower convergence rate.

    Table 1 Conversion coefficients

    Table 2 Factors and their levels

    Table 3 Orthogonal array and calculation results

    Table 3 continued

    Fig.6 Index fluctuation when k0a changes

    We drew the following diagrams after processing the data in Table 3,since it is easier to observe from a diagram the index fluctuation when each factor’s level changes.

    Figures 6-10 illustrate the mean error of each factor over the range ofcorresponding levels,respectively.Taking Fig.6 as an example,itwas drawn by calculating the mean value of the results(percentage error)at each level in the k0a column of Table 3.The above figures show that the convergence rate of a case with smaller wave numbers is significantly faster than that with larger wave numbers;the convergence rate ofthe case in surge or sway mode is faster than those in heave and rotational modes;the convergence rate of a case with a column spacing of 15a is much faster than cases with other column spacings in the defined range;the convergence rate increases significantly by increasing the draught-to-waterdepth ratio.

    Fig.7 Index fluctuation when a ii changes

    Fig.8 Index fluctuation when L/a changes

    Fig.9 Index fluctuation when d/a changes

    Fig.10 Index fluctuation when h/d changes

    Fig.11 Primary and secondary factors

    The range R of each factor can be obtained by processing the data in Figs.6-10,which helps to further analyze the primary and secondary factors.Since the wave-number level differsfromthatoftheothers,wecalculatedtherevisedrange R′to compare each factor’s effect on the index.

    Figure 11 indicates that both the oscillation mode and the ratio of the draught to water depth have greater influence on the convergence rate than wave number and column spacing do.

    Next,we designed a test case in which four cylinders are located at the vertexes of a square,with column spacing 5a,and water depth 20a,to verify the rule given above.First,we investigated the surge and pitch modes for a ratio of draught to water depth of 0.1,and then for the surge mode for ratios of 0.1 and 0.6.

    WecanobservefromFig.12a-cthattheoscillationmodes andthedraught-to-water-depthratiohavearemarkableinfluence on the convergence rate;the convergence rate of the surge mode is faster than the pitch mode;the convergence rates of cases in which the draught-to-water-depth ratio is 0.6 are faster than those with the ratio 0.1.

    Thus far,there is no effective way to properly estimate the truncated indices apart from trial and error.The above rule simplifies the investigation of convergence in the cases discussedbelow.Basedontheconclusionsdrawninthissection,cases with rotational modes or smaller ratios of draught to water depth should be given priority when investigating truncated features.If the value of some truncated items in the above cases satisfies the precision requirement,then the same value will satisfy cases with higher convergence rates,such as those with translational modes or larger ratios of the draught to water depth.As such,there is no need to perform trials for every case with different parameter combinations,which improves efficiency and saves a lot of time.

    6 Results for several cases and discussion

    6.1Hydrodynamic coefficients:the effects of evanescent modes

    As mentioned above,there are some solutions available for arrays of circular cylinders that use large-spacing approximation,neglecting the evanescent modes in the velocity potential of emanating waves.In this study,the solutions include evanescent modes.Therefore,the effects of evanescent modes on the hydrodynamic interactions can be evaluated.In this subsection,we use two different algorithms to calculate the hydrodynamic forces of an array of cylinders: one includes evanescent modes and the other does not.Then we demonstrate the effects of evanescent modes by comparing the two sets of results.

    6.1.1Two cylinders

    First,we recalculated the hydrodynamic coefficients for two cylinders(case 1 in Williams and Abul-Azm[11],where cylinder 1 oscillates and cylinder 2 is fixed).Then,we compared the two sets of results(including and without evanescent modes)with that using a modified plane-wave techniquebasedonalarge-spacingapproximation(Williams and Abul-Azm[11]),as shown in Figs.13 and 14.

    In Figs.13a,b and 14a,b,the triangular points denote the results of Williams and Abul-Azm[11],the solid lines denotetheresultsofthisstudythatincludeevanescentmodes,and the dashed lines denote the results of this study without evanescent modes.refers to the real/imaginary partof the hydrodynamic force in the k-direction on cylinder i due to oscillation in the s-mode of cylinder j.

    Fig.12 Test case results.a Oscillation modes of surge and pitch.b h/d of 0.1 and 0.6.c Oscillation modes of surge and pitch,h/d of 0.1 and 0.6

    Fig.13 Comparisons ofhydrodynamic coefficients on cylinder1 between the presentstudy(including and withoutevanescentmodes)and Williams and Abul-Azm[11],L=3a.a Added mass.b Damping coefficients

    We can see that the results without evanescent modes are in good agreement with those of Williams and Abul-Azm[11],which are based on a large-spacing approximation. For the damping coefficient,the two sets of results(including and without evanescent modes)in this study nearly coincide with each other.However,there are evident differences in the solutions of added mass for results including or without evanescent modes,especially for added mass on cylinder 2,which reveals the effects of evanescent modes.

    6.1.2Array consisting of five circular cylinders:four in a straight line,the other on a perpendicular bisector of the line

    Next,we investigated a cylinder array consisting of five truncated circularcylinders with a specialgeometry configuration,as shown in Fig.15.The configuration of the cylinder array is symmetric about the x-axis.The shaded and hollow circles in Fig.15 denote oscillating and fixed cylinders,respectively.

    Fig.14 Comparisonsofhydrodynamic coefficientson cylinder2 between the presentstudy(including and withoutevanescentmodes)and Williams and Abul-Azm[11]L=3a.a Added mass.b Damping coefficients

    Fig.15 Geometric configuration of cylinders and oscillation forms: cylinder 1 oscillates while others remain still

    In Fig.15,we see thatcylinder 1 oscillates in surge,heave or pitch mode,and the others remain still.We calculated the hydrodynamic coefficients of cylinders 1 and 3,including and without evanescent modes,with a water depth d= 20a,draught h=10a,and column spacings L=2.4a,3a,5a.The results are shown in Figs.16-21,where“cylinder 1/3-x/z/pitch-direction(cylinder 1 surge/heave/pitch and 2,3,4,5 fixed)”represents the hydrodynamic coefficients of cylinder 1/3 in the x/z/pitch-direction as cylinder 1 oscillates in surge/heave/pitch mode while othercylinders remain still.

    In Figs.16-21,the solid and hollow symbols,respectively,representthe two setsofresults(including and without evanescent modes)for the hydrodynamic coefficients of cylinder 1 or 3 for the configurations shown in Fig.15 with column spacings L=2.4a,3a,and 5a.Figures 16a,b and 17a,b,respectively,show the added mass and damping coefficients of cylinders 1 and 3 in the x-direction as cylinder 1 oscillates in surge mode.Figures 18a,b and 19a,b,respectively,show the added mass and damping coefficients of cylinders 1 and 3 in the z-direction as cylinder 1 oscillates in heave mode.Figures 20a,b and 21a,b,respectively,show the added mass and damping coefficients of cylinders 1 and 3 in the pitch-direction as cylinder 1 oscillates in pitch mode.

    Fig.16 Hydrodynamic coefficients of cylinder 1-x-direction(cylinder 1 surge and 2,3,4,5 fixed).a Added mass.b Damping coefficients

    Fig.17 Hydrodynamic coefficients of cylinder 3-x-direction(cylinder 1 surge and 2,3,4,5 fixed).a Added mass.b Damping coefficients

    Fig.18 Hydrodynamic coefficients of cylinder 1-z-direction(cylinder 1 heave and 2,3,4,5 fixed).a Added mass.b Damping coefficients

    Fig.19 Hydrodynamic coefficients of cylinder 3-z-direction(cylinder 1 heave and 2,3,4,5 fixed).a Added mass.b Damping coefficients

    We can see from Figs.16-21 that there are evident differences between the two sets of solutions for added mass(including or without evanescent modes),especially for cylinder 3.The differences increase with a decrease in column spacing.This reveals the effects of the evanescent modes.When the column spacing is L=2.4a,when the cylinder is oscillating in surge mode,the maximum difference between the two sets ofsolutions forthe added mass for cylinder1 is23%,and the maximum difference forcylinder3 is 147%.Forpitch mode(with the same column spacing),the differences are even greater than those of the surge or heave mode.It appears that the effects of the evanescent modes are more pronounced for cases in which the cylinder oscillates in rotational mode than in translationalmode;and these effects are more significant for a fixed cylinder.When the cylinders are located close together,no matter which oscillation mode is considered,the effect of evanescent modes on the added mass is significant and should not be neglected. However,even if the cylinders are located close together,there is little effect by the evanescent modes on the damping coefficients.

    Fig.20 Hydrodynamic coefficients of cylinder 1-pitch-direction(cylinder 1 pitch and 2,3,4,5 fixed).a Added mass.b Damping coefficients

    Fig.21 Hydrodynamic coefficients of cylinder 3-pitch-direction(cylinder 1 pitch and 2,3,4,5 fixed).a Added mass.b Damping coefficients

    6.2Array consisting of six circular cylinders:five arranged at the vertices of a regular pentagon inscribed in a circle and the sixth at the center of the circle Another geometric configuration investigated in this paper is an array of six truncated circular cylinders,five of which are arranged at the vertices of a regular pentagon inscribed in a circle and the sixth at the center of the circle,as shown in Fig.22,with d=20a,h=10a,and L=5a.

    For this situation,cylinder 1 oscillates in surge mode,cylinders 2 and 5 in rolland pitch modes,cylinders 3 and 4 in heave mode,and cylinder6 remainsstill.Cylinder2 oscillates in phase(in antiphase)with cylinder 5 in the pitch(roll)-direction.The prescribed amplitudes of the oscillation are

    For ζ(j)s,j represents the cylinder number and s the oscillation mode,with j=1,2,3,4,5,and s=1,2,3,4,5.

    Fig.22 Geometric configuration of six cylinders

    The hydrodynamic pressure distributions on each cylinder at different depths are shown in Fig.23a-f,where A represents wave amplitude.We can see that the pressure distributions of cylinder 2(3)and those of cylinder 5(4)are symmetric with respect to the x-axis because of the symmetry of the geometric configuration and oscillation modes. For the same reason of symmetry,the curves of the pressure distributions on cylinder 1(6)are self-symmetric about the x-axis.

    Fig.23 Pressure distributions on each cylinder for oscillation situation in subsection 6.2.a Cylinder 1.b Cylinder 2.c Cylinder 3.d Cylinder 4. e Cylinder 5.f Cylinder 6

    7 Conclusions

    In this study,we investigated the hydrodynamic interactions ofarrays oftruncated circularcylinders with relative motions(“C2 arrays”)in which each cylinder can oscillate independently in any rigid oscillation mode(including translational and rotational modes,such as surge,sway,heave,pitch,roll,and their combinations)with an arbitrarily prescribed amplitude,including evanescent modes.In addition,we used the orthogonal test method to investigate the effects of several factors influencing convergence in different combinations,and assessed the significances of these different factors.The results illustrate the effects of evanescent modes on hydrodynamic coefficients.Our conclusions are summarized as follows

    (1)The calculation results reveal that in the hydrodynamic features of an array of cylinders,as each cylinder oscillates in various modes,there are significant differences in values and variation rules.It appears that the effects of evanescent modes are more pronounced for cases in which the cylinder oscillates in rotational mode than those that oscillate in translational mode;such effects are more pronounced for the fixed cylinder in the array. When the cylinders are located in close proximity,no matter which oscillation mode is considered,the effect ofthe evanescentmodes on the added mass is significant and should not be neglected.However,even if the cylinders are located in close proximity,there is little effect of the evanescent modes on the damping coefficients.

    (2)We employed the orthogonal test method in this study,which can dramatically reduce the number of required tests,and we obtain the relationships between different influential factors and convergence behaviors for different factor combinations.The main results can be summarized as follows:the oscillation mode and draught-to-water depth ratio have more significant effects on the convergence rate than wave number and column spacing do;the convergence rates of cases in surge or sway mode are faster than those in heave and rotational modes;the convergence rate increases with a rise in the draught-to-water-depth ratio.The above resultscan play an importantrole in improving efficiency in obtaining convergentresults when analyzing arrays of circular cylinders in various arrangements and oscillation modes.

    Acknowledgments The project was supported by the National Natural Science Foundation of China(Grants 11072246,51490673)and the National Basic Research Program(973 Program)of China(Grant 2014CB046801).

    Appendix

    The unknown radiation coefficientscan be solved from the equation set

    The expressions of the known coefficients in Eq.(A1a,A1b)are

    The diffraction transfer matrices of the truncated single cylindersare obtained following the proceduregiven in Kagemoto and Yue[12],the elements of(or)are the amplitude ofthe q-th(or p-th)partialwave ofthe diffraction potential due to a single unit-amplitude incidence of mode n on cylinder j,which are listed below

    The element Tij(n,m,l)in Eq.(A7d)is as follows

    1.Falc?o,A.F.O.:Wave energy utilization:A review of the technologies.Renew.Sustain.Energy Rev.14,899-918(2010)

    2.Eriksson,M.,Waters,R.,Svensson,O.,et al.:Wave power absorption:Experiments in open sea and simulation.J.Appl.Phys.102,084910(2007)

    3.Dong,Y.-Q.:Wave Loads and Response of the Oil-extraction Platform in Deep Ocean.Tianjin University Press,Tianjin(2005)(in Chinese)

    4.Zeng,X.-H.,Li,X.-W.,Liu,Y.,etal.:Nonlineardynamic responses oftension leg platform with slack-tauttether.China Ocean Eng.23,37-48(2009)

    5.Zeng,X.-H.,Shen,X.-P.,Wu,Y.-X.:Governing equations and numerical solutions of tension leg platform with finite amplitude motion.Appl.Math.Mech.Engl.Edition 28,37-49(2007)

    6.Zeng,X.-H.,Yu,Y.,Zhang,L.,et al.:A new energy-absorbing device formotion suppression in deep-sea floating platforms.Energies 8,111-132(2015)

    7.Wang,C.-Z.,Mitra,S.,Huang,H.-C.,etal.:Finite elementanalysis of second order wave radiation by a group of cylinders in the time domain.J.Hydrodyn.25,348-361(2013)

    8.Zhou,B.-Z.,Ning,D.-Z.,Teng,B.,et al.:Fully nonlinear modeling of radiated waves generated by floating flared structures.Acta Mech.Sin.30,667-680(2014)

    9.Babarit,A.:Impact of long separating distances on the energy production of two interacting wave energy converters.Ocean Eng.37,718-729(2010)

    10.Williams,A.N.,Demirbilek,Z.:Hydrodynamic interactions in floating cylinder arrays-I:Wave scattering.Ocean Eng.15,549-583(1988)

    11.Williams,A.N.,Abul-Azm,A.G.:Hydrodynamic interactions in floating cylinder arrays-II:Wave radiation.Ocean Eng.16,217-263(1989)

    12.Kagemoto,H.,Yue,D.K.P.:Interactions among multiple threedimensional bodies in water waves:an exact algebraic method.J. Fluid Mech.166,189-209(1986)

    13.Linton,C.M.,Evans,D.V.:The interaction of waves with arrays of vertical circular cylinders.J.Fluid Mech.215,549-569(1990)

    14.Yilmaz,O.:Hydrodynamic interactions of waves with group of truncated vertical cylinders.J.Waterw.Port Coast.Ocean Eng. 124,272-279(1998)

    15.Siddorn,P.,Eatock Taylor,R.:Diffraction and independent radiation by an array of floating cylinders.Ocean Eng.35,1289-1303(2008)

    16.Chatjigeorgiou,I.K.:The hydrodynamics of arrays of truncated elliptical cylinders.Eur.J.Mech.B/Fluids.37,153-164(2013)

    17.Rao,C.R.:Factorial experiments derivable from combinatorial arrangements of arrays.R.Stat.Soc.(Suppl.)9,128-139(1947)

    18.Taguchi,G.:Performance analysis design.Int.J.Prod.Res.16,521-530(1978)

    19.The Orthogonal Test Method Authoring Group:Orthogonal Test Method.National Defence Industry Press,Beijing(1976)

    20.Fang,K.-T.,Ma,C.-X.:Orthogonal and Uniform Test Design.Sciense Press,Beijing(2001)

    21.Hedayat,A.S.,Sloane,N.J.A.,Stufken,J.:OrthogonalArrays:Theory and Application.Springer,New York(1999)

    22.Azouzi,R.,Guillot,M.:On-line prediction of surface finish and dimensionaldeviation in turning using neuralnetwork based sensor fusion.Int.J.Mach.Tools Manuf.37,1201-1217(1997)

    23.Green,P.E.,Krieger,A.M.,Wind,Y.:Thirty yearsofconjointanalysis:reflections and prospects.Interfaces 31,56-73(2001)

    24.Lee,K.H.,Yi,J.W.,Park,J.S.,et al.:An optimization algorithm using orthogonal arrays in discrete design space for structures. Finite Elem.Anal.Des.40,121-135(2003)

    25.Wang,T.,Zhou,X.-Q.,Tian,S.-B.,et al.:Numerical simulation method for rock natural stress field of a valley and its application based on orthogonal experiments.Rock Soil Mech.24,831-835(2003)

    26.Ali,M.T.,Khalil,G.M.:On hydrodynamic interaction between several freely floating vertical cylinders in waves,In:Proceeding of the 24th International Conference on Offshore Mechanics and Arctic Engineering(2005)

    27.Yilmaz,O.,Incecik,A.,Barltrop,N.:Wave enhancement due to blockage in semisubmersible and TLP structures.Ocean Eng.28,471-490(2001)

    7 September 2015/Revised:29 January 2016/Accepted:10 March 2016/Published online:20 July 2016

    ?The Chinese Society of Theoretical and Applied Mechanics;Institute of Mechanics,Chinese Academy of Sciences and Springer-Verlag Berlin Heidelberg 2016

    国产精品久久久久久精品电影| 亚洲aⅴ乱码一区二区在线播放| 夫妻午夜视频| 国产一区二区在线观看日韩| 人妻制服诱惑在线中文字幕| 波多野结衣巨乳人妻| 成人亚洲精品av一区二区| 麻豆乱淫一区二区| 国产美女午夜福利| 毛片女人毛片| 春色校园在线视频观看| 国语对白做爰xxxⅹ性视频网站| 亚洲内射少妇av| 精品国产乱码久久久久久小说| 26uuu在线亚洲综合色| 99热网站在线观看| 国产精品三级大全| 有码 亚洲区| 久久综合国产亚洲精品| 午夜老司机福利剧场| 丝袜脚勾引网站| 亚洲精品中文字幕在线视频 | 亚洲精品国产av蜜桃| 国产亚洲精品久久久com| 最近手机中文字幕大全| 一本一本综合久久| 亚洲国产色片| 色婷婷久久久亚洲欧美| 国产老妇伦熟女老妇高清| 午夜爱爱视频在线播放| 亚洲精品第二区| 婷婷色av中文字幕| 人妻一区二区av| 久久久午夜欧美精品| 高清av免费在线| 在线 av 中文字幕| 精品一区二区三卡| 午夜福利高清视频| 国产亚洲午夜精品一区二区久久 | 国产永久视频网站| 精品久久国产蜜桃| 简卡轻食公司| 成人欧美大片| 久久久亚洲精品成人影院| 九色成人免费人妻av| 国产极品天堂在线| 七月丁香在线播放| 日韩免费高清中文字幕av| 伊人久久国产一区二区| 亚洲无线观看免费| 午夜视频国产福利| 高清av免费在线| 亚洲精品国产成人久久av| 男女国产视频网站| 亚洲色图综合在线观看| 欧美成人一区二区免费高清观看| 久久精品夜色国产| 中文字幕免费在线视频6| 欧美 日韩 精品 国产| 国产成人免费观看mmmm| 久久99热这里只有精品18| 七月丁香在线播放| 狂野欧美激情性xxxx在线观看| 在现免费观看毛片| 在线亚洲精品国产二区图片欧美 | 女人被狂操c到高潮| 国产亚洲午夜精品一区二区久久 | 天堂俺去俺来也www色官网| 在线亚洲精品国产二区图片欧美 | 久久女婷五月综合色啪小说 | 免费看a级黄色片| 少妇人妻久久综合中文| h日本视频在线播放| 青春草视频在线免费观看| 成年女人在线观看亚洲视频 | 欧美日韩精品成人综合77777| 久久久久久伊人网av| 天美传媒精品一区二区| 18禁动态无遮挡网站| www.av在线官网国产| 亚洲av成人精品一二三区| 亚洲精品影视一区二区三区av| 久久精品国产亚洲av天美| 久久女婷五月综合色啪小说 | 亚洲国产成人一精品久久久| 伊人久久精品亚洲午夜| 男女下面进入的视频免费午夜| 欧美激情在线99| 亚洲不卡免费看| 老司机影院毛片| 亚洲av男天堂| 久久久色成人| 赤兔流量卡办理| 在线观看av片永久免费下载| 听说在线观看完整版免费高清| 国产视频内射| 九草在线视频观看| 中国三级夫妇交换| 国产淫片久久久久久久久| 99久久精品一区二区三区| 亚洲av国产av综合av卡| av在线天堂中文字幕| 久久精品久久久久久久性| 看免费成人av毛片| 日韩大片免费观看网站| 水蜜桃什么品种好| 国产又色又爽无遮挡免| av免费观看日本| 亚洲国产精品成人综合色| 国产日韩欧美亚洲二区| 成人毛片60女人毛片免费| 97在线视频观看| 色视频www国产| 国产成人精品福利久久| 免费不卡的大黄色大毛片视频在线观看| 日本欧美国产在线视频| 国产精品国产av在线观看| 亚洲自拍偷在线| 久久精品久久精品一区二区三区| 又爽又黄a免费视频| 大话2 男鬼变身卡| 亚洲最大成人av| 91精品国产九色| 一个人看的www免费观看视频| 最近最新中文字幕免费大全7| 国产老妇伦熟女老妇高清| 国产又色又爽无遮挡免| 99热国产这里只有精品6| 一本久久精品| 80岁老熟妇乱子伦牲交| 亚洲天堂av无毛| 丝袜美腿在线中文| videos熟女内射| 午夜福利网站1000一区二区三区| 中国国产av一级| 成人亚洲精品av一区二区| 五月开心婷婷网| 久久99精品国语久久久| 国产黄片美女视频| 一级毛片久久久久久久久女| 精品少妇久久久久久888优播| 在线观看免费高清a一片| 99热全是精品| 国产色婷婷99| 免费av不卡在线播放| av国产免费在线观看| 国产成人福利小说| 婷婷色av中文字幕| 肉色欧美久久久久久久蜜桃 | 在线观看人妻少妇| 2022亚洲国产成人精品| 王馨瑶露胸无遮挡在线观看| 麻豆乱淫一区二区| a级一级毛片免费在线观看| 成人毛片a级毛片在线播放| 激情五月婷婷亚洲| 免费黄网站久久成人精品| 夜夜看夜夜爽夜夜摸| 精品久久国产蜜桃| 在线亚洲精品国产二区图片欧美 | 精品久久久噜噜| 中文天堂在线官网| 久久99热这里只频精品6学生| 国产在线一区二区三区精| 国产精品一区二区在线观看99| 最新中文字幕久久久久| 国产精品99久久久久久久久| 80岁老熟妇乱子伦牲交| 免费av不卡在线播放| 久久精品国产自在天天线| 亚洲av电影在线观看一区二区三区 | 我要看日韩黄色一级片| 亚洲精品国产av成人精品| 欧美变态另类bdsm刘玥| 26uuu在线亚洲综合色| 男人和女人高潮做爰伦理| 国产精品不卡视频一区二区| 日韩制服骚丝袜av| 欧美xxxx性猛交bbbb| 伦理电影大哥的女人| 国产欧美另类精品又又久久亚洲欧美| 国产亚洲一区二区精品| 熟女电影av网| av在线天堂中文字幕| 精品亚洲乱码少妇综合久久| 日韩av不卡免费在线播放| 综合色丁香网| 99热这里只有是精品50| 97人妻精品一区二区三区麻豆| 免费观看的影片在线观看| 亚洲人成网站在线观看播放| 一区二区av电影网| 亚洲欧美中文字幕日韩二区| 九草在线视频观看| 97在线视频观看| 亚洲久久久久久中文字幕| 国产精品一区二区性色av| 久久人人爽av亚洲精品天堂 | 国产熟女欧美一区二区| 久久人人爽av亚洲精品天堂 | 亚洲av日韩在线播放| 99九九线精品视频在线观看视频| 中文在线观看免费www的网站| 亚洲精品乱码久久久久久按摩| 欧美zozozo另类| 男人添女人高潮全过程视频| 成人国产麻豆网| 精品人妻熟女av久视频| 久久久久久久久久人人人人人人| 97在线视频观看| 久久精品国产亚洲av天美| 亚洲精品亚洲一区二区| 亚洲精品国产av成人精品| 尾随美女入室| 国产精品国产三级国产av玫瑰| 日韩一区二区三区影片| 国产熟女欧美一区二区| 日本爱情动作片www.在线观看| 日韩视频在线欧美| 少妇人妻 视频| 精品人妻熟女av久视频| 免费看光身美女| xxx大片免费视频| 嘟嘟电影网在线观看| 国内精品美女久久久久久| 欧美日韩视频高清一区二区三区二| 伊人久久国产一区二区| 2022亚洲国产成人精品| 高清欧美精品videossex| 99re6热这里在线精品视频| 国产精品国产三级国产av玫瑰| 亚洲成人一二三区av| 麻豆久久精品国产亚洲av| 免费观看无遮挡的男女| 中文资源天堂在线| 最近最新中文字幕免费大全7| 啦啦啦啦在线视频资源| 亚洲av欧美aⅴ国产| 性插视频无遮挡在线免费观看| 国产成人午夜福利电影在线观看| 国产成年人精品一区二区| 欧美日韩亚洲高清精品| freevideosex欧美| 精品一区二区三区视频在线| 岛国毛片在线播放| 噜噜噜噜噜久久久久久91| 国产乱人偷精品视频| 黄片无遮挡物在线观看| 男人爽女人下面视频在线观看| 久久鲁丝午夜福利片| 久久精品综合一区二区三区| 男女那种视频在线观看| 国产伦精品一区二区三区视频9| 在线免费观看不下载黄p国产| 精品国产一区二区三区久久久樱花 | 97超视频在线观看视频| 麻豆精品久久久久久蜜桃| 日本免费在线观看一区| 欧美激情在线99| 亚洲精华国产精华液的使用体验| 久久久国产一区二区| 成年av动漫网址| 午夜福利高清视频| 成人特级av手机在线观看| 卡戴珊不雅视频在线播放| 亚洲国产精品专区欧美| 五月玫瑰六月丁香| 久久精品国产亚洲网站| 日韩一区二区视频免费看| 中文天堂在线官网| 欧美另类一区| 欧美xxxx性猛交bbbb| 麻豆成人av视频| 又大又黄又爽视频免费| 色播亚洲综合网| 直男gayav资源| 三级国产精品欧美在线观看| 少妇被粗大猛烈的视频| 校园人妻丝袜中文字幕| 超碰97精品在线观看| 日本一本二区三区精品| 亚洲精品自拍成人| 校园人妻丝袜中文字幕| 直男gayav资源| 交换朋友夫妻互换小说| 久久综合国产亚洲精品| 久久精品国产亚洲av天美| 国产熟女欧美一区二区| 91狼人影院| 国产老妇女一区| 久久99热这里只有精品18| 国模一区二区三区四区视频| 欧美97在线视频| 国产精品偷伦视频观看了| 69av精品久久久久久| 极品少妇高潮喷水抽搐| 亚洲精品aⅴ在线观看| 国产综合精华液| 亚洲图色成人| 永久网站在线| 99九九线精品视频在线观看视频| 午夜精品国产一区二区电影 | 国产免费一区二区三区四区乱码| 亚洲人成网站高清观看| av播播在线观看一区| 下体分泌物呈黄色| 色网站视频免费| 777米奇影视久久| 精品久久久久久久久av| 大片电影免费在线观看免费| 亚洲综合精品二区| 最近最新中文字幕大全电影3| 国产乱人视频| 一区二区三区四区激情视频| 男女那种视频在线观看| 神马国产精品三级电影在线观看| 亚洲,一卡二卡三卡| 男女那种视频在线观看| 2022亚洲国产成人精品| av卡一久久| 天堂网av新在线| 国产亚洲5aaaaa淫片| 男男h啪啪无遮挡| 黄色欧美视频在线观看| 亚洲精品,欧美精品| videossex国产| 男女无遮挡免费网站观看| 日韩大片免费观看网站| 观看免费一级毛片| 免费av毛片视频| 午夜精品国产一区二区电影 | 欧美zozozo另类| 精品少妇黑人巨大在线播放| 超碰av人人做人人爽久久| 小蜜桃在线观看免费完整版高清| 国产一区二区亚洲精品在线观看| 成人一区二区视频在线观看| 亚洲欧洲日产国产| 联通29元200g的流量卡| 大片电影免费在线观看免费| 可以在线观看毛片的网站| 在线观看三级黄色| av黄色大香蕉| 身体一侧抽搐| 欧美一区二区亚洲| 亚洲图色成人| av免费观看日本| 超碰97精品在线观看| 一区二区三区精品91| 免费黄色在线免费观看| av在线观看视频网站免费| 久久人人爽人人片av| 亚洲精品自拍成人| 久久99热6这里只有精品| 91午夜精品亚洲一区二区三区| 五月天丁香电影| 成年女人看的毛片在线观看| 亚洲精品国产av成人精品| 人体艺术视频欧美日本| 国内精品宾馆在线| 欧美精品人与动牲交sv欧美| 少妇 在线观看| 午夜日本视频在线| 麻豆乱淫一区二区| 深夜a级毛片| 狂野欧美白嫩少妇大欣赏| 亚洲色图av天堂| 噜噜噜噜噜久久久久久91| 欧美成人午夜免费资源| 国产高潮美女av| 丰满乱子伦码专区| 成年av动漫网址| 日本与韩国留学比较| 国产片特级美女逼逼视频| 1000部很黄的大片| 最近2019中文字幕mv第一页| 久久99精品国语久久久| 女人被狂操c到高潮| 国产精品偷伦视频观看了| 成年免费大片在线观看| 亚洲熟女精品中文字幕| 欧美高清性xxxxhd video| 国产精品秋霞免费鲁丝片| 成年av动漫网址| 久久精品国产自在天天线| 中文字幕免费在线视频6| 亚洲av成人精品一二三区| 国产在线一区二区三区精| 亚洲国产av新网站| 亚洲人成网站高清观看| 国产亚洲午夜精品一区二区久久 | av国产精品久久久久影院| 免费大片黄手机在线观看| 国产免费视频播放在线视频| 高清午夜精品一区二区三区| 亚洲国产色片| 亚洲真实伦在线观看| 亚洲欧美日韩另类电影网站 | 亚洲av中文字字幕乱码综合| 亚洲美女视频黄频| 国产视频内射| av在线蜜桃| 性插视频无遮挡在线免费观看| 熟女av电影| av又黄又爽大尺度在线免费看| 香蕉精品网在线| 日韩欧美 国产精品| 能在线免费看毛片的网站| 一个人看视频在线观看www免费| 色哟哟·www| 精品亚洲乱码少妇综合久久| 国产精品成人在线| 久久久久精品性色| 欧美日韩一区二区视频在线观看视频在线 | 亚洲欧美日韩东京热| www.色视频.com| 亚洲精品第二区| 制服丝袜香蕉在线| 嫩草影院入口| 久久久久久久久久成人| av一本久久久久| 国产亚洲91精品色在线| 久久精品综合一区二区三区| 免费黄频网站在线观看国产| 日韩欧美精品v在线| 丝袜脚勾引网站| 九九在线视频观看精品| 国产黄片视频在线免费观看| 免费av不卡在线播放| 免费观看的影片在线观看| 中文天堂在线官网| 干丝袜人妻中文字幕| 日本午夜av视频| 久久久久久久国产电影| 欧美日韩视频精品一区| 久久精品国产亚洲网站| 在线精品无人区一区二区三 | 亚洲天堂国产精品一区在线| 男女边摸边吃奶| 亚洲四区av| 三级男女做爰猛烈吃奶摸视频| 亚洲精品久久午夜乱码| 91精品一卡2卡3卡4卡| 久久久久久久午夜电影| 男人添女人高潮全过程视频| 一本—道久久a久久精品蜜桃钙片 精品乱码久久久久久99久播 | 在线看a的网站| 亚洲精品第二区| 哪个播放器可以免费观看大片| 99热这里只有精品一区| 我的女老师完整版在线观看| 精品久久久久久久人妻蜜臀av| 亚洲av一区综合| 五月开心婷婷网| 美女cb高潮喷水在线观看| 久久久久久九九精品二区国产| 性色av一级| 亚洲美女视频黄频| 夜夜爽夜夜爽视频| 超碰av人人做人人爽久久| 亚洲色图av天堂| 国内精品宾馆在线| 国产黄频视频在线观看| 国产乱人偷精品视频| 2018国产大陆天天弄谢| 十八禁网站网址无遮挡 | 搡老乐熟女国产| 天天一区二区日本电影三级| 美女cb高潮喷水在线观看| 亚洲精品视频女| 久久久久久久久久久丰满| 国产精品一及| 人妻 亚洲 视频| 在线看a的网站| 亚洲精品日本国产第一区| 亚洲国产日韩一区二区| 欧美一区二区亚洲| 看黄色毛片网站| 色哟哟·www| 夜夜爽夜夜爽视频| 成人国产av品久久久| 少妇人妻久久综合中文| 久热这里只有精品99| 青青草视频在线视频观看| 成年版毛片免费区| 下体分泌物呈黄色| 五月伊人婷婷丁香| 又黄又爽又刺激的免费视频.| 中文精品一卡2卡3卡4更新| 久久久欧美国产精品| 国产免费视频播放在线视频| 永久免费av网站大全| 青春草国产在线视频| 精品久久久久久电影网| 赤兔流量卡办理| 最近中文字幕高清免费大全6| 成人亚洲欧美一区二区av| 国产成人aa在线观看| 日韩中字成人| 男女无遮挡免费网站观看| 日韩成人av中文字幕在线观看| 国产熟女欧美一区二区| 久久久久性生活片| 免费av毛片视频| 亚洲人成网站在线观看播放| 国产精品久久久久久av不卡| 22中文网久久字幕| 亚洲久久久久久中文字幕| 久久久国产一区二区| 伊人久久精品亚洲午夜| 国产精品熟女久久久久浪| 国产午夜精品久久久久久一区二区三区| 尤物成人国产欧美一区二区三区| 夫妻午夜视频| 国产毛片a区久久久久| 国产老妇伦熟女老妇高清| 成人黄色视频免费在线看| videos熟女内射| 搡老乐熟女国产| 久久久久久久久久人人人人人人| 亚洲精品成人av观看孕妇| 99久久精品热视频| 国产伦精品一区二区三区四那| 女人十人毛片免费观看3o分钟| av在线蜜桃| 国产v大片淫在线免费观看| 777米奇影视久久| 国产成人一区二区在线| 国产亚洲精品久久久com| 亚洲自拍偷在线| 久久久久久久亚洲中文字幕| 欧美一区二区亚洲| 日本-黄色视频高清免费观看| 免费看av在线观看网站| 黑人高潮一二区| 国产探花极品一区二区| 国产老妇伦熟女老妇高清| 欧美日韩精品成人综合77777| 亚洲三级黄色毛片| 亚洲aⅴ乱码一区二区在线播放| 成人毛片a级毛片在线播放| 国产91av在线免费观看| 国产成人福利小说| 亚洲av二区三区四区| 亚洲成人久久爱视频| 草草在线视频免费看| 国产精品人妻久久久久久| 午夜老司机福利剧场| 亚洲精品,欧美精品| 又黄又爽又刺激的免费视频.| 观看美女的网站| 欧美三级亚洲精品| 热re99久久精品国产66热6| 大码成人一级视频| 岛国毛片在线播放| 亚洲人与动物交配视频| 日韩制服骚丝袜av| 18+在线观看网站| 国模一区二区三区四区视频| 国产精品福利在线免费观看| 天堂网av新在线| 久久6这里有精品| 久久热精品热| 极品少妇高潮喷水抽搐| 一级毛片久久久久久久久女| 三级国产精品欧美在线观看| 亚洲人成网站在线播| 精品一区在线观看国产| 高清午夜精品一区二区三区| 国产精品伦人一区二区| 精品久久久久久久久av| 夫妻性生交免费视频一级片| 久久精品国产亚洲网站| av又黄又爽大尺度在线免费看| 黄色配什么色好看| 少妇人妻精品综合一区二区| 少妇人妻久久综合中文| 99re6热这里在线精品视频| 麻豆成人午夜福利视频| 久久99热这里只有精品18| av在线亚洲专区| 中文乱码字字幕精品一区二区三区| 免费播放大片免费观看视频在线观看| 国产成人午夜福利电影在线观看| 日本一二三区视频观看| 精品久久久久久久久亚洲| 波多野结衣巨乳人妻| 国语对白做爰xxxⅹ性视频网站| 欧美日韩视频精品一区| 久久久精品免费免费高清| 性色avwww在线观看| 亚洲精品日韩av片在线观看| 免费人成在线观看视频色| 91精品伊人久久大香线蕉| 亚洲真实伦在线观看| 中文字幕久久专区| 欧美性感艳星| 三级经典国产精品| 男人和女人高潮做爰伦理| 网址你懂的国产日韩在线| 精华霜和精华液先用哪个| 中文字幕免费在线视频6| 99久久精品热视频| 极品少妇高潮喷水抽搐| 国产男女内射视频| 中国美白少妇内射xxxbb| 免费观看性生交大片5| 三级经典国产精品| 男人和女人高潮做爰伦理| 日产精品乱码卡一卡2卡三| 内地一区二区视频在线| 国产亚洲91精品色在线| 国产毛片在线视频| 亚洲婷婷狠狠爱综合网| 成人毛片a级毛片在线播放|