Simulation method for ball-end milling with coupling of ultrasonic vibration and trochoidal modulation
By combining the Z-map algorithm with ultrasonic vibration and cocycloidal modulation for ball end mill simulation, the problem of balancing large computational load and simulation accuracy in ultra-precision milling is solved, achieving efficient simulation and accurate surface quality assessment, and reducing experimental costs.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- GUANGDONG UNIV OF TECH
- Filing Date
- 2026-03-25
- Publication Date
- 2026-07-21
AI Technical Summary
Existing simulation models for ultra-precision milling have high computational demands, making it difficult to balance simulation accuracy and efficiency. Furthermore, they are costly to conduct experiments and lack effective simulation prediction and parameter optimization methods.
A simulation method for ball end mills based on Z-map algorithm coupling ultrasonic vibration and cocycloidal modulation is adopted. By introducing tool, machine tool and ultrasonic vibration parameters and combining ROI search mechanism, the simulation model is optimized to simulate the cutting process of the tool on the workpiece surface, so as to achieve high-precision and high-efficiency simulation.
It achieves high-precision surface morphology simulation, improves simulation efficiency, provides comprehensive quantitative evaluation of surface quality, supports process parameter optimization, and reduces experimental costs.
Smart Images
Figure CN122433231A_ABST
Abstract
Description
Technical Field
[0001] This application relates to the technical field of ultra-precision machining, specifically to a simulation method for ball end milling coupled with ultrasonic vibration and cocycloidal modulation. Background Technology
[0002] The continuous development and improvement of modern industrial systems and the rapid pace of technological revolution have significantly impacted the iterative upgrading of industrial technologies across various sectors. Among these advancements, the continuous progress of ultra-precision machining technology is a crucial branch of modern industrial development and technological progress. Ultra-precision machining technology encompasses a wide range of fields, including aerospace, optical instruments, electronic equipment, and medical machinery. Within the field of ultra-precision machining, milling technology is a key component. Milling involves the coordinated operation of ultra-precision machine tools, diamond cutting tools, and the workpiece. In some cutting-edge research areas, it also involves ultrasonic vibration-assisted machining modules and the integrated analysis of optics, electrochemistry, thermodynamics, and mechanics. Ultra-precision milling technology is characterized by extremely high precision, flexible and convenient machining, and a wide range of applications. However, in real-world milling operations, various experimental limitations, such as chipping of expensive diamond tools due to improper operation, unnecessary consumption of experimental consumables, trial-and-error experiments, lengthy and ineffective experiments, and the influence of external factors, all contribute to high experimental costs. Therefore, it is particularly necessary to construct an accurate simulation mathematical model that describes the interaction between the tool and the workpiece before the formal experiment begins, and to use simulation software to analyze the machining morphology, evaluate the expected surface quality, optimize process parameters, and reasonably reduce experimental costs.
[0003] Throughout history, researchers have tirelessly strived to develop large-scale simulation models that closely approximate experimental data, provide more accurate simulation predictions, adapt to complex machining processes, integrate multiple modules, and achieve high computational efficiency in the field of ultra-precision machining. These models can be broadly categorized into two schools of thought: one focusing on the advanced development of geometric modeling techniques, and the other on physical / mechanical modeling. Regarding geometric modeling techniques, researchers proposed the solid construction geometry method in the 1970s. Based on set theory and Boolean algebra, they considered the workpiece and the tool sweep volume as mathematical "solids." The cutting process was described as a Boolean difference operation (A - B) between the workpiece solid (A) and the tool sweep volume (B). However, its limitations were also apparent: the computational load increased exponentially with the number of cutting steps, leading to memory overflow and limiting its application to early, simple experimental verification. Later, with the iterative development of computer technology, simulation modeling techniques also became increasingly sophisticated. One such technique is the multi-normal depth element method based on the Z-map algorithm. Canadian researchers continuously improved upon the Z-map algorithm, creatively overturning the Z-map's argument that "a point can only store one height" based on ray tracing and discrete sampling theory. They stored multiple entry and exit points in the X, Y, Z, and even more normal directions, making it applicable to various surface data acquisitions and improving the fitting accuracy of the simulation model. Another school of thought is the physical / mechanical modeling approach, which has also made significant progress in related technologies. For example, the shear plane theory and slip line field theory proposed in the 1940s laid the foundation for modern cutting mechanics, and their related models can be used to quickly estimate cutting forces. However, the high computational and storage requirements limit its application areas. On the other hand, Japanese researchers have also achieved great success in finite element analysis (FEM). They were the first to establish an energy-based model for predicting tool wear and cutting temperature, while a team at Ohio State University made significant contributions to the application of FEM in metal forming and cutting. However, both geometric modeling and physical / mechanical modeling simulations place significant demands on computational resources. Furthermore, the computational load is closely related not only to the number of cutting passes and steps but also to simulation accuracy. Balancing and selecting among these three factors is a key consideration for current simulation modeling technologies. Therefore, using the Z-map algorithm for ultra-precision milling simulation modeling requires optimizing and reducing computational load while ensuring simulation accuracy and a high number of cutting steps, minimizing bloated calculations caused by ineffective cutting and excessive pursuit of precision. Summary of the Invention
[0004] The purpose of this invention is to provide a simulation method for ultra-precision milling of ball end mills coupled with ultrasonic vibration and cocycloidal modulation. By designing a milling simulation mathematical model based on MATLAB software programming, the method predicts and evaluates the overall surface morphology, texture, microstructure trajectory and state of the workpiece by modulating the machining parameters and trajectory morphology. At the same time, the method also performs numerical standard evaluation and comparison of the simulated surface morphology quality, in order to improve machining parameters and optimize machining process, so as to assist in actual experimental machining.
[0005] To solve the aforementioned technical problems, the basic concept of the technical solution adopted by this invention is as follows: First, the milling simulation mathematical model is implemented based on the Z-map algorithm, and its steps include introducing tool parameters, machine tool motion parameters, ultrasonic vibration parameters, and cocycloidal modulation parameters; second, kinematic calculation preprocessing needs to be completed; then, the Z-map algorithm is used to initialize the two-dimensional surface mesh, which discretizes the workpiece into a 2D mesh (XY plane). After the workpiece surface mesh initialization is completed, the milling simulation stage begins. This part is the core loop of the program, used to simulate the tool's interaction with the surface. The physical cutting process of the workpiece is based on the Z-map algorithm mentioned above. It involves multiple passes and time steps (k) loops, calculates the tool center position, defines the ROI region and searches for the cutting area, binds the contact point and the sweep arc, and updates the surface height Z_surf in real time, taking the minimum value to simulate material removal. In this step, the introduction of the ROI search region and the selection of a certain length of cutting arc segment that only removes material are mainly used to avoid invalid cutting, optimize and reduce the amount of calculation, and ensure simulation accuracy. Finally, the data post-processing outputs the topography map, XY cross-sectional height map and roughness parameter map.
[0006] A simulation method for ultra-precision milling of ball end mills coupled with ultrasonic vibration and cocycloidal modulation includes the following steps: Step S1: First, set the simulation parameters, including tool parameters, machine tool motion parameters, ultrasonic vibration parameters, cocycloidal modulation parameters, and workpiece settings. The tool parameters consist of the tool radius R. The machine tool motion parameters include spindle speed n, feed rate Vf, depth of cut ap, and step size s. The ultrasonic vibration parameters include frequency F, X-direction amplitude Ax, Z-direction amplitude Az, and cocycloidal modulation parameters A... mod The modulation frequency is f_mod, and the workpiece setting parameters include workpiece width lx, workpiece length ly, and grid precision grid_size; Step S2, next, in order to complete the linear motion plus disturbance process of the tool path and simulate planar milling, while avoiding redundant calculations within the loop and improving efficiency, kinematic calculation preprocessing needs to be completed. The parameters of this part are composed of spindle angular velocity w, total simulation time T_total, time step dt, time range vector, and cutting pass. Logically, the spindle speed n is converted into angular velocity w. At the same time, the total simulation time T_total is determined by the feed rate and workpiece length. Regarding the time step calculation part, the continuous motion is discretized into a time sequence t, combined with rotation resolution (1° / step) to capture the details of the tool sweep, while considering the coverage of the workpiece width by multiple machining passes. Step S3, then, uses the Z-map algorithm to initialize the two-dimensional surface mesh. This part discretizes the workpiece into a 2D mesh (XY plane). First, the workpiece is divided into tens of thousands of mesh points with a precision of micrometer level according to the preset mesh precision. Then, each mesh point is stored as the initial height (initially 0, indicating that the workpiece has not been cut). The motion trajectory of the tool in the time series is simulated to calculate the contact point between the tool and the workpiece, and the mesh height value is updated in time. The minimum value is taken to simulate material removal. Step S4: After initializing the workpiece surface mesh, the milling simulation stage begins. This part is the core loop of the program, used to simulate the physical cutting process of the tool on the workpiece. It is based on the Z-map algorithm mentioned above, through multiple passes and time steps (k) loops, calculates the tool center position, defines the ROI region and searches for the cutting region, binds the contact point and the sweep arc, and updates the surface height Z_surf in real time, taking the minimum value to simulate material removal. Logically, this part considers tool rotation, vibration, modulation disturbance and high-precision mesh, and optimizes efficiency based on the defined ROI region to avoid full mesh traversal. The whole process simulates Boolean cutting: the tool envelope sweeps the workpiece, leaving the minimum height. First, pre-calculation of vectors is performed, where the tool amplitude kinematics formulas for the X and Z directions are as follows:
[0007] Meanwhile, to ensure the tool simulates a cocycloidal path for sweeping cutting and feed, the kinematic formula is as follows:
[0008] The vibration displacements Vx and Vz are pre-calculated using the above kinematic formulas, and the tool motion trajectory, i.e., the modulation displacement X, is also calculated. mod and Y mod (X cosine forms cocycloidal perturbation, Y linear feed + sine perturbation), Z baseline Z base Maximum offset O max ROI search radius R search(Contact projection + 1.5 margin + step offset to ensure complete coverage), while using an angle vector with initial phase displacement to simulate the cutting edge rotation phase. The pre-calculated vectors complete the preparation work, which not only defines the tool movement trajectory but also improves the simulation speed. Then, the outer loop iterates through the passes, setting the X baseline. base To ensure the cumulative cutting across passes is calculated, the inner loop iterates through time steps (k), continuously updating the iteration count to simulate planar milling across passes. First, the current tool center Xc (pass offset X) is calculated. base +X modulation +X amplitude), Yc (Y modulation), Zc (baseline Z) base (+Z vibration), the formula for the tool coordinate center is as follows:
[0009] By specifying and adjusting the tool's coordinate position on the workpiece in real time using the tool coordinate center formula described above, we can ensure that the toolpath accurately and completely completes the cutting process, guaranteeing that the planned path meets the expected requirements and improving machining accuracy. The next step is ROI optimization, which is based on the tool center coordinates and the search radius R. search Define the tool boundary range, determine the left, right, top, and bottom rectangular index boundaries, and then judge the grid points in the potential cutting area. If a point is close to the tool or within the rectangular boundary, it will be cut; otherwise, it will be skipped. This can significantly improve the program's running speed and avoid invalid cuts. The index formulas for the X and Y directions (points within the search range) are given below. The formula for calculating the minimum X index is shown below:
[0010] X min It is the minimum mesh index calculation for the ROI in the X direction (counting from 1), which is the current tool center minus R. search Divide by the grid precision (grid_size, denoted by g here), mapping the continuous space to a discrete grid, defining the left boundary of the ROI, and then determining the formula for calculating the maximum X index as follows:
[0011] X max This is the calculation of the maximum grid index of the ROI in the X direction (Nx is the total number of grid columns in the X direction), which adds R to the current tool center. search Dividing by the grid precision (grid_size, denoted by g here), we map the continuous space to a discrete grid, define the right boundary of the ROI, and ensure that the index covers the entire search radius and is within the search range. Below is the formula for calculating the minimum Y index direction:
[0012] Y min It is the minimum mesh index calculation for the ROI in the Y direction (counting from 1), which is the current tool center minus R. search Divide by the grid precision (grid_size, denoted by g here), mapping the continuous space to a discrete grid, defining the lower boundary of the ROI, and then determining the formula for calculating the maximum Y index as follows:
[0013] Y max This is the calculation of the maximum mesh index of the ROI in the Y direction (Ny is the total number of mesh rows in the Y direction), which is the current tool center minus R. search Divide by the grid precision (grid_size, denoted by g here), mapping the continuous space to a discrete grid, defining the upper boundary of the ROI, and forming the ROI range (rectangular region) using the above formula. Then, search for potential cutting points within the ROI range. The next step is to perform angle sweep determination. The core idea and process is to simulate the tool cutting only within a finite arc segment during rotation. Check whether a point on the workpiece is located within the sector angle covered by the tool's current rotation and whether the distance is within the tool radius. If so, the surface height value is updated; otherwise, it is skipped. The process involves calculating the azimuth angle of the grid point relative to the tool center, as shown in the following formula:
[0014] Next, the current blade angle is obtained, and the formula is as follows:
[0015] Next, calculate the minimum angle between the current cutting edge angle and the azimuth angle of the cutting center (to ensure the minimum path), and then calculate the angle through which the tool rotates (sweep step length). The formula is as follows:
[0016] Then, it is determined whether the grid point is within the cutting region (based on both angle and distance conditions). If it is, it is cut, and the surface height value is updated. The formula is as follows:
[0017] Here Z surf The value is based on the real-time surface update of the Z-map algorithm, which is the principle of material removal. Its logic is that the initial workpiece height is defined as 0, and the height value is negative after the tool enters the workpiece. If the current tool height value Z... tool The height Z of the workpiece surface surf If it's even lower (deeper), it means the tool has cut to this point, meaning the Z-axis needs to be updated. surfIf the tool height is higher than the surface height (indicating the tool is moving idly, or the cut has been made deeper), then Z is not updated. surf ; Step S5, finally, involves post-processing of the data: section extraction, roughness calculation, and result visualization. First, the X and Y section height data are extracted, and the formula is shown below:
[0018] Next, the roughness values Ra, Rq, and Rt are calculated. The formula for the arithmetic mean deviation Ra is as follows:
[0019] Where N represents the sampling points, the next step is to calculate Rq. The formula for the root mean square deviation Rq is as follows:
[0020] Finally, the total height Rt of the outline is calculated, using the formula shown below:
[0021] After completing the above calculations, the values are converted to micrometers, and the surface roughness is calculated to measure and objectively reflect the surface quality of the entire plane. The formulas for the surface arithmetic mean deviation Sa and the surface root mean square deviation Sq are as follows:
[0022] Finally, the processed and calculated data are visualized, including the overall 3D surface morphology, cross-sectional profiles (X-direction and Y-direction cross-sections), and a summary of roughness parameters (Ra, Rt, Rq, Sa, Sq).
[0023] In some optional embodiments, in step S1, regarding the workpiece setting part, the lengths of the workpiece lx and ly are set to 1×1mm respectively, wherein in order to ensure computational efficiency and balance the microscopic detail accuracy requirements, the grid precision grid_size is 0.003; Simultaneously determine the cycloidal modulation parameter A. mod The value is 0.029, the modulation frequency f_mod is 28, the number of cutting edges is single-edged by default, and the tool radius is 2mm.
[0024] In some optional embodiments, the kinematic calculation preprocessing part in step S2 converts the spindle speed n into angular velocity w, as shown in the following formula:
[0025] Since the spindle speed is measured in r / min, dividing by 60 converts it to rad / s. Additionally, the total simulation time T_total is determined by the feed rate Vf and the workpiece length ly, as shown in the following formula:
[0026] Here, the unit is also uniformly s; Next is the calculation of the time step dt, the formula of which is as follows:
[0027] dt corresponds to the time step (s) of 1° rotation. The rotational speed per second is n / 60, and the degree per second is (rotational speed) × 360°. Therefore, the time for 1° is 1 / (degree / s). Finally, the cutting pass is defined, and its formula is as follows:
[0028] The above formula is used to discretize the continuous motion of the tool into a gridded point set for computation, in order to capture the microscopic sweeping details of the cutting edge.
[0029] In some optional embodiments, in step S3, the workpiece surface machining accuracy is determined by both the mesh accuracy and the workpiece size, as shown in the following formula:
[0030] In this context, grid precision is represented by g, which stands for grid_size. lx is the workpiece width, and ly is the workpiece length.
[0031] In some optional embodiments, in step S4, regarding the surface quality characteristic parameters Sa and Sq, MN is represented as a surface matrix with M×N sampling points.
[0032] The present invention has the following beneficial effects: High-precision surface morphology simulation – This method employs the Z-map algorithm combined with micron-level mesh precision to discretize the workpiece surface, and fully considers multi-dimensional factors such as tool parameters, machine tool motion parameters, ultrasonic vibration parameters, and cocycloidal modulation parameters. By simulating the precise motion trajectory of the tool over time, especially by introducing ultrasonic vibration (X / Z direction amplitude) and cocycloidal modulation disturbances, it can more realistically reproduce the complex physical cutting process in ultra-precision milling, capture microscopic tool sweep details, and thus generate a high-fidelity overall 3D surface morphology.
[0033] High computational and simulation efficiency – To address the massive computational burden of high-precision simulations, this method employs an efficient ROI search mechanism. Instead of traversing the entire mesh during cutting calculations, a rectangular boundary is dynamically defined based on the current tool center position and the search radius Rsearch. The program only evaluates and updates mesh points within this ROI that meet the cutting conditions. This local search strategy significantly reduces redundant calculations, dramatically improving program speed and simulation efficiency.
[0034] Precise physical cutting logic simulation – This method not only performs fine discretization in space but also discretizes continuous motion into a time series in time, combining rotational resolution for calculation. In particular, it introduces "angle sweep determination" logic to simulate the finite arc cutting characteristics of the tool during rotation. By calculating the difference between the azimuth angle of a grid point relative to the tool center and the current cutting edge angle, and determining whether it falls within the included angle and radius of the tool sweep's fan shape, it achieves accurate simulation of the Boolean cutting process (simulating material removal by taking the lowest height). This mechanism ensures that the simulation is not only geometrically correct but also conforms to the discontinuous cutting characteristics of actual physical cutting. Comprehensive Surface Quality Quantitative Assessment – After simulation, this method provides a complete data post-processing workflow, capable of automatically extracting cross-sectional height data in the X and Y directions from the simulation results. It not only calculates traditional line roughness parameters (… Furthermore, the surface roughness parameters were calculated. This comprehensive evaluation system, encompassing both two-dimensional contours and three-dimensional surfaces, coupled with visualized charts (3D morphology, cross-sectional contour diagrams, etc.), provides objective and quantitative data support for process parameter optimization and machining quality prediction, helping to intuitively reflect the machining quality of the entire plane. Attached Figure Description
[0035] To more clearly illustrate the technical solutions in the embodiments of this application, the accompanying drawings used in the description of the embodiments will be briefly introduced below. Obviously, the accompanying drawings described below are only some embodiments of this application. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.
[0036] Figure 1 This is a schematic diagram of the method flow of the present invention.
[0037] Figure 2 This is a schematic diagram of workpiece mesh discretization using the Z-map algorithm.
[0038] Figure 3 It is a milled surface topography diagram.
[0039] Figure 4This is a summary diagram of the XY cross-sectional height and roughness parameters. Detailed Implementation
[0040] The present application will now be described in further detail with reference to the accompanying drawings and embodiments. It should be particularly noted that the following embodiments are for illustrative purposes only and do not limit the scope of the application. Similarly, the following embodiments are only some, not all, embodiments of the present application, and all other embodiments obtained by those skilled in the art without inventive effort are within the scope of protection of the present application.
[0041] The terms "first," "second," and "third" used in the embodiments of this application are for descriptive purposes only and should not be construed as indicating or implying relative importance or implicitly specifying the number of indicated technical features. Therefore, a feature defined as "first," "second," or "third" may explicitly or implicitly include at least one of that feature. In the description of this application, "multiple" means at least two, such as two, three, etc., unless otherwise explicitly specified. All directional indications (such as up, down, left, right, front, back, etc.) in the embodiments of this application are only used to explain the relative positional relationships and movement of components in a specific posture (as shown in the figures). If the specific posture changes, the directional indication will also change accordingly. The terms "comprising" and "having," and any variations thereof, in the embodiments of this application are intended to cover non-exclusive inclusion. For example, a process, method, system, product, or device that includes a series of steps or units is not limited to the listed steps or units, but may optionally include steps or units not listed, or may optionally include other steps or components inherent to these processes, methods, products, or devices.
[0042] In this document, the term "embodiment" means that a particular feature, structure, or characteristic described in connection with an embodiment may be included in at least one embodiment of this application. The appearance of this phrase in various places throughout the specification does not necessarily refer to the same embodiment, nor is it a separate or alternative embodiment mutually exclusive with other embodiments. It will be explicitly and implicitly understood by those skilled in the art that the embodiments described herein can be combined with other embodiments.
[0043] This embodiment will be further described with reference to the accompanying drawings, and its specific measures are as follows: This application provides a simulation method for ultra-precision milling of ball end mills coupled with ultrasonic vibration and cycloidal modulation, including the following steps: Step S1: Specific parameter examples are provided below, along with accompanying drawings, to illustrate the milling model. Simulation parameters are set, including tool parameters (tool radius 2mm), machine tool motion parameters (spindle speed n = 7000 r / min, feed rate Vf = 200 mm / min, depth of cut ap = 0.05 mm, step size s = 0.2 mm), ultrasonic vibration parameters (frequency F = 5000 Hz, X-direction amplitude Ax = 0.002 mm), and cocycloidal modulation parameters (modulation amplitude A...). mod The value is 0.029, the modulation frequency f_mod is 28, and finally the workpiece settings include workpiece width lx of 1mm, workpiece length ly of 1mm, and grid precision grid_size of 0.003mm. Step S2: After setting the above tool parameters, kinematic calculation preprocessing is performed. This part consists of spindle angular velocity w, total simulation time T_total, time step dt, time range vector, and cutting passes. In this part, the spindle angular velocity w is calculated based on the value of the spindle speed n in step S1. At the same time, the total simulation time T_total is determined by the feed rate and workpiece length. Regarding the time step calculation, the continuous motion is discretized into a time sequence t, combined with the rotation resolution (1° / step) to capture the details of the tool sweep, while also considering the coverage of the workpiece width by multiple machining passes. Step S3, then, involves initializing the two-dimensional surface mesh according to the Z-map algorithm, as shown in the attached figure. Figure 2 As shown, the specific content is to discretize the workpiece into a 2D mesh (XY plane). First, according to the preset mesh accuracy grid_size of 0.003mm in step S1, the workpiece is divided into mesh points with micron-level accuracy. Then, each mesh point is stored as the initial height (initially 0, indicating that the workpiece has not been cut). The motion trajectory of the tool in the time series is simulated to calculate the contact point between the tool and the workpiece, and the mesh height value is updated in time. The minimum value is taken to simulate material removal. Step S4: After initializing the workpiece surface mesh, the milling simulation stage begins. This part is the core loop of the program, used to simulate the physical cutting process of the tool on the workpiece. It is based on the Z-map algorithm mentioned above, using multiple passes and time steps (k) to calculate the tool center position, define the ROI region and search for the cutting region, bind the contact point and the sweep arc, and update the surface height in real time, taking the minimum value to simulate material removal. Its flowchart is attached. Figure 1As shown, logically, this part considers tool rotation, vibration, modulation disturbances, and high-precision meshes. It also optimizes efficiency based on defined ROI regions to avoid full mesh traversal. The entire process simulates Boolean cutting: the tool envelope sweeps across the workpiece, leaving a minimum height. The first step is to perform pre-calculated vector work, where the tool amplitude kinematics formulas for the X and Z directions are shown below:
[0044] Meanwhile, to ensure the tool simulates a cocycloidal path for sweeping cutting and feed, the kinematic formula is as follows:
[0045] The vibration displacements Vx and Vz are pre-calculated using the above kinematic formulas, and the tool motion trajectory, i.e., the modulation displacement X, is also calculated. mod and Y mod (X cosine forms cocycloidal perturbation, Y linear feed + sine perturbation), Z baseline Z base Maximum offset O max ROI search radius R search (Contact projection + 1.5 margin + step offset to ensure complete coverage), while using an angle vector with initial phase displacement to simulate the cutting edge rotation phase. The pre-calculated vectors complete the preparation work, which not only defines the tool movement trajectory but also improves the simulation speed. Then, the outer loop iterates through the passes, setting the X baseline. base To ensure the cumulative cutting across passes is calculated, the inner loop iterates through time steps (k), continuously updating the iteration count to simulate planar milling across passes. First, the current tool center Xc (pass offset X) is calculated. base +X modulation +X amplitude), Yc (Y modulation), Zc (baseline Z) base (+Z vibration), the formula for the tool coordinate center is as follows:
[0046] By specifying and adjusting the tool's coordinate position on the workpiece in real time using the tool coordinate center formula described above, we can ensure that the toolpath accurately and completely completes the cutting process, guaranteeing that the planned path meets the expected requirements and improving machining accuracy. The next step is ROI optimization, which is based on the tool center coordinates and the search radius R. search Define the tool boundary range, determine the left, right, top, and bottom rectangular index boundaries, and then judge the grid points in the potential cutting area. If a point is close to the tool or within the rectangular boundary, it will be cut; otherwise, it will be skipped. This can significantly improve the program's running speed and avoid invalid cuts. The index formulas for the X and Y directions (points within the search range) are given below. The formula for calculating the minimum X index is shown below:
[0047] X min It is the minimum mesh index calculation for the ROI in the X direction (counting from 1), which is the current tool center minus R. search Divide by the grid precision (grid_size, denoted by g here), mapping the continuous space to a discrete grid, defining the left boundary of the ROI, and then determining the formula for calculating the maximum X index as follows:
[0048] X max This is the calculation of the maximum grid index of the ROI in the X direction (Nx is the total number of grid columns in the X direction), which adds R to the current tool center. search Dividing by the grid precision (grid_size, denoted by g here), we map the continuous space to a discrete grid, define the right boundary of the ROI, and ensure that the index covers the entire search radius and is within the search range. Below is the formula for calculating the minimum Y index direction:
[0049] Y min It is the minimum mesh index calculation for the ROI in the Y direction (counting from 1), which is the current tool center minus R. search Divide by the grid precision (grid_size, denoted by g here), mapping the continuous space to a discrete grid, defining the lower boundary of the ROI, and then determining the formula for calculating the maximum Y index as follows:
[0050] Y max This is the calculation of the maximum mesh index of the ROI in the Y direction (Ny is the total number of mesh rows in the Y direction), which is the current tool center minus R. search Divide by the grid precision (grid_size, denoted by g here), mapping the continuous space to a discrete grid, defining the upper boundary of the ROI, and forming the ROI range (rectangular region) using the above formula. Then, search for potential cutting points within the ROI range. The next step is to perform angle sweep determination. The core idea and process is to simulate the tool cutting only within a finite arc segment during rotation. Check whether a point on the workpiece is located within the sector angle covered by the tool's current rotation and whether the distance is within the tool radius. If so, the surface height value is updated; otherwise, it is skipped. The process involves calculating the azimuth angle of the grid point relative to the tool center, as shown in the following formula:
[0051] Next, the current blade angle is obtained, and the formula is as follows:
[0052] Next, calculate the minimum angle between the current cutting edge angle and the azimuth angle of the cutting center (to ensure the minimum path), and then calculate the angle through which the tool rotates (sweep step length). The formula is as follows:
[0053] Then, it is determined whether the grid point is within the cutting region (based on both angle and distance conditions). If it is, it is cut, and the surface height value is updated. The formula is as follows:
[0054] Here Z surf The value is based on the real-time surface update of the Z-map algorithm, which is the principle of material removal. Its logic is that the initial workpiece height is defined as 0, and the height value is negative after the tool enters the workpiece. If the current tool height value Z... tool The height Z of the workpiece surface surf If it's even lower (deeper), it means the tool has cut to this point, meaning the Z-axis needs to be updated. surf If the tool height is higher than the surface height (indicating the tool is moving idly, or the cut has been made deeper), then Z is not updated. surf Its 3D topographic image is attached. Figure 3 As shown; Step S5, finally, involves post-processing of the data: section extraction, roughness calculation, and result visualization. First, the X and Y section height data are extracted, and the formula is shown below:
[0055] Next, the roughness values Ra, Rq, and Rt are calculated. The formula for the arithmetic mean deviation Ra is as follows:
[0056] Where N represents the sampling points, the next step is to calculate Rq. The formula for the root mean square deviation Rq is as follows:
[0057] Finally, the total height Rt of the outline is calculated, using the formula shown below:
[0058] After completing the above calculations, the values are converted to micrometers, and the surface roughness is calculated to measure and objectively reflect the surface quality of the entire plane. The formulas for the surface arithmetic mean deviation Sa and the surface root mean square deviation Sq are as follows:
[0059] Finally, the processed and calculated data are visualized, including the overall 3D surface morphology, cross-sectional profiles (X-axis and Y-axis sections), and a summary of roughness parameters (Ra, Rt, Rq, Sa, Sq). An example of the results is shown in the attached figure. Figure 4 As shown.
[0060] Furthermore, in the kinematic calculation preprocessing part of step S2, the spindle speed n is converted into angular velocity w, and the formula is as follows:
[0061] Since the spindle speed is measured in r / min, dividing by 60 converts it to rad / s. Additionally, the total simulation time T_total is determined by the feed rate Vf and the workpiece length ly, as shown in the following formula:
[0062] Here, the unit is also uniformly s; Next is the calculation of the time step dt, the formula of which is as follows:
[0063] dt corresponds to the time step (s) of 1° rotation. The rotational speed per second is n / 60, and the degree per second is (rotational speed) × 360°. Therefore, the time for 1° is 1 / (degree / s). Finally, the cutting pass is defined, and its formula is as follows:
[0064] The above formula is used to discretize the continuous motion of the tool into a gridded point set for computation, in order to capture the microscopic sweeping details of the cutting edge.
[0065] Furthermore, in step S3, the surface machining accuracy of the workpiece is determined by both the mesh accuracy and the workpiece size, as shown in the following formula:
[0066] In this context, grid precision is represented by g, which stands for grid_size. lx is the workpiece width, and ly is the workpiece length.
[0067] Furthermore, in step S4, regarding the surface quality characteristic parameters Sa and Sq, MN represents a surface matrix with M×N sampling points.
[0068] This embodiment provides a simulation method for ultra-precision milling of ball end mills coupled with ultrasonic vibration and cocycloidal modulation. Based on the Z-map algorithm, it introduces cocycloidal modulation parameters and simulates the normal milling process of ball end mills through the above example. Compared with ordinary milling simulation models, to study the micron-level microstructure morphology and characteristic laws, this model introduces ultrasonic vibration parameter modules in the X and Z directions. By modulating the vibration amplitude and response frequency in different directions, different morphology maps can be obtained. Furthermore, while ensuring micron-level accuracy, this model introduces boundary index conditions (ROI regions) and defines cutting arc segments to search for potential cutting points, which can significantly improve computational efficiency, save computation time, and avoid invalid cutting. Finally, in terms of data visualization, the morphology variation law and the characteristics of valleys and ridges are analyzed using 3D morphology X and Y-direction cross-sectional height maps. Roughness parameters (Ra, Rt, Rq, Sa, Sq) are also summarized to illustrate the quality assessment status of different machining parameters on the morphology trajectory, providing assistance for subsequent optimization of machining processes and adjustment and improvement of machining parameters.
[0069] The above description is only a part of the embodiments of this application and does not limit the scope of protection of this application. Any equivalent device or equivalent process transformation made based on the content of this application specification and drawings, or direct or indirect application in other related technical fields, are similarly included in the patent protection scope of this application.
Claims
1. A simulation method for ball end milling coupled with ultrasonic vibration and cocycloidal modulation, characterized in that, Includes the following steps: Step S1: First, set the simulation parameters, including tool parameters, machine tool motion parameters, ultrasonic vibration parameters, cocycloidal modulation parameters, and workpiece settings. The tool parameters consist of the tool radius R. The machine tool motion parameters include spindle speed n, feed rate Vf, depth of cut ap, and step size s. The ultrasonic vibration parameters include frequency F, X-direction amplitude Ax, Z-direction amplitude Az, and cocycloidal modulation parameters A... mod The modulation frequency is f_mod, and the workpiece setting parameters include workpiece width lx, workpiece length ly, and grid precision grid_size; Step S2, next, in order to complete the linear motion plus disturbance process of the tool path and simulate planar milling, while avoiding redundant calculations within the loop and improving efficiency, kinematic calculation preprocessing needs to be completed. The parameters of this part are composed of spindle angular velocity w, total simulation time T_total, time step dt, time range vector, and cutting pass. Logically, the spindle speed n is converted into angular velocity w. At the same time, the total simulation time T_total is determined by the feed rate and workpiece length. Regarding the time step calculation part, the continuous motion is discretized into a time sequence t, combined with rotation resolution (1° / step) to capture the details of the tool sweep, while considering the coverage of the workpiece width by multiple machining passes. Step S3, then, uses the Z-map algorithm to initialize the two-dimensional surface mesh. This part discretizes the workpiece into a 2D mesh (XY plane). First, the workpiece is divided into tens of thousands of mesh points with a precision of micrometer level according to the preset mesh precision. Then, each mesh point is stored as the initial height (initially 0, indicating that the workpiece has not been cut). The motion trajectory of the tool in the time series is simulated to calculate the contact point between the tool and the workpiece, and the mesh height value is updated in time. The minimum value is taken to simulate material removal. Step S4: After initializing the workpiece surface mesh, the milling simulation stage begins. This part is the core loop of the program, used to simulate the physical cutting process of the tool on the workpiece. It is based on the Z-map algorithm mentioned above, through multiple passes and time steps (k) loops, calculates the tool center position, defines the ROI region and searches for the cutting region, binds the contact point and the sweep arc, and updates the surface height Z_surf in real time, taking the minimum value to simulate material removal. Logically, this part considers tool rotation, vibration, modulation disturbance and high-precision mesh, and optimizes efficiency based on the defined ROI region to avoid full mesh traversal. The whole process simulates Boolean cutting: the tool envelope sweeps the workpiece, leaving the minimum height. First, pre-calculation of vectors is performed, where the tool amplitude kinematics formulas for the X and Z directions are as follows: Meanwhile, to ensure the tool simulates a cocycloidal path for sweeping cutting and feed, the kinematic formula is as follows: The vibration displacements Vx and Vz are pre-calculated using the above kinematic formulas, and the tool motion trajectory, i.e., the modulation displacement X, is also calculated. mod and Y mod (X cosine forms cocycloidal perturbation, Y linear feed + sine perturbation), Z baseline Z base Maximum offset O max ROI search radius R search (Contact projection + 1.5 margin + step offset to ensure complete coverage), while using an angle vector with initial phase displacement to simulate the cutting edge rotation phase. The pre-calculated vectors complete the preparation work, which not only defines the tool movement trajectory but also improves the simulation speed. Then, the outer loop iterates through the passes, setting the X baseline. base To ensure the cumulative cutting across passes is calculated, the inner loop iterates through time steps (k), continuously updating the iteration count to simulate planar milling across passes. First, the current tool center Xc (pass offset X) is calculated. base +X modulation +X amplitude), Yc (Y modulation), Zc (baseline Z) base (+Z vibration), the formula for the tool coordinate center is as follows: By specifying and adjusting the tool's coordinate position on the workpiece in real time using the tool coordinate center formula described above, we can ensure that the toolpath accurately and completely completes the cutting process, guaranteeing that the planned path meets the expected requirements and improving machining accuracy. The next step is ROI optimization, which is based on the tool center coordinates and the search radius R. search Define the tool boundary range, determine the left, right, top, and bottom rectangular index boundaries, and then judge the grid points in the potential cutting area. If a point is close to the tool or within the rectangular boundary, it will be cut; otherwise, it will be skipped. This can significantly improve the program's running speed and avoid invalid cuts. The index formulas for the X and Y directions (points within the search range) are given below. The formula for calculating the minimum X index is shown below: X min It is the minimum mesh index calculation for the ROI in the X direction (counting from 1), which is the current tool center minus R. search Divide by the grid precision (grid_size, denoted by g here), mapping the continuous space to a discrete grid, defining the left boundary of the ROI, and then determining the formula for calculating the maximum X index as follows: X max This is the calculation of the maximum grid index of the ROI in the X direction (Nx is the total number of grid columns in the X direction), which adds R to the current tool center. search Dividing by the grid precision (grid_size, denoted by g here), we map the continuous space to a discrete grid, define the right boundary of the ROI, and ensure that the index covers the entire search radius and is within the search range. Below is the formula for calculating the minimum Y index direction: Y min It is the minimum mesh index calculation for the ROI in the Y direction (counting from 1), which is the current tool center minus R. search Divide by the grid precision (grid_size, denoted by g here), mapping the continuous space to a discrete grid, defining the lower boundary of the ROI, and then determining the formula for calculating the maximum Y index as follows: Y max This is the calculation of the maximum mesh index of the ROI in the Y direction (Ny is the total number of mesh rows in the Y direction), which is the current tool center minus R. search Divide by the grid precision (grid_size, denoted by g here), mapping the continuous space to a discrete grid, defining the upper boundary of the ROI, and forming the ROI range (rectangular region) using the above formula. Then, search for potential cutting points within the ROI range. The next step is to perform angle sweep determination. The core idea and process is to simulate the tool cutting only within a finite arc segment during rotation. Check whether a point on the workpiece is located within the sector angle covered by the tool's current rotation and whether the distance is within the tool radius. If so, the surface height value is updated; otherwise, it is skipped. The process involves calculating the azimuth angle of the grid point relative to the tool center, as shown in the following formula: Next, the current blade angle is obtained, and the formula is as follows: Next, calculate the minimum angle between the current cutting edge angle and the azimuth angle of the cutting center (to ensure the minimum path), and then calculate the angle through which the tool rotates (sweep step length). The formula is as follows: Then, it is determined whether the grid point is within the cutting region (based on both angle and distance conditions). If it is, it is cut, and the surface height value is updated. The formula is as follows: Here Z surf The value is based on the real-time surface update of the Z-map algorithm, which is the principle of material removal. Its logic is that the initial workpiece height is defined as 0, and the height value is negative after the tool enters the workpiece. If the current tool height value Z... tool The height Z of the workpiece surface surf If it's even lower (deeper), it means the tool has cut to this point, meaning the Z-axis needs to be updated. surf If the tool height is higher than the surface height (indicating the tool is moving idly, or the cut has been made deeper), then Z is not updated. surf ; Step S5, finally, involves post-processing of the data: section extraction, roughness calculation, and result visualization. First, the X and Y section height data are extracted, and the formula is shown below: Next, the roughness values Ra, Rq, and Rt are calculated. The formula for the arithmetic mean deviation Ra is as follows: Where N represents the sampling points, the next step is to calculate Rq. The formula for the root mean square deviation Rq is as follows: Finally, the total height Rt of the outline is calculated, using the formula shown below: After completing the above calculations, the values are converted to micrometers, and the surface roughness is calculated to measure and objectively reflect the surface quality of the entire plane. The formulas for the surface arithmetic mean deviation Sa and the surface root mean square deviation Sq are as follows: Finally, the processed and calculated data are visualized, including the overall 3D surface morphology, cross-sectional profiles (X-direction and Y-direction cross-sections), and a summary of roughness parameters (Ra, Rt, Rq, Sa, Sq).
2. The simulation method for ball end milling coupled with ultrasonic vibration and cocycloidal modulation according to claim 1, characterized in that: In step S1, regarding the workpiece setting, the lengths of the workpiece lx and ly are set to 1×1mm respectively. To ensure computational efficiency and balance the accuracy requirements of microscopic details, the grid precision grid_size is set to 0.
003. Simultaneously, the cocycloidal modulation parameter A is determined. mod The value is 0.029, the modulation frequency f_mod is 28, the number of cutting edges is single-edged by default, and the tool radius is 2mm.
3. The simulation method for ball end milling coupled with ultrasonic vibration and cocycloidal modulation according to claim 1, characterized in that: In step S2, the kinematic calculation preprocessing part converts the spindle speed n into angular velocity w, and the formula is as follows: As shown: Since the spindle speed is measured in r / min, dividing by 60 converts it to rad / s. Additionally, the total simulation time T_total is determined by the feed rate Vf and the workpiece length ly, as shown in the following formula: Here, the unit is also unified as seconds (s). Then, the calculation of the time step dt is shown in the following formula: dt corresponds to the time step (s) of 1° rotation. The rotational speed per second is n / 60, and the degree per second is (rotational speed) × 360°. Therefore, the time for 1° is 1 / (degree / s). Finally, the cutting pass is defined, and its formula is as follows: The above formula is used to discretize the continuous motion of the tool into a gridded point set for computation, in order to capture the microscopic sweeping details of the cutting edge.
4. The simulation method for ball end milling coupled with ultrasonic vibration and cocycloidal modulation according to claim 1, characterized in that: In step S3, the surface machining accuracy of the workpiece is determined by both the mesh accuracy and the workpiece size, as shown in the following formula: In this context, grid precision is represented by g, grid_size is represented by lx, and workpiece width is represented by ly.
5. The simulation method for ball end milling coupled with ultrasonic vibration and cocycloidal modulation according to claim 1, characterized in that: In step S4, regarding the surface quality characteristic parameters Sa and Sq, MN represents the surface matrix with M×N sampling points.