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

    Dynam ic Pre-ultim ate Strength Evaluation of Containership based on a 2D M odified Hydroelasticity M ethod in Extrem e W aves

    2016-05-15 13:24:13LIUWeiqinPEIZhiyongWUWeiguoSUZUKIKatsuyuki
    船舶力學(xué) 2016年9期
    關(guān)鍵詞:武漢理工大學(xué)船體彎矩

    LIU Wei-qin,PEI Zhi-yong,WU Wei-guo,SUZUKI Katsuyuki

    (1.Departments of Navel Architecture,Ocean and Structural Engineering,School of Transportation,Wuhan University of Technology,Wuhan 430063,China;2.Department of Systems Innovation,Graduate School of Engineering,the University of Tokyo,Tokyo 113-8656,Japan)

    Dynam ic Pre-ultim ate Strength Evaluation of Containership based on a 2D M odified Hydroelasticity M ethod in Extrem e W aves

    LIU Wei-qin1,PEI Zhi-yong1,WU Wei-guo1,SUZUKI Katsuyuki2

    (1.Departments of Navel Architecture,Ocean and Structural Engineering,School of Transportation,Wuhan University of Technology,Wuhan 430063,China;2.Department of Systems Innovation,Graduate School of Engineering,the University of Tokyo,Tokyo 113-8656,Japan)

    Extreme waves have led to many accidents and losses of ships at sea.In this paper,a 2D modified hydroelasticity method is proposed as a means of studying the dynamic pre-ultimate strength of a containership instead of post-ultimate strength when traversing extreme waves,while considering the ultimate strength of the ship.Traditional ultimate strength evaluations are undertaken by a quasi-static assumption and the dynamic wave effect is not considered.And,the dynamic response of a ship as induced by a wave is studied on the basis of the hydroelasticity theory so that the nonlinear structural response of the ship cannot be obtained for large waves.So,a 2D modified hydroelasticity method,which takes the coupling between time-domain waves and the nonlinear ship beam into account,is proposed.This method is based on a hydroelasticity method and a Smith method that combine the wave load and the structural nonlinearity.A dynamic reduction in rigidity related to deformation could influence the strength and curvature of a ship’s beam;So,it is input into a dynamic hydrodynamic formula rather than being regarded as a constant structural rigidity in a hydroelastic equation.A number of numerical extreme wave models are selected for computing the modified hydroelasticity,such that large deformations occur and nonlinear dynamic VBM is generated when the ship traverses these extreme waves.Four important extreme wave patterns such as maximum wave height, regular wave height,wave speed,and wave length ratio are changed to study their influence on modified hydroelastic results including nonlinear dynamic VBM and deformational angle at midship,and the modified hydroelastic results are computed and compared with results obtained with the hydroelasticity method,some differences are observed and conclusions are drawn.

    dynamic pre-ultimate strength;2D hydroelasticity;Smith method; extreme wave;nonlinear VBM;angle

    0 Introduction

    The conventional tendency in ship engineering is to evaluate structures based on theirultimate strength,requiring that the ultimate strength of hull girders be evaluated according to the Common Structural Rules(CSR)for Tankers and Bulk Carriers,as specified by the International Association of Classification Societies,Ltd.,(IACS,2006).Many researchers have been pursuing a means of developing an efficient and accurate method of evaluating the ultimate strength of hull girders.Most approaches to computing ultimate strength have been based on static methods.Once ultimate strength of ship is reached,most accidents involving ships have resulted from with a dynamic extreme wave.Pre-ultimate strength denotes the strength still unreached the ultimate strength point while post-ultimate strength means that the strength beyond the ultimate strength point.It is general practice to study pre-ultimate strength of ship for the society of naval architecture engineering and its related rules.So this paper just calculates dynamic pre-ultimate strength of a ship instead of post-ultimate strength in extreme waves or other giant waves by using a 2D modified hydroelastic method.

    ‘Extreme waves’are more precisely defined as waves with a height more than twice the significant wave height.Furthermore,an extreme wave differs from regular head waves in that it is a single wave.Peliniovsky et al[1]cataloged major marine accidents around the world,and found that at least 22 super-carriers had been lost because of encounters to extreme waves between 1969 and 1994,causing 525 fatalities.He performed several computations of wave loads encountered in the Pacific and the Atlantic.Haver et al[2]recorded an extreme wave at the Draupner platform in the North Sea,providing indisputable evidence that such waves exist.While the significant wave height was 11.92 m,the maximum wave height reached 25.63 m. To understand the physical mechanisms related to such extreme waves,Waseda[3]and Minami[4]have been studying how such waves are formed by simulating extreme waves in a laboratory tank,waves were generated while referring to wave simulations performed with a nonlinear numerical wave tank program,NWT2D.

    Traditional hydroelastic theory cannot model the nonlinear structural response of a ship to extreme waves,and the most popular static methods for calculating the nonlinear strength of a ship do not address the dynamic wave problem.Yamamoto et al[5]proposed a time-domain 2D hydroelastic method to calculate the motions and longitudinal strength of a ship as it traverses in waves,while taking the specially nonlinear slamming force was into consideration. However,this method cannot be used to compute the dynamic nonlinear structural response of a ship when faced with in an extreme wave.Masaoka et al[6]proposed a numerical approach to analyze the collapse of a ship’s hull girder when traversing an extreme wave,which was based on a structural analysis system that they had developed.The dynamic FEM system can consider a ship’s longitudinal collapse behavior,whereby the nonlinear fluid forces acting on a ship’s hull were calculated using the Ursell-Tasai method.This paper considers fluid forces acting on ship’s hull by using finite boundary method.Iijima et al[7]proposed a modified hydroelasticity method for a model test,whereby two rigid segments are connected by plastic hinges of different rigidities,rotational angle at a midship section.Different from scale model of above modified hydroelasticity method of Iijima,this paper uses an actual containership to consider the ship collapse due to an extreme wave.Interaction and coupling between in ship structureand waves is also addressed.Xu et al[8]conducted a number of experimental investigations into the collapse behavior of a box-shape hull girder subjected to extreme wave-induced loads, some sacrificial specimens with circular pillar and trough shapes which respectively show different bending moment-displacement characteristics were mounted to compare the dynamic collapse characteristics of the hull girder in waves.And then the post-ultimate strength behavior of the hull girder including the global deformation and motions under extreme wave induced loads are studied by Xu et al[9].Masaoka and Okada[10]used a FEM beam analyzing system to investigate longitudinal collapse behavior by considering the interaction between structural and fluid non-linearity,the collapse behavior of the cross sections of the hull girders is expressed as non-linear bending moment-curvature relationship,an important research direction is guided to modified hydroelasticity.Pei et al[11]developed a total simulation system combining load calculation and structural collapse analysis applied to simulate progressive collapse behaviour of a single-hull Kamsarmax type bulk carrier,the applicability to structure system,high accuracy and sufficient efficiency of ISUM had been demonstrated,but coupling between ship structure and waves are not considered.Liu et al[12]used a strip code to calculate the wave loads of a containership when it transverses extreme waves,and then,these wave loads were applied to a nonlinear dynamic FEM model to obtain the nonlinear dynamic VBM and angles at midship.Liu et al[13]proposed a modified hydroelastic method to evaluate ship nonlinear dynamic strength of ship beam in extreme waves,regarding ship structure as a pure elasto-plastic beam with sulciform section,no buckling effect of ship girder is considered while Liu et al[14]considers the interaction and coupling between in ship structure and waves,a simplified progressive collapse method taking the buckling effect of ship structure into account is combined into the 2D modified hydroelasticity method.

    Although FEM was used to calculate nonlinear dynamic strength of ship induced by large waves,we have to face some problems due to FEM.Computational efficiency is too slow and computational cost is too large for shell/solid FE model when ship’s nonlinearity is considered,so the FEM beam model is only choice as Ref.[10]for the modified hydroelastic problem.FEM needs a fixed boundary condition to calculate strength of ship so that the large ship motions induced by extreme waves are not obtained,even dynamic problems to balanced loads, the boundary conditions are added only to suppress the rigid body motion(fix only 6DOFs), which does not lose generality,present commercial codes do not supply a function to account for the wave-structure interaction.This method that is capable of considering 2D modified hydroelasticity was developed by combining the hydroelasticity and Smith method in this paper,the computational efficiency is improved,the interaction and coupling between the ship structure and large waves were studied using a simplified progressive collapse method that accounts for the buckling effect of ship structure,and a FEM beam module is used to calculate the vary ship beam modal shape by varying rigidity of ship beam in this method.At last, extensive extreme wave patterns are analyzed for their influence on nonlinear dynamic strength and large deformation.

    1 Theoretical background

    1.1 Overview

    Hydroelasticity method has been solving elastic response of ship structure in waves.However,this method can not compute dynamic nonlinear structural response of ship in extreme waves.In addition,simplified progressive collapse method(Smith method) has always been simple but efficient method to analyze the progressive collapse behavior of ship girder under longitudinal bending,static ultimate strength can be obtained,but the static method can not analyze the nonlinear dynamic collapse behavior of ship girder in time-domain extreme waves.In this paper,Smith method and 2D hydroelasticity method are integrated to develop a 2D modified hydroelasticity method,this method considers plasticity and buckling nonlinearities of ship section,rigidity of ship beam which is dependent on deformational curvature is alterable rather than constant value in hydroelastic theory,so nonlinear dynamic strength and large deformation can be obtained.Fig.1 shows methodology of modified hydroelasticity method.

    1.2 2D Hydroelasticity theory

    Fig.1 Methodology of hydroelasticity-plasticity method

    A ship can be modeled as a non-uniform free-free beam,and is often assumed to be a Bernoulli-Euler beam.A right-handed coordinate system O( x0,y0,z0)is fixed in space.The(x0,y0)plane corresponds to the surface of the water,and Z0is directed downwards.Another right-handed coordinate system O x,y,()z moves forward at a constant speed U,equal to the speed of the ship,and the origin o coincides with a point on the surface of the water,dropped perpendicularly from the aft of the ship.Here,x corresponds to the forward motion of the ship and z is vertically downwards.At the initial time t=0,the two coordinate systems coincide.In terms of the incident waves,it was first considered that a regular wave with an amplitude of ζ0propagates in the negative x direction.Fig.2 shows the coordinate system of the 2D modified hydroelasticity method.For an infinitesimal segment,d x,we can apply the following equilibrium formula.

    Fig.2 Coordinate system of the 2D modified hydroelasticity method

    where M=vertical bending moment;F=shear force;fe=external force;m=mass per unit length of ship;and w=vertical deflection including elastic-plastic and rigid deflection.Moreover,the vertical bending moment M can be expressed as

    where η=equivalent structural damping coefficient;E=Young’s modulus;and I=section moment of inertia.The product term EI is called the‘rigidity of the ship’.Substituting Eq.(2)into Eq. (1),gives Eq.(3),as follows:

    Eq.(3)defines the hydroelastic formula equilibrium.Here,the external force feincludes the wave damping force fr,the weight of the ship fg,the restoring force fs,and the hydrodynamic force fm′.

    where the hydrodynamic force fm′takes the hydrodynamic effect of a wave on the ship into account.Even extreme wave,this force remains dominant.Yamamoto et al[5]derived the hydrodynamic force equation presented in Eq.(5).This takes the slamming effect and wave nonlinearity into account.

    where MHis added mass derived by the boundary element method,Eq.(5)includes the slamming force component,fslam.

    The second part of the hydrodynamic force is the damping force,which is proportional to the speed V:

    where N( x,t)is wave-making damping coefficient per strip which is different by section changing,it is also obtained by using the boundary element method.

    The restoring force is given by the following Equation:

    The inertial force is given by the following equation:

    The force of the ship’s weight due to gravity is expressed as follows:

    In addition,the modified hydroelasticity method combines the progressive collapse method with hydroelastic theory.The quasi-static method is used to determine the structural section bending moment and rigidity.The rigidity is determined from the following equation:

    where EIpis the nonlinear rigidity,Mpis the nonlinear elastic-plastic bending moment,and w″is the curvature of the section.The relationship between the rigidity and curvature is determined using the progressive collapse method.Moreover,as the rigidity changes,the natural modal shape can be calculated to obtain the renewed vertical deflection wp.The vertical deflection wpshould be satisfied with a boundary condition of free-free.This is expanded by superposing a series of the structural natural modal shape.

    where Wj(x)which is structural natural dry modal shape.So the deflection up can be expressed from the mode shapes of a free-free beam.Eq.(12)corresponds to modal superposition,which is not suitable for a nonlinear 3D model.Because buckling may occur with a 3D ship model as it traverses an extreme wave,the deformed shape will be abnormal as long as buckling occurs,making modal superposition inapplicable to the 3D model.However,as this study considers a beam model,it is assumed that the beam model does not incur any buckling,and its deformed shape is normal.

    The angle deformation can be addressed by using differential equation of global deflection wp.

    The beam rigidity is reduced to consider its capacity to resist deformation under the influence of an extreme wave,so that modal superposition can be applied to a beam model with 2D modified hydroelasticity.The dynamic rigidity,EIpwhich is related to the curvature and modified modal shape,is substituted into Eq.(3).

    Eq.(12)is substituted into Eq.(3)to give a dynamic modified hydroelastic equilibrium formula,as given by Eq.(14).This gives the plastic effect,and the added mass and wave-damping coefficient are obtained by using the 2D finite boundary element method.

    Eq.(15)involves the term of inertia force matrix[M],the term of damping force matrix[C],the term of stiffness force matrix[K]and the term of external force matrix[F],all matrix of Eq.(15)are expressed as from Eq.(16)to Eq.(19).

    The elements of the matrix are as above,the subscript of i or j denotes the number of natural modals,with the first two modals being the rigid ship motions,and the subsequent modals being the flexibility modals.Once the initial conditions are set,the increments can be obtained for a time series by using a numerical integration scheme such as the Newmark-β method.

    As long as the acceleration term satisfies the iteration condition,the calculation continues to next time step.ε is iterative permissible error.

    In this paper,the modified hydroelasticity method is verified for its convergence to determine the number of the integral points.Convergence is verified by using the relative acceleration convergence rate as Eq.(25)whenis just satisfied the iteration condition as Eq. (24)at the iteration step of n+1.

    2 Containership model and com putation of ultimate strength

    2.1 Containership model

    The final objective of this study is to compute nonlinear dynamic strength for a containership in a simulated extreme ocean wave. For this reason,a 500TEU containership is introduced to access its modified hydroelastic computation and wave load.The main principal dimensions of 500TEU container ship are summarized in Tab.1.And the body plan of 500TEU containership is presented in Fig.3.

    In this paper,a module for calculating modal shape of FEM beam is used.The wave lengths of all extreme waves would be defined as being equal to the ship length because VBM is the largest so that the condition is the most dangerous.The deformational shape of the ship beam is highly dependent on the wave length.Even if the first flexuous natural mode was assumed to be dominant when the wave length was set such that it was equal to the ship length,former ten flexuous natural modes are applied to superpose the deformational shape of ship in the extreme waves in the study.Therefore,it is sufficient to present the deformational shape of ship as the former ten flexuous natural modes.Fig.4 shows ten flexuous modal shapes of ship and their natural frequencies.

    Fig.3 Body plan of 500TEU container ship

    Tab.1 Principal dimensions of ship model

    Continue Tab.1

    Fig.4 Former ten flexuous modal shapes of ship

    Fig.5 Ship weight distribution

    This computation is carried out based on an assumption that ship is under design full loading whose draft line is designed loading line.Fig.5 shows the weight distribution of 500-TEU containership.Generally,maximum weight is distributed around the middle part of ship due to the loading of cargo while less weight is appeared at two ends of ship.For the same reason,weight distribution requires the corresponding strength and rigidity distribution.In fact,no matter rigidity or the moment-curvature relation is different along ship length,maximum rigidi-ty appears at midship.Fig.6 presents the elastic rigidity(EI)distribution of 500-TEU containership in elastic stage.

    Fig.6 Ship elastic rigidity distribution

    2.2 Ultimate strength of 500TEU containership

    For reasons of high efficiency and low computational cost,the progressive collapse method has been applied to the computation of the ultimate strength of a ship’s girder.In this study, structural nonlinearity was considered based on the computation of the progressive collapse method,with a flat structural section of a ship being divided into different stiffened elements according to their collapse failure modes,plasticity,and buckling.

    The code for the simplified progressive collapse method was programmed to calculate the ultimate strength and determine the relationship between the bending-curvature and rigiditycurvature for a 500TEU containership.Ideal elastic-plastic steel material was used,with a yield stress of 2.35e+008 Pa.Plasticity was taken into account,and the buckling collapse behavior of the stiffened plate elements was considered under compression.In this calculation,the relations of bending-curvature and rigidity-curvature are different along the ship length.Each relation of rigidity-curvature dominates the strength at each strip though it is a great of possible that ship collapse is appeared at midship.Fig.7 shows a structural section of a midship sec-tion of a 500TEU containership and the stress distribution of the structural elements under ultimate sagging strength.The Smith method is used to calculate the ultimate strength,and the three typical stiffened elements(A,B,and C)shown in Fig.7 are assumed to buckle under pressure.Element A represents plate element,element B represents the ordinary stiffened element with angle steel,and element C represents the stiffened element with a T section bar.

    Fig.7 Stress distribution of structural section

    Fig.8 Relationship between average strain and average stress of three stiffened elements

    The Fig.8 shows the relationship between the average strain and average stress of the three stiffened elements.In Fig.8,relative stress σ/σYmeans stress divided by yield stress and relative strain ε/εYdenotes strain divided by yield strain.This relationship is as stated in the Common Structural Rules for Tankers and Bulk Carriers.

    The ultimate bending moments are presented as follows:

    where Ms is the ultimate sagging bending moment,which occurs when the curvature is equal to 6.842 0e-04(1/m).Mh is the ultimate hogging bending moment.The curvature for each bending moment is plotted in Fig.9.It can be seen that the ultimate sagging bending moment is smaller than the ultimate hogging bending moment.It is known that the deck is further from the neutral axis than the keel of the ship,such that buckling of the deck occurs when the ship is experiencing sagging bending,but bottom buckling cannot occur easily in the hogging bending case.

    Fig.9 Bending moment-curvature of 500TEU containership at midship

    Fig.10 Rigidity-curvature of 500TEU containership at midship

    The sectional rigidity is calculated through the application of Eq.(7).The sectional rigidity in the elastic stage is EI=2.26e+012(Nm2),falls in the elasto-plastic stage,becomes infinitely close to zero at the point corresponding to the ultimate bending moment,and then becomes negative.In the case of sagging bending,when the curvature of a section exceeds 2. 013 0e-04,that section of the ship’s beam enters the elasto-plastic stage.Fig.10 shows the rigidity-curvature curve,from which it can be determined that the rigidity remains constant while the ship is in the elastic stage but falls rapidly when the ship enters the elasto-plasticstage,ultimately falling to zero at the point corresponding to the ultimate bending moment.It is known that the load falls beyond the ultimate strength point while the sectional rigidity becomes negative.This negative rigidity results in the reduction of ship strength further.However,this paper just studies the dynamic pre-ultimate strength,the post-ultimate strength is not studied in this paper.So the negative rigidity of ship is not integrated in the modified hydroelastic code.

    2.3 Extreme wave models

    A numerical wave tank code NWT2D(New Wave Tanker 2D)performed in some previous studies is used to compose numerical extreme wave models.The program NWT2D was developed by National Maritime Research Institute in Japan.Extreme wave generation is performed by wave superposition of linear numerical regular wave and numerical focusing wave,focusing wave is generated in specific position and time,focusing wave amplitude and energy superposition are coincided with regular wave amplitude and energy to compose extreme wave which has a higher wave height.Fig.11 shows the numerical extreme wave by superposition of numerical regular wave and focusing wave,(a)is numerical regular wave elevation,(b)is numerical focusing wave elevation,(a)is superposition with(b)to generate(c),(c)is numerical extreme wave.In this study,numerical extreme wave model is also linear wave.

    Fig.11 Numerical extreme wave by superposition between numerical regular wave and focusing wave

    Fig.12 Wave elevation profile of containership at midship

    Three extreme wave models(Case 1,Case 2 and Case 3)are presented in Tab.2 as typical extreme wave model,their regular waves are kept similarly to make sure computational comparability,regular wave height is 2.8 m while the ratio of wave length to ship length(Lw/ L)is equal to 1.0.Extreme wave height is paid a close attention,three different focusing wave heights are used to change extreme wave heights such as 6.0,11.4 and 11.6.The ratios ex-treme wave height divided by D(molded depth)are 0.7,1.33 and 1.35.These three extreme wave Cases are defined to indicate three typical strength condition of elastic strength,elastoplastic strength and ultimate strength.Time-domain profiles of three extreme wave models are presented in Fig.12.

    Tab.2 Typical extreme wave model

    3 Calculation of modified hydroelasticity

    In this paper,a method for calculating modified hydroelasticity is proposed to enable the study of the nonlinear dynamic response of a containership to the extreme waves described above.Firstly,the convergence of modified hydroelasticity is discussed by remeshing ship beam to different number of elements along ship length.And then a number of computations including deformational angle and VBM at midship were performed using the TPSLAM modified hy-droelastic code,they are plotted and compared with the calculations of hydroelasticity,so that the difference between the nonlinear and linear dynamic results can be obtained transparently.In addition,the course of any collapse and the reasons for the collapse are studied for extreme waves.

    Fig.13 Midship angle obtained with modified hydroelasticity and hydroelasticity

    A parameter of angle derived from Eq.(9)is used to manifest itself as deformation of the ship’s beam,it is solved from differential equation of global deflection wp by length.The angle at the midship section warrants close attention.The low rigidity of the beam generates a large deformational angle.Fig.13 shows the angle curves obtained with the modified hydroelasticity and hydroelasticity methods at the midship section when subjected to the three extreme waves described above.Fig.13,(a)shows that there is no difference between the two methods,(b)indicates that the modified hydroelasticity method generates a larger angle deformation than the hydroelasticity method,with a residual deformational angle(rs2=0.075°)being generated after the first large peak.Next,(c)shows a much greater deformation,with the maximum angle obtained with the modified hydroelasticity method being 1.103 6 degrees with a residual angle(rs3=0.665 9°),while the maximum angle obtained with the hydroelasticity method alone has a very small value of 0.345 degrees.

    Fig.14 Sagging bending moment and elasto-plastic unloading

    It is observed from Fig.13(b)and(c)that residual deformational angles are resulted in. The modified hydroelasticity method considers residual deformation due to elasto-plastic unloading.As to above three extreme waves,sagging bending results in large deformation and residual deformation.As long as ship enters elasto-plasticity,unrecoverable residual deformation is induced.Fig.14 shows the elasto-plastic unloading routes of Case 2 and Case 3.Residual angle rs2 induced in Case 2 and residual angle rs3 induced in Case 3 are generated due toelasto-plastic unloading.Dynamic modified hydroelastic analysis shown in Fig.13 reveals the residual angles rs2 and rs3 after the top of dynamic angle curves.

    Fig.15 VBM as obtained with modified hydroelasticity and hydroelasticity

    In this study,the nonlinear vertical bending moment(VBM)is obtained from the calculation of the relationship between the curvature and bending moment,as shown in Fig.9.In previous research,VBM was calculated by using linear hydroelastic code since the structural rigidity is assumed to be constant when using the hydroelasticity method.In this study,VBM was computed by using a modified hydroelasticity method,with the structural nonlinear rigidity of the ship’s beam being regarded as a dynamic time-domain,as determined from the deformational curvature.A large angle occurs when an extreme wave impacts a ship’s beam as shown in Fig.13(c),is calculated by using the modified hydroelasticity method.Therefore,a larger VBM peak is obtained when using the modified hydroelasticity method,in comparison with the hydroelasticity method.Fig.15 presents the VBM curves as obtained with the modified hydroelasticity and hydroelasticity methods at midship points with the above three extreme waves.Fig.15(a)shows the VBM for the Case of CASE 1 at S.S.5.0.It can be seen that results obtained with the hydroelastic-plastic and hydroelastic methods are in very good agreement,meaning that the CASE 1 extreme wave generates only an elastic response from the ship. Next,(b)shows that there is very little difference at the first large peak,but that the VBM as determined using the modified hydroelasticity method is larger than that obtained with the hydroelasticity method,meaning that the extreme Case of CASE 2 produces an elastic-plastic response.Finally,(c)shows a considerable difference for the first large peak,with the maxi-mum VBM being attained at the ultimate sagging bending moment Ms.This implies that in the extreme wave Case of CASE 3,the response is fatal.Both(b)and(c)indicate that the modified hydroelasticity method produces a larger VBM result than the hydroelasticity method when ship beam enters the stage of large deformation.

    Fig.16 Hydrodynamic VBM and hydrostatic VBM of CASE 3

    Fig.17 Models obtained with hydroelasticity and modified hydroelasticity models(t=72 s)

    It is obvious from Fig.15(c)that a sharp dynamic VBM peak occurs at time=72 seconds when CASE 3 extreme wave is encountering the ship.It is necessary to analyze the component of the dynamic VBM peak.The wave forces are classified by hydrodynamic and hydrostatic force as Eq.(4)explains,the hydrodynamic component includes wave damping load,inertia load and slamming load,the hydrostatic component is decided by interaction between ship weight and buoyancy.CASE 3 is used to analyze dynamic VBM components,modified hydroelastic VBM and hydroelastic VBM have the same hydrostatic components,and Fig.16 reveals modified hydroelastic VBM,hydroelastic VBM and hydrostatic VBM of CASE 3.It is obvious that no matter modified hydroelastic VBM peak or hydroelastic VBM peak is much larger than the peak of hydrostatic VBM,it indicates that hydrodynamic component plays a dominant role to result in the top of dynamic VBM when extreme wave of CASE 3 is encountering to the ship model.Fig.17 shows models of the hydroelasticity and modified hydroelasticity forextreme wave CASE 3 when time is equal to 72 seconds.By comparing the results obtained with the modified hydroelastic model,a large deformation is generated for hydroelastic-plastic model while the hydroelastic model produces an inconspicuous deformation.In addition,a significant slamming phenomenon happens at the bow of ship,it means a giant impacting force would be generated in upper-forward section,so that the ultimate sagging bending moment is generated at this instant,slamming force is dominant to induce ultimate sagging bending moment,and fatal sagging bending is earlier than hogging bending.

    Fig.18 Rigidity distribution of three extreme wave

    Fig.19 Time-domain modified hydroelastic model in CASE 1

    Fig.18 shows rigidity distribution of three extreme wave when time is equal to 72.0 seconds for above three extreme wave Cases.It is obvious from Fig.18 that the ship rigidity changes greatly along ship length for Cases of Hf/D=1.35 and Hf/D=1.33,the smallest rigidity appears around the midship due to ship collapse is happened,while Case of Hf/D=0.7 does not change the rigidity of ship so that only elastic response of ship is resulted.

    The Case of CASE 1 was studied carefully to determine the deformational course of the ship model when struck by the extreme wave.As we know that the CASE 3 wave model consists of a regular wave in its beginning and end parts and an extreme wave in the middle part, the responses for both a regular and extreme wave are obtained.Fig.19 shows a time-domain modified hydroelastic model for the CASE 3 wave.Fig.19 presents six images of the modified hydroelastic model at different times(t=67 s,71 s 72 s,76 s,and 78 s).The images for t=67 s show the model for the regular wave,while the extreme wave begins to pass the ship model, which is elevated at the bow when t=71 s,such that a large deformation is generated when t= 72 s,and then the extreme wave arrives at the stern of the ship,which falls at t=76 s.At t=78 s,the extreme wave has passed the ship model.

    4 Discussion on convergence of the 2D method of m odified hydroelasticity

    4.1 Discussion on convergence of the method regarding to element number

    The calculation results are verified for the method of modified hydroelasticity by discussing its convergence.The way verifying the convergence of the modified hydroelasticity is to change the number of subdivisions of ship beam. Fig.20 shows the number of ship subdivisions (N=10,20,30,40 and 50).Relative acceleration convergence rate as Eq.(25)is used to ac-count for the convergence of modified hydroelasticity.Fig.21 shows the relative acceleration convergence rates as the number of ship subdivisions is changed from 10 to 50 while extreme wave model is CASE 3 which wave parameters are revealed in Tab.2.Fig.21 reveals three curves of relative acceleration convergence rates when time is equal to 67.8 s,72.0 s and 91.9 s.It is observed from Fig.21 that relative acceleration convergence rates are reducing as the number of ship subdivisions is increasing,and the relative acceleration convergence rate is in close proximity to zero when the number of ship subdivision is 50.According to the tendency of three curves shown in Fig.21,it is indicated that the calculation model with large number of ship subdivision has smaller relative acceleration convergence rate.It means that the convergence of the modified hydroelasticity is affirmative.Fig.22 shows the ratio of VBM at midship to sagging ultimate bending moment Ms when the number of ship subdivisions is changed.Also three time steps(67.8 s,72.0 s and 91.9 s)are chosen to be shown.Nearly zero VBM,regular VBM peak and extreme maximum VBM peak are induced at these three time steps.It is observed that nearly zero values are resulted when the number of ship subdivisions is changed at time step of 67.8 s,VBM of time step of 72.0 s has a significant increasing since the number of ship subdivisions from 10 to 20 but the ratio of VBM to Ms almost keeps to 1.0 after the number of ship subdivisions is 20.The VBM of time step of 91.9 s has a very tiny growth since the number of ship subdivisions is increased to 20.It is indicated that the number of ship subdivisions has a tiny effect on VBM results since it is increased to 20,and extreme wave model of CASE 3 leads to sagging ultimate bending moment after N is equal to 20.Therefore, the number of ship subdivisions is determined as 20 in this paper based on the consideration of computational cost of the code.

    Fig.20 Number of ship subdivisions

    Fig.21 Relative convergence rate as number of ship elements is changed

    Fig.22 VBM/Ms as number of ship elements is changed

    4.2 Discussion on convergence of the method regarding to mode number

    In addition,the convergence of the 2D modified hydroelasticity regarding to the number of mode is discussed.The way discussing the convergence of the modified hydroelasticity is to change the number of ship mode.In this paper,the numbers of ship modal involving Nm=1,2, 3,4,5,7,10,15 are used to discuss the applicability of modal method while extreme waveCase is fixed at CASE 3.The Nm-th modal shapes are superposed to obtain ship deformation and ship’s rigidity is varying.Each relative acceleration convergence rate is recorded when the number of ship modal involving Nmis equal to 1,2,3,4,5,7,10 and 15.Fig.23 shows the relative acceleration convergence rates as the number of ship mode is changed to 1,2,3,4, 5,7,10 and 15.It is observed from Fig.23 that the relative convergence rate is notched up small increases when the number of mode is increasing,the relative convergence rate is not sensitive if the number of modal is larger than 10.On the whole,the relative convergence rate increases with the number of modal.Similarly,the ratio of VBM at midship to sagging ultimate bending moment Ms regarding to the number of modal is discussed.Fig.24 shows that the ratio of VBM at midship to sagging ultimate bending moment Ms when the number of mode is changed.It is obvious from Fig.24 that the Max.VBM/Ms keeps constantly as long as the number of mode is larger than 3 though it is increasing slightly when number of mode is between 1 and 3.It is indicated that the number of mode has no effect on maximum VBM as long as it is larger than 3.It means that former three modal shapes are dominant for ship deformation regarding to extreme wave CASE 3,because the CASE3 in which wave length is close to ship length is significant for global ship deformation instead of local deformation simulated by larger modes.

    Fig.23 Relative convergence rate as number of mode is changed

    Fig.24 VBM/Ms as number of mode is changed

    5 Structural evaluation of 500TEU containership in extreme waves

    In order to be applicable to ship designing,as many as wave parameters are studied as well to get general conclusions.The main selection for extreme wave case is supposed to cover most extreme wave patterns,four important extreme wave patterns are considered such as extreme wave height,regular wave height,wave speed and wave length.The four extreme wave patterns are changed to study their influence on nonlinear dynamic strength and deformation of ship.Extreme wave height means the maximum wave height as the mid-part of Fig.12 shows, extreme wave height ratio(Hf/D)which is the ratio of extreme wave height(Hf)to moldeddepth(D)is used to denote the extreme height level.

    5.1 Extreme wave pattern for extreme wave height

    To study the dynamic strength and evaluate the dynamic pre-ultimate strength,the influence of the wave/height ratio(Hf/D)on the ship’s strength and deformation was studied.In this computation,the regular wave height,Hr,has been determined to be 2.8 m and the Froude number of the encountered wave is determined to be 0.3.According to the computational results obtained with the modified hydroelasticity and hydroelasticity methods,Fig.25 shows the maximum VBM/Ms curve as Hf/D is increasing.The ratio of the maximum VBM to Ms indicates the strength of the ship,Ms is sagging ultimate bending moment.All of the points on the blue curve were calculated using the hydroelasticity method,while the points on the red curve were computed using the modified hydroelasticity method with structural nonlinearity taken into account.It is obvious that there is a great difference between them.It indicates that an actual nonlinear dynamic VBM has an ultimate value,and does not increase infinitely,as the height of an extreme wave height increases.It can be seen that,as long as Hf/D is equal to or greater than 1.35,the modified hydroelastic model reaches an ultimate strength and does not increase infinitely.When Hf/D is around 1.55,the maximum VBM of the hydroelastic model exceeds the value of the ultimate bending moment.In the range of Hf/D from 1.3 to 1.55,there is a difference in spread between the two methods,nonlinear VBM calculated by modified hydroelasticity is larger than linear VBM calculated hydroelasticity.

    Fig.25 Max.VBM/Mu curves between hydroelastcity and hydroelasticity-plasticity method

    Fig.26 Maximum angles by changing extreme wave height ratio

    The difference between the two methods is compared based on another viewpoint of deformational angle,that is taken the ship’s deformation into account.The deformational angle at the midship section(S.S.5.0)was computed.Therefore,the influence of Hf/D on the angle at the midship section was also examined.Fig.26 shows the angle curve at S.S.5.0 as Hf/D is increasing.The blue points denote the hydroelastic computations of angle,and the red points represent the modified hydroelastic computations of the angle.It is obvious that there is a difference in the angle since Hf/D exceeds 1.32.It can be seen that a large deformation is generated with the modified hydroelastic model but that the maximum angle of the hydroelastic model increases almost linearly with small amplitude.

    5.2 Extreme wave pattern for regular wave height

    An extreme wave pattern for regular wave height(Hr)is studied to investigate its influence on ship ultimate strength and deformation.The regular wave height,whose interval is 0.5 m,is increased between 0.5 m and 5.5 m while other wave patterns are determined as Fr=0.3, and Lw/L=1.In order to investigate the difference of maximum extreme wave height on ship strength,three extreme wave heights(Hf=11.8,10.32,and 6.0 m)are arranged.Fig.27 shows Max.VBM/Ms curves by changing regular wave height while wave Froude number is 0.3 and wave length ratio Lw/L=1,curves of a1 and a2 mean that their extreme wave height is 11.8 m, curves of b1 and b2 denote that their extreme wave height is 10.32 m,curves of c1 and c2 have extreme wave height as 6.0m,red curves denote that their results are calculated by modified hydroelasticity method while blue curves are calculated by hydroelastic method.It is observed that the extreme waves whose maximum extreme wave is 11.8 m and 10.32 m induces the difference between two methods,it implies that large deformation and nonlinear dynamic strength are resulted in,while the extreme wave(Hf=6.0 m)manifests that no difference between two methods,it means that no large deformation is occurred.Although a small increasing tendency of maximum VBM when regular is increased,maximum extreme wave height plays a dominant role on the large deformation and nonlinear strength.As Fig.27 shows,the effect of regular wave height on ship strength and deformation is rather finite compared with maximum wave height.If it can play a significant effect on ship strength,the maximum wave height needs rather enough amplitude.

    Fig.27 Maximum VBM s by changing regular wave height

    Fig.28 Maximum angles by changing regular wave height

    Fig.28 reveals maximum angle at midship while regular wave height is increasing from 0. 5 m to 3.0 m.By comparison between modified hydroelastic and hydroelastic results,the large deformation is generated since regular wave height is 2.0 m for modified hydroelastic computations while the deformational angle is relatively small manifested by hydroelasticity.

    5.3 Extreme wave pattern for encountering wave speed

    The extreme wave pattern for the encountered wave speed was used to study the influence of the encountered wave speed.Froude numbers are changed to study the influence on shipstrength while other wave patterns are determined as(Hf=11.6 m,Hr=2.78 m,Lw/L=1).Fig. 29 shows the max.VBM/Ms points obtained with the hydroelasticity and modified hydroelasticity methods versus the Froude number,which increases from 0.05 to 0.35 while the wave height ratio Hf/D is fixed at 1.35.It can also be observed that the modified hydroelasticity method produces an ultimate bending moment while the hydroelasticity method does not.Furthermore,the modified hydroelasticity method has a smaller Froude number(Fr=0.3)to make structural entering ultimate VBM while the hydroelasticity method has a larger value of around Fr=0.325 to exceed 1.The VBM difference between the two method occurrs since Fr=0.25.This means that the encountered wave speed has a significant influence on ship strength.Particularly,the use of the modified hydroelasticity method incurs a different effect regarding the encountered wave speed,in comparison with the hydroelasticity method,given that the structural nonlinearity is considered.

    The effect of the wave speed on the ship deformation is also taken into account.Fig.30 shows the maximum angle points obtained with the hydroelasticity and modified hydroelasticity methods as the Froude number changes from 0.05 to 0.35.It can be seen that the two methods produce different results,with the deformation increasing gradually up to a large value when the Froude number is 0.25.The ultimate strength is generated when the Froude number is 0.3.

    Fig.29 Maximum VBM s by changing Froude number

    Fig.30 Maximum angles by changing Froude number

    5.4 Extreme wave pattern for wave length ratio

    Fig.31 Maximum VBM s by changing wave length ratio

    The extreme wave pattern of wave length ratio is discussed to study the strength difference when is changed.The wave length ratio Lw/L denotes the wave length divided by ship length,and wave length ratio is changing from 0.4 to 2.0 while maximum wave height is 13 m,regular wave height is 2.78 m,and Froude number is 0. 3.Fig.31 is maximum VBM/Ms curves by changing wave length ratio.It is seen thatmaximum VBM is occurred around Lw/L=1, maximum nonlinear VBM calculated by modified hydroelasticity arrives ultimate sagging bending moment while maximum linear VBM calculated by hydroelasticity does not arrive the value of ultimate sagging bending moment,the great difference of VBM s occurs in the span around Lw/L=1 where it is from the point(Lw/L=0.9)to the point(Lw/L=1. 35),and ultimate bending moment is generated from the point(Lw/L=0.95)to the point (Lw/L=1.1),other range is the same so that only elastic strength is induced.

    Fig.32 shows the maximum angle curves by changing wave length ratio.The wave length ratio is plotted on the horizontal axis,and maximum angle is plotted on the vertical axis,the red curve denotes the maximum angle calculated by modified hydroelasticity while the blue curve means the maximum angle calculated by hydroelasticity,it is obvious that maximum values of angles occur around the Lw/L=1,maximum angles calculated by modified hydroelasticity and hydroelasticity are the same when the wave length ratio is smaller than 0.9 or larger than 1.45,but the angle difference happens when wave length ratio is between 0.9 and 1.45 because the large deformation calculated by modified hydroelasticity appears,and the angles calculated modified hydroelasticity exceed to 1.137°in which ultimate bending moment occurred when the wave length ratio is between 0.95 and 1.1,it also means that the ship is disable to survive between 0.95 and 1.1.

    Fig.32 Maximum angles by changing wave length ratio

    6 Conclusions

    This paper proposed a 2D modified hydroelasticity method which combines Smith method and hydroelasticity method,a number of extreme wave models were selected for modified hydroelastic computations,and a number of computational results were obtained and compared with those obtained with the hydroelastic method.So we can draw the following conclusions:

    (1)The modified hydroelasticity method is able to calculate the dynamic nonlinear strength and structural deformation of a ship beam subject to extreme waves.

    (2)Nonlinear dynamic strength calculated by modified hydroelasticity has a smaller extreme wave height ratio,regular wave height,and wave speed to make ship arrive the value of ultimate strength of ship than the linear static strength calculated by hydroelasticity,and modified hydroelasticity is easier to result in ultimate strength than hydroelasticity when wave length ratio is around 1.0.

    (3)Extreme wave height plays an important role to make ship collapse,while other waveparameters have secondary effects on ship safety because they should combine extreme wave height to determine the ship safety.

    The current method considers nonlinear ship structure considering plasticity and buckling effect,beam model is applied in a 2D modified hydroelastic code.It is future work to take into account of 3D effect of modified hydroelasticity,every realistic structural member of ship will be computed in future.In fact,the different ship model has different relation of rigidity(EIp) and curvature which has a important effect on nonlinear dynamic response of ship in the same extreme wave model,it also is another perspective that ship model(nonlinear moment-curvature relation)is changed when extreme wave is fixed.There must be a lot of interesting results to be discussed for further studies.

    [1]Peliniovsky E,Kharif C.Physical mechanisms of the rogue wave phenomenon[J].European Journal of Mech.B/Fluids, 2003,22:603-634.

    [2]Haver S,Andersen O J.Freak waves:Rare realizations of a typical population or a typical realization of a rare population [C]//Proc.10th International Society of Offshore and Polar Eng Conf.Seattle,USA,1990,3:123-130.

    [3]Waseda T,Rheem C K,Sawamura J,Yuhara T,et al.Extreme wave generation in laboratory wave tank[C].Proc.15th International Society of Offshore and Polar Eng Conf.,2005:1-9.

    [4]Minam i M,Sawada H,Tanizawa K.Study of ship responses and wave loads in freak wave[J].The Japan Society of Naval Architects and Ocean Engineers,2006.

    [5]Yamamoto,Yoshiyuki,Fujino,Masataka,Fukaswa,Toichi.Motion and longitudinal strength of a ship in head sea and the effects of non-linearities[C].Conference of the society of naval architects of Japan.Nov.1978.

    [6]Masaoka K,Okada H.A numerical approach for ship hull girder collapse behavior in wave[C]//ISOPE2003.Honolulu, Hawaii,USA,2003:369-375.

    [7]Iijima K,Kimura K,Xu W,Fujikubo M.Modified hydroelasticity approach to predicting the post-ultimate strength behavior of ship’s hull girder in waves[J].Journal of Marine Science and Technology,2011,16(4):379-389.

    [8]Xu W,Iijima K,Wada R,Fujikubo M.Experimental evaluation of the post-ultimate strength behavior of a ship’s hull girder in waves[J].J Marine Sci.Appl.,2012:34-43.

    [9]Xu W,Iijima K,Wada R,Fujikubo M.Parametric dependencies of the post-ultimate strength behavior of a ship’s hull girder in waves[J].J Mar Sci Technol,2012,17:203-215.

    [10]Masaoka K,Okada H.A numerical approach for ship hull girder collape behavior in waves[C].Proceedings of ISOPE2003, 2003:369-375.

    [11]Pei Z,Iijima K,Wada R,Fujikubo M,et al.Simulation on progressive collapse behaviour of whole ship model under extreme waves using idealized structural unit method[J].Marine Structures,2015:104-133.

    [12]Liu W,Suzuki K,Shibanuma K.Nonlinear dynamic response and structural evaluation of container ship in large extreme waves[J].Journal of Offshore Mechanics and Arctic Engineering,2015:64-73.

    [13]Liu W,Suzuki K,Shibanuma K.Nonlinear dynamic response and strength evaluation of a containership beam in extreme waves based on hydroelasticity-plasticity method[C].Conference Proceedings of the International Society of Offshore and Polar Engineers,2014,4:652-657.

    [14]Liu W,Suzuki K,Shibanuma K.A two-dimensional hydroelastoplasticity method of a container ship in extreme waves[J]. Journal of Offshore Mechanics and Arctic Engineering,2015:84-93.

    基于一個二維修正水彈性方法的集裝箱船體梁動態(tài)前極限強度評估研究

    劉維勤1,裴志勇1,吳衛(wèi)國1,鈴木克幸2

    (1.武漢理工大學(xué)交通學(xué)院海洋工程系,武漢430063;2.東京大學(xué)工學(xué)系研究科系統(tǒng)創(chuàng)成專攻,東京113-8656,日本)

    海上極端波因其巨大的波高常常導(dǎo)致船體的極限破壞。該文提出了一個二維的修正水彈性方法來研究一個集裝箱船船體梁在極端波中的動態(tài)前極限強度。傳統(tǒng)的極限強度評估基于準(zhǔn)靜態(tài)方法,沒有動態(tài)效應(yīng)被考慮。而船體在波浪下的動態(tài)結(jié)構(gòu)響應(yīng)是基于水彈性方法,傳統(tǒng)的水彈性方法并不能計算船體梁的動態(tài)非線性強度。該二維修正的水彈性方法考慮時域波浪和非線性船體梁之間的耦合,將水彈性方法和Smith方法結(jié)合,用Smith方法計算船體梁的剛度,而其剛度與船體梁的強度和變形曲率有關(guān)。所以該時域的非線性剛度被用于修改水彈性方法里的常數(shù)項的結(jié)構(gòu)梁剛度。幾組極端波模型被用以產(chǎn)生船體梁的大變形和非線性動態(tài)垂向彎矩。文中分別采用修正水彈性方法和普通水彈性方法,通過改變四個重要的極端波參數(shù)如極端波最大波高、規(guī)則波的波高、波速和波長等來研究其對船體梁船中處的大變形轉(zhuǎn)角和非線性垂向彎矩的影響,通過采用修正的水彈性方法計算得來的結(jié)果與水彈性方法計算得來的結(jié)果相比較,得到了一些差異和結(jié)論。

    動態(tài)極限強度;二維修正水彈性方法;Smith方法;極端波;非線性垂向彎矩;轉(zhuǎn)角

    TV131.2

    A

    劉維勤(1985-),男,武漢理工大學(xué)交通學(xué)院講師;裴志勇(1974-),男,武漢理工大學(xué)交通學(xué)院副教授;吳衛(wèi)國(1960-),男,武漢理工大學(xué)交通學(xué)院教授;鈴木克幸(1964-),男,日本東京大學(xué)工學(xué)系研究科教授。

    TV131.2

    A

    10.3969/j.issn.1007-7294.2016.09.005

    1007-7294(2016)09-1121-26

    Received date:2016-07-22

    Foundation item:Supported by National Natural Science Foundation of China:the study of Hydroelasto-plasticity Method of Ship Structure(51509197)

    Biography:LIU Wei-qin(1985-),male,lecturer,E-mail:liuweiqin_123@sina.com; PEI Zhi-yong(1974-),female,associate professor.

    猜你喜歡
    武漢理工大學(xué)船體彎矩
    船體行駛過程中的壓力監(jiān)測方法
    零彎矩設(shè)計理論在連續(xù)梁橋中的應(yīng)用研究
    《武漢理工大學(xué)學(xué)報(交通科學(xué)與工程版)》征稿簡則
    《武漢理工大學(xué)學(xué)報(交通科學(xué)與工程版)》征稿簡則
    CFRP-PCPs復(fù)合筋連續(xù)梁開裂截面彎矩計算方法研究
    Lanterne-volant
    鋼-混疊合連續(xù)梁負(fù)彎矩區(qū)計算分析
    板孔式有彎矩平衡梁應(yīng)用技術(shù)及研究
    焊接殘余應(yīng)力對船體結(jié)構(gòu)疲勞強度的影響分析
    焊接(2015年9期)2015-07-18 11:03:51
    幾何形態(tài)和視覺感知的探討
    国产精品综合久久久久久久免费| 欧美日韩精品网址| 在线a可以看的网站| 亚洲色图av天堂| 此物有八面人人有两片| 久久99热这里只有精品18| 香蕉丝袜av| 精品欧美国产一区二区三| 欧美日韩亚洲国产一区二区在线观看| 国产精品99久久久久久久久| 男女午夜视频在线观看| 欧美日韩中文字幕国产精品一区二区三区| 亚洲欧美日韩高清专用| 国产欧美日韩精品一区二区| 欧美不卡视频在线免费观看| 久久亚洲精品不卡| 这个男人来自地球电影免费观看| 一本精品99久久精品77| 精品国产超薄肉色丝袜足j| 国产精品久久久人人做人人爽| 午夜亚洲福利在线播放| 国产v大片淫在线免费观看| 搡老熟女国产l中国老女人| 男人舔女人的私密视频| 天天添夜夜摸| 久久久久性生活片| 香蕉久久夜色| 国产一区二区在线观看日韩 | 日韩精品青青久久久久久| 国产欧美日韩一区二区三| 国产精品,欧美在线| 国产成人精品无人区| 精品无人区乱码1区二区| 看片在线看免费视频| 国产99白浆流出| 国产精品一及| 高清毛片免费观看视频网站| 欧美国产日韩亚洲一区| 亚洲成人免费电影在线观看| 男女床上黄色一级片免费看| 久久午夜亚洲精品久久| 国产激情偷乱视频一区二区| 久久久久亚洲av毛片大全| 久久精品综合一区二区三区| 国产视频内射| 国产毛片a区久久久久| 亚洲自拍偷在线| 亚洲性夜色夜夜综合| 亚洲人成电影免费在线| 看片在线看免费视频| 一区二区三区高清视频在线| 精品乱码久久久久久99久播| 后天国语完整版免费观看| 国产高清三级在线| ponron亚洲| 亚洲美女黄片视频| 舔av片在线| 午夜激情福利司机影院| 国产伦在线观看视频一区| 变态另类成人亚洲欧美熟女| 成人三级黄色视频| 男女之事视频高清在线观看| 午夜精品在线福利| 国产精品野战在线观看| 亚洲乱码一区二区免费版| 一个人免费在线观看的高清视频| 精品无人区乱码1区二区| 国产亚洲精品综合一区在线观看| 国产野战对白在线观看| 免费搜索国产男女视频| 国产精品一区二区免费欧美| 午夜福利高清视频| 国产精品av久久久久免费| 久久久久性生活片| 白带黄色成豆腐渣| 亚洲av免费在线观看| 男人的好看免费观看在线视频| 国产 一区 欧美 日韩| 变态另类成人亚洲欧美熟女| 免费高清视频大片| 国产真实乱freesex| 美女被艹到高潮喷水动态| 精品电影一区二区在线| 听说在线观看完整版免费高清| 中文字幕精品亚洲无线码一区| 国产1区2区3区精品| 国内精品美女久久久久久| 天堂动漫精品| 精品无人区乱码1区二区| 亚洲性久久影院| 成人三级黄色视频| 插阴视频在线观看视频| 日韩一区二区视频免费看| 91午夜精品亚洲一区二区三区| 国产亚洲午夜精品一区二区久久 | 国模一区二区三区四区视频| 国产成人精品一,二区| 亚洲电影在线观看av| 婷婷六月久久综合丁香| 国产片特级美女逼逼视频| 啦啦啦韩国在线观看视频| 国产爱豆传媒在线观看| 淫秽高清视频在线观看| 午夜福利成人在线免费观看| av.在线天堂| 五月玫瑰六月丁香| 国产视频内射| 伦理电影大哥的女人| 日韩制服骚丝袜av| 成人av在线播放网站| 午夜视频国产福利| 久久精品久久久久久噜噜老黄 | 永久网站在线| 嫩草影院精品99| 成人鲁丝片一二三区免费| 国产成人freesex在线| 日韩精品有码人妻一区| 大香蕉97超碰在线| 91精品一卡2卡3卡4卡| 日韩欧美精品v在线| 99国产精品一区二区蜜桃av| 亚洲在线自拍视频| 又粗又爽又猛毛片免费看| 噜噜噜噜噜久久久久久91| 欧美3d第一页| 最近中文字幕2019免费版| 五月玫瑰六月丁香| 中文字幕精品亚洲无线码一区| 国产精品99久久久久久久久| 97超视频在线观看视频| 国产成人aa在线观看| kizo精华| 少妇丰满av| 亚洲精品国产av成人精品| 大又大粗又爽又黄少妇毛片口| 久久精品人妻少妇| 日日摸夜夜添夜夜爱| 国产av一区在线观看免费| 欧美日韩一区二区视频在线观看视频在线 | 人妻少妇偷人精品九色| 欧美激情久久久久久爽电影| av又黄又爽大尺度在线免费看 | 欧美最新免费一区二区三区| 久久精品夜夜夜夜夜久久蜜豆| 国产精品久久视频播放| 久久久久久久久中文| 国产午夜精品久久久久久一区二区三区| 日韩成人av中文字幕在线观看| 久久久国产成人精品二区| 国语自产精品视频在线第100页| 亚洲一级一片aⅴ在线观看| 麻豆成人午夜福利视频| 国内精品一区二区在线观看| 亚洲熟妇中文字幕五十中出| 99久国产av精品| 一级黄色大片毛片| 九九热线精品视视频播放| 欧美成人精品欧美一级黄| 三级男女做爰猛烈吃奶摸视频| 精品午夜福利在线看| 午夜激情福利司机影院| 欧美日韩综合久久久久久| 国产成人免费观看mmmm| 99热全是精品| 中文亚洲av片在线观看爽| 九色成人免费人妻av| 国产亚洲一区二区精品| 天堂中文最新版在线下载 | 久久欧美精品欧美久久欧美| 18禁裸乳无遮挡免费网站照片| 国产精品野战在线观看| 亚洲美女视频黄频| 国产真实乱freesex| 久久久久久久久久成人| 国产亚洲91精品色在线| 特大巨黑吊av在线直播| 91午夜精品亚洲一区二区三区| 国产高清有码在线观看视频| 黄色欧美视频在线观看| 26uuu在线亚洲综合色| 午夜福利视频1000在线观看| 日韩国内少妇激情av| 最近最新中文字幕免费大全7| 水蜜桃什么品种好| .国产精品久久| 日本猛色少妇xxxxx猛交久久| 干丝袜人妻中文字幕| 国产老妇女一区| 又粗又硬又长又爽又黄的视频| 国产精品一二三区在线看| 久久久久久久久久久免费av| 国产亚洲午夜精品一区二区久久 | 久久久久性生活片| 亚洲欧美精品自产自拍| 免费大片18禁| 日本一本二区三区精品| 精华霜和精华液先用哪个| 亚洲国产最新在线播放| 国产精品99久久久久久久久| 天堂√8在线中文| 成人二区视频| 在线观看av片永久免费下载| 岛国在线免费视频观看| 国产免费又黄又爽又色| 草草在线视频免费看| 久久久国产成人免费| 97超碰精品成人国产| 99热网站在线观看| 特大巨黑吊av在线直播| 日本五十路高清| 午夜福利在线观看吧| 人人妻人人看人人澡| 婷婷六月久久综合丁香| 男人狂女人下面高潮的视频| 久久久亚洲精品成人影院| 丝袜美腿在线中文| 国产亚洲精品av在线| 久久久成人免费电影| 一级黄片播放器| 在线观看66精品国产| 91久久精品国产一区二区三区| 精品99又大又爽又粗少妇毛片| 边亲边吃奶的免费视频| 亚洲国产最新在线播放| 岛国毛片在线播放| 国产高清国产精品国产三级 | av又黄又爽大尺度在线免费看 | 91午夜精品亚洲一区二区三区| 毛片一级片免费看久久久久| 日日摸夜夜添夜夜添av毛片| 久久人人爽人人爽人人片va| 亚洲欧美精品综合久久99| 我的老师免费观看完整版| 精品免费久久久久久久清纯| 波多野结衣巨乳人妻| 精品一区二区免费观看| 91精品国产九色| 国产女主播在线喷水免费视频网站 | 国语对白做爰xxxⅹ性视频网站| 免费av观看视频| 色综合亚洲欧美另类图片| 最近的中文字幕免费完整| 免费av不卡在线播放| 亚洲最大成人av| 国产大屁股一区二区在线视频| 国产一区二区三区av在线| 好男人在线观看高清免费视频| 黄色配什么色好看| 偷拍熟女少妇极品色| 国产 一区 欧美 日韩| 少妇人妻一区二区三区视频| 国产精品一区www在线观看| 免费av不卡在线播放| 18禁在线播放成人免费| 看黄色毛片网站| 色综合亚洲欧美另类图片| 久久精品久久久久久久性| 国产免费福利视频在线观看| 亚洲国产最新在线播放| 国产黄片美女视频| 校园人妻丝袜中文字幕| 日日摸夜夜添夜夜添av毛片| 丰满乱子伦码专区| 男女那种视频在线观看| 日韩精品有码人妻一区| 18禁动态无遮挡网站| 日本黄大片高清| 少妇猛男粗大的猛烈进出视频 | 国产精品麻豆人妻色哟哟久久 | 国产伦精品一区二区三区视频9| 国产精品久久久久久精品电影小说 | 欧美一区二区亚洲| 国产精品野战在线观看| 国产国拍精品亚洲av在线观看| 99热这里只有精品一区| 亚洲成人中文字幕在线播放| 99热全是精品| 亚洲精品日韩在线中文字幕| 国产精品久久久久久精品电影小说 | 少妇裸体淫交视频免费看高清| 精品久久久噜噜| 国产伦精品一区二区三区视频9| 成人鲁丝片一二三区免费| 国内少妇人妻偷人精品xxx网站| 性色avwww在线观看| 天天一区二区日本电影三级| 超碰av人人做人人爽久久| 欧美潮喷喷水| 成人国产麻豆网| 婷婷色av中文字幕| 亚洲乱码一区二区免费版| 中文字幕av在线有码专区| 亚洲av免费在线观看| 欧美成人午夜免费资源| av免费观看日本| 人妻系列 视频| 免费观看在线日韩| 国产精华一区二区三区| 白带黄色成豆腐渣| 色播亚洲综合网| 国产精品人妻久久久久久| 99久久中文字幕三级久久日本| 午夜福利成人在线免费观看| 日韩中字成人| 久久草成人影院| 婷婷色麻豆天堂久久 | 成人三级黄色视频| 久久精品国产99精品国产亚洲性色| 成人毛片a级毛片在线播放| 日韩欧美精品免费久久| 国产精品1区2区在线观看.| 日韩成人av中文字幕在线观看| 久久久久久久久中文| 国产国拍精品亚洲av在线观看| 国产亚洲av片在线观看秒播厂 | 九九久久精品国产亚洲av麻豆| 亚洲一区高清亚洲精品| 只有这里有精品99| 久久99热6这里只有精品| 久久精品人妻少妇| 麻豆成人午夜福利视频| 久久久久久九九精品二区国产| 国产亚洲av嫩草精品影院| 国产黄片美女视频| 好男人视频免费观看在线| 九草在线视频观看| 免费搜索国产男女视频| 国语自产精品视频在线第100页| 免费观看性生交大片5| 日韩成人伦理影院| 一级毛片久久久久久久久女| 国产免费男女视频| 精品免费久久久久久久清纯| 成人高潮视频无遮挡免费网站| 听说在线观看完整版免费高清| 少妇的逼水好多| 男女边吃奶边做爰视频| 成人午夜精彩视频在线观看| 精品久久久久久久久av| 美女被艹到高潮喷水动态| 国产亚洲91精品色在线| 一级毛片aaaaaa免费看小| 日韩制服骚丝袜av| 老女人水多毛片| 亚洲av.av天堂| 亚洲精品一区蜜桃| 日日撸夜夜添| 99热这里只有是精品50| 国产男人的电影天堂91| 九九爱精品视频在线观看| 亚洲第一区二区三区不卡| 亚洲五月天丁香| 色综合色国产| 岛国在线免费视频观看| 看十八女毛片水多多多| av视频在线观看入口| 欧美一级a爱片免费观看看| 国产不卡一卡二| 干丝袜人妻中文字幕| 老司机影院毛片| 成人漫画全彩无遮挡| 日韩欧美在线乱码| 乱人视频在线观看| 尾随美女入室| 精品免费久久久久久久清纯| 精品酒店卫生间| 日日摸夜夜添夜夜爱| 国产毛片a区久久久久| av视频在线观看入口| 人妻制服诱惑在线中文字幕| av卡一久久| 女人久久www免费人成看片 | 日本免费一区二区三区高清不卡| 三级毛片av免费| 日韩av在线免费看完整版不卡| 午夜爱爱视频在线播放| 久久久午夜欧美精品| 秋霞在线观看毛片| 深夜a级毛片| 国产v大片淫在线免费观看| 97人妻精品一区二区三区麻豆| 国产一级毛片在线| 久久久久精品久久久久真实原创| 欧美成人午夜免费资源| 我的女老师完整版在线观看| 久久久久久久亚洲中文字幕| 国产成年人精品一区二区| 天堂网av新在线| 国产成年人精品一区二区| 亚洲在线观看片| АⅤ资源中文在线天堂| 国产成人午夜福利电影在线观看| 在线观看美女被高潮喷水网站| 尾随美女入室| 日日摸夜夜添夜夜爱| 三级经典国产精品| 国产在视频线精品| 精品一区二区三区人妻视频| 在线免费十八禁| 久久久久久九九精品二区国产| 国产黄色视频一区二区在线观看 | 亚洲精品乱码久久久v下载方式| 热99在线观看视频| 男女那种视频在线观看| 亚洲欧洲日产国产| 国产成人一区二区在线| 国产免费视频播放在线视频 | 色吧在线观看| 亚洲在线观看片| 一边亲一边摸免费视频| 亚洲欧美成人精品一区二区| 日本免费a在线| 国产毛片a区久久久久| 性插视频无遮挡在线免费观看| 精品不卡国产一区二区三区| 国产精品嫩草影院av在线观看| 精品一区二区三区人妻视频| 非洲黑人性xxxx精品又粗又长| 婷婷色麻豆天堂久久 | 日韩大片免费观看网站 | 色综合亚洲欧美另类图片| 老师上课跳d突然被开到最大视频| 99热网站在线观看| 国产淫语在线视频| or卡值多少钱| 岛国毛片在线播放| 男人的好看免费观看在线视频| 国内精品美女久久久久久| 国产精品美女特级片免费视频播放器| 村上凉子中文字幕在线| 亚洲av电影不卡..在线观看| 亚洲一级一片aⅴ在线观看| 国产探花在线观看一区二区| 亚洲精品国产成人久久av| 国产伦理片在线播放av一区| 亚洲国产精品成人综合色| 中文字幕亚洲精品专区| 草草在线视频免费看| 国产 一区 欧美 日韩| 天堂√8在线中文| 久久久欧美国产精品| 少妇裸体淫交视频免费看高清| 国产精品一区二区三区四区免费观看| 春色校园在线视频观看| 久久久亚洲精品成人影院| 久久精品久久久久久久性| 久久精品国产自在天天线| 99久久成人亚洲精品观看| 欧美精品国产亚洲| 久久久欧美国产精品| 嫩草影院新地址| 国产精品1区2区在线观看.| 国产亚洲av嫩草精品影院| 久久久亚洲精品成人影院| 久久久国产成人精品二区| 国产精品av视频在线免费观看| 又黄又爽又刺激的免费视频.| 成年女人永久免费观看视频| 国产成人freesex在线| av国产久精品久网站免费入址| 亚洲中文字幕一区二区三区有码在线看| 乱系列少妇在线播放| 久久99蜜桃精品久久| 欧美3d第一页| 只有这里有精品99| 少妇人妻精品综合一区二区| 中文字幕av成人在线电影| 亚洲人与动物交配视频| 国产高潮美女av| 日本欧美国产在线视频| 久久久久国产网址| 久久精品久久精品一区二区三区| 亚洲性久久影院| 久久久精品欧美日韩精品| 欧美变态另类bdsm刘玥| 长腿黑丝高跟| 高清在线视频一区二区三区 | 精品久久久久久久久av| 尤物成人国产欧美一区二区三区| 村上凉子中文字幕在线| 国产黄色视频一区二区在线观看 | 国产亚洲av嫩草精品影院| 亚洲无线观看免费| 少妇的逼好多水| videos熟女内射| 国产伦精品一区二区三区视频9| av女优亚洲男人天堂| 国产精品久久久久久久电影| 乱人视频在线观看| 真实男女啪啪啪动态图| 国产精品av视频在线免费观看| 成人二区视频| 亚洲真实伦在线观看| 熟女人妻精品中文字幕| 欧美日韩综合久久久久久| 欧美性感艳星| a级一级毛片免费在线观看| 日韩大片免费观看网站 | 精品不卡国产一区二区三区| 最近中文字幕高清免费大全6| 精品酒店卫生间| 精品欧美国产一区二区三| 黄片无遮挡物在线观看| 高清毛片免费看| 国产av不卡久久| 搞女人的毛片| 国产精华一区二区三区| 欧美日本亚洲视频在线播放| 乱码一卡2卡4卡精品| 狂野欧美激情性xxxx在线观看| 大香蕉久久网| 午夜精品在线福利| 97人妻精品一区二区三区麻豆| av在线老鸭窝| 亚洲av一区综合| 全区人妻精品视频| 建设人人有责人人尽责人人享有的 | 亚洲av熟女| 国产免费福利视频在线观看| 久久久国产成人精品二区| 高清日韩中文字幕在线| 国产精品不卡视频一区二区| 最后的刺客免费高清国语| 久久久午夜欧美精品| 久久久久久久久久久丰满| 97在线视频观看| 黄色一级大片看看| 久久久色成人| 3wmmmm亚洲av在线观看| 联通29元200g的流量卡| 丝袜美腿在线中文| 在线观看av片永久免费下载| 亚洲不卡免费看| 久久人人爽人人爽人人片va| 看免费成人av毛片| 亚洲国产色片| 日韩精品有码人妻一区| 国产在线男女| 久久午夜福利片| 久久久亚洲精品成人影院| 亚洲不卡免费看| 久久久久久大精品| 黄色日韩在线| 非洲黑人性xxxx精品又粗又长| 久久99热这里只有精品18| 国产成人午夜福利电影在线观看| 亚洲国产最新在线播放| 久久99蜜桃精品久久| 国产一区有黄有色的免费视频 | 能在线免费看毛片的网站| 国产精品国产三级专区第一集| 日韩欧美国产在线观看| 我要搜黄色片| av在线蜜桃| 精品99又大又爽又粗少妇毛片| 男女边吃奶边做爰视频| .国产精品久久| 麻豆精品久久久久久蜜桃| 丰满人妻一区二区三区视频av| 国产精品99久久久久久久久| 日韩亚洲欧美综合| 久久精品久久久久久噜噜老黄 | 一二三四中文在线观看免费高清| 国产精品综合久久久久久久免费| 高清毛片免费看| 中文字幕熟女人妻在线| 高清午夜精品一区二区三区| 男人狂女人下面高潮的视频| 一边亲一边摸免费视频| 人人妻人人澡人人爽人人夜夜 | 亚洲久久久久久中文字幕| 久久久精品94久久精品| 美女脱内裤让男人舔精品视频| 成年免费大片在线观看| 中文精品一卡2卡3卡4更新| 一个人看视频在线观看www免费| 在线免费十八禁| 日本wwww免费看| 久久久久免费精品人妻一区二区| 大香蕉97超碰在线| 26uuu在线亚洲综合色| 久久久久久久久中文| 国产在视频线精品| 五月伊人婷婷丁香| 小蜜桃在线观看免费完整版高清| 亚洲综合精品二区| 婷婷色综合大香蕉| 波多野结衣巨乳人妻| 国产在视频线精品| www日本黄色视频网| 国产在线男女| 春色校园在线视频观看| 少妇被粗大猛烈的视频| 18+在线观看网站| 男人的好看免费观看在线视频| 亚洲成人中文字幕在线播放| 欧美变态另类bdsm刘玥| 欧美激情在线99| 免费观看在线日韩| 久久精品国产亚洲网站| 亚洲乱码一区二区免费版| 日韩大片免费观看网站 | 在线a可以看的网站| 亚洲欧美中文字幕日韩二区| av天堂中文字幕网| 国产成人aa在线观看| av在线蜜桃| 国产单亲对白刺激| 国产精品国产三级国产av玫瑰| 亚洲av成人精品一区久久| 麻豆国产97在线/欧美| 伦精品一区二区三区| 国产精品女同一区二区软件| 国产熟女欧美一区二区| 日本午夜av视频|