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

    Statistical higher-order multi-scale method for nonlinear thermo-mechanical simulation of composite structures with periodically random configurations

    2024-02-28 11:53:44DONGHaoCUIJunzhi

    DONG Hao, CUI Jun-zhi

    (1.School of Mathematics and Statistics,Xidian University,Xi’an 710071,China;2.Academy of Mathematics and Systems Science,Chinese Academy of Sciences,Beijing 100190,China)

    Abstract: Stochastic multi-scale modeling and simulation for nonlinear thermo-mechanical problems of composite structures with complicated random microstructures remains a challenging issue.In this paper,we develop a novel statistical higher-order multi-scale (SHOMS) method for nonlinear thermo-mechanical simulation of composite structures with periodically random configurations,which is designed to overcome limitations of prohibitive computation involving the macro-scale and micro-scale.By virtue of statistical multi-scale asymptotic analysis and Taylor series method,the SHOMS computational model is rigorously derived for accurately analyzing nonlinear thermo-mechanical responses of random composite structures both in the macro-scale and micro-scale.Moreover,the local error analysis of SHOMS solutions in the point-wise sense clearly illustrates the crucial indispensability of establishing the higher-order asymptotic corrected terms in SHOMS computational model for keeping the conservation of local energy and momentum.Then,the corresponding space-time multi-scale numerical algorithm with off-line and on-line stages is designed to efficiently simulate nonlinear thermo-mechanical behaviors of random composite structures.Finally,extensive numerical experiments are presented to gauge the efficiency and accuracy of the proposed SHOMS approach.

    Key words: random composite structures;nonlinear thermo-mechanical simulation;SHOMS computational model;space-time multi-scale algorithm;local error analysis

    1 Introduction

    Random composite materials have been extensively applied in a variety of engineering sectors,such as aviation,aerospace,civil construction and smart structures,etc.By randomly distributing high-performance fibrous or particulate materials into ordinary matrix material,these synthetic composite materials exhibit high temperature resistance,high fatigue resistance and high fracture resistance,etc[1,2].Especially in aviation and aerospace industries,engineering structures often served under extreme heat environment while the thermal and mechanical properties of component materials exhibit obviously temperature-dependent feature.These complicated nonlinear behaviors and randomly geometric heterogeneities of the composite structures raise a grand challenge to effective numerical simulation[3].

    To the best of our knowledge,traditional numerical methods including the finite element method(FEM)[4],boundary element method[5]and meshless method[6]have been adopted to the simulation of nonlinear thermo-mechanical problems.Moreover,Abdoun et al.[7]used homotopy and asymptotic numerical method to simulate and analyze the thermal buckling and vibration of laminated composite plates with temperature-dependent properties.Najibi et al.[8]employed higher-order graded finite element method to conduct thermal stress analysis for a hollow FGM cylinder with nonlinear temperature-dependent material properties.In Ref.[9],the state space method and transfer-matrix method are adopted to obtain the displacements and stresses for the thick beams with temperature-dependent material properties under thermo-mechanical loads.However,the direct numerical simulation for composite materials needs a tremendous amount of computational resources or even ineffective to capture their microscopic behaviors due to highly heterogeneous components.

    To accomplish effective modeling and efficient simulation for inhomogeneous materials,scientists and engineers presented a variety of multi-scale methods,such as asymptotic homogenization method(AHM)[10],multi-scale finite element method (MsFEM)[11],heterogeneous multi-scale method(HMM)[12],variational multi-scale method(VMS)[13],multi-scale eigenelement method (MEM)[14],localized orthogonal decomposition method (LOD)[15]and finite volume based asymptotic homogenization theory(FVBAHT)[16],etc.However,numerical computation and theoretical analysis find that most of above-mentioned multi-scale methods are lower-order(first-order) multi-scale method in essence[17,18],which can only capture macroscopic and inadequate microscopic information of heterogeneous materials,especially for high-contrast composite materials.To improve inadequate numerical accuracy of classical lower-order multi-scale approaches,Cui and his research team systematically developed a class of higher-order (second-order) multi-scale methods,whose numerical accuracy is significantly improved for simulating authentic composite materials.Hence,these higher-order multi-scale approaches are extensively used in multi-physics coupling problem,stochastic problem,structural mechanics problem and nonlinear multi-scale problem of heterogeneous materials,etc[19-23].The reviews of above-mentioned multi-scale approaches show that these methods have a strong potential to encourage important advances in modeling and simulating a certain range of composites’ behaviors.However,they still need to be improved for random composite materials with complex non-deterministic microstructure.The uncertainties in the microstructure prominently affect the mechanical properties of the composite materials.Some stochastic multi-scale computational schemes have been established in recent years based on perturbation-based stochastic finite element method[24-26],spectral stochastic finite element method[27-30]and stochastic collocation method[31,32]for specific problems.Especially in Ref.[30],the spectral stochastic FEM is applied to analysis multiple-random field RC structures.Furthermore,combining with Monte Carlo method,the higher-order multi-scale methods proposed by Cui and his research team have been applied to simulate a wide range of physical behaviors of random composite materials[19,20,23,33,34].However,there are few works about multi-scale thermo-mechanical simulation of random composite materials with temperature-dependent properties.Hence,it is of great theoretical and engineering values to develop effective multi-scale approaches for nonlinear thermo-mechanical simulation of random composite materials.

    2 Statistical higher-order multi-scale computational model

    2.1 Microscopic representation of random composite

    The primary challenge for solving random multi-scale problems pertains to their auxiliary cell problems defined on the entire spaceΝ(Ν=2,3).To tackle this challenge,by using “periodization” and “cutoff” techniques in previous studies[19,20,23],the unit cell problems defined on macroscopic composite structureΩare approximated by transforming them into unit cell problems on a finite domain so-called microscopic unit cellYs(s=1,2,3,…denotes the index of random samples) with spatial sizeεand infinite random sampling,see Fig.1 for a schematic explanation,which corresponds to random sampleωsobeying a given probability distribution modelP(ω) with random variableω.

    In this study,the investigated composite structures are comprised of matrix and randomly distributed reinforced (or soften) particles or fibers,or porous media[35,36],as shown in Fig.2.Based on computer representation idea and its improved algorithm devised by Li et al.[37,38],we employ open-source Freefem++ software to establish the detailed computer representation algorithm for generating microscopic configurations of random composite materials as follows.

    Fig.1 Composite structures with periodically random configurations

    Fig.2 Random composite structures with different microscopic configurations

    (S1) Regarding random particulate and fibrous materials in Fig.2,the probability distribution model is first employed to generate the random geometric parameters (x1,x2,a,b,θ1) or (x1,x2,x3,a,b,c,θ1,θ2,θ3) for 2D or 3D randomly distributed configurations.

    (S2) Then,judge whether the newly generated configuration is located inside RVE and whether the newly generated configuration intersects with other previously generated configurations.We use whether the distance between the central points of the previously generated configurations and newly generated configuration is greater than the sum of the radii of both configurations as criterion.Additionally,to enhance the packing ratio of microscopic inclusions,we use the following rule as revised discriminate criterion:There exist such intersection points on previous configuration those connect the centers of previous ones with the points on the surfaces of new generated configuration.

    (S3) When generating a sufficient amount of microscopic configurations,mesh generation algorithm based on Delaunay Refinement method (Freefem++ command:“buildmesh” or “tetg” for 2D or 3D geometrical configurations,respectively) is adopted to create the microscopic configurations of the investigated random composites[39].

    For geometric parameters(x1,x2,a,b,θ1)in 2D case,x1andx2represent the central coordinates of the elliptical inclusion for thex-axis andy-axis.aandbdenote the lengths of the long half-axis and short half-axis of the elliptical inclusion andθ1represents the intersection angle between the long half-axis of the elliptical inclusion and thex-axis.For geometric parameters (x1,x2,x3,a,b,c,θ1,θ2,θ3) in 3D case,x1,x2andx3are the central coordinates of the ellipsoidal inclusion for thex-axis,y-axis andz-axis.a,bandcdenote the lengths of the long half-axis,middle half-axis and short half-axis of the ellipsoidal inclusion.θ1,θ2andθ3are the three Euler angles of the ellipsoidal inclusion,respectively.Moreover,by increasing the ratio of their long half-axis to short half-axis,the elliptical or ellipsoidal particles can be changed as fibrous inclusions.To sum up,the above methodologies accomplish the effective generation of finite element mesh for random composites.

    2.2 Stochastic multi-scale nonlinear thermo-mechanical problems

    Based on the classical thermo-mechanical model[3],the stochastic governing equations for describing the nonlinear thermo-mechanical problems of composite structures with periodically random configurations are set up,whose material parameters all possess the temperature-dependent properties.

    (1)

    where notationΩrepresents a bounded convex domain inΝ(Ν=2,3) with a boundary ?ΩT∪?Ωu∪?Ωq∪?Ωσ.In the micro-scale,the domainΩcan be defined by a statistical periodic layout of microscopic unit cellYscorresponding to random sampleωs.The characteristic size of microscopic cellYsis characterized by a parameterε.TheTε(x,t,ω) andε(x,Tε,ω) is the mass density;cε(x,Tε,ω) is specificthe second order thermal modulus tensor.Furthermore,we assume that all material parameters satisfy Lipschitz continuous condition with respect to temperature variableTεand are statistical periodicis the prescribed temperature on the boundaryis the prescribed heat flux normal to the boundary ?Ωqwith the normal vectorni,andis the prescribed traction on the boundary ?Ωσwith the normal vectoris the initial temperature;The internal heat source and body forces are represented byh(x,t) andfi(x,t),respectively.To begin with,let us sety=x/εas microscopic coordinates of statistical periodic unit cellYs=[0,1]N.With this notation,we have the chain rule for the spatial scales as follows.

    (2)

    which will be extensively used in the sequel.

    2.3 Statistical higher-order multi-scale computational model

    (3)

    By aid of the above Taylor’s formula(4) and multi-index notation(5),the material parameters depending on temperatureTεthus can be expanded as follows[22]

    kij(y,T0+εT1+ε2T2+O(ε3),ω)=

    kij(y,T0,ω)+D(0,1,0)kij(y,T0,ω)[εT1+ε2T2+

    O(ε3)]2+O([εT1+ε2T2+O(ε3)]3)=

    kij(y,T0,ω)+εT1D(0,1,0)kij+

    (4)

    cε(x,Tε,ω)=c(0)+εc(1)+ε2c(2)+O(ε3)

    (5)

    By substituting equations (3)-(5) into the multi-scale coupling problem (1) and utilizing the chain rule (2),we can expand the derivatives and match terms with the same order of the small parameterε.Then,according to the similar higher-order multi-scale analysis in Ref.[22],O(ε-2)-order equations allow us to reasonably obtain the following results

    (6)

    And then,according to O(ε-1)-order equations,the first-order correctorsT1andui1can be individually constructed as follows

    (7)

    (8)

    (9)

    (10)

    Subsequently,we perform a volume integral on both sides ofO(ε0)-order equations on microscopic unit cellYsand apply the Gauss theorem onO(ε-1)-order equations.By further applying the Kolmogorov’s strong law of large numbers,the macroscopic homogenized equations associated with multi-scale problem (1) are derived as follows

    (11)

    where the macroscopic homogenized material parameters in (11) are defined as follows

    (12)

    The effective material coefficients in (12) are evaluated using the following formulas,which correspond to microscopic unit cellYswith random sampleωs

    (13)

    After that,by substituting (6,7) intoO(ε0)-order equations,and then subtractingO(ε0)-order equations from (11),the concrete expressions forT2andui2can be defined as follows

    (14)

    (15)

    (16)

    (17)

    (18)

    (19)

    (20)

    (21)

    (22)

    (23)

    (24)

    Summing up,the macro-micro coupled SHOMS asymptotic solutions of the multi-scale nonlinear dynamic thermo-mechanical problem (1) are established.

    2.4 Local error analysis of statistical multi-scale solutions in point-wise sense

    To begin with,two kinds of the multi-scale approximate solutions of stochastic multi-scale problem (1) are defined as follows

    (25)

    Next,we introduce the following residual functions for thelocal numerical accuracy analysis.

    (26)

    (27)

    (28)

    Noting the residual equations (27),it can be concluded that the residual error of SLOMS solutions is of order-O (1) in the point-wise sense,primarily due to the presence of terms F0and S0i.Comparing with the residual equations (28),the residual error of SHOMS solutions is of order-O(ε) in the point-wise sense.This implies that the SHOMS solutions can satisfy the local energy conservation of thermal equation and local momentum conservation of mechanical equations of the original stochastic multi-scale equations (1) in the point-wise sense,which can still provide the required accuracy for engineering computation and accurately capture the microscopic oscillating behaviors exhibited by random composite materials.

    3 Space-time multi-scale numerical algorithm

    In this section,a new space-time multi-scale numerical algorithm with off-line and on-line stages is presented for the stochastic nonlinear governing equations (1) based on the FEM in spatial region and the finite difference method (FDM) in time direction.The detailed algorithm is listed as follows.

    3.1 Off-line computation for microscopic cell problems

    3.2 On-line computation for macroscopic homogenized problem and statistical higher-order multi-scale solutions

    (3) It should be highlighted that,decoupled thermal equation is a nonlinear system which can not be computed directly.Herein,the direct iterative method is employed for simulating the nonlinear system.

    (4) For arbitrary point (x,t)∈Ω×(0,T*),to apply the interpolation method to acquire the corresponding values of first-order auxiliary cell functions,second-order auxiliary cell functions and homogenized solutions.The spatial derivatives in formulas (7,14) are evaluated by the average technique on relative elements[40],and the temporal derivatives in formula (14) are evaluated by using the difference scheme at every time steps.Then,the temperature fieldT(2ε)(x,t,ω) and displacement fieldu(2ε)(x,t,ω) are computed.Moreover,we can further employ the higher-order interpolation method and post-processing technique in[41,42]to obtain the high-accuracy SHOMS solutions.

    4 Numerical experiments

    In this section,the numerical accuracy andefficiency of the proposed SHOMS method are demonstrated by several numerical examples.Since obtaining the exact solutions for the multi-scale nonlinear problem (1) is extremely difficult or even impossible,direct numerical simulation (DNS) solutions on the fine grid for the multi-scale nonlinear problem (1) are taken as the reference solutions denoted asTDNS(x,t) anduDNS(x,t).In the following numerical experiments,some error notations are introduced,Terr0,Terr1 and Terr2 represent theL2norm errors,TErr0,TErr1 and TErr2 represent theH1semi-norm errors for macroscopic homogenized solution,SLOMS solutions and SHOMS solutions of temperature field respectively;Uerr0,Uerr1 and Uerr2 represent theL2norm errors,UErr0,UErr1 and UErr2 represent theH1semi-norm errors for macroscopic homogenized solution,SLOMS solutions and SHOMS solutions of displacement field respectively.

    Example 1Validation of the SHOMS me-thod for nonlinear thermo-mechanical simulation in periodic composite structure.

    A 2D composite structure with periodically microscopic configurations is investigated here,whose macrostructureΩand unit cellY^sare shown in Fig.3,whereΩ=(x1,x2)=[0,1]×[0,1] cm2and periodic unit cell sizeε=1/5.This 2D composite structure is clamped on its four boundaries,and the temperature at the boundaries is kept at 373.15 K.In addition,the internal heat source and body forces are set ash=5000 J/(cm2·s) and (f1,f2)=(-2000,-2000) N/cm2.And,the material parameters of 2D composite structure are demonstrated in Table 1.

    Fig.3 Illustration of investigated 2D composite structure Implementing the SHOMS method to multi-scale nonlinear coupling equations (1) in time interval t∈[0,1] s with temporal step Δt=0.002 s,we define service temperature range [273.15,873.15] K of investigated composite structure and 60 representative macroscopic parameters in temperature range.In this example,the total auxiliary cell problems need to be solved off-line 4380 times,in which the 13 first-order cell functions and 60 second-order cell functions are solved on 60 macroscopic temperature interpolation points.The information of FEM meshes is listed in Tab.2.After off-line computation for microscopic cell problems and on-line computation for macroscopic homogenized problems and higher-order multi-scale solutions,we depict the computational results in Fig.4~Fig.8.

    Tab.1 Material parameters of 2D composite structure

    Tab.2 Summary of computational cost

    Fig.4 Numerical results of temperature field at t=0.2 s

    Fig.5 Numerical results of displacement field at t=0.2 s

    Fig.6 Numerical results of temperature field at t=1.0 s

    Fig.7 Numerical results of displacement field at t=1.0 s

    Fig.8 Evolutive relative errors of temperature and displacement fields According to the computational resource cost in Tab.2,the SHOMS method can greatly economize computer memory without losing precision.Actually,both the SHOMS method and direct numerical simulation are performed on a HP desktop workstation equipped with an Intel(R) Core(TM) i7-8750H processor (2.20 GHz) and 16.0 GB of internal memory.As illustrated in Fig.4~Fig.7,we can conclude that the higher-order multi-scale solutions can accurately capture the microscopic oscillatory behaviors and preferably approximate the exact solutions of the investigated 2D composite structure compared with macroscopic homogenized solutions and lower-order multi-scale solutions,especially for temperature field.From the evolutive relative errors in Fig.8,it can clearly demonstrate that the two-stages space-time multi-scale numerical algorithm is accurate and stable in the long-time numerical simulation.Furthermore,it is worth emphasizing that the presented SHOMS approach remains effective even for a relatively small parameter ε,namely existing a great number of microscopic unit cells in inhomogeneous structures.At this time,the high-resolution DNS simulation can not guarantee the convergence for the investigated large-scale problems.This obvious advantage of the SHOMS approach is of great application values for engineering computation.

    Example 2Application of the SHOMS method for equivalent material parameters computation of random composite structure.

    In this example,two kinds of composite materials with matrix Ti-6Al-4V and random inclusion ZrO2,and matrix SiC and random inclusion C are investigated by the SHOMS method,as exhibited in Fig.9.The detailed material parameters for the investigated composite materials are presented in the following Tab.3 and Tab.4.

    Fig.9 Several random RVEs of the composites employed for predicting equivalent material parameters

    Tab.3 Material property parameters of Ti-6Al-4V/ZrO2 composite

    Tab.4 Material property parameters of SiC/C composite

    By using the SHOMS method,the equivalent material parameters at macro-scale are obtained by the mean value of 50 randomly microscopic samples.The corresponding results are depicted in Fig.10.According to the numerical results in Fig.10,we can conclude that the predictive values of Ti-6Al-4V/ZrO2composite and SiC/C composite fall between lower and upper bounds of Voigt-Reuss method,Hashin-Shtrikman method and also approximate the predicted values of Hobbs method.Hence,the proposed SHOMS can be employed to predict the temperature-dependent equivalent material properties of Ti-6Al-4V/ZrO2composite and SiC/C composite.

    Fig.10 A comparison of the predictive results of equivalent material parameters

    Example 3Application of the SHOMS method fornonlinear thermo-mechanical simulation in random composite structure.

    This example mainly study the nonlinear thermo-mechanical simulation of 2D composite structure with randomly microscopic configurations,as depicted in Fig.11.In addition,the setting of initial-boundary conditions,heat source,body forces and material parameters in this example is the same as those of Example 1.

    Fig.11 Illustration of investigated 2D composite structure Applying thenovel SHOMS method to multi-scale nonlinear coupling equations (1) within time interval t∈[0,1] s with temporal step Δt=0.002 s,we establish the service temperature range of the investigated composite structure as [273.15,873.15] K.In this service temperature range,we distribute 60 representative macroscopic parameters detailed information of FEM meshes is listed in Tab.5.After off-line computation for microscopic cell problems and on-line computation for macroscopic homogenized problems and higher-order multi-scale solutions,we present the computational results in Fig.12~Fig.16.

    Tab.5 Summary of computational cost

    Fig.12 Numerical results of temperature field at t=0.2 s

    Fig.13 Numerical results of displacement field at t=0.2 s

    Fig.14 Numerical results of temperature field at t=1.0 s

    Fig.15 Numerical results of displacement field at t=1.0 s As indicated Tab.5,thenovel SHOMS method can significantly reduce computer memory without losing precision.Moreover,the numerical results in Fig.12~Fig.15 reveal that the higher-order multi-scale solutions can accurately capture the microscopic oscillatory behaviors and provide preferable approximations to the exact solutions compared with macroscopic homogenized solutions

    Fig.16 Evolutive relative errors of temperature and displacement fields and lower-order multi-scale solutions,especially for temperature field.The evolutive relative errors shown in Fig.16 clearly demonstrate the accuracy and stability of the two-stages space-time multi-scale numerical algorithm in the long-time numerical simulation.Furthermore,it is important to highlight that the SHOMS approach remains effective even for a relatively small parameter ε,which corresponds to a large number of microscopic unit cells in heterogeneous structures.In contrast,the high-resolution DNS simulation fails to converge for the investigated large-scale problems.This prominent computational advantage is of significant practical value of the SHOMS approach in engineering computations.

    5 Conclusions and outlook

    In the present work,a novel statistical higher-order multi-scale method is developed for effectively simulating nonlinear thermo-mechanical problems of composite structures with periodically random configurations and temperature-dependent properties,which serve under extreme heat environment.The main contributions of this work are threefold:First,the statistical multi-scale formulations with the higher-order correction terms are established for composite structures with periodically random configurations.Second,the local error estimations for the statistical multi-scale solutions of nonlinear thermo-mechanical systems are derived in detail.Third,a space-time numerical algorithm with off-line and on-line stages is designed to overcome the prohibitive computation of direct numerical simulation.Furthermore,numerical results demonstrate that the presented SHOMS approach can effectively simulate nonlinear thermo-mechanical coupling behaviors with less computational cost and accurately capture the microscopic oscillatory information caused by randomly heterogeneous configurations.Besides,the proposed SHOMS approach can accurately predict the equivalent material parameters of random composite structures compared with the predictive results of some theoretical models,which illustrate that high temperature field has a remarkable effect on macroscopic thermo-mechanical properties.

    In the future,the SHOMS method will be extended to more complex nonlinear problems including thermal convection and radiation effects under extreme thermal environment.Additionally,machine learning approaches and parallel algorithm will be introduced in the off-line stage of SHOMS framework,in order to avoid repetitive statistical computation and improve computational efficiency.

    午夜福利乱码中文字幕| 亚洲伊人久久精品综合| 国产人伦9x9x在线观看| 高清欧美精品videossex| 超碰97精品在线观看| 亚洲国产毛片av蜜桃av| 老司机影院毛片| 后天国语完整版免费观看| 伦理电影免费视频| 中文字幕精品免费在线观看视频| 欧美+亚洲+日韩+国产| 他把我摸到了高潮在线观看 | 两个人看的免费小视频| 欧美精品亚洲一区二区| 亚洲欧美色中文字幕在线| 午夜福利在线免费观看网站| 午夜精品久久久久久毛片777| 超碰97精品在线观看| avwww免费| 90打野战视频偷拍视频| 日本一区二区免费在线视频| 人人妻人人添人人爽欧美一区卜| 亚洲五月色婷婷综合| 国产一区有黄有色的免费视频| 极品教师在线免费播放| 亚洲色图 男人天堂 中文字幕| 亚洲精品国产色婷婷电影| 亚洲欧美一区二区三区久久| 久久人人97超碰香蕉20202| 亚洲精品美女久久久久99蜜臀| 午夜福利在线观看吧| 老司机午夜福利在线观看视频 | 欧美老熟妇乱子伦牲交| 亚洲精华国产精华精| 欧美成人免费av一区二区三区 | 国产精品98久久久久久宅男小说| 亚洲第一av免费看| 中文字幕人妻丝袜制服| 99热网站在线观看| 美国免费a级毛片| 久久久久久久国产电影| 日韩三级视频一区二区三区| 麻豆国产av国片精品| 欧美精品av麻豆av| 亚洲精品一卡2卡三卡4卡5卡| 久久国产精品男人的天堂亚洲| 99国产综合亚洲精品| 国产极品粉嫩免费观看在线| 久久青草综合色| 黄网站色视频无遮挡免费观看| 欧美午夜高清在线| 久久久久国内视频| 午夜激情久久久久久久| 欧美国产精品va在线观看不卡| 男女午夜视频在线观看| 日韩欧美免费精品| 久久久久久亚洲精品国产蜜桃av| 91九色精品人成在线观看| 亚洲 国产 在线| 性高湖久久久久久久久免费观看| 国产精品美女特级片免费视频播放器 | 国产在线一区二区三区精| 在线天堂中文资源库| 美女主播在线视频| 窝窝影院91人妻| 妹子高潮喷水视频| 十八禁网站免费在线| 欧美人与性动交α欧美精品济南到| 国产日韩一区二区三区精品不卡| 国产精品亚洲一级av第二区| 91麻豆精品激情在线观看国产 | 满18在线观看网站| 自拍欧美九色日韩亚洲蝌蚪91| 高清欧美精品videossex| 精品亚洲成国产av| 国产有黄有色有爽视频| 一级片'在线观看视频| 老司机深夜福利视频在线观看| 国产亚洲精品第一综合不卡| 香蕉丝袜av| 最新在线观看一区二区三区| 高清欧美精品videossex| 欧美在线黄色| av网站免费在线观看视频| 午夜激情久久久久久久| 精品少妇一区二区三区视频日本电影| 18禁国产床啪视频网站| 久久九九热精品免费| 精品第一国产精品| 飞空精品影院首页| 国产亚洲欧美在线一区二区| 女警被强在线播放| 亚洲欧美一区二区三区久久| 欧美av亚洲av综合av国产av| 母亲3免费完整高清在线观看| 免费久久久久久久精品成人欧美视频| 亚洲欧美精品综合一区二区三区| 亚洲视频免费观看视频| a在线观看视频网站| 50天的宝宝边吃奶边哭怎么回事| 999精品在线视频| 天堂动漫精品| 三级毛片av免费| 男男h啪啪无遮挡| av天堂在线播放| 久久中文字幕人妻熟女| 可以免费在线观看a视频的电影网站| 夜夜骑夜夜射夜夜干| 免费少妇av软件| 亚洲国产精品一区二区三区在线| 国产精品久久电影中文字幕 | 韩国精品一区二区三区| 日韩欧美免费精品| svipshipincom国产片| 侵犯人妻中文字幕一二三四区| 他把我摸到了高潮在线观看 | 日韩欧美国产一区二区入口| 亚洲精品中文字幕一二三四区 | 亚洲熟妇熟女久久| 久久久精品国产亚洲av高清涩受| 精品国产国语对白av| 国产精品久久久av美女十八| 一级片免费观看大全| 亚洲欧洲精品一区二区精品久久久| 精品少妇一区二区三区视频日本电影| 一本一本久久a久久精品综合妖精| 亚洲天堂av无毛| 国产精品电影一区二区三区 | 一个人免费在线观看的高清视频| 成年女人毛片免费观看观看9 | 久久久精品区二区三区| 久久久久精品人妻al黑| 国产麻豆69| 又黄又粗又硬又大视频| 久久亚洲精品不卡| 免费人妻精品一区二区三区视频| 亚洲 欧美一区二区三区| 国产欧美日韩一区二区三区在线| 欧美激情高清一区二区三区| 国产成人啪精品午夜网站| 精品福利观看| 久久av网站| 操出白浆在线播放| 欧美日韩国产mv在线观看视频| 亚洲国产毛片av蜜桃av| 欧美日韩成人在线一区二区| 男男h啪啪无遮挡| 无遮挡黄片免费观看| 狠狠狠狠99中文字幕| 国产激情久久老熟女| 精品国产亚洲在线| 国产一区二区三区在线臀色熟女 | 久久香蕉激情| 香蕉久久夜色| 午夜福利乱码中文字幕| 欧美日韩福利视频一区二区| 99久久99久久久精品蜜桃| 丁香欧美五月| a在线观看视频网站| 99精品在免费线老司机午夜| 色在线成人网| 18在线观看网站| 日日夜夜操网爽| tube8黄色片| 中亚洲国语对白在线视频| 久久人人爽av亚洲精品天堂| 99国产极品粉嫩在线观看| 亚洲国产精品一区二区三区在线| 黑人欧美特级aaaaaa片| 99国产极品粉嫩在线观看| 国产成人精品无人区| 欧美大码av| 女人被躁到高潮嗷嗷叫费观| 性色av乱码一区二区三区2| 欧美国产精品一级二级三级| 日韩制服丝袜自拍偷拍| 狠狠狠狠99中文字幕| 搡老乐熟女国产| 丝袜在线中文字幕| 天堂中文最新版在线下载| 高清在线国产一区| 韩国精品一区二区三区| 久久亚洲真实| 看免费av毛片| 亚洲国产毛片av蜜桃av| 久久九九热精品免费| 一区二区三区乱码不卡18| 国产高清激情床上av| 成年人黄色毛片网站| 国产亚洲精品一区二区www | 亚洲av国产av综合av卡| 国产av国产精品国产| 精品亚洲成a人片在线观看| 中亚洲国语对白在线视频| 国产99久久九九免费精品| 亚洲少妇的诱惑av| 十八禁人妻一区二区| 啦啦啦视频在线资源免费观看| 国产成人系列免费观看| 久9热在线精品视频| 在线观看舔阴道视频| 欧美黑人精品巨大| 欧美日韩一级在线毛片| 一边摸一边抽搐一进一出视频| 亚洲精品美女久久av网站| 夜夜夜夜夜久久久久| 激情在线观看视频在线高清 | 日韩有码中文字幕| 亚洲色图av天堂| 国产在线视频一区二区| 欧美一级毛片孕妇| 999久久久精品免费观看国产| 成在线人永久免费视频| 91大片在线观看| 亚洲精品在线观看二区| 久久天躁狠狠躁夜夜2o2o| 脱女人内裤的视频| 两个人免费观看高清视频| 亚洲色图 男人天堂 中文字幕| 夜夜骑夜夜射夜夜干| a级片在线免费高清观看视频| 热99re8久久精品国产| 老司机在亚洲福利影院| 亚洲第一青青草原| 少妇的丰满在线观看| 在线观看www视频免费| 免费观看人在逋| 国产人伦9x9x在线观看| 亚洲美女黄片视频| 后天国语完整版免费观看| 成年人午夜在线观看视频| 99精品欧美一区二区三区四区| 99热网站在线观看| 精品少妇一区二区三区视频日本电影| 久久中文字幕一级| 亚洲欧美日韩高清在线视频 | 亚洲情色 制服丝袜| 精品亚洲成a人片在线观看| 十八禁人妻一区二区| 一区二区三区乱码不卡18| 美女高潮喷水抽搐中文字幕| 狂野欧美激情性xxxx| 欧美日韩黄片免| 一个人免费在线观看的高清视频| 黑人欧美特级aaaaaa片| 在线观看免费视频网站a站| 久久午夜亚洲精品久久| 五月天丁香电影| 日本黄色日本黄色录像| 国产精品久久电影中文字幕 | 久久亚洲精品不卡| 老司机午夜十八禁免费视频| 热re99久久精品国产66热6| 精品免费久久久久久久清纯 | 每晚都被弄得嗷嗷叫到高潮| 麻豆av在线久日| 窝窝影院91人妻| 国产野战对白在线观看| 黄片小视频在线播放| 黑人巨大精品欧美一区二区蜜桃| 国产精品av久久久久免费| www.精华液| 黑人巨大精品欧美一区二区mp4| 久久久水蜜桃国产精品网| 免费观看a级毛片全部| 一区二区三区国产精品乱码| 老熟妇仑乱视频hdxx| 国产不卡av网站在线观看| 亚洲人成电影免费在线| 欧美日韩国产mv在线观看视频| 丁香六月天网| 久久中文字幕人妻熟女| √禁漫天堂资源中文www| 国产精品免费视频内射| 久久亚洲精品不卡| 欧美精品一区二区免费开放| 免费在线观看日本一区| 91字幕亚洲| 国产精品av久久久久免费| 老司机福利观看| 国产精品1区2区在线观看. | 成人特级黄色片久久久久久久 | av片东京热男人的天堂| 一二三四在线观看免费中文在| 亚洲全国av大片| 国产一区二区激情短视频| 久久久精品国产亚洲av高清涩受| 午夜福利影视在线免费观看| 国产不卡一卡二| 99在线人妻在线中文字幕 | 午夜激情久久久久久久| 亚洲伊人色综图| 老鸭窝网址在线观看| netflix在线观看网站| 亚洲av美国av| 多毛熟女@视频| 天天影视国产精品| 久久九九热精品免费| 亚洲成人手机| 91成人精品电影| 欧美在线一区亚洲| 国产人伦9x9x在线观看| 久久国产精品影院| 欧美黑人精品巨大| 成在线人永久免费视频| 高清av免费在线| 国产精品国产高清国产av | 国产真人三级小视频在线观看| 成在线人永久免费视频| 狂野欧美激情性xxxx| 两性午夜刺激爽爽歪歪视频在线观看 | 久久中文看片网| 99riav亚洲国产免费| 国产精品一区二区免费欧美| 飞空精品影院首页| 黄色视频,在线免费观看| 亚洲成a人片在线一区二区| 老司机在亚洲福利影院| 亚洲男人天堂网一区| 精品国产乱码久久久久久小说| 黑人欧美特级aaaaaa片| 韩国精品一区二区三区| 欧美激情高清一区二区三区| 亚洲情色 制服丝袜| 一区福利在线观看| 色综合欧美亚洲国产小说| 免费看a级黄色片| 国产精品美女特级片免费视频播放器 | 天天添夜夜摸| 亚洲成人手机| 国产免费现黄频在线看| 国产免费视频播放在线视频| 国产在线视频一区二区| 高潮久久久久久久久久久不卡| 1024香蕉在线观看| svipshipincom国产片| 亚洲第一av免费看| 亚洲第一欧美日韩一区二区三区 | 日本av免费视频播放| 侵犯人妻中文字幕一二三四区| 满18在线观看网站| 成年版毛片免费区| 人成视频在线观看免费观看| 五月天丁香电影| 18禁美女被吸乳视频| 亚洲午夜理论影院| 日韩视频一区二区在线观看| 精品国产国语对白av| 90打野战视频偷拍视频| 亚洲av片天天在线观看| 久久国产精品大桥未久av| 日日摸夜夜添夜夜添小说| 99re6热这里在线精品视频| 在线播放国产精品三级| 人妻 亚洲 视频| 自拍欧美九色日韩亚洲蝌蚪91| 丁香欧美五月| 婷婷丁香在线五月| 亚洲国产毛片av蜜桃av| 久久久国产一区二区| 国产精品麻豆人妻色哟哟久久| 精品久久久久久电影网| 免费黄频网站在线观看国产| 最新在线观看一区二区三区| 亚洲美女黄片视频| 老司机福利观看| 丁香欧美五月| av网站在线播放免费| 在线观看人妻少妇| 777米奇影视久久| 亚洲国产毛片av蜜桃av| 欧美成人免费av一区二区三区 | 欧美黄色片欧美黄色片| 丰满迷人的少妇在线观看| 91麻豆精品激情在线观看国产 | 老司机福利观看| 国产精品麻豆人妻色哟哟久久| 国产区一区二久久| 国产成人系列免费观看| 成人三级做爰电影| 国产av一区二区精品久久| 亚洲欧洲日产国产| 热99久久久久精品小说推荐| 国产精品一区二区免费欧美| 久久av网站| 欧美久久黑人一区二区| 国产一区二区三区综合在线观看| av国产精品久久久久影院| 人人妻人人爽人人添夜夜欢视频| 一级毛片女人18水好多| 亚洲成av片中文字幕在线观看| 91老司机精品| 岛国毛片在线播放| 热re99久久国产66热| 国产精品自产拍在线观看55亚洲 | 精品人妻1区二区| 欧美日韩福利视频一区二区| 日韩欧美国产一区二区入口| 狂野欧美激情性xxxx| 少妇精品久久久久久久| 老司机福利观看| 久久精品熟女亚洲av麻豆精品| 成人18禁在线播放| 国产91精品成人一区二区三区 | 成人亚洲精品一区在线观看| 亚洲av成人一区二区三| 久久中文字幕人妻熟女| 一边摸一边抽搐一进一出视频| 天天操日日干夜夜撸| 美女国产高潮福利片在线看| 高清毛片免费观看视频网站 | 9色porny在线观看| 成年人免费黄色播放视频| 免费一级毛片在线播放高清视频 | 午夜激情av网站| 日韩人妻精品一区2区三区| 国产精品99久久99久久久不卡| 中亚洲国语对白在线视频| 黄网站色视频无遮挡免费观看| 成人手机av| 老汉色∧v一级毛片| 女警被强在线播放| videos熟女内射| 久久精品国产99精品国产亚洲性色 | 性色av乱码一区二区三区2| 热99国产精品久久久久久7| 国产成人精品无人区| 变态另类成人亚洲欧美熟女 | 精品亚洲乱码少妇综合久久| 一本—道久久a久久精品蜜桃钙片| 免费黄频网站在线观看国产| √禁漫天堂资源中文www| 国产精品国产av在线观看| 久久ye,这里只有精品| 国产91精品成人一区二区三区 | 久久久精品94久久精品| 亚洲,欧美精品.| 啦啦啦中文免费视频观看日本| 久久久水蜜桃国产精品网| 欧美激情极品国产一区二区三区| √禁漫天堂资源中文www| 午夜福利视频在线观看免费| 亚洲一卡2卡3卡4卡5卡精品中文| 9热在线视频观看99| 在线天堂中文资源库| 国产精品亚洲av一区麻豆| 大型黄色视频在线免费观看| 桃花免费在线播放| 国产精品香港三级国产av潘金莲| 国产精品一区二区在线不卡| 国产精品二区激情视频| 精品少妇内射三级| 国产伦人伦偷精品视频| 亚洲精品一卡2卡三卡4卡5卡| 高清毛片免费观看视频网站 | xxxhd国产人妻xxx| 老司机影院毛片| 欧美性长视频在线观看| 这个男人来自地球电影免费观看| 老司机靠b影院| 免费不卡黄色视频| 日韩 欧美 亚洲 中文字幕| www.熟女人妻精品国产| 国产在线视频一区二区| 夜夜骑夜夜射夜夜干| 欧美精品一区二区大全| 999久久久精品免费观看国产| bbb黄色大片| 一个人免费看片子| 在线播放国产精品三级| 大香蕉久久网| 夜夜骑夜夜射夜夜干| 交换朋友夫妻互换小说| 蜜桃国产av成人99| 99久久精品国产亚洲精品| 国产xxxxx性猛交| 成人18禁高潮啪啪吃奶动态图| 19禁男女啪啪无遮挡网站| 欧美日韩国产mv在线观看视频| 在线看a的网站| 国产视频一区二区在线看| 成年人黄色毛片网站| 成人国产一区最新在线观看| 婷婷丁香在线五月| 亚洲成av片中文字幕在线观看| 一区二区三区国产精品乱码| 满18在线观看网站| 久久久精品区二区三区| e午夜精品久久久久久久| 欧美在线黄色| 国产精品一区二区免费欧美| 俄罗斯特黄特色一大片| 亚洲色图av天堂| 成年人黄色毛片网站| 国产日韩欧美亚洲二区| 少妇的丰满在线观看| 考比视频在线观看| 精品国产乱子伦一区二区三区| 香蕉丝袜av| 亚洲欧美精品综合一区二区三区| 69精品国产乱码久久久| 久久国产精品男人的天堂亚洲| av天堂久久9| 国产三级黄色录像| 色婷婷久久久亚洲欧美| tocl精华| 一个人免费看片子| 国产精品一区二区在线不卡| 国产精品自产拍在线观看55亚洲 | 国产免费av片在线观看野外av| 亚洲精品久久午夜乱码| 丝袜美腿诱惑在线| 国产精品98久久久久久宅男小说| 视频区图区小说| 欧美在线黄色| 黄色a级毛片大全视频| 久久人妻福利社区极品人妻图片| 亚洲精品成人av观看孕妇| 欧美日韩福利视频一区二区| 欧美在线黄色| 日日夜夜操网爽| 一夜夜www| 黄网站色视频无遮挡免费观看| 性少妇av在线| 国产精品国产av在线观看| 久久 成人 亚洲| 国产精品久久久久久人妻精品电影 | 18禁美女被吸乳视频| 亚洲中文字幕日韩| 日韩大片免费观看网站| 亚洲精品一二三| 午夜日韩欧美国产| 18在线观看网站| 老熟女久久久| 久久久精品94久久精品| 一个人免费在线观看的高清视频| 国产精品99久久99久久久不卡| 别揉我奶头~嗯~啊~动态视频| 少妇精品久久久久久久| 日本av手机在线免费观看| 香蕉丝袜av| 日韩一区二区三区影片| 国产区一区二久久| a级片在线免费高清观看视频| 99久久99久久久精品蜜桃| 一边摸一边抽搐一进一小说 | 露出奶头的视频| 亚洲一区中文字幕在线| 99精品久久久久人妻精品| 亚洲国产看品久久| 日日夜夜操网爽| 久久久久久久精品吃奶| 亚洲 国产 在线| 欧美中文综合在线视频| 日韩三级视频一区二区三区| 黄色视频在线播放观看不卡| 别揉我奶头~嗯~啊~动态视频| 天天影视国产精品| 国精品久久久久久国模美| 51午夜福利影视在线观看| 999久久久精品免费观看国产| 欧美一级毛片孕妇| 一二三四社区在线视频社区8| 高清欧美精品videossex| 波多野结衣一区麻豆| 丁香欧美五月| 国产99久久九九免费精品| 国产精品久久久久久人妻精品电影 | 最近最新免费中文字幕在线| 韩国精品一区二区三区| 搡老乐熟女国产| 国产精品av久久久久免费| 精品一区二区三区视频在线观看免费 | 精品一品国产午夜福利视频| 亚洲 欧美一区二区三区| 大型av网站在线播放| 丝袜在线中文字幕| 多毛熟女@视频| 国产成人精品在线电影| 黄色成人免费大全| 亚洲免费av在线视频| 美女高潮喷水抽搐中文字幕| 亚洲va日本ⅴa欧美va伊人久久| 国产欧美日韩一区二区精品| 久久狼人影院| 一区二区三区国产精品乱码| 欧美亚洲 丝袜 人妻 在线| 午夜福利视频在线观看免费| 免费在线观看日本一区| 最近最新中文字幕大全免费视频| 纯流量卡能插随身wifi吗| 黄色丝袜av网址大全| 99精国产麻豆久久婷婷| 久久久久久免费高清国产稀缺| 大片电影免费在线观看免费| 国产亚洲欧美在线一区二区| 少妇 在线观看| 国产成人欧美在线观看 | 操出白浆在线播放| 两个人免费观看高清视频| 午夜91福利影院| 制服诱惑二区| 久久狼人影院| 大片免费播放器 马上看| 久久九九热精品免费| 国产xxxxx性猛交| 中文字幕高清在线视频| 视频在线观看一区二区三区| 丁香六月天网| 中文字幕av电影在线播放| 一级片免费观看大全| 极品少妇高潮喷水抽搐| 脱女人内裤的视频| 久久免费观看电影| 无遮挡黄片免费观看|