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

    Optimal Analysis for Shakedown of Functionally Graded(FG)Bree Plate with Genetic Algorithm

    2014-04-17 06:17:53ZhengPeng4andHu
    Computers Materials&Continua 2014年7期

    H.Zheng,X.Peng,3,4and N.Hu

    1 Introduction

    The rapid development of aircraft and space technologies requires materials to work in more and more severe environments,such as high and varying temperature with large temperature gradients.With the distributions of the thermal and mechanical properties properly designed according to the requirements in practical applications,functionally graded materials(FGMs)play an irreplaceable role in improving the performance of the composite materials.The gradient distribution of the thermal and mechanical properties in composites,induced by the change in the distribution of the reinforcement phase,endues FGMs with excellent adaptability to severe loading and environment conditions[Lee(1994)].

    The development achieved in recent years implies the great potential of FGMs in a wide range of thermal,biomedical and structural applications[Sureshet al.(1998);Watariet al.(2004);Reddy(2011);Qiuet al.(1999);Taniet al.(2001)],and the plate structures made of FGMs have also been widely used in practical engineering.Great progress has been made in the research in the mechanical properties of Functionally graded plates(FGPs),such as buckling,bending and vibration[Feldomanet al.(1997);Na and Kim(2006);Aghabaei and Reddy(2009);Altenbach and Eremeyev(2009);Wu and Huang(2009)],crack and thermal fracture[Wanget al(2000);Guoet al.(2005)],transient thermal stress[Chenget al.(2000);Velet al.(2002,2003);Andrewet al.(2010)],static and dynamic responses[Zhonget al.(2003);Elishakoffet al.(2005);Gilhooleyet al.(2007)],shakedown analysis[Penget al.(2009a,b)],etc.

    On the other hand,the concept of FGMs also provides possibilities for optimizing the microstructure parameters of composites to achieve high performance and better utilization of materials[Noda(1990)].The design of the distributions of the thermal and mechanical properties in FGMs has always been an important topic of research since the concept of FGMs was proposed.It is largely due to the excellent and unique advantages of FGMs,for instance,the combination of the advantages of different materials,the compatibilities between different materials and the good adaptability to various environments[Parashkevova(2004)].In the past decade,considerable attention has been devoted to FGM representations(Computer-Aided Design),design validation(Computer-Aided Engineering),fabrication(Computer-Aided Manufacturing)and material heterogeneity optimization[Kou(2012)].

    In the development of modern structures,optimization is one of the most essential topics in the development of FGMs.Markworth and Saunders(1995)optimized the ceramic/metal composition with certain constraints to maximize or minimize the heat fl ux through the materials,where the normal thermal stress pro files were calculated and some unusual behavior was found.Ootaoet al.(1999)optimized the composition in an inhomogeneous hollow sphere with arbitrarily distributed and continuously varying material properties with a neural network approach.Cho and Ha(2002)optimized the volume fraction of Al2O3particles in a Ni/Al2O3composite to minimize the steady state thermal stress,in which the interiorpenalty-function method and the golden section method were employed,together with finite differential method for the sensitivity analysis and an appropriate material property estimation for calculating the thermo-mechanical properties in the graded layer.Chen and Tong(2005)presented a systematic numerical technique to perform sensitivity analysis of coupled thermo-mechanical problems.General formulations were presented based on finite element model by making use of the direct and the adjoint methods.

    The initial progress in the analysis of the shakedown of the FG Bree plate,which is subjected to the coupled constant mechanical loading and cyclically varying temperature loading,has been reported in the authors’previous papers[Penget al.(2009a,b);Zhenget al.(2012)].Different from our previous work,in this paper we focus on the optimal distribution of the reinforcement particles in the FG Bree plate,for the purpose to enhance the capability of the FG plate to bear cyclic thermal loading.In order to achieve a more accurate result,we characterize the distribution of the material properties in the thickness of a FG Bree plate with a piecewise exponential distribution[Guoet al.(2007)],and the effective mechanical property of a material element of the FG plate is evaluated with a double-inclusion mean field approach[Juet al.(1994)].Recent developments of"Computational Grains"also enable a direct simulation of a large number of inclusions for composite or FGM,without FEM meshing of inclusions/matrix[Dong,et al.(2012a,b;2013)].In order to achieve the best shakedown capability of the FG plate,the distribution of reinforcement particles is optimized with the Genetic Algorithm(GA)[Cavalcantiet al.(1997);Chenet al.(2000);Choet al,(2004);Huanget al.(2002);Khalilet al.(2004);Kouet al.(2006,2009);Praveenet al.(1998);Wadleyet al.(2003);Williamset al.(2005)].The optimization model for the shakedown of the FG plate is programmed in MATLAB.Two numerical examples are presented to demonstrate the validity of the proposed approach.

    2 Constitutive model,static and kinematic shakedown theorems

    Assuming small deformation,for initially isotropic and plastically incompressible materials,the constitutive model adopted can be expressed as

    where εijis strain,are its elastic,plastic and thermal components,respectively,which are determined with

    whereE,ν and α are the Young’s modulus,the Poisson’s ratio and the coefficient of linear expansion,respectively,σijandsijare stress and its deviatoric component,θ andθ0are temperature and reference temperature,respectively,δijis the Kronecker delta,and[Penget al.,(1996)]

    s0yis a material constant related to initial yield andf(λ)>0 is a function describing the hardening of material.It can be seen thatdλ is non-negative in any plastic deformation process.The following hardening function is adopted in analysis,

    It can be seen in Eqs.(3)and(4)thatf(λ)increases with the development of plastic deformation and tends to an asymptotic valuedcorresponding to the ultimate strength σy=d·σ0yas plastic deformation fully develops when λ →∞,indicating a saturated state of hardening.Substituting Eq.(3)into Eq.(2)3yields the following Mises-type yield condition

    Introducing the following loading function

    It can be seen that any state of stress should satisfyF(sij,k)≤0,andF(sij,k)=0 de fi nes a loading surface.

    Given a set of actual stress and hardening states,sijandk,and a set of allowable stress and hardening states,s?ij,k?,which satisfy

    And making use of the following inequality

    One can prove that[Penget al.(1993)]

    2.1 Static Shakedown Theorem[Peng et al.(2009a)]

    If there exist a time-independent residual stress fieldρˉijand a time-independent fieldk?such that for all the load variations within a given load domain ?,the following condition holds

    then the total energy dissipated in any allowable load path is bounded.

    In Ineq.(10),sEijis the purely elastic solution of the deviatoric stress determined by external loads and 1≤k?≤d.

    2.2 Kinematic Shakedown Theorem[Peng et al.(2009a)]

    If there exist over a certain time interval(t1,t2),a history of load resulting in a history of purely elastic stresssEij(x,t),and a history of plastic straineˉij(x,t)resulting in a kinematic ally admissible increment such that

    with?ˉui=0 onSu(the boundary where displacement is prescribed),and if shakedown occurs in the given structure,the condition as Ineq.(12)should be satisfied for all kinematically admissible plastic strain cycles.

    is a dissipation function.

    The following relationship can be obtained for practical application by substituting Eq.(2)3and Eq.(3)into Ineq.(12)by following the definition by K?nig(1987),

    and a set of inequalities βk-≤βk≤βk+(k=1,2,...,m)de fi nes the domain ? of loads.

    3 Optimal model

    3.1 The Bree plate

    A plate(Figure 1)of thicknesshis subjected to loads(Px,Py)per unit length in two mutually orthogonal directions.The surfaces of the plate are subjected to temperatures θ2and θ1which vary cyclically,as shown in Figure 1.The cycle time ?tis assumed large compared with characteristic heat conduction time,and the change,between θ0(a reference temperature)and θ0+?ˉθ,is assumed to take place sufficiently slowly for steady state conditions of prevail.The strain εxand εyare assumed to be uniform throughout the thickness of the plate.This problem is a simulation of the behavior of a thin walled tube,in the context of a nuclear fuel can design problem by Bree(1967)for homogeneous material and perfect plasticity.

    Figure 1:The Bree plate.

    Figure 2:Piecewise-exponential distribution model for FG plate with arbitrarily distributed properties.

    3.2 PWED model and material properties

    The piece-wise exponential distribution model(PWED model)[Guoet al.(2007)]is adopted for the distribution of the material properties in an FG plate,as shown in Figure 2.The plate is assumed to be divided intoLparallel layers in the direction of thickness,in each of which the thermal-mechanical properties vary exponentially in the direction of thickness.By this way,the actual thermal-material properties of the plate can be approximated with a set of exponential functions.At both surfaces of a layer the properties are identical with the actual properties of the material.Therefore,the properties of an FG plate can be approximated with sufficient accuracy if the thickness of each layer is sufficiently small.

    The plate of thicknesshis assumed to be divided intoLparallel layers,each of which is marked with subscripti(i=1,2,...,L)counting from the bottom surface of the plate,thei-th layer is located betweenz=hi-1andz=hiand at the bottom surfaceh0=-h/2,while at the top surfacehL=h/2.SupposefM(z)is the real distribution of a thermal/mechanical property in a layer,we can approximatefM(z)withM(z),i.e.,

    Given a set of actual values of a thermal/mechanical property,fM(hi),i=1,2,...,L,for each segmentAMiandBMican be solved from Eqs.(17a,b)as

    whereMcan be replaced respectively with the Young’s modulusE,the coefficient of thermal expansion α ,thermal conductivity λ ,and yield stress σ0y,etc.,therefore,

    Considering that the variation of the Poisson’s ratio in each layer is insignificant and in order to simplify the analysis,we assume the Poisson’s ratio in each layer is a constant taking the mean of the Poisson’s ratio over the thickness,i.e.,

    If the volume fraction of particles athiis ξ (hi),thenfE(hi),fα(hi)andfλ(hi)can be obtained as[Juet al.(1994);Shen(1998)]

    Here the subscript“m”and “c”denote matrix and inclusion,respectively;Km,KcandKare the bulk moduli of the matrix,the inclusion,and the composite,respectively;μm,μcand μ are shear moduli of the matrix,the inclusion,and the composite,respectively.

    Combining the mean field scheme[Juet al.(1994)]and the constitutive model mentioned in section 2,one can obtainfy(hi).

    3.3 Distribution of temperature

    Assuming steady-state heat transfer,the distribution of the temperature in a structure can be described with the following equation of heat conduction,without considering source heat,

    For uniaxial heat transfer,as shown in Figure 3,Eq.(32)can be simplified as

    Figure 3:Heat tran ree plate

    Making the boundary condition shown in Figure 3,we obtain

    Letting??θ=θ1-θ2,?θ(z)=θ(z)-θ2,Eq.(34)can rewritten as

    For the distribution of the temperature in the Bree plate shown in Figure 1,

    Substituting Eq.(21)into Eq.(35a),there have been

    And then?θi(z)can be obtained as

    Here

    3.4 Purely elastic solution

    For the FG Bree plate,we will be concerned with the solution when the displacement iny-direction is fixed,i.e.,εy=0,for simplicity.In this case,we need not be concerned with the yield condition iny-direction[Bree(1967)].The components of the strain in the plate subjected to stress σx,σycan be expressed as

    where ?θi(z)is the temperature change atz.Keeping in mind that εy=0,it can be solved from Eq.(38)that

    (1)Purely elastic solution of the stress distribution σPcaused byPx

    WhenPxis applied individually,i.e.,?θ(z)=0,the equilibrium inx-direction gives

    Substituting Eq.(39)and?θ(z)=0 into Eq.(40a),one obtains

    Keeping in mind that ν is treated as a constant(Eq.(23))and using Eq.(39),we obtain the purely elastic solution of the stress distribution σpcaused byPxas

    Making use of the PWED model,the denominator on the RHS of Eq.(41a)can further be expressed as

    Thus,Eq.(41a)can be expressed in the following discrete form

    (2)Purely elastic solution of the stress distribution σθcaused by?ˉθ

    When?ˉθ is applied individually,one has the following condition of equilibrium

    Keeping in mind that

    And substituting Eq.(39)into Eq.(43)yields

    Combining Eq.(45)with Eq.(39),the distribution of the thermal stress can be determined as

    Making use of the PWED model,we obtain the following discrete form

    Substituting Eqs.(19)-(21)and Eq.(37a)into Eq.(46b),we can obtain

    Compared with Eq.(43),there isHere,the fi rst term in the square bracket on the right side of Eq.(47b)can be expressed as

    3.5 Static shakedown analysis

    The analysis for the shakedown of the Bree plate includes the determination of the boundary of initial yield,the boundary between the areas of shakedown and incremental collapse,and the boundary between the areas of shakedown and reversed plasticity.

    The plate contains two parts,in one partfi(z)≤0 and in the other partfi(z)≥0.For the thermal and mechanical loading shown in Figure 1,we can obtain the following shakedown boundaries.

    (1)Initial yield

    No plastic deformation takes place on the condition that

    (2)Static shakedown boundaries

    Suppose there are a time-independent residual stress field ρˉxi(z)and a time-independentfieldk?(1 ≤k?≤d),the shakedown will occur to the plate if the following condition is satisfied in one temperature change cycle

    where σyi(z)=k*σ0yi(z)andk*=d.

    The shakedown boundaries of the plate contain two parts:the boundary between the area of shakedown and that of incremental collapse,and the boundary between the area of shakedown and that of reversed plasticity.

    (a)Shakedown boundary corresponding to reversed plasticity

    It can be obtained fromIneq.(49)that

    The equality of both sides of anyone in Ineq.(50)at any point in the cross section implies the equality of both sides of Ineqs.(49a)and(49c),or the equality of Ineqs.(49b)and(49d),indicating that reversed plasticity occurs at this point in the cross section.

    (b)Shakedown boundary corresponding to incremental collapse

    Making use of the given time-independent residual stress fieldˉρxi(z)and noticing that??θ=?ˉθ asn?t≤t≤(n+0.5)?t,Ineqs.(49a,b)can be rewritten in this duration as

    Since??θ=0 as(n+0.5)?t≤t≤(n+1)?t,in this duration Ineqs.(49c,d)will be reduced to

    Ineq.(51)indicates that at each point in the region wherefi(z)≥0 of the cross section,the stress σ+=σPi(z)+σθi(z,?ˉθ)+ˉρxi(z)reaches σyi(z);while at each point in the region wherefi(z)≤ 0 of the cross section,the stress σ-= σPi(z)+ˉρxi(z)reaches σyi(z).That is,duringn?t≤t≤(n+0.5)?t,σ+in a part of the Bree plate reaches σyi(z),and during(n+0.5)?t≤t≤ (n+1)?t,σ-at the in the rest part of the Bree plate reaches σyi(z).The two parts of the cross section may flow forward alternatively.

    3.6 Objective of optimization

    To theorize objective function is the key in GA approach.In this paper the objective of the optimization is the maximization of?ˉθ in the FG Bree plate under the condition of shakedown.The temperature?ˉθ can be solved from Ineqs.(50a,b)and Eq.(47a)as

    The objective of the optimization is to find?ˉθmax.

    3.7 Variables and constraint of optimization

    The objective to be optimized is the distribution of volume fraction of the reinforcement particles in a FG Bree plate.The distribution function is assumed as

    whereh=60mm is the thickness of the plate,and the variables to be optimized are

    a1,a2,P1andP2.The boundaries of these variables are assumed as

    Opt_0.3:a1,a2∈[0.0,0.3];P1,P2∈[0.0,2]

    Opt_0.4:a1,a2∈[0.0,0.4];P1,P2∈[0.0,3]

    Opt_0.5:a1,a2∈[0.0,0.5];P1,P2∈[0.0,3]

    The constraint to the optimization is that the gross volume fraction of the reinforcement particles in the plate is 15%,i.e.,

    The shakedown of the plate is determined with Ineqs.(48)-(51),and the best distribution of the volume fraction of the reinforcement particles is obtained by optimization calculation.

    3.8 Computing process of GA approach

    In this paper,the fi rst important thing is to program subroutines in matlab language about the constraint function and the fitting function which is related to the objective function.The computer process of the optimization is depicted in Figure 4.

    The parameters in GA approach are given as follows.The population type is double vector and the population size is 100.The numerical value of generations is 100.The crossover fraction is 0.8 and mutation fraction is 0.01.The constraint function is solved by penalty function method and the penalty factor is 100.

    Figure 4:Computer process of the optimization.

    4 Numerical examples

    4.1 Optimal design of Al/SiC FG Bree plate

    The first example is the optimization of an Al/SiC FG Bree plate.The reinforcement phase is SiC,and the total volume fraction of the SiC particles keeps 15%in the optimization process.The material properties of Al and SiC[Shen(1998)]are listed in Table 1.

    Table 1:Material properties of Al and SiC.

    The optimal design results of the Al/SiC FG plate are

    The shakedown boundaries of Al/SiC FG Bree plates,with the volume fraction distribution functions of the SiC particles given above,are shown in Figure 4 with the curves marked with curve Opt_0.3,curve Opt_0.4 and curve Opt_0.5,respectively.For comparison,the shakedown boundaries of the Al/SiC FG plates with the volume fraction distribution functions given in Eqs.(59)-(61)are also shown in Figure 4,respectively,denoted as Linear 1,Linear 2 and Homo_0.15.

    The distributions of the stress in the plates corresponding to points E1through E6are shown in Figures 6-11,respectively,whereˉρxdenotes the residual stress,σθthe thermal stress,σPthe mechanical stress, σythe yield strength,σ+the stress σθ+σP+ˉρxand σ-the stress σP+ˉρxin the plates.

    Figure 5:Shakedown area of different Al/SiC FG Bree plate after optimization.

    Figure 6:Stress distributions in Al/SiC FG plate corresponding to point E1.

    Figures 6 and 7 show respectively the stress distributions in the Al/SiC FG plates corresponding to the points E1(withVcdetermined by Eq.(59))and E2(withVcdetermined by Eq.(60))in Figure 5.It can be seen in Figure 5 that the?ˉθ at E2is larger than that at E1.The comparison between the stress distributions shown in Figures 6 and 7 indicates that the difference in?ˉθ can be attributed to the difference between the distributions of stress,and the difference between the material properties,for example,the yield strength.One can find in Figure 6 that,the weakest point where both σ+and σ-reach the yield strength,is located at the lower surface of the plate,where the yield strength of the material is the lowest.While in Figure 7 the weakest point, where bothσ+andσ-reach the yield strength,is located at the upper surface of the plate,where the yield strength of the material is the lowest.

    Figure 7:Stress distributions in Al/SiC FG plate corresponding to point E2.

    Figure 8:Stress distributions in Al/SiC effective homogeneous plate corresponding to point E3.

    The stress distributions in the plate,corresponding to the point E3on the curve marked with Homo_0.15(Figure 5)are shown in Figure 8,in which the SiC particles distribute uniformly withVc=0.15(Eq.(61)).It can be seen that for the residual stress fi eld shown in Figure 8(a),the obtained σ+and σ-reach the yield strength at both the upper and the lower surfaces.The corresponding?ˉθ(Figure 5)is a little larger than that corresponding to E2.The comparison between the results corresponding to points E1,E2and E3implies the possibility to achieve a better shakedown boundary by optimizing the distribution of the material properties with a proper distribution of particle volume fraction.

    The optimized shakedown boundary corresponding to the distribution of Eq.(56)is given in Figure 5 with the curve marked with Opt_0.3.Given the residual stress field shown in Figure 9(a),the stress distributions at point E4are shown in Figure 9(b),where it can be seen that reversed plastic deformation takes place at both the upper and the lower surfaces.Compared with the result at either E1,E2or E3,the?θˉ at E4of Opt_0.3 increases to some extent.

    Figure 9:Stress distributions in Al/SiC FG plate corresponding to point E4.

    The optimized shakedown boundary corresponding to the distribution of Eq.(57)is given in Figure 5 with the curve marked with Opt_0.4.Given the residual stress field shown in Figure 10(a),the stress distributions at point E5are shown in Figure 10(b),where it can be seen that reversed plastic deformation takes place at both the upper and the lower surfaces.Compared with the result at E4,the?θˉ at E5in the curve marked with Opt_0.4 increases,which could be attributed to that the maximum particle volume fraction at the upper surface is increased from 0.3 to 0.4,which enhances the capability of the plate to bear the thermal loading.

    The optimized shakedown boundary corresponding to the distribution of the particle volume fraction of Eq.(58)is given in Figure 5 with the curve marked with Opt_0.5.Given the residual stress field shown in Figure 11(a),the stress distributions at point E6are shown in Figure 11(b),where it can be seen that reversed plastic deformation takes place at both two surface(one is the upper surface,the other one is not the lower surface but a little higher than it).Compared with the result at E5,the?ˉθ at E6in the curve marked with Opt_0.5 further increases,which could be attributed to that the maximum particle volume fraction at the upper surface is increased from 0.4 to 0.5,which further enhances the capability of the plate to bear the thermal loading.

    The comparison between the?ˉθ at the points E1through E6in Figure 5 shows that the?ˉθ at E6,corresponding to the distribution of the particle volume fraction,Eq.(58),is the largest.If the maximum particle volume fraction is limited to 0.5,this?ˉθ can be regarded as?ˉθmax.

    Figure 10:Stress distributions in Al/SiC FG plate corresponding to point E5.

    Figure 11:Stress distributions in Al/SiC FG plate corresponding to point E6.

    4.2 Optimal design of Ti/Si3N4FG Bree plate

    In order to further verify the optimal design model,the optimization for shakedown capability of a Ti/Si3N4FG Bree plate is performed,in which the Si3N4particles are the reinforcement phase.Material constants of Ti and Si3N4[Cho et al.(2004)]are listed in Table 2.

    We also assume three classes of the Ti/Si3N4FG Bree plates,of which the distribution of the Si3N4particle volume function is described with Eq.(54),wherea1,P1anda2,P2are the parameters to be optimized.If the average particle volume is limited to 0.15,and the maximum local particle volume fractions of the three classes are 0.3(Class 1),0.4(Class 2)and 0.5(Class 3),respectively,the corresponding shakedown boundaries are shown in Figure 12 with the curves marked with Opt_0.3,o Opt_0.4 and Opt_0.5,respectively.With the optimizeda1,P1anda2,P2,the three classes distribution functions can be expressed as

    Table 2:Material constants of Ti and Si3N4.

    For comparison,the shakedown boundaries of Ti/Si3N4FG plates with the particle volume fraction distribution functions of Eqs.(59)-(61)are also given in Figure 12 with the curves marked with Linear 1,Linear 2 and Homo_0.15 respectively.

    Figure 13 show the stress distributions in the FG plate corresponding to point E1on the curve marked as Linear 1 in Figure 13.Figure 14 shows the stress distributions in the FG plate corresponding to point E2on the curve marked as Linear 2 in Figure 11.The results shown in Figure 13 and Figure 14 are respectively similar to that in Figure 6 and Figure 7.

    Then,the shakedown of the equivalent homogenous plate is analyzed in the same way as numerical example 1,and the boundary which is marked as Homo_0.15 is shown in Figure 12.Stress distributions corresponding to point E3are shown in Figure 15.

    Be similar to the numerical example of Al/SiCFG plate,the optimization of Ti/Si3N4FG plate is given.The results of optimal volume fraction distribution functions of Si3N4have been given as Eqs.(62)-(64),and the shakedown boundaries of corresponding FG plates have been shown as Opt_0.3,Opt_0.4 and Opt_0.5 in Figure 12.As shown in Figure 12,the plate of which the reinforcement phase volume fraction function is optimized can be endured larger?ˉθmax,i.e.,?ˉθmaxof whichever Opt_0.3,Opt_0.4 and Opt_0.5 is larger than that of Homo_0.15.

    Figure 12:Shakedown area of different Ti/Si3N4FG Bree plate after optimization.

    Figure 13:Stress distributions in Ti/Si3N4FG plate corresponding to point E1.

    Figure 14:Stress distributions in Ti/Si3N4FG plate corresponding to point E2.

    Figure 15:Stress distributions in Ti/Si3N4effective homogenous plate corresponding to point E3.

    Figure 17:Stress distributions in Ti/Si3N4FG plate corresponding to point E5.

    Figure 18:Stress distributions in Ti/Si3N4FG plate corresponding to point E6.

    The stress distributions corresponding to point E4on Opt_0.3,point E5on Opt_0.4 and point E6on Opt_0.5 have been respectively shown in Figures 16-18.In these figures the phenomena are similar with those in Figures 9-11,although the constituents of the FG plates are different.

    5 Conclusion and discussion

    The optimization of the volume fraction distribution function of reinforcement particles in FG Bree plates was performed with the GA approach.The objective of optimization was to gain the maximum of the?ˉθ applied to the plates under the constraint that the average volume fraction of the reinforcement particles keeps constant.An optimal model was proposed and two numerical examples were presented,from which the following conclusion may be drawn:

    1.An approach was developed for the optimization of the distribution of the reinforcement particles in the plate for the purpose to enhance the shakedown capability of the plate.The distribution of the material properties in the thickness of the plate was characterized with a piecewise exponential distribution,and the effective mechanical property of the material of the plate was evaluated with a double-inclusion mean fi eld scheme,and the distribution of the volume fraction of the reinforcement particles was optimized with the Genetic Algorithm.The validity of the approach was demonstrated by the two numerical examples.

    2.The shakedown capability of a FG plate is determined by both the applied thermal-mechanical loading and the thermal-mechanical properties of the plate.For a FG plate,the latter can be designed by properly distributing the volume fraction of the reinforcement particles,which provides the possibility to achieve desired or optimal thermal-mechanical properties of the plate to substantially enhance the shakedown capability of the plate.It was found in this paper that the distribution of the particle volume fraction affects markedly the shakedown capability of the FG plate,and a proper distribution of the particle volume fraction can substantially enhance the shakedown capability of the plate.

    3.The concept of the piecewise distribution of the material properties is in accordancewith that of the piecewise approach used in the shakedown analysis;therefore,it is of particular advantage for the optimization of the shakedown capability of the FG plate,because it may bring great convenience to the optimization without introducing additional difficulties.

    Acknowledgement:The authors gratefully acknowledge the financial support to this work from NSFC under Grant Numbers 11272364 and 11332013.

    Andrew,J.G.;Senthil,S.V.(2010):Transient multiscale thermoelastic analysis of functionally graded materials.Compos.Struct.,vol.92,pp.1372-1390.

    Aghababaei,R.;Reddy,J.N.(2009):Nonlocal third-order shear deformation plate theory with application to bending and vibration of plates.J.Sound Vib.,vol.326,pp.277–289.

    Altenbach,H.;Eremeyev,V.A.(2009):Eigen-vibrations of Plates made of Functionally Graded Material.CMC:Computers,Materials,&Continua,vol.9,pp.153-178.

    Bree,J.(1967):Elastic plastic behavior of thin tubes subjected to internal pressure and intermittent high heat fluxes with application to fast nuclear reactor fuel elements.J.Strain Anal.,vol.2,pp.226-238.

    Cavalcanti,P.R;Carvalho,P.C.P;Martha,L.F.(1997):Non-manifold modeling:an approach based on spatial subdivision.Computer-Aided Design,vol.29,pp.209-220.

    Chen,B.;Tong,L.(2005):Thermo-mechanically coupled sensitivity analysis and design optimization of functionally graded materials.Comput.Methods Appl.Mech.Eng.,vol.194,pp.1891-1911.

    Chen,M.;Tucker,J.V.(2000):Constructive volume geometry.Comput.Graph.Forum,vol.19,pp.281-293.

    Cheng,Z.Q.;Batra,R.C.(2000):Three-dimensional thermoelastic deformations of a functionally graded elliptic plate.Composites Part B,vol.31,pp.97-106.

    Cho,J.R.;Choi,J.H.(2004):A yield-criteria tailoring of the volume fraction in metal-ceramic functionally graded material.Euro.J.Mech./A Solids,vol.23,pp.271-281.

    Cho,J.R.;Ha,D.Y.(2002):Volume fraction optimization for minimizing thermal stress in Ni–Al2O3functionally graded materials.Mater.Sci.Eng.A,vol.334,pp.147–155.

    Cho,J.R.;Shin,S.W.(2004):Material composition optimization for heat-resistingFGMs by artificial neural network.Composites Part A:Applied Science and Manufacturing,vol.35,pp.585–594.

    Dong,L.;Atluri,S.N.(2012a):T-Trefftz voronoi cell fi nite elements with elastic/rigid inclusions or voids for micro mechanical analysis of composite and porous materials.CMES:Computers,Materials,&Continua,vol.83,pp.183-219.

    Dong,L.;Atluri,S.N.(2012b):Development of 3D Trefftz voronoi cells with ellipsoidal voids&/or elastic/rigid inclusions for micro mechanical modeling of heterogeneous materials.CMC:Computers Materials and Continua,vol.30,pp.39-81.

    Dong,L.;Atluri,S.N.(2013):SGBEM voronoi cells(SVCs),with embedded arbitrary-shaped inclusions,voids,and/or cracks,for micro mechanical modeling of heterogeneous materials.CMC:Computers,Materials&Continua,vol.33,pp.111-154.

    Elishakoff,I.;Gentilini,C.(2005):Three-dimensional flexure of rectangular plates made of functionally graded materials.J.Appl.Mech.,vol.72,pp.788-791.

    Feldoman,E.;Aboudi,J.(1997):Buckling analysis of functionally graded plates subjected to uniaxial loading.Compos.Struct.,vol.38,pp.29-36.

    Guo,L.C.;Noda,N.(2007):Modeling method for a crack problem of functionally graded materials with arbitrary properties piecewise exponential model.Int.J.Solids Struct.,vol.44,pp.6768-6790.

    Guo,L.;Wu,L.;Sun,Y.;Ma,L.(2005):The transient fracture behavior for a functionally graded layered structure subjected to an in-plane impact load.Acta Mech.Sinica,vol.21,pp.257-266.

    Gilhooley,D.F.;Batra,R.C.;Xiao,J.R.;Mccarthy,M.A.;Gillespie,J.W.(2007):Analysis of thick functionally graded plates by using higher-order shear and normal deformable plate theory and MLPG method with radial basis functions.Compos.Struct.,vol.80,pp.39-552.

    Huang,J.;Fadel,G.M;Blouin,V.Y;Grujicic,M.(2002):Bi-objective optimization design of functionally gradient materials.Materials and Design,vol.23,pp.657-666.

    Ju,J.W.;Chen,T.M.(1994):Micro mechanics and effective moduli of elastic composites containing randomly dispersed ellipsoidal inhomogeneities.Acta Mech.,vol.103,pp.103-121.

    Khalil,S.;Nam,J.;Darling,A.;Sun,W.(2004):Multi-nozzle biopolymer deposition for freeform fabrication of tissue construct.In:Proc.15th solid freeform fabrication symposium,pp.826-837.

    Kou,X.Y.;Parks,G.T.;Tan,S.T.(2012):Optimal design of functionally graded materials using a procedural model and particle swarm optimization.Computer-Aided Design,vol.44,pp.300-310.

    Kou,X.Y.(2006):Computer-aided design of heterogeneous objects.Ph.D.thesis.University of Hong Kong.

    Kou,X.Y.;Tan,S.T.;Sze,W.S.(2006):Modeling complex heterogeneous objects with nonmanifold heterogeneous cells.Computer-Aided Design,vol.38,pp.457-474.

    Kou,X.Y;Tan,S.T.(2009):Robust and efficient algorithms for rapid prototyping of heterogeneous objects.Rapid Prototyping J.,vol.15,pp.5-18.

    K?nig,J.A.(1987):Shakedown of elastic-plastic Structures.PWN-Polish Scienti fi c Publishers.

    Lee,Y.D.;Erdogan,F.(1994):Residual/thermal stress in FGM and laminated thermal barrier coating.Int.J.Fracture,vol.69,pp.145-165.

    Na,K.S.;Kim,J.H.(2006):Three-dimensional theCompos.Struct.,vol.73,pp.413-422.

    Noda,N.;Tsuji,T.(1990):Steady thermal stresses in a plate of functionally gradient material.Proc.1st Int.Symp.Functionally Gradient Materials,pp.339–344.Mark worth,A.J.;Saunders,J.H.(1995):A model of structure optimization for functionally graded materials.Mater.Lett.,vol.22,pp.103–107.

    Ootao,Y.;Kawamura,R.;Tanigawa,Y.;Imamura,R.(1999):Optimization of material composition of nonhomogeneous hollow sphere for thermal stress relaxation making use of neural network.Comput.Methods Appl.Mech.Eng,vol.180,pp.185-201.

    Parashkevova,L.;Ivanova,J.;Bontcheva,N.(2004):Optimal design of func-tionally graded plates with thermo-elastic-plastic behavior.Comptes Rendus Mecaniquvol.332,pp.493-498.

    Peng,X.;Fan,J.;Zeng,X.(1996):Analysis for plastic buckling of thin-walled cylinders via non-classical constitutive theory of plasticity.In.J.Solids and Struct.,vol.33,pp.4495-4509.

    Peng,X.;Hu,N.;Zheng,H.;Fang,C.(2009a):Analysis of shakedown of FG Bree plate subjected to coupled thermal-mechanical loadings.Acta Mech.Solida Sinica,vol.22,pp.95-108.

    Peng,X.;Ponter,A.R.S.(1993):Extremal properties of Endochronic plasticity,Part II:Extremal path of the Endochronic constitutive equation with a yield surface and application.Int.J.Plasticity,vol.9,pp.567-581.

    Peng,X.;Zheng,H.;Hu,N.;Fang,C.(2009b):Static and kinematic shakedown analysis of FG plate subjected to constant mechanical load and cyclically varying temperature change.Composite Structures,vol.91,pp.212-221.

    Praveen,G.N.;Reddy,J.N.(1998):Nonlinear transient thermoelastic analysis of functionally graded ceramic-metal plates.Int.J.Solids and Struct.,vol.35,pp.4457-4476.

    Qiu,J.;Tani,J.;Warkentin,D.J.(1999):Stress analysis of RAINBOW actuators and relief of stress by gradation of material composition.Japan.Soc.Appl.Electromagn.Mech.,vol.7,pp.185192

    Reddy,J.N.(2011):Microstructure-dependent couple stress theories of functionally graded beams.J.Mech.Phys.Solids,vol.59,pp.2382-2399.

    Shen,Y.(1998):Thermal expansion of metal–ceramic composites: a three-dimensionalanalysisMater.Sci.Eng.A,vol.252,pp.269-275.

    Suresh,S.;Mortensen,A.(1998):Fundamentals of functionally graded materials:processing and thermomechanical behaviour of graded metals and metal-ceramic composites.IOM Communications Ltd,London.

    Suresh,S.;Mortensen,A.(1998):Fundamentals of functionally graded materials.London,UK:Institute of Materials.

    Tani,J.;Qiu,J.;Morit,A.T.(2001):Functionally gradient piezoelectric actuatorsTrans.Mater.Res.Soc.Japan,vol.26,pp.283286

    Vel,S.S.;Batra,R.C.(2002):Exact solution for thermoelastic deformations of functionally graded thick rectangular plates.AIAA J.,vol.40,pp.1421-1433.

    Vel,S.S.;Batra,R.C.(2003):Three-dimensional analysis of transient thermal stresses in functionally graded plates.Int.J.Solids.Struct.,vol.40,pp.7181-7196.

    Wadley,H.N.G.;Fleck,N.A.;Evans,A.G.(2003):Fabrication and structural performance of periodic cellular metal sandwich structures.Compos.Sci.Technol.,vol.63,pp.2331-2343.

    Wang,B.L.;Han,J.C.;Du,S.Y.(2000):Crack problems for functionally graded materials under transient thermal loading.J.Thermal Stresses,vol.23,pp.143-168.

    Watari,F.;Yokoyama,A.;Omori,M.;Hirai,T.;Kondo,H.;Uo,M.;Kawasaki,T.(2004):Biocompatibility of materials and development to functionally graded implant for biomedical application.Compos.Sci.Technol.,vol.64,pp.893-908.

    Williams,C.B.;Mistree,F.;Rosen,D.W.(2005):Investigation of additive manufacturing processes for the manufacture of parts with designed mesostructure.In:Proc.ASME 2005 design engineering technical conferences,pp.353-364.

    Wu,C.P.;Huang,S.E.(2009):Three-Dimensional Solutions of Functionally Graded Piezo-Thermo-Elastic Shells and Plates Using a Modi fi ed Pagano Method.CMC:Computers,Materials,&Continua,vol.12,pp.251-282.

    Zheng,H.;Peng,X.;Hu,N.(2012):Analysis for Shakedown of Functionally Graded Plate Subjected to Thermal-Mechanical Loading with Piecewise-ExponentialDistribution of Material Properties.CMES:Comp.Model.Eng.Sci.,vol.86,pp.505-531.

    Zhong,Z.;Shang,E.T.(2003):Three-dimensional exact analysis of simply supported functionally gradient piezoelectric plates.Int.J.Solids Struct.,vol.40,pp.5335-5352.

    亚洲伊人久久精品综合| 在线a可以看的网站| 国产午夜福利久久久久久| av网站免费在线观看视频| eeuss影院久久| av一本久久久久| av国产久精品久网站免费入址| 亚洲精品第二区| 国产成人a区在线观看| 免费黄频网站在线观看国产| 可以在线观看毛片的网站| 欧美激情久久久久久爽电影| 免费看日本二区| 日本一本二区三区精品| 国产成人精品福利久久| 91久久精品电影网| 美女国产视频在线观看| 嫩草影院新地址| 高清欧美精品videossex| 久热久热在线精品观看| 亚洲国产日韩一区二区| 中文欧美无线码| 亚洲国产欧美人成| 水蜜桃什么品种好| 精品久久久久久久人妻蜜臀av| 美女脱内裤让男人舔精品视频| 精华霜和精华液先用哪个| 国产成人福利小说| 亚洲va在线va天堂va国产| 中文字幕制服av| 亚洲无线观看免费| 午夜福利在线在线| 听说在线观看完整版免费高清| 精品一区二区三区视频在线| 美女被艹到高潮喷水动态| 国产精品久久久久久久电影| 97超视频在线观看视频| 男人和女人高潮做爰伦理| 欧美成人午夜免费资源| 波野结衣二区三区在线| 久久久久精品久久久久真实原创| 色视频在线一区二区三区| 国产中年淑女户外野战色| 亚洲三级黄色毛片| 国产爱豆传媒在线观看| 成人特级av手机在线观看| 好男人在线观看高清免费视频| 精品人妻一区二区三区麻豆| 熟女电影av网| 国产精品人妻久久久影院| 一本久久精品| 精品视频人人做人人爽| 亚洲精品久久久久久婷婷小说| 精品久久久久久久久av| 少妇人妻精品综合一区二区| 尾随美女入室| 亚洲精品日韩av片在线观看| 嫩草影院新地址| 亚洲国产成人一精品久久久| 成人亚洲欧美一区二区av| 国产精品秋霞免费鲁丝片| 亚洲高清免费不卡视频| 久久久久久久亚洲中文字幕| 中文字幕av成人在线电影| av又黄又爽大尺度在线免费看| 高清欧美精品videossex| 国产成人精品一,二区| 成年女人看的毛片在线观看| xxx大片免费视频| 热99国产精品久久久久久7| 亚洲人成网站高清观看| 亚洲一区二区三区欧美精品 | 久久精品久久久久久久性| 免费av观看视频| 2022亚洲国产成人精品| 一区二区三区四区激情视频| 新久久久久国产一级毛片| 精品国产一区二区三区久久久樱花 | 午夜亚洲福利在线播放| 国产精品福利在线免费观看| 18禁裸乳无遮挡动漫免费视频 | 真实男女啪啪啪动态图| 看黄色毛片网站| 亚洲成色77777| 久久精品人妻少妇| 在线观看三级黄色| 一区二区三区免费毛片| 嘟嘟电影网在线观看| 特大巨黑吊av在线直播| 国产免费又黄又爽又色| 国产毛片在线视频| 久久国产乱子免费精品| 久久午夜福利片| 中文字幕人妻熟人妻熟丝袜美| 99视频精品全部免费 在线| 最近最新中文字幕免费大全7| 日本一二三区视频观看| 人人妻人人看人人澡| av在线观看视频网站免费| 99九九线精品视频在线观看视频| 国产在线一区二区三区精| 欧美三级亚洲精品| av一本久久久久| 女人被狂操c到高潮| 国产一区二区三区综合在线观看 | 91在线精品国自产拍蜜月| 成人亚洲精品av一区二区| 女人久久www免费人成看片| 国产精品秋霞免费鲁丝片| 亚洲va在线va天堂va国产| 男人添女人高潮全过程视频| 国产精品成人在线| 六月丁香七月| 联通29元200g的流量卡| 亚洲精品自拍成人| av在线播放精品| 亚洲精品一区蜜桃| 国产免费一级a男人的天堂| 欧美另类一区| 精品99又大又爽又粗少妇毛片| 久久女婷五月综合色啪小说 | 精品99又大又爽又粗少妇毛片| 欧美区成人在线视频| 97超视频在线观看视频| 欧美丝袜亚洲另类| 九草在线视频观看| 大码成人一级视频| 性色av一级| 黄色一级大片看看| 又粗又硬又长又爽又黄的视频| 欧美xxxx性猛交bbbb| 日日摸夜夜添夜夜添av毛片| 久久久成人免费电影| 国产亚洲精品久久久com| 成人欧美大片| 久久精品国产自在天天线| 日韩免费高清中文字幕av| 联通29元200g的流量卡| 日韩精品有码人妻一区| 日本欧美国产在线视频| a级一级毛片免费在线观看| av在线观看视频网站免费| 国产成人免费无遮挡视频| 国产日韩欧美亚洲二区| 男人舔奶头视频| 精品少妇久久久久久888优播| 高清在线视频一区二区三区| 91久久精品电影网| 男男h啪啪无遮挡| 丝袜喷水一区| 91aial.com中文字幕在线观看| 国产免费视频播放在线视频| 欧美潮喷喷水| 26uuu在线亚洲综合色| 午夜老司机福利剧场| 美女被艹到高潮喷水动态| 亚洲欧美日韩无卡精品| 日本欧美国产在线视频| 国产欧美亚洲国产| 久久久久久久午夜电影| 欧美成人午夜免费资源| av在线观看视频网站免费| 大片免费播放器 马上看| 成人欧美大片| 少妇丰满av| 亚洲欧美清纯卡通| 99re6热这里在线精品视频| 久久久久久久久大av| 下体分泌物呈黄色| 成人免费观看视频高清| 精品午夜福利在线看| 一级黄片播放器| 欧美激情久久久久久爽电影| 国产亚洲最大av| 欧美日韩亚洲高清精品| 午夜福利在线观看免费完整高清在| av一本久久久久| 精品99又大又爽又粗少妇毛片| 国产免费一级a男人的天堂| 日韩成人av中文字幕在线观看| 丝袜喷水一区| 少妇丰满av| 麻豆成人午夜福利视频| 国产黄片美女视频| 舔av片在线| 亚洲精品乱码久久久v下载方式| 边亲边吃奶的免费视频| 免费观看性生交大片5| 在线观看人妻少妇| 日本色播在线视频| 熟女人妻精品中文字幕| 大陆偷拍与自拍| 嫩草影院精品99| av线在线观看网站| 国产精品福利在线免费观看| 亚洲国产色片| 国产精品一区二区在线观看99| 日本猛色少妇xxxxx猛交久久| 国产一级毛片在线| 纵有疾风起免费观看全集完整版| 少妇猛男粗大的猛烈进出视频 | 国产黄片视频在线免费观看| 久久精品熟女亚洲av麻豆精品| 亚洲成人中文字幕在线播放| 美女cb高潮喷水在线观看| 免费av观看视频| 日日撸夜夜添| 特大巨黑吊av在线直播| 男的添女的下面高潮视频| 国产高清有码在线观看视频| 国产成人精品一,二区| 少妇人妻一区二区三区视频| 午夜福利视频1000在线观看| 国产毛片a区久久久久| 大码成人一级视频| 尤物成人国产欧美一区二区三区| 国产真实伦视频高清在线观看| 午夜福利视频1000在线观看| 十八禁网站网址无遮挡 | 菩萨蛮人人尽说江南好唐韦庄| 免费av观看视频| 精品久久久久久久久亚洲| 日本欧美国产在线视频| 大陆偷拍与自拍| 18禁在线播放成人免费| 免费少妇av软件| 亚洲熟女精品中文字幕| 一级黄片播放器| 国产精品久久久久久久电影| 久久久久久久精品精品| 欧美变态另类bdsm刘玥| 少妇人妻 视频| 熟女人妻精品中文字幕| 国产免费福利视频在线观看| 王馨瑶露胸无遮挡在线观看| 欧美日本视频| 亚洲自偷自拍三级| 久久精品国产鲁丝片午夜精品| 久久99精品国语久久久| 大片免费播放器 马上看| 另类亚洲欧美激情| 99久久中文字幕三级久久日本| 亚洲精品中文字幕在线视频 | 日韩一区二区视频免费看| 一区二区av电影网| 国产69精品久久久久777片| 精品熟女少妇av免费看| 大片免费播放器 马上看| 日韩成人伦理影院| 国产亚洲最大av| 国产精品麻豆人妻色哟哟久久| 精品人妻视频免费看| 国产老妇伦熟女老妇高清| 久久久久久久久久久免费av| 观看美女的网站| 国产淫语在线视频| 一本色道久久久久久精品综合| 女的被弄到高潮叫床怎么办| 男人舔奶头视频| 久久国内精品自在自线图片| 久久女婷五月综合色啪小说 | 小蜜桃在线观看免费完整版高清| 2018国产大陆天天弄谢| 观看免费一级毛片| 黄片wwwwww| 久久精品夜色国产| 三级国产精品欧美在线观看| 97热精品久久久久久| 色5月婷婷丁香| xxx大片免费视频| 熟妇人妻不卡中文字幕| 日韩欧美精品v在线| 亚洲人与动物交配视频| 久久久久久九九精品二区国产| 午夜老司机福利剧场| 日本熟妇午夜| 青青草视频在线视频观看| 尤物成人国产欧美一区二区三区| 熟女电影av网| eeuss影院久久| av网站免费在线观看视频| 免费电影在线观看免费观看| 亚洲国产精品成人久久小说| 亚洲性久久影院| 国产成人福利小说| 大又大粗又爽又黄少妇毛片口| 少妇人妻久久综合中文| 精品99又大又爽又粗少妇毛片| 国产精品偷伦视频观看了| 久久久久精品性色| 日韩欧美 国产精品| 久久午夜福利片| 看免费成人av毛片| 午夜亚洲福利在线播放| 欧美+日韩+精品| 国产视频首页在线观看| 天堂网av新在线| 久久久久久伊人网av| 一级av片app| 午夜福利在线在线| 国产中年淑女户外野战色| 麻豆久久精品国产亚洲av| 一区二区三区免费毛片| 在线a可以看的网站| 欧美zozozo另类| 日韩成人伦理影院| 国产精品女同一区二区软件| 高清视频免费观看一区二区| 成人漫画全彩无遮挡| 久久久精品免费免费高清| 18禁裸乳无遮挡动漫免费视频 | 久久精品久久久久久噜噜老黄| 久久99热6这里只有精品| 91久久精品电影网| 亚洲国产精品成人久久小说| 欧美精品一区二区大全| 亚洲成人久久爱视频| 成年女人看的毛片在线观看| 一本一本综合久久| 永久网站在线| 国内精品宾馆在线| 国产av码专区亚洲av| 欧美zozozo另类| 不卡视频在线观看欧美| 18禁在线播放成人免费| 国产精品不卡视频一区二区| 男女边吃奶边做爰视频| 秋霞伦理黄片| 国产亚洲91精品色在线| 精品一区二区三卡| 又粗又硬又长又爽又黄的视频| av在线老鸭窝| 97超视频在线观看视频| 国产成人a区在线观看| 日韩大片免费观看网站| 观看免费一级毛片| 亚洲国产欧美在线一区| 黄色欧美视频在线观看| 又粗又硬又长又爽又黄的视频| 日本一本二区三区精品| 又粗又硬又长又爽又黄的视频| 在线观看一区二区三区激情| 黑人高潮一二区| 精品人妻一区二区三区麻豆| 97超视频在线观看视频| 国产日韩欧美在线精品| 自拍欧美九色日韩亚洲蝌蚪91 | 极品少妇高潮喷水抽搐| 日韩精品有码人妻一区| 又爽又黄a免费视频| 精品久久久久久电影网| 亚洲久久久久久中文字幕| 男女下面进入的视频免费午夜| 香蕉精品网在线| 亚洲欧美精品专区久久| 欧美3d第一页| 看非洲黑人一级黄片| 嘟嘟电影网在线观看| 成人黄色视频免费在线看| 成年免费大片在线观看| 91aial.com中文字幕在线观看| 国产69精品久久久久777片| av播播在线观看一区| 欧美丝袜亚洲另类| 亚洲天堂av无毛| 美女高潮的动态| 国产亚洲5aaaaa淫片| 最近2019中文字幕mv第一页| 亚洲精品自拍成人| 女人被狂操c到高潮| 少妇人妻一区二区三区视频| 免费看av在线观看网站| av免费观看日本| 美女视频免费永久观看网站| 日韩不卡一区二区三区视频在线| 亚洲天堂国产精品一区在线| 18禁裸乳无遮挡动漫免费视频 | 亚洲av一区综合| 亚洲国产精品999| 在线观看国产h片| 在线播放无遮挡| 欧美一区二区亚洲| www.av在线官网国产| av在线播放精品| 国产精品.久久久| 国产成人免费无遮挡视频| h日本视频在线播放| 天天躁日日操中文字幕| 亚洲经典国产精华液单| 国模一区二区三区四区视频| 日本黄大片高清| av女优亚洲男人天堂| 中文天堂在线官网| 熟女人妻精品中文字幕| av播播在线观看一区| 91精品国产九色| 男女下面进入的视频免费午夜| 亚洲图色成人| 18禁裸乳无遮挡动漫免费视频 | 精品少妇久久久久久888优播| 亚洲av一区综合| 亚洲最大成人中文| av卡一久久| 波野结衣二区三区在线| 在现免费观看毛片| 国产高潮美女av| 水蜜桃什么品种好| 国产真实伦视频高清在线观看| 亚洲av不卡在线观看| 99热国产这里只有精品6| 久久久久久久久久久丰满| 美女国产视频在线观看| 免费人成在线观看视频色| 精品一区在线观看国产| 欧美+日韩+精品| 国产成人免费无遮挡视频| 老司机影院毛片| 午夜福利视频1000在线观看| 久久久亚洲精品成人影院| 午夜激情福利司机影院| 婷婷色综合大香蕉| 国产亚洲av片在线观看秒播厂| 七月丁香在线播放| 联通29元200g的流量卡| 深爱激情五月婷婷| 人人妻人人澡人人爽人人夜夜| 欧美xxxx黑人xx丫x性爽| 三级经典国产精品| 国产精品嫩草影院av在线观看| 午夜福利网站1000一区二区三区| 久久99热6这里只有精品| eeuss影院久久| 精品人妻一区二区三区麻豆| 麻豆国产97在线/欧美| 欧美日本视频| 男男h啪啪无遮挡| 亚洲国产高清在线一区二区三| 色婷婷久久久亚洲欧美| 亚洲av中文av极速乱| 性色av一级| 视频区图区小说| 亚洲不卡免费看| 国产成年人精品一区二区| 日韩欧美一区视频在线观看 | 国产成人福利小说| 久久久久久久久久成人| 嫩草影院入口| 日本黄大片高清| 秋霞在线观看毛片| 亚洲精品成人久久久久久| 又爽又黄a免费视频| 高清av免费在线| 亚洲人成网站在线观看播放| 精品少妇久久久久久888优播| 国产高清不卡午夜福利| 黄色视频在线播放观看不卡| 免费观看a级毛片全部| 别揉我奶头 嗯啊视频| 18禁裸乳无遮挡动漫免费视频 | 91久久精品国产一区二区成人| 高清午夜精品一区二区三区| 乱系列少妇在线播放| 国产真实伦视频高清在线观看| 亚洲精品乱码久久久v下载方式| 免费观看性生交大片5| 日本熟妇午夜| 美女cb高潮喷水在线观看| 亚洲激情五月婷婷啪啪| 亚洲人成网站在线播| 色吧在线观看| 欧美日本视频| av国产久精品久网站免费入址| 91久久精品电影网| 国产一级毛片在线| 国产欧美另类精品又又久久亚洲欧美| 久久亚洲国产成人精品v| 毛片一级片免费看久久久久| 身体一侧抽搐| 天堂俺去俺来也www色官网| 久久精品久久久久久噜噜老黄| 久久久久久久亚洲中文字幕| 偷拍熟女少妇极品色| 亚洲色图av天堂| 久久久久网色| av国产免费在线观看| 亚洲天堂国产精品一区在线| 嫩草影院精品99| 黄色日韩在线| 亚洲欧美一区二区三区国产| 视频区图区小说| 男女无遮挡免费网站观看| 18禁在线播放成人免费| 熟女电影av网| 久久99热6这里只有精品| 精品一区二区免费观看| 国产精品国产av在线观看| 啦啦啦中文免费视频观看日本| 一级毛片黄色毛片免费观看视频| 青春草国产在线视频| 日韩不卡一区二区三区视频在线| 免费人成在线观看视频色| 亚洲精品久久午夜乱码| 日本黄大片高清| 热re99久久精品国产66热6| 国产精品一及| 国产伦理片在线播放av一区| 久久久久九九精品影院| 国产精品成人在线| 丝袜喷水一区| 视频中文字幕在线观看| 全区人妻精品视频| 人妻系列 视频| 日本wwww免费看| 国产欧美日韩精品一区二区| 亚洲,一卡二卡三卡| 1000部很黄的大片| 亚洲国产成人一精品久久久| 九色成人免费人妻av| 十八禁网站网址无遮挡 | 欧美日韩国产mv在线观看视频 | 激情五月婷婷亚洲| 国产白丝娇喘喷水9色精品| 亚洲最大成人手机在线| 中文欧美无线码| 日本一本二区三区精品| 嫩草影院新地址| 国内精品美女久久久久久| 91精品国产九色| 一级片'在线观看视频| 在线精品无人区一区二区三 | 成人一区二区视频在线观看| 久久精品久久精品一区二区三区| 午夜福利网站1000一区二区三区| 国产欧美日韩一区二区三区在线 | 亚洲国产日韩一区二区| 下体分泌物呈黄色| 亚洲av成人精品一二三区| 日韩一区二区视频免费看| 91精品伊人久久大香线蕉| 久久精品人妻少妇| av黄色大香蕉| 亚洲av一区综合| 制服丝袜香蕉在线| 中文字幕亚洲精品专区| 国产av不卡久久| 久久久久国产网址| 水蜜桃什么品种好| www.色视频.com| 国产又色又爽无遮挡免| 日韩免费高清中文字幕av| 日韩欧美精品v在线| 夜夜爽夜夜爽视频| av在线播放精品| 在线观看av片永久免费下载| 内地一区二区视频在线| 国国产精品蜜臀av免费| 人人妻人人澡人人爽人人夜夜| 亚洲精品成人av观看孕妇| 亚洲伊人久久精品综合| 国产成人一区二区在线| 2021少妇久久久久久久久久久| av免费观看日本| 午夜免费观看性视频| 午夜福利视频精品| 国产极品天堂在线| 亚洲av在线观看美女高潮| 狂野欧美激情性xxxx在线观看| 在线a可以看的网站| 美女高潮的动态| 大片电影免费在线观看免费| 亚洲av二区三区四区| 麻豆乱淫一区二区| 国产欧美另类精品又又久久亚洲欧美| 欧美激情久久久久久爽电影| 亚洲精品第二区| 免费不卡的大黄色大毛片视频在线观看| 国产成人freesex在线| 国产美女午夜福利| 99久久中文字幕三级久久日本| 国产探花极品一区二区| 亚洲最大成人手机在线| 日韩精品有码人妻一区| 国产精品不卡视频一区二区| 男人和女人高潮做爰伦理| 免费播放大片免费观看视频在线观看| kizo精华| 久久久久网色| www.色视频.com| 亚洲在久久综合| 在线免费十八禁| 国产精品av视频在线免费观看| 人人妻人人看人人澡| 亚洲四区av| 亚洲欧洲国产日韩| 日韩在线高清观看一区二区三区| 国产精品熟女久久久久浪| 亚洲欧美成人综合另类久久久| 亚洲国产精品999| 国产精品久久久久久精品电影小说 | 最近最新中文字幕免费大全7| 亚洲国产日韩一区二区| 99热6这里只有精品| 中国三级夫妇交换| 99久久九九国产精品国产免费| 国产乱来视频区| 国产欧美另类精品又又久久亚洲欧美| 白带黄色成豆腐渣| 国产69精品久久久久777片| 有码 亚洲区| 夜夜看夜夜爽夜夜摸| 国产爽快片一区二区三区| 国产精品精品国产色婷婷| 伦理电影大哥的女人| 建设人人有责人人尽责人人享有的 | 久久精品国产自在天天线|