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

    Stress analysis and damage evolution in individual plies of notched composite laminates subjected to in-plane loads

    2017-11-21 12:54:43HuJunshnZhngKifuChengHuiLiuPingZouPengSongDnlong
    CHINESE JOURNAL OF AERONAUTICS 2017年1期

    Hu Junshn,Zhng Kifu,*,Cheng Hui,Liu Ping,Zou Peng,Song Dnlong

    aThe Ministry of Education Key Lab of Contemporary Design and Integrated Manufacturing Technology,Northwestern Polytechnical University,Xi’an 710072,China

    bSchool of Mechanical Engineering,Northwestern Polytechnical University,Xi’an 710072,China

    Stress analysis and damage evolution in individual plies of notched composite laminates subjected to in-plane loads

    Hu Junshana,Zhang Kaifua,*,Cheng Huib,Liu Pinga,Zou Penga,Song Danlonga

    aThe Ministry of Education Key Lab of Contemporary Design and Integrated Manufacturing Technology,Northwestern Polytechnical University,Xi’an 710072,China

    bSchool of Mechanical Engineering,Northwestern Polytechnical University,Xi’an 710072,China

    Damage evolution;Failure mode;Failure strength;Layer-by-layer stresses;User-defined subroutine

    This work aims to investigate local stress distribution,damage evolution and failure of notched composite laminates under in-plane loads.An analytic method containing uniformed boundary equations using a complex variable approach is developed to present layer-by-layer stresses around the notch.The uniformed boundary equations established in series together with conformal mapping functions are capable of dealing with irregular boundary issues around the notch and at infinity.Stress results are employed to evaluate the damage initiation and propagation of notched composites by progressive damage analysis(PDA).A user-de fined subroutine is developed in the finite element(FE)model based on coupling theories for mixed failure criteria and damage mechanics to ef ficiently investigate damage evolution as well as failure modes.Carbon/epoxy laminates with a stacking sequence of[45°/0°/-60°/90°]sare used to investigate surface strains,in-plane load capacity and microstructure of failure zones to provide analytic and FE methods with strong validation.Good agreement is observed between the analytic method,the FE model and experiments in terms of the stress(strain)distributions,damage evaluation and ultimate strength,and the layerby-layer stress components vary according to a combination effect of fiber orientation and loading type,causing diverse failure modes in individuals.

    1.Introduction

    Carbon fiber reinforced polymer(CFRP)composites are widely used in all fields of aerospace,automobile,electronic power,and mechanical engineering for several advantages offered over metals,ceramics,and plastics.These include low density,robust specific strength,low thermal expansion,corrosion resistance and designable characteristics for lightweight,efficient structures.Cutouts,especially circular and elliptical shapes,commonly appear in composite structural components due to requirements for stability,maneuverability,low weight optimization and accessibility of other systems.This is often the case in aircraft where composite structures such as wing spars and thin walls are drilled for electronic wires and hydraulic pipes1or to facilitate assembly operations.2Other instances occur not as part of the initial design of essential structures but due to material defects or unexpected damage during a service cycle.3These are serious and dangerous cases since they are not supposed to be there.Whatever their reasons for being,the integrity and continuity of fiber and matrix in the composite are destroyed;this makes the notch region the weakest part of the structure and causes serious local stress concentrations in the vicinity.As a consequence,structural capacity is reduced and uncertain reactions to external loads may occur and lead to unexpected fractures in service.In order to avoid potential safety hazards,accurate analysis of local stress levels and damage evolution of notched laminates are of great significance to the utilization of the material and lay the foundation for engineering applications.

    The determination of stress fields in notched anisotropic plates has been the focus of many scholars for a long time.Muskhelishvili4first introduced complex potential theory to the isotropic elastic plates and successfully obtained an accurate solution of stress distribution.Analytical solutions for the stress distribution around holes of different shapes in anisotropic plates were given by Lekhnitskii5using series methods.Savin6presented a much simpler approach by conformal mapping of Cauchy integrals.Gao7used a biaxial loading factor together with an arbitrary orientation angle previously used to solve the problem of a plate with biaxial loading at infinity,avoiding the superposition of a solution with two uniaxial loading problems.The stress analysis research above established the base for later study of notched composite laminates.Ukadgaonker et al.developed a general solution for stress around oval holes,8triangular holes9with rounded corners10and irregular shaped holes11in cross-ply and angleply orthotropic plates under in-plane loading by superposition of two stage solutions for boundary conditions.The stresses were employed together with a series of failure criteria,which were an extension of the Von Mises criterion for quadratic interaction,to calculate the first ply failure(FPF)strength.Similar solutions around a rectangular hole in an infinite isotropic and anisotropic plate were given by Pan et al.12and Rao et al.,13respectively,using the complex variable method.Sharma14suggested general solutions for determining the stress distribution around polygonal holes and investigated the effect of hole geometry and loading pattern on the stress concentration factor.Batista15used a modified solution to solve problems of stress distribution around polygonal holes of complex geometry in an infinite plate subjected to uniform loading at infinity.Remeepazhand and Jafari16also studied the central polygonal hole problem in composite plates using two simple equations with parameters λ,c,nandwcontrolling the size,shape and bluntness of corners.Toubal et al.17experimentally investigated the tensile strain field of composite plates in the presence of stress concentration caused by geometrical cutouts consisting of circular holes by Electronic Speckle Pattern Interferometer(ESPI).The stress obtained in experiments is consistently lower than the analytical and numerical models.

    This literature mostly concentrated on stress distribution in anisotropic plates as well as laminated composites,which are equivalent anisotropic plates,thus the stresses around cutouts were simplified to be uniform throughout,regardless of thickness.It should be noted that composite laminates contain several plies which possess dissimilar properties due to different fiber orientations,so the stress distribution in individual plies should be thoroughly investigated.

    The above stress research methodologies and conclusions establish a theoretical basis for damage and failure of anisotropic plates.As for failure analysis,the damage propagation process is often depicted based on progressive damage analysis(PDA)utilizing numerical and experimental methods.18Lapczyk and Hurtado19proposed an anisotropic damage model suitable for predicating failure and post-failure behavior in fiber-reinforced materials.The plane stress formulation is used and the response of the undamaged material is assumed to be linearly elastic.The evaluation law is based on fracture energy dissipation and implemented in a finite element code.Zahari and EI-Zafrany20developed a progressive analysis algorithm based on Tsai-Hill failure to model the non-linear material behavior and capture the compressive response of woven glass/epoxy composite plates via non-linear finite element analysis.Rakesh et al.21introduced a generic finite element model to investigate the failure of unidirectional glass fiber reinforced plastic(UD-GFRP)composite laminates with drilled holes under tensile testing and compared the results with experimental work done earlier.Effects of joint geometry and stacking sequence on the bearing strength and damage mode were investigated by Ondurucu et al.22and specimens were examined for failure modes using a scanning electron microscope(SEM).Kim et al.23employed rosette strain gauges to measure the strain around the joint hole during insertion of stainless steel pins into glass fiber reinforced plastic(GFRP)specimens,and results were compared with the finite element method.Satapathy et al.24presented a modified fiber failure fatigue model to characterize the behavior of laminated composites with a central circular hole under in-plane fatigue loading.Martins et al.25studied the influence of diameter and thickness on failure pressure during tube burst tests and employed the progressive failure analysis using a damage model by a user subroutine(UMAT)implemented in ABAQUS software to understand the behavior of composite tubes under internal pressure.Lee et al.26proposed an evaluation method for the progressive failure of composite laminates built upon Puck failure criterion by implanted UMAT to efficiently analyze the progressive failure phenomenon in glass/carbon fiber-reinforced composite laminates.Compression tests were performed by Aljibori et al.27on 16 fiber-glass laminated plates with and without circular cut-outs to investigate the effects of varying the centrally located circular cut-out sizes and fiber angle-ply orientations on the ultimate load.Similar experiments were also conducted by Abu et al.28to investigate the influence of cut-outs on multilayer Kevlar-29/epoxy composite laminated plates.

    These remarkable works contribute a lot to the failure analysis of notched anisotropic composites whose layups are quasiisotropic,cross-ply,and angle-ply.At the same time,numerical and advanced experimental methods for progressive failure have greatly developed.However,the failure criteria employed in the literature were linear expressions;whereas composite laminatesin specialapplication are notquasi-isotropic lay-up,composite materials accumulate damage in shear,leading to nonlinear stress-strain relations.The relationship between stress distribution and damage evolution thus needs further investigation.

    In the present study,an analytical model for layer-by-layer stress around notches in individual plies of general composite laminates is investigated using a complex variable approach that contains uniformed boundary functions.The functions established in series together with conformal mapping functions are capable of dealing with irregular boundary issues both around the notch and at infinity.A stress solution is employed to evaluate the progressive damage and failure of notched laminates based on coupling theories for mixed failure criteria and damage mechanics.A user-defined subroutine is then developed in ABAQUS to simulate the damage evolution under tensile and compressive loading.Experiments for strain and failure are conducted using a 3D digital speckle measurement system and microscopy study to validate the proposed analysis method.

    2.Analysis for local stress distribution

    2.1.Problem configuration of notched laminates

    A composite laminate with a traction free notch under inplane loading is investigated in Fig.1.The inplane loading is in the form of stress resultants in an arbitrary coordinate systemx′oy′rotated clockwise by an angle β from the global coordinate systemxoy,as shown in Fig.1(a).Dis the diameter of the circle that encloses the notch.It can be written in terms of Gao’s loading conditions7using stress loading σ∞x′,y′as follows:

    The notched plate is projected to a unit circle on ζ plane under normalized loading in Fig.1(b)with an angle θ to denote hole position.The equivalent loading σ∞x,yalongxoyaxes corresponding to original loading alongx′oy′axes is obtained by transformation of axis:

    By substituting Eqs.(3)and(4)into Eq.(7),a forth order differential equation is obtained:

    Fig.1 Problem configuration and solution of a notched composite under in-plane loading.

    The general solution of bi-harmonic Eq.(6)is determined by roots of the following characteristic equation:

    In general,it can be proved that Eq.(7)has four distinct roots,which are always complex conjugates in pairs:

    where α1,α2,β1,β2are real numbers and β1> 0,β2> 0.μ1,μ2and their conjugates are complex parameters of a plane stress state.They are determined by effective compliance coefficients aijand can be used to describe the degree of material anisotropy.

    Through the integral operation of Eq.(6),the stress functionU(x,y)can be expressed by the real part of two complex functionsF1(z1)andF2(z2):

    wherezp=x+ μpy(p=1,2)are anisotropic complex coordinates,andShrepresents the boundary conditions of the laminate in Fig.1(a).So,the problem is reduced to finding two complex analytical stress functions that satisfy the boundary conditions at the notch edge as well as infinity.

    2.2.Determination of analytic stress functions

    Here,in the issue of plane stress state,the uniformed boundary equations φ(z1)and ψ(z2)are developed in the form of derivatives ofF1andF2:

    where parameters γ and γ′are complex variables of the notch boundary.ParametersB,B′andC′are loading condition contents.Parametersanandbnare series coefficients determined by the boundary conditionsSharound the cutouts and at infinity.

    2.2.1.Remote boundary conditions

    At the infinite plate edge(|z|→∞),the mean stress components reach the applied stress loading of the normalized boundary condition:

    Through substituting Eqs.(9)and(11)into Eq.(4),the parametersB,B′andC′are obtained by the following equation group:

    2.2.2.Boundary conditions at the notch edge

    At the notch edge,which is traction free,the notch area on thezplane is projected onto the exterior of a unit circle on ζ complex plane using the conformal mapping method(Fig.1).The unit circle is defined by the Euler formula ζ=eiθin Fig.1(b).This opens up the possibility of a simple general approach to determine the stress functions of the actual notch contour independently.

    The mapping function for the elliptic notch in anisotropic cases fromzplane to ζ plane is defined by:

    whereR=(a+b)/2 is a constant for size of the hole,andm1=(a-b)/(a+b)is defined to control the ellipticity.Parametersaandbare major and minor axes of the elliptic hole,respectively.These are typically:(1)m1=0 ora=b:circular holes;(2)0<m1<1 ora≠b:elliptical holes;(3)m1→ 1 ora?b:blunted cracks.

    By the projection process and taking Eq.(13)into the uniformed boundary Eq.(10)on ζ plane,the boundary equations in ζ plane are written as:

    It is seen that the boundary of the notch is free of traction,which indicates that the mean resultant forces on the notch are zero:

    Parameters γ and γ′are confirmed to be zero,and the logarithmic terms drop out.The rest of the parameters in series terms are obtained by substituting Eqs.(14)and(16)into Eq.(15):

    whereK1,K2,K3andK4are complex constants.They are determined by:

    2.3.Local stresses around notches of individual plies

    By substituting the analytic stress functions into Eq.(4),the mean stress components around the hole in global coordinates(Cartesian coordinates in Fig.1)are obtained through the following equations:

    According to the classical laminate theory29,30,laminate strains can be written in terms of the midplane strains ε0x,yand midplane curvatures κx,yas follows:

    where εx, εyand γxyarex-direction,y-direction and shear strain components inzthickness layer,and ε0x, ε0yand γ0xyarex-direction,y-direction and shear strain components in the midplane.κx,κyand κxyarex-direction,y-direction and torsional curvatures in midplane,respectively.

    For the symmetric laminates studied here,it can be proven that the coupling stiffness coefficients Bij(i,j=1,2,6)are zero,which indicates the resultant and moment terms are uncoupled.Hence the curvature of midplane κx,yis equal to zero.So,the strain at an arbitrary point of the notch edge throughout its thickness is the same as midplane strain and the mean strain:

    The stress components in thekth ply in principal direction σ(k)are calculated by the following equation modi fied from

    L,TStress-Strain Equation and the stress transformation of axis:

    3.Analysis for damage and failure

    As mentioned before,one of the main purposes of determining local stress distribution in individual plies around the notch is to investigate the damage and failure of the perforated laminate accurately.In progressive damage analysis(PDA)theory,31composites are materials characterized by complicated failure phenomena that interact.The onset of failure in an individual ply generally does not lead to the structure’s collapse,and this may still not occur until the failure has spread to each of the multiple plies,indicating that some plies fail before others.Since the notch region is the weakest part of a composite plate,the damage initiates around the notch.

    3.1.Failure criterion

    Most failure criteria are extended from the von Mises criterion for quadratic interaction.32They allow the failure of a whole ply to be checked,but ignore the distinct failure modes of composites.Given this,the mixed failure criteria,which combines Hanshin criteria and maximum criteria,differs greatly from them.It not only allows distinguishing between failure in tension compression as well as fiber and matrix,but also takes nonlinear shear stress-strain behavior into consideration.The mixed failure criteria33are given by:

    whereemt,emc,eft,efcandefsare failure indices that denote matrix tensile,matrix compression,fiber tensile,fiber compression and the fiber-matrix shear-out,respectively.The parameter α is a material constant that is determined experimentally according to Chang and Lessard34and commonly takes the value of 2.44 × 10-8MPa-3in the User’s Manual of Abaqus 6.11.Once an index exceeds 1.0,the corresponding failure occurs.The termsXt,Xcare tensile and compressive strengths in the fiber direction andYt,Ycare tensile and compressive strengths in the traverse direction.Scis the in-plane shear strength.

    3.2.Material property degradation rules

    From the phenomenological viewpoint,the defect initiation/propagation phenomenon for a material is strongly related to the degradation of the material’s capacity,and the most visible aspect is the degradation of the material’s stiffness.Therefore,it is possible to describe the progressive failure of a composite laminate by numerically representing the correlation between an increase in the material’s internal damage and a decrease in the material’s stiffness.The Camanho and Matthews degradation rules35(Table 1)are adopted because they are closest to the real conditions according to the research of Chang.34ParametersE1,E2,G12and υ12are the longitudinal Young’s modulus,transverse Young’s modulus,shear modulus and Poison’s ratios of the lamina.E′1,E′2,G′

    12and υ′12are material properties after degradation.

    3.3.Progressive failures of notched laminates

    After using Eq.(23)to obtain the stress distribution around the notch in each individual ply,the mixed failure criteria is employed to check the failure state around the notch edge.On the occurrence of some sort of failure,the failing ply should be degraded according to Table 1,which causes redistribution of stresses around the notch.Thus the failure state of the notch edge will be rechecked.The procedure of stress redistribution and failure recheck will continue until all plies fail or fiber failure occurs,causing sharp decline of carrying capacity.So,the location of damage initiation and final failure strength can be obtained during the analysis.Failure strength here refers to the maximum permissible remote stress loading applied.30

    4.Numerical implementation

    The stress distribution,damage and strength of notched specimens were also investigated using FE analysis implemented in the commercial software ABAQUS 6.11/standard.A composite laminate with a general stacking sequence of[45°/0°/-60°/90°]swas chosen as the investigated object,and the material properties are given in Table 2.The specimen configuration for the open hole tensile(OHT)test is illustrated in Fig.2 according to ASTM D5766M-02a Standard.36The open hole compressive(OHC)test followed ASTM D 6484/6484M-04.37The specimen geometry was similar to that used for the OHC test.The aspect ratio met the conditionW/D≥4 so that the external boundary effect was negligible.38

    In the FE model,the 3D eight-node layered continuum shell element was employed for the laminate.The mesh was structured using linear elements with reduced-integration(SC8R),and the notch region was highly refined by partition to capture the high near-hole stress gradients.A planar mesh size of 0.188 mm×0.275 mm in(inx-yplane)in the vicinity of the hole was used while a coarser mesh was used away from the zone of the interest to computation cycles.In the thickness direction,a single layer of elements was used to model each ply.The approximate element size was 0.206 mm×0.188 mm×0.375 mm.Successive space discretization was compared,varying the element size surrounding the hole to verify that the final mesh size was enough to provide accurate results without expensive computational cost.A total of 92,345 eight-node brick elements were employed to model the specimen.The specimen was encastre on one side and applied tensile or compressive stress loading on the other side.Two lines,between which the length wasL-2S,were marked to define the elongation corresponding to the part of the specimen that was not griped.The mesh and boundary conditions(BCs)are shown in Fig.3.

    The progressive model is implemented in ABAQUS by applying a separate subroutine called user-define field(USDFLD).39In the analysis,the material properties depend on six field variables(FVs).The first five FVs from FV1 to FV5 represent failure indexes corresponding toemt,emc,eft,efcandefs.The last one,FV6,denotes the shear damage parameter.The algorithm is illustrated in Fig.4.In the beginning,the FVs and the solution-dependent state variables(SDVs)are set to 0 and the material properties are equal to their initial values.Then the load is increased gradually;For each load increment,several iterations are necessary before the analysis converges to an equilibrium state.At the end of each increment,stresses,failure indices and shear damage parameter are computed at the integration points.The value of the shear damage parameter is directly assigned to the corresponding FVs and used to degrade the shear modulus.The values of the failure indices are stored as SDVs.Once an SDV exceeds 1.0,the corresponding FV is set to 1 and maintains this value until the end of the analysis.Meanwhile,the material properties are automatically reduced according to the degradation rules in Table 1.With each increasing applied load,the procedure is repeated until the program terminates due to excessive element distortion.Since equilibrium will not be reestablished after degradation of properties,small load increments should be set to obtain an accurate solution.40

    Table1 Property degradation rulesofCamanho and Matthews.

    5.Experimental investigation

    The damage and failure behavior of notched composite laminates is also investigated on a Bairoe electronic universal testing machine with a load capacity of 100 kN.The tests were run in displacement control at a rate of 2 mm/min.Loading stopped after dropping 30%from the peak load.A non-contact measurement method,3D digital speckle measurement system(DSMS),was used to investigate the strain field in the vicinity of the hole in the notched composite plate.The system is based on a 3D digital image correlation(DIC)system with binocular stereo vision technique that together perform fullfield 3D measurements of surface strain and out-of-plane deformation of the joints,as shown in Fig.5(a).The principle is to track an applied surface pattern during loading bycontinuously taking digital images of the surface.41Two charge coupled device(CCD)cameras and two light emitting diode(LED)sources are required for measurement.To ensure data synchronization,DSMS was connected with the testing machine by a synchronous data transmission device.

    Table 2 Material properties of carbon-epoxy T300/3526.

    Fig.2 Specimen configuration of composite for open hole test.

    Fig.3 Mesh and BCs of the FE model.

    Fig.4 Algorithm for ABAQUS user-defined subroutine.

    Fig.5 Test machine with digital speckle measurement system and specimen processing.

    The CFRP laminate was fabricated through hand lay-up method using resin prepregs,and specimens were processed with a water jet cutter to avoid defects caused by cutting heat.The material properties and specimen geometry are the same as the FE model in Section 5.Holes were drilled and reamed with carbide-tipped tools to attain fine surface quality,.then surfaces of the specimens were polished with sandpaper(first using a grain size of 800 then a grain size of 2000)to ensure a smooth surface.All the specimens were cleaned with alcohol to remove grease and dust before the test.A white background was painted on the testing surface of the specimens using spray paint,and later small black spray paint droplets were subtly applied with an appropriate density to make a high contrastpattern,as shown in Fig.5(b).Images of speckle pattern were recorded at a frequency of 5 Hz.The VIC 3D software calculated the accumulated strain by comparing the subsequent images with the first image,which was taken as a reference.

    6.Results and discussion

    6.1.Experimental veri fication of the strain distribution around the notch

    The strain distributions around the notch on the surface of the specimen under tensile and compressive loading are presented in Figs.6 and 7,respectively.The test was repeated three times,and the strains in the vicinity of the open hole were extracted by DSMS(Fig.5(b))to verify the validation of the theory and FE model.Since the strains around the notch are center symmetric,only those whose angular position θ ranges from 0°to 180°are revealed.

    The experimental results show relatively good trending alongside the proposed analytical solution and FE method with acceptable amplitude errors.There are two possible reasons for the errors:Either the open hole is a singularity at the surface of the specimen,or the fragmentary grid on the hole boundary cannot be used to calculate the strain field.42In fact,the strains obtained by DSMS are the same as for the ring belt,which is a small offset from the hole edge.The strains fade much more sharply along the radial direction at the points where the strain concentrations are the most serious at the edge of the hole,so the obtained strains from ring belt near these points are apparently smaller than those from the hole edge.

    The small error between analytical and FE results is due to the algorithm in the FE method,which is a numerical and approximate solution as compared to an analytical solution.Actually,the strains under tension and compression possess the opposite distribution law due to the tensile and compressive deformation state of specimen in either thexoryaxial direction.For strains in thexaxial direction,maximum strain occurs at the 90°point where by pass load goes through net tension plane and causes larger strain than any other part of the specimen.As for strains in theyaxial direction,extreme points occur at regions around 0°and 90°due to transverse contracting under tension(Fig.6(b))or transverse stretching under compression(Fig.7(b)).Extreme points at 45°and 135°are caused by shear effect.

    6.2.Stress distribution around the hole in individual layers

    The layer-by-layer stress distributions under tensile and compressive loading are provided separately in Figs.8 and 9 to give clear insight into the mechanical behavior of notched composite laminates.Good agreement is observed between theory and FEA,thus providing confidence in the accuracy of the present results.Stress components in individual plies are shown to be highly non-uniform throughout the laminate thickness.The normalized longitudinal stresses in an individual layer have a much larger magnitude than the transverse and tangential stresses,reflecting the specific characteristic of anisotropic material in which carbon fibers bear most of the load.

    It is seen that the stress concentration level in a layer depends on the relative angle between the fiber orientation and loading direction.The same dependence exists for concentration regions.For longitudinal stress,the smaller the angle is,the higher concentration level found near the region where the fiber direction is tangential to the hole perimeter.This can explain why the 0°layer exhibit the highest stress peaks of±8.65 near 91.8°and is the most effective in carrying the axial load while the 90°layer has the smallest concentration level of±2.25 at 0°.The transverse stresses in each layer exhibit peaks at locations where the fiber direction deviates most from the loading direction or where the fiber direction is normal to the notch boundary.So the 90°layer shows the largest concentration of±0.62 near 90°while the 0°layer has a level of±0.18 near 0°.Actually,there is adverse distribution rule between longitudinal and transverse stresses.For tangential stress influenced by the shear effect,higher degrees of concentration are found in layers whose fiber orientation possesses an angle in the loading direction.

    Fig.6 Strain distribution around circular notch of composite laminates under tensile loading.

    Fig.7 Strain distribution around circular notch of composite laminates under compressive loading.

    Fig.8 Stress distribution around circular notch in individual layers of composite laminates under tensile loading.

    Fig.9 Stress distribution around circular notch in individual layers of composite laminates under compressive loading.

    6.3.Damage initiation and propagation of notched composite laminates

    Damage initiation and propagation was evaluated by comparing screenshots of various failure modes in different increment steps(IS).Some of the important damage initiation locations were also predicted using the analytical method in section 3.3.A microscopy study that investigated micro-scale failures occurring in the vicinity of the hole using scanning electron microscope(SEM)is shown in Fig.12 to validate the damage mode simulated in the FE model.

    Fig.10 Damage initiation and propagation of notched composite laminates under tensile loading.

    Fig.11 Damage initiation and propagation of notched composite laminates under compressive loading.

    Fig.10 shows the failure mode(i.e.,matrix tensile,matrix compress,fiber tensile,fiber compress,and fiber-matrix shear)in different ply orientations under tensile loading.Matrix tensile failure,which onset first at 92.72°in the 90°layer,was the dominant failure mode in early increment steps.This damage,though possibly not affecting tensile strength much,was caused by transverse stress concentration and spread to other layers with increasing loading.Fiber tensile failure occurred at 88.71°in the 0°layer and was accompanied by matrix compressive failure in later increment steps.As loading increased,fiber damage propagated across the net tension plane causing dramatic decrease in load capacity.Matrix- fiber shear failure was observed due to excessive transverse shrinkage in late loading increments.The failure mode under compressive loading where matrix compressive failure initiation was observed at 88.4°in the 90°layer and propagates in early increment steps is presented in Fig.11.Fiber compressive failure was found along with matrix- fiber shear in the 0°layer before specimen buckling deformation.The initiation location of fiber compressive damage occurred where fiber direction was normal to the hole perimeter.Tensile failure rarely appeared in the process.

    Microscopy investigation of damaged specimen form tests was conducted to back up the progressive failure results from FE model.In the microscopic morphology,matrix crack(Fig.12(a))and fiber breakage(Fig.12(b)),which corresponded to matrix tensile failure and fiber tensile failure from FE results in Fig.10,were observed as common micro-scale failure modes in OHT specimens.Matrix crush(Fig.12(c)),fiber buckling(Fig.12(d))and matrix- fiber shear(Fig.12(e))were dominant failure modes in compressive specimens,just as the FE model simulated in Fig.11.Slight delamination is also found in the experiment because excessive load destroyed the compatibility of deformation among layers.This is beyond the elastic range and not considered in the analytical solution.

    To further support the progressive failure results from the FE model,strain concentration and crack growth path on the specimen surface were captured by DSMS to compare with the damage of the 45°layer(surface layer)in the composites,and the results are listed in Figs.13 and 14.For OHT testing,the initial damage points induced by matrix tensile failure(Fig.13(c))matched very well with the zones of strain concentration(Fig.13(a)),and the damage propagation path follows the crack profile.Fiber tensile failure(fiber tensile rupture)is also observed in both test results(Fig.13(b))and the FE results(Fig.13(d)).As for OHC testing,the initial damage points of failure(Fig.14(c))were in accordance with the zones of strain concentration around the hole(Fig.14(a)).The fibermatrix shear-out failure in the 45°layer(Fig.14(d))resulted in surface cracking of the specimen(Fig.14(b))because they possessed the same growth path.The above analysis provided sufficient credibility to the validation of the FE model for progressive failure analysis.

    Fig.12 Scanning electron microscope images of damage types in the vicinity of the notch.

    Fig.13 Comparison of failure results between experimental test and FE model under tensile loading.

    Fig.14 Comparison of failure results between experimental test and FE model under compressive loading.

    6.4.Failure strength

    The load-displacement plots of experimental and FE simulation results for the notched composite laminates are shown in Figs.15 and 16.The slope of the linearity in the experimental tests(i.e.elastic specimen stiffness)is seen to be very well captured by the numerical models,indicating that the composite can be regarded as linear-elastic material despite its brittleness.The stiffness of the specimen(slope of the curve)from the test is smaller than that from FE results due to the defects of the composite laminates and the degradation rules for the material properties.The OHC test shows slight nonlinear behavior at the start of the curve because the specimen slightly buckled under compressive loading.Since the OHC specimens were not ruptured,compressive failure is much more moderate than in the OHT test.The nonlinear behavior observed near the peak load is associated with severe damage in the specimens leading to fracture.

    Fig.15 Load-displacement curves of notched composite laminates under tensile test.

    Fig.16 Load-displacement curves of notched composite laminates under compressive test.

    In the tensile test,there is a protuberant point at a load of 247 MPa,which is mainly caused by the initiation of fiber tensile damage around the notch in the 0°layer along with serious matrix tensile failure in all plies.With each load increment,the fiber tensile damage spreads across the tension plane at 0°through thickness into 90°,-60°,+45°layers in the vicinity of the hole,resulting in drastic fluctuations in the curve.The experimental ultimate strength was 484.6 MPa with an error of 7.7%compared with the FE result of 524.9 MPa.In the compressive test,the main damage mode in an earlier state was matrix compressive failure.Fiber buckling failure occurs,especially in the 0°layer in the net tension plan,causing a quick shock in the load capacity curve.As the load increased,buckling failure extended beyond the surroundings of the notch and caused bending of the specimen.Since the specimens were not ruptured,the compressive failure is much more moderate than in the tensile test.There was a difference of 5.3%between the test strength of 206.1 MPa and FE result of 217.5 MPa.

    7.Conclusions

    This research investigates the layer-by layer stress components,and the initiation and propagation of damage and failure in notched composite laminates by analytic method and FE model,both of which are validated by experimental tests.There is good agreement between the analytic method and the FE model in terms of the stress(strain)distributions,damage initiation and location of failure modes,and experiments on strain by DSMS and microscopy study by SEM provide them with strong suitability and validation.The proposed analytic fundamentals present a dependable and manageable alternative to evaluate stress level and damage initiation quickly,while the FE model simulates the damage evolution process thoroughly.Several conclusions can be drawn from the research:

    (1)The stress components are all non-uniform and periodic in different plies and caused by the anisotropy of the laminates.The magnitude of longitudinal stress is much larger than that of transverse stress and tangential stress due to the high elastic modulus of the fiber contrasting with the matrix.

    (2)The layer-by-layer stress distributions around cutouts are affected by fiber orientation and loading direction.The high stress level of σLmainly occurs where the fiber direction is tangential to the hole perimeter.The concentration of σTis observed where the fiber direction is normal to the notch boundary.The τLTis more influenced by shear effect and possesses more extreme points than σLand σT.

    (3)Matrix dominant failure is common in early states of inplane loading.It initiates at the hole edge in the ply whose fiber orientation is perpendicular to the load direction and then propagates to other plies.Fiber damage occurs and builds until matrix damage causes ruptures.While fiber buckling occurs in most plies,fiber breakage is found in the ply where the fiber orientation is parallel with loading direction.

    (4)Load-displacement plots exhibit a long linear state,indicating that the composites can be regarded as linearelastic material despite of their brittleness.Matrix damage evolution does not affect tensile strength much while fiber damage initiation may cause sudden shocks in plots.CFRP laminates possess better tensile load capacity than compressive ones.

    Acknowledgements

    This work was sponsored by the National Natural Science Foundation of China,with three different programs(No.51275410,No.51305349 and No.51305352)that supports the present work financially.The authors would like to acknowledge the editors and the anonymous reviewers for their insightful comments.

    1.Ostergaard MG,Ibbotson AR,Roux OL,Prior AM.Virtual testing ofaircraftstructures.CEASAeronautJ2011;1(1–4):83–103.

    2.Singh I,Bhatnagar N.Drilling-induced damage in uni-directional glass fiber reinforced plastic(UD-GFRP)composite laminate.Int J Adv Manuf Tech2006;27(9–10):877–82.

    3.Senthil K,Arockiarajan A,Palaninathan R,Santhosh B,Usha KM.Defects in composite structures:Its effects and prediction methods–acomprehensivereview.ComposStruct2013;106:139–49.

    4.Muskhelishvili NI.Some basic problems of the mathematical theory of elasticity.Noordhoff P.Ltd.,1953,123–129.

    5.Lekhnitskii SG.Anisotropic plates.New York:Gordon and Breach;1968.p.41–5.

    6.Savin GN.Stress concentration around holes.New York:Pergamon Press;1961.p.169–74.

    7.Gao XL.A general solution of an infinite elastic plate with an elliptic hole under biaxial loading.Int J Pres Ves Pip1996;67(1):95–104.

    8.Ukadgaonker VG,Awasare PJ,Gade SV.Stress analysis of door and window of a passenger aircraft.Aeronaut J2000;104(1039):409–14.

    9.Ukadgaonker VG,Rao DKN.Stress distribution around triangularholesin anisotropic plates.ComposStruct1999;45(3):171–83.

    10.Ukadgaonker VG,Rao DKN.A general solution for stresses around holes in symmetric laminates under inplane loading.Compos Struct2000;49(3):339–54.

    11.Ukadgaonker VG,Kakhandki V.Stress analysis for an orthotropic plate with an irregular shaped hole for different in-plane loading conditions– Part 1.Compos Struct2005;70(3):225–74.

    12.Pan ZX,Cheng YS,Liu J.Stress analysis of a finite plate with a rectangular hole subjected to uniaxial tension using modified stress functions.Int J Mech Sci2013;75:265–77.

    13.Rao DKN,Babu MR,Reddy KRN,Sunil D.Stress around square and rectangular cutouts in symmetric laminates.Compos Struct2010;92(12):2845–59.

    14.Sharma DS.Stress distribution around polygonal holes.Int J Mech Sci2012;65(1):115–24.

    15.Batista M.On the stress concentration around a hole in an infinite plate subject to a uniform load at infinity.Int J Mech Sci2011;53(4):254–61.

    16.Remeepazhand J,Jafari M.Stress analysis of perforated composite plates.Compos Struct2005;71(3–4):463–8.

    17.Toubal L,Karama M,Lorrain B.Stress concentration in a circular hole in composite plate.Compos Struct2005;68(1):31–6.

    18.Liu PF,Zheng JY.Recent developments on damage modeling and finite element analysis for composite laminates:a review.Mater Des2010;31(8):3825–34.

    19.Lapczyk I,Hurtado JA.Progressive damage modeling in fiberreinforcedmaterials.ComposPartA–ApplSci2007;38(11):2333–41.

    20.Zahari R,EI-Zafrany A.Progressive failure analysis of composite laminated stiffened plates using the finite strip method.Compos Struct2009;87(1):63–70.

    21.Rakesh PK,Singh I,Kumar D.Failure prediction in glass fiber reinforced plastics laminates with drilled hole under uni-axial loading.Mater Des2010;31(6):3002–7.

    22.Ondurucu A,Esendemir U,Tunay RF.Progressive failure analysis of glass-epoxy laminated composite pinned-joints.Mater Des2012;36:617–25.

    23.Kim SY,He B,Shim CS,Kim D.An experimental and numerical study on the interference-fit pin installation process for cross-ply glass fiber reinforced plastics(GFRP).Compos Part B-Eng2013;54:153–62.

    24.Satapathy MR,Vinayak BG,Jayaprakash K,Naik NK.Fatigue behavior of laminated composites with a circular hole under inplane multiaxial loading.Mater Des2013;51:347–56.

    25.Martins LAL,Bastian FL,Netto TA.Reviewing some design issues for filamentwound composite tubes.MaterDes2014;55:242–9.

    26.Lee CS,Kim JH,Kim SK,Ryu DM,Lee JM.Initial and progressive failure analyses for composite laminates using Puck failure criterion and damage-coupled finite element method.Compos Struct2015;121:406–19.

    27.Aljibori HSS,Chong WP,Mahlia TMI,Chong WT,Edi P,Alqrimli H,et al.Load-displacement behavior of glass fiber/epoxy composite plates with circular cut-outs subjected to compressive load.Mater Des2010;31(1):466–74.

    28.Talib AR Abu,Ramadhan AA,Rafie ASM,Zahari R.Influence of cut-out hole on multi-layer Kevlar-29/epoxy composite laminated plates.Mater Des2013;43:89–98.

    29.Jones RM.Mechanics of composite materials.2nd ed.Philadelphia:Taylor&Francis;1999.p.109–21.

    30.Muc A,Rys J,Latas W.Limit load-carrying capacity for spherical laminated shells under external-pressure.Compos Struct1993;25(1–4):295–303.

    31.Zhang JY,Liu FR,Zhao LB,Chen YL,Fei BJ.A progressive damage analysis based characteristic length method for multi-bolt composite joints.Compos Struct2014;108:915–23.

    32.Pietropaoli E.Progressive failure analysis of composite structures using a constitutive material model(USERMAT)developed and implemented in ANSYS?.ApplCompos Mater2012;19(3–4):657–68.

    33.Ribeiro ML,Vandepitte D,Tita V.Damage model and progressive failure analyses for filament wound composite laminates.Appl Compos Mater2013;20(5):975–92.

    34.Chang FK,Lessard LB.Damage tolerance of laminated composites containing an open hole and subjected to compressive loadings.1.Analysis.J Compos Mater1991;25(1):2–43.

    35.Camanho PP,Matthews FL.A progressive damage model for mechanically fastened joints in composite laminates.J Compos Mater1999;33(24):2248–80.

    36.ASTM International.D5766/D5766M-02a Standard test method for open hole tensile strength of polymer–matrix composite laminates.2007.

    37.ASTM International.D6484/D6484M-04 Standard test method for open-hole compressive strength of polymer–matrix composite laminates.2004.

    38.Wu HC,Mu B.On stress concentrations for isotropic/orthotropic plates and cylinders with a circular hole.Compos Part B-Eng2003;34(2):127–34.

    39.Hibbit,Karlssonand Sorensen Inc.ABAQUS/Standard version 6.11,user’s manual.2011.

    40.Huang ZM.Inelastic and failure analysis of laminate structures by ABAQUS incorporated with a general constitutive relationship.J Reinf Plast Comp2007;26(11):1135–81.

    41.Ekh J,Schon J,Melin LG.Secondary bending in multi fastener,composite-to-aluminium single shear lap joints.Compos Part BEng2005;36(3):195–208.

    42.Vieille B,Aucher J,Taleb L.Overstress accommodation in notched woven-ply thermoplastic laminates at high-temperature:numerical modeling and validation by digital image correlation.Compos Part B-Eng2013;45(1):290–302.

    23 February 2016;revised 5 April 2016;accepted 5 August 2016

    Available online 21 December 2016

    ?2016 Chinese Society of Aeronautics and Astronautics.Production and hosting by Elsevier Ltd.This is anopenaccessarticleundertheCCBY-NC-NDlicense(http://creativecommons.org/licenses/by-nc-nd/4.0/).

    *Corresponding author.

    E-mail address:zhangkf@nwpu.edu.cn(K.Zhang).

    Peer review under responsibility of Editorial Committee of CJA.

    秋霞伦理黄片| kizo精华| 免费观看av网站的网址| 91国产中文字幕| 五月玫瑰六月丁香| 免费观看无遮挡的男女| 亚洲成人一二三区av| 欧美精品高潮呻吟av久久| 国产精品99久久久久久久久| 最新的欧美精品一区二区| 18禁裸乳无遮挡动漫免费视频| 免费大片18禁| 三级国产精品片| 免费人成在线观看视频色| 国产精品人妻久久久久久| 欧美三级亚洲精品| 一本色道久久久久久精品综合| 亚洲,一卡二卡三卡| a 毛片基地| 欧美 日韩 精品 国产| av视频免费观看在线观看| 色网站视频免费| 18在线观看网站| 欧美日韩精品成人综合77777| 国产午夜精品久久久久久一区二区三区| 91午夜精品亚洲一区二区三区| 成人毛片60女人毛片免费| 国产男女超爽视频在线观看| 免费观看在线日韩| 久久免费观看电影| 街头女战士在线观看网站| 18在线观看网站| av在线观看视频网站免费| 天天躁夜夜躁狠狠久久av| 午夜福利网站1000一区二区三区| 欧美日韩成人在线一区二区| 亚洲色图 男人天堂 中文字幕 | 亚洲五月色婷婷综合| 精品一区二区三区视频在线| 午夜激情福利司机影院| 久久精品国产亚洲网站| 久热久热在线精品观看| 久久这里有精品视频免费| 在线观看人妻少妇| 日韩一区二区视频免费看| 欧美一级a爱片免费观看看| 精品酒店卫生间| 天天躁夜夜躁狠狠久久av| 大码成人一级视频| 老司机亚洲免费影院| 久久精品国产a三级三级三级| 欧美老熟妇乱子伦牲交| 国产黄频视频在线观看| 男女国产视频网站| 一级黄片播放器| 久久狼人影院| 久久久久久久久久久久大奶| 亚洲欧美成人综合另类久久久| 欧美性感艳星| 热re99久久精品国产66热6| 午夜免费观看性视频| 精品熟女少妇av免费看| 国产男人的电影天堂91| 老司机影院成人| av在线app专区| 伦理电影大哥的女人| 亚洲熟女精品中文字幕| 91久久精品国产一区二区成人| 久久久精品区二区三区| 免费大片黄手机在线观看| 久久99蜜桃精品久久| av一本久久久久| 欧美日韩av久久| 午夜影院在线不卡| 三级国产精品欧美在线观看| 日韩av免费高清视频| 春色校园在线视频观看| 国产免费又黄又爽又色| 久久热精品热| 亚洲欧洲精品一区二区精品久久久 | a 毛片基地| 欧美一级a爱片免费观看看| 夜夜看夜夜爽夜夜摸| 精品久久久久久电影网| 99久久精品国产国产毛片| 丝袜喷水一区| 男的添女的下面高潮视频| 久久人人爽人人爽人人片va| 男人爽女人下面视频在线观看| 天天影视国产精品| 国产永久视频网站| 国产男女内射视频| 久久人人爽人人片av| 国产又色又爽无遮挡免| av播播在线观看一区| 在线观看国产h片| 建设人人有责人人尽责人人享有的| 免费看不卡的av| 男女啪啪激烈高潮av片| 欧美日韩视频高清一区二区三区二| 建设人人有责人人尽责人人享有的| 国产精品不卡视频一区二区| 欧美日韩精品成人综合77777| 啦啦啦中文免费视频观看日本| 亚洲精品中文字幕在线视频| 少妇高潮的动态图| 成人国产麻豆网| 亚洲高清免费不卡视频| 午夜福利视频在线观看免费| 精品人妻偷拍中文字幕| 欧美另类一区| 一边摸一边做爽爽视频免费| 国产精品麻豆人妻色哟哟久久| 97超碰精品成人国产| 91精品一卡2卡3卡4卡| 色94色欧美一区二区| 夫妻性生交免费视频一级片| 寂寞人妻少妇视频99o| 91国产中文字幕| 中文天堂在线官网| av.在线天堂| 亚洲国产日韩一区二区| 考比视频在线观看| 精品亚洲成a人片在线观看| 日本黄色片子视频| 丝袜脚勾引网站| 亚洲成人手机| 日韩成人av中文字幕在线观看| 搡老乐熟女国产| 啦啦啦中文免费视频观看日本| av卡一久久| 视频中文字幕在线观看| 国产精品熟女久久久久浪| 建设人人有责人人尽责人人享有的| 乱码一卡2卡4卡精品| 国产伦理片在线播放av一区| 亚洲精品乱码久久久久久按摩| 免费播放大片免费观看视频在线观看| 国产成人精品婷婷| 久久人人爽人人片av| 不卡视频在线观看欧美| 一级毛片 在线播放| 中文乱码字字幕精品一区二区三区| 人妻 亚洲 视频| 午夜日本视频在线| 精品卡一卡二卡四卡免费| 久久久国产欧美日韩av| 成人手机av| 国产免费视频播放在线视频| 大又大粗又爽又黄少妇毛片口| 一级,二级,三级黄色视频| 精品少妇久久久久久888优播| 99热这里只有精品一区| 亚洲精品aⅴ在线观看| 国产视频首页在线观看| 亚洲图色成人| 欧美三级亚洲精品| 激情五月婷婷亚洲| 亚洲成人一二三区av| 老司机影院毛片| 亚洲精品乱久久久久久| 狂野欧美白嫩少妇大欣赏| 日日摸夜夜添夜夜添av毛片| 美女主播在线视频| 亚洲国产最新在线播放| 国产淫语在线视频| 91精品三级在线观看| 这个男人来自地球电影免费观看 | 国产精品一区二区在线观看99| 国产精品国产三级国产av玫瑰| av女优亚洲男人天堂| a级片在线免费高清观看视频| 看非洲黑人一级黄片| 插阴视频在线观看视频| 另类亚洲欧美激情| 一区在线观看完整版| 精品国产国语对白av| 色婷婷久久久亚洲欧美| 最黄视频免费看| 国产成人精品婷婷| 亚洲成人手机| 午夜精品国产一区二区电影| 久久鲁丝午夜福利片| 免费播放大片免费观看视频在线观看| 免费观看性生交大片5| 99热这里只有精品一区| 国产男女超爽视频在线观看| av又黄又爽大尺度在线免费看| 国产亚洲精品久久久com| 飞空精品影院首页| 精品一区二区三卡| 大陆偷拍与自拍| 一区二区av电影网| 久久人人爽av亚洲精品天堂| 久久久久网色| 在线观看人妻少妇| 大陆偷拍与自拍| av福利片在线| 久久综合国产亚洲精品| 欧美精品亚洲一区二区| 一级毛片aaaaaa免费看小| 日本-黄色视频高清免费观看| 国产欧美日韩综合在线一区二区| 91在线精品国自产拍蜜月| 99热这里只有是精品在线观看| 少妇的逼水好多| a级毛片在线看网站| 2018国产大陆天天弄谢| 丝袜在线中文字幕| 黑丝袜美女国产一区| 搡老乐熟女国产| 高清在线视频一区二区三区| 一级,二级,三级黄色视频| 91久久精品国产一区二区成人| 亚洲丝袜综合中文字幕| 欧美bdsm另类| 免费日韩欧美在线观看| 一本一本久久a久久精品综合妖精 国产伦在线观看视频一区 | 中文字幕人妻熟人妻熟丝袜美| 男女啪啪激烈高潮av片| 99久久综合免费| 午夜免费观看性视频| 少妇人妻 视频| 亚洲欧美中文字幕日韩二区| 十八禁网站网址无遮挡| 日韩三级伦理在线观看| 国产无遮挡羞羞视频在线观看| 亚洲久久久国产精品| 亚洲av欧美aⅴ国产| 欧美bdsm另类| 免费观看无遮挡的男女| 最近2019中文字幕mv第一页| 成年av动漫网址| 欧美激情国产日韩精品一区| 国产熟女午夜一区二区三区 | 美女脱内裤让男人舔精品视频| 2018国产大陆天天弄谢| 老司机亚洲免费影院| 香蕉精品网在线| 80岁老熟妇乱子伦牲交| 大话2 男鬼变身卡| 九九久久精品国产亚洲av麻豆| 2021少妇久久久久久久久久久| 18+在线观看网站| 男人操女人黄网站| 嫩草影院入口| 国产色婷婷99| 日韩av不卡免费在线播放| 久久久久视频综合| 国产精品久久久久久精品电影小说| 国产不卡av网站在线观看| 又大又黄又爽视频免费| 精品少妇内射三级| 乱码一卡2卡4卡精品| 免费观看av网站的网址| 黑人猛操日本美女一级片| 在线观看www视频免费| 日韩成人伦理影院| 母亲3免费完整高清在线观看 | 美女福利国产在线| 亚洲精品第二区| 精品人妻一区二区三区麻豆| 人人妻人人澡人人爽人人夜夜| 看非洲黑人一级黄片| 亚洲欧美日韩卡通动漫| 99九九在线精品视频| 久久99热6这里只有精品| 国产在线视频一区二区| 欧美激情国产日韩精品一区| 99热国产这里只有精品6| 国产在视频线精品| 成年女人在线观看亚洲视频| 国产精品成人在线| 日日摸夜夜添夜夜爱| 久久国产精品男人的天堂亚洲 | 亚洲av日韩在线播放| 九九久久精品国产亚洲av麻豆| 王馨瑶露胸无遮挡在线观看| 亚洲av成人精品一区久久| 免费观看在线日韩| 久久久久久久久大av| 久久 成人 亚洲| 久久久a久久爽久久v久久| 成人毛片a级毛片在线播放| 国产欧美另类精品又又久久亚洲欧美| 精品一区在线观看国产| 人妻人人澡人人爽人人| 99热这里只有精品一区| 久久女婷五月综合色啪小说| 欧美日韩在线观看h| 精品少妇久久久久久888优播| 国产一级毛片在线| 日本黄色日本黄色录像| 亚洲欧美成人综合另类久久久| 韩国高清视频一区二区三区| 十八禁网站网址无遮挡| 亚洲美女视频黄频| 免费日韩欧美在线观看| 国产不卡av网站在线观看| 日本午夜av视频| 日韩欧美精品免费久久| 少妇精品久久久久久久| freevideosex欧美| 亚洲欧美日韩卡通动漫| 十八禁高潮呻吟视频| av网站免费在线观看视频| 国产亚洲午夜精品一区二区久久| 在线观看www视频免费| 国产精品久久久久久久电影| 啦啦啦视频在线资源免费观看| 少妇的逼水好多| 国产成人精品婷婷| 久久久国产欧美日韩av| 亚洲怡红院男人天堂| 国产精品一区www在线观看| 夜夜看夜夜爽夜夜摸| 99久久综合免费| 日韩,欧美,国产一区二区三区| 在线观看国产h片| 午夜免费观看性视频| 精品午夜福利在线看| 精品人妻一区二区三区麻豆| 国产极品粉嫩免费观看在线 | 日本免费在线观看一区| 精品少妇内射三级| 日本黄色日本黄色录像| 成人漫画全彩无遮挡| 80岁老熟妇乱子伦牲交| 大片电影免费在线观看免费| av在线老鸭窝| 国产成人91sexporn| 女的被弄到高潮叫床怎么办| a级毛片在线看网站| 国产片特级美女逼逼视频| 夜夜看夜夜爽夜夜摸| videossex国产| 性高湖久久久久久久久免费观看| 日本av手机在线免费观看| 另类亚洲欧美激情| 国产精品三级大全| 美女中出高潮动态图| 亚洲av.av天堂| 久久精品国产a三级三级三级| 最黄视频免费看| 亚洲精品中文字幕在线视频| 国产精品偷伦视频观看了| 99热6这里只有精品| 久久久亚洲精品成人影院| 美女国产高潮福利片在线看| 永久网站在线| 插阴视频在线观看视频| 欧美日韩综合久久久久久| 毛片一级片免费看久久久久| 国产精品一区二区三区四区免费观看| 日本爱情动作片www.在线观看| 色5月婷婷丁香| 亚洲国产精品一区三区| 国产又色又爽无遮挡免| av一本久久久久| 赤兔流量卡办理| 免费久久久久久久精品成人欧美视频 | 亚洲精品乱码久久久久久按摩| 综合色丁香网| 乱人伦中国视频| 久久鲁丝午夜福利片| 免费播放大片免费观看视频在线观看| 曰老女人黄片| 人妻制服诱惑在线中文字幕| 自拍欧美九色日韩亚洲蝌蚪91| 韩国av在线不卡| 美女脱内裤让男人舔精品视频| 亚洲av电影在线观看一区二区三区| 精品人妻一区二区三区麻豆| 2022亚洲国产成人精品| 18禁在线无遮挡免费观看视频| av视频免费观看在线观看| 一区二区av电影网| 国产 精品1| 成人亚洲精品一区在线观看| 国国产精品蜜臀av免费| 国产在线一区二区三区精| 欧美精品一区二区免费开放| 久久99热6这里只有精品| 国产精品偷伦视频观看了| 国产精品一二三区在线看| 在线观看免费日韩欧美大片 | 欧美bdsm另类| 女的被弄到高潮叫床怎么办| 日韩伦理黄色片| 自线自在国产av| 亚洲欧美成人综合另类久久久| 国产深夜福利视频在线观看| 亚洲伊人久久精品综合| 乱码一卡2卡4卡精品| 国产精品国产三级国产av玫瑰| 又黄又爽又刺激的免费视频.| 色吧在线观看| 在线观看一区二区三区激情| 黑人高潮一二区| 欧美97在线视频| 一级a做视频免费观看| 久久这里有精品视频免费| 欧美另类一区| 欧美激情 高清一区二区三区| 成人18禁高潮啪啪吃奶动态图 | 精品久久国产蜜桃| 中国国产av一级| 亚洲欧美成人精品一区二区| 免费大片18禁| 看免费成人av毛片| 久久 成人 亚洲| 亚洲人成77777在线视频| 亚洲高清免费不卡视频| 熟女电影av网| 一本大道久久a久久精品| 亚洲av日韩在线播放| 熟女人妻精品中文字幕| 91午夜精品亚洲一区二区三区| 一边摸一边做爽爽视频免费| 欧美精品一区二区免费开放| 亚洲国产成人一精品久久久| 国产伦精品一区二区三区视频9| 新久久久久国产一级毛片| 亚州av有码| 成人综合一区亚洲| 午夜久久久在线观看| 国产不卡av网站在线观看| 久久人人爽av亚洲精品天堂| www.av在线官网国产| 99久久综合免费| 中文字幕人妻丝袜制服| a级毛片黄视频| 一级,二级,三级黄色视频| 丰满饥渴人妻一区二区三| 欧美日韩视频精品一区| 熟妇人妻不卡中文字幕| 成人国产av品久久久| 亚洲熟女精品中文字幕| 熟妇人妻不卡中文字幕| 日韩大片免费观看网站| 成年女人在线观看亚洲视频| 国产高清国产精品国产三级| 久久精品国产亚洲av涩爱| 午夜福利视频精品| 中文字幕人妻丝袜制服| 黄色视频在线播放观看不卡| 18在线观看网站| 91精品国产九色| 日韩制服骚丝袜av| 精品久久久噜噜| 夫妻午夜视频| 国产毛片在线视频| 午夜福利视频精品| 成年美女黄网站色视频大全免费 | 美女主播在线视频| 亚洲少妇的诱惑av| 少妇的逼好多水| 欧美激情 高清一区二区三区| av国产精品久久久久影院| 老司机影院毛片| 老司机影院成人| 汤姆久久久久久久影院中文字幕| 欧美日韩亚洲高清精品| 草草在线视频免费看| 国产成人精品一,二区| 欧美少妇被猛烈插入视频| 在线精品无人区一区二区三| 在线观看免费高清a一片| 两个人的视频大全免费| av不卡在线播放| 视频在线观看一区二区三区| 久久久久久久国产电影| 五月玫瑰六月丁香| 亚洲精华国产精华液的使用体验| 卡戴珊不雅视频在线播放| 精品熟女少妇av免费看| 麻豆精品久久久久久蜜桃| 最新中文字幕久久久久| 国产一区有黄有色的免费视频| 一级毛片电影观看| 亚洲精品美女久久av网站| 午夜老司机福利剧场| av网站免费在线观看视频| 欧美3d第一页| 色吧在线观看| 97超碰精品成人国产| 一二三四中文在线观看免费高清| 日本黄大片高清| 国产精品一二三区在线看| 精品国产一区二区久久| 日韩视频在线欧美| 国产精品女同一区二区软件| 精品国产一区二区三区久久久樱花| 人妻少妇偷人精品九色| 3wmmmm亚洲av在线观看| 少妇精品久久久久久久| 国产日韩一区二区三区精品不卡 | 欧美性感艳星| 成人亚洲精品一区在线观看| 欧美日韩av久久| 中文欧美无线码| av黄色大香蕉| 国产高清国产精品国产三级| .国产精品久久| 国产片特级美女逼逼视频| 国产精品无大码| 精品亚洲成国产av| 在线观看国产h片| 波野结衣二区三区在线| a 毛片基地| 最黄视频免费看| 午夜影院在线不卡| 综合色丁香网| 久久人妻熟女aⅴ| 国产视频内射| 一级毛片 在线播放| 亚洲精华国产精华液的使用体验| 王馨瑶露胸无遮挡在线观看| 日韩av在线免费看完整版不卡| 午夜精品国产一区二区电影| 免费久久久久久久精品成人欧美视频 | 国产免费又黄又爽又色| 人人妻人人添人人爽欧美一区卜| 精品熟女少妇av免费看| 妹子高潮喷水视频| av免费在线看不卡| 中国美白少妇内射xxxbb| 午夜老司机福利剧场| 久久久午夜欧美精品| 99热这里只有精品一区| 国国产精品蜜臀av免费| 亚洲av二区三区四区| 男人爽女人下面视频在线观看| av免费观看日本| 黄色欧美视频在线观看| 一区二区三区四区激情视频| 国产女主播在线喷水免费视频网站| 亚洲国产欧美日韩在线播放| 一级片'在线观看视频| 国产一区二区三区综合在线观看 | 精品久久久噜噜| 校园人妻丝袜中文字幕| 国产成人精品无人区| 国产精品蜜桃在线观看| 国产男人的电影天堂91| 欧美激情极品国产一区二区三区 | 久久久精品免费免费高清| 亚洲欧洲日产国产| 国产午夜精品一二区理论片| 91久久精品国产一区二区三区| 久久久国产欧美日韩av| 最近最新中文字幕免费大全7| 极品少妇高潮喷水抽搐| 26uuu在线亚洲综合色| 久久精品国产鲁丝片午夜精品| 亚洲一级一片aⅴ在线观看| 久久国产精品男人的天堂亚洲 | 日本wwww免费看| 国产日韩欧美亚洲二区| 在线观看免费日韩欧美大片 | 一区在线观看完整版| 蜜桃久久精品国产亚洲av| 免费播放大片免费观看视频在线观看| 水蜜桃什么品种好| 制服丝袜香蕉在线| 日本欧美国产在线视频| 成年美女黄网站色视频大全免费 | 99国产精品免费福利视频| 精品少妇久久久久久888优播| 亚洲精品美女久久av网站| 少妇猛男粗大的猛烈进出视频| 国产精品久久久久久久电影| 美女国产高潮福利片在线看| 亚洲av国产av综合av卡| 99热这里只有是精品在线观看| 国产精品一国产av| 18禁观看日本| 久久久久久伊人网av| 国产日韩欧美在线精品| 男女边吃奶边做爰视频| 国产精品不卡视频一区二区| 最近中文字幕高清免费大全6| 久久精品人人爽人人爽视色| 中文字幕精品免费在线观看视频 | 国产黄色免费在线视频| 少妇人妻精品综合一区二区| 肉色欧美久久久久久久蜜桃| 国产日韩一区二区三区精品不卡 | 精品久久久久久久久亚洲| 久久精品国产亚洲av天美| 一本色道久久久久久精品综合| 日韩精品免费视频一区二区三区 | 久久久久久久久大av| 久久这里有精品视频免费| 久久午夜福利片| 中文欧美无线码| 亚洲国产精品国产精品| 有码 亚洲区| 亚洲久久久国产精品| 国产精品久久久久成人av| 欧美bdsm另类| 在线播放无遮挡| 日本欧美国产在线视频| 亚洲欧洲国产日韩| av视频免费观看在线观看| 午夜福利影视在线免费观看| 五月伊人婷婷丁香| 国产精品国产av在线观看| 我要看黄色一级片免费的| 男人操女人黄网站| 99久久精品一区二区三区| 人妻人人澡人人爽人人| 婷婷色综合大香蕉| 80岁老熟妇乱子伦牲交| 美女cb高潮喷水在线观看|