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

    A novel method for solving shortest tool length based on compressing 3D check surfaces relative to tool postures

    2021-04-06 10:25:36XiangyuLIJunxueRENXimingLVYukeZHOUCongleLIU
    CHINESE JOURNAL OF AERONAUTICS 2021年2期

    Xiangyu LI, Junxue REN, Ximing LV, Yuke ZHOU, Congle LIU

    Key Laboratory of High Performance Manufacturing for Aero Engine (Northwestern Polytechnical University), Ministry of Industry and Information Technology, Xi’an 710072, China

    KEYWORDS CNC;Compressed check surface;Multi-axis machining;Shortest tool length;Silhouette curve;Tool holding system

    Abstract Solving the shortest tool length quickly under a known tool trajectory in multi-axis machining of complex channel parts is an urgent problem in industrial production.To solve this problem, a novel and efficient method is proposed which is featured by extracting only a few necessary curves from the check surface instead of sampling the entire surface. By rotating and compressing the 3D check surface relative to all tool postures,the boundaries of the area occupied by the 2D compressed surfaces are the essential elements for determining the shortest tool length.A tracking-based numerical algorithm is introduced to efficiently solve the silhouette curves which are formed in compressing.To define the multi-taper shaped tool holding system(THS)which is commonly used in production,a characterization model for THS profile is established.A model for solving the shortest tool length is finally constructed based on the critical interference relationship between the THS profile and all compressed boundary curves.For acceleration,the boundary splines are segmented according to their knot vectors.Then a new concept called the axis-aligned tool length box(AATB)is introduced,which can provide a conservative range of tool length for a spline segment.By scanning the AATBs of all spline segments,the very few effective spline segments that may ultimately determine the shortest tool length are filtered out.This acceleration method makes the solution for the shortest tool length more focused and efficient.The results of experimental examples are also reported to validate the efficiency and accuracy of the proposed algorithm.

    1. Introduction

    With the advanced design concepts, more and more complex channel parts enclosed by free-form surfaces are adopted in modern turbomachinery, such as the blisk in aero-engine and the closed impeller in gas-turbine. Multi-axis CNC machining has been the exclusive means for manufacturing this kind of part owing to its unique capability of continuously changing the tool orientation. However, the interference-free cutter which can reach the deepest part of the channel has to be thin and long,1as illustrated in Fig.1(left).As a result,the cutter is too weak in stiffness and thus extremely susceptible to deflection,2,3vibration4,5and wear,6,7which is particularly worrisome as it often severely reduces the cutting efficiency8–10and undermines the geometric accuracy of the final machined part.11,12To strengthen the cutter stiffness in multi-axis machining,especially for hard-to-cut materials,in general,current research adopts two natural strategies:shortening the tool length of the entire tool path or optimizing the tool shank shape.13–15Compared with the latter, the former strategy is more economical and easier to implement in practice. Hence,the shortest tool length is preferred and widely adopted in practice. Besides, the shortest tool length is also considered as an important factor in scheduling the tool sequence,as well as optimizing tool axes.16

    Quickly solving the shortest tool length with a reference tool trajectory has always been an urgent problem in industrial production. Generally, an ordinary multi-axis machining program contains hundreds or even thousands of tool postures,if we try to calculate the exact shortest tool length of each posture one by one, it will require a large number of interference checking operations. Such a heavy calculation load often makes the task of solving shortest tool length very timeconsuming (even up to hundreds of minutes), which far exceeds the expectation of programmer (preferably in a few minutes). An easier way in practice is, based on the programmer’s experience about the structure of the channel, choosing those tool postures located in deepest part of the channel to participate in calculation, as they are considered most likely to ultimately determine the shortest tool length. Then the initially obtained shortest tool length will be updated until all postures are just interference-free.This approach is less robust because of its heavy reliance on human experience. Even so,the computing efficiency is still unacceptable, especially for the complicated tool holding system.To further reduce the calculation load,one conservative approach is designing a simple safe geometry(e.g.,the safe cylinder in Fig.1(left)),or bounding box) to replace the complex check surfaces in interference detection. However, the conservative tool length is often too long which may severely reduce the tool stiffness,just as shown in Fig.1(left).Conceivably,if a smaller and more complicated tool handle(e.g.,the conical tool handle which is gaining popularity in manufacturing) is adopted and allowed to enter the channel,the tool length may be significantly reduced,as shown in Fig. 1(right). In fact, besides tool handle, the spindle or other related geometries of machine tool may also restrict the shortest tool length in some special cases. Therefore, to ensure the calculated shortest tool length more reliable, the actual shape of the tool handle, spindle, and other related geometries have to be fully considered. Certainly, it will require much more computing resources.

    Fig. 1 A long cutter with cylindrical tool handle (left) and a short cutter with conical tool handle (right).

    Regarding efficiently calculating the shortest tool length itself, some key works, though not too much, focusing on this problem is perhaps worthwhile to review. Inui et al.17–19proposed a method to calculate the shortest tool length in 3-axis and 5-axis machining. The length was derived by computing the distance between the cutter location points and an offset shape of the workpiece by the inverted zero-length cutter(i.e., the Minkowski sum shape). As the operations of offsetting workpiece and measuring minimum distance need a large number of calculations, the method is inefficient. Besides, the method is limited in dealing with the simplest cylindrical tool handle.Also,Kaneko and Horio20implemented a similar algorithm. Lee et al.21used a hybrid model which is composed of the Z-map model and the Triangulation model to inspect the collision caused by tool holder,then they applied an increment method to find the optimal tool length. Different from these studies, Ahmed and Chen22used the technology of Kd-tree to gradually remove most of the irrelevant geometric elements,i.e., the discretized triangular patches. As a result, the efficiency of calculating the shortest tool length was improved.Actually, in essence, solving the shortest tool length is still the problem of interference detection. Detection of the minimum distance23is the most commonly used method.To generalize the method to different complex geometries, both the obstacle surface and the tool handle have to be discretized into simple geometric elements, such as points, line segments, and triangular patches. By judging the minimum distance between the tool handle and the obstacle surface against the safe value,the shortest tool length can be obtained. Because the topological relationship between two sets of discrete geometries no longer needs to be concerned, hence,the simplicity and strong robustness are the significant benefits of this method. On the negative side, the calculation amount is typically huge (e.g.,millions of points),as a result,the calculation speed and accuracy are heavily dependent on discrete methods. This type of operation is most suitable for running on a graphics processing unit (GPU) with many simple processors. For example, Bi et al.24,25and Wang et al.26,27adopted GPU to accelerate the solution of the shortest tool length. However, testing all discretized points on the check surface in the above methods is not the best choice. Conceivably, most of the sampling points are radially away from the tool axial line and are impossible to determine the shortest tool length.

    Fig.2 Revolving space formed by rotating check surface around tool axial line.

    A complete tool holding system (THS) contains tool handle, spindle, and other related geometry bodies. Generally,THS is a revolving body with respect to the tool axial line. It is noticed that the tool length is minimum when THS is just contacting with check surfaces at a critical point p, as shown in Fig.2.Let h be the axial height of p from the cutter tip point pct. The plane crossing p and perpendicular to tool axis T is denoted as Λ. The intersection curve of plane Λ and channel surfaces S is denoted as c. Actually, the critical point p is the closest point on curve c to the tool axial line and the minimum distance is exactly the radius r of THS at p. Actually, r and h are more important than the specific position of p on the curve c for solving the shortest tool length. On the other hand, all points on c except p are far away from the tool axial line.Theoretically, these points are useless for solving the shortest tool length and should be abandoned to reduce the computational burden. Hence, searching the critical point is a breakthrough for efficiently determining the shortest tool length, which is also the innovation of this paper.

    Imagine that if we rotate the check surface S around the tool axial line,then the formed revolving space Ω is the design space of tool and THS.Obviously,the meridian ζ of Ω is made up of all critical points. Therefore, our idea about solving the shortest tool length is: adjusting the tool length to make the THS profile is just contacting with the meridian of the revolving space of check surface, the resultant tool length is the shortest. Some research works also mentioned the revolving space. Kim28introduced a concept of safe space between the tool axis and the check surface under a known tool path.Chaves-Jacob et al.29proposed the concept of Adapted Tool Shape and used it to optimize the cutter shape in flank milling to reduce the machining error. However, their description of solving the revolving space is simplified. To achieve the best conical cutter with maximum stiffness, a similar concept of Cutter Profile Design Space (CPDS) was also proposed in our previous research.13Unfortunately, none of the above three works concern about the problem of solving the shortest tool length.

    In a different note, considering the shortest tool length is the global minimum for all tool postures in a tool path. The bottleneck tool postures that restrict the global shortest tool length are rarely, especially for those huge NC programs.Chiou et al.16also noticed this phenomenon. They shortened the global shortest tool length by adjusting the tool axes of those bottleneck tool postures. Hence, calculating the exact shortest tool lengths of every tool posture one by one is inefficient and meaningless.Because most of them will eventually be replaced by the minimum one. Inspired by this idea, filtering out the few effective meridians of all revolving spaces will accelerate the process of solving the global shortest tool length(see Section 3 for the details).Motivated by the lack of an efficient algorithm for calculating the shortest tool length in a complex channel enclosed by freeform surfaces, to fulfill the need, we present a novel solution with fully considering the actual shape of THS.The main components,as well as the contributions of our solution,are as follows:

    ●Given a tool posture, by rotating and compressing a check surface with respect to the tool axial line, the boundary curves of the 2D compressed surface are the necessary elements to determine the shortest tool length.Then,a numerical tracking algorithm is presented to efficiently trace out the silhouette curves, which are newly formed boundaries in compressing.

    ●To define the multi-taper shaped THS which is most commonly used in industry,a generic model of tool holding system is proposed. Based on the critical interference relationship between the THS profile and all compressed boundary curves, the model for solving the shortest tool length is established.

    ●To accelerate the solution, these compressed boundary splines are divided into spline segments according to their knots.Then,we propose a new concept of axis-aligned tool length box (AATB) for a single spline segment, which can be quickly obtained by taking advantage of the convex hull property of NURBS spline. Based on all AATBs, the very few effective spline segments which might ultimately determine the shortest tool length can be quickly filtered out.

    The remainder of this paper is organized as follows. In Section 2, we present the model for solving the shortest tool length based on the THS profile and all boundary curves of the compressed check surface. Also, an efficient algorithm for tracking the silhouette curves is introduced herein. Next,in Section 3, by filtering out the very few effective boundary spline segments, a novel acceleration method for solving the shortest tool length is proposed. Section 4 reports the detailed results and discussions of examples. Finally, we conclude the paper in Section 5.

    2. Model for solving the shortest tool length

    As already alluded,when the THS is just contacting the check surfaces for all postures, the tool length is minimum. According to this strategy, we will build the model of solving the shortest tool length in this section.Specifically,a generic model of tool holding system and the method for extracting the boundaries of compressed check surfaces are needed and will be described as follows.

    2.1. Characterization model of generic tool holding system

    Although there are many types of THS in industrial production, most of them are multi-taper shaped and their profiles are composed of line segments (Note: those very special profiles are not considered in this article, e.g., spline and conic curve), as shown in Fig. 3 (left). Taking the cutter tip point pctas the origin of a two-dimensional cutter coordinate system,the tool axial line as the h-axis and the radius of THS as the raxis,then the THS profile can be characterized by multiple line segments which make the THS radius is monotonically increasing along the tool axis, as illustrated in Fig. 3(right).Obviously,the height of the THS profile is the function of tool length L and can be abstractly expressed as Γ(L). Let us describe the characterization of Γ(L) as follows.

    Fig. 3 An actual THS (left) and characterization of THS profile in cutter coordinate system (right).

    The distribution of r-coordinate for all endpoints of Γ(L)segments is {R1,R2,...,Rk,...,RK,RK+1}, where 1 ≤k ≤K, Rk

    Since Γ(L)can only move along the h-axis,to facilitate the description of the filtering method in Section 3, the line r=Rk(i=1,2,...,K) is called a transboundary line of THS profile, and the area of Rk

    2.2. Compressing 3D check surface relative to tool posture

    As already explained in the introduction, the meridian of the revolving space formed by rotating the check surface relative to the tool axial line is the key element for solving the shortest tool length.It is noticed that the meridian is formed from a few special curves on the check surface.This property can be guaranteed for the channel parts in aero-engines because their enclosed surfaces are always free-from and typically very smooth. Let us describe the compressing method and a tracking-based algorithm for efficiently solving the meridian.

    2.2.1. Rotating and compressing check surface relative to tool posture

    Referring Fig. 4, the smooth surface S is one of the original check surfaces. Normally, S can be represented as a NURBS surface and can be written as Eq. (2) (see the Eq. (4.13) in30),{ pi,j} are its control points. A known cutter posture (pct,T)is defined by the cutter tip point pctand the unit tool axis vector T.

    In fact,we are most concerned about the boundaries of the area occupied by the compressed surface. Because only these boundaries participate in determining the shortest tool length,while the interior points of the area are not. That is the main reason why the method of sampling the entire check surface is inefficient. A compressed surface usually has two types of boundaries, i.e., the four edges and some silhouettes which are newly formed in the compressing process. The edge curves can be easily extracted based on the surface equation and omitted here.As for the silhouettes,any point on it is called a characteristic point (CP). On the compressed surface C, CP is the closest perpendicular point to the h-axis with current axial height. Then, the silhouette curves (i.e., all CPs) need to be searched out and will be presented in the next section.

    Fig. 4 A check surface is rotated about tool axial line and compressed into a 2D folded surface.

    2.2.2. Tracking silhouette curves

    Generally,a silhouette curve starts from and ends at the edge of the compressed surface,which is denoted as a type I silhouette.However,it could happen(though extremely rare)that a silhouette curve entirely lies in the interior of the compressed surface,and this kind of silhouette is denoted as a type II silhouette.The details of efficiently tracking for silhouettes on parametric surface C(u,v),i.e.,C(u,v)=[Cr(u,v)Ch(u,v)]T,will be presented as follows. To ensure the calculated CPs are uniformly distributed (within the given tolerance) on the silhouette curve,let Ls(e.g., 0.01 mm) be the minimum Euclidean distance between any two neighboring CPs on the silhouette.

    (1) Geometric property of the silhouette curve

    For any CP P⊥(Cr(u,v),Ch(u,v)) on silhouette with tool axial height h⊥, it should satisfy the condition of ‘‘the closest perpendicular point”. To describe it, we construct a Lagrange-form function, as in Eq. (6), where variable λ is the Lagrange multiplier. With this combined form, the CP should satisfy the extreme condition of Eq. (6), as given in Eq. (7).

    Obviously, for any potential CP on surface C, if it lies on the silhouette curve,its parameter(u,v)must satisfy the condition W(u,v) = 0.

    (2) Tracking for type I silhouette

    Since any type I silhouette begins and ends at some boundary edges of C(u,v), we will numerically construct the type I silhouette curve by the following procedure:

    Step 1: Find the starting CPs on the edges of the compressed surface.

    All starting CPs can be found by the following strategy.Divide each edge of surface C(u,v) into enough intervals.The interval size is no more than Ls. If the signs of W(u,v)are different at the two ends of one specific interval, then this interval is considered containing a starting CP(here we made a reasonable assumption:there will be at most one starting CP in an interval owing to the smallness of Ls). The accurate target CP will be found out by the binary searching method. In contrast,the remaining intervals need to be further tested.The test method is:randomly sample Ls/ε points in the test interval(ε is the required displacement accuracy of the linear axes of any typical multi-axis CNC machine tool, e.g., 10-3mm); if the signs of W(u,v)at two neighboring sample points are different,then their midpoint is a starting CP within the tolerance.

    Step 2: Track an entire type I silhouette.

    Randomly pick a starting CP as an initial point for iteration, the tracking procedure towards the interior of C(u,v) is then executed to trace out the entire type I silhouette. In general, the tracking will stop at another starting CP which also lies on one edge of C(u,v). The two starting CPs then are marked as ‘‘searched”. The iterative tracking process from an obtained CP Ck(u,v) to the next CP Ck+1(u+Δu,v+Δv) is performed in two steps:

    Step 2.1: Initial search along the contour of W=0 by step Ls.

    To restrict the search length, the Euclidean distance between the current CP Ckand the next one Ck+1is set to be Ls, thus we have Eq. (9). Expanding it into Taylor series and omitting the second-order term and after, yields Eq.(10):

    Substituting Eq. (12) into Eq. (11), then Eq. (13) is obtained.

    After Step 2.1,in general,due to the omission of the second and all the higher-order terms in Eq. (10), the actual W(u+Δu(k),v+Δv(k)) does not equal to zero in general, as shown in Fig.6.To ensure the correctness and quality of tracking, we need to further search from (u+Δu(k),v+Δv(k)) to a more accurate Ck+1, which is Step 2.2.

    Step 2.2: Accurate search along the gradient of W by an iterative method.

    Taking (u+Δu(k),v+Δv(k))(which is obtained in Step 2.1)as the initial point, we employ the well-known twodimensional damped Newton method to accurately track a neighboring CP that meets the condition W=0, as shown in Fig.6.This Newton method is well known for its quadratic convergence.The iterative formula for the propagation of surface parameter (u, v) is given in Eq. (16). Using this formula,the exact parameters of the next CP can be obtained after a few iterations under a given threshold of accuracy:

    Fig. 6 Tracking type I silhouette on compressed surface.

    Step 3: Construct a NURBS curve through all CPs which are searched in Step 2.The simple chord length method which is described in30is adopted here (the details are omitted here).As a result, this curve is an identified type I silhouette curve.

    Step 4: Repeat Step 1 to Step 3 until all starting CPs are marked as ‘‘searched”. It indicates that all type I silhouette curves have been found.

    (3) Tracking for type II silhouette

    Normally,type II silhouettes rarely appear in smooth channel surfaces, such as the aero-engine parts. We numerically search and construct type II silhouettes by the following simple procedure:

    Fig. 5 Determination of searching step: initial tracking (left) and intermediate tracking (right).

    a. Sample the compressed surface into a dense set of points based on some controlling tolerance.Then surface C and its uv domain are divided into many patches based on these sample points. We reasonably assume that there is at most one CP in each patch.

    b. Identify the patches containing a type II CP. If the W signs of the four corner points of a patch are not the same with each other,which illustrates that the existence of a type I CP or a type II CP in this patch, then the patch is marked. Remove those patches crossed by the type I silhouettes from all the marked patches. As a result,the remaining patches are identified as containing type II CP.

    c. Taking the midpoint of any identified patch as the initial point, we still employ the damped Newton method to accurately search for an interior type II CP that meets the condition W=0. Starting from this initial CP, the tracking process of Step 2 in Section 2.2.2 will be executed and eventually stop once a loop (i.e., a type II silhouette) is formed. Similarly, execute Step 3 in Section 2.2.2 and construct a type II silhouette curve.Next, remove the patches crossed by the obtained silhouette.

    d. Repeat (c) until there is no marked patch containing type II CP,which means all type II silhouettes have been found.

    So far, all the boundary curves of the compressed surface are obtained, including the edges and the newly formed type I and type II silhouettes.

    2.3. Model of shortest tool length

    As explained previously,a combination of a tool posture and a check surface corresponds to a 2D compressed surface.Referring Fig. 7(left), given a channel enclosed by M freeform surfaces (i.e., the check surfaces) and a toolpath containing N non-interference tool postures, by the means in Section 2.2, a total of M×N compressed surfaces are obtained and accompanied by s boundary curves Bi(μ), i = 1,2,...,s,s ≥4×M×N, μ is the curve parameter. In the example shown in Fig. 7, M=1, N=34, s=170, there are 34×4=136 compressed edges and 34 compressed silhouettes(all are type Ⅰ).

    If the THS profile is just contacting with Bi(μ) (see Fig. 7(right)),it means the THS is also just contacting with the check surfaces. As a result, the current tool length is exactly the shortest tool length Lshortest, as given in Eq. (17):

    Eq. (17) is an abstract expression, let us present a clearer form which is easier to implement in below. Sampling Bi(μ)into a dense set of points based on some controlling tolerance,let bi,p(Bi,p,r,Bi,p,h) be the pth sample point on the boundary curve Bi(μ), p = 1,2,...,ni, niis the number of samples. If Rk

    Eq. (18) is the mapping from any boundary point to the tool length restricted by itself. Obviously, it is an important basis for determining the shortest tool length. Substitute Eq.(18) into Eq. (17), Lshortestcan be intuitively expressed as:

    Eq. (19) is exactly the model for solving the shortest tool length with a given tool trajectory in a complex channel,which is based on the THS profile is just contacting with all boundaries of compressed surfaces.

    3. Acceleration for shortest tool length

    Although we have found out all necessary boundaries of compressed surfaces in advance,the operation of sampling boundaries in Eq. (19) still needs a lot of calculation resources.Actually, referring Fig. 7(right), we can draw an interesting observation:only part of the outermost envelope of all boundaries are possible to ultimately determine the shortest tool length while most of the internal boundaries are not. Trying to obtain the complete envelope curve is not wise, because many time-consuming calculations still need to be processed,such as solving the exact intersection points between splines.To accelerate,we turn to filter out the very few effective spline segments from all boundaries. Fortunately, the strong convex hull property of NURBS spline can be used for just that. Let us describe the details.

    Fig. 7 Solving shortest tool length of a known tool path (left) based on all compressed boundaries (right). (Compressed edges are colored by black and compressed silhouettes are colored by pink).

    3.1. AATB of boundary spline segment

    3.1.1. Boundary spline segment

    For a compressed boundary B(μ), it can be defined by a 3thdegree NURBS curve and expressed as:

    where μ is the curve parameter, {Pi} are the n + 1 control points, {wi} are the weights, {Ni,3(μ)} are the 3th-degree B-spline basis functions defined on the nonperiodic (and nonuniform) knot vector U = {0,0,0,0,u4,u5,...,un,1,1,1,1}. Then B(u) will be split into n + 1 spline segments according to its knot vector. As a result, the spline segment is the basic unit for filtering. According to the strong convex hull property of NURBS curve,30that is, if u ∈[ui,ui+1), then B(μ)lies within the convex hull of the four control points Pi-3,Pi-2,...,Pi(i.e., the control polygon), as shown in Fig. 8.

    For a spline segment, let rminand rmaxbe the left and right limits of its control polygon, respectively. Given a specific THS, if both Rkk) are met, then the control polygon crosses the influenced ranges of the kth to the qth THS profile segments. Similarly, it means this control polygon affects tool length by restricting the kth to the qth segment of THS profile. For the convenience of presentation, the former spline segment is marked as type I while the latter is marked as type II.

    Fig. 8 A spline segment and its control polygon.

    3.1.2. AATb

    If we directly use the irregular control polygons to filter the effective spline segments,it will require a large number of polygon intersection operations,which is too complicated and inefficient.The Axis-aligned Bounding Box(AABB)of the control polygon is a useful and efficient filtering tool although it only works for the multi-cylindrical THS(such as the one in Fig.15(left)).To apply to the more generic THS which contains conical tool handles, inspired by the idea of AABB, here we propose a new concept which can give a conservative range of tool length for a spline segment: the axis-aligned tool length box(AATB) of control polygon.

    For an AATB,rminand rmaxare also its left and right limits,respectively. Here, let lmaxand lminbe the maximum and minimum values of the shortest tool length determined by the control polygon of spline segment, respectively. Referring Fig. 9(left), lmaxand lminare also the upper and lower limits of the AATB, respectively. To avoid visual confusion, an AATB is displayed by a four-color box (the left, right, upper and lower lines are colored by pink,cyan,orange,and blue,respectively.)Obviously, AATB is determined by the positions of control polygon relative to Γ(L). Noticeably, without a specific THS, AATB cannot exist independently. For a known spline segment, rminand rmaxare easy to obtain while lminand lmaxneed to be calculated by the following steps:

    Step 1: According to Eq. (18), calculate the shortest tool lengths respectively determined by the four control points of the control polygon.

    Step 2: If the spline segment is type II, then perform this step; otherwise, go to Step 3. Since the convex hull of type II spline segment crosses the multiple influenced ranges of the THS profile segment, it is necessary to additionally solve all the intersection points between the control polygon and the transboundary lines of THS profile (i.e., r=Rk,k=1,2,...,K).Referring Fig.9(left),the two upward-pointing triangles on the control polygon are the intersection points.Next,calculate the shortest tool lengths determined by the left and right limit points of intersection points respectively. The hcoordinates of four downward-pointing triangles in Fig.9(left)reflect the mapped tool lengths.

    Fig.9 AATB of control polygon(left)and filtering the effective AATBs (right).

    Step 3: The minimum and maximum values of all shortest tool lengths solved in Step 1 and Step 2 are lminand lmax,respectively.

    The h-coordinate of AATB can be used to directly and consistently reflect the shortest tool lengths determined by different spline segments. In other words, relative to all AATBs, it can be considered that the complicated THS profile(no matter the multi-cylinders,the multi-tapers,or the mixture of the first two) has been transformed to be a horizontal line h=Lshortest,as shown in Fig.9(right).As a result,AATB can be utilized to quickly filter out the effective spline segments.

    By carefully observing Fig.9(left),we can get such a conclusion: the height range of AATB of a type II spline segment is conservative. Specially, if the spline segment is located within the (k+1)th to the qth influenced ranges, then the lmaxof the current AATB is ‘‘false high” because the ‘‘real” upper limit of the spline segment should be lower than lmax.However,to further carefully distinguish this special case is unnecessary and skipped here,as it does not miss the correct effective spline segments.

    For the example shown in Fig. 7, the AABBs and AATBs of all boundary spline segments are illustrated in Fig. 10(left)and Fig.10(right),respectively.Although both of them are displayed by the four-color boxes, their shapes may be different.Specifically, only for the type I spline segment in the special case, i.e., the spline segment is completely located in a single influenced range of cylindrical THS profile segment, dose its AATB and its AABB have the same shape, for example, the rightmost spline segment in Fig. 9(right). Certainly, this conclusion can also be inferred from Eq. (18). Compared with AABB, AATB is more applicable because it is applicable for multi-taper shaped THS (including multi-cylinder, multitaper, and cylinder-taper mixed).

    3.2. Filtering effective spline segments based on AATBs

    The filtering criterion for effective boundary spline segments is:let Lmaxbe the maximum lminof all AATBs,if an AATB satisfies lmax>Lmax,then the corresponding spline segment is effective, that is, it is possible to ultimately determine the shortest tool length. According to the above criteria, all effective AATBs and spline segments for both type I and type II will be quickly filtered out. Note: it may conservatively filter out the extremely few ‘‘false” effective type II spline segments.

    By applying the above criterion,the effective AATBs of the example shown in Fig. 7 are illustrated in Fig. 11. From the zoomed figure, it can be observed that all horizontal orange lines (i.e., lmaxof effective AATBs) are higher than the highest horizontal blue line h=Lmax(i.e., the maximum of lminfor all effective AATBs). Remarkably, the number of effective AATBs is very small (i.e., 28) compared to that one (i.e.,5560)in Fig.10.By tracing the forming process of these effective spline segments,we can find they come from six boundary curves(i.e.,three edges and three silhouettes),and the number only accounts for 3.53% of the total number of all boundary curves. Further tracing, these six curves come from only three tool postures. Evidently, only very few tool postures and very small areas on the check surface are effective to ultimately determine the shortest tool length. Besides that, it is conceivable that the number of effective AATBs and spline segments will stay at a small level even though the number of tool postures may increase dramatically(the results of several tests can confirm it). In other words, this acceleration method is welltargeted and robust.

    Fig. 10 AABBs (left) and AATBs (right) of all control polygons.

    Fig. 11 All effective spline segments and AATBs.

    3.3. Calculating shortest tool length based on effective spline segments

    Fortunately, the required accuracy of the shortest tool length in industrial production is relatively low, e.g., 0.1 mm. As the identified effective spline segments are very few and short,hence, we adopted the simple sampling strategy already mentioned in Section 2.3 to solve the exact shortest tool length according to Eq. (18). Certainly, the distance between two neighboring sample points should not exceed the required accuracy. The mapped tool lengths of all sample points are shown by the black solid curves in Fig. 12. As a result, the h-coordinate of the highest mapped point (e.g., the red star in Fig. 12) is the shortest tool length Lshortest.

    We can draw several interesting observations from Fig.12.One is the tool length has a jump near the line r=R2. The main reason is,for the given THS profile in Fig.7(right),there is a height difference H1between the end-point of the non-vertical part of 1th segment and the start-point of the non-vertical part of the 2th segment. Another interesting phenomenon is that the shape of mapped tool lengths who satisfies R2

    Fig. 12 All mapped tool lengths of sample points on effective spline segments.

    3.4.Tracing critical interference information in machining scene

    Noticeably, the highest mapped point in the 2D rh plane can reflect the critical interference information in the 3D machining scene with the shortest tool length. By carefully tracing the formation process of the highest point,the following items are clear:the riskiest tool posture,the critical interference position, which check surface and which area (edge or the interior silhouette)participate in forming the critical interference point,which part of the THS restricts the shortest tool length, etc.Tracing and analyzing the example in Fig. 7, it is clear that the 1th tool posture is the most dangerous.Besides,the bottom disc of tool handle is just touching the unique check surface at one edge, which restricts the shortest tool length. These judgments are consistent with the actual state in the 3D machining scene, as illustrated in Fig. 13. Of course, the traced interference information can be applied to adjust those bottleneck tool axes (just like the approach in16) or modify the shape of the tool handle to shorten the global shortest tool length.

    4. Results of experimental examples and discussion

    Fig. 13 Tracing riskiest tool posture and critical interference position (red star).

    To verify and evaluate the accuracy and applicability of the proposed numerical algorithm for efficiently solving the shortest tool length,based on the software platforms of Unigraphics V18.0 and MATLAB 2017b and the integrated development environment of Microsoft Visual Studio 2015, we developed a prototype software system written by C++ language and applied it to a typical open blisk, which is a complex channel part used in aero-engine. Referring Fig. 14, the length of the blade is about 160 mm, a conical ball-end cutter BR2C2 (the ball radius is 2 mm and the taper is 2°) is adopted for finish machining of the pressure surface of the blade in a five-axis machine tool. To improve machining efficiency and quality,it is desirable to use the shortest tool length to enhance the tool stiffness as much as possible. The required accuracy of the shortest tool length is 0.01 mm.The tool axis is changed along the tool path under the rule of the constant lead and title angles, i.e., 18° and 0°, respectively. There are five tool paths colored by cyan near the root of the blade (For the convenience of display, the width of cutting line is set to be larger than normal value). Without counting the non-cutting motions, a total of 196 interference-free tool postures participate in the calculation of the shortest tool length.In this example,two check surfaces are considered,i.e.,the pressure surface Spand the suction surface Ss. Both of them have 27×34=918 control points. When searching the silhouette curves, the step Lsis set to 0.01 mm.

    Fig. 15 shows two typical THSs. The left one adopted a multi-cylinder hydraulic tool handle, denoted as #1-THS.According to the characterization model of THS in Section 2.1,the geometric parameters of its profile are {(9.5 mm, 54 mm,90 deg); (20 mm, 49 mm, 90 deg); (32 mm, 16 mm, 90 deg);(87.5 mm, 45 mm, 90 deg); (109 mm, 41 mm, 90 deg);(1000 mm, 1000 mm, 90 deg)}. Similarly, the right one used a multi-taper heat-shrinkable too handle (there is a conical segment in THS profile), denoted as #2-THS. Its parameters are{(9 mm, 0 mm, 90 deg); (13 mm, 94 mm, 3 deg); (32 mm,26 mm, 90 deg); (48 mm, 30 mm, 90 deg); (1000 mm,1000 mm, 90 deg)}. Note: the height of the 1st segment of#2-THS is 0 mm. It can be seen that both THS profiles are complicated and can reflect the actual THS tightly.Our prototype software is implemented on an average computer with an Intel(R) Core (TM) i5-2300 M CPU with a 2.80 GHz processor.

    4.1. Calculation results

    4.1.1. Boundaries of all compressed surfaces

    Fig. 14 A five-axis tool trajectory for machining open blisk.

    Fig. 15 Two THSs. #1-THS (left); #2-THS (right).

    The compressed surfaces formed by Spand Ssrelative to all tool postures are shown in Fig. 16(left) (or Fig. 16(right)).For all boundaries, there are 2×196×4=1568 edges and 313 silhouettes (all are type I). By adjusting the tool length,the two THS profiles are just contacting with these boundaries,as shown in Fig.16.Intuitively,only very few boundary curves are close to the THS profiles.

    4.1.2. Effective boundary spline segments and AATBs

    According to the geometric parameters of two THSs, the AATBs of all boundary spline segments are illustrated in Fig. 17. The identified effective AATBs and spline segments are shown in Fig. 18. For #1-THS, only 9/62627=0.014%of AATBs are effective; likewise, the percentage of #2-THS is 20/62627=0.032%. The proportion of effective AATBs is also the proportion of effective spline segments. Such a small proportion validates the motivation and advantages of the AATB proposed by us. In other words, by segmenting the THS profile and the compressed boundaries and overall scanning all AATBs, the proposed method can accurately and quickly discover all effective boundary spline segments, which make the next sampling approach for solving the shortest tool length more focused.This acceleration method avoids wasting precious computing resources on sampling the majority boundaries which are impossible to contribute to the shortest tool length. Therefore, the complete solution presented can be easily applied in industrial production.

    By tracing all effective spline segments of #1-THS, they come from only 4 compressed boundaries,which only account for 4/1568=0.026% of all boundaries. Similarly, all effective spline segments of #2-THS are from 7 original compressed boundaries and the proportion is 7/1568=0.045%. When the tool length is minimum, the relationships between THS profile and the effective compressed boundaries are illustrated in Fig. 19. The results in this figure are consistent with that in Fig. 16 but with fewer compressed boundaries.

    4.1.3. Shortest tool length and running time

    The mapped tool lengths at the sample points on the effective boundary spline segments are shown by the black solid curves in Fig. 20. The solved shortest tool length of #1-THS is 71.1194 mm. By tracing the most risk point which ultimately determines the global shortest tool length,it can be found that the 171th tool posture which is located on the deepest tool path is most dangerous. When #1-THS move to that posture with the shortest tool length, the tool handle has partially entered the blisk channel, meanwhile, the bottom edge of the second cylinder in tool handle is just interfering with the blade tip(edge) of Ss,as shown in Fig. 21 (left). Similarly, the obtained shortest tool length of #2-THS is 59.1993 mm. Different from#1-THS, the 9th tool posture which is located on the shallowest tool path is the most dangerous.At this time,the side of the first cone in tool handle is just contacting with a corner point of Ss, as illustrated in Fig. 21 (right). For the two different THSs in this example, the most dangerous tool posture and the critical interference position are greatly different, which is mainly related to the tool axes and the specific geometric parameters of tool holding system.

    Fig.16 THS profiles are just contacting with all compressed boundaries.#1-THS(left);#2-THS(right).(Compressed edges are colored by black and compressed silhouettes are colored by pink).

    Fig. 17 All AATBs. #1-THS (left); #2-THS (right).

    Fig. 18 Effective AATBs (four-color boxes) and spline segments (black solid curves). #1-THS (left); #2-THS (right). (Dark red solid lines are transboundary lines of THS profile).

    Fig. 19 All effective compressed boundaries. #1-THS (left); #2-THS (right). (Compressed edges are colored by black and compressed silhouettes are colored by pink).

    Fig.20 Tool lengths(black solid curves)mapped by all effective boundary spline segments and global shortest tool length(red star).#1-THS (left); #2-THS (right).

    The components of running time for two THSs are listed in Table 1. It can be observed that the procedure of solving the compressed boundaries takes the most calculation time for both THSs, about 25.991/28.354=91.67%. Reducing this part of time is definitely our focus for improving the algorithm in the future,e.g.,using the parallel algorithms to track the silhouette curves.

    4.2. Comparison with other methods

    Fig. 21 Critical interference position (red star) under the riskiest tool posture with the shortest tool length: #1-THS (left); #2-THS(right).

    Table 1 The components of running time for two THSs.

    To further verify the calculation quality and efficiency of the proposed algorithm, we also compared it with two commonly used methods for calculating the shortest tool length, which are the method of calculating the minimum distance (Method 1) and the method of fully-discretization of the check surfaces(Method 2). The procedure of Method 1 is: determining whether THS body interferes with the channel by calculating the minimum distance between the two entities; using the dichotomy method to get the precise minimum tool length of each tool posture (the range of tool length is (0 mm,300 mm), and the iteration accuracy is 0.01 mm); as a result,the global minimum tool length is the goal. The program is implemented in C ++ language and serial mode. The procedure of Method 2 is: sampling the entire check surfaces and determining the minimum tool length of each tool posture based on these sample points, and then the global minimum tool length of all postures is the target. To evaluate the algorithm more objectively, this method is implemented in MATLAB language and parallel mode. As the size of both check surfaces is about 159.76 mm×89.38 mm and the discrete density is 0.1 mm×0.1 mm, then the total number of sample points is 2×1598×894=2.86×106(Note: if the density becomes 0.01 mm×0.01 mm, the number of samples is as high as 2.86×108, an ordinary computer can not handle them well). Considering that the parallel Method 2 performed not well for a THS with multiple tapers (because it has to frequently judge which taper of THS is affected by the sample point,which seriously affects the efficiency of parallel computing), hence a simple and conservative THS (recorded as #3-THS) which is composed of two cylinders was adopted in this test, as shown in Fig. 22(left). Its parameters are {(25 mm,50 mm, 90 deg); (109 mm, 102 mm, 90 deg); (1000 mm,1000 mm, 90 deg)}. The critical interference state under the riskiest tool posture with the shortest tool length is illustrated in Fig. 22(right).

    For Method 1 and Method 2, both of them need to calculate the exact shortest tool length of each tool posture one by one. In Fig. 23, we only present the results calculated by Method 2. The calculation results and running time of three methods are compared in Table 2.

    From Fig.23 and Table 2,it can be seen that the results of the proposed algorithm and the two comparison methods are consistent, such as the shortest tool length and the critical interference position.However,under the same hardware conditions, it is gratifying that the running time of the proposed method is only 27.293/139.552=19.56% of Method 1(though their accuracies are the same), which is 27.293/110.235=24.76% of the Method 2 (though the accuracy of the former is higher than the latter).

    Method 1 costs most of the time to calculate the shortest tool length of each tool posture, as a result, solving the global shortest tool length is inefficient. Although the proposed method also needs to obtain the compressed boundaries of each tool posture one by one, overall filtering the effective AATBs can quickly pick out the few effective spline segments.Compared with Method 1, our method can not only improve the calculation efficiency but also guarantee the accuracy.

    As can be seen that the calculation results of three methods in this example are very close, theoretically, the accuracy in Method 2 should be lower than that of the other two methods.We thought this is just a coincidence for this test. Because the critical interference position is located at the corner of the check surface (see Fig. 22(right)). Even though Method 2 has a low discrete accuracy (i.e., 0.01 mm×0.01 mm), the corners of the check surface are always sampled in general. If the critical interference position was not sampled,the low accuracy of Method 2 will be leaked. Due to using the strategy of fullydiscretization, Method 2 requires the high-performance computer hardware.When the required precision is extremely high,the computational cost of Method 2 will increase by the order of square. However, for the proposed algorithm, the number of CP points on a silhouette curve will increase linearly,which will result in a linear increase in calculation load. Compared with Method 2, our method has lower requirements for computer hardware and is easier to apply in industrial production.

    In short,it can be concluded that the proposed method for calculating the shortest tool length can not only improve computing efficiency but also can ensure accuracy.

    Fig. 22 #3-THS (left) and critical interference position (red star) under riskiest tool posture (right).

    Fig.23 Required tool lengths for each tool posture(blue and orange bars represent critical interference with Sp and Ss,respectively)and global shortest tool length in update (green stepped curve).

    Table 2 Calculation results and running time of three methods.

    4.3. Application

    4.3.1. Optimizing tool handle for shortest tool length

    For the tool paths in Fig. 14, here the conical ball-end cutter BR2C2 was respectively clamped by eight different tool handles (#1-THS, #2-THS, and #3-THS are the same with those in Section 4.1 and 4.2; the physical entity of #5 tool handle is given in Fig. 3(left)) with their shortest tool length Lshortest,as illustrated in Fig. 24, the other geometries (such as spindle)in THS are the same as #1-THS; then we compared their tool stiffnesses. The tool material is the common K44 cemented carbide, whose Young’s modulus is 490 GPa. The calculation formula for the stiffness of a generic cutter in13(which has been physically verified) was used herein and the results are shown in Fig. 25.

    Fig. 24 Eight different tool handles with shortest tool length.(Red stars are critical interference positions).

    Fig. 25 Comparison of shortest tool lengths and tool stiffnesses of eight different THSs.

    It can be seen that different tool handles may have different shortest tool lengths and may result in huge differences in tool stiffness. For example, the ratio of the shortest tool length of#4-THS to that one of #3-THS is 40.12 mm/127.29 mm=31.52%, while the ratio of stiffness is(12.7986×105N·m-1)/(2.8489×105N·m-1)=4.49. #4-THS is the optimal one that can significantly improve the tool stiffness and will be beneficial for the CNC machining efficiency and quality in practical (This idea has been applied to industrial manufacturing by authors). Besides, another interesting phenomenon is that the shortest tool lengths of #4 and #8 are the same while their tool handle shapes are different. Certainly, obtaining the maximum combined stiffness of cutter and tool handle may be a better option, it will be explored in the future.

    Fig. 26 Transitional-taper cutter with double cylinders.

    4.3.2. Selecting appropriate transitional-taper cutter

    In the multi-axis machining of complex channel parts, as described previously, using a full-conical cutter (the tapered shank is assumed to be long enough)to replace a cylinder cutter can enhance the stiffness.However,compared with the latter, the applicability of the former is poor for other channels due to its fixed taper,which may easily cause the tapered shank near the ball head to interfere with the machined surface.Besides, the conical cutter is not easy to be resharpened once its cutting edges are worn, because the entire tapered shank has to be ground. A compromise choice that may slightly reduce the stiffness is the transitional-taper cutter with double cylinders,as shown in Fig.26.The cylinder with diameter D1is the only clamping part of the tool,and the other cylinder with diameter D0is allowed to be ground while the transitional cone(the taper is βcand the cone height is hc) is not. In industrial production, using this type of cutter is more economic thanks to its better applicability and improvement of stiffness.

    Only constrained by the current THS, the shortest tool length L1can be calculated by the proposed algorithm. However,it still cannot guarantee that all resharpened transitionaltaper cutters in the tool storehouse are available due to their non-standard h0.For a specific transitional-taper tool,whether it interferes with the check geometries still needs further judgment.Specifically,the judgment method is as follows:only taking the clamping cylinder (with the diameter D0) and the transition cone as a new THS, a new shortest tool length L2can be calculated by our algorithm. Let Δhcbe the height difference between the critical interference point (such as the red star in Fig. 26) and the bottom disc of the transition cone.Then h0≥L2-Δhccan guarantee this transitional-taper cutter is interference-free with the check geometries. Therefore,the shortest clamping length of this cutter, i.e., the L in Fig. 26, is equal to max{L1,h0+hc|h0≥L2-Δhc}, which can guarantee neither the THS nor the transitional-taper cutter interferes with the check geometries.

    This method of selecting and clamping the transitionaltaper cutter has been successfully applied in practice to choose the strongest cutter.For example,here we take#4-THS(whose shortest tool length for the full conical cutter is minimum in Section 4.3.1) for the tool paths in Fig. 14, then L1is 40.12 mm. Two different transitional-taper cutters (i.e.,BR2C45D8 and BR2C4D7), the cylindrical cutter (BR2) and the full-conical cutter(BR2C2) in the tool storehouse are four alternative cutters.Fig.27 shows the geometric parameters and the shortest clamping lengths of the candidate cutters.

    Fig. 27 Four different candidate cutters with their shortest clamping lengths.

    Clearly,among the four candidate cutters,the transitionaltaper cutter BR2C45D8 is the most suitable for the current tool path because its stiffness is the maximum, i.e.,13.03×105N·m-1. Comparing the stiffness of the cylindrical cutter BR2 and the transitional-taper cutter BR24C45D8,the latter is 5.57 times higher than the former. Hence, using a transitional-taper cutter is indeed an available solution to improve the cutter stiffness in practice although a careful calculation is required.

    5. Conclusion and future work

    This paper proposed a novel method to quickly solve the shortest tool length of a known tool trajectory in multiaxis machining. The main idea is not to sample the entire check surfaces but to find out the few necessary curves which may determine the shortest tool length. By rotating and compressing the 3D check surface relative to all tool postures,the boundary curves (i.e., the edges and the newly formed silhouette curves in compressing) of the 2D compressed surfaces are the necessary elements. To efficiently solve silhouette curves, a tracking-based numerical algorithm is given. Based on the proposed characterization model for the multi-taper shaped THS, a model for solving the shortest tool length is constructed according to the critical interference relationship between the THS profile and all compressed boundaries. To speed up, a new concept of AATB is proposed, which can give a conservative range of tool length for a boundary spline segment. By scanning all AATBs,the very few effective boundary spline segments that may ultimately determine the shortest tool length were then filtered out. This method makes the solution for solving the shortest tool length more focused and efficient. The verification examples confirmed it.

    About the future research on the subject, we plan to identify the few bottleneck tool postures for the shortest tool length(e.g., the tool postures near the 174th one in Fig. 23) in advance based on the relationship between the channel features and tool postures. We believe it will further improve the efficiency of solving the shortest tool length. Moreover,extracting the local effective regions on check surface may be another acceleration method worth trying, especially for large parts,such as the leading edge protection cap of the composite blade. Besides, another interesting work is to improve the acceleration method to apply to the spline-shaped THS profile.Certainly, upgrading the prototype with advanced programming techniques (e.g., multithreading and GPU acceleration)will make the proposed algorithm more practical in industrial production, which is also another task for authors in future.

    Declaration of Competing Interest

    The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

    Acknowledgment

    The authors gratefully acknowledge the support of National Science and Technology Major Project of China (No. JPPTKF2016).We also greatly appreciate Dr.Luo Ming for the linguistic assistance.

    a级片在线免费高清观看视频| 美女午夜性视频免费| 亚洲欧美色中文字幕在线| www.av在线官网国产| 男女免费视频国产| 一本色道久久久久久精品综合| 日本av手机在线免费观看| 少妇 在线观看| 亚洲美女搞黄在线观看| 国产精品久久久久久久久免| 国产片特级美女逼逼视频| 欧美久久黑人一区二区| 我的亚洲天堂| 青春草亚洲视频在线观看| 日韩精品有码人妻一区| 少妇人妻 视频| 最新的欧美精品一区二区| 国产精品久久久久久人妻精品电影 | 欧美精品亚洲一区二区| 欧美另类一区| 制服诱惑二区| 亚洲一码二码三码区别大吗| 精品酒店卫生间| 熟女av电影| 国产亚洲精品第一综合不卡| 午夜精品国产一区二区电影| 女人被躁到高潮嗷嗷叫费观| 国产精品久久久久久久久免| 建设人人有责人人尽责人人享有的| 麻豆精品久久久久久蜜桃| 亚洲三区欧美一区| 亚洲精品国产av蜜桃| 亚洲国产欧美网| 最近最新中文字幕免费大全7| 男人操女人黄网站| 亚洲av中文av极速乱| 欧美在线一区亚洲| 精品国产乱码久久久久久男人| 午夜免费观看性视频| 亚洲第一av免费看| 国产在线视频一区二区| 一级毛片黄色毛片免费观看视频| 777米奇影视久久| 无限看片的www在线观看| 乱人伦中国视频| 大香蕉久久网| 日日撸夜夜添| 亚洲自偷自拍图片 自拍| 一级黄片播放器| 另类亚洲欧美激情| 日韩av不卡免费在线播放| 久久久国产一区二区| 一本大道久久a久久精品| 亚洲伊人色综图| 午夜福利在线免费观看网站| 爱豆传媒免费全集在线观看| 最近2019中文字幕mv第一页| 亚洲精品久久成人aⅴ小说| 满18在线观看网站| 亚洲精品乱久久久久久| 国产精品一区二区精品视频观看| 国产精品香港三级国产av潘金莲 | 超色免费av| 伦理电影免费视频| 亚洲欧美一区二区三区国产| kizo精华| 99精品久久久久人妻精品| 日本爱情动作片www.在线观看| 丰满饥渴人妻一区二区三| 欧美老熟妇乱子伦牲交| 免费观看a级毛片全部| 精品少妇久久久久久888优播| 免费看av在线观看网站| 亚洲人成77777在线视频| av片东京热男人的天堂| 久久久久网色| 久久精品aⅴ一区二区三区四区| 交换朋友夫妻互换小说| av一本久久久久| 亚洲精品第二区| 亚洲视频免费观看视频| 一区在线观看完整版| 一二三四中文在线观看免费高清| 新久久久久国产一级毛片| 大陆偷拍与自拍| 国产麻豆69| 精品国产国语对白av| av在线app专区| 国产精品一区二区在线观看99| 日韩制服骚丝袜av| 在线 av 中文字幕| 欧美老熟妇乱子伦牲交| 精品一区二区三区四区五区乱码 | 九色亚洲精品在线播放| 亚洲精品美女久久久久99蜜臀 | 久久97久久精品| 国产精品久久久久久久久免| 香蕉丝袜av| 大码成人一级视频| 韩国av在线不卡| 久久亚洲国产成人精品v| 在线观看免费日韩欧美大片| 午夜免费鲁丝| 精品亚洲成a人片在线观看| 欧美日韩一区二区视频在线观看视频在线| 欧美国产精品一级二级三级| 亚洲专区中文字幕在线 | 欧美日韩精品网址| 欧美激情 高清一区二区三区| 亚洲av日韩在线播放| 啦啦啦在线观看免费高清www| 亚洲一区中文字幕在线| 亚洲成av片中文字幕在线观看| 美女主播在线视频| 精品卡一卡二卡四卡免费| 亚洲成人免费av在线播放| 51午夜福利影视在线观看| 2018国产大陆天天弄谢| 国产精品久久久久久精品古装| 大香蕉久久网| 午夜激情久久久久久久| 99热全是精品| 亚洲av电影在线进入| 国产精品亚洲av一区麻豆 | 国产精品成人在线| 精品国产国语对白av| 日本欧美视频一区| 国产精品一国产av| 女人爽到高潮嗷嗷叫在线视频| 久久国产精品男人的天堂亚洲| 欧美精品人与动牲交sv欧美| 国产一区二区 视频在线| 老司机靠b影院| videosex国产| 欧美少妇被猛烈插入视频| 天堂俺去俺来也www色官网| 欧美日本中文国产一区发布| 亚洲国产av新网站| 国产黄色视频一区二区在线观看| 岛国毛片在线播放| 最近最新中文字幕免费大全7| 欧美日韩精品网址| 1024香蕉在线观看| 69精品国产乱码久久久| 久久影院123| 伦理电影免费视频| 两个人免费观看高清视频| 精品一品国产午夜福利视频| 人人妻人人爽人人添夜夜欢视频| 黄色视频不卡| 国产精品av久久久久免费| 黄频高清免费视频| 精品一区二区免费观看| 日韩中文字幕视频在线看片| 成人18禁高潮啪啪吃奶动态图| 人人澡人人妻人| 国产人伦9x9x在线观看| 国产亚洲最大av| 精品少妇内射三级| 精品少妇内射三级| 人妻人人澡人人爽人人| 亚洲欧美激情在线| 久久久久久免费高清国产稀缺| 久久久久久人人人人人| 女人高潮潮喷娇喘18禁视频| 久久天堂一区二区三区四区| 青草久久国产| 亚洲自偷自拍图片 自拍| 蜜桃在线观看..| av福利片在线| 日韩 亚洲 欧美在线| 在线精品无人区一区二区三| av又黄又爽大尺度在线免费看| 99热全是精品| 久久精品国产亚洲av高清一级| 青春草亚洲视频在线观看| 丝瓜视频免费看黄片| 亚洲七黄色美女视频| 国产精品一区二区在线不卡| 熟妇人妻不卡中文字幕| 亚洲精品视频女| 亚洲精品久久午夜乱码| 国产淫语在线视频| 亚洲欧洲国产日韩| 99热全是精品| 亚洲av日韩在线播放| 欧美日韩视频高清一区二区三区二| 久久精品国产a三级三级三级| 男女下面插进去视频免费观看| 人人澡人人妻人| 在线观看人妻少妇| 亚洲人成电影观看| 99久久综合免费| 一区福利在线观看| 极品人妻少妇av视频| 97精品久久久久久久久久精品| 看免费成人av毛片| 国产精品免费大片| 国产乱人偷精品视频| 亚洲国产欧美一区二区综合| 久久韩国三级中文字幕| 蜜桃在线观看..| 少妇人妻久久综合中文| 悠悠久久av| 国产成人精品久久久久久| 综合色丁香网| 国产色婷婷99| 亚洲欧洲国产日韩| 亚洲成人一二三区av| 你懂的网址亚洲精品在线观看| av一本久久久久| 国产乱人偷精品视频| 亚洲人成电影观看| 在线 av 中文字幕| 最近中文字幕2019免费版| 久久久久精品国产欧美久久久 | www.精华液| 国产爽快片一区二区三区| 永久免费av网站大全| 日韩成人av中文字幕在线观看| 免费黄频网站在线观看国产| av福利片在线| 在线观看www视频免费| 精品国产露脸久久av麻豆| 人人澡人人妻人| 免费少妇av软件| 最黄视频免费看| 老司机影院成人| 欧美激情极品国产一区二区三区| 中文字幕最新亚洲高清| 久久婷婷青草| 可以免费在线观看a视频的电影网站 | 天天操日日干夜夜撸| 黑人欧美特级aaaaaa片| 一区福利在线观看| 成年女人毛片免费观看观看9 | 日韩电影二区| 国产97色在线日韩免费| 纯流量卡能插随身wifi吗| 天天操日日干夜夜撸| 毛片一级片免费看久久久久| 日本黄色日本黄色录像| 搡老乐熟女国产| av在线老鸭窝| 久久久精品国产亚洲av高清涩受| 国产无遮挡羞羞视频在线观看| 亚洲七黄色美女视频| 卡戴珊不雅视频在线播放| 免费看不卡的av| av片东京热男人的天堂| 韩国精品一区二区三区| 欧美日韩福利视频一区二区| 精品国产超薄肉色丝袜足j| 免费观看av网站的网址| 热re99久久国产66热| 亚洲,欧美精品.| 成年女人毛片免费观看观看9 | 秋霞伦理黄片| 国产成人免费观看mmmm| 欧美 日韩 精品 国产| 日韩一本色道免费dvd| 丝袜美足系列| 日韩,欧美,国产一区二区三区| 黑人猛操日本美女一级片| 日日啪夜夜爽| 亚洲国产av影院在线观看| 国产日韩欧美视频二区| 美女大奶头黄色视频| 久久精品国产综合久久久| 深夜精品福利| 两个人看的免费小视频| 国产成人午夜福利电影在线观看| 男女边吃奶边做爰视频| 一区福利在线观看| 高清av免费在线| 男女午夜视频在线观看| 国产欧美日韩综合在线一区二区| 免费看av在线观看网站| 国产精品成人在线| 亚洲欧美日韩另类电影网站| 欧美在线一区亚洲| 五月天丁香电影| 老司机影院毛片| 青青草视频在线视频观看| 欧美黑人精品巨大| 国产激情久久老熟女| 精品亚洲乱码少妇综合久久| 国产国语露脸激情在线看| 啦啦啦啦在线视频资源| 国产极品天堂在线| 亚洲欧美一区二区三区国产| 母亲3免费完整高清在线观看| 国产在视频线精品| 精品国产露脸久久av麻豆| 啦啦啦在线免费观看视频4| 国产精品国产av在线观看| 狂野欧美激情性xxxx| 狠狠精品人妻久久久久久综合| 精品免费久久久久久久清纯 | 蜜桃在线观看..| 中文天堂在线官网| 女人精品久久久久毛片| 国产99久久九九免费精品| 另类精品久久| 精品少妇黑人巨大在线播放| 免费观看a级毛片全部| 久久久久人妻精品一区果冻| 国产精品久久久久久精品电影小说| 亚洲国产毛片av蜜桃av| av.在线天堂| 18禁裸乳无遮挡动漫免费视频| 欧美 亚洲 国产 日韩一| 午夜免费观看性视频| 国产 一区精品| 亚洲精品国产av成人精品| 黄色视频不卡| 国产成人免费无遮挡视频| 午夜福利一区二区在线看| 亚洲av福利一区| 十八禁高潮呻吟视频| 别揉我奶头~嗯~啊~动态视频 | 男人爽女人下面视频在线观看| 精品酒店卫生间| 久久久久久久精品精品| 欧美日本中文国产一区发布| 国精品久久久久久国模美| 高清视频免费观看一区二区| 天堂8中文在线网| 国产av一区二区精品久久| 国产精品人妻久久久影院| 老熟女久久久| 欧美人与性动交α欧美精品济南到| 韩国高清视频一区二区三区| 只有这里有精品99| 久久久久久免费高清国产稀缺| 国产极品天堂在线| 激情视频va一区二区三区| 日韩一本色道免费dvd| 曰老女人黄片| 亚洲精品aⅴ在线观看| 欧美人与性动交α欧美软件| 中文字幕人妻丝袜制服| 欧美国产精品va在线观看不卡| 香蕉丝袜av| 女人爽到高潮嗷嗷叫在线视频| 国产精品一区二区在线观看99| 一二三四中文在线观看免费高清| 国产一区二区三区综合在线观看| 国产亚洲欧美精品永久| 欧美精品亚洲一区二区| 中文字幕人妻丝袜制服| 国产精品亚洲av一区麻豆 | avwww免费| 在线观看免费日韩欧美大片| 91精品伊人久久大香线蕉| 亚洲精品中文字幕在线视频| 叶爱在线成人免费视频播放| 日本一区二区免费在线视频| av在线播放精品| 三上悠亚av全集在线观看| 九草在线视频观看| 国产国语露脸激情在线看| 99九九在线精品视频| 丝袜喷水一区| 成人毛片60女人毛片免费| 啦啦啦在线观看免费高清www| 国产精品免费视频内射| 国产色婷婷99| 久久99精品国语久久久| 少妇猛男粗大的猛烈进出视频| 美女午夜性视频免费| 亚洲色图 男人天堂 中文字幕| 激情视频va一区二区三区| 成人手机av| 韩国高清视频一区二区三区| 观看av在线不卡| 欧美亚洲日本最大视频资源| 最黄视频免费看| 久久久久久久国产电影| 亚洲美女视频黄频| 国产乱人偷精品视频| 亚洲av男天堂| 亚洲国产看品久久| 激情五月婷婷亚洲| 人妻一区二区av| 欧美久久黑人一区二区| 19禁男女啪啪无遮挡网站| 亚洲 欧美一区二区三区| 捣出白浆h1v1| 国产免费一区二区三区四区乱码| 免费女性裸体啪啪无遮挡网站| 99久国产av精品国产电影| 免费观看av网站的网址| 97人妻天天添夜夜摸| 男女边摸边吃奶| 免费黄色在线免费观看| 国产精品偷伦视频观看了| 亚洲精品中文字幕在线视频| 日本wwww免费看| 国产在线一区二区三区精| 99热网站在线观看| 美女中出高潮动态图| 午夜久久久在线观看| 久久久精品免费免费高清| 国产人伦9x9x在线观看| 欧美黄色片欧美黄色片| 久久国产亚洲av麻豆专区| 成人免费观看视频高清| 人妻人人澡人人爽人人| 精品福利永久在线观看| 亚洲精品成人av观看孕妇| 丝袜人妻中文字幕| 飞空精品影院首页| 亚洲av男天堂| 一区福利在线观看| 久久久久久久久免费视频了| 亚洲精品美女久久久久99蜜臀 | 久久人人97超碰香蕉20202| 菩萨蛮人人尽说江南好唐韦庄| 无遮挡黄片免费观看| xxx大片免费视频| 国产淫语在线视频| 9色porny在线观看| 国产免费现黄频在线看| 婷婷色综合www| 亚洲国产精品一区二区三区在线| 亚洲精品成人av观看孕妇| 如日韩欧美国产精品一区二区三区| 在线亚洲精品国产二区图片欧美| 中国三级夫妇交换| 久热这里只有精品99| 精品第一国产精品| 成人亚洲精品一区在线观看| www日本在线高清视频| 欧美精品人与动牲交sv欧美| 亚洲精品中文字幕在线视频| 国产成人91sexporn| av在线老鸭窝| 久久ye,这里只有精品| 十八禁网站网址无遮挡| 女人爽到高潮嗷嗷叫在线视频| 成人毛片60女人毛片免费| 久久人人爽人人片av| 日本色播在线视频| 乱人伦中国视频| xxx大片免费视频| 久久天堂一区二区三区四区| 丁香六月欧美| 中文字幕另类日韩欧美亚洲嫩草| 黄色视频在线播放观看不卡| 日本色播在线视频| 亚洲成人免费av在线播放| 欧美中文综合在线视频| 人妻一区二区av| 中文字幕高清在线视频| 51午夜福利影视在线观看| 在线免费观看不下载黄p国产| 91国产中文字幕| 熟妇人妻不卡中文字幕| 成人三级做爰电影| 欧美黑人欧美精品刺激| 嫩草影院入口| 最近最新中文字幕大全免费视频 | 国产精品二区激情视频| 国产成人免费无遮挡视频| 狠狠婷婷综合久久久久久88av| 下体分泌物呈黄色| 欧美精品一区二区免费开放| 乱人伦中国视频| 欧美日韩亚洲国产一区二区在线观看 | 人人妻人人澡人人爽人人夜夜| 亚洲精品久久成人aⅴ小说| 七月丁香在线播放| 久久久久久久大尺度免费视频| 午夜福利乱码中文字幕| 亚洲av男天堂| 高清在线视频一区二区三区| 视频区图区小说| 婷婷成人精品国产| av免费观看日本| 嫩草影视91久久| 色婷婷av一区二区三区视频| 国产乱人偷精品视频| 精品久久久久久电影网| 999久久久国产精品视频| 一级毛片 在线播放| 精品一品国产午夜福利视频| 久久97久久精品| 亚洲第一区二区三区不卡| 麻豆乱淫一区二区| 人人妻人人澡人人看| 午夜影院在线不卡| 天堂8中文在线网| 99国产综合亚洲精品| a 毛片基地| 别揉我奶头~嗯~啊~动态视频 | 国产 精品1| 哪个播放器可以免费观看大片| 高清在线视频一区二区三区| av一本久久久久| 成人免费观看视频高清| 国产 精品1| 看非洲黑人一级黄片| 国产日韩欧美在线精品| 午夜免费鲁丝| 免费女性裸体啪啪无遮挡网站| 国产日韩欧美视频二区| 色婷婷av一区二区三区视频| 国产不卡av网站在线观看| 午夜日韩欧美国产| 日韩av在线免费看完整版不卡| 男人舔女人的私密视频| 夫妻午夜视频| 90打野战视频偷拍视频| 亚洲欧美精品综合一区二区三区| 视频区图区小说| 日韩一区二区视频免费看| xxxhd国产人妻xxx| 黄频高清免费视频| 999精品在线视频| 日韩中文字幕视频在线看片| 国产精品二区激情视频| 热99国产精品久久久久久7| 免费人妻精品一区二区三区视频| 国产野战对白在线观看| 久久人人爽av亚洲精品天堂| 精品国产露脸久久av麻豆| 国产精品国产三级国产专区5o| 99久久精品国产亚洲精品| 国产av一区二区精品久久| 国产成人a∨麻豆精品| 天天躁夜夜躁狠狠久久av| 成年av动漫网址| 免费av中文字幕在线| 亚洲成色77777| 侵犯人妻中文字幕一二三四区| 久久久久精品久久久久真实原创| 午夜免费观看性视频| 午夜福利视频精品| 男女边摸边吃奶| 日日啪夜夜爽| 欧美久久黑人一区二区| 欧美xxⅹ黑人| 中文字幕制服av| 男女边吃奶边做爰视频| 欧美最新免费一区二区三区| 狂野欧美激情性bbbbbb| 国产成人91sexporn| 国产伦理片在线播放av一区| 一本大道久久a久久精品| 香蕉国产在线看| 国产xxxxx性猛交| 校园人妻丝袜中文字幕| 男女下面插进去视频免费观看| 国产欧美日韩一区二区三区在线| 国产精品三级大全| 免费高清在线观看日韩| 天堂俺去俺来也www色官网| 中文字幕色久视频| 亚洲在久久综合| 热99国产精品久久久久久7| 99精品久久久久人妻精品| 午夜免费观看性视频| 99久久99久久久精品蜜桃| 这个男人来自地球电影免费观看 | 久久久久人妻精品一区果冻| 免费观看人在逋| 老司机靠b影院| www.精华液| 久久国产精品大桥未久av| 丰满饥渴人妻一区二区三| 成年女人毛片免费观看观看9 | 男的添女的下面高潮视频| 老熟女久久久| 欧美久久黑人一区二区| 亚洲国产精品一区二区三区在线| 国产爽快片一区二区三区| 成人影院久久| 精品国产乱码久久久久久男人| 国产成人欧美| 免费久久久久久久精品成人欧美视频| 国产成人精品在线电影| 性色av一级| 高清av免费在线| 亚洲人成77777在线视频| 国产精品一二三区在线看| 中文精品一卡2卡3卡4更新| 国产欧美亚洲国产| 精品少妇一区二区三区视频日本电影 | 嫩草影院入口| 91国产中文字幕| 国产淫语在线视频| 久久精品国产综合久久久| 如日韩欧美国产精品一区二区三区| 精品视频人人做人人爽| netflix在线观看网站| 一区二区三区精品91| 美女中出高潮动态图| 9191精品国产免费久久| 青草久久国产| 精品一品国产午夜福利视频| 日韩欧美精品免费久久| 日本色播在线视频| 国产精品久久久人人做人人爽| 又黄又粗又硬又大视频| 美女主播在线视频| 蜜桃在线观看..| 亚洲精品日本国产第一区| 精品一区二区免费观看| 亚洲精品美女久久久久99蜜臀 | 黄片小视频在线播放| a级毛片黄视频| 国产成人免费观看mmmm| 黄色一级大片看看| 青春草国产在线视频|