Milling stability analysis method based on Bernoulli distribution and hybrid drive method

Abstract

The present invention discloses a milling stability analysis method based on Bernoulli distribution and hybrid drive method. The method establishes a milling dynamics model to construct a milling stability judgment model, establishes a probability function of Bernoulli distribution to judge the stability of the milling process, constructs a correction function, establishes a data-driven model of the spectral radius, and uses a genetic algorithm to optimize other attribute parameters. Finally, the optimized attribute parameters are obtained and an SLD graph is drawn. The present invention uses the spectral radius of the state transfer matrix to judge milling stability through theoretical analysis. By analyzing the properties of the SLD graph, two different optimization methods are used to first quickly determine the natural circular frequency to draw an SLD graph with higher accuracy, and then other attribute parameters are determined by comprehensively considering accuracy and stability.

Description

Technical Field

[0001] The present invention relates to a milling stability lobe diagram optimization analysis technology, and in particular to a milling stability analysis method based on Bernoulli distribution and a hybrid drive method. Background Art

[0002] Due to the large size, high material removal rates, numerous thin-walled areas, and poor rigidity of large aerospace structural components, high-speed milling not only easily causes deformation but also vibration, making it difficult to control machining accuracy and surface quality. The dynamic variation in chip thickness caused by the regenerative effect is the most significant factor in milling chatter. Milling stability analysis using stability lobe diagrams is a widely used method to guide actual machining, addressing the inevitable chatter that occurs during milling.

[0003] Currently, related research focuses on improving the accuracy of SLD maps. Due to the inaccurate input parameters of physical model-based methods and the lack of versatility and physical interpretability of data-driven methods, few studies have combined the characteristics of these two methods to analyze milling stability. Summary of the Invention

[0004] In order to accurately perform milling stability analysis and find reasonable vibration and stability boundary conditions, the present invention combines the advantages of physical model method and data-driven method to propose a hybrid-driven analysis method to predict milling stability, and finally redraws the stability lobe diagram to find the vibration stability and stability boundary.

[0005] The present invention achieves the above-mentioned purpose by adopting the following technical solutions: A milling stability analysis method based on Bernoulli distribution and hybrid drive method, the steps of which are as follows:

[0006] Step 1: Establish a milling dynamics model to construct a milling stability judgment model, and then use the spectrum radius to solve the chatter stability.

[0007] S01: The dynamic milling process considering the regenerative effect is described by a second-order time-delay differential equation:

[0008]

[0009] Where M is the modal mass matrix, q(t) is the tool vibration displacement vector, C is the modal damping matrix, K is the modal stiffness matrix, and a p is the axial cutting depth, D(t) is the dynamic milling force matrix, and T is the cutting cycle of a cutter tooth;

[0010] S02: Use state-space formal equations to describe the dynamic milling process:

[0011]

[0012] Where, x(t)=[q(t),p(t)] T , A is the constant coefficient matrix related only to the modal parameters, and B(t) is the periodic coefficient matrix related only to the dynamic milling force;

[0013] S03: Describe the vibration equation in any time interval:

[0014]

[0015] Where, t i is a discrete time, and the vibration equation when not cutting is

[0016] S04: represents the state items x1, x2, x3, x4 at time t1, t2, t3, t4 respectively:

[0017]

[0018] Where, t f It is the time period from the initial moment to the start of cutting by the cutter teeth;

[0019] S05: Separate the state term and the time delay term:

[0020] Discrete point t i+4 The state item x at i+4 express:

[0021]

[0022] Separate the state term and the time delay term of Equation (5):

[0023]

[0024] Where:

[0025] S06: Discrete form of continuous form state equation:

[0026] The simultaneous equations (4, 6) represent the discrete form of the continuous state equation in equation (2):

[0027] GX m+1 =HX m+1-T (7)

[0028] Representation of the state transfer matrix within a milling cycle:

[0029] Ψ=G -1 H (8)

[0030] Formula (8) is the milling stability judgment model; stability is judged according to the spectral radius of the state transfer matrix Ψ;

[0031] Table 1 Relationship between the state transfer matrix spectrum radius and system state

[0032] State transfer matrix spectral radius System Status λ(Ψ)<1 The system is in a stable state λ(Ψ)>1 The system is in an unstable state λ(Ψ)=1 The system is in a critical state

[0033] Step 2: Establish a probability function of Bernoulli distribution to judge the stability of the milling process: Since the experimental parameter collection is not completely accurate in actual situations, the milling state judgment calculated by the spectrum radius in step 1 is inaccurate. Therefore, the experimental parameters are corrected by combining the experiment with the stability judgment model to improve the prediction accuracy. The specific process is as follows:

[0034] S01: Determine the milling stability results by spectrum radius:

[0035] Solving the spectrum radius of the state transfer matrix Ψ through experiments requires the milling process parameters x=(Ω,a p ) and attribute parameters θ=(ω x ,ω y ,c x ,c y ,k x ,k y ,D t ,D r ), where: Ω is the spindle speed, a p is the milling depth, ω x ,ω y is the natural circular frequency, c x ,c y is the damping ratio, k x ,k y is the modal stiffness, D t ,D r is the cutting force coefficient;

[0036] Assuming that according to formula (8), the milling stability results are as follows:

[0037]

[0038] Where y is the stable state, λ(Ψ) is the spectral radius of the state transition matrix Ψ;

[0039] S02: Convert the discrete quantity of stability results into a continuous quantity and express the milling stability in a probabilistic form:

[0040]

[0041] Where p(y|θ,x) is the result of milling stability judgment, is the chatter probability. The value of the chatter probability is a logical function constructed based on the spectrum radius property as shown in formula (11):

[0042]

[0043] Where η is a hyperparameter;

[0044] S03: Determine the prediction results based on the Bernoulli distribution law:

[0045] When λ(θ,x)=1 Taking 0.5 as the critical value, we make the judgment as shown in Table 2: Regardless of the actual situation, if the stability judgment is correct, the value of p(y|θ,x) is always greater than 0.5;

[0046] Table 2 Bernoulli distribution law

[0047]

[0048] Step 3: Construct a correction function: Conduct n subsequent flutter experiments, update the stability lobe diagram, and obtain more reasonable attribute parameters. The specific steps are as follows:

[0049] S01: Construct pre-correction function f(θ,x):

[0050] Perform n vibration experiments and construct a pre-correction function through Bernoulli distribution:

[0051]

[0052]

[0053] Where u(ζ) is the judgment function, f(θ,x) is the pre-correction function;

[0054] If the predicted result is the same as the experimental result, the correct judgment f(θ,x) increases by 1, that is, the value of the pre-correction function is the same as the correct experimental group value;

[0055] S02: Construct the re-correction function g(θ,x):

[0056] The shortcoming of the pre-correction function is that there may be chatter points corresponding to different milling conditions between pre-correction curves 1 and 2. In order to obtain a safer SLD, it is necessary to ensure that it is as stable as possible below the SLD curve from the perspective of machining accuracy and efficiency;

[0057] Milling condition distribution points The larger the value, the more stable the SLD graph is, and the correction function is constructed:

[0058]

[0059]

[0060] Where v(ζ) is the judgment function and g(θ,x) is the correction function;

[0061] S03: Construct a security correction function h(θ,x):

[0062] f(θ,x) represents the number of correct judgments, g(θ,x) represents the security of correct judgments, and constructs a security correction function:

[0063] h(θ,x)=f(θ,x)+g(θ,x) (16)

[0064] Where h(θ,x) is the security correction function;

[0065] Step 4: Establish a data-driven model for the spectral radius: Improve the calculation speed and generalize it to the range of all milling process parameters. The specific steps are as follows:

[0066] S01: Randomly sample attribute parameters within the error range to conduct simulation experiments and establish a connection between them and the spectral radius;

[0067] S02: Determine the BP neural network topology;

[0068] S03: After establishing the topological structure of the BP neural network, set relevant parameters to train the network model;

[0069] Step 5: Optimize the natural circular frequency using the golden section method:

[0070] The error range is set to the maximum for f(θ,x), and the golden section method is used to optimize and adjust the natural circular frequency to the ideal state.

[0071]

[0072] Where σ ω is the allowable error of the natural circular frequency;

[0073] S01: Set the upper and lower boundaries of the golden section ω L 、ω R ;

[0074] S02: Set the upper and lower boundaries ω′ that may be used as the new error range after each update by the golden section method L ,ω′ R ;

[0075] S03: Set the update direction of the golden section method:

[0076] The possible direction function direct(·) for each update is defined as follows: L ,(ω L +ω R ) / 2] and perform 10 equidistant samplings, and calculate an f(θ i ,x) to get Similarly, in [(ω L +ω R ) / 2,ω R ] Perform 10 equal-distance samplings to obtain Select the one with the larger value as the update direction;

[0077] S04: Set the threshold for update optimization, (ω L +ω R ) / 2 is the optimal natural circular frequency after optimization;

[0078] Step 6: Genetic algorithm optimizes other attribute parameters:

[0079] The error range is set to maximize h(θ,x), and the genetic algorithm is used to optimize and adjust the damping ratio, modal stiffness and cutting force coefficient to the ideal state;

[0080]

[0081] Where π is other attribute parameters, σ π is the allowable error of other attribute parameters;

[0082] S01: Set the parameters of the genetic algorithm:

[0083] S02: Optimization is performed using a genetic algorithm, and the fitness is defined as:

[0084]

[0085] Where, e(U) is the fitness;

[0086] Finally, the optimized attribute parameter θ is obtained and the SLD diagram is drawn.

[0087] The present invention uses the spectrum radius of the state transfer matrix to judge the milling stability through theoretical analysis. The discrete variables of the stability analysis are converted into continuous probability variables in a probabilistic form for easy optimization, and the correct judgment result is made for the actual situation through the Bernoulli distribution law. While constructing the correction function, not only the accuracy of the prediction is considered, but also the security of the prediction is guaranteed. By establishing a data-driven model of the spectrum radius, the input sample of the prediction model is expanded using error construction simulation experiments, the calculation speed is improved compared with the general physical model, and it is extended to the continuous interval range of all attribute parameters, that is, after determining the milling process parameters and conducting experiments, the corresponding SLD graph under the conditions can be quickly solved to make stability judgments. By analyzing the properties of the SLD image, two different optimization methods are used to first quickly determine the inherent circular frequency to draw an SLD graph with higher accuracy, and then determine other attribute parameters with comprehensive consideration of accuracy and stability. BRIEF DESCRIPTION OF THE DRAWINGS

[0088] Figure 1 It is a milling dynamics model diagram of the present invention; Figure 2 is a graph of the chatter probability function of the present invention;

[0089] Figure 3 is a schematic diagram of the same pre-correction function value in the present invention;

[0090] Figure 4a It is the SLD curve with a chatter probability greater than 0.5; Figure 4b It is the SLD curve with a chatter probability greater than 0.5;

[0091] Figure 5 It is a network topology diagram;

[0092] Figure 6a is the effect of natural circular frequency on SLD; Figure 6b is the effect of damping ratio on SLD;

[0093] Figure 6c is the effect of modal stiffness on SLD; Figure 6d is the effect of cutting force coefficient on SLD;

[0094] Figure 7 It is the data point for dynamic optimization;

[0095] Figure 8a It is the SLD diagram of whether there is chatter before optimization; Figure 8b is the prediction accuracy of the SLD graph before optimization;

[0096] Figure 9a It is the SLD diagram of whether there is chatter after the golden section optimization;

[0097] Figure 9b is the prediction accuracy of the SLD graph after optimization by the golden section method;

[0098] Figure 10a It is the SLD diagram of whether there is chatter after genetic algorithm optimization;

[0099] Figure 10b is the prediction accuracy of the SLD graph after genetic algorithm optimization; DETAILED DESCRIPTION

[0100] The present invention will be further described below with reference to the accompanying drawings and embodiments.

[0101] A milling stability analysis method based on Bernoulli distribution and hybrid drive method, the steps are as follows:

[0102] Step 1: Establish milling dynamics model (such as Figure 1 As shown), F is the milling force, f z is the feed rate, Ω is the spindle speed, k is the modal stiffness, and c is the damping ratio.

[0103] S01: Dynamic milling process considering regenerative effects:

[0104] It is described by the second-order delay differential equation as shown in formula (1):

[0105]

[0106] In formula (1), M is the modal mass matrix, C is the modal damping matrix, K is the modal stiffness matrix, q(t) is the tool vibration displacement vector, f(t) is the static milling force matrix, D(t) is the dynamic milling force matrix, and a p is the axial cutting depth, T=60 / NΩ is the cutting cycle of one tooth, and N is the number of teeth.

[0107] Since the static force does not affect the stability of the system, Equation (1) can be degenerated into a form that omits the static milling force matrix:

[0108]

[0109] In formula (2), the relationship between the modal mass matrix M, the modal stiffness matrix K and the natural circular frequency ω is:

[0110] |K-ω 2 M|=0 (3)

[0111] S02: Convert the second-order delay differential equation into a state-space form equation:

[0112] If the order x(t)=[q(t),p(t)] T , then formula (2) can be expressed as:

[0113]

[0114] In formula (4), is a constant coefficient matrix related only to the modal parameters, It is a periodic coefficient matrix related only to the dynamic milling force and satisfies the condition B(t)=B(t+T).

[0115] S03: Convert the continuous time equation into discrete time equation:

[0116] The moment when the cutter tooth leaves the workpiece after the previous cutting cycle is the initial moment t0 of the current cycle, and the time period from t0 to the start of cutting is t f , at t = t0 + t f The time period from when the workpiece is cut into until the end of the current cycle is Tt f At t f The blade teeth vibrate freely during the time period Tt fThe forced vibration of the cutter teeth in the time period. A cutter tooth period T can be divided into a non-cutting time period and a cutting time period. The forced vibration time period is discretized into m time periods with an interval length of τ, and the continuous time is expressed as discrete time as shown in formula (5):

[0117] t i =t0+t f +(i-1)τ (5)

[0118] Where, i=1,2,...,m,m+1.

[0119] In formula (4), p B(t)[x(t)-x(tT)] is considered as a homogeneous equation The non-homogeneous terms of are transformed into formula (6):

[0120]

[0121] If t∈[t i ,t i+1 ], the vibration equation in any time interval is described by formula (6) as formula (7):

[0122]

[0123] Since in [t0,t0+t f ] The tool inside does not cut, so B(s) = 0, and equation (6) degenerates into equation (8):

[0124]

[0125] S04: represents the state item x1 at time t1:

[0126] Note x i =x(t i ), x i-T =x(t i -T), B i =B(t i ), B i-T =B(t i -T), and express equation (8) as equation (9):

[0127]

[0128] S05: represents the state item x2 at time t2:

[0129] On the interval [t1, t2], the state term x2 at the discrete point t2 is expressed as (10) using equation (7):

[0130]

[0131] In formula (10) The trapezoidal quadrature formula is used to approximate the formula (11):

[0132]

[0133] Formula (11) is classified and organized according to the state term and time lag term as shown in Formula (12):

[0134]

[0135] S06: represents the state item x3 at time t3:

[0136] Similarly, the state term x3 at discrete point t3 is expressed as (13) from (7):

[0137]

[0138] In formula (13) The Simpson quadrature formula is approximated as formula (14):

[0139]

[0140] Similarly, formula (14) is classified and organized according to the state term and time lag term as formula (15):

[0141]

[0142] S07: represents the state item x4 at time t4:

[0143] Similarly, the state term x4 at discrete point t4 is expressed as (16) from equation (7):

[0144]

[0145] Formula (16) is further described by Newton's quadrature formula as formula (17):

[0146]

[0147] In the interval [t i ,t i+4 ], the discrete point t is i+4 The state item x at i+4 It is expressed as formula (18):

[0148]

[0149] S08: Separate the state term and the time delay term by Cotes quadrature formula:

[0150] The state term and time delay term of equation (18) are separated as follows:

[0151]

[0152] In formula (19):

[0153]

[0154]

[0155] S09: Solve the discrete form of the final continuous form of the state equation:

[0156] Simultaneously combining equations (9, 12, 15, 17, 19), the discrete form of the continuous form state equation of equation (4) is equation (22):

[0157] GX m+1 =HX m+1-T (twenty two)

[0158] In formula (22)

[0159]

[0160]

[0161]

[0162] Therefore, the state transition matrix within a milling cycle can be expressed as:

[0163] Ψ=G -1 H (26)

[0164] According to Floquet theory, the system stability can be judged based on the spectral radius of the state transfer matrix Ψ (as shown in Table 1). Equation (26) is called the judgment model of milling stability.

[0165] Table 1 Relationship between the state transfer matrix spectrum radius and system state

[0166] State transfer matrix spectral radius System Status λ(Ψ)<1 The system is in a stable state λ(Ψ)>1 The system is in an unstable state λ(Ψ)=1 The system is in a critical state

[0167] Step 2: Establish a probability function based on Bernoulli distribution. Analyze the stability of the milling process. At the same time, the milling process parameters x=(Ω,a p ) and attribute parameters θ=(ω x ,ω y ,c x ,c y ,k x ,k y ,D t ,D rSince the milling process is subject to the constraints of multiple factors and the measurement data is inaccurate, it is necessary to correct the attribute parameters. Assume that the attribute parameter error σ=(σ ωx ,σ ωy ,σ cx ,σ cy ,σ kx ,σ ky ,σ Dt ,σ Dr ). The probability function of the Bernoulli distribution can be applied to the correction of attribute parameters within the error range. The specific process is as follows:

[0168] S01: Determine the milling stability results by spectrum radius:

[0169] The milling stability results are expressed as (27) using (26):

[0170]

[0171] In Equation (27), y is the stability state, where 0 indicates stability and 1 indicates instability, and λ(Ψ) is the spectral radius of the state transition matrix Ψ.

[0172] S02: Express the milling stability results in a probabilistic form:

[0173] Continuous quantity optimization has better results in dynamic optimization problems. The Bernoulli distribution is used to convert the discrete quantity of 0 or 1 in the stability result y into a probability value. A continuous quantity. Definition is the probability of chatter, is the probability of stability, and the milling stability result in Equation (27) is expressed in probability form as Equation (28):

[0174]

[0175] S03: Construct the functional relationship between spectrum radius and chatter probability:

[0176] Flutter probability The appropriate representation method is the key factor to achieve the equivalent representation between the stability result y and the probability p. Based on the properties of the spectrum radius, the logic function is constructed as shown in Equation (29). The value of the hyperparameter η will affect the speed of the target approach (such as the attached Figure 2 As shown), here we take η = 4:

[0177]

[0178] S04: Determine the prediction results based on the Bernoulli distribution law:

[0179] When λ(θ,x)=1 Taking 0.5 as the critical value, the following judgment is made: no matter what the actual experimental results are, as long as the stability judgment is correct, the value of p(y|θ,x) is always greater than 0.5 (as shown in Table 2). Table 2 Bernoulli distribution law

[0180]

[0181] Step 3: Construct a correction function and conduct n subsequent flutter experiments to update and correct the stability lobe diagram, thereby obtaining more reasonable attribute parameters. The specific steps are as follows:

[0182] S01: Construct pre-correction function f(θ,x):

[0183] Set the step function u(ζ) as p(y i |θ,x i ) activation function, retaining the correct judgment results and outputting 1, filtering out the wrong judgment results. Perform n vibration experiments to obtain Bernoulli distribution values, process them through the activation function u(ζ) and then sum them to construct the pre-correction function f(θ,x), as shown in Equations (30, 31):

[0184]

[0185]

[0186] If there are n sets of experiments used for dynamic correction and m sets of experiments with correct stability judgments, the accuracy of the model judgment is m / n. With each correct judgment, f(θ, x) increases by 1. The value of the pre-correction function is the same as the value of the correct experimental set. The accuracy of the model judgment can also be expressed as f(θ, x) / n.

[0187] S02: Construct the re-correction function g(θ,x):

[0188] The pre-correction function has the following shortcomings (such as Figure 3 The figure contains two pre-corrected SLD curves 1 and 2. The circle indicates stability and the box indicates vibration. Both curves make correct judgments. The corresponding pre-correction function and The values ​​of are equal, making it impossible to determine the effectiveness of the correction. Between pre-correction curves 1 and 2, there may be chatter points corresponding to different milling conditions. To achieve a safer SLD, from the perspective of machining accuracy and efficiency, it is necessary to ensure that the surface below the SLD curve remains as stable as possible. Selecting pre-correction curve 1 as the final corrected SLD curve is more effective than pre-correction curve 2.

[0189] The stability state corresponding to each milling condition distribution point has been converted to a chatter probability between 0 and 1 The value of the probability > 0.5 (such as Figure 4a As shown), images with probability < 0.5 (as shown Figure 4b As shown). For the milling condition distribution points, The larger the value of , the more stable the SLD graph result. The re-correction function g(θ,x) is further constructed as follows:

[0190]

[0191]

[0192] Where v(ζ) is the ReLU function, which can only retain the specific values ​​of the correct part. i =1, Take the larger value, v(ζ)p(y i |θ,x i ) also takes a larger value, then g(θ,x) takes a larger value; if the chatter is y i =0, Take the larger value, v(ζ)p(y i |θ,x i ) takes a smaller value, 1-v(ζ)p(y i |θ,x i ) takes a larger value, then g(θ,x) still takes a larger value.

[0193] S03: Construct a security correction function h(θ,x):

[0194] From the perspective of the correction function construction process, f(θ,x) represents the number of correct judgments, g(θ,x) represents the safety of correct judgments, and the value of g(θ,x) does not affect the number of correct judgments. The two do not affect each other. A more reasonable SLD safety correction function h(θ,x) can be constructed as follows:

[0195] h(θ,x)=f(θ,x)+g(θ,x) (34)

[0196] Step 4: During the update and calculation of the correction function, the spectral radius λ needs to be recalculated for each iteration. If Equation (26) is used for each update, it will take a long time to calculate. Therefore, a data-driven model of the spectral radius is established, which has a much faster calculation speed than the physical model method. The specific steps are as follows:

[0197] S01: Simulation experiment through error:

[0198] For a given milling condition, within the error range, attribute parameters are randomly sampled. After performing a spectrum radius simulation experiment using Equation (26), a BP neural network is selected to establish the relationship between attribute parameters and spectrum radius, that is, a spectrum radius neural network prediction model under certain milling conditions is established.

[0199] S02: Determine the BP neural network topology:

[0200] The network structure consists of 1 input layer, n hidden layers and 1 output layer. The input layer has 8 neurons θ i (1≤i≤8), the output layer has only one neuron, that is, the spectral radius Λ that needs to be predicted (such as Figure 5 shown).

[0201] The number of neurons in the first hidden layer, m1, is determined according to a common empirical formula to be m1 = 17. Starting from the second hidden layer, the number of neurons is determined using the network structure growth method, that is, the number of neurons is increased sequentially until there is no significant increase in the predicted error and the actual error, and the final network topology is determined.

[0202] S03: Train the network model:

[0203] After establishing the topological structure of the BP neural network, relevant parameters are set to train the network model.

[0204] Step 5: Optimize the natural circular frequency using the golden section method:

[0205] Within the error range, adjust the natural circular frequency parameters to the ideal state to maximize f(θ,x) and correct the SLD. Since the left and right positions of the SLD graph are only related to the natural circular frequency (such as Figure 6a As shown in (35), the golden section method is used to optimize the natural circular frequency.

[0206]

[0207] S01: Set the upper and lower boundaries of the golden section:

[0208] Set the natural circular frequency ω at the initial error σ ω The upper and lower boundaries of the range ω L 、ω R .

[0209] S02: Set the update boundary of the golden section method:

[0210] Set the upper and lower boundaries ω′ that may be used as the new error range after each update by the golden section method L ,ω′ R .

[0211] S03: Set the update direction of the golden section method:

[0212] The possible direction function direct(·) for each update is defined as follows: L ,(ω L +ωR ) / 2] interval, and calculate an f(θ i ,x), we get Similarly in [(ω L +ω R ) / 2,ω R ] interval is also sampled 10 times with equal intervals, and we get Select the larger value as the update direction.

[0213] S04: Set thresholds for update optimization:

[0214] The iterative process is carried out until the current circular frequency error does not exceed the threshold ε ω This termination condition. (ω L +ω R ) / 2 is taken as the optimal value of the natural circular frequency after the golden section method optimization.

[0215] Step 6: Genetic algorithm optimizes other attribute parameters:

[0216] Within the error range, adjust the damping ratio, modal stiffness and cutting force coefficient parameters to the ideal state to maximize h(θ,x) and thus correct the SLD. Figure 6b 、 Figure 6c and Figure 6d As shown in (36), the genetic algorithm is used to optimize the damping ratio, modal stiffness and cutting force coefficient.

[0217]

[0218] S01: Set the parameters of the genetic algorithm:

[0219] Set the population size, upper and lower limits of variables, fitness function, selection function, crossover probability and mutation probability.

[0220] S02: Optimization using genetic algorithms:

[0221] The genetic algorithm is used to solve Equation (37) to find the damping ratio, modal stiffness, and cutting force coefficient that maximize the objective function h(θ,x). The fitness can be defined as:

[0222]

[0223] Example: A milling experiment was carried out on a milling system consisting of a milling machine NCT EmR-610Ms, a two-tooth end mill with a helix angle β = 30° and a diameter of 16 mm, and an aluminum alloy 2024-T351 workpiece. In order to facilitate the observation of the effect, only a single lobe diagram of the SLD diagram was optimized to verify the proposed milling stability analysis method based on Bernoulli distribution and hybrid drive method.

[0224] Step 1: Attribute parameter analysis:

[0225] The cutting force was collected by a Kistler 9129AA multi-component dynamometer and an NI-9234 acquisition card, and the modal parameters were obtained by an Endevco hammer 2302-10 and a PCB accelerometer 352C32. The experimentally measured property parameters are shown in Table 3.

[0226] Table 3 Experimental values ​​of property parameters

[0227]

[0228] S01: Use the attribute parameters to obtain the SLD map through the Cotes integration method:

[0229] Using the milling system attribute parameters shown in Table 3, the SLD diagram is obtained by the Cotes integration method of formula (26), see Figure 6a 、 Figure 6b 、 Figure 6c and Figure 6d The solid black line in .

[0230] S02: Analyze the influence of natural circular frequency on SLD diagram:

[0231] Natural circular frequency ω x ,ω y Impact on SLD graph (such as Figure 6a As shown), 90% ω in the legend x =4426rad / s, 90%ω y =4257rad / s, and so on. As can be seen from the figure, the natural circular frequency mainly changes the left and right positions of the SLD graph.

[0232] S03: Analyze the influence of other factors on the SLD graph besides the natural circular frequency:

[0233] The influence of damping ratio, modal stiffness and cutting force coefficient on SLD diagram (e.g. Figure 6b 、 Figure 6c 、 Figure 6d As shown in the figure, it can be seen that the numerical changes of other factors except the natural circular frequency mainly change the height and width of the SLD graph.

[0234] Step 2: Spectral radius calculation:

[0235] S01: Simulation experiment using attribute parameters and their errors:

[0236] Within the error range, 1000 simulation experiments were conducted on the attribute parameters using a uniformly distributed random number method.

[0237] S02: Calculate the spectrum radius corresponding to different milling conditions under each set of attribute parameters:

[0238] Figure 7 The figure in the center shows the scattered data points for dynamic optimization. Circles represent stability, crosses represent chatter, and triangles indicate indeterminate results. The grid node step size for spindle speed is set to 100, and the grid node step size for depth of cut is set to 0.125. The horizontal axis represents the spindle speed discretized into 85 nodes, and the vertical axis represents the depth of cut discretized into 27 nodes. The coordinate plane contains 2295 grid nodes, each corresponding to a milling condition. Through 1000 simulation experiments, the spectral radius corresponding to 2295 milling conditions under 1000 sets of attribute parameters was solved.

[0239] S03: Establish the prediction model topology of spectrum radius through BP neural network:

[0240] The SLD graph for each attribute parameter in the 1000 simulations contained 2295 grid nodes, and a 1000 x 2295 neural network prediction model was established. The input neurons were 8 neurons, representing the attribute parameters, and the output neurons were 1 neuron, representing the spectral radius. The network structure was set to 8-17-8-1.

[0241] S04: Set the BP neural network related parameters to train the model:

[0242] The network activation functions were set to F1 = logsig, F2 = logsig, and F3 = purelin, with a maximum iteration cycle of 100 and a learning rate of 0.01. The input parameters were preprocessed and standardized, and the training algorithm, trainlm, was selected. 950 groups randomly selected from the 1000 simulated groups were used as the training set, and the remaining 50 groups were used as the test set.

[0243] Step 3: Dynamically optimize the single leaflet graph using the golden section method:

[0244] See Figure 7 ,Here, only the stable and chattering data points are used to update and ,correct the SLD graph.,The network distribution of the milling working condition is selected to ,just cover one lobe of the SLD for 43 milling processing experiments.,The SLD graph before optimization (e.g. Figure 8a As shown), the circle represents stability and the cross represents vibration. The accuracy is shown in Figure 8bAs shown in the figure, the lower triangle represents the correct prediction, the upper triangle represents the wrong prediction, and the other figures are analogous. The accuracy of the SLD before optimization is 70%. Figure 9a As shown), set the threshold ε ω 1e -6 , its accuracy is 93% (as Figure 9b shown).

[0245] Step 4: Genetic algorithm dynamically optimizes the single leaflet graph:

[0246] Keep the optimal value of the natural circular frequency unchanged and further optimize the other six attribute parameters. The SLD diagram after genetic algorithm optimization (such as Figure 10a As shown), its accuracy (as Figure 10b (as shown in Figure 2), the population size was set to 50, the upper and lower limits of the variables were set to ±20%, and the input parameters were the damping ratio, modal stiffness, and cutting force coefficient in both directions. The optimization objective was to maximize the value of the fitness function, i.e., the correction function. The selection function was random uniform selection, with a crossover probability of 0.8 and a mutation probability of 0.2. The SLD accuracy after genetic algorithm optimization was 98%. The values ​​of all model parameters after two optimizations are shown in Table 4.

[0247] Table 4 Attribute parameters after optimizing single petal

[0248]

Claims

1. The milling stability analysis method based on Bernoulli distribution and hybrid drive method is characterized by: Here are the steps: Step 1: Establish a milling dynamics model to construct a milling stability judgment model, and then use the spectrum radius to solve the chatter stability. S01: The dynamic milling process considering the regenerative effect is described by a second-order time-delay differential equation: Where M is the modal mass matrix, q(t) is the tool vibration displacement vector, C is the modal damping matrix, K is the modal stiffness matrix, and a p is the axial cutting depth, D(t) is the dynamic milling force matrix, and T is the cutting cycle of a cutter tooth; S02: Use state-space formal equations to describe the dynamic milling process: Where, x(t)=[q(t),p(t)] T , A is the constant coefficient matrix related only to the modal parameters, and B(t) is the periodic coefficient matrix related only to the dynamic milling force; S03: Describe the vibration equation in any time interval: Where, t i is a discrete time, and the vibration equation when not cutting is S04: represents the state items x1, x2, x3, x4 at time t1, t2, t3, t4 respectively: Where, t f It is the time period from the initial moment to the start of cutting by the cutter teeth; S05: Separate the state term and the time delay term: Discrete point t i+4 The state item x at i+4 express: Separate the state term and the time delay term of Equation (5): Where: S06: Discrete form of continuous form state equation: The simultaneous equations (4, 6) represent the discrete form of the continuous state equation in equation (2): GX m+1 =HX m+1-T (7) Representation of the state transfer matrix within a milling cycle: Ψ=G -1 H (8) Formula (8) is the milling stability judgment model; stability is judged according to the spectral radius of the state transfer matrix Ψ; Table 1 Relationship between the state transfer matrix spectrum radius and system state Step 2: Establish a probability function of Bernoulli distribution to judge the stability of the milling process: Since the experimental parameter collection is not completely accurate in actual situations, the milling state judgment calculated by the spectrum radius in step 1 is not accurate. Therefore, the experimental parameters are corrected by combining the experiment with the stability judgment model to improve the prediction accuracy. The specific process is as follows: S01: Determine the milling stability results by spectrum radius: Solving the spectrum radius of the state transfer matrix Ψ through experiments requires the milling process parameters x=(Ω,a p ) and attribute parameters θ=(ω x ,ω y ,c x ,c y ,k x ,k y ,D t ,D r ), where: Ω is the spindle speed, a p is the milling depth, ω x ,ω y is the natural circular frequency, c x ,c y is the damping ratio, k x ,k y is the modal stiffness, D t ,D r is the cutting force coefficient; Assuming that according to formula (8), the milling stability results are as follows: Where y is the stable state, λ(Ψ) is the spectral radius of the state transition matrix Ψ; S02: Convert the discrete quantity of stability results into a continuous quantity and express the milling stability in a probabilistic form: Where p(y|θ,x) is the result of milling stability judgment, is the chatter probability. The value of the chatter probability is a logical function constructed based on the spectrum radius property as shown in formula (11): Where η is a hyperparameter; S03: Determine the prediction results based on the Bernoulli distribution law: When λ(θ,x)=1 Taking 0.5 as the critical value, we make the judgment as shown in Table 2: Regardless of the actual situation, if the stability judgment is correct, the value of p(y|θ,x) is always greater than 0.5; Table 2 Bernoulli distribution law Step 3: Construct a correction function; conduct n subsequent flutter experiments, update the stability lobe diagram, and obtain more reasonable attribute parameters. The specific steps are as follows: S01: Construct pre-correction function f(θ,x): Perform n vibration experiments and construct a pre-correction function through Bernoulli distribution: Where u(ζ) is the judgment function, f(θ,x) is the pre-correction function; If the predicted result is the same as the experimental result, the correct judgment f(θ,x) increases by 1, that is, the value of the pre-correction function is the same as the correct experimental group value; S02: Construct the re-correction function g(θ,x): The shortcoming of the pre-correction function is that there may be chatter points corresponding to different milling conditions between pre-correction curves 1 and 2. In order to obtain a safer SLD, it is necessary to ensure that it is as stable as possible below the SLD curve from the perspective of machining accuracy and efficiency; Milling condition distribution points The larger the value, the more stable the SLD graph is, and the correction function is constructed: Where v(ζ) is the judgment function and g(θ,x) is the correction function; S03: Construct a security correction function h(θ,x): f(θ,x) represents the number of correct judgments, g(θ,x) represents the security of correct judgments, and constructs a security correction function: h(θ,x)=f(θ,x)+g(θ,x) (16) Where h(θ,x) is the security correction function; Step 4: Establish a data-driven model for the spectral radius: Improve the calculation speed and generalize it to the range of all milling process parameters. The specific steps are as follows: S01: Randomly sample attribute parameters within the error range to conduct simulation experiments and establish a connection between them and the spectral radius; S02: Determine the BP neural network topology structure; S03: After establishing the topological structure of the BP neural network, set relevant parameters to train the network model; Step 5: Optimize the natural circular frequency using the golden section method: The error range is set to the maximum for f(θ,x), and the golden section method is used to optimize and adjust the natural circular frequency to the ideal state. Where σ ω is the allowable error of the natural circular frequency; S01: Set the upper and lower boundaries of the golden section ω L 、ω R ; S02: Set the upper and lower boundaries ω′ that may be used as the new error range after each update by the golden section method L ,ω′ R ; S03: Set the update direction of the golden section method: The possible direction function direct(·) for each update is defined as follows: L ,(ω L +ω R ) / 2] and perform 10 equidistant samplings, and calculate an f(θ i ,x) to get Similarly, in [(ω L +ω R ) / 2,ω R ] Perform 10 equal-distance samplings to obtain Select the one with the larger value as the update direction; S04: Set the threshold for update optimization, (ω L +ω R ) / 2 is the optimal natural circular frequency after optimization; Step 6: Genetic algorithm optimizes other attribute parameters: The error range is set to maximize h(θ,x), and the genetic algorithm is used to optimize and adjust the damping ratio, modal stiffness and cutting force coefficient to the ideal state; Where π is other attribute parameters, σ π is the allowable error of other attribute parameters; S01: Set the parameters of the genetic algorithm: S02: Optimization is performed using a genetic algorithm, and the fitness is defined as: Where, e(U) is the fitness; Finally, the optimized attribute parameter θ is obtained and the SLD diagram is drawn.

Images