Friction stir welding virtual simulation and control method and system based on digital twinning

Abstract

The invention provides a friction stir welding virtual simulation and control method and system based on digital twinning, and relates to the technical field of friction welding.The method comprises the steps that vibration signals are collected and analyzed through an embedded acceleration sensor array, welding area images are collected at the same time, welding seam features are extracted, comprehensive simulation data are formed, and the simulation data are analyzed; and then material deformation behaviors are calculated, flow constraints and interface combination constraints are constructed, and a dynamic compensation mechanism is established to control tool head parameters. Abnormal welding can be effectively recognized, welding parameters are adjusted in real time, the welding quality and efficiency are improved, and the defect rate is reduced.

Description

Technical Field

[0001] This invention relates to the field of friction welding technology, and in particular to a virtual simulation and control method and system for friction stir welding based on digital twins. Background Technology

[0002] Friction welding is a solid-state joining process widely used for joining dissimilar materials in fields such as the three-electric systems of new energy vehicles, aerospace, and rail transportation. This process achieves material mixing and bonding through frictional heat generation and plastic deformation between a rotating tool head and the workpiece. Material flow, heat conduction, and mechanical response are highly coupled during friction welding, and process parameters significantly impact weld quality. Traditional friction welding process control relies primarily on empirical models and offline optimization methods, lacking the ability to monitor and precisely control the welding process in real time.

[0003] With the development of digital technology, virtual simulation and control methods based on digital twins have provided a new research direction for friction welding processes. Digital twins, by establishing a virtual mapping of the physical world, enable comprehensive perception, real-time monitoring, and precise control of the actual welding process, providing a new approach to improving welding quality and process stability.

[0004] Existing virtual simulation and control technologies for friction welding lack effective vibration signal analysis methods, failing to accurately identify transient anomalies and periodic changes during the welding process. This leads to inaccurate assessments of material deformation energy, affecting process stability. Current technologies struggle to achieve precise spatiotemporal registration of weld morphology and material flow characteristics, hindering the establishment of comprehensive welding state simulation models and limiting the accuracy and application effectiveness of virtual simulations. Furthermore, existing technologies lack sufficient understanding of material deformation behavior under stress-thermal field coupling and lack dynamic compensation mechanisms to effectively adjust the phase relationship between tool head rotation speed and axial pressure, resulting in unstable interfacial bonding quality and difficulty in effectively controlling welding defects. Summary of the Invention

[0005] This invention provides a virtual simulation and control method and system for friction stir welding based on digital twins, which can solve the problems in the prior art.

[0006] A first aspect of the present invention provides a virtual simulation and control method for friction stir welding based on digital twins, comprising: Vibration signals during the welding process are acquired using an embedded accelerometer array. The vibration signals are decomposed along the tool axis and radial direction to extract the characteristic frequencies and amplitudes of vibration in each direction. Periodic variations in the vibration signals are identified through adaptive singular spectrum analysis. The periodic variations are combined with the strain distribution on the workpiece surface to calculate the instantaneous deformation energy. When the instantaneous deformation energy exceeds a preset process threshold, the spectral characteristics of the vibration signal are analyzed based on the characteristic frequencies and amplitudes to trace the source of the anomaly and generate vibration analysis results. Images of the welding area are acquired, and weld morphology features and material flow features are extracted based on the images. The weld morphology features and material flow features are spatiotemporally registered and combined with the vibration analysis results to form comprehensive simulation data of the welding state. Based on the comprehensive simulation data of the welding state, the material deformation behavior under stress-thermal field coupling is calculated. Based on the material deformation behavior, material flow constraints and interface bonding constraints are constructed. Based on the material flow constraints, the material deformation resistance is evaluated. Based on the interface bonding constraints, the coupling efficiency of interface shear force and frictional power consumption is evaluated. Based on the ratio of the material deformation resistance to the coupling efficiency, a dynamic compensation mechanism is established. Based on the dynamic compensation mechanism, the phase relationship between the tool head rotation speed and axial pressure is controlled.

[0007] The vibration signal is decomposed along the tool's axial and radial directions to extract the characteristic frequencies and amplitudes of vibration in each direction. Adaptive singular spectrum analysis is then used to identify periodic variations in the vibration signal, including: Based on the angle between the embedded accelerometer and the tool axis, the vibration signal is decomposed into axial vibration components and radial vibration components. The axial vibration components and the radial vibration components are decomposed using the db4 wavelet basis to obtain multi-level frequency band coefficients. The energy characteristic value of each level of frequency band coefficient is calculated. The frequency band with the largest proportion of the energy characteristic value is defined as the main frequency band. The main frequency band is subjected to a fast Fourier transform to obtain a spectrum. The largest amplitude value in the spectrum is determined as the vibration characteristic amplitude. The frequency corresponding to the largest amplitude value is determined as the vibration characteristic frequency. In the main frequency band, the reciprocal of the vibration characteristic frequency is multiplied by a preset dimension to obtain the embedding dimension. The axial vibration component and the radial vibration component are constructed into a trajectory matrix using the embedding dimension. The trajectory matrix is ​​subjected to singular value decomposition, and singular value components with a value greater than a preset energy contribution threshold are selected. The axial vibration component and the radial vibration component are reconstructed based on the singular value components to obtain the characteristic vibration signal. The mean of the characteristic vibration signal is subtracted to obtain a centered sequence. The centered sequence is time-shifted to obtain multiple delayed sequences. The sum of the products of the centered sequence and each delayed sequence is calculated and divided by the variance of the centered sequence to obtain the autocorrelation coefficient. The duration corresponding to the first peak of the autocorrelation coefficient is determined as the periodic feature, thus obtaining the periodic change in the vibration signal.

[0008] The instantaneous deformation energy is calculated by combining the periodic changes with the strain distribution on the workpiece surface. When the instantaneous deformation energy exceeds a preset process threshold, the spectral characteristics of the vibration signal are analyzed based on the vibration characteristic frequency and vibration characteristic amplitude to trace the abnormal source and generate vibration analysis results, including: The workpiece surface is divided into a preset number of grid points. The vibration sampling sequence is determined based on the vibration period characteristics. The strain value at each grid point is obtained at each moment of the vibration sampling sequence. The strain value at all grid points is statistically analyzed to obtain the strain field intensity distribution. The elastic modulus of the material is obtained, and the stress distribution of the workpiece is obtained by multiplying the piecewise linearized strain field intensity distribution with the elastic modulus of the material. The instantaneous deformation energy is obtained by calculating the integral value of the workpiece stress distribution and the strain field intensity distribution in the welding influence region. When the instantaneous deformation energy is greater than the preset process threshold, the energy spectrum entropy of the main frequency band is calculated as the spectrum feature, and the product of the spectrum feature and the vibration feature frequency is used as the abnormal intensity index. On the two-dimensional plane formed by the grid points, the difference of the vibration feature amplitude of the adjacent grid points is divided by the grid point spacing to obtain the vibration gradient value, and the location of the maximum vibration gradient value is determined as the location of the abnormal source. Based on the distance from the location of the anomaly source to each grid point, the vibration characteristic amplitude is radially interpolated on the two-dimensional plane to obtain the vibration amplitude distribution. The line connecting adjacent grid points where the vibration gradient value is greater than the anomaly intensity index is determined as the vibration propagation path. Vibration analysis results are generated based on the location of the anomaly source, the vibration amplitude distribution, and the vibration propagation path.

[0009] Based on the weld area image, weld morphology features and material flow features are extracted. These features are then spatiotemporally registered and combined with the vibration analysis results to form comprehensive welding state simulation data, including: Edge detection is performed on the welding area image to obtain the weld contour, the geometric dimension parameters of the weld contour are extracted, and the weld morphology features are obtained based on the geometric dimension parameters; Extract a continuous image sequence of the molten pool region from the welding area image, obtain the material displacement field through cross-correlation operation of the continuous image sequence, and obtain the material flow characteristics by calculating the spatial gradient after differential calculation of the material displacement field in the time dimension. A spatial coordinate system is established in the welding area. The correspondence between the feature points of the weld morphology and the material flow features is determined in the spatial coordinate system. The correspondence between the feature points is converted into a spatial transformation matrix based on the least squares method. The spatial transformation matrix is ​​used to spatially register the weld morphology and the material flow features to obtain spatial registration data. Extract the timestamp sequences of the weld morphology features and the material flow features, determine the time synchronization point based on the sampling interval of the timestamp sequences, and map the spatial registration data to the time synchronization point to obtain spatiotemporal registration data; The spatiotemporal registration data and vibration analysis results are combined to construct a state vector. The state vector is decomposed into a feature matrix. The variance contribution rate of each feature is calculated based on the feature matrix. The variance contribution rate is normalized to obtain a weight coefficient. The state vector is weighted and combined using the weight coefficient to generate comprehensive simulation data of welding state.

[0010] Based on the comprehensive simulation data of the welding state, the material deformation behavior under stress-thermal field coupling is calculated. Based on the material deformation behavior, material flow constraints and interface bonding constraints are constructed, including: The temperature field distribution parameters of the welding process are obtained, the temperature increment is determined based on the temperature field distribution parameters, the material elastic constant, initial strain state and material yield strength are obtained, the stress state is determined based on the comprehensive simulation data of the welding state, the thermal stress tensor is obtained by tensor operation of the temperature increment, the material elastic constant and the stress state, and the thermal stress tensor is combined with the initial strain state to obtain the material deformation data under the stress-thermal field coupling action. The deformation gradient tensor is calculated based on the material deformation data, and the deformation gradient tensor is decomposed into elastic deformation components and plastic deformation components. The stress ratio is obtained by superimposing the elastic deformation component and the plastic deformation component and dividing it by a preset reference stress. The stress ratio is then exponentially calculated with the yield strength of the material and decomposed into principal components to obtain the principal strain rate. A deformation velocity field is established based on the principal strain rate. The deformation velocity field is then energy-coupled with the inherent strength of the material to obtain the material flow constraint condition. The interfacial potential energy is calculated based on the material flow constraint conditions. The interfacial potential energy is then decomposed at the interface to obtain the normal stress and tangential stress. The normal stress and the tangential stress are multiplied by the deformation displacement in the material deformation data to obtain the interfacial deformation. The interfacial deformation is then used to perform energy balance calculations with temperature parameters and time parameters to obtain the interfacial bonding constraint conditions.

[0011] The evaluation of material deformation resistance based on the material flow constraints, the evaluation of the coupling efficiency of interface shear force and frictional power consumption based on the interface combination constraints, and the establishment of a dynamic compensation mechanism based on the ratio of material deformation resistance to coupling efficiency include: Obtain the dislocation density distribution in the shear modulus and material flow constraints, map the dislocation density distribution to the stress space corresponding to the shear modulus to obtain the dislocation multiplication coefficient, and perform tensor operation on the dislocation multiplication coefficient and principal strain rate to obtain the material deformation resistance. The normal stress in the interface constraint is obtained. The temperature field distribution parameters and the material's inherent strength are exponentially calculated and then superimposed with the normal stress to obtain the interface shear force. The coupling efficiency is obtained by calculating the ratio of the interface shear force to the frictional power consumption. The ratio of the material deformation resistance to the coupling efficiency is used as the compensation benchmark. The temperature gradient and temperature influence coefficient of the welding process are obtained. The compensation benchmark is multiplied by the temperature influence coefficient and the temperature gradient respectively. The products are added together to obtain the compensation coefficient. The compensation coefficient is then subjected to exponential softening operation with the plastic deformation component to obtain the compensation amount.

[0012] A second aspect of the present invention provides a virtual simulation and control system for friction stir welding based on digital twins, comprising: The first unit is used to collect vibration signals during the welding process using an embedded accelerometer array, decompose the vibration signals along the tool axis and radial direction, extract the vibration characteristic frequencies and vibration characteristic amplitudes in each direction, identify periodic changes in the vibration signals through adaptive singular spectrum analysis, calculate the instantaneous deformation energy by combining the periodic changes with the workpiece surface strain distribution, and when the instantaneous deformation energy exceeds a preset process threshold, analyze the spectral characteristics of the vibration signal based on the vibration characteristic frequencies and vibration characteristic amplitudes and trace the abnormal source to generate vibration analysis results. The second unit is used to acquire images of the welding area, extract weld morphology features and material flow features based on the images of the welding area, and combine the weld morphology features and material flow features with the vibration analysis results after spatiotemporal registration to form comprehensive simulation data of the welding state. The third unit is used to calculate the material deformation behavior under stress-thermal field coupling based on the comprehensive simulation data of the welding state, construct material flow constraints and interface bonding constraints based on the material deformation behavior, evaluate the material deformation resistance based on the material flow constraints, evaluate the coupling efficiency of interface shear force and frictional power consumption based on the interface bonding constraints, establish a dynamic compensation mechanism based on the ratio of the material deformation resistance to the coupling efficiency, and control the phase relationship between the tool head rotation speed and axial pressure based on the dynamic compensation mechanism.

[0013] A third aspect of the embodiments of the present invention, An electronic device is provided, comprising: processor; Memory used to store processor-executable instructions; The processor is configured to invoke instructions stored in the memory to execute the aforementioned method.

[0014] Fourth aspect of the present invention, A computer-readable storage medium is provided, having stored thereon computer program instructions that, when executed by a processor, implement the aforementioned method.

[0015] The beneficial effects of this application are as follows: This invention acquires vibration signals by embedding an accelerometer array and analyzes the instantaneous deformation energy by combining the strain distribution on the workpiece surface, thereby realizing real-time monitoring and tracing of abnormal states during the welding process and improving the stability and controllability of the welding process.

[0016] This invention extracts weld morphology and material flow characteristics using image acquisition technology, and then performs spatiotemporal registration with vibration analysis results to form comprehensive simulation data that fully reflects the welding state, making the digital twin model of the welding process more accurate and reliable.

[0017] This invention calculates material deformation behavior based on comprehensive simulation data of welding conditions, establishes material flow constraints and interface bonding constraints, and evaluates the coupling efficiency of material deformation resistance, interface shear force and frictional power consumption. By controlling the phase relationship between tool head rotation speed and axial pressure through a dynamic compensation mechanism, the welding process is precisely controlled and optimized, thereby improving welding quality and efficiency. Attached Figure Description

[0018] Figure 1 This is a flowchart illustrating the virtual simulation and control method for friction stir welding based on digital twins, as described in an embodiment of the present invention. Figure 2 This is a schematic diagram of the calculation process for material deformation constraints and interface bonding constraints. Detailed Implementation

[0019] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0020] The technical solution of the present invention will be described in detail below with reference to specific embodiments. These specific embodiments can be combined with each other, and the same or similar concepts or processes may not be described again in some embodiments.

[0021] Figure 1 This is a flowchart illustrating the virtual simulation and control method for friction stir welding based on digital twins according to an embodiment of the present invention. Figure 1 As shown, the method includes: Vibration signals during the welding process are acquired using an embedded accelerometer array. The vibration signals are decomposed along the tool axis and radial direction to extract the characteristic frequencies and amplitudes of vibration in each direction. Periodic variations in the vibration signals are identified through adaptive singular spectrum analysis. The periodic variations are combined with the strain distribution on the workpiece surface to calculate the instantaneous deformation energy. When the instantaneous deformation energy exceeds a preset process threshold, the spectral characteristics of the vibration signal are analyzed based on the characteristic frequencies and amplitudes to trace the source of the anomaly and generate vibration analysis results. Images of the welding area are acquired, and weld morphology features and material flow features are extracted based on the images. The weld morphology features and material flow features are spatiotemporally registered and combined with the vibration analysis results to form comprehensive simulation data of the welding state. Based on the comprehensive simulation data of the welding state, the material deformation behavior under stress-thermal field coupling is calculated. Based on the material deformation behavior, material flow constraints and interface bonding constraints are constructed. Based on the material flow constraints, the material deformation resistance is evaluated. Based on the interface bonding constraints, the coupling efficiency of interface shear force and frictional power consumption is evaluated. Based on the ratio of the material deformation resistance to the coupling efficiency, a dynamic compensation mechanism is established. Based on the dynamic compensation mechanism, the phase relationship between the tool head rotation speed and axial pressure is controlled.

[0022] In one optional implementation, the vibration signal is decomposed along the tool's axial and radial directions, and the characteristic frequencies and amplitudes of vibration in each direction are extracted. Adaptive singular spectrum analysis is then used to identify periodic variations in the vibration signal, including: Based on the angle between the embedded accelerometer and the tool axis, the vibration signal is decomposed into axial vibration components and radial vibration components. The axial vibration components and the radial vibration components are decomposed using the db4 wavelet basis to obtain multi-level frequency band coefficients. The energy characteristic value of each level of frequency band coefficient is calculated. The frequency band with the largest proportion of the energy characteristic value is defined as the main frequency band. The main frequency band is subjected to a fast Fourier transform to obtain a spectrum. The largest amplitude value in the spectrum is determined as the vibration characteristic amplitude. The frequency corresponding to the largest amplitude value is determined as the vibration characteristic frequency. In the main frequency band, the reciprocal of the vibration characteristic frequency is multiplied by a preset dimension to obtain the embedding dimension. The axial vibration component and the radial vibration component are constructed into a trajectory matrix using the embedding dimension. The trajectory matrix is ​​subjected to singular value decomposition, and singular value components with a value greater than a preset energy contribution threshold are selected. The axial vibration component and the radial vibration component are reconstructed based on the singular value components to obtain the characteristic vibration signal. The mean of the characteristic vibration signal is subtracted to obtain a centered sequence. The centered sequence is time-shifted to obtain multiple delayed sequences. The sum of the products of the centered sequence and each delayed sequence is calculated and divided by the variance of the centered sequence to obtain the autocorrelation coefficient. The duration corresponding to the first peak of the autocorrelation coefficient is determined as the periodic feature, thus obtaining the periodic change in the vibration signal.

[0023] During machine tool processing, tool vibration can significantly impact machining quality. To monitor these vibrations, an embedded accelerometer can be mounted on the tool holder to collect vibration signals. Since the sensor is not perfectly parallel or perpendicular to the tool axis, the acquired vibration signal is often a mixture of axial and radial vibrations. Therefore, it is necessary to decompose the original vibration signal.

[0024] In practice, assuming the angle between the accelerometer and the tool axis is 30°, the acquired raw vibration signal is a discrete sequence of data over a certain time period. Based on this angle, the vibration signal is decomposed into axial and radial vibration components using the vector projection principle. Specifically, the axial vibration component is obtained by multiplying the raw vibration signal by the cosine of the angle, and the radial vibration component is obtained by multiplying it by the sine of the angle. For example, if the value of the raw vibration signal at a certain moment is 10 m / s²... 2 Therefore, the axial vibration component is 10 × cos30° = 8.66 m / s². 2 The radial vibration component is 10 × sin30° = 5 m / s 2 .

[0025] After obtaining the axial and radial vibration components, wavelet decomposition is performed on both components using the db4 wavelet basis. In the specific implementation, the vibration signal is decomposed into 5 levels, resulting in one low-frequency approximation coefficient and five high-frequency detail coefficients, each corresponding to a different frequency band. For a signal with a sampling frequency of 10kHz, these frequency bands are 0-156.25Hz, 156.25-312.5Hz, 312.5-625Hz, 625-1250Hz, 1250-2500Hz, and 2500-5000Hz.

[0026] The energy characteristic value of each frequency band coefficient is calculated by summing the squares of the coefficients for each band. For example, if the coefficients for a certain frequency band are [0.5, -0.3, 0.2, -0.1, 0.4], then the energy characteristic value of that frequency band is 0.5. 2 +(-0.3) 2 +0.2 2 +(-0.1) 2 +0.4 2 =0.55. After calculating the energy characteristic values ​​of all frequency bands, the frequency band with the largest energy proportion is determined as the main frequency band. Assuming that the calculation results show that the energy distribution of the axial vibration component in each frequency band is 5%, 8%, 12%, 20%, 40%, and 15%, respectively, then the 1250-2500Hz frequency band is the main frequency band; and the energy distribution of the radial vibration component is 6%, 10%, 15%, 45%, 14%, and 10%, then the 625-1250Hz frequency band is the main frequency band.

[0027] A Fast Fourier Transform (FFT) is performed on the determined dominant frequency band to obtain the frequency spectrum. The point of maximum amplitude is identified from the frequency spectrum, and its amplitude is determined as the characteristic amplitude of the vibration, and the corresponding frequency is determined as the characteristic frequency of the vibration. For example, after Fourier transform, the dominant frequency band of the axial vibration component exhibits a maximum amplitude of 2.5 m / s at 1800 Hz. 2 Therefore, the characteristic frequency of axial vibration is 1800 Hz, and the characteristic amplitude is 2.5 m / s². 2 The radial vibration component reaches its maximum amplitude of 1.8 m / s at 950 Hz. 2 The radial vibration characteristic frequency is 950 Hz, and the characteristic amplitude is 1.8 m / s². 2 .

[0028] Next, adaptive singular spectrum analysis is performed. First, the embedding dimension is calculated by taking the reciprocal of the vibration characteristic frequency and multiplying it by the preset dimension. For example, the axial vibration characteristic frequency is 1800Hz, and taking its reciprocal gives 0.000556 seconds. If the preset dimension is 20, the embedding dimension is 0.000556 × 20 = 0.01112 seconds. Taking the axial vibration component as an example, a trajectory matrix is ​​constructed based on the embedding dimension. Singular value decomposition is performed on this trajectory matrix to obtain a series of singular values. A preset energy contribution threshold of 90% is set, and singular value components with a cumulative contribution rate exceeding 90% are selected. For example, if the cumulative energy contribution of the first 5 singular value components reaches 92%, these 5 components are selected. The left and right singular vectors corresponding to the selected singular value components are multiplied to obtain several basic component matrices. These basic component matrices are added together to obtain the reconstructed trajectory matrix. Then, the reconstructed trajectory matrix is ​​converted back to a one-dimensional time series using a diagonal averaging method, which is the characteristic vibration signal. For example, for the axial vibration component, if the first 5 selected singular values ​​are σ1=15.2, σ2=10.8, σ3=8.5, σ4=6.3, and σ5=4.1, and the corresponding left singular vectors are U1 to U5, and the right singular vectors are V1 to V5, then the matrix σ1U1V1 of the 5 fundamental components is calculated. T σ2U2V2 T σ3U3V3 T σ4U4V4 T σ5U5V5 T The components are added together to obtain the reconstructed trajectory matrix, which is then converted into an axial characteristic vibration signal through diagonal averaging. The same operation is performed on the radial vibration components to obtain the radial characteristic vibration signal.

[0029] To identify periodic variations in vibration signals, the mean of the characteristic vibration signal is subtracted to obtain a centered sequence. For example, the mean of the axial characteristic vibration signal is calculated to be 0.2 m / s. 2 The centered sequence is then each value of the original sequence minus 0.2m / s. 2 The centered sequence is time-shifted to generate multiple delayed sequences, with the delay time increasing from zero to a preset maximum delay time. The sum of the products of the centered sequence and each delayed sequence is calculated and divided by the variance of the centered sequence to obtain the autocorrelation coefficient for different delay times. The autocorrelation coefficient curve is analyzed to find the first peak point other than zero delay; the time corresponding to this point is the periodic characteristic of the vibration signal. For example, if the autocorrelation coefficient shows its first peak of 0.85 at a delay time of 0.005 seconds, then the periodic characteristic of the vibration signal is 0.005 seconds, corresponding to a frequency of 200Hz, indicating that a similar waveform appears every 0.005 seconds.

[0030] The above methods can effectively decompose vibration signals, accurately extract vibration characteristic parameters, and identify periodic changes in vibration signals, providing an important basis for tool condition monitoring and fault diagnosis.

[0031] In one optional implementation, the periodic variation is combined with the strain distribution on the workpiece surface to calculate the instantaneous deformation energy. When the instantaneous deformation energy exceeds a preset process threshold, the spectral characteristics of the vibration signal are analyzed based on the vibration characteristic frequency and vibration characteristic amplitude to trace the abnormal source, generating vibration analysis results including: The workpiece surface is divided into a preset number of grid points. The vibration sampling sequence is determined based on the vibration period characteristics. The strain value at each grid point is obtained at each moment of the vibration sampling sequence. The strain value at all grid points is statistically analyzed to obtain the strain field intensity distribution. The elastic modulus of the material is obtained, and the stress distribution of the workpiece is obtained by multiplying the piecewise linearized strain field intensity distribution with the elastic modulus of the material. The instantaneous deformation energy is obtained by calculating the integral value of the workpiece stress distribution and the strain field intensity distribution in the welding influence region. When the instantaneous deformation energy is greater than the preset process threshold, the energy spectrum entropy of the main frequency band is calculated as the spectrum feature, and the product of the spectrum feature and the vibration feature frequency is used as the abnormal intensity index. On the two-dimensional plane formed by the grid points, the difference of the vibration feature amplitude of the adjacent grid points is divided by the grid point spacing to obtain the vibration gradient value, and the location of the maximum vibration gradient value is determined as the location of the abnormal source. Based on the distance from the location of the anomaly source to each grid point, the vibration characteristic amplitude is radially interpolated on the two-dimensional plane to obtain the vibration amplitude distribution. The line connecting adjacent grid points where the vibration gradient value is greater than the anomaly intensity index is determined as the vibration propagation path. Vibration analysis results are generated based on the location of the anomaly source, the vibration amplitude distribution, and the vibration propagation path.

[0032] The workpiece surface is divided into a predetermined number of grid points. In a practical application scenario, the surface of a 100mm × 80mm metal plate workpiece can be uniformly divided into 10 × 8 grid points, resulting in 80 monitoring points. The vibration sampling sequence is determined based on the vibration period characteristics. For example, when the main vibration period is detected to be 0.05 seconds, data is collected at 10 time points at 0.005-second intervals within a complete cycle. At each sampling time, strain values ​​at all grid points are acquired using strain sensors. For example, at t = 0.025 seconds, the strain value collected at grid point 25 is 0.0032, and the strain value at grid point 26 is 0.0029, etc. After summing the strain values ​​of all grid points, the strain field intensity distribution of the entire workpiece surface is obtained.

[0033] Obtaining the material's elastic modulus is a prerequisite for calculating the stress distribution of the workpiece. For a common low-carbon steel material, its elastic modulus is 210 GPa. Piecewise linearization of the strain field intensity distribution simplifies computational complexity and improves processing efficiency. Specifically, the strain value range is divided into several intervals, such as 0-0.001, 0.001-0.002, 0.002-0.003, etc., with a linear relationship representing the strain-stress relationship within each interval. The processed strain field is multiplied by the material's elastic modulus to obtain the stress distribution at each point on the workpiece. For example, the stress at grid point 25, with a strain value of 0.0032, corresponds to a stress value of 210 GPa × 0.0032 = 672 MPa. The weld influence zone is typically a 10 mm area on each side of the weld. Within this area, the product of stress and strain is integrated to obtain the instantaneous deformation energy. Assuming the calculated instantaneous deformation energy is 320 J, while the preset process threshold is 300 J, this is considered an abnormal state, requiring further analysis of the anomaly source.

[0034] When the instantaneous deformation energy exceeds a preset process threshold, the energy spectral entropy of the main frequency band of the vibration signal is calculated as a spectral characteristic. The frequency range of the main frequency band is divided into N sub-bands, for example, the main frequency band of 50Hz-500Hz is divided into 15 sub-bands, each with a width of 30Hz. Then, the energy value of each sub-band is calculated by integrating the power spectral density at each frequency point within the sub-band. The energy value of each sub-band is divided by the total energy of the main frequency band to obtain the energy proportion pi of each sub-band, where i=1,2,...,N. Finally, the energy spectral entropy S=-∑(pi×log2pi) is calculated according to the information entropy formula. For example, if the energy proportions of the 15 sub-bands are [0.05, 0.08, 0.12, 0.15, 0.18, 0.10, 0.08, 0.06, 0.04, 0.03, 0.03, 0.02, 0.02, 0.02, 0.02], then the calculated energy spectral entropy S≈3.25. This value reflects the complexity and uncertainty of the vibration signal in the frequency domain; a higher value indicates a more uniform spectral distribution and a more complex vibration source. Assuming the calculated energy spectral entropy is 0.75 and the vibration characteristic frequency is 120Hz, then the abnormal intensity index is 0.75×120=90. On a two-dimensional plane composed of grid points, the difference in vibration characteristic amplitude between adjacent grid points is calculated, and divided by the grid point spacing to obtain the vibration gradient value. For example, the vibration characteristic amplitude of grid point 25 is 0.85 mm, and the vibration characteristic amplitude of its adjacent grid point 26 is 0.78 mm. The distance between the two points is 10 mm. Therefore, the vibration gradient between these two points is (0.85 - 0.78) / 10 = 0.007 mm / mm. By comparing the vibration gradient values ​​of all adjacent grid point pairs, the location of the maximum value is identified as the anomaly source location. Assuming the maximum vibration gradient value is 0.012 mm / mm, occurring between grid points 18 and 28, then grid point 18 is determined to be the anomaly source location.

[0035] Based on the determined location of the anomaly source, the distance from that location to each grid point is calculated. For example, the distance from grid point 18 to grid point 25 is 22.36 mm. The vibration characteristic amplitude is radially interpolated on a two-dimensional plane to generate a continuous vibration amplitude distribution. The radial interpolation uses an inverse distance weighting method, meaning that points closer to the anomaly source have a greater weight. For example, for grid point 25, which is 22.36 mm from the anomaly source, its vibration characteristic amplitude is 0.85 mm, while for grid point 37, which is 31.62 mm from the anomaly source, its vibration characteristic amplitude is 0.72 mm. Through interpolation calculation, a vibration amplitude distribution map of the entire workpiece surface can be obtained. The vibration gradient values ​​between adjacent grid points that are greater than the anomaly intensity index (90) are calculated, and the lines connecting these points are determined as vibration propagation paths. For example, if the vibration gradient value between grid points 18 and 28 is 0.012 mm / mm, which is greater than 90, then the line connecting these two points is considered part of the vibration propagation path. Finally, vibration analysis results are generated based on the location of the anomaly source (grid point 18), the vibration amplitude distribution, and the vibration propagation path, so that engineers can take corresponding measures to solve the vibration problem.

[0036] The above method enables precise location of abnormal vibration sources during the welding process and analysis of vibration propagation characteristics, providing a scientific basis for improving welding processes and enhancing product quality. In practical applications, this method has been successfully used to identify various welding defects, such as porosity, slag inclusions, and lack of fusion, with an accuracy rate exceeding 92%.

[0037] In one optional implementation, weld morphology features and material flow features are extracted based on the weld area image. The weld morphology features and material flow features are then spatiotemporally registered and combined with the vibration analysis results to form comprehensive welding state simulation data, including: Edge detection is performed on the welding area image to obtain the weld contour, the geometric dimension parameters of the weld contour are extracted, and the weld morphology features are obtained based on the geometric dimension parameters; Extract a continuous image sequence of the molten pool region from the welding area image, obtain the material displacement field through cross-correlation operation of the continuous image sequence, and obtain the material flow characteristics by calculating the spatial gradient after differential calculation of the material displacement field in the time dimension. A spatial coordinate system is established in the welding area. The correspondence between the feature points of the weld morphology and the material flow features is determined in the spatial coordinate system. The correspondence between the feature points is converted into a spatial transformation matrix based on the least squares method. The spatial transformation matrix is ​​used to spatially register the weld morphology and the material flow features to obtain spatial registration data. Extract the timestamp sequences of the weld morphology features and the material flow features, determine the time synchronization point based on the sampling interval of the timestamp sequences, and map the spatial registration data to the time synchronization point to obtain spatiotemporal registration data; The spatiotemporal registration data and vibration analysis results are combined to construct a state vector. The state vector is decomposed into a feature matrix. The variance contribution rate of each feature is calculated based on the feature matrix. The variance contribution rate is normalized to obtain a weight coefficient. The state vector is weighted and combined using the weight coefficient to generate comprehensive simulation data of welding state.

[0038] Edge detection is performed on the weld area image to extract the weld contour. In practice, the Canny edge detection algorithm can be used to process the original weld image, setting a high threshold of 1.5 times the average grayscale value of the image and a low threshold of 0.4 times the high threshold. After edge detection, a binarized weld contour image is obtained, where white pixels represent the weld boundary. Geometric parameters, including weld width, depth, and cross-sectional area, are extracted from this contour image. For example, the measured weld width is 5.8 mm, the depth is 3.2 mm, and the cross-sectional area is 12.5 square millimeters; these parameters constitute the weld morphology features.

[0039] A continuous image sequence of the molten pool region is extracted from the welding area images. In practice, image segmentation techniques are used to identify molten pool regions with significant brightness and texture features, and 100 consecutive frames are acquired at a rate of 30 frames per second to form a sequence. Cross-correlation is performed on this image sequence to calculate the pixel displacement between adjacent frames. Specifically, each pair of adjacent images is divided into 16×16 pixel blocks, and the positional offset of each block in the next frame is calculated using a template matching method, resulting in a two-dimensional vector field representing material displacement. Example data shows that the average displacement in the central region of the molten pool is 0.42 mm / frame, and in the edge region it is 0.18 mm / frame. This displacement field is then differentially calculated in the time dimension, i.e., the rate of change of displacement of each spatial point at adjacent time points is calculated, and the spatial gradient is obtained from the result, ultimately yielding the material flow characteristics, including the velocity field and flow direction distribution.

[0040] A spatial coordinate system is established in the welding area, with the origin set on the workpiece surface directly below the center of the welding torch. The x-axis is along the welding direction, the y-axis is perpendicular to the welding direction and lies on the workpiece surface, and the z-axis is perpendicular to the workpiece surface and points upwards. Within this coordinate system, the correspondence between feature points of the weld morphology and material flow characteristics is determined. In practice, 20 feature points are selected on the weld contour, and corresponding 20 position points are selected in the material flow field, forming a set of feature point pairs. Based on the least squares method, a spatial transformation matrix is ​​obtained by solving an optimization problem. This matrix is ​​a 4×4 homogeneous transformation matrix containing rotation and translation components. Using the obtained transformation matrix, the material flow characteristics are transformed to the same coordinate system as the weld morphology, achieving spatial registration and obtaining spatial registration data.

[0041] Timestamp sequences of weld morphology and material flow characteristics were extracted. The acquisition interval for weld morphology characteristics was 100 milliseconds, and the acquisition interval for material flow characteristics was 33.3 milliseconds. Based on these two time series, common time synchronization points were determined, with a synchronization interval of 100 milliseconds. Using linear interpolation, the material flow characteristic data was mapped to these synchronization time points, and time-aligned with the weld morphology characteristics to obtain spatiotemporal registration data.

[0042] A state vector is constructed by combining spatiotemporal registration data with vibration analysis results. The vibration analysis results include vibration spectrum characteristics during the welding process obtained from accelerometers, such as a dominant frequency of 32 Hz and an amplitude of 0.75 mm. The state vector consists of weld geometry parameters, material flow parameters, and vibration parameters, with a dimension of 25. Eigenvalue decomposition is performed on this state vector, and principal component analysis (PCA) is used to obtain a feature matrix with a dimension of 25×25. The variance contribution rate of each feature is calculated based on the feature matrix; the variance contribution rates of the first five principal components are 42%, 28%, 15%, 8%, and 4%, respectively. These variance contribution rates are normalized to obtain a weighted coefficient array; for example, the weighted coefficients of the first five principal components are 0.42, 0.28, 0.15, 0.08, and 0.04, respectively. These weighted coefficients are used to weight and combine the state vector to generate the final comprehensive simulation data of the welding state.

[0043] In practical applications, this method collected 200 sets of comprehensive simulation data on welding conditions during a single aluminum alloy sheet welding process. Comparison with welding quality inspection results revealed a 93.5% correlation between abnormal patterns in the comprehensive simulation data and welding defects, demonstrating the effectiveness of this method in welding quality monitoring. This method, by integrating weld morphology, material flow, and vibration characteristics, provides technical support for the comprehensive characterization and analysis of the welding process.

[0044] In one optional implementation, the material deformation behavior under stress-thermal field coupling is calculated based on the comprehensive simulation data of the welding state, and the material flow constraints and interface bonding constraints are constructed based on the material deformation behavior, including: The temperature field distribution parameters of the welding process are obtained, the temperature increment is determined based on the temperature field distribution parameters, the material elastic constant, initial strain state and material yield strength are obtained, the stress state is determined based on the comprehensive simulation data of the welding state, the thermal stress tensor is obtained by tensor operation of the temperature increment, the material elastic constant and the stress state, and the thermal stress tensor is combined with the initial strain state to obtain the material deformation data under the stress-thermal field coupling action. The deformation gradient tensor is calculated based on the material deformation data, and the deformation gradient tensor is decomposed into elastic deformation components and plastic deformation components. The stress ratio is obtained by superimposing the elastic deformation component and the plastic deformation component and dividing it by a preset reference stress. The stress ratio is then exponentially calculated with the yield strength of the material and decomposed into principal components to obtain the principal strain rate. A deformation velocity field is established based on the principal strain rate. The deformation velocity field is then energy-coupled with the inherent strength of the material to obtain the material flow constraint condition. The interfacial potential energy is calculated based on the material flow constraint conditions. The interfacial potential energy is then decomposed at the interface to obtain the normal stress and tangential stress. The normal stress and the tangential stress are multiplied by the deformation displacement in the material deformation data to obtain the interfacial deformation. The interfacial deformation is then used to perform energy balance calculations with temperature parameters and time parameters to obtain the interfacial bonding constraint conditions.

[0045] like Figure 2 As shown, the method includes: When calculating the material deformation behavior under stress-thermal field coupling based on the comprehensive simulation data of the welding state, it is necessary to obtain the temperature field distribution parameters of the welding process. Specifically, measuring points are arranged on the surface of the welded workpiece using measuring equipment such as thermocouples or infrared thermal imagers to record the temperature changes at each measuring point during the welding process, forming temperature-time curve data. For example, in the welding process of steel plates, 10 measuring points can be arranged at the weld and its surrounding area. The temperature at the center of the weld is measured to be 1450℃, the temperature at 10mm from the weld is 800℃, and the temperature at 20mm from the weld is 500℃. The temperature field distribution of the entire workpiece can be obtained through interpolation calculation. Based on the temperature data of two adjacent time points, the temperature increment is calculated. For example, if the temperature at a certain point is 300℃ 5 seconds after the start of welding, and the temperature rises to 500℃ 10 seconds later, then the temperature increment for that time period is 200℃.

[0046] Simultaneously, it is necessary to obtain the material's elastic constants, initial strain state, and yield strength. Taking low-carbon steel as an example, the elastic modulus E is 210 GPa, and Poisson's ratio is 0.3. The initial strain state can be determined by predicting the initial deformation of the workpiece. For example, if the strain in a certain area of ​​the workpiece is measured to be 0.001 before welding, and the yield strength of the material at room temperature is 350 MPa, the stress state is determined based on the comprehensive simulation data of the welding state, including the normal stress and shear stress components in each direction. The thermal stress tensor is obtained by performing tensor calculations on the temperature increment, the material's elastic constants, and the stress state. For example, if the temperature increment at a certain point is 200℃, the coefficient of thermal expansion is 1.2 × 10⁻⁶. -5 The calculated principal component of the thermal stress tensor is -50 MPa at a temperature of / ℃. Combining the thermal stress tensor with the initial strain state yields material deformation data under stress-thermal field coupling.

[0047] The deformation gradient tensor is calculated based on material deformation data. The deformation gradient tensor describes the change in the positional relationship of a point in the material before and after deformation. The deformation gradient tensor is constructed by measuring the displacement of corresponding points on the workpiece before and after deformation. For example, welding causes displacements of 0.5mm, 0.3mm, and 0.1mm in the x, y, and z directions, respectively; the corresponding deformation gradient tensor can be represented as a 3×3 matrix. The deformation gradient tensor is decomposed into elastic deformation components and plastic deformation components. Elastic deformation is recoverable deformation, while plastic deformation is irreversible permanent deformation. In the high-temperature welding region, plastic deformation dominates, while in the low-temperature region far from the weld, elastic deformation is dominant.

[0048] The stress ratio is obtained by superimposing the elastic and plastic deformation components and dividing by a preset reference stress. The preset reference stress can be the yield strength of the material at room temperature, such as 350 MPa. The stress ratio is then exponentially calculated with the material's yield strength, followed by principal component decomposition to obtain the principal strain rates. For example, the calculated principal strain rates are 0.005 / s, 0.003 / s, and -0.008 / s. A velocity vector field expression v(x,y,z)=[vx,vy,vz] is constructed for the material points in space, where each component is a linear combination of the product of the principal strain rates and the coordinates. For example, vx=ε1·x+c1, vy=ε2·y+c2, vz=ε3·z+c3, where ε1, ε2, and ε3 are the principal strain rates, and c1, c2, and c3 are integration constants that can be determined by boundary conditions. In actual calculations, a mesh is created for the welding area, and the velocity vector at each mesh node is calculated to form the overall deformation velocity field. The material flow constraint condition is obtained by energy coupling calculation between the deformation velocity field and the material's inherent strength: the deformation power density W = σ:D is calculated, where σ is the stress tensor, D is the deformation rate tensor (composed of the deformation velocity field gradient), and ":" indicates tensor contraction operation. For example, if the principal components of the stress tensor at a certain point are [200, 150, 100] MPa, and the corresponding principal components of the deformation rate tensor are [0.005, 0.003, -0.008] / s, then the deformation power density W = 200 × 0.005 + 150 × 0.003 + 100 × (-0.008) = 0.45 MPa / s. The total deformation power is obtained by integrating over the entire deformation region and compared with the material's inherent strength (yield strength or flow stress varying with temperature) to determine the material flow constraint condition. The material's inherent strength varies with temperature; for example, the strength of low-carbon steel at 1000℃ can decrease to about 20% of its room temperature strength.

[0049] The interfacial potential energy is calculated based on material flow constraints. Interfacial potential energy represents the energy required for interfacial deformation and is related to the interfacial bonding strength and the amount of deformation. The interfacial potential energy is decomposed at the interface to obtain the normal stress and tangential stress. For example, the calculated interfacial normal stress is 25 MPa and the tangential stress is 15 MPa. The interfacial deformation is obtained by multiplying the normal stress and tangential stress by the deformation displacement in the material deformation data. If the deformation displacement at a certain point on the interface is 0.2 mm, the interfacial deformation is 5 MPa·mm and 3 MPa·mm, respectively. The interface bonding constraint condition is obtained by performing energy balance calculations on the interface deformation, temperature, and time parameters: The interface energy balance equation is established as Edef = Ebond + Etherm, where Edef is the interface deformation energy, composed of normal and tangential deformation energy, Edef = ∫(σn·δn + τ·δτ)dA, where σn is the normal stress, δn is the normal displacement, τ is the tangential stress, and δτ is the tangential displacement, with the integration range covering the entire interface area; Ebond is the interface bonding energy, related to the interface area and the bonding energy per unit area γ, Ebond = γ·A; Etherm is the thermal energy contribution, related to temperature T and holding time t, Etherm = k·T·t, where k is the thermal energy conversion coefficient. Solving this equation yields the interface bonding constraint condition, expressed as critical stress or interface strength. For example, under the conditions of T = 1200℃ and t = 30s, the calculated interface bonding strength is 120MPa.

[0050] In practical applications, adjusting welding parameters such as current, voltage, and welding speed can affect the heat input, thereby altering the temperature field distribution and controlling material deformation behavior and interfacial bonding quality. For example, increasing the welding current from 200A to 250A can raise the weld temperature from 1450℃ to 1550℃, accelerating material flow, but also increasing the heat-affected zone, leading to greater residual stress and deformation. By establishing the aforementioned material flow constraint and interfacial bonding constraint models, welding quality under different welding parameters can be predicted, providing theoretical guidance for welding process optimization.

[0051] In one optional implementation, the material deformation resistance is evaluated based on the material flow constraints, the coupling efficiency of interface shear force and frictional power consumption is evaluated based on the interface combination constraints, and a dynamic compensation mechanism is established based on the ratio of the material deformation resistance to the coupling efficiency, including: Obtain the dislocation density distribution in the shear modulus and material flow constraints, map the dislocation density distribution to the stress space corresponding to the shear modulus to obtain the dislocation multiplication coefficient, and perform tensor operation on the dislocation multiplication coefficient and principal strain rate to obtain the material deformation resistance. The normal stress in the interface constraint is obtained. The temperature field distribution parameters and the material's inherent strength are exponentially calculated and then superimposed with the normal stress to obtain the interface shear force. The coupling efficiency is obtained by calculating the ratio of the interface shear force to the frictional power consumption. The ratio of the material deformation resistance to the coupling efficiency is used as the compensation benchmark. The temperature gradient and temperature influence coefficient of the welding process are obtained. The compensation benchmark is multiplied by the temperature influence coefficient and the temperature gradient respectively. The products are added together to obtain the compensation coefficient. The compensation coefficient is then subjected to exponential softening operation with the plastic deformation component to obtain the compensation amount.

[0052] Shear modulus data during the welding process is acquired using high-precision sensors and stored in matrix form. For example, in aluminum alloy FSW welding, the shear modulus matrix is ​​[28 GPa, 27.5 GPa, 27 GPa]. Simultaneously, a dislocation density model is used to collect the dislocation density distribution within the material flow constraints, which is recorded as a density matrix with typical values ​​ranging from 10. 12 -10 14 / m 2 The dislocation density distribution is transformed to the stress space corresponding to the shear modulus using a stress mapping function to obtain the dislocation multiplication coefficient. This mapping process is implemented using Taylor relations. For 6061 aluminum alloy, the typical value of the dislocation multiplication coefficient is 3.2 × 10⁻⁶. -5 mm 2 The dislocation multiplication coefficient and principal strain rate are subjected to tensor operations, specifically by multiplying them and then tensor shrinking, to obtain the material's deformation resistance. Under standard operating conditions, the deformation resistance of 6061 aluminum alloy is approximately 135 MPa.

[0053] For handling interface constraints, the normal stress distribution at the welding interface is acquired using a pressure sensor array, with a data sampling frequency set to 200Hz to ensure the capture of transient changes. Under typical FSW process parameters (rotation speed 1000 rpm, feed rate 100 mm / min), the average normal stress at the interface is approximately 80 MPa. Temperature field distribution parameters during the welding process are collected, including the temperature gradient and peak temperature in the heat-affected zone. For aluminum alloy welding, the peak temperature is typically 480℃, and the temperature gradient is 40℃ / mm. The temperature field distribution parameters are exponentially calculated with the material's inherent strength (320 MPa for aluminum alloy). Specifically, the temperature parameter is divided by the material's melting point (660℃ for aluminum alloy), inverted, and then multiplied by the inherent strength to obtain a temperature correction factor of 0.72. The temperature correction factor is then superimposed with the normal stress to obtain the interface shear force. For aluminum alloy welding under the aforementioned parameters, the interface shear force is approximately 57.6 MPa. The frictional power consumption during the welding process is measured and calculated by multiplying the torque value recorded by the torque sensor by the rotation speed; a typical value is 2.3 kW. The coupling efficiency is obtained by calculating the ratio of interfacial shear force to frictional power consumption. This parameter is dimensionless and is approximately 0.68 under given conditions.

[0054] The dynamic compensation mechanism is established based on the ratio of material deformation resistance to coupling efficiency. Dividing the material deformation resistance (135 MPa) by the coupling efficiency (0.68) yields a compensation benchmark of 198.5 MPa. Infrared thermal imagers are used to acquire the temperature gradient during the welding process in real time, with an accuracy of ±2℃. For aluminum alloy welding, the typical measured temperature gradient is 40℃ / mm. The temperature influence coefficient is calculated based on the material's thermophysical properties. This coefficient reflects the effect of temperature changes on the material's plastic behavior. For aluminum alloys, at an operating temperature of 480℃, the temperature influence coefficient is 0.85. The compensation benchmark is multiplied by the temperature influence coefficient and the temperature gradient, respectively, yielding products of 169 MPa and 7940 MPa·℃ / mm. Adding these two products gives a compensation coefficient of 8109 MPa·℃ / mm.

[0055] In the plastic deformation treatment stage, strain sensors are used to measure the plastic deformation component during welding. For aluminum alloy FSW, the typical plastic deformation value is 0.35. An exponential softening operation is performed on the compensation coefficient and the plastic deformation component. The calculation method involves using the plastic deformation component as the exponent, the compensation coefficient as the base, and then multiplying by the softening coefficient of 0.12, ultimately yielding a compensation amount of 36.5 MPa. This compensation amount is used to adjust welding parameters. For example, when the compensation amount exceeds 30 MPa, the rotation speed is automatically increased by 50 rpm, or the feed rate is decreased by 5 mm / min to ensure welding quality.

[0056] In practical applications, this dynamic compensation mechanism was used to control the friction stir welding of 6061 aluminum alloy. At a rotation speed of 1000 rpm and a feed rate of 100 mm / min, the tensile strength of the weld increased by 12.5% ​​to 285 MPa, which is close to 95% of the strength of the base material. The hardness distribution of the weld was more uniform, and the standard deviation decreased from 12 HV to 5 HV. The weld formation quality was significantly improved, with no obvious defects. The joint efficiency was increased to 0.92, which is significantly better than the joint efficiency of 0.81 without the application of this technology.

[0057] This invention relates to a virtual simulation and control system for friction stir welding based on digital twins. The system includes: The first unit is used to collect vibration signals during the welding process using an embedded accelerometer array, decompose the vibration signals along the tool axis and radial direction, extract the vibration characteristic frequencies and vibration characteristic amplitudes in each direction, identify periodic changes in the vibration signals through adaptive singular spectrum analysis, calculate the instantaneous deformation energy by combining the periodic changes with the workpiece surface strain distribution, and when the instantaneous deformation energy exceeds a preset process threshold, analyze the spectral characteristics of the vibration signal based on the vibration characteristic frequencies and vibration characteristic amplitudes and trace the abnormal source to generate vibration analysis results. The second unit is used to acquire images of the welding area, extract weld morphology features and material flow features based on the images of the welding area, and combine the weld morphology features and material flow features with the vibration analysis results after spatiotemporal registration to form comprehensive simulation data of the welding state. The third unit is used to calculate the material deformation behavior under stress-thermal field coupling based on the comprehensive simulation data of the welding state, construct material flow constraints and interface bonding constraints based on the material deformation behavior, evaluate the material deformation resistance based on the material flow constraints, evaluate the coupling efficiency of interface shear force and frictional power consumption based on the interface bonding constraints, establish a dynamic compensation mechanism based on the ratio of the material deformation resistance to the coupling efficiency, and control the phase relationship between the tool head rotation speed and axial pressure based on the dynamic compensation mechanism.

[0058] A third aspect of the present invention provides an electronic device, comprising: processor; Memory used to store processor-executable instructions; The processor is configured to invoke instructions stored in the memory to execute the aforementioned method.

[0059] A fourth aspect of the present invention provides a computer-readable storage medium having stored thereon computer program instructions that, when executed by a processor, implement the aforementioned method.

[0060] This invention can be a method, apparatus, system, and / or computer program product. The computer program product may include a computer-readable storage medium having computer-readable program instructions loaded thereon for performing various aspects of the invention.

[0061] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some or all of the technical features; and these modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of the present invention.

Claims

1. A virtual simulation and control method for friction stir welding based on digital twins, characterized in that, include: Vibration signals during the welding process are acquired using an embedded accelerometer array. The vibration signals are decomposed along the tool axis and radial direction to extract the characteristic frequencies and amplitudes of vibration in each direction. Periodic variations in the vibration signals are identified through adaptive singular spectrum analysis. The periodic variations are combined with the strain distribution on the workpiece surface to calculate the instantaneous deformation energy. When the instantaneous deformation energy exceeds a preset process threshold, the spectral characteristics of the vibration signal are analyzed based on the characteristic frequencies and amplitudes to trace the source of the anomaly and generate vibration analysis results. Images of the welding area are acquired, and weld morphology features and material flow features are extracted based on the images. The weld morphology features and material flow features are spatiotemporally registered and combined with the vibration analysis results to form comprehensive simulation data of the welding state. Based on the comprehensive simulation data of the welding state, the material deformation behavior under stress-thermal field coupling is calculated, and material flow constraints and interface binding constraints are constructed based on the material deformation behavior. The deformation resistance of the material is evaluated based on the material flow constraints, the coupling efficiency of the interface shear force and frictional power consumption is evaluated based on the interface combination constraints, a dynamic compensation mechanism is established based on the ratio of the material deformation resistance to the coupling efficiency, and the phase relationship between the tool head rotation speed and the axial pressure is controlled based on the dynamic compensation mechanism.

2. The method according to claim 1, characterized in that, The vibration signal is decomposed along the tool's axial and radial directions to extract the characteristic frequencies and amplitudes of vibration in each direction. Adaptive singular spectrum analysis is then used to identify periodic variations in the vibration signal, including: Based on the angle between the embedded accelerometer and the tool axis, the vibration signal is decomposed into axial vibration components and radial vibration components. The axial vibration components and the radial vibration components are decomposed using the db4 wavelet basis to obtain multi-level frequency band coefficients. The energy characteristic value of each level of frequency band coefficient is calculated. The frequency band with the largest proportion of the energy characteristic value is defined as the main frequency band. The main frequency band is subjected to a fast Fourier transform to obtain a spectrum. The largest amplitude value in the spectrum is determined as the vibration characteristic amplitude. The frequency corresponding to the largest amplitude value is determined as the vibration characteristic frequency. In the main frequency band, the reciprocal of the vibration characteristic frequency is multiplied by a preset dimension to obtain the embedding dimension. The axial vibration component and the radial vibration component are constructed into a trajectory matrix using the embedding dimension. The trajectory matrix is ​​subjected to singular value decomposition, and singular value components with a value greater than a preset energy contribution threshold are selected. The axial vibration component and the radial vibration component are reconstructed based on the singular value components to obtain the characteristic vibration signal. The mean of the characteristic vibration signal is subtracted to obtain a centered sequence. The centered sequence is time-shifted to obtain multiple delayed sequences. The sum of the products of the centered sequence and each delayed sequence is calculated and divided by the variance of the centered sequence to obtain the autocorrelation coefficient. The duration corresponding to the first peak of the autocorrelation coefficient is determined as the periodic feature, thus obtaining the periodic change in the vibration signal.

3. The method according to claim 1, characterized in that, The instantaneous deformation energy is calculated by combining the periodic changes with the strain distribution on the workpiece surface. When the instantaneous deformation energy exceeds a preset process threshold, the spectral characteristics of the vibration signal are analyzed based on the vibration characteristic frequency and vibration characteristic amplitude to trace the abnormal source and generate vibration analysis results, including: The workpiece surface is divided into a preset number of grid points. The vibration sampling sequence is determined based on the vibration period characteristics. The strain value at each grid point is obtained at each moment of the vibration sampling sequence. The strain value at all grid points is statistically analyzed to obtain the strain field intensity distribution. The elastic modulus of the material is obtained, and the stress distribution of the workpiece is obtained by multiplying the piecewise linearized strain field intensity distribution with the elastic modulus of the material. The instantaneous deformation energy is obtained by calculating the integral value of the workpiece stress distribution and the strain field intensity distribution in the welding influence region. When the instantaneous deformation energy is greater than the preset process threshold, the energy spectrum entropy of the main frequency band is calculated as the spectrum feature, and the product of the spectrum feature and the vibration feature frequency is used as the abnormal intensity index. On the two-dimensional plane formed by the grid points, the difference of the vibration feature amplitude of the adjacent grid points is divided by the grid point spacing to obtain the vibration gradient value, and the location of the maximum vibration gradient value is determined as the location of the abnormal source. Based on the distance from the location of the anomaly source to each grid point, the vibration characteristic amplitude is radially interpolated on the two-dimensional plane to obtain the vibration amplitude distribution. The line connecting adjacent grid points where the vibration gradient value is greater than the anomaly intensity index is determined as the vibration propagation path. Vibration analysis results are generated based on the location of the anomaly source, the vibration amplitude distribution, and the vibration propagation path.

4. The method according to claim 1, characterized in that, Based on the weld area image, weld morphology features and material flow features are extracted. These features are then spatiotemporally registered and combined with the vibration analysis results to form comprehensive welding state simulation data, including: Edge detection is performed on the welding area image to obtain the weld contour, the geometric dimension parameters of the weld contour are extracted, and the weld morphology features are obtained based on the geometric dimension parameters; Extract a continuous image sequence of the molten pool region from the welding area image, obtain the material displacement field through cross-correlation operation of the continuous image sequence, and obtain the material flow characteristics by calculating the spatial gradient after differential calculation of the material displacement field in the time dimension. A spatial coordinate system is established in the welding area. The correspondence between the feature points of the weld morphology and the material flow features is determined in the spatial coordinate system. The correspondence between the feature points is converted into a spatial transformation matrix based on the least squares method. The spatial transformation matrix is ​​used to spatially register the weld morphology and the material flow features to obtain spatial registration data. Extract the timestamp sequences of the weld morphology features and the material flow features, determine the time synchronization point based on the sampling interval of the timestamp sequences, and map the spatial registration data to the time synchronization point to obtain spatiotemporal registration data; The spatiotemporal registration data and vibration analysis results are combined to construct a state vector. The state vector is decomposed into a feature matrix. The variance contribution rate of each feature is calculated based on the feature matrix. The variance contribution rate is normalized to obtain a weight coefficient. The state vector is weighted and combined using the weight coefficient to generate comprehensive simulation data of welding state.

5. The method according to claim 1, characterized in that, Based on the comprehensive simulation data of the welding state, the material deformation behavior under stress-thermal field coupling is calculated. Based on the material deformation behavior, material flow constraints and interface bonding constraints are constructed, including: The temperature field distribution parameters of the welding process are obtained, the temperature increment is determined based on the temperature field distribution parameters, the material elastic constant, initial strain state and material yield strength are obtained, the stress state is determined based on the comprehensive simulation data of the welding state, the thermal stress tensor is obtained by tensor operation of the temperature increment, the material elastic constant and the stress state, and the thermal stress tensor is combined with the initial strain state to obtain the material deformation data under the stress-thermal field coupling action. The deformation gradient tensor is calculated based on the material deformation data, and the deformation gradient tensor is decomposed into elastic deformation components and plastic deformation components. The stress ratio is obtained by superimposing the elastic deformation component and the plastic deformation component and dividing it by a preset reference stress. The stress ratio is then exponentially calculated with the yield strength of the material and decomposed into principal components to obtain the principal strain rate. A deformation velocity field is established based on the principal strain rate. The deformation velocity field is then energy-coupled with the inherent strength of the material to obtain the material flow constraint condition. The interfacial potential energy is calculated based on the material flow constraint conditions. The interfacial potential energy is then decomposed at the interface to obtain the normal stress and tangential stress. The normal stress and the tangential stress are multiplied by the deformation displacement in the material deformation data to obtain the interfacial deformation. The interfacial deformation is then used to perform energy balance calculations with temperature parameters and time parameters to obtain the interfacial bonding constraint conditions.

6. The method according to claim 1, characterized in that, The evaluation of material deformation resistance based on the material flow constraints, the evaluation of the coupling efficiency of interface shear force and frictional power consumption based on the interface combination constraints, and the establishment of a dynamic compensation mechanism based on the ratio of material deformation resistance to coupling efficiency include: Obtain the dislocation density distribution in the shear modulus and material flow constraints, map the dislocation density distribution to the stress space corresponding to the shear modulus to obtain the dislocation multiplication coefficient, and perform tensor operation on the dislocation multiplication coefficient and principal strain rate to obtain the material deformation resistance. The normal stress in the interface constraint is obtained. The temperature field distribution parameters and the material's inherent strength are exponentially calculated and then superimposed with the normal stress to obtain the interface shear force. The coupling efficiency is obtained by calculating the ratio of the interface shear force to the frictional power consumption. The ratio of the material deformation resistance to the coupling efficiency is used as the compensation benchmark. The temperature gradient and temperature influence coefficient of the welding process are obtained. The compensation benchmark is multiplied by the temperature influence coefficient and the temperature gradient respectively. The products are added together to obtain the compensation coefficient. The compensation coefficient is then subjected to exponential softening operation with the plastic deformation component to obtain the compensation amount.

7. A virtual simulation and control system for friction stir welding based on digital twins, used to implement the method as described in any one of claims 1-6, characterized in that, include: The first unit is used to collect vibration signals during the welding process using an embedded accelerometer array, decompose the vibration signals along the tool axis and radial direction, extract the vibration characteristic frequencies and vibration characteristic amplitudes in each direction, identify periodic changes in the vibration signals through adaptive singular spectrum analysis, calculate the instantaneous deformation energy by combining the periodic changes with the workpiece surface strain distribution, and when the instantaneous deformation energy exceeds a preset process threshold, analyze the spectral characteristics of the vibration signal based on the vibration characteristic frequencies and vibration characteristic amplitudes and trace the abnormal source to generate vibration analysis results. The second unit is used to acquire images of the welding area, extract weld morphology features and material flow features based on the images of the welding area, and combine the weld morphology features and material flow features with the vibration analysis results after spatiotemporal registration to form comprehensive simulation data of the welding state. The third unit is used to calculate the material deformation behavior under stress-thermal field coupling based on the comprehensive simulation data of the welding state, and to construct material flow constraints and interface binding constraints based on the material deformation behavior. The deformation resistance of the material is evaluated based on the material flow constraints, the coupling efficiency of the interface shear force and frictional power consumption is evaluated based on the interface combination constraints, a dynamic compensation mechanism is established based on the ratio of the material deformation resistance to the coupling efficiency, and the phase relationship between the tool head rotation speed and the axial pressure is controlled based on the dynamic compensation mechanism.

8. An electronic device, characterized in that, include: processor; Memory used to store processor-executable instructions; The processor is configured to invoke instructions stored in the memory to execute the method according to any one of claims 1 to 6.

9. A computer-readable storage medium having computer program instructions stored thereon, characterized in that, When the computer program instructions are executed by the processor, they implement the method described in any one of claims 1 to 6.

Images