Propulsion system vibration active compensation control method based on RBF load mapping

Abstract

The application discloses a kind of based on RBF load mapping's propulsion system vibration active compensation control method, comprising: the fluctuating pressure time domain data of unsteady flow field around propeller is collected to construct flow field fluctuating pressure time domain matrix;Based on radial basis function, construct fluid-solid coupling mapping kernel function matrix;The above-mentioned matrix is input to RBF neural network mapping model, and the wet surface node load time domain sequence of propeller is obtained;Based on Lasso regularization constraint, the key node with the greatest contribution to vibration is screened;The sparse load time domain sequence is carried out frequency domain transformation, and the spatial distribution characteristics of sparse key node set are combined to generate structure surface load spectrum matrix;Structure surface load spectrum matrix is input to anti-phase control network, and anti-phase vibration compensation driving signal is applied to piezoelectric actuator array, to realize propulsion system vibration active compensation control.The application can significantly improve the response speed, suppression precision and engineering applicability of propulsion system vibration control.

Description

Technical Field

[0001] This application relates to the field of vibration control technology for ship propulsion systems, and in particular to an active compensation control method for propulsion system vibration based on RBF load mapping. Background Technology

[0002] During the propulsion process of underwater vehicles, the thruster interacts with a complex, unsteady flow field, generating intense pulsating pressure loads. These loads are transmitted to the thruster structure and hull through the fluid-structure interaction interface, inducing high-frequency vibrations and structural acoustic radiation, severely impacting stealth and crew comfort. Traditional vibration control methods mostly rely on passive vibration isolation or active control strategies based on PID feedback. The former has limited effectiveness in suppressing mid-to-high frequency vibrations, while the latter suffers from insufficient compensation accuracy due to control lag and difficulty in adapting to the rapid time-varying characteristics of the flow field loads.

[0003] In recent years, progress has been made in fluid-structure interaction load prediction methods based on computational fluid dynamics (CFD), which can interpolate and map the fluctuating pressure at fluid nodes to the loads at the wetted surface nodes of the structure through radial basis functions (RBF). However, existing technologies mostly employ static, offline mapping methods, loading the CFD time-domain data into the structural finite element model in one go after Fourier transform for frequency domain analysis. This approach cannot adapt to the continuous changes in the flow field during actual ship navigation. At the same time, existing mapping methods typically reconstruct the loads at all nodes of the wetted surface in all dimensions, resulting in high computational redundancy, which is difficult to meet the stringent real-time requirements of active control systems. Furthermore, traditional methods lack a closed-loop mechanism that couples accurate load mapping with active actuator drive, leading to a disconnect between load identification and vibration compensation, and failing to achieve true adaptive vibration suppression. Summary of the Invention

[0004] The main objective of this application is to provide a propulsion system vibration active compensation control method based on RBF load mapping, which aims to solve the technical problems of existing propulsion system vibration control relying on lag-type PID strategies, inaccurate fluid-structure interaction load transfer and excessive computational load, and inability to achieve real-time active suppression of high-frequency vibrations caused by unsteady flow fields.

[0005] To achieve the above objectives, this application proposes an active vibration compensation control method for propulsion systems based on RBF load mapping. The active vibration compensation control method for propulsion systems based on RBF load mapping includes: The time-domain data of pulsating pressure in the unsteady flow field around the thruster are collected, and the flow field pulsating pressure time-domain matrix is ​​constructed by extracting flow field features and aligning them with time and space. Based on the Wendland C2 type radial basis function with local curvature adaptive compactly supported domain, and combined with the geometric nonlinear characteristics of the fluid-structure interaction interface, a fluid-structure interaction mapping kernel function matrix is ​​constructed. The time-domain matrix of the flow field pulsating pressure and the fluid-structure interaction mapping kernel function matrix are input into the RBF neural network mapping model to obtain the time-domain sequence of the nodal load on the wet surface of the thruster. The time-domain sequence of the load at the wet surface nodes is processed by a sparse identification algorithm. Based on Lasso regularization constraints, the key nodes that contribute the most to the vibration are selected to obtain a sparse key node set and a sparse load time-domain sequence. The sparse load time-domain sequence is transformed in the frequency domain, and combined with the spatial distribution characteristics of the sparsified key node set, a structural surface load spectrum matrix is ​​generated. The surface load spectrum matrix of the structure is input into the anti-phase control network of the piezoelectric actuator array to generate an anti-phase vibration compensation drive signal. The anti-phase vibration compensation drive signal is then applied to the piezoelectric actuator array of the thruster structure to achieve active vibration compensation control of the propulsion system.

[0006] Furthermore, to achieve the above objectives, this application also proposes an active vibration compensation control device for a propulsion system based on RBF load mapping. The active vibration compensation control device for a propulsion system based on RBF load mapping includes: The module is used to collect the time-domain data of the pulsating pressure of the unsteady flow field around the thruster, and construct the time-domain matrix of the pulsating pressure of the flow field by extracting flow field features and aligning them with time and space. The construction module is also used to construct a fluid-structure interaction mapping kernel function matrix based on the Wendland C2 type radial basis function of the locally curvature adaptive compactly supported domain, combined with the geometric nonlinear characteristics of the fluid-structure interaction interface. The mapping module is used to input the time-domain matrix of the flow field pulsating pressure and the fluid-structure interaction mapping kernel function matrix into the RBF neural network mapping model to obtain the time-domain sequence of the nodal load on the wet surface of the propeller. The filtering module is used to process the time-domain sequence of the wet surface node load through a sparse identification algorithm, and to filter the key nodes that contribute the most to the vibration based on Lasso regularization constraints, so as to obtain a sparse key node set and a sparse load time-domain sequence. The transformation module is used to perform frequency domain transformation on the sparse load time domain sequence and generate a structural surface load spectrum matrix by combining the spatial distribution characteristics of the sparsified key node set. The control module is used to input the surface load spectrum matrix of the structure into the anti-phase control network of the piezoelectric actuator array, generate an anti-phase vibration compensation drive signal, and apply the anti-phase vibration compensation drive signal to the piezoelectric actuator array of the thruster structure to realize active vibration compensation control of the propulsion system.

[0007] In addition, to achieve the above objectives, this application also proposes a non-transitory computer-readable storage medium storing a computer program thereon, which, when executed by a processor, implements the active vibration compensation control method for propulsion systems based on RBF load mapping as described above.

[0008] This application proposes one or more technical solutions that, by real-time acquisition of the pulsating pressure time-domain data of the unsteady flow field around the thruster and combining it with flow field feature extraction and spatiotemporal alignment mechanisms, achieves an efficient conversion of the pulsating pressure from physical acquisition to structured matrix representation. This is based on the Wendland model of locally curvature adaptive compactly supported domains. By incorporating the C2-type radial basis function and the geometric nonlinear characteristics of the fluid-structure interaction interface, a fluid-structure interaction mapping kernel function matrix is ​​constructed, possessing both geometrically adaptive compact support properties and local coupling strength sensing capabilities. This allows the pulsating pressure at discrete nodes in the fluid domain to be transmitted to the structural surface along the true geodesic path of the wetted surface, significantly suppressing load transfer distortion caused by the Euclidean distance assumption in high-curvature regions. The time-domain matrix of the pulsating pressure in the flow field and the fluid-structure interaction mapping kernel function matrix are input into an RBF neural network mapping model. Utilizing the local interpolation characteristics of the radial basis function and the nonlinear fitting capability of the neural network, real-time and accurate mapping of the pulsating pressure at fluid nodes to the dynamic load at the nodes on the wetted surface of the structure is achieved. A sparse identification algorithm based on Lasso regularization constraints is applied to the time-domain sequence of the load at the wetted surface nodes, automatically selecting the key nodes with the largest vibration contribution and eliminating redundant nodes, thus reducing the dimensionality of the full-dimensional load reconstruction problem to a smaller scale. The key node subspace significantly reduces the real-time computational load of the active control system while ensuring control accuracy. It performs frequency domain transformation on the sparse load time-domain sequence and, combined with the spatial distribution characteristics of the sparse key node set, generates a structural surface load spectrum matrix. This matrix is ​​input into the piezoelectric actuator array anti-phase control network to generate an anti-phase vibration compensation drive signal, which is directly applied to the piezoelectric actuator array of the thruster structure. Using a feedforward approach, controlled micro-vibrations with opposite phases and matched amplitudes are generated before the vibration response occurs, achieving real-time active suppression of high-frequency vibrations caused by unsteady flow fields. This solves the technical problems of existing propulsion system vibration control relying on lag-type PID strategies, inaccurate fluid-structure interaction load transfer, and excessive computational load. It achieves real-time active suppression of high-frequency thruster vibrations while ensuring computational accuracy, significantly improving the response speed, suppression accuracy, and engineering applicability of propulsion system vibration control. Attached Figure Description

[0009] The accompanying drawings, which are incorporated in and form part of this specification, illustrate embodiments consistent with this application and, together with the description, serve to explain the principles of this application.

[0010] To more clearly illustrate the technical solutions in the embodiments of this application or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, for those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0011] Figure 1 This is a flowchart illustrating an embodiment of the active vibration compensation control method for propulsion systems based on RBF load mapping in this application. Figure 2 This is a schematic diagram of the module structure of the active vibration compensation control device for a propulsion system based on RBF load mapping, according to an embodiment of this application.

[0012] The purpose, features, and advantages of this application will be further explained in conjunction with the embodiments and with reference to the accompanying drawings. Detailed Implementation

[0013] It should be understood that the specific embodiments described herein are merely illustrative of the technical solutions of this application and are not intended to limit this application.

[0014] To better understand the technical solution of this application, a detailed description will be provided below in conjunction with the accompanying drawings and specific implementation methods.

[0015] It should be noted that the executing entity in this embodiment can be a computing service device with data processing, network communication, and program execution functions, such as a tablet computer, personal computer, or mobile phone, or an electronic device capable of performing the above functions, such as an active vibration compensation control device for a propulsion system based on RBF load mapping. The following description uses an active vibration compensation control device for a propulsion system based on RBF load mapping as an example to illustrate this embodiment and the subsequent embodiments.

[0016] Based on this, embodiments of this application provide a propulsion system vibration active compensation control method based on RBF load mapping, referring to... Figure 1 , Figure 1 This is a flowchart illustrating the first embodiment of the active vibration compensation control method for propulsion systems based on RBF load mapping in this application.

[0017] In this embodiment, the active vibration compensation control method for propulsion systems based on RBF load mapping includes steps S10~S60: Step S10: Collect time-domain data of pulsating pressure in the unsteady flow field around the thruster, and construct the time-domain matrix of pulsating pressure in the flow field by extracting flow field features and aligning them with time and space.

[0018] It should be noted that unsteady flow field refers to the flow state in which physical quantities such as velocity and pressure change significantly over time in the flow field around the propeller when the propeller operates under complex incoming flow conditions such as the wake of the ship, the free surface, and maneuvering motions. Pulsating pressure time-domain data refers to a discrete-time series data set characterizing the time-varying characteristics of flow field pressure, collected by an array of pressure sensors arranged in the flow field around the propeller. Flow field feature extraction is a signal processing procedure that extracts the dominant pressure pulsating modes related to structural excitation from the raw pulsating pressure data and suppresses flow noise interference components. Spatiotemporal alignment is a data organization process that matches and normalizes the flow field pressure data with the propeller structural nodes in terms of time and space dimensions. The flow field pulsating pressure time-domain matrix is ​​a structured flow field pressure data organized in matrix form, output after flow field feature extraction and spatiotemporal alignment, where the row dimension corresponds to spatial sampling points and the column dimension corresponds to time sampling points.

[0019] In actual operation, a multiphysics sensor cluster is deployed in the flow field around the propeller using a circumferential array and an axial gradient. Specifically, a circumferential pressure sensor array is arranged within a range of 0.2D to 0.5D upstream of the propeller disk, where D is the propeller diameter, with a circumferential resolution of at least 16 points. The array is densely packed in the blade tip gaps and along the leading and trailing edges. An axial pressure sensor array is deployed along the axial direction with a 0.1D spacing gradient within a range of 0.1D to 1.0D downstream of the propeller. In addition to pressure sensors, the sensor cluster integrates turbulent pulsation velocity sensors and fluid temperature sensors to support the discrimination of local flow field conditions. Pulsation pressure time-domain data from each sensor channel are synchronously acquired at a sampling frequency of at least 20kHz, with each channel's data length covering at least 100 propeller rotation cycles.

[0020] In one feasible implementation, step S10 may include: deploying a cluster of multiphysics sensors in a circumferential array and axial gradient in the flow field around the thruster to collect pulsating pressure time-domain data of the unsteady flow field around the thruster; determining a local state discrimination factor of the flow field based on the turbulent pulsating velocity and fluid temperature, and removing outlier data from the pulsating pressure time-domain data based on the local state discrimination factor to obtain cleaned pulsating pressure time-domain data; performing multi-scale wavelet packet decomposition on the cleaned pulsating pressure time-domain data to obtain the dominant pressure pulsation component within a preset frequency band; and rotating the dominant pressure pulsation component based on the thruster shaft frequency encoded signal. Periodic synchronization alignment is performed to generate a phase-locked pressure sequence that is phase-locked with the blade rotation. The phase-locked pressure sequence is then normalized in spatiotemporal coordinates to construct a standardized pulsating pressure time-domain signal matrix. Singular value decomposition is performed on the standardized pulsating pressure time-domain signal matrix to extract the first K dominant mode pressure components, generating a compressed dominant mode pressure time-domain matrix, where K is a positive integer. Based on the wet surface node topology numbering of the propeller structure finite element model, a spatiotemporal mapping index is established between the compressed dominant mode pressure time-domain matrix and the structural nodes. The flow field pulsating pressure time-domain matrix corresponding to the spatiotemporal of the structural wet surface nodes is determined according to the spatiotemporal mapping index.

[0021] It should be noted that, as is understandable, a multiphysics sensor cluster refers to a composite sensor array unit that integrates the sensing functions of multiple physical quantities such as pressure, velocity, and temperature. These are typically manufactured using MEMS technology and are characterized by their small size, fast response, and ease of array deployment. A circumferential array refers to an array configuration where sensors are arranged uniformly or at specific angular intervals along the circumference of the propeller. Axial gradient deployment refers to a spatial configuration strategy where sensors are arranged along the propeller axis with non-equidistant or equidistant spacing, typically with denser placement in regions of rapid flow field changes.

[0022] It is understandable that turbulent pulsating velocity refers to the random pulsating component of instantaneous velocity relative to time-averaged velocity in a flow field, and is a key physical quantity characterizing turbulence intensity. Fluid temperature refers to the temperature of the medium in a local region of the flow field; temperature changes affect fluid density, viscosity, and pressure sensor calibration characteristics. The flow field local state discrimination factor is a comprehensive quantitative index used to evaluate the stability of the flow state and the reliability of data in a local region of the flow field. Outlier removal refers to the data cleaning operation that identifies and removes abnormal sampling points affected by strong interference or sensor malfunction based on the discrimination factor.

[0023] Understandably, wavelet packet decomposition is a more refined time-frequency analysis method than wavelet decomposition, recursively decomposing both the high-frequency and low-frequency components of the signal to form frequency bands of equal bandwidth. The pressure pulsation dominant component refers to the set of frequency band sub-signals containing the main propeller excitation frequency components, i.e., the shaft frequency and its harmonics, after wavelet packet decomposition. The flow noise interference component refers to the set of frequency band sub-signals containing incoherent interference components such as broadband random turbulent pulsations and hydrodynamic cavitation noise.

[0024] Understandably, the shaft frequency encoded signal refers to a pulse sequence signal generated by an encoder installed on the propeller shaft system, synchronized with the shaft rotation, with each pulse corresponding to a fixed angular increment of the shaft. Rotational cycle synchronization alignment is a signal processing technique that uses the shaft rotation period as a time reference and performs phase-locked superposition and averaging of pressure data collected at different times. The phase-locked pressure sequence is a pressure pulsation cycle sequence output after rotational cycle synchronization alignment, which strictly corresponds to the blade rotation phase.

[0025] Understandably, spatiotemporal coordinate normalization is a data preprocessing operation that unifies data from different sensor channels to the same spatial coordinate system and time scale. The standardized pulsating pressure time-domain signal matrix is ​​a structured pressure data matrix with unified dimensions and coordinate reference, output after spatiotemporal coordinate normalization. Singular value decomposition (SVD) is a matrix factorization method that decomposes a matrix into the product of left singular vectors, singular values, and right singular vectors; it is widely used in data dimensionality reduction and principal component extraction. The dominant modal pressure component is the pressure space-time coupled mode corresponding to the larger singular value after SVD, characterizing the main energy distribution pattern of flow field pressure pulsations.

[0026] It is understandable that the wet surface node topology number is the global set of numbers and adjacency information of all wet surface nodes in contact with the fluid domain in the finite element model of the propeller structure. The spatiotemporal mapping index is an index data structure that establishes the correspondence between spatial sampling points in the fluid domain and wet surface nodes in the structural domain, supporting data mapping and interpolation operations from the fluid domain to the structural domain.

[0027] In the specific implementation, the circumferential array of the multiphysics sensor cluster is arranged 0.3D upstream of the propeller disk, with 24 measuring points evenly distributed circumferentially, covering the entire 360° circumference. In the blade tip gap region, 4 additional dense measuring points are added at the corresponding position of each blade, improving the circumferential resolution to 96 points. The axial gradient array extends from 0.5D upstream to 1.0D downstream, with a total of 15 axial stations arranged, and the spacing gradually changes from 0.05D upstream to 0.1D downstream. The pressure sensors in the sensor cluster are piezoelectric or fiber optic, with a measurement range covering -50kPa to +200kPa and a frequency response range of 0.1Hz to 50kHz, meeting the wideband requirements for unsteady flow field pulsating pressure measurement.

[0028] Turbulent pulsating velocities are measured using a hot-wire anemometer or laser Doppler velocimeter integrated into a sensor cluster, with the sampling frequency synchronized with the pressure sensor. Fluid temperature is measured using thermocouples. The local state discrimination factor λ of the flow field is defined as: in, This represents the root mean square value of the turbulent fluctuation velocity. The local time-averaged flow velocity is obtained by moving average of the instantaneous flow velocity time series, with an averaging window of 10 rotation cycles. T is the measured fluid temperature. For reference temperature, the inlet flow temperature of the thruster or the design operating temperature is taken. α and β are preset weighting coefficients, for example, α=0.7, β=0.3.

[0029] A state discrimination threshold λ_th is set. When λ > λ_th, the data at that sampling point is determined to be affected by strong turbulent fluctuations or temperature drift, marked as abnormal data, and removed, resulting in cleaned pulsating pressure time-domain data. The state discrimination threshold λ_th is calibrated based on the flow field conditions. For example, if λ_th = 0.15, when the λ value of a certain sensor channel continuously exceeds the threshold for 5 rotation cycles, the data of that channel is determined to be abnormal, and the removal procedure is initiated, using spatial interpolation results from adjacent normal channel data to fill the gap.

[0030] Multi-scale wavelet packet decomposition was performed on the cleaned pulsating pressure time-domain data using the Daubechies 8th-order wavelet basis function with a decomposition level of J=5, resulting in 2^J=32 frequency band subspaces. Based on the characteristic frequency range of the thruster excitation (typically 1 to 10 times the shaft frequency), a subset of frequency bands containing the dominant excitation frequency was selected as the dominant pressure pulsation component, while the remaining terms were suppressed as flow noise interference components.

[0031] Wavelet packet decomposition is implemented using the Mallat pyramid algorithm. Let the original signal be s(t), and the decomposition level J=5, to obtain the wavelet packet coefficients W. j,n (k), where j=0,1,...,J is the decomposition level, and n=0,1,...,2 j-1 The nodes are numbered, and k is the discrete-time index. The dominant excitation frequency range is determined based on the thruster design parameters, such as the shaft frequency f. shaft If the frequency is 5Hz, then the dominant excitation frequency range is 5~50Hz. The set of wavelet packet nodes containing this frequency range is selected as the dominant component of the pressure pulsation. The coefficients of the remaining nodes are set to zero before reconstruction to obtain the denoised pressure signal.

[0032] The shaft frequency encoding signal is provided by a photoelectric encoder or magnetoelectric speed sensor mounted on the propeller shaft system, outputting N per revolution. p pulses, for example, N p=1024. Using the zero-phase moment of the shaft frequency encoded signal as the time reference, the dominant pressure pulsation components of each sensor channel are divided according to the rotation period and superimposed and averaged to generate a phase-locked pressure sequence that is phase-locked with the blade rotation, effectively eliminating the incoherent components of random turbulent pulsation.

[0033] Spatiotemporal coordinate normalization includes time coordinate normalization and spatial coordinate normalization. Spatial coordinate normalization transforms the physical coordinates of each sensor to a cylindrical coordinate system with the propeller center as the origin and the axis as the z-axis, and performs dimensionless processing. Time coordinate normalization maps the phase angle of the phase-locked sequence to the interval [0, 2π]. A standardized pulsating pressure time-domain signal matrix P∈R is constructed. (M×N) Where M is the number of spatial sampling points and N is the number of time sampling points per cycle.

[0034] Singular value decomposition is performed on the standardized pulsating pressure time-domain signal matrix to extract the top K dominant mode pressure components, i.e., retaining the components corresponding to the top K largest singular values, to generate the compressed dominant mode pressure time-domain matrix P. K ∈R (M×K) The selection of K is based on the cumulative singular value energy ratio criterion; for example, to retain more than 99% of the signal energy, K is a positive integer. After singular value decomposition, the cumulative singular value energy ratio is defined as: in, The percentage of energy for the Kth cumulative singular value. Let M be the k-th singular value, and min(M,N) be the smaller of the number of rows and columns of matrix P.

[0035] When η(K) first exceeds a preset energy threshold, such as 99%, extraction stops and the current K value is determined, achieving data dimensionality reduction and compression while retaining the core energy information of pressure pulsations in the flow field. For unsteady flow fields in the thruster, a K value in the range of 5 to 20 can capture more than 99% of the pressure pulsation energy.

[0036] The propeller structure finite element model is a discretized mesh model containing the propeller body, hub, and shaft connection structure, created using professional finite element software such as ANSYS and ABAQUS. All wetted surface nodes in contact with the fluid are extracted, and a unified topological numbering system is applied to these nodes in a global coordinate system, forming a node number list. For each fluid domain spatial sampling point in the compressed dominant modal pressure time-domain matrix, distance-weighted matching is performed between its normalized cylindrical coordinates and the coordinates of the structural wetted surface nodes. A corresponding flow field sampling point index relationship is established for each structural node, thus creating a spatiotemporal mapping index between the compressed dominant modal pressure time-domain matrix and the structural nodes. Finally, based on the established spatiotemporal mapping index, the spatiotemporal flow field pulsating pressure time-domain matrix corresponding to each structural wetted surface node is obtained through inverse distance-weighted interpolation.

[0037] Step S20: Based on the Wendland C2 type radial basis function of the locally curvature adaptive compactly supported domain, and combined with the geometric nonlinear characteristics of the fluid-structure interaction interface, construct the fluid-structure interaction mapping kernel function matrix.

[0038] It should be noted that the local curvature adaptive compactly supported domain is a localized approximation strategy that dynamically adjusts the range of the radial basis function support domain based on the local geometric curvature of the fluid-structure interaction interface. The Wendland C2-type radial basis function is a positive definite radial basis function with compact support characteristics. This function possesses C2 continuity and compact support properties, making it suitable for interpolation of large-scale scattered data. The geometric nonlinearity of the fluid-structure interaction interface refers to the characteristic of the load transfer path deviating from the linear assumption due to local curvature changes, surface distortion, and edge effects at the fluid-structure interaction interface, i.e., the wetted surface of the propeller. The fluid-structure interaction mapping kernel function matrix describes the mapping relationship between nodal pressure in the fluid domain and nodal load on the wetted surface of the structure; its elements are composed of the values ​​of the radial basis functions at the fluid-structure interaction interface.

[0039] In one feasible implementation, step S20 may include: extracting the set of wet surface element center coordinates and the set of element normal vectors from the finite element model of the propeller structure, and obtaining the set of fluid domain node coordinates; determining the three-dimensional Euclidean distance matrix between the wet surface element center and the fluid domain node based on the set of wet surface element center coordinates and the set of fluid domain node coordinates, and determining the local Gaussian curvature distribution and the average curvature distribution of the wet surface through local surface fitting based on the set of element normal vectors; determining the curvature adaptive compactly supported domain adjustment factor according to the local Gaussian curvature distribution and the average curvature distribution of the wet surface, and correcting the three-dimensional Euclidean distance matrix into an equivalent distance matrix along the geodesic of the wet surface based on the curvature adaptive compactly supported domain adjustment factor, wherein the curvature adaptive compactly supported domain adjustment factor characterizes the nonlinear distortion effect of the local curvature of the fluid-structure interaction interface on the load transfer path; determining the local strength index of fluid-structure interaction based on the turbulent kinetic energy density distribution and the local vorticity distribution at the fluid domain node, and adjusting the Wendland coefficient of the local curvature adaptive compactly supported domain according to the local strength index of fluid-structure interaction. The support domain radius distribution of the C2-type radial basis function is obtained, and the adjusted support domain radius distribution is obtained. Based on the equivalent distance matrix and the adjusted support domain radius distribution, a radial basis function kernel matrix is ​​constructed, and a set of linear equations with weight coefficients at the center of the wetted surface unit as the constraint point is established. The set of linear equations with weight coefficients is solved by the constrained least squares method to obtain the mapping weight coefficient vector from the fluid domain node to the center of the wetted surface unit. The mapping weight coefficient vector is fused with the radial basis function kernel matrix to obtain the fluid-structure interaction mapping kernel function matrix.

[0040] It can be understood that the set of coordinates of the wetted surface element centers is the set of coordinates of the geometric center points of all wetted surface elements in contact with the fluid domain in the global coordinate system within the finite element model of the propeller structure. The set of element normal vectors is the set of the outward normal unit vectors of each wetted surface element, representing the local orientation of the element surface. The set of coordinates of the fluid domain nodes is the set of coordinates of the discrete nodes of the fluid domain in the computational fluid dynamics (CFD) mesh in the global coordinate system.

[0041] Understandably, the three-dimensional Euclidean distance matrix is ​​a distance matrix that records the straight-line distance between the center of the wetted surface unit and the nodes of the fluid domain. Local surface fitting is a method that approximates the true geometry of a point cloud or mesh data using a low-order polynomial surface within a local neighborhood. Gaussian curvature is a differential geometric invariant describing the local bending characteristics of a surface, reflecting its elliptical, hyperbolic, or parabolic properties. Mean curvature is a differential geometric quantity describing the local average degree of bending of a surface. The curvature adaptive compact support domain adjustment factor is a correction coefficient that dynamically adjusts the support range of the radial basis function based on the magnitude of the local curvature; the greater the curvature, the smaller the adjustment factor, and the more contracted the compact support domain.

[0042] Understandably, the equivalent distance matrix is ​​the distance matrix after curvature adaptive correction, and its elements reflect the effective load transfer distance after considering geometric nonlinear effects. Turbulent kinetic energy density is a scalar physical quantity characterizing the intensity of turbulent fluctuations, defined as half the mean square value of the turbulent velocity fluctuations. Local vorticity distribution is a vector physical quantity distribution characterizing the local rotational characteristics of the flow field. The fluid-structure interaction local intensity index is a dimensionless index that comprehensively quantifies the local load transfer intensity at the fluid-structure interaction interface.

[0043] It is understandable that the support domain radius distribution is the set of parameters representing the effective range of the radial basis functions at each fluid domain node. Constrained least squares is an optimization method that introduces regularization constraints into the least squares objective function to handle ill-conditioned or underdetermined linear equation systems. The mapped weight coefficient vector is a set of weight coefficients describing the contribution of each node in the fluid domain to the load at the center of the wetted surface element.

[0044] In the specific implementation, wet surface information is extracted from the finite element model of the thruster structure. Wet surface elements are typically composed of surface elements of shell or solid elements, with element types including triangular and quadrilateral elements. For each wet surface element, its geometric center coordinates are calculated: in, Let be the geometric center coordinates of the i-th wetted surface element. The number of unit nodes, Let be the coordinates of the kth node of the i-th wet surface element.

[0045] The geometric center coordinates of all wetted surface elements form the wetted surface element center coordinate set. Similarly, the outward normal unit vector of each wetted surface element is extracted. The element normal vector is calculated through the cross product of the element node coordinates. The outward normal unit vectors of all elements form the element normal vector set. At the same time, the global coordinates of all discrete nodes in the fluid domain are read from the CFD calculation result file to form the fluid domain node coordinate set. Based on the wetted surface element center coordinate set and the fluid domain node coordinate set, the three-dimensional straight-line distance between any wetted surface element center and any fluid domain node is calculated. All distances are then sequentially filled into the corresponding positions in the matrix to obtain the three-dimensional Euclidean distance matrix.

[0046] Then, taking the center of each wet surface unit as the center, select the units in its predetermined neighborhood, and perform second-order polynomial local surface fitting on the geometric center coordinates of all units in the neighborhood. Based on the surface parameters obtained by fitting, solve for the first and second fundamental form coefficients of the unit position. Then, calculate the Gaussian curvature and mean curvature of the position through differential geometry formulas. After traversing all wet surface units, the complete local Gaussian curvature distribution and mean curvature distribution of the wet surface can be obtained.

[0047] The curvature-adaptive compact support region adjustment factor reflects the nonlinear distortion effect of local curvature on the load transfer path. The greater the curvature, the shorter the effective distance of load transfer, and the compact support region should be reduced accordingly. The curvature-adaptive compact support region adjustment factor is as follows: in, This is the curvature adaptive compact support domain adjustment factor corresponding to the i-th wetted surface unit. For the preset curvature sensitivity coefficient, such as =10, Let be the Gaussian curvature at the location of the i-th wetted surface unit. Let be the average curvature at the location of the i-th wetted surface unit.

[0048] Curvature adaptive compact support region adjustment factor The value range of is (0,1). When the local surface is a plane, that is... = When =0, =1, the equivalent distance is equal to the Euclidean distance, when the local curvature increases. A decrease in curvature shortens the equivalent distance and causes the compactly supported region to shrink. (Curvature sensitivity coefficient) Calibrated based on the typical curvature range of the wetted surface of the propeller, for example, for propeller blades. Take a value between 5 and 15.

[0049] Each distance element in the three-dimensional Euclidean distance matrix is ​​corrected according to the curvature adaptive compact support domain adjustment factor. The original straight-line Euclidean distance is multiplied by the adjustment factor to obtain the corrected equivalent distance. All the corrected distances are arranged in the original manner to form an equivalent distance matrix along the geodesic line of the wet surface.

[0050] Next, the turbulent kinetic energy density and vorticity vector modulus at each fluid domain node are read from the CFD calculation results. Based on these two physical quantities, the local strength index of fluid-structure interaction is calculated, which can be specifically expressed as: in, Let be the local strength index of the fluid-structure interaction at the j-th fluid domain node. The turbulent kinetic energy density at this node. This represents the average turbulent kinetic energy density at all fluid domain nodes. The vorticity modulus of this node. This represents the average vorticity modulus of all fluid domain nodes.

[0051] The initial support domain radius is adjusted based on the local strength index of the fluid-structure interaction. A larger local strength index indicates a stronger fluid-structure interaction load at that location, requiring an increased support domain radius to ensure interpolation accuracy. Therefore, the adjusted support domain radius can be expressed as: in, The adjusted support domain radius for the j-th fluid domain node. The initial support domain radius is preset. This is the strength adjustment factor. This is the average value of the local strength index of fluid-structure interaction at all nodes.

[0052] Substituting the adjusted support domain radius distribution into the Wendland C2 type radial basis function, and combining it with the obtained equivalent distance matrix, the radial basis function values ​​between each constraint point and the fluid domain node are determined. The values ​​are then filled into the matrix to obtain the radial basis function kernel matrix. Furthermore, a linear equation system of weight coefficients is established based on the load interpolation constraint at the center of the wet surface element. A constraint least squares objective function is constructed by introducing a Tikhonov regularization term. The stable mapping weight coefficient vector is obtained by solving the QR decomposition. Finally, the mapping weights are combined with the kernel matrix to obtain the complete fluid-structure interaction mapping kernel function matrix.

[0053] Step S30: Input the time-domain matrix of the flow field pulsating pressure and the fluid-structure interaction mapping kernel function matrix into the RBF neural network mapping model to obtain the time-domain sequence of the nodal load on the wet surface of the propeller.

[0054] It should be noted that the RBF neural network mapping model is a feedforward neural network that uses radial basis functions as the activation functions of the hidden layers. Its mapping from the hidden layers to the output layers is a linear weighted combination, which has global approximation capability and fast learning characteristics. The wet surface nodal load time-domain sequence is a discrete-time series data set that characterizes the dynamic load on each node of the wet surface of the thruster structure over time after being mapped by the RBF neural network.

[0055] In actual implementation, the time-domain matrix of the fluctuating pressure in the flow field is windowed with an overlap rate of 50% according to a preset time window length, generating multiple short time series segments. The fluid-structure interaction mapping kernel function matrix is ​​used as the initial hidden layer weights to initialize the RBF neural network mapping model. The network topology is as follows: the input layer dimension equals the number of fluid domain nodes N. f The hidden layer is activated using radial basis functions, and the number of hidden layer nodes is equal to the number of wet surface elements N. w The output layer dimension is equal to the number of wet surface nodes N. s N s ≥N w This is because a surface element may correspond to multiple structural nodes. The center coordinates of the wetted surface element are taken as the center of the hidden layer, and the initial width parameter is determined by the radius distribution of the support domain. Then, short-time series fragments and the fluid-structure interaction mapping kernel function matrix are input in parallel into the initialized RBF neural network mapping model.

[0056] In one feasible implementation, step S30 may include: performing overlapping windowing processing on the flow field pulsating pressure time-domain matrix according to a preset time window length to generate multiple short time sequence segments; initializing the network topology of the RBF neural network mapping model using the fluid-structure interaction mapping kernel function matrix as the initial hidden layer weights to obtain the initialized RBF neural network mapping model; inputting the short time sequence segments and the fluid-structure interaction mapping kernel function matrix in parallel into the initialized RBF neural network mapping model, and updating the network output layer weights in real time using an online recursive least squares algorithm to obtain the updated network output layer weights; based on the updated... The new network output layer weights and the fluid-structure interaction mapping kernel function matrix are used to perform radial basis interpolation on the fluctuating pressure of the fluid nodes within each time window to determine the dynamic pressure distribution vector at the center of the wet surface unit. Based on the dynamic pressure distribution vector and the area attribute of the corresponding wet surface unit, the concentrated force vector of the surface unit is obtained. The concentrated force vector of the surface unit is weighted and distributed to all structural nodes belonging to the corresponding wet surface unit according to the shape function or area coordinate of the corresponding wet surface unit to obtain the equivalent dynamic node force on each structural node. The equivalent dynamic node forces of the structural nodes corresponding to all wet surface units are aggregated and stacked in time sequence to generate the time domain sequence of the thruster wet surface node load.

[0057] It should be noted that the time window length is the duration of a short data segment during online learning of the RBF neural network mapping model, typically ranging from 10 to 50 rotation cycles. Overlapping windowing is a data segmentation technique that divides a continuous time series into segments of fixed length with partial overlap between adjacent segments, used to improve temporal resolution and smooth window boundary effects.

[0058] It is understandable that the RBF neural network mapping model is a three-layer feedforward neural network architecture, including an input layer, hidden layers, and an output layer. Unlike traditional RBF networks, the hidden layer weights in this implementation are not randomly initialized or fixed through offline training, but are directly used as initial weights based on the fluid-structure interaction (FSI) mapping kernel function matrix. This allows the network to possess physically interpretable FSI priors during the initialization phase, significantly shortening the online convergence time. Network topology refers to the connection methods and dimensional configurations between the layers of the neural network. The online recursive least squares algorithm is a recursive least squares method that updates model parameters sample by sample, suitable for adaptive learning of real-time data streams.

[0059] It is understandable that radial basis interpolation is a scattered data interpolation method based on radial basis functions, which estimates the function values ​​of unknown points by weighted combination of function values ​​at known points. The dynamic pressure distribution vector is the pressure vector at the center of the wetted surface element that changes over time. The concentrated force vector of the surface element is the equivalent concentrated force vector obtained by multiplying the pressure at the element center by the element area. The equivalent dynamic nodal force is the nodal force vector obtained by distributing the concentrated force of the surface element to each node of the element according to the shape function.

[0060] In the specific implementation, the preset time window length T w Take 20 rotational cycles. For the shaft frequency f shaft =5Hz thruster, T w =4s. An overlap rate of 50% means that adjacent time windows overlap by 10 rotation cycles, or 2s. The total number of time windows equals the total time length divided by the time window step size. The time window step size is the time length of the non-overlapping part, i.e., Tw / 2. The total time length corresponds to the total physical duration of the CFD unsteady computation.

[0061] During the initialization of the RBF neural network mapping model, the number of hidden layer nodes is N. w This refers to the number of wet surface elements. Each hidden layer node corresponds to a radial basis function center, and the center coordinates are taken from the coordinates of the wet surface element center. The initial width parameter σ... i Take the radius R of the support domain i 0.5 times, that is, σ i =R i / 2.

[0062] Short-time sequence fragments and the fluid-structure interaction (FSI) mapping kernel function matrix are input in parallel into the initialized RBF neural network mapping model. The network output layer is then updated in real time using an online recursive least squares algorithm. The core of the online recursive least squares algorithm is to use a forgetting factor to weight the old data while retaining the update weights of the new data, ensuring that the mapping model can quickly track vibration deviations caused by changes in the propulsion system load.

[0063] Short-time series segments and the fluid-structure interaction mapping kernel function matrix are input in parallel into the initialized RBF neural network mapping model. For each time t in the short-time series segment, the hidden layer output is first calculated. Simultaneously, calibration pressure sensors deployed on the propeller structure surface, located at key locations such as the hub and blade root, acquire measured calibration values ​​of the wetted surface unit centers. The recursive update objective of the output layer weights is defined as minimizing the calibration pressure estimation error. The weight correction direction and step size are dynamically adjusted using a forgetting factor, and the output layer weights are updated in real-time within each time window, resulting in the updated output layer weights. These updated output layer weights incorporate the flow field state evolution information within the current time window, enabling the mapping model to adaptively track unsteady flow fields.

[0064] The concentrated force vector of the surface element acts at the element center, while the degrees of freedom of the structural finite element model are defined at the element nodes. Therefore, the concentrated force of the element needs to be equivalently distributed to each node to form a nodal force vector. The distribution method must follow the principle of virtual work equivalence, that is, the work done by the concentrated force of the element on any virtual displacement is equal to the sum of the work done by the equivalent nodal forces at each node on the corresponding nodal virtual displacement. Depending on the element type, there are two specific methods: for triangular elements, area coordinates (i.e., centroid coordinates) are used for weighted distribution; for quadrilateral elements, bilinear shape function weighted distribution is used to obtain the equivalent dynamic nodal force corresponding to each structural node. After completing the interpolation calculation and nodal force distribution for all time windows, the nodal forces at all times are stacked in temporal order to obtain the complete time-domain sequence of the nodal loads on the wetted surface of the thruster.

[0065] Step S40: Process the time-domain sequence of the wet surface node load using a sparse identification algorithm, and select the key nodes that contribute the most to vibration based on Lasso regularization constraints to obtain a sparse key node set and a sparse load time-domain sequence.

[0066] It should be noted that sparse identification algorithms are compressed sensing algorithms that identify and extract the few key elements that contribute most to the dynamic response of a system from a large-scale dataset. Lasso regularization is a regularization technique that introduces an L1 norm penalty term in regression analysis to induce sparse solutions, compressing unimportant regression coefficients to zero through soft thresholding. The sparse key node set is the set of wet surface node numbers that contribute most to the vibration response of the thruster structure, selected after sparse identification. The sparse load time-domain sequence is a dimensionality-reduced time-domain data sequence containing only the load components corresponding to the sparse key node set.

[0067] In one feasible implementation, step S40 may include: extracting the first N order main wet surface mode shape vectors of the finite element model of the thruster structure, where N is a positive integer; determining the modal participation factor vectors of each wet surface node based on the main wet surface mode shape vectors, and constructing a spatial feature matrix weighted by the modal participation factors; determining the temporal energy density of each node load component in the time-domain sequence of the wet surface node loads, and constructing a regression target vector with the time-domain energy density as the response variable; using the spatial feature matrix as the prediction variable and the regression target vector as the response variable, establishing an L1 positive... A Lasso sparse regression model with regularization constraints is constructed, and an initial intensity coefficient for regularization is set. The Lasso sparse regression model is iteratively optimized using the coordinate descent method. During the iteration process, a soft threshold operation is applied to the mapping weights. After the iteration converges, the node indices corresponding to the non-zero elements of the mapping weights are extracted to generate a sparse key node set. The vibration contribution ranking of each node in the sparse key node set is determined. Based on the vibration contribution ranking and the sparse key node set, the load components of the corresponding key nodes are extracted from the time-domain sequence of the wet surface node loads to obtain the sparse load time-domain sequence.

[0068] It should be noted that the main wet surface mode shape vector is the shape vector of the first N modes corresponding to the nodal degrees of freedom of the wet surface in the finite element modal analysis of the structure, characterizing the vibration mode distribution of the structure at each natural frequency. The modal participation factor vector is a normalized vector describing the relative participation of each node in each mode, reflecting the modal coupling characteristics of the nodal vibration response.

[0069] Understandably, time-domain energy density is the cumulative energy of the load signal per unit time, reflecting the intensity level of the dynamic load on that node. The regression objective vector is a response vector composed of the normalized values ​​of the time-domain energy density of each node, serving as the fitting target for Lasso regression. L1 regularization constraint is a constraint form that introduces a penalty term for the L1 norm of the parameter to be solved into the optimization objective function, possessing the characteristic of inducing sparse solutions. Soft thresholding is a thresholding operation that compresses parameters with absolute values ​​less than a threshold to zero and shrinks the remaining parameters towards zero.

[0070] Understandably, coordinate descent is an iterative optimization algorithm that performs a one-dimensional search along each coordinate direction, suitable for convex optimization problems with L1 penalties, such as Lasso. The sparsified key node set is a subset of wet surface node numbers that significantly contribute to the structural vibration response, retained after sparse identification. Vibration contribution ranking is a priority ranking based on a comprehensive evaluation of the modal participation degree and load energy level of each key node.

[0071] In practical implementation, modal analysis employs the Block Lanczos method or the Subspace iteration method to extract the first N=30 wetted surface modes of the thruster structure. The mode shape vectors need to be mass-normalized. For large and complex thruster structures, modal synthesis or automatic multi-stage substructures are used to improve modal calculation efficiency. Based on the main wetted surface mode shape vectors, the modal participation factor vector (MPF) of each wetted surface node is determined. For the j-th wetted surface node, its MPF is... j ∈R N The nth element is: in, This represents the mode shape value of the nth mode at the j-th node. This represents the total number of nodes on the wet surface. This represents the mode shape value of the nth mode at the kth node.

[0072] Then, a spatial feature matrix is ​​constructed by weighting the modal participation factors, where the spatial feature matrix X∈R Ns×N The j-th row and n-th column is X j,n =MPF j,n Where j corresponds to the node number on the wetted surface, and n corresponds to the modal order. The temporal energy density of each node load component in the time-domain sequence of the wetted surface node loads is determined, and the temporal energy density E of the j-th node is given. j The target vector can be obtained by integrating the square of the time series of the node load and dividing it by the total duration. Then, the time-domain energy density of all nodes is normalized and a regression target vector is constructed.

[0073] By using the spatial feature matrix as the predictor variable and the regression target vector as the response variable, a Lasso sparse regression model with L1 regularization can be established. The initial regularization strength coefficient λ is set. Lasso =0.1, this value can be adaptively adjusted according to the total number of nodes and the desired sparsity.

[0074] Lasso regularization strength coefficient λ Lasso The selection adopts a cross-validation strategy. The data is divided into training and validation sets according to time, with a ratio of 8:2, in λ∈

[10] . -4 10 0100 candidate values ​​are logarithmically uniformly sampled within the range. A Lasso model is trained for each candidate value, and the mean square error of prediction (MSE) is calculated on the validation set. The λ that minimizes the validation MSE is selected as the optimal regularization strength.

[0075] An adaptive threshold strategy is used to control the size of the sparsified key node set. A target sparsity ρ is set. target =10%, if the initial λ Lasso The obtained sparsity ρ≠ρ target Then adjust λ using a binary search. Lasso Until |ρ-ρ target |<1%. This adaptive strategy ensures that the sparsification process achieves the expected computational dimensionality reduction effect under different operating conditions and structural configurations.

[0076] The vibration contribution of each node in the sparsified critical node set is ranked, and the vibration contribution is determined by comprehensively considering the modal participation factor and load energy density. After ranking the vibration contribution, the top 20% of critical nodes are subject to key monitoring and high-precision control, while the remaining 80% of nodes are estimated using low-frequency sampling or interpolation, thereby optimizing the allocation of computing resources.

[0077] Step S50: Perform frequency domain transformation on the sparse load time domain sequence, and generate a structural surface load spectrum matrix by combining the spatial distribution characteristics of the sparse key node set.

[0078] It should be noted that frequency domain transformation is a mathematical transformation method that converts a time-domain signal into a frequency-domain representation, revealing the frequency components and energy distribution of the signal. The structural surface load spectrum matrix is ​​a structured frequency domain data organized in matrix form, generated after frequency domain transformation and spatial feature fusion, and its dimensions cover three levels: nodes, frequency points, and frequency bands.

[0079] In one feasible implementation, step S50 may include: performing a windowed short-time Fourier transform on each key node component in the sparse load time-domain sequence to obtain a complex spectrum sequence of each key node at different times; extracting amplitude spectrum components and phase spectrum components from the complex spectrum sequence, and dividing the frequency band based on the thruster excitation characteristic frequency to generate the spectral energy density distribution within each frequency band; constructing a frequency domain-space joint index based on the three-dimensional spatial coordinates of each node in the sparse key node set, wherein the index key of the frequency domain-space joint index is a spatial location code generated based on the three-dimensional spatial coordinates of the node. The index value is the spectral identifier of the corresponding node; the amplitude spectral component and the phase spectral component are bound to the corresponding frequency domain-space joint index to establish a frequency domain-space joint index table; according to the spectral energy density distribution, the dominant excitation frequency and the corresponding harmonic components are identified, and the spectral components of each key node at the dominant excitation frequency are topologically aggregated to obtain the topological aggregation result; based on the topological aggregation result and the frequency domain-space joint index table, a frequency domain-space joint representation tensor is constructed; the frequency domain-space joint representation tensor is reorganized according to the node-frequency-band dimension to obtain the structural surface load spectrum matrix.

[0080] Understandably, the windowed short-time Fourier transform (STFT) is a time-frequency analysis method that performs Fourier transform on a segmented, windowed time-domain signal, balancing time and frequency resolution. The complex spectrum sequence is the complex numerical time-frequency representation output by the STFT, containing amplitude and phase information. The amplitude spectral component is the magnitude of the complex spectrum, reflecting the energy intensity of each frequency component. The phase spectral component is the argument of the complex spectrum, reflecting the phase relationship between the frequency components.

[0081] Understandably, frequency band division is the operation of discretizing a continuous spectrum into several characteristic frequency bands based on the thruster excitation characteristics. The spectral energy density distribution is the spatial distribution of accumulated energy within each frequency band, reflecting the distribution characteristics of vibration energy across different frequency ranges and spatial locations. The frequency domain-spatial joint index is a multi-dimensional index structure that simultaneously encodes frequency features and spatial locations, supporting efficient multi-condition queries.

[0082] It is understandable that spatial location encoding is a dimensionality reduction encoding method that maps three-dimensional spatial coordinates to one-dimensional index keys. Spectral identifiers are tag information containing metadata about the spectral characteristics of nodes. Topological aggregation is an operation that clusters and merges spectral data based on the spatial adjacency relationships of nodes, aiming to reduce data dimensionality while preserving spatial distribution characteristics. The frequency-space joint representation tensor is a multi-dimensional data structure that integrates frequency, spatial, and temporal information.

[0083] In actual implementation, a windowed short-time Fourier transform is performed column-by-column on each critical node component in the sparse load time-domain sequence. Let the load time sequence of the j-th critical node be f. j(t), a Hanning window w(t) is used for windowing, with a window length of 2 to 4 rotation cycles, i.e., T. window The time interval is 0.4~0.8s, the overlap rate is 75%, and the frequency resolution is Δf=1 / T. window The complex spectrum sequences of each key node at different times are obtained.

[0084] Extracting the amplitude spectral component |F from a complex spectrum sequence j (τ,ω)| and phase spectral component arg(F) j (τ,ω)). Frequency band division is based on the thruster excitation characteristic frequency. The thruster excitation characteristic frequency includes the shaft frequency f. shaft and its harmonics, such as leaf frequency f blade =Z·f shaft Z represents the number of blades and the low-frequency modulation component caused by the uneven flow.

[0085] The frequency band division boundaries are dynamically adjusted based on the actual operating conditions of the propeller. Under mooring conditions, the main excitations are the shaft frequency and blade frequency, along with their lower harmonics. Under cruise conditions, low-frequency modulation components caused by wake inhomogeneities need to be considered, i.e., 0.2~0.8 times the shaft frequency. Under maneuvering conditions, a transient impact frequency band needs to be added, i.e., >5 times the blade frequency. The dynamic adjustment of the frequency band division is automatically triggered by the operating condition identification module based on the shaft frequency encoded signal and the flow field state discrimination factor. In this embodiment, the frequency band division scheme is set as follows: low-frequency band (0.5~2 times the shaft frequency), shaft frequency band (2~5 times the shaft frequency), blade frequency band (Z-2~Z+2 times the shaft frequency), and high-frequency band (>Z+2 times the shaft frequency). This generates the spectral energy density distribution within each frequency band. For the b-th frequency band, the spectral energy density of the j-th critical node in that frequency band is the ratio of the sum of squares of the amplitude spectral components of all frequency points in that frequency band to the bandwidth. By traversing all critical nodes, the spectral energy density distribution of the entire sparse critical node set in each frequency band can be obtained.

[0086] Based on the three-dimensional spatial coordinates (x, y) of each node in the sparse key node set j ,y j ,z j The frequency-space joint index is constructed by encoding the three-dimensional coordinates using a Z-order space-filling curve to generate a one-dimensional spatial location code, which serves as the index key. The index value stores the amplitude and phase metadata of each frequency band of the corresponding node, i.e., the spectrum identifier. The Z-order curve has good spatial locality preservation, and the encoded values ​​of adjacent spatial locations are usually similar, supporting spatial proximity retrieval based on range queries. By sequentially binding the amplitude spectral components and phase spectral components to the corresponding frequency-space joint index, the frequency-space joint index table can be constructed.

[0087] Based on the spectral energy density distribution of each frequency band, frequency points with an energy proportion exceeding 5% of the total energy are selected as dominant excitation frequencies. Simultaneously, their corresponding harmonic components are extracted. Based on the spatial adjacency grid of key nodes, topological aggregation is performed on the spectral components of nodes at the same dominant frequency point. Spectral components of adjacent nodes with a spatial distance less than a set threshold and a spectral feature similarity greater than a set threshold are merged into a single topological aggregation unit. The central spatial coordinates and average spectral features of each unit are retained, yielding the topological aggregation result. Based on the topological aggregation result and the constructed frequency-space joint index table, a three-dimensional frequency-space joint representation tensor is constructed according to the three dimensions of frequency point, spatial node, and frequency band. This tensor is then rearranged according to the node-frequency-frequency band dimensions, ultimately resulting in a tensor of dimension N. s '×N f ×N b The structural surface load spectrum matrix, where N s 'N' represents the total number of nodes in the sparse key node set. f N represents the total number of dominant incentive frequencies. b This represents the total number of frequency bands divided. The spatial distance threshold for topological aggregation is determined based on the thruster structural characteristics and actuator arrangement spacing. For large propeller thrusters, the threshold is set to 0.1 to 0.2 times the blade radius; for waterjet thrusters, the threshold is set to 0.05 to 0.1 times the nozzle diameter.

[0088] Step S60: Input the surface load spectrum matrix of the structure into the anti-phase control network of the piezoelectric actuator array to generate an anti-phase vibration compensation drive signal, and apply the anti-phase vibration compensation drive signal to the piezoelectric actuator array of the thruster structure to realize active vibration compensation control of the propulsion system.

[0089] It should be noted that the piezoelectric actuator array anti-phase control network is a frequency-domain feedforward actuator drive signal generation network. Its core principle is to generate drive signals at each frequency point that are 180 degrees out of phase with the load on the structural surface and proportional in amplitude, so that the controlled vibration generated by the piezoelectric actuator cancels out the structural vibration caused by the flow field excitation. Unlike the ex-post error correction mode of traditional PID feedback control, the feedforward control of this implementation can generate compensation action before the vibration response occurs, fundamentally eliminating control lag.

[0090] The anti-phase vibration compensation drive signal is a control drive signal whose amplitude matches the target vibration signal in the frequency domain but whose phase is opposite. It generates a compensating force that cancels out the structural vibration through the inverse piezoelectric effect of the piezoelectric actuator. The residual vibration response feedback sequence is the sequence of actual vibration response signals collected by an array of acceleration sensors arranged on the surface of the structure after the control drive is applied.

[0091] This embodiment provides a propulsion system vibration active compensation control method based on RBF load mapping. By real-time acquisition of the time-domain data of pulsating pressure in the unsteady flow field around the propeller and combining flow field feature extraction and spatiotemporal alignment mechanisms, it achieves efficient conversion of the pulsating pressure from physical acquisition to structured matrix representation. This is based on the Wendland model of locally curvature adaptive compactly supported domains. By incorporating the C2-type radial basis function and the geometric nonlinear characteristics of the fluid-structure interaction interface, a fluid-structure interaction mapping kernel function matrix is ​​constructed, possessing both geometrically adaptive compact support properties and local coupling strength sensing capabilities. This allows the pulsating pressure at discrete nodes in the fluid domain to be transmitted to the structural surface along the true geodesic path of the wetted surface, significantly suppressing load transfer distortion caused by the Euclidean distance assumption in high-curvature regions. The time-domain matrix of the pulsating pressure in the flow field and the fluid-structure interaction mapping kernel function matrix are input into an RBF neural network mapping model. Utilizing the local interpolation characteristics of the radial basis function and the nonlinear fitting capability of the neural network, real-time and accurate mapping of the pulsating pressure at fluid nodes to the dynamic load at the nodes on the wetted surface of the structure is achieved. A sparse identification algorithm based on Lasso regularization constraints is applied to the time-domain sequence of the load at the wetted surface nodes, automatically selecting the key nodes with the largest vibration contribution and eliminating redundant nodes, thus reducing the dimensionality of the full-dimensional load reconstruction problem to a smaller scale. The key node subspace significantly reduces the real-time computational load of the active control system while ensuring control accuracy. It performs frequency domain transformation on the sparse load time-domain sequence and, combined with the spatial distribution characteristics of the sparse key node set, generates a structural surface load spectrum matrix. This matrix is ​​input into the piezoelectric actuator array anti-phase control network to generate an anti-phase vibration compensation drive signal, which is directly applied to the piezoelectric actuator array of the thruster structure. Using a feedforward approach, controlled micro-vibrations with opposite phases and matched amplitudes are generated before the vibration response occurs, achieving real-time active suppression of high-frequency vibrations caused by unsteady flow fields. This solves the technical problems of existing propulsion system vibration control relying on lag-type PID strategies, inaccurate fluid-structure interaction load transfer, and excessive computational load. It achieves real-time active suppression of high-frequency thruster vibrations while ensuring computational accuracy, significantly improving the response speed, suppression accuracy, and engineering applicability of propulsion system vibration control.

[0092] Based on the first embodiment of this application, in the second embodiment of this application, the content that is the same as or similar to that in the first embodiment described above can be referred to the above description and will not be repeated hereafter. Based on this, step S60 includes steps S601 to S608: Step S601: Normalize the amplitude of the surface load spectrum matrix of the structure, and perform phase unwinding processing on the phase spectrum of each key node component channel by channel, encoding it into a frequency domain conditional tensor containing amplitude and phase channels.

[0093] Understandably, amplitude normalization is a preprocessing operation that divides each element of the spectrum matrix by its global maximum absolute value, mapping the data range to [-1, 1]. Phase unwrapping is a signal processing technique that eliminates 2π phase angle jumps and restores continuous phase changes. A frequency domain conditional tensor is a structured frequency domain data tensor that simultaneously encodes amplitude and phase information.

[0094] In the specific implementation, for each element X(i,j,k) in the surface load spectrum matrix of the structure, where i is the node index, j is the frequency point index, and k is the frequency band index, amplitude normalization is performed to obtain the normalized amplitude X'(i,j,k)=X(i,j,k) / max(|X|), where max(|X|) is the maximum absolute value in the entire spectrum matrix. Then, the phase component is extracted for each key node and each frequency band. Since the principal phase value obtained by Fourier transform is restricted to the interval (-π,π], discontinuous jumps will occur when the actual phase change exceeds 2π. Therefore, the least squares phase unwrapping algorithm is used to restore the continuous phase change relationship channel by channel to eliminate phase jump errors. Finally, the normalized amplitude and the unwrapped phase are treated as two independent channels and concatenated to obtain a dimension N. s '×N f ×N b The frequency domain conditional tensor is 2 × 2, where the last dimension corresponds to the amplitude channel and the phase channel.

[0095] Step S602: Obtain the spatial layout coordinates and installation direction angle of each actuator channel in the piezoelectric actuator array on the surface of the thruster structure, and determine the spatial coupling strength between each actuator and each node in the sparse key node set.

[0096] It should be noted that the piezoelectric actuator array is a distributed drive system composed of multiple piezoelectric ceramic actuators arranged in a specific spatial layout, with each channel controlled independently. The installation orientation angle is the angular relationship between the actuator's main driving direction and the local coordinate system of the structure. Spatial coupling strength is an indicator that quantifies the actuator's effectiveness in controlling the vibration of structural nodes, taking into account both spatial distance and the degree of directional matching.

[0097] In the specific implementation, based on the finite element modal analysis results of the thruster structure, the structural local stiffness matrix and vibration transfer function of each actuator's operating position are extracted. For the m-th piezoelectric actuator and the j-th key node, their spatial coupling strength C mj The spatial distance attenuation term is obtained by multiplying the spatial distance attenuation term and the orientation matching term. The spatial distance attenuation term uses the geodesic distance between the actuator and the node as a variable and is calculated using a Gaussian function; the greater the distance, the lower the coupling strength. The orientation matching term is calculated by the cosine of the angle between the actuator's installation direction and the node's principal vibration direction; the closer the orientation match, the higher the coupling strength. By traversing all actuator-key node combinations, the final result is a term with dimension N.a ×N s The spatial coupling strength matrix of ', where N a This represents the total number of channels in the piezoelectric actuator array.

[0098] Step S603: Based on the spatial coupling strength, construct the actuator-node association weight matrix. The dimension of the actuator-node association weight matrix is ​​determined according to the number of actuator channels and the number of key nodes.

[0099] It should be noted that the actuator-node association weight matrix is ​​a normalized weight matrix that describes the proportion of each actuator's contribution to the vibration compensation of each node.

[0100] In the specific implementation, the spatial coupling strength between all actuators and nodes is normalized. For each critical node, the sum of the coupling strengths of all actuators is normalized to 1, as shown in the following formula: in, The compensation contribution weight for the m-th actuator to the j-th critical node.

[0101] After normalization, the actuator-node association weight matrix is ​​obtained. The weight is positively correlated with the spatial coupling strength. The higher the spatial coupling strength, the greater the weight allocated to the corresponding actuator, ensuring that the vibration compensation energy is concentrated on the control path with high contribution.

[0102] Step S604: Input the frequency domain conditional tensor and the actuator-node association weight matrix into the inverse phase control network, and map the node frequency domain features into the preliminary driving spectrum of each actuator channel through a fully connected mapping layer.

[0103] It should be noted that the anti-phase control network is a neural network or adaptive filter that realizes the mapping from frequency domain features to the driving spectrum. The anti-phase control network includes a fully connected mapping layer, a frequency domain amplitude inversion layer, and an output normalization layer. The fully connected mapping layer takes the normalized amplitude and continuous phase of the key nodes as input, and performs weighted aggregation of the features of each node according to the actuator-node association weights to generate the initial frequency domain feature output corresponding to each actuator channel.

[0104] In the specific implementation, the frequency domain condition tensor is expanded along the node dimension and multiplied with the actuator-node association weight matrix of dimension Na×Ns' to obtain the weighted frequency domain features of each actuator channel corresponding to each frequency point and frequency band. After retaining the dimensions of the amplitude channel and the phase channel, it is input into two nonlinearly activated fully connected layers to output the preliminary driving spectrum. The preliminary driving spectrum completes the spatial dimension transformation from node load features to actuator driving features.

[0105] Step S605: In the frequency domain amplitude inversion layer of the anti-phase control network, a phase inversion operation is performed on the preliminary driving spectrum, and amplitude matching correction is performed according to the electromechanical coupling efficiency and structural frequency response function of the piezoelectric actuator to obtain the corrected frequency domain driving spectrum.

[0106] Understandably, the frequency domain amplitude inversion layer is the core functional layer in the network that performs phase inversion and amplitude correction. Electromechanical coupling efficiency is the efficiency coefficient of a piezoelectric actuator in converting electrical energy into mechanical energy, and it is frequency-dependent. The structural frequency response function is a frequency domain function that describes the transmission characteristics of a structure between the excitation point and the response point.

[0107] In the specific implementation, the phase channel component in the preliminary drive spectrum is extracted, and π is directly added to the phase component of all frequency points and all frequency bands to complete the phase reversal, so that the phase of the drive signal is 180 degrees out of phase with the original structural load, which meets the core requirement of anti-phase cancellation. Then, the amplitude channel component of the preliminary drive spectrum is extracted, and combined with the electromechanical coupling efficiency coefficient of the piezoelectric actuator at the corresponding frequency point and the amplitude of the structural frequency response function corresponding to the actuator's position, the amplitude component is corrected for each frequency point. The corrected amplitude is equal to the original amplitude multiplied by the reciprocal of the structural frequency response function amplitude, and then divided by the electromechanical coupling efficiency of the piezoelectric actuator to compensate for the electromechanical conversion loss of the actuator and the frequency dependence distortion of the structural transmission characteristics, ensuring that the compensated vibration amplitude generated by the drive at each frequency point is accurately matched with the original flow field excitation vibration amplitude. Finally, the corrected amplitude and the reversed phase are spliced ​​together to obtain the corrected frequency domain drive spectrum.

[0108] Step S606: Convert the corrected frequency domain drive spectrum into a time domain voltage signal through inverse fast Fourier transform, and introduce actuator saturation constraint and voltage limiting nonlinear correction to generate a preliminary time domain drive sequence.

[0109] It should be noted that the inverse fast Fourier transform is a mathematical operation that converts a frequency domain signal back to a time domain signal. Actuator saturation constraint is the maximum output constraint condition of a piezoelectric actuator due to material physical limits and drive power supply limitations. Voltage limiting nonlinear correction is a nonlinear processing method that applies soft or hard limiting to signals exceeding the actuator's voltage tolerance range.

[0110] It is worth noting that when there are components in the corrected frequency domain drive spectrum that exceed the rated voltage range of the piezoelectric actuator, if the inverse transformation is directly performed to obtain the time domain signal, the actuator will be unable to output the corresponding drive amplitude, or even damage the actuator or drive power supply. Therefore, saturation constraint correction must be introduced before generating the final drive sequence.

[0111] In the specific implementation, the voltage amplitude of the original time-domain voltage sequence obtained by the inverse fast Fourier transform is judged point by point to determine whether it exceeds the preset rated voltage range: for signal points that do not exceed the range, the original value is kept unchanged; for signal points that exceed the range, soft limiting is used, that is, the excess part is nonlinearly compressed proportionally to avoid high-frequency harmonic distortion caused by hard limiting, and finally a preliminary time-domain drive sequence that meets the physical constraints of the actuator is obtained.

[0112] Step S607: Perform actuator frequency response function compensation and group delay equalization on the preliminary time-domain drive sequence to obtain an anti-phase vibration compensation drive signal.

[0113] It should be noted that actuator frequency response function compensation is a pre-compensation for the influence of the actuator's own dynamic response characteristics on control signal distortion. Group delay equalization is a signal processing operation that corrects time alignment errors caused by differences in system group delays among different frequency components.

[0114] It is understandable that the electromechanical response of a piezoelectric actuator is not completely flat, and there are inherent differences in gain and phase response at different frequencies. These differences can cause the output compensation vibration to deviate from the expected design, so the drive signal needs to be frequency response corrected in advance.

[0115] In practice, the frequency response function of each piezoelectric actuator channel is obtained in advance through experimental measurement. After performing a Fourier transform on the initial time-domain drive sequence to the frequency domain, the frequency point is multiplied by the reciprocal of the amplitude of the actuator's frequency response function to compensate for the inherent phase delay of the phase component. Then, the inverse Fourier transform is used to transform it back to the time domain. Finally, the sliding window minimum mean square algorithm is used to equalize the group delay difference across the entire frequency band, eliminate the time offset of different frequency components, and ensure that the compensated vibration of each frequency component reaches the target control position simultaneously. Finally, a time-aligned, amplitude- and phase-accurate anti-phase vibration compensation drive signal is obtained.

[0116] Step S608: Apply the anti-phase vibration compensation drive signal to the piezoelectric actuator array of the thruster structure to generate a residual vibration response feedback sequence, and perform active vibration compensation control of the propulsion system based on the residual vibration response feedback sequence.

[0117] Understandably, the residual vibration response feedback sequence is the actual remaining vibration response signal of the structure after active control is applied, used to evaluate the control effect and guide adaptive parameter adjustment. The frequency band energy proportion is the ratio of the target frequency band energy to the total spectral energy, reflecting the relative intensity of the vibration in that frequency band. The local residual vibration weighted value is a local control effect evaluation index that comprehensively considers the amplitude and phase consistency of the neighboring vibrations. Closed-loop iterative optimization is an adaptive optimization process that continuously adjusts control parameters based on feedback information, gradually approaching the optimal control effect.

[0118] It is worth noting that active vibration compensation control cannot achieve optimal results with a single open-loop output. Affected by flow field load fluctuations, structural state changes, and actuator performance drift, local residual vibrations may still exist after a single control. Therefore, it is necessary to perform closed-loop iterative adjustments based on the residual response.

[0119] In practice, real-time vibration response signals are collected by piezoelectric sensors placed at key locations in the structure. After analog-to-digital conversion and detrending preprocessing, a residual vibration response feedback sequence is obtained. The feedback sequence is then subjected to frequency band energy calculation to extract the energy proportion of the target control frequency band and the local residual vibration weighted value of each key node. When the local residual vibration weighted value is higher than a preset threshold, the residual vibration response is re-input into the anti-phase control network to adjust the amplitude and phase parameters of the drive spectrum of each actuator channel. Multiple rounds of closed-loop iterative optimization are performed until the local residual vibration weighted value of all key nodes is lower than the threshold, ultimately achieving high-precision active compensation for the vibration of the propulsion system.

[0120] In one feasible implementation, step S608 may include: applying the anti-phase vibration compensation drive signal to the piezoelectric actuator array of the thruster structure after power amplification; acquiring the compensated time-domain vibration acceleration signal through an acceleration sensor array arranged on the surface of the thruster structure; performing bandpass filtering and double integration on the time-domain vibration acceleration signal to generate a residual vibration response feedback sequence containing acceleration, velocity, and displacement components; determining the frequency band energy ratio within the target frequency band based on the residual vibration response feedback sequence, and determining the local residual vibration weighting value of the key node neighborhood based on the root mean square value and phase dispersion of the vibration response within a preset neighborhood range of each key node in the residual vibration response feedback sequence; when the frequency band energy ratio exceeds a preset proportion, correcting the piezoelectric actuator drive spectrum gain coefficient and the output layer weight of the anti-phase control network according to the frequency band energy ratio, and correcting the local residual vibration weighting value of the key node neighborhood. The regularization intensity coefficient of the sparse identification algorithm is used to obtain the corrected driving spectrum gain coefficient, the corrected output layer weight, and the corrected regularization intensity coefficient. The corrected output layer weight is injected into the anti-phase control network to obtain the parameter-updated anti-phase control network. The corrected regularization intensity coefficient is injected into the sparse identification algorithm to drive the adaptive update of the sparse key node set to obtain the updated sparse key node set. Based on the updated sparse key node set and the parameter-updated anti-phase control network, the generation of the structural surface load spectrum and the generation of the anti-phase vibration compensation driving signal are re-executed, and the amplitude of the driving signal is adjusted based on the corrected driving spectrum gain coefficient during the generation process to obtain the updated anti-phase vibration compensation driving signal. The updated anti-phase vibration compensation driving signal is applied to the piezoelectric actuator array of the thruster structure to achieve closed-loop iterative optimization of the active vibration compensation control of the propulsion system.

[0121] It should be noted that power amplification is a signal conditioning process that amplifies the low-voltage signal output by the controller to the high voltage required by the piezoelectric actuator. The accelerometer array is a distributed measurement system composed of multiple piezoelectric or MEMS accelerometers, arranged at key vibration response locations on the thruster structure surface. In this embodiment, the power amplifier adopts a bipolar high-voltage output design, with an output voltage range of ±150V, a peak current of 3A, and a bandwidth of DC~20kHz, meeting the power requirements of wideband vibration control. The accelerometers are IEPE type piezoelectric accelerometers with a sensitivity of 100mV / g, a range of ±50g, and a frequency response of 0.5~10kHz, fixed to the structural surface by magnetic bases or adhesive bonding. The sensor array arrangement follows the modal shape coverage principle, ensuring that the responses of each major mode can be effectively observed.

[0122] Understandably, bandpass filtering is a digital filtering operation that removes DC drift, high-frequency electrical noise, and higher-order structural mode components from the measurement signal. Double integration is the process of successively converting the acceleration signal into velocity and displacement signals through frequency domain integration or time domain numerical integration.

[0123] The bandpass filter employs a zero-phase digital filter, such as a bidirectional filter implemented using the `filtfilt` function, to avoid phase distortion. The lower passband limit is 0.5Hz to eliminate DC drift and extremely low-frequency structural drift, while the upper passband limit is determined based on the target control frequency band and set at 500Hz, covering the range from the shaft frequency to five times the blade frequency. The quadratic integration uses a frequency domain integration method: after performing an FFT on the acceleration signal, it is divided by (iω). 2 The displacement spectrum is obtained and then IFFT is performed back to the time domain; the velocity signal is divided by (iω). The frequency domain integral suffers from noise amplification in the low-frequency band, which is compensated for by a high-frequency attenuation factor.

[0124] It is worth noting that the frequency band energy ratio is a global indicator for evaluating the anti-phase compensation effect, reflecting the concentration of residual vibration energy within the target frequency band in the total residual energy. The target frequency band is usually defined as the union of the propeller blade frequency BPF and its 2nd to 4th harmonics, i.e., [f BPF ,2f BPF ,3f BPF ,4f BPF The frequency bands are formed by taking ±5% bandwidth from each. The acceleration component a in the residual vibration response feedback sequence... filt Performing a Fast Fourier Transform on (t,s) yields the spectrum A(f,s). The frequency band energy percentage is then calculated. The calculation formula is: in, For the target frequency band set, Sampling frequency, Let f be the power spectral density of the s-th channel at frequency f.

[0125] ∈[0,1], the closer the value is to 0, the better the residual vibration suppression effect of the target frequency band is, and the larger the value is, the more it indicates insufficient compensation or resonance amplification.

[0126] The weighted value of local residual vibration in the neighborhood of key nodes is a spatial localization index for evaluating the compensation effect. It is used to identify which nodes in the sparsified key node set still have significant residual vibration in their neighborhoods, thus guiding the adjustment of the regularization intensity in subsequent sparse identification algorithms. For each node m in the sparsified key node set, its preset neighborhood range is defined as the distance from the node coordinate y to the local residual vibration. m Centered on, radius r n The spherical region N(m) = {s:||y = 0.05m~0.15m} s -ym ||2≤r n}, where y s Let be the spatial coordinates of the s-th accelerometer. If a sensor is located within the intersection of the neighborhoods of multiple nodes, it is assigned according to the nearest neighbor principle. The root mean square value of the vibration response within the neighborhood and the phase dispersion reflect the overall intensity and phase consistency of the residual vibration within the neighborhood, respectively. The weighted fusion of these two values ​​yields the weighted value of the local residual vibration corresponding to that node. Specifically, by extracting the vibration spectral components of the target frequency band from all sensor channels within the neighborhood N(m), and calculating the square root of the squared average of the vibration amplitude at each frequency point, the root mean square value σ of the vibration response of the target frequency band corresponding to that neighborhood is obtained. rms (m); then, for the vibration phase components of all sensors at the same frequency point within the neighborhood, calculate the standard deviation of the phase samples to obtain the phase dispersion at that frequency point, and then take the average of the phase dispersion of all target frequency points to obtain the average phase dispersion σ of the neighborhood. phase (m); final local residual vibration weighted value w res ( m)=ασ rms (m)+βσ phase (m), where α and β are pre-set weighting coefficients that satisfy α+β=1 and can be adjusted according to control requirements. Usually, α=0.6 and β=0.4 are taken.

[0127] The preset ratio is usually set to 0.1~0.2. That is, when the residual energy in the target frequency band accounts for more than 10%~20% of the total energy, the current compensation effect is deemed insufficient, and the closed-loop correction mechanism is triggered. This threshold needs to be set according to the type of thruster and the stealth index. For the thrusters of underwater vehicles with high stealth requirements, a lower value, such as 0.08, is advisable to ensure a more stringent vibration suppression level.

[0128] The drive spectrum gain coefficient is a global amplitude adjustment factor used to quickly compensate for drive signal amplitude mismatch caused by abrupt changes in flow field state or actuator efficiency drift. It is adaptively adjusted based on the frequency band energy ratio: the higher the frequency band energy ratio, the larger the corrected drive spectrum gain coefficient, achieving a rapid increase in global compensation intensity. The output layer weights of the anti-phase control network are adjusted node by node according to the weighted value of local residual vibration. The greater the residual vibration, the greater the adjustment range of the corresponding channel weights, ensuring that higher compensation intensity is obtained at local strong residual locations.

[0129] The regularization intensity coefficient is positively correlated with the maximum local residual vibration weight value of all key nodes: when the local residual vibration weight values ​​are generally large, increasing the regularization intensity prompts the sparse identification algorithm to select fewer but more critical control nodes, avoiding over-control; when there are only a few local high residual nodes, decreasing the regularization intensity allows the sparse identification algorithm to expand the number of nodes and specifically enhance the local compensation capability.

[0130] It is worth noting that the corrected output layer weights are sent to the FPGA implementation module of the anti-phase control network via the control bus. In the FPGA, the update of the output layer weights employs a ping-pong buffer mechanism: new weights are first written to a spare weight storage area, and after the current control cycle ends, the new weights are activated via hardware signal switching, ensuring that the weight update process does not interrupt the drive signal output. The updated anti-phase control network performs frequency domain amplitude inversion with the corrected output layer weights in the next control cycle, generating an updated preliminary drive spectrum. Simultaneously, the corrected drive spectrum gain coefficients are applied to the amplitude correction stage, optimizing both the global amplitude and local mapping relationship of the drive signal.

[0131] After the modified regularization intensity coefficient is injected into the sparse identification algorithm, it triggers a re-solution of the Lasso sparse regression model. Due to real-time requirements, instead of starting the coordinate descent iteration from scratch, it uses the current sparse weights as the initial values ​​for a warm start, performing only a small number of coordinate descent iterations, such as 5-10, to allow the weights to quickly converge to the optimal solution under the new regularization intensity. After convergence, the node indices corresponding to non-zero weights are extracted to generate an updated set of sparse key nodes. Simultaneously, the vibration contribution ranking of each node is updated for use in calculating the local residual vibration weighted values ​​for the next control cycle.

[0132] It is worth noting that, based on the updated sparse key node set and the parameter-updated anti-phase control network, the sparse RBF load mapping is re-executed to obtain the updated structural surface load spectrum. Then, the updated drive spectrum is generated through frequency band decomposition and anti-phase projection. After the amplitude is corrected by the drive spectrum gain coefficient, the actuator frequency response function compensation and group delay equalization are completed sequentially, finally obtaining the updated anti-phase vibration compensation drive signal. The updated drive signal is amplified and applied to the piezoelectric actuator array to achieve a new round of compensation effect feedback adjustment. Through multiple iterations, the residual vibration is gradually suppressed to the design requirement range, completing the high-precision active compensation control of the propulsion system vibration.

[0133] In this embodiment, amplitude normalization and phase unwinding processes eliminate control ambiguities caused by differences in spectral amplitude magnitudes and phase jumps. Distance-direction dual-factor modeling of spatial coupling strength achieves precise spatial matching between actuator driving energy and the dominant vibration region. Row normalization of the actuator-node association weight matrix establishes a weighted mapping channel from the node load space to the actuator driving space. Nonlinear characteristic transformation of the fully connected mapping layer enables a high-dimensional abstract mapping from the frequency domain conditional tensor to the preliminary driving spectrum. Phase reversal and amplitude matching correction of the frequency domain amplitude inversion layer generate a frequency domain driving spectrum that is phase-opposite to the structural load and compensated for by electromechanical efficiency and structural transmission loss. Soft saturation voltage limiting and actuator frequency response function compensation ensure accurate restoration of the time-domain driving signal within physical constraints. Group delay equalization eliminates waveform distortion caused by time delay differences in different frequency band components. Closed-loop iterative correction driven by residual vibration feedback achieves adaptive dynamic optimization of the driving spectrum gain coefficient, network weights, and sparse parameters.

[0134] It should be noted that the above examples are only for understanding this application and do not constitute a limitation on the active vibration compensation control method for propulsion systems based on RBF load mapping in this application. Any simple modifications based on this technical concept are within the protection scope of this application.

[0135] This application also provides a propulsion system vibration active compensation control device based on RBF load mapping, please refer to... Figure 2 The active vibration compensation control device for propulsion systems based on RBF load mapping includes: Module 10 is used to collect time-domain data of pulsating pressure in the unsteady flow field around the thruster, and constructs a time-domain matrix of pulsating pressure in the flow field by extracting flow field features and aligning them with time and space.

[0136] The construction module 10 is also used to construct a fluid-structure interaction mapping kernel function matrix based on the Wendland C2 type radial basis function of the locally curvature adaptive compactly supported domain, combined with the geometric nonlinear characteristics of the fluid-structure interaction interface.

[0137] The mapping module 20 is used to input the time-domain matrix of the flow field pulsating pressure and the fluid-structure interaction mapping kernel function matrix into the RBF neural network mapping model to obtain the time-domain sequence of the nodal load on the wet surface of the propeller.

[0138] The filtering module 30 is used to process the time-domain sequence of the wet surface node load through a sparse identification algorithm, and to filter the key nodes that contribute the most to the vibration based on Lasso regularization constraints, so as to obtain a sparse key node set and a sparse load time-domain sequence.

[0139] Transformation module 40 is used to perform frequency domain transformation on the sparse load time domain sequence and generate a structural surface load spectrum matrix by combining the spatial distribution characteristics of the sparsified key node set.

[0140] The control module 50 is used to input the surface load spectrum matrix of the structure into the anti-phase control network of the piezoelectric actuator array, generate an anti-phase vibration compensation drive signal, and apply the anti-phase vibration compensation drive signal to the piezoelectric actuator array of the thruster structure to realize active vibration compensation control of the propulsion system.

[0141] The RBF load mapping-based active vibration compensation control device for propulsion systems provided in this application employs the RBF load mapping-based active vibration compensation control method for propulsion systems described in the above embodiments. It solves the technical problems of existing propulsion system vibration control relying on lag-type PID strategies, inaccurate fluid-structure interaction load transfer, and excessive computational load, thus failing to achieve real-time active suppression of high-frequency vibrations caused by unsteady flow fields. Compared with the prior art, the beneficial effects of the RBF load mapping-based active vibration compensation control device for propulsion systems provided in this application are the same as those of the RBF load mapping-based active vibration compensation control method for propulsion systems provided in the above embodiments. Furthermore, other technical features of the RBF load mapping-based active vibration compensation control device for propulsion systems are the same as those disclosed in the methods of the above embodiments, and will not be repeated here.

[0142] This application also provides a non-transitory computer-readable storage medium storing a computer program thereon, which, when executed by a processor, implements the active vibration compensation control method for propulsion systems based on RBF load mapping as described above.

[0143] The above are only some embodiments of this application and do not limit the patent scope of this application. All equivalent structural transformations made under the technical concept of this application and using the contents of the specification and drawings of this application, or direct / indirect applications in other related technical fields, are included in the patent protection scope of this application.

Claims

1. A method for active vibration compensation control of a propulsion system based on RBF load mapping, characterized in that, The method includes: The time-domain data of pulsating pressure in the unsteady flow field around the thruster are collected, and the flow field pulsating pressure time-domain matrix is ​​constructed by extracting flow field features and aligning them with time and space. Based on the Wendland C2 type radial basis function with local curvature adaptive compactly supported domain, and combined with the geometric nonlinear characteristics of the fluid-structure interaction interface, a fluid-structure interaction mapping kernel function matrix is ​​constructed. The time-domain matrix of the flow field pulsating pressure and the fluid-structure interaction mapping kernel function matrix are input into the RBF neural network mapping model to obtain the time-domain sequence of the nodal load on the wet surface of the thruster. The time-domain sequence of the load at the wet surface nodes is processed by a sparse identification algorithm. Based on Lasso regularization constraints, the key nodes that contribute the most to the vibration are selected to obtain a sparse key node set and a sparse load time-domain sequence. The sparse load time-domain sequence is transformed in the frequency domain, and combined with the spatial distribution characteristics of the sparsified key node set, a structural surface load spectrum matrix is ​​generated. The surface load spectrum matrix of the structure is input into the anti-phase control network of the piezoelectric actuator array to generate an anti-phase vibration compensation drive signal. The anti-phase vibration compensation drive signal is then applied to the piezoelectric actuator array of the thruster structure to achieve active vibration compensation control of the propulsion system.

2. The method as described in claim 1, characterized in that, The time-domain data of the fluctuating pressure in the unsteady flow field around the thruster are collected. Through flow field feature extraction and spatiotemporal alignment, a flow field fluctuating pressure time-domain matrix is ​​constructed, including: A cluster of multiphysics sensors is deployed in a circumferential array and an axial gradient in the flow field around the thruster to collect time-domain data of pulsating pressure in the unsteady flow field around the thruster. The local state discrimination factor of the flow field is determined based on the turbulent pulsating velocity and fluid temperature, and the abnormal data of the pulsating pressure time domain data is removed based on the local state discrimination factor of the flow field to obtain the cleaned pulsating pressure time domain data. Multi-scale wavelet packet decomposition is performed on the cleaned pulsating pressure time-domain data to obtain the dominant pressure pulsation component within a preset frequency band. Based on the propeller shaft frequency encoded signal, the rotational cycle of the dominant component of the pressure pulsation is synchronized to generate a phase-locked pressure sequence that is phase-locked with the blade rotation. The phase-locked pressure sequence is normalized in time and space to construct a standardized pulsating pressure time-domain signal matrix; Singular value decomposition is performed on the standardized pulsating pressure time-domain signal matrix to extract the first K dominant mode pressure components and generate a compressed dominant mode pressure time-domain matrix, where K is a positive integer; Based on the wet surface node topology numbering of the finite element model of the thruster structure, a spatiotemporal mapping index between the time domain matrix of the dominant mode pressure after compression and the structural nodes is established. The time-domain matrix of flow field pulsation pressure corresponding to the wet surface node of the structure is determined based on the spatiotemporal mapping index.

3. The method as described in claim 1, characterized in that, The Wendland C2 type radial basis function based on the locally curvature adaptive compactly supported domain, combined with the geometric nonlinear characteristics of the fluid-structure interaction interface, constructs a fluid-structure interaction mapping kernel function matrix, including: Extract the set of center coordinates of wet surface elements and the set of element normal vectors from the finite element model of the thruster structure, and obtain the set of node coordinates in the fluid domain; The three-dimensional Euclidean distance matrix between the center of the wet surface unit and the fluid domain node is determined based on the set of coordinates of the wet surface unit center and the set of coordinates of the fluid domain node. The local Gaussian curvature distribution and the average curvature distribution of the wet surface are determined by local surface fitting based on the set of unit normal vectors. Based on the local Gaussian curvature distribution and the average curvature distribution of the wet surface, an adaptive compact support domain adjustment factor for curvature is determined, and the three-dimensional Euclidean distance matrix is ​​corrected to an equivalent distance matrix along the geodesic of the wet surface based on the adaptive compact support domain adjustment factor for curvature. The adaptive compact support domain adjustment factor for curvature characterizes the nonlinear tortuous effect of the local curvature of the fluid-structure interaction interface on the load transfer path. Based on the turbulent kinetic energy density distribution and local vorticity distribution at the fluid domain nodes, the local strength index of fluid-structure interaction is determined, and the support domain radius distribution of the Wendland C2 type radial basis function of the locally curvature adaptive compactly supported domain is adjusted according to the local strength index of fluid-structure interaction to obtain the adjusted support domain radius distribution. Based on the equivalent distance matrix and the adjusted support domain radius distribution, a radial basis function kernel matrix is ​​constructed, and a system of linear equations with weight coefficients at the center of the wet surface unit as the constraint point is established. The linear equations for the weighting coefficients are solved by constrained least squares method to obtain the mapping weighting coefficient vector from the fluid domain node to the center of the wetted surface unit. The mapping weight coefficient vector is fused with the radial basis function kernel matrix to obtain the fluid-structure interaction mapping kernel function matrix.

4. The method as described in claim 1, characterized in that, The step of inputting the time-domain matrix of the flow field pulsating pressure and the fluid-structure interaction mapping kernel function matrix into the RBF neural network mapping model to obtain the time-domain sequence of the nodal load on the wetted surface of the thruster includes: The flow field pulsating pressure time domain matrix is ​​subjected to overlapping windowing processing according to a preset time window length to generate multiple short time sequence segments; Using the fluid-structure interaction mapping kernel function matrix as the initial hidden layer weights, the network topology of the RBF neural network mapping model is initialized to obtain the initialized RBF neural network mapping model. The short time sequence segment and the fluid-structure interaction mapping kernel function matrix are input in parallel into the initialized RBF neural network mapping model. The weights of the network output layer are updated in real time by the online recursive least squares algorithm to obtain the updated network output layer weights. Based on the updated network output layer weights and the fluid-structure interaction mapping kernel function matrix, radial basis interpolation is performed on the fluid node pulsating pressure in each time window to determine the dynamic pressure distribution vector at the center of the wet surface unit. Based on the dynamic pressure distribution vector and the area properties of the corresponding wetted surface unit, the concentrated force vector of the surface unit is obtained. The concentrated force vector of the surface element is weighted and distributed to all structural nodes belonging to the corresponding wet surface element according to the shape function or area coordinate of the corresponding wet surface element, so as to obtain the equivalent dynamic nodal force on each structural node. The equivalent dynamic nodal forces of all structural nodes corresponding to the wet surface units are aggregated and stacked in time sequence to generate the time-domain sequence of wet surface nodal loads of the thruster.

5. The method as described in claim 1, characterized in that, The process involves processing the time-domain sequence of nodal loads on the wet surface using a sparse identification algorithm, and selecting the key nodes with the greatest contribution to vibration based on Lasso regularization constraints. This yields a sparse key node set and a sparse load time-domain sequence, including: Extract the first N main wet surface mode shape vectors from the finite element model of the thruster structure, where N is a positive integer; Based on the main wet surface mode shape vectors, the modal participation factor vectors of each wet surface node are determined, and a spatial feature matrix weighted by the modal participation factors is constructed. Determine the time-domain energy density of each nodal load component in the time-domain sequence of the wet surface nodal load, and construct a regression target vector with time-domain energy density as the response variable; Using the spatial feature matrix as the prediction variable and the regression target vector as the response variable, a Lasso sparse regression model with L1 regularization constraint is established, and the initial regularization strength coefficient is set. The Lasso sparse regression model is iteratively optimized using the coordinate descent method. During the iteration process, a soft threshold operation is applied to the mapping weights. After the iteration converges, the node indices corresponding to the non-zero elements of the mapping weights are extracted to generate a sparsified key node set. The vibration contribution ranking of each node in the sparsified key node set is determined, and based on the vibration contribution ranking and the sparsified key node set, the load components of the corresponding key nodes are extracted from the wet surface node load time-domain sequence to obtain the sparse load time-domain sequence.

6. The method as described in claim 1, characterized in that, The step of performing a frequency domain transformation on the sparse load time-domain sequence, combined with the spatial distribution characteristics of the sparsified key node set, to generate a structural surface load spectrum matrix includes: Windowed short-time Fourier transform is performed on each key node component in the sparse load time-domain sequence to obtain the complex spectrum sequence of each key node at different times. The amplitude spectrum component and phase spectrum component are extracted from the complex spectrum sequence, and the frequency band is divided based on the thruster excitation characteristic frequency to generate the spectral energy density distribution in each frequency band. Based on the three-dimensional spatial coordinates of each node in the sparse key node set, a frequency domain-space joint index is constructed. The index key of the frequency domain-space joint index is a spatial location code generated based on the three-dimensional spatial coordinates of the node, and the index value is the spectral identifier of the corresponding node. Bind the amplitude spectrum component and the phase spectrum component to the corresponding frequency domain-space joint index to establish a frequency domain-space joint index table; Based on the spectral energy density distribution, the dominant excitation frequency and its corresponding harmonic components are identified, and the spectral components of each key node at the dominant excitation frequency are topologically aggregated to obtain the topological aggregation result. Based on the topological aggregation results and the frequency domain-space joint index table, a frequency domain-space joint representation tensor is constructed. The frequency-space joint characterization tensor is reorganized according to the node-frequency point-frequency band dimension to obtain the structural surface load spectrum matrix.

7. The method as described in claim 1, characterized in that, The step of inputting the surface load spectrum matrix of the structure into the anti-phase control network of the piezoelectric actuator array to generate an anti-phase vibration compensation drive signal, and applying the anti-phase vibration compensation drive signal to the piezoelectric actuator array of the thruster structure to achieve active vibration compensation control of the propulsion system includes: The load spectrum matrix of the structure surface is normalized in magnitude, and the phase spectrum of each key node component is unwrapped channel by channel, and encoded into a frequency domain conditional tensor containing both magnitude and phase channels. Obtain the spatial layout coordinates and installation orientation angles of each actuator channel in the piezoelectric actuator array on the surface of the thruster structure, and determine the spatial coupling strength between each actuator and each node in the sparse key node set. Based on the spatial coupling strength, an actuator-node association weight matrix is ​​constructed. The dimension of the actuator-node association weight matrix is ​​determined according to the number of channels and key nodes of each actuator. The frequency domain conditional tensor and the actuator-node association weight matrix are input into the inverse phase control network, and the node frequency domain features are mapped to the preliminary driving spectrum of each actuator channel through a fully connected mapping layer. In the frequency domain amplitude inversion layer of the anti-phase control network, a phase inversion operation is performed on the preliminary driving spectrum, and amplitude matching correction is performed based on the electromechanical coupling efficiency and structural frequency response function of the piezoelectric actuator to obtain the corrected frequency domain driving spectrum. The corrected frequency domain drive spectrum is converted into a time domain voltage signal by inverse fast Fourier transform, and actuator saturation constraint and voltage limiting nonlinear correction are introduced to generate a preliminary time domain drive sequence. The initial time-domain drive sequence is subjected to actuator frequency response function compensation and group delay equalization to obtain an anti-phase vibration compensation drive signal. The anti-phase vibration compensation drive signal is applied to the piezoelectric actuator array of the thruster structure to generate a residual vibration response feedback sequence, and the active vibration compensation control of the propulsion system is performed based on the residual vibration response feedback sequence.

8. The method as described in claim 7, characterized in that, The step of applying the anti-phase vibration compensation drive signal to the piezoelectric actuator array of the thruster structure to generate a residual vibration response feedback sequence, and performing active vibration compensation control of the propulsion system based on the residual vibration response feedback sequence, includes: The anti-phase vibration compensation drive signal is amplified and applied to the piezoelectric actuator array of the thruster structure. The compensated time-domain vibration acceleration signal is collected by an acceleration sensor array arranged on the surface of the thruster structure. The time-domain vibration acceleration signal is subjected to bandpass filtering and double integration to generate a residual vibration response feedback sequence containing acceleration, velocity and displacement components. The frequency band energy ratio within the target frequency band is determined based on the residual vibration response feedback sequence, and the local residual vibration weighting value of the neighborhood of the key node is determined based on the root mean square value and phase dispersion of the vibration response within the preset neighborhood range of each key node in the residual vibration response feedback sequence. When the frequency band energy ratio exceeds a preset ratio, the piezoelectric actuator drive spectrum gain coefficient and the output layer weight of the anti-phase control network are corrected according to the frequency band energy ratio, and the regularization intensity coefficient of the sparse identification algorithm is corrected according to the local residual vibration weight value of the neighborhood of the key node, so as to obtain the corrected drive spectrum gain coefficient, the corrected output layer weight and the corrected regularization intensity coefficient. The corrected output layer weights are injected into the anti-phase control network to obtain the parameter-updated anti-phase control network. The modified regularization strength coefficient is injected into the sparse identification algorithm to drive the sparse key node set to be updated adaptively, thus obtaining the updated sparse key node set. Based on the updated sparse key node set and the updated anti-phase control network, the generation of the structural surface load spectrum and the generation of the anti-phase vibration compensation drive signal are re-executed. During the generation process, the amplitude of the drive signal is adjusted based on the corrected drive spectrum gain coefficient to obtain the updated anti-phase vibration compensation drive signal. The updated anti-phase vibration compensation drive signal is applied to the piezoelectric actuator array of the thruster structure to achieve closed-loop iterative optimization of the active vibration compensation control of the propulsion system.

9. A propulsion system vibration active compensation control device based on RBF load mapping, characterized in that, The device includes: The module is used to collect the time-domain data of the pulsating pressure of the unsteady flow field around the thruster, and construct the time-domain matrix of the pulsating pressure of the flow field by extracting flow field features and aligning them with time and space. The construction module is also used to construct a fluid-structure interaction mapping kernel function matrix based on the Wendland C2 type radial basis function of the locally curvature adaptive compactly supported domain, combined with the geometric nonlinear characteristics of the fluid-structure interaction interface. The mapping module is used to input the time-domain matrix of the flow field pulsating pressure and the fluid-structure interaction mapping kernel function matrix into the RBF neural network mapping model to obtain the time-domain sequence of the nodal load on the wet surface of the propeller. The filtering module is used to process the time-domain sequence of the wet surface node load through a sparse identification algorithm, and to filter the key nodes that contribute the most to the vibration based on Lasso regularization constraints, so as to obtain a sparse key node set and a sparse load time-domain sequence. The transformation module is used to perform frequency domain transformation on the sparse load time domain sequence and generate a structural surface load spectrum matrix by combining the spatial distribution characteristics of the sparsified key node set. The control module is used to input the surface load spectrum matrix of the structure into the anti-phase control network of the piezoelectric actuator array, generate an anti-phase vibration compensation drive signal, and apply the anti-phase vibration compensation drive signal to the piezoelectric actuator array of the thruster structure to realize active vibration compensation control of the propulsion system.

10. A non-transitory computer-readable storage medium having a computer program stored thereon, characterized in that, When the computer program is executed by the processor, it implements the active vibration compensation control method for propulsion systems based on RBF load mapping as described in any one of claims 1 to 8.

Images