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

    Analysis of Local Fracture Strain and Damage Limit of Advanced High Strength Steels using Measured Displacement Fields and FEM

    2015-12-13 03:32:53MaSatoandTakada
    Computers Materials&Continua 2015年6期

    N.Ma,K.Satoand K.Takada

    Analysis of Local Fracture Strain and Damage Limit of Advanced High Strength Steels using Measured Displacement Fields and FEM

    N.Ma1,2,K.Sato3and K.Takada4

    The local mechanical behaviors of advanced high strength steels undergoing a very large strain from uniform plastic deformation to fracture were investigated with the aid of a measured displacement field and a measurement based FEM.As a measurement method,a digital image grid method(DIGM)was developed and the three-direction transient displacement field on uniaxial tensile test pieces was measured.Combining the measured transient displacement field with the finite element method,a measurement based FEM(M-FEM)was developed for the computation of distribution of the local strains,local stresses and ductile damage accumulation in a tensile test piece.Furthermore,the local fracture strain and damage limit of several advanced high strength steel sheets(980MPa/t1.2mm,980MPa/t1.6mm,1180MPa/t1.6mm)were identified by uniaxial tensile tests and the measurement based FEM.The identified damage limit of materials agreed very well compared with that measured by a conventional press test,and the validity of the measuring method and measurement based FEM was verified.

    Digital image grid method,Measured displacement field,Measurement based FEM,Local fracture strain,Damage limit,Advanced high strength steels.

    1 Introduction

    The finite element method(FEM)has been widely used in the analysis of mechanical behaviors of various materials and structures.To know the detail distribution of the local strains and stresses by FEM,the analyzing object is generally divided into very small elements.However,the local strains and stresses in a small element were difficult to be measured by conventional experiments.

    Recently,with the aid of advanced measurement methods such as a digital image correlation method,the local mechanical deformation behaviors in materials can be directly observed by experiments.These directly measuring technologies are very helpful in developing some new computational methods and also in verifying the numerical simulation models.

    A digital image correlation method(DIC)was early developed by Peters et al.(1981),Chu et al.(1985),and Sutton et al.(1986).Pan et al.(2009)reviewed the measuring technology of DIC for strain field in details.Huang et al.(2011)evaluated the anisotropic parameter R-value;Zeng and Xia(2010)determined stressstrain curve beyond uniform elongation.Based on DIC based measuring methods,Ramazani et al.(2013)investigated deformation characterization of failure initiation in a dual phase steel;Paul(2012)predicted failure modes of dual phase steels.Coppieters et al.(2014)measured the large strains using the digital image correlation method and identified post-necking strain hardening phenomena.In this year,Hopmann et al.(2015)measured the strain rate dependent material properties of polymers and Sato et al.(2015)measured the strain rate dependent fracture behaviors of advanced high strength steels,Jian et al.(2015)investigated the damage initiation of dual phase steels under different stress triaxiality.

    Up to now, experimental measurements and numerical simulations using FEM were performed independently.Their seamless combination was not deeply discussed.An early application of their combination was tried by Ueda et al.(1989)who computed the distributions of three dimensional internal welding residual stresses by elastic FEM using measured released strains.Recently,great attentions are being paid to the applications of the displacement field measured by digital image correlation method.Besnard et al.(2006)developed a method to estimate the displacement from a pair of images and evaluated the displacement field for finiteelement analysis.Dupuy et al.(2010)mapped the measured displacement field to finite element model and analyzed elastic bending behaviors on single shear lap joint.Tarigopula et al.(2008)extended DIC to the study of large plastic deformation problem.Furthermore,Roux and Hild(2008)estimated the damage using measured results by DIC.Yoneyama(2011)developed a smoothing scheme for measured displacements and then calculated the strain concentration around a hole in a plate specimen.

    Generally,the digital image correlation method is based on a random pattern painted on the surface of structure samples. The painted random pattern is easy to be broken and accompanied measuring errors may increase when the deformation becomes large.Ma et al.(2014)developed a digital image grid method(DIGM)using a regular pattern with grids etched electrochemically on the surface of a test piece.The grids etched electrochemically on the surface of test pieces can suffer a very large deformation until fracture.The regular patterns used in DIGM can be arranged as a rectangle shape or a triangle shape which is similar to the element shape of FEM.The grids of regular patterns can be used as the nodes of FEM.These grids can also be modelled by particles if the meshless local petrov-galerkin method(MLPG)[Liu et al.(2006),Sladeki et al.(2013)]is employed.Recently,an advanced MLPGEshelby Method was proposed by Han and Atluri(2014)which has advantages to deal with material separation.If the measured displacement fields are combined with MLPG-Eshelby Method in the future work,the detailed fracture behaviors such as crack initiation and propagation will be investigated.

    In the present article,the digital image grid method(DIGM)with regular patterns was used to measure the transient displacement field of tensile test pieces during loading process. Then, the measured grid displacements are considered as the nodal displacements in widely employed FEM,and the strains and stresses can be directly estimated.This combination is here named as a measurement based finite element method(M-FEM).Then applications of M-FEM to the computation of distribution of the local strains,the local stresses and the local ductile damage accumulation in steel materials during tensile tests are represented.Furthermore,the ductile damage limit of several advanced high strength steels was identified by the proposed M-FEM.The identified damage limit agreed well with the value measured by a conventional method and the validity of M-FEM was verified.

    2 Displacement field measurement by digital image grid method

    2.1 Digital image grid method(DIGM)

    The digital image grid method used in the uniaxial tensile test is schematically shown in Fig.1(a).It consists of a two-camera system for static tensile testing which is controlled by a computer.The specifications of the digital cameras used in the experiments are shown in Fig.1(b).The digital images were taken every 2.0 sec in the tensile test.The grids with a diameter of 0.25 mm or 0.5 mm are electrically etched with a certain pitch length named as grid pitch length.The grid pitch length can be 0.5mm or 1.0mm or 2.0mm as shown in Fig.1(c)depending on the measuring requests.The size of grids and grid pitch length are designed by considering the limitation of the degree of etching precision and accuracy for the local strain measurement.The center of the grid was determined by the fitting of elliptic curves from each image.Using the digital image correlation method,the center of the distorted grids was determined accurately even a large distortion occurs after the strain localization.The 3D coordinates of the grids on the specimen were calculated using the center of the fitted ellipses.The three displacement components ux,uy,uzat the grids were calculated from the coordinate changes in the x,y and z directions,which were continuously measured within a certain time interval.

    Figure 1:Measuring system of digital image grid method:(a)Static tensile test machine with two digital cameras,(b)Specifications of the digital camera,(c)conditions of the grid pitch length and grid size.

    In the digital image grid method,the deformation of the etched grids on the specimen is analyzed using a computer vision algorithm(2006).The computer vision algorithm is able to reconstruct a three-dimensional(3D)shape from multiple two-dimensional(2D)images.The relationship between a 3D-coordinate point,M(x,y,z)in the three dimensional global coordinate system and a point m(u,v)in the two dimensional camera coordinate system is represented by Eq.(1).

    where λ is a scale parameter,A is a camera matrix including the focal length(au,av),the center of the image(u0,v0),and a share factor scale,s.The R and T matrixes are the transformation and rotation matrixes defined by Eq.(2)and shown in Figure 2 respectively

    Figure 2:Rotation and translation of the stereo camera model.

    In the Eq.(1),the part of A[R|T]is condensed into a 3×4 matrix,P,as follows

    The parameters of the camera matrix P can be estimated using a pattern with known 3D coordinate points and the projected points in the camera coordinate system shown in Fig.2.Once P is determined,the transformation and rotation of the camera are calculated by the decomposition of matrix P.When the positions and rotations of all cameras are calculated,the 3D position of the object can be calculated from multiple images captured by the cameras.The calculation problem is a linear simultaneous equation with 3D-coordinate parameters x,y and z.The equation could be solved using the inverse matrix of P.

    2.2 Measuring object

    Figure 3 shows a uniaxial tensile test piece of a high strength steel sheet 980MPa with 1.6mm in thickness and the displacement field on the highlight area was measured by the digital image grid method.The marks?are the grids used in DIGM where the displacements are to be continuously measured.

    2.3 Displacement field measured by DIGM

    When the length of the measuring area of a tensile test piece of a steel sheet 980MPa/t1.6mm is pulled from the original 52.7711mm to 61.289mm with a con-trolled velocity 2.5mm/min,the displacement field of component Ux in the tensile direction(x)directly measured by DIGM with grid pitch length 1.0mm is shown in Figure 4.The points marked by?,3 and 2 are the focused points located at edges and center of the tensile test piece.The distribution of the measured displacement component Uy in the transverse direction(y)is shown in Figure 5.The obvious necking in the width direction of a tensile test piece can be observed.The measured necking displacement Uy at the two edges of the width direction(y)has a little difference.The distribution of displacement component Uz in the thickness direction(z)is shown in Figure 6.The through thickness necking can be obviously observed at the center of a test piece in the width direction.

    Figure 7 shows the history of measured displacements Ux,Uy and Uz at the three points marked by○(upEdge),◇(lowEdge)and □(Center),respectively.Among the three displacement components, the tensile directional displacement component Ux is largest.The value of the measured displacement Ux at three points has a little difference.The deformation directions of the displacement Uy between the two points located at the up edge and low edge,respectively,are opposite,and the value at the up edge point is larger than that at the low edge point.The displacement Uz at the two edge points is almost the same and that at the center point is quite different.

    Figure 3:A tensile test piece and the grids for displacement measurement by DIGM.

    From the directly measured displacement fields of components Ux,Uy and Uz shown in Figures 4-7,it can be seen that the distribution of displacements along the width direction of the tensile test piece is not symmetric because the experimental conditions are not as simple as the conditions used in the standard FEM.Therefore,the combination of the directly measured displacement fields with FEM should play a very important role in investigating the practical deformation behaviors occurred in experiments.

    Figure 4:Distribution of measured displacement component Ux(mm)

    Figure 5:Distribution of measured displacement component Uy(mm)

    Figure 6:Distribution of measured displacement component Uz(mm)

    Figure 7:History of measured displacements Ux,Uy and Uz at two edges and center

    3 Measurement based FEM for local strains,stresses and damage computation

    The standard FEM for nonlinear problems is based on Eq.(4)using the stiffness matrix[K]and nodal displacement increment{?u(t)}at timetif external nodal force increment{?F(t)}is applied to an analyzed model.

    If the nodal displacement increment{?u(t)}is obtained by solving Eq.(4),the total displacement{u(t)}at nodes,the strain increment{?ε(t)},the total strain{ε(t)},the stress increment{?σ(t)}and the total stress{σ(t)}in elements can be computed in sequence by Eqs.(5-9).Furthermore,the strength and damage accumulation of materials can be predicted.

    Where,[B]is the strain-displacement matrix depending on the element types used in FEM and[De]is the elastic matrix of material defined by Young’s modules and Poisson’s ratio.

    With the progress of measurement technologies,the nodal displacement increment{?u(t)}in Eq.(4)can be partially or fully measured.If the nodal displacement at partial nodes is measured and expressed by{?ˉu},the nodal displacement vector{?u(t)}of all nodes used in FEM can be classified into the unknown nodal displacement increment{?u}and known nodal displacement{?ˉu}.Therefore,the basic Eq.(4)of the standard FEM can be transferred into Eq.(10)which can be called as a partial measurement based FEM.

    If the displacements{?u}at all the nodes of a FE model are measured,it is unnecessary to solve Eq.(10).The total displacement at nodes,strains and stresses in elements of FE models can be computed by Eqs.(5)-(9)using measured results.This method is here named as a measurement based FEM which can be used for detail evaluation of the local strains and stresses existing in materials.

    Generally,the computation for strains and stresses by standard FEM is performed following the flow shown in Figure 8.If the displacement field at the interesting zones was measured,the computational flow based on the measurement based FEM for local strains and stresses in materials can be simply summarized by Figure 9.The main differences of the measurement based FEM from the standard FEM are highlighted in Fig.9.One is that the mesh used in M-FEM must be corresponding to the measuring points and another is that the nodal displacement increment{?u}used in M-FEM is just the measured value.In the measurement based FEM,the isotropic Mises yield function or the anisotropic Hill yield function can be employed.After the strains and stresses are computed,their damage variables in elements for strength evaluation can be also computed.In this study,a ductile damage originally proposed by Cockcroft-Latham(1968)was mainly discussed,because it has been successfully applied to prediction of the crack occurrence during metal forming by Takuda et al.(2009)and to the impact strength evaluation by Takada et al.(2015)after implemented into commercial software LS-DYNA by Ma(2013).

    Figure 8:Flow of standard FEM.

    Figure 9:Flow of measurement based FEM.

    4 Strain distributions and local fracture strain computed by M-FEM

    4.1 Mesh division

    To use the measurement based FEM,the measuring grids shown in Fig.3 are defined as nodes and rectangle shell elements corresponding the measuring zone are generated as shown in Figure 10.The mesh size is equal to the grid pitch length 1.0mm used in the measurement.A plane stress state in the uniaxial tensile test piece is assumed and a membrane shell element with one integration point through the thickness direction was employed.Therefore,the measured in-plane components Ux and Uy were only applied to nodes of the analyzing model.The strain component εzin the thickness direction of shell elements is computed from in-plane strain components εxand εy.

    Figure 10:Shell elements corresponding to the measuring grids in DIGM.

    4.2 Local strain distributions

    Figure 11 and Figure 12 represent the distribution of the maximum in-plane principal strain(Major strain)and its deformation direction just before crack occurrence.A large tensile deformation band can be obviously observed.The vector direction of the major strain is approximately parallel to the tensile direction.The element located at the center of the large deformation band has a largest principal strain.The maximum local strain is about 0.62 just before cracking.In the other word,the local fracture strain 0.62 measured with 1.0mm gauge length is much larger than the elongation of 980MPa/t1.6mm steel sheet which is evaluated using 50mm gauge length.

    Figure 13 and Figure 14 show the distribution of the minimum in-plane principal strain(Minor strain)and its deformation direction,respectively.A compressive deformation direction of the minor strain is approximately parallel to the width direction of the tensile test piece.

    Figure 15 shows the thickness distribution just before crack occurrence.The center in the width direction of the tensile test piece has a thinnest thickness 1.14mm and the local thickness strain is about-33.4%.

    Figure 11:Major strain distribution.

    Figure 12:Major strain direction.

    Figure 13:Minor strain distribution.

    4.3 Local fracture strain

    From the local strain distributions shown in Fig.11 and Fig.15,it has been known that the crack initiated at the center of the uniaxial tensile test piece.Referring to the FLD expression for crack estimation in sheet metal forming,the local inplane principal strains(ε1,ε2)at the crack initiation element and their historical path up to fracture are shown in Figure 16.When the in-plane principal strainsε1andε2are small,their relation i.e.the strain path is linear and its ratioβis about 0.5 representing the uniaxial stress state.When the strains become larger,the strain path curve changes its direction parallel to the vertical axis and the ratioβof the major strain increment dε1and minor strain increment dε2changes to near zero,which represents the plane strain state.The local fracture strain marked by×measured with the 1.0mm grid pitch length for steel 980MPa/t1.6mm is 0.62.

    4.4 Local necking deformation characteristics

    To investigate the local deformation characteristics,the local strains(εx,εy,εz)at the center of the transverse section where crack initiated in the uniaxial tensile tests and their change with strain-x over 50mm in gauge length which is equal to the tensile stroke if multiplied by 50mm,are shown in Figure 17.The broken red line in

    Figure 14:Minor strain direction.

    Figure 15:Thickness distribution before crack occurrence.

    Figure 16:Local strain path and local fracture strain of 980MPa/t1.6mm.

    Fig.17 represents an assumed relation if the local strain εxwas always equal to the average strain-x.When the average strain-x is less than about 0.08,the local strain εxis approximately equal to the average strain-x and it can be understood that the uniform deformation occurred.Below this value,the three local strain components εx,εyand εzhave a linear relation with the average strain-x,and the difference between εyand εzis small.If a local strain ratio dεy/dεzis given by Eq.(11)and named as a local anisotropic parameter Rlocalto evaluate the deformation characteristics,its value is close to 1.0 and the material has an isotropic deformation feature when the strain is small.

    When the average strain-x is larger than 0.08,the changes of the local strains εx,εyand εzwith the average strain-x become nonlinear.If looking at the changes of the local strains in details,the changing tendency of the local strains εx,εyand εzis different when the average strain-x is between 0.08~0.13 and larger than 0.13.When the average strain-x is between 0.08 and 0.13,the local thickness strain εzchanges a little bit larger than the local strain εy.From another view,the local strain ratio dεy/dεz,i.e.the local anisotropic parameter Rlocalbecomes less than 1.0 and an anisotropic deformation occurred.When the average strain-x is larger than 0.13,the local thickness strain εzchanges much faster than the local strain εy.This is because the through thickness necking occurred.This phenomenon can also be described by the decreasing change of the local anisotropic value Rlocalwith local strain as shown in Figure 18.The Ravin Fig.18 which was measured from 50mm gauge length,is the conventional average anisotropic R-value for description of material anisotropic deformation,and is generally considered as a constant.

    To clearly distinguish diffusion necking in the width direction and local necking in the thickness direction,Figure 19 shows the changes of local thinning strain εzat three locations(Center,upEdge,lowEdge)and average shrinkage strain εyin the width direction of the tensile test piece as drawn by purple broken line(ey:average of GL25mm).The red dotted line(ey:assumed uniform)in this figure was drawn to express the assumed isotropic uniform shrinkage.When the average strain-x with the gauge length GL50mm is less than 0.08,the thinning strain εzat all three locations and average shrinkage strain εyin the transverse direction(GL25mm)overlapped with the assumed isotropic uniform strain εy.When the average strain-x is between 0.8-0.13,the thinning strain εzat the three locations have a little difference.When the average strain-x is larger than 0.13,the thinning strain εzat the center increased much faster than that at the up edge and low edge.This means that the local necking and crack occurred at the center of the tensile test piece.

    Figure 17:Changes of local strains.

    Figure 18:Changes of local anisotropic value.

    Figure 19:Historical change of thinning strain ez and shrinkage strain ey.

    5 Stress distribution computed by M-FEM

    To compute the stresses in the tensile test piece by M-FEM,the nonlinear stress and strain relation under plastic deformation state must be given or identified.Generally,the strain hardening phenomena are coupled with yield function and stress cycles.In this uniaxial tensile test,kinematic hardening is not considered because there is no unloading process.For steel materials,the Hill yield function is often selected for the description of stresses and strains under the plasticity state.Therefore,the Hill yield function with averaged anisotropic value r=0.94 for a 980MPa/t1.6mm steel sheet is here employed.Furthermore,a uniaxial tensile test piece as shown in Figure 20(a)was modeled by finite element with 1.0mm in mesh size as shown in Figure 20(b),and then stress-strain curve was identified by fitting the computed force-stroke curve with the measured force-stoke curve.The stroke is de fined by the change of the 50mm gauge length(GL50mm)as shown in both Fig.20(a)and Fig.20(b).Figure 19(c)shows the plastic strain distribution computed by FEM and a very large strain concentrated at the center of the tensile specimen.Figure 21 represents the identified stress-strain curve based on well used swift hardening rule in which the initial yield stress,the strength coefficient and the strain hardening exponent are 524MPa,1459MPa and 0.104,respectively.Using the identified stress-strain curve,the force-stroke curve(F_fem)reproduced by FEM was compared with the experimental curve(F_test)as shown in Figure 22.It can be seen that a good agreement between FEM and experimental measurement was obtained.Therefore,the identified stress-strain curve can be accepted for the computation of local stresses.

    Figure 20:Identification method of stress-strain relation.

    Figure 21:Identified stress-strain relation of 980MPa/t1.6mm steel sheet.

    Figure 22:Comparison of force-stroke curves between FEM and experimental test.

    Figure 23:Equivalent plastic strain.

    Figure 24:The maximum principal stress.

    Figure 23 and Figure 24 show the distributions of equivalent plastic strain and the maximum principal stress computed by M-FEM.From computed results,it can be easily observed that a large plastic strain and principal tensile stress exist at the center of the width direction of tensile test piece.Although the distributions of stress and strain in the low stress zones are not smooth enough due to the measuring errors,their effect is little in the large strain and high stress zone.The results can be considered reliable in general.

    6 Damage distribution and damage limit identified by M-FEM

    When both the stresses and the plastic strain are computed by M-FEM,various ductile damage values can be calculated following their damage rules.Here,the ductile damage C-value proposed by Cockcroft and Latham given by Eq.(12)is computed in which εf,σ1and dˉεpare the local fracture strain,the 1stprincipal stress and the equivalent plastic strain increment,respectively.

    Figure25shows the distribution of the ductile damage value in the tensile test piece.The maximum damage C-value is produced at the center of the large deformation band.This can be understood that the crack may be initiated at the center position more than the edges of the tensile test piece.The crack initiation position was also observed by experiments as shown in Figure 26,which verified the computed results of measurement based FEM.If the fracture section is observed in details,it can be seen the fracture face at the crack initiation position is perpendicular to the tensile direction.

    When both the damage accumulation in the center element where the crack initiated shown in Fig.25 and the local tensile strain in the same element shown in Fig.11,are drawn by a graph,the historical change of the damage accumulation with the local strain during the tensile test can be represented by Figure 27.It can be seen that the damage limit of high strength steel 980MPa/t1.6mm and the local fracture strain are 870MPa,0.62,respectively.

    Figure 25:Ductile damage value.

    Figure 26:Crack initiation position.

    Figure 27:Ductile damage accumulation at crack initiation element and its limit.

    7 Local fracture strain and damage limit of various advanced high strength steels

    As shown in Figure 28,the local strain path and local fracture strain for three advanced high strength steel sheets(1180MPa/t1.6mm, 980MPa/t1.2mm,980MPa/t1.6mm)were measured with different grid pitch lengths(GL=2.0mm,1.0mm,0.5mm)in DIGM.Then,the damage limit was identified by the measurement based FEM.The mark“+”in this figure represents the local fracture strains when the crack initiated at the center of the tensile test pieces.The measured local fracture strain and identified damage limit for the three advanced high strength steel sheets(1180MPa/t1.6mm,980MPa/t1.2mm, 980MPa/t1.6mm) are shown in Table 1.With the decreasing of grid pitch length,the local fracture strain and ductile damage limit increased.The damage limit measured based on a conventionally complicated Marciniak press test(1967)was also represented in Tab.1.The damage limit identified by the measurement based FEM agreed well with that measured by a conventional press test.The measuring method DIGM and developed measurement based FEM are accurate and reliable in analysis of the local mechanical behaviors such as local strains,local stresses and local damage of materials.

    Figure 28:Local strain path and local fracture strain of advanced high strength steels.

    Table 1:Local fracture strain and damage limit of advanced high strength steel sheets.

    8 Summaries

    (1)Displacement field in the uniaxial tensile test piece of advanced high strength steel sheets was successfully measured by digital image grid method(DIGM).

    (2)A measurement based finite element method(M-FEM)for investigation of local strains,stresses and local damage of materials was developed.

    (3)Local strains,stresses and damage distribution in uniaxial tensile test piece of 980MPa steel were computed by developed measurement based FEM.

    (4)The maximum value of local principal strain,principal stress and ductile damage accumulation existed at the center of the width direction of the tensile test piece which was consistent with the crack occurrence position observed in experiments.

    (5)The local fracture strain and damage limit of advanced high strength steel sheets(980MPa/t1.6mm,980MPa/t1.2mm,1180MPa/t1.6mm)were accurately measured by a uniaxial tensile test and identified by measurement based FEM compared with those measured by a conventional press test.

    (6)The effect of grid pitch length on the ductile damage limit was investigated if different mesh sizes of FEM are employed in strength prediction.With the decreasing of grid pitch length,the ductile damage limit increased.

    Chu,T.C.;Rasson,W.F.;Sutton,M.A.;Peters,W.H.(1985):Applications of digital-image-correlation techniques to experimental mechanics.Experimental Mechanics,vol.25,pp.232–244.

    Cockcroft,M.G.;Latham,D.J.(1968):Ductility and the workability of metals.Journal of the Institute of Metals,vol.96,pp.33-39.

    Coppieters,S.;Ichikawa,K.;Kuwabara,T.(2014):Identification of Strain

    Hardening Phenomena in Sheet Metal at Large Plastic Strains..Procedia Engineering,vol.81,pp.1288–1293.

    Dupuy,J.S.;Lachuad,F.;Piquet,R.;Huet,J.(2010):Finite element model matching based on optical measurement fields on single shear lap joint.Annual conference proceedings of society for experimental mechanics,June 70-10,Indianapolis,Indiana,USA,53-62.

    Han,Z.D.;Atluri,S.N.(2014):Eshelby Stress Tensor T:a Variety of Conservation Laws for T in Finite Deformation Anisotropic Hyperelastic Solid&Defect Mechanics,and the MLPG-Eshelby Method in Computational Finite Deformation Solid Mechanics-Part I.CMES,vol.97,no.1,pp.1-34.

    Han,Z.D.;Atluri,S.N.(2014):Onthe(Meshless Local Petrov-Galerkin)MLPGEshelby Method in Computational Finite Deformation Solid Mechanics-Part II.CMES,vol.97,no.3,pp.199-237.

    Hopmann,C.;Klein,J.(2015):Determination of strain rate dependent material data for FEA crash simulation of polymers using digital image correlation.Computational materials science,100,Part B,April,pp.181–190.

    Huang,G.;Yan,B.;Xia,Z.(2011):Measurement of r-values of High Strength Steels Using Digital Image Correlation.SAE Int.J.Mater.Manuf.,vol.4,no.1,pp.385-395.

    Lian,J.;Yang,H.;Vajragupta,N.;Münstermann,S.;Bleck,W.(2014):A method to quantitatively upscale the damage initiation of dual-phase steels under various stress states from microscale to macroscale.Computational materials science,vol.94,November,pp.245–257.

    Liu,H.T.;Han,Z.D.;Rajendran,A.M.;Atluri,S.N.(2006):Computational Modeling of Impact Response with the RG Damage Model and the Meshless Local Petrov-Galerkin(MLPG)Approaches.Computers,Materials&Continua,vol.4,Issue 1,pp.43-54.

    Ma,N.(2013):Accurate simulation on failure and springback of sheet metal forming.Sokeizai,vol.54,no.4,pp.21-26(in Japanese).

    Ma,N.;Takada,K;Sato,K.(2014):Measurement of local strain path and identification of ductile damage limit based on simple tensile test..Procedia Engineering,vol.81,pp.1402-1407.

    Marciniak,Z.;Kuczynski,K.(1967):Limit strains in the processes of stretchforming sheet metal International.Journal of Mechanical Science,vol.9,pp.609–620.

    Pan,B.;Qian,K.;Xie,H.;Asundi,A.(2009):Two-dimensional digital image correlation for in-plane displacement and strain measurement:a review.Measure-ment Science and Technology,vol.20.

    Paul,S.K.(2012):Micromechanics based modeling of Dual Phase steels:Prediction of ductility and failure modes.Computational Materials Science,vol.56,April,pp.34–42.

    Peters,W.H.;Ranson,W.F.(1981):Digital imaging techniques in experimental stress analysis.Optical Engineering,vol.21,pp.27–431.

    Ramazani,A.;Schwedt,A.;Aretz,A.;Prahl,U.;Bleck,W.(2013):Characterization and modelling of failure initiation in DP steel..Computational materials science,vol.75,July,pp.35-44.

    Roux,S.;Hild,F.(2008):Ditigal image mechanical indentification,Experimental Mechanics,Doi:101007/s11340-007-9103-03.

    Sato,J.(2006):Recovering multiplye view geometry from manual projections of multiply cameras.Int.J.Comput.,vol.66,pp.123-140.

    Sato,K.;Yu,Q.;Hiramoto,J.;Urabe,T.;Yoshitake,A.(2015):A method to investigate strain rate effects on necking and fracture behaviors of advanced highstrength steels using digital imaging strain analysis.Int.J.Imp.Eng.,vol.75,pp.11-26.

    Sladek,J.;Stanak,P.;Han,Z-D.;Sladek,V.;Atluri,S.N.(2013):Applications of the MLPG Method in Engineering&Sciences:A Review.Computer Modeling in Engineering&Sciences,vol.92,Issue 5,pp.423-475

    Sutton,M.A.(1986):Application of an optimized digital correlation method to planar deformation analysis,Computer Vision/Computer Graphics Collaboration Techniques.Computer Science,vol.4,pp.143-150.

    Takada,K.;Sato,K.;Ma,N.(2015):Fracture Prediction for Automotive Bodies using a Ductile Fracture Criterion and a Strain-Dependent Anisotropy Model.SAE Int.J.Mater.Manuf.,vol.8,Issue 3.

    Takuda,H.;Hama,T.;Nishida,K.;Yoshida,T.;Nitta,J.(2009):Prediction of forming limit in stretch flanging by finite element simulation combined with ductile fracture criterion.Computer Methods in Materials Science,vol.9,no.1,pp.137-142.

    Tarigopula,V.;Hopperstad,O.S.;Langseth,M.;Clusen,A.H.;Hild,F.;Lademo,O.G.;Eriksson,M.(2008):A study of large plastic deformations in dual phase steel using digital image correlation and FE analysis..Experimental Mechanics,vol.48,no.2,pp.181-196.

    Ueda,Y.;Fukuda,K.(1989):New Measuring Method of Three-Dimensional,Residual Stresses in Long Welded Joints Using Inherent Strains as Parameters.Trans.of the ASME,J.Eng.Materials and Technology,vol.111,pp.1-8.

    Yoneyama,S.(2011):Smoothing Measured Displacements and Computing Strains Utilising Finite Element Method.Strain,vol.47,no.2,pp.258–266.

    Zeng,D.;Xia,Z.(2010):Extending Tensile Curves beyond Uniform Elongation Using Digital Image Correlation:Capability Analysis.SAE Int.J.Mater.Manuf.,vol.3,no.1,pp.702-710

    1Osaka University,Osaka,Japan.

    2JSOL Corporation,Osaka,Japan.

    3JFE Steel Corporation,Hiroshima,Japan.

    4Honda R&D Co.Ltd,Tochigi,Japan.

    精品熟女少妇八av免费久了| 亚洲国产精品合色在线| 国产午夜精品久久久久久| 亚洲熟妇熟女久久| 国产精品爽爽va在线观看网站| 国产精品久久久久久人妻精品电影| 久久久精品国产亚洲av高清涩受| 亚洲av片天天在线观看| 人人妻人人澡欧美一区二区| 俄罗斯特黄特色一大片| 夜夜夜夜夜久久久久| 免费一级毛片在线播放高清视频| 国产一区二区三区视频了| 国产免费av片在线观看野外av| 久久精品国产亚洲av高清一级| aaaaa片日本免费| 啦啦啦韩国在线观看视频| 国产亚洲av高清不卡| 日韩精品免费视频一区二区三区| 一个人免费在线观看电影 | av福利片在线| 精品国内亚洲2022精品成人| 国产高清视频在线播放一区| 琪琪午夜伦伦电影理论片6080| 香蕉国产在线看| 欧美色欧美亚洲另类二区| 一进一出好大好爽视频| 十八禁人妻一区二区| 亚洲av电影在线进入| 亚洲七黄色美女视频| 一边摸一边抽搐一进一小说| 欧美zozozo另类| 久久午夜亚洲精品久久| 欧美日韩亚洲国产一区二区在线观看| 中文字幕av在线有码专区| 黄色 视频免费看| 免费在线观看完整版高清| 国产精品久久视频播放| 中文字幕最新亚洲高清| 国产精品免费视频内射| 99热只有精品国产| 久9热在线精品视频| 亚洲在线自拍视频| 亚洲成人免费电影在线观看| 亚洲av美国av| 国产成人aa在线观看| 日日夜夜操网爽| 精品一区二区三区视频在线观看免费| 1024香蕉在线观看| 高清毛片免费观看视频网站| 亚洲天堂国产精品一区在线| 久久婷婷成人综合色麻豆| 久久久久久久久久黄片| 免费在线观看黄色视频的| 人妻丰满熟妇av一区二区三区| 在线观看www视频免费| 色综合亚洲欧美另类图片| 18美女黄网站色大片免费观看| 欧洲精品卡2卡3卡4卡5卡区| 一区二区三区激情视频| 国产精品九九99| 久99久视频精品免费| 啦啦啦韩国在线观看视频| 免费在线观看影片大全网站| 亚洲人成77777在线视频| 亚洲成a人片在线一区二区| 黑人操中国人逼视频| 久久精品91蜜桃| 国内少妇人妻偷人精品xxx网站 | 欧美乱妇无乱码| 亚洲性夜色夜夜综合| 日韩欧美免费精品| 久久精品成人免费网站| 99热6这里只有精品| 成人18禁在线播放| 特级一级黄色大片| 午夜福利欧美成人| 18禁裸乳无遮挡免费网站照片| 国产片内射在线| 精品少妇一区二区三区视频日本电影| 一区福利在线观看| 午夜免费激情av| 观看免费一级毛片| 身体一侧抽搐| 成在线人永久免费视频| www国产在线视频色| 亚洲中文日韩欧美视频| 国产av麻豆久久久久久久| 亚洲欧美日韩无卡精品| av中文乱码字幕在线| 久久精品夜夜夜夜夜久久蜜豆 | 欧美日韩乱码在线| 欧美日韩福利视频一区二区| 精品久久蜜臀av无| 女同久久另类99精品国产91| 日本精品一区二区三区蜜桃| 欧美日韩一级在线毛片| 看黄色毛片网站| 色哟哟哟哟哟哟| 久久精品影院6| 国产亚洲av嫩草精品影院| 亚洲av日韩精品久久久久久密| 国产精品一区二区精品视频观看| 99精品欧美一区二区三区四区| 俺也久久电影网| 国产精品永久免费网站| 免费在线观看影片大全网站| 51午夜福利影视在线观看| 欧美成人免费av一区二区三区| 亚洲午夜精品一区,二区,三区| 欧美黑人精品巨大| 欧美色视频一区免费| 一级毛片女人18水好多| 久久久国产欧美日韩av| 精品久久久久久久末码| 脱女人内裤的视频| 少妇人妻一区二区三区视频| 亚洲aⅴ乱码一区二区在线播放 | 久久这里只有精品中国| 桃红色精品国产亚洲av| 高潮久久久久久久久久久不卡| 天天躁夜夜躁狠狠躁躁| 香蕉丝袜av| 色av中文字幕| 不卡一级毛片| 制服丝袜大香蕉在线| 久久久久久久久免费视频了| 国产精品免费视频内射| 麻豆国产97在线/欧美 | 国产爱豆传媒在线观看 | 国产精品1区2区在线观看.| 午夜免费观看网址| 夜夜夜夜夜久久久久| 成年版毛片免费区| 亚洲av日韩精品久久久久久密| 首页视频小说图片口味搜索| 亚洲av五月六月丁香网| 香蕉av资源在线| 久久午夜亚洲精品久久| 亚洲最大成人中文| 美女高潮喷水抽搐中文字幕| 又爽又黄无遮挡网站| 三级男女做爰猛烈吃奶摸视频| 亚洲精品色激情综合| 757午夜福利合集在线观看| 97人妻精品一区二区三区麻豆| 男女视频在线观看网站免费 | a级毛片a级免费在线| 久久久久久亚洲精品国产蜜桃av| 又紧又爽又黄一区二区| 黄色视频,在线免费观看| 蜜桃久久精品国产亚洲av| 日韩中文字幕欧美一区二区| 露出奶头的视频| 国产v大片淫在线免费观看| 日韩欧美在线乱码| 欧美 亚洲 国产 日韩一| 亚洲精品国产一区二区精华液| 欧美黄色淫秽网站| 国产探花在线观看一区二区| 老司机在亚洲福利影院| 丁香欧美五月| 亚洲av日韩精品久久久久久密| 精品久久蜜臀av无| 亚洲中文字幕一区二区三区有码在线看 | 午夜福利高清视频| 免费人成视频x8x8入口观看| 观看免费一级毛片| 一级a爱片免费观看的视频| 亚洲精华国产精华精| 国产激情欧美一区二区| 亚洲美女视频黄频| 精品熟女少妇八av免费久了| 黄色成人免费大全| 国内毛片毛片毛片毛片毛片| 日本一本二区三区精品| 国产亚洲精品第一综合不卡| 香蕉国产在线看| 一本综合久久免费| 久久久水蜜桃国产精品网| 欧美色欧美亚洲另类二区| 亚洲专区国产一区二区| 床上黄色一级片| 又爽又黄无遮挡网站| 亚洲熟女毛片儿| 久久久国产精品麻豆| 欧美成人免费av一区二区三区| 亚洲精品在线美女| 一个人观看的视频www高清免费观看 | 欧美国产日韩亚洲一区| 婷婷精品国产亚洲av在线| 国产精品久久久久久人妻精品电影| 看免费av毛片| 丰满的人妻完整版| 成年免费大片在线观看| 日本成人三级电影网站| 国产伦人伦偷精品视频| 日日爽夜夜爽网站| 国产一区二区在线观看日韩 | 亚洲av成人av| 老司机午夜福利在线观看视频| 午夜两性在线视频| 亚洲aⅴ乱码一区二区在线播放 | 中国美女看黄片| 欧美高清成人免费视频www| 国产蜜桃级精品一区二区三区| 亚洲一卡2卡3卡4卡5卡精品中文| 亚洲欧美精品综合一区二区三区| 日韩欧美一区二区三区在线观看| 国产精品一区二区三区四区久久| 变态另类丝袜制服| 丰满人妻熟妇乱又伦精品不卡| 1024香蕉在线观看| 99国产精品一区二区蜜桃av| 欧美3d第一页| 无限看片的www在线观看| 免费av毛片视频| 国内精品久久久久久久电影| 久久精品91蜜桃| 欧美日韩乱码在线| 99国产精品99久久久久| 国内揄拍国产精品人妻在线| 久久久久性生活片| 精品不卡国产一区二区三区| 一夜夜www| 黄色a级毛片大全视频| 精品久久久久久,| 给我免费播放毛片高清在线观看| 久久久久国产精品人妻aⅴ院| 麻豆国产97在线/欧美 | 亚洲国产精品sss在线观看| 国产成人影院久久av| 女人爽到高潮嗷嗷叫在线视频| 岛国视频午夜一区免费看| 成年女人毛片免费观看观看9| 妹子高潮喷水视频| 欧美黑人巨大hd| 9191精品国产免费久久| 长腿黑丝高跟| 亚洲欧洲精品一区二区精品久久久| 波多野结衣高清无吗| 久久天躁狠狠躁夜夜2o2o| 麻豆一二三区av精品| 精品国产美女av久久久久小说| 很黄的视频免费| 动漫黄色视频在线观看| 久久天躁狠狠躁夜夜2o2o| 香蕉久久夜色| 琪琪午夜伦伦电影理论片6080| 正在播放国产对白刺激| 国产又色又爽无遮挡免费看| 操出白浆在线播放| 在线观看www视频免费| 国产成人精品久久二区二区91| 精品国产乱子伦一区二区三区| 97人妻精品一区二区三区麻豆| 精品熟女少妇八av免费久了| 不卡一级毛片| 男人的好看免费观看在线视频 | 日韩高清综合在线| 婷婷亚洲欧美| 国产伦人伦偷精品视频| 搡老熟女国产l中国老女人| 欧美黄色片欧美黄色片| 日韩欧美国产在线观看| avwww免费| 老司机午夜十八禁免费视频| 老熟妇仑乱视频hdxx| 色精品久久人妻99蜜桃| 亚洲欧美日韩高清在线视频| 日本免费一区二区三区高清不卡| 999久久久精品免费观看国产| 在线观看www视频免费| 亚洲,欧美精品.| 激情在线观看视频在线高清| 精品国产超薄肉色丝袜足j| 亚洲av成人av| 国产成人av教育| 亚洲中文日韩欧美视频| 一进一出好大好爽视频| 亚洲人与动物交配视频| 日韩欧美 国产精品| 老鸭窝网址在线观看| 成人永久免费在线观看视频| 91成年电影在线观看| 两个人看的免费小视频| 国产精品,欧美在线| 成人三级做爰电影| 国产在线观看jvid| 久久精品国产99精品国产亚洲性色| 日本免费一区二区三区高清不卡| 亚洲色图 男人天堂 中文字幕| 亚洲专区国产一区二区| 亚洲五月天丁香| 一本综合久久免费| 丝袜美腿诱惑在线| 亚洲av片天天在线观看| 国产精品爽爽va在线观看网站| 日本一本二区三区精品| 丰满人妻一区二区三区视频av | 亚洲精品美女久久久久99蜜臀| 色尼玛亚洲综合影院| 亚洲成人精品中文字幕电影| e午夜精品久久久久久久| 老熟妇仑乱视频hdxx| 成人手机av| 在线观看66精品国产| 69av精品久久久久久| 熟妇人妻久久中文字幕3abv| 亚洲九九香蕉| 亚洲国产日韩欧美精品在线观看 | 亚洲成人国产一区在线观看| 中文字幕最新亚洲高清| 1024香蕉在线观看| 色尼玛亚洲综合影院| 欧美在线黄色| 欧美成人免费av一区二区三区| 女人高潮潮喷娇喘18禁视频| 成年免费大片在线观看| 午夜福利18| 很黄的视频免费| 日韩免费av在线播放| 又大又爽又粗| 真人做人爱边吃奶动态| 欧美三级亚洲精品| avwww免费| 久久久久亚洲av毛片大全| 熟女少妇亚洲综合色aaa.| 国产视频一区二区在线看| 国产黄a三级三级三级人| 精品久久蜜臀av无| 搡老岳熟女国产| 在线看三级毛片| 伊人久久大香线蕉亚洲五| 亚洲一卡2卡3卡4卡5卡精品中文| 国产精品综合久久久久久久免费| 麻豆国产av国片精品| 亚洲中文字幕日韩| 日本一二三区视频观看| 久久久久久免费高清国产稀缺| 丝袜人妻中文字幕| 久久精品91蜜桃| 嫁个100分男人电影在线观看| 亚洲专区国产一区二区| 高清毛片免费观看视频网站| 精品国产亚洲在线| 麻豆国产av国片精品| 欧美日韩黄片免| 俺也久久电影网| 小说图片视频综合网站| 国模一区二区三区四区视频 | 日本 av在线| 成人国语在线视频| 欧美在线黄色| 国产精品 欧美亚洲| 在线观看一区二区三区| 国产精品永久免费网站| 美女黄网站色视频| 黑人欧美特级aaaaaa片| 亚洲精品美女久久av网站| 啦啦啦观看免费观看视频高清| 亚洲国产欧洲综合997久久,| 国产探花在线观看一区二区| 亚洲av日韩精品久久久久久密| 看片在线看免费视频| 国产精品九九99| 久久 成人 亚洲| 久久久水蜜桃国产精品网| 亚洲 欧美一区二区三区| 长腿黑丝高跟| 精品久久久久久久人妻蜜臀av| 亚洲精品粉嫩美女一区| 亚洲国产欧美网| 最近视频中文字幕2019在线8| 啦啦啦观看免费观看视频高清| 久久精品影院6| 成人欧美大片| 国产精品自产拍在线观看55亚洲| 欧美绝顶高潮抽搐喷水| 欧美在线一区亚洲| 一边摸一边做爽爽视频免费| 亚洲色图 男人天堂 中文字幕| 麻豆久久精品国产亚洲av| 午夜福利在线观看吧| 黄色毛片三级朝国网站| 国产精品久久久av美女十八| 香蕉丝袜av| 熟女少妇亚洲综合色aaa.| 日韩精品青青久久久久久| 国产精品九九99| 久久九九热精品免费| 十八禁人妻一区二区| 夜夜躁狠狠躁天天躁| 国产成人系列免费观看| 操出白浆在线播放| 婷婷亚洲欧美| 免费高清视频大片| 国产成人欧美在线观看| 亚洲专区国产一区二区| 国产精品精品国产色婷婷| xxx96com| 美女大奶头视频| 日韩有码中文字幕| 国产亚洲av高清不卡| 色播亚洲综合网| 日日夜夜操网爽| 欧美午夜高清在线| 午夜精品一区二区三区免费看| 九色成人免费人妻av| 国产av麻豆久久久久久久| 国产精品爽爽va在线观看网站| 亚洲国产日韩欧美精品在线观看 | 亚洲自拍偷在线| 久久精品国产99精品国产亚洲性色| 国产真人三级小视频在线观看| 法律面前人人平等表现在哪些方面| 久久天堂一区二区三区四区| 成人一区二区视频在线观看| 国产精品自产拍在线观看55亚洲| 欧美激情久久久久久爽电影| 嫁个100分男人电影在线观看| 久久婷婷成人综合色麻豆| bbb黄色大片| 国产野战对白在线观看| 人妻久久中文字幕网| 日韩欧美一区二区三区在线观看| 免费观看精品视频网站| 99国产综合亚洲精品| 欧美绝顶高潮抽搐喷水| 18美女黄网站色大片免费观看| 国产午夜福利久久久久久| 在线观看免费午夜福利视频| 日韩成人在线观看一区二区三区| 国内久久婷婷六月综合欲色啪| 久久久久久久久久黄片| 精品高清国产在线一区| 国产高清有码在线观看视频 | а√天堂www在线а√下载| 国产精品1区2区在线观看.| 高清毛片免费观看视频网站| 男人舔女人的私密视频| 欧美一级毛片孕妇| 成人国语在线视频| 国产精品爽爽va在线观看网站| 免费观看精品视频网站| 露出奶头的视频| 不卡av一区二区三区| a级毛片在线看网站| 一区二区三区激情视频| 日韩有码中文字幕| 美女免费视频网站| 最新美女视频免费是黄的| 久久久水蜜桃国产精品网| 亚洲成av人片在线播放无| 日本a在线网址| 在线观看66精品国产| 国产av又大| av天堂在线播放| 婷婷丁香在线五月| 岛国视频午夜一区免费看| 国产av麻豆久久久久久久| 国内毛片毛片毛片毛片毛片| 一二三四在线观看免费中文在| 青草久久国产| 人人妻,人人澡人人爽秒播| 亚洲第一欧美日韩一区二区三区| 国产精品,欧美在线| 国产1区2区3区精品| 丝袜美腿诱惑在线| 日本成人三级电影网站| bbb黄色大片| 搡老岳熟女国产| 国产人伦9x9x在线观看| 欧美性猛交黑人性爽| 国产免费男女视频| 嫁个100分男人电影在线观看| 男人舔女人的私密视频| 成人亚洲精品av一区二区| 亚洲激情在线av| 免费在线观看日本一区| 免费看美女性在线毛片视频| av免费在线观看网站| 亚洲欧美日韩高清专用| 欧美日本亚洲视频在线播放| 国产久久久一区二区三区| 免费高清视频大片| 免费在线观看视频国产中文字幕亚洲| 欧美成人一区二区免费高清观看 | 亚洲成人精品中文字幕电影| 国产高清videossex| 国产av不卡久久| 国产熟女xx| 午夜福利在线观看吧| 日韩免费av在线播放| 亚洲欧美精品综合久久99| 国产av一区二区精品久久| 啦啦啦观看免费观看视频高清| 香蕉丝袜av| 欧美成人午夜精品| aaaaa片日本免费| 久久精品国产亚洲av香蕉五月| 99精品在免费线老司机午夜| 18禁美女被吸乳视频| 精品国产亚洲在线| 首页视频小说图片口味搜索| 宅男免费午夜| 亚洲国产日韩欧美精品在线观看 | 亚洲成av人片在线播放无| 国内揄拍国产精品人妻在线| 精品乱码久久久久久99久播| 午夜免费观看网址| 欧美另类亚洲清纯唯美| av欧美777| 欧洲精品卡2卡3卡4卡5卡区| 国产精品av视频在线免费观看| 一级作爱视频免费观看| 精品久久久久久,| 观看免费一级毛片| 舔av片在线| 久久国产精品影院| 久久久久性生活片| 18禁国产床啪视频网站| 国产熟女xx| 一本久久中文字幕| 午夜影院日韩av| 观看免费一级毛片| 91av网站免费观看| 国产99白浆流出| 欧美中文日本在线观看视频| 免费电影在线观看免费观看| 欧美乱码精品一区二区三区| 又大又爽又粗| 欧美日本亚洲视频在线播放| 黑人欧美特级aaaaaa片| 国产黄色小视频在线观看| 国产99白浆流出| 可以免费在线观看a视频的电影网站| 美女高潮喷水抽搐中文字幕| 日韩国内少妇激情av| 好男人在线观看高清免费视频| 亚洲成人中文字幕在线播放| 成人一区二区视频在线观看| 99久久无色码亚洲精品果冻| 欧美一区二区精品小视频在线| 99热这里只有精品一区 | 免费观看人在逋| 亚洲中文字幕日韩| 狂野欧美激情性xxxx| 中文字幕人妻丝袜一区二区| 超碰成人久久| 国产av又大| 色噜噜av男人的天堂激情| 丁香六月欧美| 母亲3免费完整高清在线观看| 97碰自拍视频| 制服人妻中文乱码| 午夜亚洲福利在线播放| 亚洲成人久久性| 色av中文字幕| 视频区欧美日本亚洲| 麻豆成人av在线观看| 成年免费大片在线观看| 日本精品一区二区三区蜜桃| 又黄又爽又免费观看的视频| 亚洲人与动物交配视频| 香蕉av资源在线| 无限看片的www在线观看| 欧美色视频一区免费| 男人舔女人下体高潮全视频| 精品久久久久久久人妻蜜臀av| 18美女黄网站色大片免费观看| 最好的美女福利视频网| 久久久精品大字幕| 亚洲九九香蕉| 不卡一级毛片| 欧美中文综合在线视频| 午夜福利免费观看在线| 中文亚洲av片在线观看爽| 国产69精品久久久久777片 | 久99久视频精品免费| 国产免费av片在线观看野外av| 久久久久免费精品人妻一区二区| 男人的好看免费观看在线视频 | 精品一区二区三区av网在线观看| 国产成人精品久久二区二区91| 欧美+亚洲+日韩+国产| 欧美大码av| а√天堂www在线а√下载| 精品电影一区二区在线| 欧美中文日本在线观看视频| 婷婷六月久久综合丁香| 最近最新中文字幕大全电影3| 国产蜜桃级精品一区二区三区| 法律面前人人平等表现在哪些方面| 伦理电影免费视频| 麻豆成人av在线观看| 欧美成狂野欧美在线观看| 成人18禁高潮啪啪吃奶动态图| 777久久人妻少妇嫩草av网站| 中亚洲国语对白在线视频| 亚洲人成网站在线播放欧美日韩| 免费看十八禁软件| 亚洲精品中文字幕在线视频| 亚洲成a人片在线一区二区| 一a级毛片在线观看| 91老司机精品| 欧美性猛交黑人性爽| 欧美乱色亚洲激情| 男女视频在线观看网站免费 | 亚洲无线在线观看| 日韩国内少妇激情av| 久久久久国内视频| 好男人电影高清在线观看| 婷婷丁香在线五月| 欧美黑人欧美精品刺激|