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

    Explicit Isogeometric Topology Optimization Method with Suitably Graded Truncated Hierarchical B-Spline

    2023-02-26 10:18:02HaoranZhuXinhaoGaoAodiYangShutingWangXiandaXieandTifanXiong

    Haoran Zhu,Xinhao Gao,Aodi Yang,Shuting Wang,Xianda Xie and Tifan Xiong

    School of Mechanical Science and Engineering,Huazhong University of Science and Technology,Wuhan,430074,China

    ABSTRACT This work puts forward an explicit isogeometric topology optimization(ITO)method using moving morphable components(MMC),which takes the suitably graded truncated hierarchical B-Spline based isogeometric analysis as the solver of physical unknown (SGTHB-ITO-MMC).By applying properly basis graded constraints to the hierarchical mesh of truncated hierarchical B-splines(THB),the convergence and robustness of the SGTHB-ITOMMC are simultaneously improved and the tiny holes occurred in optimized structure are eliminated,due to the improved accuracy around the explicit structural boundaries.Moreover,an efficient computational method is developed for the topological description functions(TDF)of MMC under the admissible hierarchical mesh,which consists of reducing the dimensionality strategy for design space and the locally computing strategy for hierarchical mesh.We apply the above SGTHB-ITO-MMC with improved efficiency to a series of 2D and 3D compliance design problems.The numerical results show that the proposed SGTHB-ITO-MMC method outperforms the traditional THB-ITO-MMC method in terms of convergence rate and efficiency.Therefore,the proposed SGTHB-ITO-MMC is an effective way of solving topology optimization(TO)problems.

    KEYWORDS Isogeometric topology optimization; moving morphable components; truncated hierarchical B-spline; suitably graded hierarchical mesh

    1 Introduction

    TO is an engineering optimization method that seeks the optimal material distribution in a prescribed design domain with specified conditions.In the last three decades,a series of TO methods have been proposed and evolved,including solid isotropic material with penalization [1,2],evolutionary structural optimization[3,4],and level set method(LSM)[5–7].Most of these traditional TO methods use the finite element method (FEM) as the solver for the unknown physical field over the design domain,which suffers from numerical instability and geometric discretization errors,due to the low discontinuity between adjacent elements and the disconnection between computer aided design(CAD)and computer aided engineering(CAE).

    Hughes et al.[8]proposed isogeometric analysis(IGA)to substitute the FEM method,which uses the CAD mathematical primitive to represent the field unknowns of the partial differential equations(PDEs).Since IGA has the advantages of high boundary continuity between adjacent elements and elimination of geometric discretization errors,Qian [9] used the relative density of control points as design variables and embedded the physical design domain into a rectangular parametric design domain parametrized by the tensor product B-spline.Gao et al.[10] used ITO for the continuum structure using the enhanced density distribution function.Wang et al.[11] designed a new highefficiency ITO,which improves the efficiency in three aspects: mesh scale reduction,acceleration strategy for the solver,and design variables reduction.Yu et al.[12]presented a multiscale ITO method,where the configuration and layout of microstructures are optimized synchronously.Zhao et al.[13]developed an ITO method for the design problems with an arbitrarily shaped design domain in the framework of T-spline based IGA.Chen et al.[14,15] applied the isogeometric boundary element based TO method to the design problems of acoustic structures for enhancing the sound-absorption performance.Due to the element-wise discreteness inherited from the variable density method,the ITO methods mentioned above easily get into trouble with the curse of dimensionality for generating accurate structural boundaries.

    To resolve the issues that occurred in the variable density-based ITO methods,one alternative is the explicit TO method by Guo et al.[16] using MMC,which originated from the fact that the TO results can be viewed as the optimal distribution of a limited number of components.In the explicit TO method,the geometric characteristic parameters of a limited number of discrete components,e.g.,the coordinates of the centroid and the inclination angle of a component,are treated as the design variables,by which the explicit geometric information of structural boundaries can be generated directly from the converged design variables.Zhang et al.[17]devised the finite-circle method to avoid the overlapping of components for the explicit TO method.Hou et al.[18] adopted the NURBSbased IGA method to obtain the structural response analysis and sensitivity analysis for the MMCbased explicit TO method.Then,an explicit ITO method is proposed with the geometric parameters of moving morphable voids (MMV) [19] treated as the design variables,which utilizes two mesh resolution schemes at different discretization levels.Gai et al.[20]proposed an explicit ITO method that uses the closed B-splines curves to describe the MMV.Zhang et al.[21] integrated the MMV method and the NURBS-based IGA to optimize the 3D shell structure.Since the sensitivity analysis of the explicit ITO depends on the structural boundaries,the accuracy of sensitivities around the structural boundaries determines the optimization accuracy of the explicit ITO method.However,due to the tensor product structure of B-splines,the analytical mesh of explicit ITO cannot perform local refinement of the boundary through the adaptive refinement technique,which leads to a severe contradiction between the analytical efficiency and optimization accuracy.

    To resolve the contradiction between the efficiency and accuracy of the explicit ITO method,the problem of enforcement of global refinement for IGA mesh must be solved first.Currently,the main solutions for this issue are listed as: hierarchical B-splines [22],T-splines [14],LR-splines [23] and hierarchical box-splines[24].Among these methods,the hierarchical B-spline is the most widely used one because of its provable linear independence and subdivision property [25].No?l et al.[26] used truncated hierarchical B-splines to describe the parametrization of the level set function for LSM,which improved the accuracy of the structural boundaries with limited increased computational effort for the updating of design variables.Xie et al.[27]proposed an adaptive explicit ITO method based on hierarchical B-splines,which refines the isogeometric mesh associated with the regions near the optimized structural boundaries locally and achieves a significant balance between the computational efficiency and optimization accuracy.It has been stated that the THB is superior to the hierarchical Bsplines in adaptive IGA because it improves the conditionality of the stiffness matrix and reduces the bandwidth of the stiffness matrix [28].Subsequently,Xie et al.[29] put forward the fully adaptive explicit ITO method with more advanced computing efficiency in terms of truncated hierarchical B-splines,where the isogeometric mesh is locally refined and coarsened simultaneously and the corresponding basis function space recovers the property of partition of unity.However,in the current adaptive explicit ITO method,tiny holes are easily observed in the optimized structures,which should be eliminated for the benefits of both the numerical robustness of the optimization process and the manufacturing of optimized structures.The underlying reason for the aforementioned weaknesses of the adaptive explicit ITO method is that the analytical accuracy of the THB-based IGA method is deteriorated in the regions of structural boundaries,resulting from the non-constrained hierarchical differences between the non-vanishing THB basis functions.

    As pointed out in[30],the admissible mesh can effectively improve the accuracy of the THB-based IGA method.To overcome the shortcoming of the current adaptive explicit ITO method,the present work proposes an advanced adaptive explicit ITO method using MMC in terms of suitably graded truncated hierarchical B-spline (SGTHB-ITO-MMC).Based on the marking strategy proposed in[31],this work develops a constrained marking strategy for imposing the suitably graded constraint on the hierarchical mesh,which enlarges the set of active elements marked to be refined and shrinks the set of deactivated elements marked to be coarsened.With the imposed suitably graded constraint,the accuracy and robustness of the adaptive explicit ITO method are improved.Besides,to further improve the efficiency in computing the TDF values for the hierarchical mesh,we designed an improved TDF calculation strategy for the deletion of the micro components and the local TDF computing strategy for the Gaussian quadrature points.

    The rest of this paper is organized as follows:Section 2“Suitably graded THB based ITO using MMC” introduces the basic theory of the SGTHB-ITO-MMC method,including the definitions corresponding to THB and admissible hierarchical meshes.The optimization model and marking strategy as well as the improved TDF calculation strategy for SGTHB-ITO-MMC are introduced in detail in Section 3“Optimization model of SGTHB-ITO-MMC”.Section 4“Overview scheme for SGTHB-ITO-MMC”presents an overall procedure for SGTHB-ITO-MMC.In Section 5“Numerical examples”,the effectiveness of SGTHB-ITO-MMC is validated by several 2D and 3D benchmarks.Finally,this work is concluded in Section 6“Conclusions”.

    2 Suitably Graded THB Based ITO Using MMC

    This section provides the theoretical foundation for the proposed SGTHB-ITO-MMC,which mainly includes the following two aspects: the construction of THB and the related concepts of admissible hierarchical meshes.

    2.1 Truncated Hierarchical B-Splines

    2.1.1 Hierarchical B-Spline Spaces

    For a given knot vectorΞ={ξ1,ξ2,···,ξn+p+1},nandprepresent the number of control points and the order of basis function,respectively.Supposing that several function spaces of B-spline are nested:V0?V1?...?Vn?1,these nested spaces are obtained by a sequence of knot vectors{Ξ0,Ξ1,...,ΞN?1}and defined on an initial domain Ω0.For multivariate case,Vl(l=1,2,...,N?1)is constructed by the tensor product structure of the univariate basis function space with the knot vectorΞlbisected fromΞ0l-times.Then,the knot intervals constituting a meshQland an arbitrary knot interval are referred to as an element of levell.Define Ω0?Ω1?...?ΩN?1as a nested sequence of domains,where Ωlis the domain of levelland represents the union of the elements onl?1 level marked to be subdivided.Subsequently,a hierarchical mesh HE is made up of active elements on different levels and is defined as:

    where Γlrepresents the union of active cells of levellwithQ?Ωl∧Q/?Ωl+1.

    Once the hierarchical mesh HE is determined [32],the corresponding hierarchical B-spline function space can be constructed by the following steps:

    ? Initialization:H(HE)0={β∈B0: suppβ/=?};

    ? A recursive formula:H(HE)l+1=

    whereHAl+1={β{∈H(HE)l: suppβ/?}Ωl+1}andHBl+1={β∈Bl+1: suppβ?Ωl+1},with suppβ=x:β(x)/=0 ∧x∈Ω0;

    ?H(HE)=H(HE)N?1.

    2.1.2 Truncated Operation

    The hierarchical B-splines basis functions are inevitably spanning the active elements belonging to different levels,which results in the lacking of the essential property of partition of unity for numerical analysis.To recover this important property,a truncated operation should be applied to the two-scale relationship existing in the basis functions of hierarchical B-splines on two consecutive levels,which is formulated as:

    wheretruncl(s)is a truncated active basis function ofsof levell,(s)is termed as the coefficient ofβbelonging to levell+1.Therefore,the THB function spaceH Tcan be generated by the extended truncated operation for active basis functions on non-consecutive levels,as follows:

    An illustration of hierarchical B-splines basis function space and its corresponding THB basis function space are shown in Fig.1.

    Figure 1:Illustration of hierarchical B-splines space H (HE)and THB space H T (HE)

    2.2 Suitably Graded Admissible Meshes

    To ensure that the hierarchical computational mesh of the proposed SGTHB-ITO-MMC model satisfies the specified suitably graded constraints,this section reviews the key definitions related to the admissible hierarchical mesh for THB[30,33],which are elaborated as follows.

    Definition 2.1.For an admissible mesh HE of classm,the truncated basis functions in THB basis function spaceH T(HE)taking nonzero values over any active element belong to no more thanmsuccessive levels,as shown in Fig.2.

    Figure 2:Admissible hierarchical meshes and their associated THB basis function with three different suitably graded constraints

    Definition 2.2.The multilevel support extension of any an elementQ∈Qlon levelk(0≤k≤l)is defined as:

    Definition 2.3.For a levellelementQ,the refining neighborhoodNr(HE,Q,m)with admissibilitymis defined as:

    Definition 2.4.For a levellelementQ,the coarsening neighborhoodNc(HE,Q,m)with admissibilitymis defined as:

    Fig.3 illustrates the definitions shown in Eqs.(4)and(5)as well as(6),which is vital to implement the suitably graded refinement and coarsening for the THB hierarchical mesh.The green elements in Fig.3 represent adjacent elements of level 0 for the active element Q (marked as yellow),which is an essential parameter for the admissible refinement and coarsening algorithms to check whether the active element Q violates the suitably graded constraint or not.The combination of Eqs.(4)and(5)can determine the grid cells adjacent to the marked refining active cell that should be refined locally to satisfy the specified admissible elements,and the combination of Eqs.(4)and(6)can determine the coarsening cells during the admissible coarsening of the hierarchical mesh.

    Figure 3: Illustration of the multilevel support extended domain S,refined neighboring domain Nr,and coarsened neighboring domain Nc for an active element Q

    3 Optimization Model of SGTHB-ITO-MMC

    The mathematical optimization model of the SGTHB-ITO-MMC method is firstly proposed in this section.Then,the constrained marking strategy is devised by integrating the fully adaptive marking strategy with the suitable constraint imposed on the hierarchical mesh,which controls the local refinement and coarsening of the hierarchical mesh for SGTHB-ITO-MMC.Finally,the improved TDF calculation strategy is described to improve the efficiency of SGTHB-ITO-MMC.

    3.1 Mathematical Optimization Model

    According to the optimization model presented in[27],the optimization model of SGTHB-ITOMMC is formulated as:

    In the above formula,drepresents explicit geometric design variables of MMC,which are similar to the shape explicitly expressed by seven parameters,u(x) is the displacement vector of the active control points of the admissible hierarchical mesh,Cdenotes the objective function of the optimization model representing the structural compliance,N0?1 is the maximum level of the hierarchical mesh,denotes the number of active elements of levell,uleis the local displacement vector of thee-th active element belonging to levell,andKrepresent the stiffness matrix of thee-th active element of levell,and the global stiffness matrix,respectively,which are related by the transformation matrixwithcollecting the components ofnon-vanishing on thee-th active element of levell,Frepresents the vector of external force,ρelis the relative material density of thee-th element of levell,andrepresent the volume constraint limit and the elemental volume of thee-th active element of levell,andUddenotes the admissible space for the explicit geometric design variables.

    It should be noted that the ersatz material model and the sensitivity analysis are rather straightforward for SGTHB-ITO-MMC,which can be referred to[29]and are not provided in this work for brevity.In this work,the concept of suitably graded constraint is implemented by altering the marking strategy of hierarchical mesh,which can be presented in the next section.

    3.2 Constrained Marking Strategy

    To implement the adaptivity of hierarchical mesh under the specified suitably graded constraint,a constrained marking strategy is developed.The proposed marking strategy in this work consists of three aspects:triggering strategy,local refining,and local coarsening strategies under suitably graded constraints,where the triggering strategy is calculated in the same way as the one proposed in[29].The local refining strategy and the local coarsening strategy are divided into two stages:(1)initial marking based on TDF values at Gaussian quadrature points;(2)updating the initial marked set under suitably graded constraints.

    3.2.1 Local Refined Marking Strategy

    Once the marking strategy is triggered,the initial local refining of the marking strategy formulated as Eq.(8)finds the elements to be refined by resorting to the topological description function(TDF)values of the active element Gaussian points,as indicated by the green elements shown in Fig.4.

    whererepresents the flag determining whether thee-th element of thelev-th level is to be refined or not,is the TDF value of thes-th Gaussian quadrature point of thee-th active element of thel-th level,is the tolerance limit describing the neighborhood of the structural boundaries.

    If the initial refining element set marked to be refined is obtained,the updating procedures for the initial marked set are described in Algorithms 1 and 2 to fulfill the specified suitably graded constraint during the local refinement of the hierarchical mesh.Algorithm 1 takes the original hierarchical mesh HE,THB basis function spaceT,and the set of initial refining for marking elementsMrefas well as the specified suitably graded constraintmas the input parameters,which plays a role in updating the marked elements in a manner of level-by-level and the entry point for the Algorithm 2.Then,Algorithm 2 refines the elements that do not satisfy the hierarchical constraints by checking the refinement neighborhood of each initial refining element belonging toMref,which obtains from Eq.(8).These elements are added to theMref,thus ensuring that the hierarchical mesh still satisfies the suitably graded requirements of the hierarchical mesh.The newly added elements marked to be refined are shown in the green striped element in Fig.4.

    Algorithm 1 updating_refining_set Input:HE,T,Mref,m Output:Mnew_ref 1:for lev=1:numel(Mref)do 2:Mnew_ref←mark_recursive(HE,T,Mref,lev,m)3:end for

    Algorithm 2 mark_recursive Input:HE,T,Mref,lev,m Output:M′ref 1:neighbors ←hspace_get_H_neighborhood(T,HE,lev,Mlev ref,m)2:if neighb ors!=?t m+1 hen 3:k=lev–4:Mk ref ←Mref ∪neighbors 5:M′ref ←mark_recursive(HE,T,Mkref,lev,m)6:end if

    Figure 4: Illustration of the differences in the variations of the hierarchical mesh between different suitably graded constraints with m=2 and m=∞:(a)Initial hierarchical mesh;(b1)Mesh to refined;(c1)Mesh to coarsened;(d1)Final hierarchical mesh;(b2)Mesh to refined;(c2)Mesh to coarsened;(d2)Final hierarchical mesh

    3.2.2 Local Coarsened Marking Strategy

    When the local refinement of the hierarchical mesh is accomplished,the elements away from the structural boundaries should be marked to be coarsened for SGTHB-ITO-MMC.Similar to the refining strategy,the local coarsening strategy is also divided into two stages:(1)obtaining the initial marking set to be coarsened; (2) updating the coarsened element set.According to the TDF values of the Gaussian quadrature points of the children for each deactivated element,the status of that deactivated element to be refined or not is determined as follows:

    wheredenotes the flag that determines theb-th deactivated element of thelev-th level in the set of admissible deactivated elementsto be coarsened or not withdefined by Eq.(10),denotes the TDF value of theu-th Gaussian quadrature point of thev-th child for theb-th deactivated element oflev-th level,is the limit value used to define the regions far away from the structural boundaries.

    whereis the union of deactivated elements oflev-th level,andare the union of active elements of levellevandlev+1,children(e)denotes all children of a deactivated element which belongs to the subtraction betweenand.

    Then,the marked coarsened set is updated by the procedures presented in Algorithm 3 to guarantee the admissible requirement resulting from the suitably graded constraint.In Algorithm 3,the deactivated element should be reactivated and its children should remain unchanged in the hierarchical mesh,if its coarsening neighborhood with admissibilitymis null;otherwise,that deactivated element should be removed from the reactivated element setMcoa,which is indicated by the blue stripes in Fig.4.

    Algorithm 3 hmsh_coarsen_admissible Input:HE,T,Mcoa,m Output:HE,new_deactivated_elements,M′coa 1:for Q ∈Rc do 2:if Nc(HE,Q,m)==?then 3:Qc←get_children(Q)4:update HE by activating Q and removing its children Qc 5:M′coa ←Mcoa ∪Q 6:end if 7:end for

    Fig.4 shows the variations of the hierarchical mesh during the complete mesh adaptivity under the suitably graded constraints withm=2 andm=∞,from which can be observed that the constrained marking strategy proposed in this work generates smaller and controllable differences in hierarchical level between adjacent active elements.

    3.3 Improved TDF Calculation Strategy for SGTHB-ITO-MMC

    To improve the computational efficiency of SGTHB-ITO-MMC method,this work proposes the following strategy for the benefit of accelerating the calculation of TDF values for the hierarchical mesh.On the one hand,the computing complexity of TDF is reduced by removing the unnecessary TDF computing associated with micro components,with the micro components removed from the design space.On the other hand,the TDF generation is calculated locally inconsistent with the hierarchy feature of the hierarchical mesh.

    3.3.1 Dimensionality Reduction Strategy

    To avoid the numerical singularity,the lower limit of the geometric design variables of MMC is set to an extremely small positive value.It can lead to the numerical phenomenon that the islanding micro component exists in the optimized structure,as shown in Fig.5.

    Figure 5:Illustration of the occurrence of islanding micro component in the optimized result

    These islanding micro components result in the continuously refined mesh around them and that leads to the increasing deteriorative mesh quality.It also results from an excessive number of degrees of freedom,which decreases the computational efficiency of IGA.Therefore,it is necessary to reduce the dimensionality of the geometric design space by removing the design variables associated with micro components in the optimization process.If a component satisfies the requirements of Eq.(11)in the specified consecutive iterations,the associated micro-component is removed from the design space of SGTHB-ITO-MMC.

    whereLiis the half length of thei-th component,t1iandt2ias well ast3iis the half thickness of thei-th component at locations of two ends and middle,εis the user-defined size threshold for determining the micro components.

    3.3.2 Local Computing Strategy

    Apart from the micro components,the redundant TDF calculation also resulted from the way of calculating the TDF values for determining the ersatz material model.In this model,the Gaussian points of the hierarchical mesh inside the active elements are properly graded,and the Gaussian quadrature points do not span different levels of active elements,which enables the local updating of TDF for the hierarchical mesh.Besides,the active elements away from the structural boundaries do not require calculating the TDF values of the Gaussian quadrature points,since SGTHB-ITO-MMC is an essentially boundary-driven TO method.Therefore,the values of TDF of a structural component can be calculated only by Gaussian quadrature points near the boundary of the component,as shown in Fig.6.

    Figure 6: Active Gaussian points (red) and inactive Gaussian quadrature points (black) when ?n setting 0

    In the proposed local TDF computing strategy,the set of active Gaussian quadrature points is defined in Eq.(12)for the hierarchical mesh.

    wherePactiveis then-th level mesh,?ijis the TDF value of thei-th component at thej-th Gaussian quadrature point,?nis the TDF value threshold for then-th level,andPis the point set of Gaussian quadrature points,which are dynamically altered with the adaptivity of the hierarchical mesh.

    4 Overview Scheme for SGTHB-ITO-MMC

    This section mainly describes the execution process of the SGTHB-ITO-MMC method.The execution flow of the SGTHB-ITO-MMC method is shown in Fig.7,which mainly consists of six modules: input module,mesh adaptivity module,ersatz material model module,IGA module,sensitivity analysis module,and output module.

    Figure 7:Illustration of the flowchart of the SGTHB-ITO-MMC algorithm

    Among them,the input module takes the mathematical model parameters in Section 3.1“Mathematical optimization model” as the input parameters and initializes these parameters.The mesh adaptivity module performs the local refinement and coarsening of hierarchical mesh based on the constrained marking strategy as described in Section 3.2“Constrained marking strategy”.In the ersatz material model module,the improved TDF calculation strategy presented in Section 3.3 “Improved TDF calculation strategy for SGTHB-ITO-MMC”is used to improve the efficiency of determining the material physical properties for all active elements.IGA module is represented as the IGA,which is used to calculate the structural responses under the hierarchical mesh.For the sensitivity analysis module,SGTHB-ITO-MMC uses the method of moving asymptotes (MMA) as the optimizer for updating the geometric design variables to achieve the minimization of the structural compliance.Finally,the output module is used to output the optimized explicit geometric parameters of all MMC.

    5 Numerical Examples

    Four numerical examples are used to demonstrate the effectiveness of the SGTHB-ITO-MMC framework with an improved TDF calculation strategy,which is running on MATLAB R2021a with Windows 10 as the software operating system,Intel (R) Core (TM) i5-10210U CPU @ 1.60 GHz–2.11 GHz and 16 GB of RAM as the hardware system.These examples adopted MMA as the optimizer.Young’s modulusEand Poisson’s ratioυare set to 1 and 0.3 as the constant physical values for solid material,respectively.Besides,SGTHB-ITO-MMC withm=∞is identical to the THITO-MMC method presented in[29].

    Through the Messerschmidt-Bolkow-Blohm beam (MBB) problem,Section 5.1 “MBB beam problem” verifies the validity of the suitably graded constraint in improving the optimized results and convergence,and the effectiveness of the improved TDF calculation strategy in enhancing the efficiency for SGTHB-ITO-MMC.In Section 5.2“Short beam problem”,the effectiveness of SGTHBITO-MMC and the improved TDF calculation strategy are validated at two different scales of geometric design space.Finally,the SGTHB-ITO-MMC method is extended to the 3D problem in Section 5.3“3D cantilever beam problem”.

    5.1 MBB Beam Problem

    This section chooses the MBB beam as the first example to verify the effectiveness of SGTHBITO-MMC.Fig.8 depicts the problem configurations and initial design of MMC,where the aspect ratio of the rectangular design domain is 3:1 and the design domain is initially divided into 45×15 uniform IGA mesh.Moreover,the upper limit of the volume fraction is 0.4 for solid material.Maximum iterations equaling 300 and the convergence criterion are used to control the termination of the proposed optimization algorithm.

    Figure 8:Problem model and the initial design for MBB beam problem

    Fig.9 presents the optimized results between different SGTHB-ITO-MMC withm=2,m=3,andm=∞,where the tiny holes that occurred in the converged design of SGTHB-ITO-MMC withm=∞are eliminated in these results by SGTHB-ITO-MMC withm=2,m=3 and the iterative steps is increasing reduced as the stricter suitably graded constraint imposed on the hierarchical mesh.

    Figure 9:Optimized results by SGTHB-ITO-MMC with(a)m=2,(b)m=3,and(c)m=∞

    The convergence histories and the variations of the number of DOFs are illustrated in Fig.10,whereIterationsandObjare referred to as the number of iteration steps and the value of the objective function,respectively.According to the curves presented in Fig.10a,the objective functions are simultaneously reduced as the increase of iterative steps for SGTHB-ITO-MMC,and the SGTHBITO-MMC withm=3 andm=2 reduce the iteration steps by 7.33%and 36.00%than the one withm=∞,respectively.Moreover,the variations in the number of DOFs are depicted in Fig.10b,which indicates that stricter suitably constraint leads to a significant increase in DOFs.However,due to the significant reduction in iterations at the stage with the hierarchies equaling four levels for SGHTBITO-MMC withm=2,the increase in computational effort is very limited for the increased DOFs.Therefore,it is concluded that imposing suitably graded constraints is very necessary to implement the adaptive explicit ITO method in terms of both the convergence rate and numerical robustness.

    Figure 10:Convergence curve of the objective function and the variations of the number of DOFs for SGTHB-ITO-MMC with m=2,m=3,and m=∞:(a)Convergence curve of the objective function;(b)Number of DOFs

    For validating the effectiveness of the improved TDF calculation strategy described in Section 3.3,this work applies it to SGTHB-ITO-MMC withm=2.The optimized structures are depicted in Fig.11 for the SGTHB-ITO-MMC with or without using the proposed improved TDF calculation strategy.Compared with SGTHB-ITO-MMC without using the improved TDF calculation strategy,the number of iteration steps and DOFs are reduced by 9.38%and 12.36%by applying the improved TDF calculation strategy to SGTHB-ITO-MMC,without altering the objective function and the final converged structure.It can be observed that with the use of the improved TDF calculation strategy,the dense mesh disappears away from the structural boundaries as shown in the red box in Fig.11,which contributes to reducing the number of DOFs and improving the optimization efficiency.Fig.12 illustrates the variations of the objective function and the TDF calculation time during the optimization process of SGTHB-ITO-MMC using an improved TDF calculation strategy or not.According to the results presented in Fig.12a,the number of iterative steps for SGTHB-ITO-MMC is reduced by 9.38%,and the total TDF calculation time is decreased by 76.29%when the improved TDF calculation strategy is used.In Fig.12b,the calculation time has abruptly changed at the 80-th,110-th,124-th,and 146-th,which is caused by the fact that the TDF values on Gaussian quadrature of all active cells should be calculated to determine the marked cells when the hierarchical mesh is required to be altered,and the improved TDF calculation strategy is triggered to reduce the calculation time once the adaptivity of the hierarchical mesh is finished.Therefore,the proposed improved TDF calculation strategy has a positive impact on SGTHB-ITO-MMC in both convergence rate and optimization efficiency.

    Figure 11: Comparisons of the optimized results without using and using the improved TDF calculation strategy for SGTHB-ITO-MMC:(a)Without using the improved TDF strategy;(b)Using the improved TDF strategy

    Figure 12:Optimization process of between without using and using the improved TDF strategy:(a)Convergence histories of the objective function;(b)Variations of TDF calculation time

    5.2 Short Beam Problem

    This section discusses the effectiveness of the SGTHB-ITO-MMC for the short beam and validates the improved TDF calculation strategy under two different dimensionalities of the geometric design space.

    5.2.1 Regular Geometric Design Space Case

    The problem setting and initial design of MMC are described in Fig.13.In this problem,the aspect ratio of the design domain is 2:1,which are initially divided into 28×14 uniform IGA mesh.The upper limit of the volume fraction is 0.4 for solid material.The maximum number of steps and convergence criteria are identical to those in Section 5.1“MBB beam problem”.

    Figure 13:Problem model and initial design for short beam problem

    Fig.14 illustrates the optimized results generated by SGTHB-ITO-MMC withm=2,m=3 andm=∞.It is not difficult to find that the convergence rate of SGTHB-ITO-MMC withm=3 andm=2 are respectively improved by 13.99% and 27.16% than SGTHB-ITO-MMC withm=∞.At the same time,the optimized results exist a tiny hole for SGTHB-ITO-MMC withm=∞andm=3,which vanishes on the optimal design by SGTHB-ITO-MMC withm=2.Fig.15 depicts the convergence curves of the objective function for the aforementioned three cases,and it shows the numerical stability of the proposed SGTHB-ITO-MMC method.Finally,the effectiveness of SGTHB-ITO-MMC is verified for the 2D cantilever in a regular geometric design space.

    Figure 14:Optimized results by SGTHB-ITO-MMC with(a)m=2,(b)m=3 and(c)m=∞

    Fig.16 illustrates the effect of the improved TDF calculation strategy on the converged results.Using the improved TDF calculation strategy,the number of iteration steps is reduced by 5.65%.The variation curves of the objective function and TDF calculation time are presented in Fig.17,which show the advantages of the proposed improved TDF calculation strategy in enhancing the convergence rate and the computational efficiency for SGTHB-ITO-MMC.It can be concluded that the improved TDF calculation strategy is valid for the short beam optimized by the SGTHB-ITO-MMC method.

    Figure 15: Convergence curve of objective function for short beam for SGTHB-ITO-MMC with m=2,m=3 and m=∞

    Figure 16:Optimization results of between without using and using the improved TDF strategy:(a)Without using the improved TDF strategy;(b)Using the improved TDF strategy

    Figure 17:Optimization process of between without using and using the improved TDF strategy:(a)Convergence histories of the objective function;(b)Variations of TDF calculation time

    5.2.2 Enlarged Geometric Design Space Case

    To verify the effectiveness of the proposed SGTHB-ITO-MMC in an enlarged geometric design space,the initial designs are shown in Fig.18 for the MMC,with the design domain divided into 4×2 equal sub-regions.During the optimization process,the central coordinates of the MMC must stay in their original sub-regions.The design domain is discretized into 56×28 uniform IGA mesh initially,and the upper limit of volume fraction is 0.4 for solid material.As for the convergence criteria,the one used here is identical to the one presented in Subsection 5.2.1“Regular geometric design space case”,except that the maximum iteration steps are set to 1000.

    Figure 18:Initial design for the short beam with enlarged geometric design space

    Based on the optimized results illustrated in Fig.19 for SGTHB-ITO-MMC withm=2 andm=∞,the boundaries smoothness and the connections of components are significantly improved by applying the suitably graded constraint withm=2 than these obtained by SGTHB-ITO-MMC without considering suitably graded constraint for hierarchical mesh.Moreover,the convergence histories are compared betweenm=2 andm=∞,where the number of iterative steps can be reduced from 1000 to 809 by applying the suitably graded constraint.Meanwhile,Fig.20 expresses the convergence histories are compared betweenm=2 andm=∞,where the number of iterative steps can be reduced from 1000 to 809.Finally,it finds that the proposed SGTHB-ITO-MMC withm=2 is superior to the one withm=∞for the enlarged design space,which shows the robustness of the proposed method concerning the dimensionality of the geometric design space.

    5.3 3D Cantilever Beam Problem

    In this section,the SGTHB-ITO-MMC method is extended to the 3D design problem by taking the 3D cantilever as the last numerical example.As shown in Fig.21,the left facet of the design domain is fixed and the right bottom side is subjected to a uniform vertical downward line load.To obtain the displacement field,the design domain is initially discretized into 15×15×2 tri-quadratic IGA elements with the suitably graded constraints taken asm=2 imposed on the hierarchical mesh.The maximum number of steps is 300 and the convergence criterion

    The initial structural components layout for the 3D cantilever,and two intermediate optimized structural layouts,as well as the final convergent result,are shown in Fig.22,which also includes variations of hierarchical mesh and the convergence process.Under the suitably graded constraints,the boundaries of the 3D cantilever become smoother and the hierarchical mesh is properly graded in the neighborhood of structural boundaries,which verifies the effectiveness of the proposed SGTHB-ITOMMC method in 3D structural TO problems.In addition,the objective function tends to converge with the increase of the iteration step,and there is no obvious oscillation in the objective function value during the duration of the TO process,which reflects the robustness of the SGTHB-ITO-MMC method in solving the 3D design problems.These results show that the SGTHB-ITO-MMC method proposed in this work can effectively solve the 3D compliance TO problem.

    Figure 19:Optimized results by SGTHB-ITO-MMC with(a)m=2 and(b)m=∞in large geometric design space

    Figure 20:Convergence curve of the objective function for SGTHB-ITO-MMC with m=2 and m=∞

    Figure 21:Problem model for 3D cantilever beam problem

    Figure 22:3D cantilever beam convergence curve in the objective function and optimization process for MMC as well as the variations of hierarchical mesh

    6 Conclusions

    This work proposes an explicit ITO method in the framework of suitably graded truncated hierarchical B-splines,which is established based on taking the suitably graded constraint into the consideration of the marking strategy.Furthermore,an improved TDF calculation strategy is put forward by reducing the dimensionality of the geometric design space and generating the TDF values locally for the active elements of the hierarchical mesh.By incorporating the suitably graded constraint into the explicit adaptive ITO method,the drawbacks existing in the optimal MMC designs,such as component discontinuities,tiny holes inside the structure,and zigzag boundaries,can be either eliminated or mitigated to a large extent.Besides,with the stricter suitably graded imposed,the convergence rate of SGTHB-ITO-MMC is increasingly accelerated.With the aid of the proposed improved TDF calculation strategy,the computational efficiency in generating the TDF values and convergence rate are simultaneously improved for SGTHB-ITO-MMC without affecting the optimal designs.Moreover,the proposed method can be applied to 3D design problems.

    In the current explicit adaptive ITO method,the versatility is limited by the basic geometric configuration of MMC.To overcome the aforementioned issue,we will extend the current explicit adaptive ITO method to the MMV framework in the future.Furthermore,the stress-constrained problem solved by the adaptive explicit ITO method is also taken as one of our aims.

    Funding Statement: This work has been supported by the National Key R&D Program of China (2020YFB1708300) and the Project funded by the China Postdoctoral Science Foundation(2021M701310).The authors would like to thank the anonymous reviewers,whose suggestions helped to substantially improve this manuscript.

    Conflicts of Interest:The authors declare that they have no conflicts of interest to report regarding the present study.

    国产精品一区二区在线不卡| 一本一本久久a久久精品综合妖精| 成人国产一区最新在线观看| 国产精品综合久久久久久久免费 | 变态另类成人亚洲欧美熟女 | 成熟少妇高潮喷水视频| 成人国语在线视频| 欧美人与性动交α欧美精品济南到| 国产成人精品久久二区二区免费| 大码成人一级视频| 国产亚洲精品第一综合不卡| 国产精品一区二区在线观看99| 精品一品国产午夜福利视频| 91麻豆精品激情在线观看国产 | av中文乱码字幕在线| 久久国产乱子伦精品免费另类| 国产高清激情床上av| 国产91精品成人一区二区三区| 国产精品免费一区二区三区在线 | bbb黄色大片| 久久久久国产一级毛片高清牌| 最新美女视频免费是黄的| xxx96com| 美女福利国产在线| 久久精品亚洲熟妇少妇任你| 少妇 在线观看| 国产一区二区三区视频了| 身体一侧抽搐| 亚洲伊人色综图| 亚洲国产看品久久| 国产蜜桃级精品一区二区三区 | 怎么达到女性高潮| 在线观看免费高清a一片| 高清欧美精品videossex| 午夜激情av网站| 日韩中文字幕欧美一区二区| 亚洲精品在线美女| 91字幕亚洲| 欧美在线黄色| 亚洲成av片中文字幕在线观看| 国产一区二区三区在线臀色熟女 | 最近最新中文字幕大全免费视频| 色精品久久人妻99蜜桃| 香蕉国产在线看| aaaaa片日本免费| 国产97色在线日韩免费| 人妻久久中文字幕网| 欧美久久黑人一区二区| 亚洲国产欧美网| 国产精品久久久久成人av| 国产精品一区二区精品视频观看| 建设人人有责人人尽责人人享有的| 三上悠亚av全集在线观看| 在线视频色国产色| 午夜免费观看网址| 成人免费观看视频高清| 亚洲aⅴ乱码一区二区在线播放 | 女人久久www免费人成看片| 亚洲精品中文字幕一二三四区| 超色免费av| 日本wwww免费看| 欧美另类亚洲清纯唯美| 人人妻人人添人人爽欧美一区卜| 欧美乱色亚洲激情| 精品一区二区三区视频在线观看免费 | 午夜福利在线免费观看网站| 男女免费视频国产| 嫁个100分男人电影在线观看| 国产国语露脸激情在线看| 精品午夜福利视频在线观看一区| 精品第一国产精品| 一本综合久久免费| 亚洲欧美日韩另类电影网站| 精品电影一区二区在线| 欧美在线黄色| 在线天堂中文资源库| 日本a在线网址| 成年人黄色毛片网站| 十分钟在线观看高清视频www| 午夜成年电影在线免费观看| 窝窝影院91人妻| 老熟妇乱子伦视频在线观看| 亚洲欧美精品综合一区二区三区| 一级a爱视频在线免费观看| 美女视频免费永久观看网站| av在线播放免费不卡| 大香蕉久久网| 国产成人精品无人区| 男女免费视频国产| 在线观看免费日韩欧美大片| 国产欧美日韩综合在线一区二区| 国产精品亚洲av一区麻豆| 国产日韩一区二区三区精品不卡| 成人精品一区二区免费| 变态另类成人亚洲欧美熟女 | 十八禁人妻一区二区| 波多野结衣一区麻豆| 日韩一卡2卡3卡4卡2021年| 欧美成人免费av一区二区三区 | 好男人电影高清在线观看| 精品人妻在线不人妻| 五月开心婷婷网| 国产精品 国内视频| 国产国语露脸激情在线看| 国产精品欧美亚洲77777| 99国产精品99久久久久| 好看av亚洲va欧美ⅴa在| 久久精品国产综合久久久| 一二三四在线观看免费中文在| 男人舔女人的私密视频| 少妇被粗大的猛进出69影院| videosex国产| 国产成人影院久久av| 亚洲成a人片在线一区二区| 黑人欧美特级aaaaaa片| 高清欧美精品videossex| 12—13女人毛片做爰片一| 麻豆乱淫一区二区| 嫩草影视91久久| 国产蜜桃级精品一区二区三区 | 天天躁日日躁夜夜躁夜夜| 亚洲五月色婷婷综合| 乱人伦中国视频| 男女床上黄色一级片免费看| 久久人妻熟女aⅴ| 欧美乱妇无乱码| 亚洲成人手机| 国产在视频线精品| 丰满饥渴人妻一区二区三| 一级a爱片免费观看的视频| 天堂动漫精品| 久久久久久久久久久久大奶| 亚洲成国产人片在线观看| 精品少妇久久久久久888优播| 美女国产高潮福利片在线看| 国产精品久久久久久精品古装| 国产又色又爽无遮挡免费看| 五月开心婷婷网| 国产黄色免费在线视频| 欧美黄色淫秽网站| 久久国产亚洲av麻豆专区| 亚洲专区国产一区二区| 欧美av亚洲av综合av国产av| 亚洲中文av在线| 国产黄色免费在线视频| 欧美一级毛片孕妇| 亚洲中文日韩欧美视频| 少妇 在线观看| 亚洲第一av免费看| 久久久久久久午夜电影 | 久久99一区二区三区| 老司机午夜十八禁免费视频| 精品国产乱子伦一区二区三区| 激情在线观看视频在线高清 | 免费在线观看亚洲国产| 亚洲情色 制服丝袜| 激情在线观看视频在线高清 | 国产单亲对白刺激| 一区二区三区国产精品乱码| 动漫黄色视频在线观看| 叶爱在线成人免费视频播放| 国产精华一区二区三区| 正在播放国产对白刺激| 757午夜福利合集在线观看| 日本五十路高清| 亚洲精品美女久久久久99蜜臀| 亚洲五月天丁香| 成人三级做爰电影| 一二三四社区在线视频社区8| 天天操日日干夜夜撸| 久久久国产成人免费| 国内久久婷婷六月综合欲色啪| 亚洲自偷自拍图片 自拍| 999久久久国产精品视频| www.999成人在线观看| 狠狠婷婷综合久久久久久88av| 一区二区三区精品91| 欧美一级毛片孕妇| 91九色精品人成在线观看| 久久久久久人人人人人| 99精品欧美一区二区三区四区| 亚洲av欧美aⅴ国产| 亚洲第一欧美日韩一区二区三区| 色播在线永久视频| 12—13女人毛片做爰片一| 在线免费观看的www视频| 亚洲五月天丁香| 法律面前人人平等表现在哪些方面| 中文字幕精品免费在线观看视频| 韩国精品一区二区三区| 啦啦啦 在线观看视频| 如日韩欧美国产精品一区二区三区| 丰满人妻熟妇乱又伦精品不卡| 成年版毛片免费区| 久久国产精品男人的天堂亚洲| 天天添夜夜摸| 国产麻豆69| 热re99久久国产66热| 18禁观看日本| 亚洲精品国产精品久久久不卡| 欧美成狂野欧美在线观看| 9色porny在线观看| 不卡一级毛片| 人人妻人人爽人人添夜夜欢视频| 大型av网站在线播放| 国产成人av教育| av国产精品久久久久影院| 欧美精品亚洲一区二区| 国产日韩欧美亚洲二区| 高清视频免费观看一区二区| 亚洲色图 男人天堂 中文字幕| 精品国产一区二区久久| 9热在线视频观看99| 精品人妻熟女毛片av久久网站| 大码成人一级视频| 国产在线一区二区三区精| 一级,二级,三级黄色视频| 日日夜夜操网爽| 在线视频色国产色| 美女 人体艺术 gogo| 欧美日韩精品网址| 久久香蕉国产精品| 亚洲avbb在线观看| 中文字幕精品免费在线观看视频| 午夜免费成人在线视频| 国产主播在线观看一区二区| 久久99一区二区三区| 美女 人体艺术 gogo| 国产人伦9x9x在线观看| 欧美精品啪啪一区二区三区| 波多野结衣av一区二区av| 久久精品亚洲av国产电影网| 国产精品99久久99久久久不卡| 啦啦啦 在线观看视频| 国产激情久久老熟女| 男人舔女人的私密视频| 欧美日韩视频精品一区| 亚洲欧美激情在线| 亚洲国产欧美一区二区综合| 欧美国产精品va在线观看不卡| www.999成人在线观看| 大片电影免费在线观看免费| 国产日韩欧美亚洲二区| bbb黄色大片| 亚洲欧美日韩另类电影网站| cao死你这个sao货| 我的亚洲天堂| 亚洲伊人色综图| 叶爱在线成人免费视频播放| tocl精华| av天堂久久9| 亚洲精品国产区一区二| 国产不卡av网站在线观看| av超薄肉色丝袜交足视频| 777米奇影视久久| 亚洲自偷自拍图片 自拍| 中文亚洲av片在线观看爽 | 国产亚洲欧美98| 欧美日韩国产mv在线观看视频| 妹子高潮喷水视频| 国产精品久久久人人做人人爽| 色精品久久人妻99蜜桃| 在线av久久热| 中文字幕精品免费在线观看视频| 国产xxxxx性猛交| 交换朋友夫妻互换小说| 欧美另类亚洲清纯唯美| 久热爱精品视频在线9| 另类亚洲欧美激情| 午夜精品在线福利| 国产激情欧美一区二区| 久久精品国产a三级三级三级| 欧美激情高清一区二区三区| 亚洲第一欧美日韩一区二区三区| 久久久久国产精品人妻aⅴ院 | 亚洲综合色网址| 免费日韩欧美在线观看| 精品视频人人做人人爽| 最新的欧美精品一区二区| 久久人妻av系列| 国产亚洲一区二区精品| 亚洲av美国av| 在线观看日韩欧美| 亚洲一区二区三区不卡视频| 国产单亲对白刺激| 免费在线观看视频国产中文字幕亚洲| 天天影视国产精品| 国产成人精品无人区| 国产精品一区二区在线不卡| 欧美日韩黄片免| 欧美激情极品国产一区二区三区| 99久久国产精品久久久| 香蕉久久夜色| 成年动漫av网址| 十八禁人妻一区二区| 精品国产乱子伦一区二区三区| 嫩草影视91久久| 91av网站免费观看| 欧美中文综合在线视频| 国产亚洲一区二区精品| x7x7x7水蜜桃| 校园春色视频在线观看| 国产精品98久久久久久宅男小说| 在线永久观看黄色视频| 制服诱惑二区| 少妇猛男粗大的猛烈进出视频| 亚洲少妇的诱惑av| 亚洲一码二码三码区别大吗| 夜夜爽天天搞| 亚洲欧洲精品一区二区精品久久久| 在线观看免费视频日本深夜| 国产精品二区激情视频| 人人妻人人澡人人看| 精品久久久精品久久久| 亚洲国产看品久久| 最新的欧美精品一区二区| 午夜免费鲁丝| 黑人巨大精品欧美一区二区蜜桃| 久久精品人人爽人人爽视色| 成人黄色视频免费在线看| 夜夜躁狠狠躁天天躁| 免费日韩欧美在线观看| 久久精品国产a三级三级三级| av有码第一页| 99riav亚洲国产免费| 精品国内亚洲2022精品成人 | 少妇被粗大的猛进出69影院| 成年人黄色毛片网站| 精品高清国产在线一区| 91成人精品电影| 高清黄色对白视频在线免费看| 亚洲国产中文字幕在线视频| 黑人巨大精品欧美一区二区mp4| 亚洲av美国av| 国产麻豆69| 又黄又爽又免费观看的视频| 欧美黑人精品巨大| 一区在线观看完整版| 亚洲欧美一区二区三区久久| 搡老岳熟女国产| 1024视频免费在线观看| 国产精品一区二区在线不卡| 三上悠亚av全集在线观看| 叶爱在线成人免费视频播放| 黑人猛操日本美女一级片| 美女视频免费永久观看网站| 国产精品av久久久久免费| 国产成人啪精品午夜网站| 悠悠久久av| 99在线人妻在线中文字幕 | 天天躁狠狠躁夜夜躁狠狠躁| 色精品久久人妻99蜜桃| 国产1区2区3区精品| 91成人精品电影| 又大又爽又粗| 午夜久久久在线观看| 亚洲欧美日韩另类电影网站| 亚洲久久久国产精品| 午夜精品久久久久久毛片777| av电影中文网址| 80岁老熟妇乱子伦牲交| 在线天堂中文资源库| 国产精品久久久av美女十八| 久久人妻熟女aⅴ| 老司机午夜十八禁免费视频| 精品久久久精品久久久| 在线观看免费视频网站a站| 免费女性裸体啪啪无遮挡网站| 日本vs欧美在线观看视频| 久久精品国产综合久久久| 村上凉子中文字幕在线| 亚洲av第一区精品v没综合| 国产激情久久老熟女| 身体一侧抽搐| 90打野战视频偷拍视频| 精品久久久久久久久久免费视频 | 亚洲 国产 在线| 久9热在线精品视频| 悠悠久久av| 久久精品亚洲av国产电影网| 少妇粗大呻吟视频| 视频区图区小说| 亚洲专区中文字幕在线| 极品人妻少妇av视频| 丝瓜视频免费看黄片| 黄色片一级片一级黄色片| 欧美日韩亚洲高清精品| 很黄的视频免费| 女人被躁到高潮嗷嗷叫费观| 免费在线观看黄色视频的| 纯流量卡能插随身wifi吗| 99香蕉大伊视频| 飞空精品影院首页| 99久久国产精品久久久| 国产精品影院久久| 精品久久蜜臀av无| av视频免费观看在线观看| 欧美国产精品一级二级三级| a级片在线免费高清观看视频| 男女高潮啪啪啪动态图| 国产不卡av网站在线观看| 丰满迷人的少妇在线观看| 免费在线观看完整版高清| 50天的宝宝边吃奶边哭怎么回事| 日本精品一区二区三区蜜桃| 亚洲情色 制服丝袜| 捣出白浆h1v1| 91精品国产国语对白视频| 成人免费观看视频高清| 热re99久久国产66热| 国产精品av久久久久免费| 国产在视频线精品| 人成视频在线观看免费观看| 变态另类成人亚洲欧美熟女 | 国产精品av久久久久免费| 精品卡一卡二卡四卡免费| 黑丝袜美女国产一区| 看黄色毛片网站| 亚洲九九香蕉| 国产精品国产高清国产av | 成人av一区二区三区在线看| 日韩有码中文字幕| 亚洲国产欧美日韩在线播放| 亚洲精华国产精华精| 国产亚洲精品一区二区www | 黑人猛操日本美女一级片| 国产一区在线观看成人免费| 伦理电影免费视频| 日韩欧美三级三区| 亚洲欧美日韩高清在线视频| 一个人免费在线观看的高清视频| 国产在视频线精品| 啪啪无遮挡十八禁网站| 亚洲三区欧美一区| 三级毛片av免费| 国产男靠女视频免费网站| 久久久久国产一级毛片高清牌| 国产不卡一卡二| 欧美国产精品一级二级三级| 人妻一区二区av| 亚洲精品av麻豆狂野| 女性生殖器流出的白浆| 天堂中文最新版在线下载| 黑人操中国人逼视频| 欧美激情高清一区二区三区| www日本在线高清视频| 免费人成视频x8x8入口观看| 一级黄色大片毛片| 在线观看免费午夜福利视频| 男女午夜视频在线观看| 多毛熟女@视频| a级片在线免费高清观看视频| 免费在线观看日本一区| 国产精品自产拍在线观看55亚洲 | 国产精品久久视频播放| 国产99白浆流出| 免费观看a级毛片全部| 乱人伦中国视频| 国产成人欧美在线观看 | 一级毛片女人18水好多| 久久国产精品大桥未久av| 亚洲自偷自拍图片 自拍| 色精品久久人妻99蜜桃| 12—13女人毛片做爰片一| 国产亚洲欧美98| 好男人电影高清在线观看| 女人精品久久久久毛片| avwww免费| 夫妻午夜视频| 捣出白浆h1v1| 黑人操中国人逼视频| 欧美在线一区亚洲| 精品国产超薄肉色丝袜足j| 午夜老司机福利片| 无遮挡黄片免费观看| 成人18禁在线播放| 亚洲精品美女久久av网站| 亚洲av片天天在线观看| 亚洲七黄色美女视频| 亚洲欧美激情在线| 国产成人av教育| 国产真人三级小视频在线观看| 国产片内射在线| 99国产精品免费福利视频| 日日摸夜夜添夜夜添小说| 国产有黄有色有爽视频| 久热爱精品视频在线9| 中文字幕人妻熟女乱码| 午夜影院日韩av| 久久久久精品国产欧美久久久| 99国产精品免费福利视频| 国产亚洲精品一区二区www | 午夜视频精品福利| 久久精品国产99精品国产亚洲性色 | 很黄的视频免费| 一级a爱片免费观看的视频| 久久久久精品国产欧美久久久| 国产欧美亚洲国产| 精品久久久精品久久久| 国产精品 国内视频| 少妇粗大呻吟视频| 婷婷丁香在线五月| 亚洲精品美女久久久久99蜜臀| 在线播放无遮挡| 69av精品久久久久久| 精品一区二区三区视频在线观看免费| 老司机午夜十八禁免费视频| 久久久成人免费电影| 窝窝影院91人妻| 亚洲av一区综合| 日本成人三级电影网站| 亚洲人与动物交配视频| 午夜a级毛片| 国内精品久久久久精免费| 乱人视频在线观看| 网址你懂的国产日韩在线| 一本一本综合久久| 国产午夜精品久久久久久一区二区三区 | 国产成+人综合+亚洲专区| 母亲3免费完整高清在线观看| 精品一区二区三区人妻视频| 免费在线观看亚洲国产| 丰满的人妻完整版| 小蜜桃在线观看免费完整版高清| 欧美成人a在线观看| 日韩欧美 国产精品| 老司机福利观看| 久久中文看片网| 久久久久免费精品人妻一区二区| 在线观看午夜福利视频| 国产精品久久久久久亚洲av鲁大| 国产久久久一区二区三区| 久久精品综合一区二区三区| 日韩精品中文字幕看吧| 麻豆国产av国片精品| 国产一区在线观看成人免费| 怎么达到女性高潮| 欧美激情在线99| 不卡一级毛片| 精品日产1卡2卡| 国产精品一区二区三区四区久久| 在线十欧美十亚洲十日本专区| 国产亚洲av嫩草精品影院| 午夜影院日韩av| 国产高清三级在线| 国产精品女同一区二区软件 | 欧美乱色亚洲激情| 亚洲成av人片免费观看| 亚洲乱码一区二区免费版| 男女那种视频在线观看| 香蕉丝袜av| 全区人妻精品视频| 夜夜夜夜夜久久久久| 午夜免费观看网址| 欧美乱色亚洲激情| 国产亚洲精品久久久com| 97人妻精品一区二区三区麻豆| 国产精品一区二区免费欧美| 俄罗斯特黄特色一大片| 国产一区二区三区视频了| www国产在线视频色| 欧美黄色片欧美黄色片| 国产午夜精品论理片| 男人和女人高潮做爰伦理| 日韩成人在线观看一区二区三区| 观看免费一级毛片| 观看美女的网站| 亚洲中文字幕一区二区三区有码在线看| www国产在线视频色| 国产高清视频在线观看网站| 亚洲国产精品合色在线| 欧美日韩一级在线毛片| 一区福利在线观看| 在线播放无遮挡| 俺也久久电影网| 久久人人精品亚洲av| 一本精品99久久精品77| 欧美日韩乱码在线| 一进一出抽搐gif免费好疼| 国产精品久久久久久久久免 | 亚洲av一区综合| 成年女人看的毛片在线观看| 精品日产1卡2卡| 国产成人福利小说| 又紧又爽又黄一区二区| 两人在一起打扑克的视频| 欧美午夜高清在线| 韩国av一区二区三区四区| 免费观看的影片在线观看| 亚洲国产精品sss在线观看| 人妻丰满熟妇av一区二区三区| 亚洲黑人精品在线| 精品人妻偷拍中文字幕| 制服人妻中文乱码| 国产伦人伦偷精品视频| 人妻夜夜爽99麻豆av| 精品99又大又爽又粗少妇毛片 | 免费电影在线观看免费观看| 国产成+人综合+亚洲专区| 别揉我奶头~嗯~啊~动态视频| 色噜噜av男人的天堂激情| 在线播放无遮挡| 黄片大片在线免费观看| 亚洲成人中文字幕在线播放| 99久久精品一区二区三区| 嫩草影院入口| 日韩欧美精品v在线| 精品国内亚洲2022精品成人| 中国美女看黄片| 久久人人精品亚洲av| 国产激情欧美一区二区| 亚洲精华国产精华精| 在线视频色国产色| 99热精品在线国产| 18禁裸乳无遮挡免费网站照片|