Modelling method for underwater propeller ring propeller based on simulation

Abstract

The application relates to the technical field of underwater propeller modeling, and discloses a modeling method for an annular propeller of an underwater propeller based on simulation emulation. The method is based on three-dimensional geometric parameters and material attribute data of the annular propeller, and an initial propeller blade fluid domain grid model is constructed by using a discrete element simulation algorithm; after a multiphase flow coupling boundary condition is applied to the initial model, fluid dynamic load distribution data are obtained through a transient flow field solver; propeller blade deformation characteristic values are generated according to the load data, and if the deformation characteristic values exceed a preset tolerance interval, a dynamic grid reconstruction instruction is activated; based on the instruction, fluid-structure interaction interface pressure pulsation data are extracted, a stress concentration coefficient is calculated in combination with material attributes, and an abnormal signal is generated if the stress concentration coefficient exceeds a safety threshold; when the signal is triggered, an eddy current field tracking module is called to update the fluid domain grid model.

Description

Technical Field

[0001] This invention relates to the field of underwater propeller modeling technology, specifically to a simulation-based method for modeling annular underwater propellers. Background Technology

[0002] In the field of underwater equipment development, propellers, as core power components, directly affect the equipment's navigation efficiency, stability, and energy consumption. Annular propellers, due to their characteristics of low water flow disturbance, high propulsion efficiency, and low noise, are widely used in underwater robots, deep-sea exploration equipment, and submarines. However, the complexity of the underwater environment presents numerous challenges to the design and performance analysis of annular propellers. Accurate modeling and simulation technologies are crucial prerequisites for achieving optimized propeller design.

[0003] Traditional propeller modeling methods are mostly based on empirical formulas or static mesh simulations, which fail to accurately reflect the dynamic interaction between the complex underwater flow field and the propeller structure. During the modeling process, the initial mesh model often relies on simplified geometric parameters, neglecting the influence of material properties on structural deformation, leading to significant deviations between simulation results and actual operating conditions. Furthermore, the coupling effects of multiphase flows in the underwater environment generate complex hydrodynamic loads. Traditional methods handle these boundary conditions rather coarsely, only acquiring local flow field data and failing to comprehensively capture the dynamic load distribution during propeller operation.

[0004] In existing simulation techniques, the mesh model remains static once it is determined, making it difficult to address the deformation of propellers under hydrodynamic forces. As the blades deform, the deviation between the static mesh and the actual structure gradually increases, leading to distortion in subsequent flow field solutions and stress analysis results. Furthermore, traditional methods lack an effective feedback mechanism for monitoring the stress state of the blades, making it impossible to adjust simulation parameters in real time based on stress concentration, thus hindering the early identification of potential structural failure risks during the propeller design phase.

[0005] In handling eddy current fields, traditional modeling methods often employ mesh refinement strategies that pre-define fixed regions, failing to dynamically adjust the refinement range and accuracy based on changes in the flow field. During underwater annular propeller operation, the eddy current field distribution around the blades is complex and dynamically changes with operating conditions. Fixed mesh refinement methods either result in insufficient refinement leading to the loss of eddy current details or excessive refinement increasing the computational load and reducing simulation efficiency. These issues collectively result in significant limitations of traditional modeling methods in the dynamic performance analysis and structural reliability assessment of annular propellers, making it difficult to meet the development requirements of high-precision, high-reliability underwater propulsion systems. Summary of the Invention

[0006] The purpose of this invention is to provide a simulation-based modeling method for annular propellers of underwater thrusters to solve the problems mentioned in the background art.

[0007] To achieve the above objectives, this invention provides a simulation-based method for modeling annular propellers of underwater thrusters, the method comprising:

[0008] Based on the three-dimensional geometric parameters and material property data of the annular propeller of the underwater thruster, an initial blade fluid domain mesh model was constructed using the discrete element simulation algorithm.

[0009] After applying multiphase flow coupled boundary conditions to the initial blade fluid domain mesh model, the hydrodynamic load distribution data of the annular propeller during operation is obtained through a transient flow field solver.

[0010] Based on the hydrodynamic load distribution data, blade deformation characteristic values ​​are generated. If the blade deformation characteristic values ​​exceed the preset deformation tolerance range, a dynamic mesh reconstruction command is activated.

[0011] Based on the dynamic mesh reconstruction command, pressure pulsation data of the fluid-structure interaction interface is extracted, and the stress concentration factor of the blade structure is calculated in combination with the material property data. If the stress concentration factor exceeds the preset stress concentration factor safety threshold set based on the material property data, an abnormal stress state signal of the blade is generated.

[0012] When the abnormal stress state signal of the blade is triggered, the eddy current field tracking module is invoked to update the fluid domain mesh model.

[0013] Preferably, the three-dimensional geometric parameters include the blade tilt angle distribution curve, hub diameter, and blade tip clearance of the annular propeller, and the location data of fluid separation points on the blade surface are acquired by particle image velocimetry technology.

[0014] The blade tilt angle distribution curve is discretized and the rate of change of tilt angle with radial position is calculated to obtain the blade angle of attack gradient parameter;

[0015] By mapping the blade tilt angle distribution curve with the fluid separation point location data in spatial coordinates, a reference value for the blade aerodynamic efficiency is obtained.

[0016] The blade angle of attack gradient parameters of the initial blade fluid domain grid model are corrected based on the aerodynamic efficiency reference value.

[0017] Preferably, the process of constructing the multiphase flow coupled boundary conditions includes:

[0018] The characteristic frequencies of cavitation effect within the rotation cycle of the annular propeller are extracted, and the characteristic frequencies of cavitation effect are superimposed with the natural frequencies of the propeller structure for spectral analysis.

[0019] If the result of the spectral superposition analysis exceeds the resonance threshold range, then the eddy current field suppression control signal is output.

[0020] The turbulent viscosity coefficient of the transient flow field solver is adjusted based on the eddy field suppression control signal.

[0021] Preferably, the method for calculating the blade deformation characteristic value is as follows:

[0022] The blade structure is divided into several deformation monitoring sub-regions, and the real-time load fluctuation amplitude of each deformation monitoring sub-region is obtained;

[0023] The deviation between the real-time load fluctuation amplitude and the material yield limit value of the corresponding sub-region is calculated to obtain the local deformation risk coefficient.

[0024] The blade deformation characteristic value is output by performing a weighted summation operation on all local deformation risk coefficients.

[0025] Preferably, the execution process of the dynamic mesh reconstruction instruction includes:

[0026] The coordinates of high-risk deformation monitoring sub-regions are located based on abnormal blade stress state signals.

[0027] Acquire fluid velocity gradient data and pressure pulsation peak data in the high-risk deformation monitoring sub-region;

[0028] Fluid velocity gradient data and pressure pulsation peak data are input into an adaptive mesh optimization algorithm to generate local mesh refinement level parameters.

[0029] The node distribution density of the fluid domain mesh model is updated based on the local mesh refinement level parameters.

[0030] Preferably, the operation flow of the eddy current field tracking module includes:

[0031] The vortex coreline was extracted based on the updated fluid domain grid model to obtain the spatial evolution trajectory of the blade wake vortex structure.

[0032] The spatial evolution trajectory is matched with the preset eddy current decay model. If the similarity matching result is lower than the eddy current stability threshold, the tip eddy correction factor is output.

[0033] Adjust the eddy viscosity equation coefficients of the transient flow field solver based on the tip eddy correction factor.

[0034] Preferably, the weighted summation method for the local deformation risk coefficient is as follows:

[0035] Obtain the cumulative material fatigue damage value for each deformation monitoring sub-region and label it as the sub-region damage weight coefficient;

[0036] The local deformation risk coefficient is multiplied by the corresponding sub-region damage weight coefficient to obtain the weighted deformation risk value.

[0037] The arithmetic mean of all weighted deformation risk values ​​is processed to output the overall deformation characteristic value of the blade.

[0038] Preferably, the execution steps of the adaptive mesh optimization algorithm include:

[0039] Construct a correlation matrix between the location data of fluid separation points on the blade surface and the peak pressure pulsation data;

[0040] The priority ranking of grid encryption regions is determined based on the distribution of eigenvalues ​​in the correlation matrix.

[0041] The boundary layer thickness parameters of the fluid domain mesh model are coupled and iteratively calculated with the blade tip clearance value according to priority to generate local mesh refinement level parameters.

[0042] Preferably, when the eddy field tracking module generates the blade tip vortex correction factor during execution, it simultaneously collects the turbulent kinetic energy attenuation rate of the blade tail vortex structure; and performs a difference calculation between the turbulent kinetic energy attenuation rate and the eddy stability threshold to obtain the vortex core dissipation compensation coefficient.

[0043] The iteration step size parameter of the vortex viscosity equation coefficients of the transient flow field solver is corrected based on the vortex core dissipation compensation coefficient.

[0044] Preferably, the method for obtaining the sub-region damage weight coefficient includes:

[0045] Extract the stress amplitude spectral density distribution data for each deformation monitoring sub-region, and calculate the equivalent number of alternating stress cycles based on the stress amplitude spectral density distribution data;

[0046] The equivalent alternating stress cycle number is matched with the material SN curve, and the sub-region damage weight coefficient update command is output.

[0047] Compared with the prior art, the beneficial effects of the present invention are:

[0048] This simulation-based modeling method for underwater propeller rings significantly improves the accuracy and adaptability of the modeling process through multi-stage collaborative optimization. In the initial mesh model construction stage, relying on 3D geometric parameters and material property data, combined with discrete element method (DEM) simulation algorithms, it can more realistically reproduce the geometric features and material properties of the propeller blades, laying a practical foundation for subsequent simulations. This construction method avoids the distortion of geometric and material properties in traditional simplified modeling, making the initial model closer to the physical essence of the propeller.

[0049] The application of multiphase flow coupled boundary conditions makes the simulation environment closer to the complex actual underwater working conditions. The hydrodynamic load distribution data obtained by the transient flow field solver can comprehensively reflect the dynamic load on the blades under different navigation conditions, breaking the limitation of traditional methods that can only capture local flow field information, and providing more comprehensive basic data for subsequent deformation analysis and stress assessment.

[0050] Monitoring deformation eigenvalues ​​and activating dynamic mesh reconstruction commands form an effective dynamic adjustment mechanism. When the blade deformation exceeds the preset range, dynamic mesh reconstruction can promptly correct the deviation between the mesh model and the actual structure, avoiding the simulation accuracy degradation caused by the accumulation of deformation in the static mesh. This dynamic response mechanism allows the simulation process to follow the structural changes of the blade in real time, ensuring that the flow field analysis and structural state always maintain a correspondence, thus improving the dynamic fit of the simulation.

[0051] In the fluid-structure interaction analysis phase, by extracting pressure pulsation data and calculating the stress concentration factor in conjunction with material properties, the structural stress state of the blade under complex flow fields can be accurately captured. The generation of anomaly signals in the stress state provides a direct basis for identifying potential structural risks, enabling the modeling process to not only reflect flow field characteristics but also delve into the mechanical response of the associated structure, thus achieving a collaborative analysis of fluid dynamics and structural performance.

[0052] The combination of the eddy current field tracking module and local mesh refinement parameter updates optimizes the efficiency of mesh resource allocation. Dynamically adjusting local mesh refinement parameters based on stress anomaly signals improves mesh accuracy in complex eddy current regions to capture detailed flow fields, while maintaining a reasonable mesh density in stable regions to control computational load, achieving a balance between simulation accuracy and efficiency. This adaptive mesh optimization approach avoids the problems of insufficient accuracy or low efficiency in traditional fixed mesh refinement strategies, ensuring that simulations remain highly efficient and accurate in complex underwater environments.

[0053] This method constructs a simulation system that is closer to the actual operating state through the whole process of geometric modeling, boundary conditions, dynamic adjustment, stress analysis and mesh optimization. It can more comprehensively and accurately reflect the dynamic characteristics and structural response of the annular propeller in the underwater environment, and meet the needs of underwater propulsion research and development for high-precision modeling. Attached Figure Description

[0054] Figure 1 This is a schematic diagram illustrating the working principle of the underwater propeller modeling method based on simulation described in this invention.

[0055] Figure 2 A flowchart for correcting the three-dimensional geometric parameters and mesh model of a ring propeller;

[0056] Figure 3A flowchart for calculating the characteristic values ​​of blade deformation;

[0057] Figure 4 This is a flowchart illustrating the operation of the eddy current field tracking module. Detailed Implementation

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

[0059] Please see Figure 1 This invention provides a simulation-based method for modeling annular propellers of underwater thrusters, the method comprising:

[0060] Based on the three-dimensional geometric parameters and material property data of the underwater propeller's annular propeller, an initial blade fluid domain mesh model was constructed using a discrete element method (DEM) simulation algorithm. The three-dimensional geometric parameters include the blade tilt angle distribution curve, hub diameter, and tip clearance. Material property data includes parameters such as elastic modulus, density, and yield strength. The DEM simulation algorithm employs a particle-based discretization method, dividing the blade surface into discrete elements to generate a fluid domain mesh containing nodes and surface elements. The initial mesh resolution is preset to a uniform distribution based on the propeller size and operating conditions. After the initial blade fluid domain mesh model is generated, multiphase flow coupling boundary conditions are applied. These boundary conditions include an interfacial interaction model of the fluid phases (such as water, air, or cavitation), describing the fluid interaction between the liquid-gas interface and the blade surface. The transient flow field solver uses a numerical method based on the Navier-Stokes equations to calculate the hydrodynamic load distribution data under the boundary conditions. This data includes field quantities such as pressure, velocity, and shear force, covering all mesh points on the blade surface. Based on the hydrodynamic load distribution data, blade deformation characteristic values ​​are generated. The blade deformation characteristic value is a quantitative index calculated by real-time monitoring of load distribution and compared with a preset deformation tolerance range. The preset deformation tolerance range is set according to material properties and propeller design standards, such as 80% to 90% of the material's yield strength. If the blade deformation characteristic value exceeds this range, a command-based dynamic mesh reconstruction mechanism is activated. After the dynamic mesh reconstruction mechanism is triggered, pressure pulsation data from the fluid-structure interface is extracted, including pressure amplitude and frequency components. Combined with material property data, the stress concentration factor of the blade structure is calculated. The stress concentration factor is analyzed using the finite element method to identify local stress peak regions. Based on the calculation results, an abnormal stress state signal for the blade is generated and stored as a Boolean variable in the control system. When the abnormal stress state signal is triggered, the eddy current field tracking module is invoked. The eddy current field tracking module updates the local refinement parameters of the fluid domain mesh model based on eddy current theory. These local refinement parameters include mesh node density and element size adjustment coefficients, aiming to enhance simulation accuracy. The local refinement parameters are the core parameters used to control the mesh fineness in different regions of the fluid domain mesh model. They mainly include mesh node density, cell size adjustment coefficient, and local mesh refinement level parameters.The generation process is triggered by anomaly signals in the stress state: the coordinates of the high-risk deformation monitoring sub-region of the blade are located based on the anomaly signals, and the fluid velocity gradient data and pressure pulsation peak data of the region are acquired simultaneously; these data are input into an adaptive mesh optimization algorithm, which determines the priority ranking of the mesh refinement region by constructing a correlation matrix between the fluid separation point position on the blade surface and the pressure pulsation peak value, and performs coupled iteration with the boundary layer thickness parameter and the blade tip clearance value to generate local mesh refinement level parameters; then, based on these level parameters, the distribution rules of the mesh node density and the adjustment coefficient of the element size are further determined to form local refinement parameters that can be directly used to update the mesh model.

[0061] Example 1: See Figure 2 The processing of three-dimensional geometric parameters is based on the blade tilt angle distribution curve, the hub diameter, and the blade tip clearance. The blade tilt angle distribution curve is represented by a parametric equation, which includes radial position variables from the hub to the blade tip and outputs the tilt angle value at the corresponding position point. The parametric equation is fitted to discrete design points using cubic spline interpolation to form a continuous spatial curve. The hub diameter is used as a basic structural parameter input to the geometric model to determine the fixed diameter range of the propeller's central axis region. The blade tip clearance directly constrains the minimum distance parameter between the blade tip and the fluid domain shell, and this parameter is integrated into the initial model through geometric Boolean operations. The application of particle image velocimetry technology uses a dual-pulse laser light source and a high-speed camera system. The laser light source is projected onto the flow field region to form a sheet light source plane, and the high-speed camera captures images of the tracer particle motion at microsecond intervals. The tracer particles are barium titanate microparticles with a particle size distribution controlled within the range of 10-15 micrometers to adapt to underwater flow field tracking. Through continuous frame image analysis, the particle displacement vector field is calculated using a cross-correlation algorithm to identify the location of the boundary layer separation point on the blade surface. The separation point location data is stored as a dataset containing three-dimensional coordinates, timestamps, and local velocity vectors.

[0062] The spatial coordinate mapping between the blade inclination distribution curve and the fluid separation point location data is achieved based on Euclidean transformation. The inclination curve is discretized into a sequence of points, with each point associated with a radial position attribute. For each separation point coordinate, the perpendicular foot of its projection onto the radial position plane is calculated, establishing a mapping relationship with the discrete inclination points. The mapping relationship is optimized and registered using the least squares method, outputting the deviation matrix between the design inclination angle and the actual inclination angle for each separation point. Based on this deviation matrix, the weighted rate of change of the lift coefficient and drag coefficient is calculated, and this rate of change is dimensionless to output an aerodynamic efficiency reference value. Specifically, the momentum-blade element theory is used to derive the benchmark efficiency model, converting the offset between the actual separation point location and the theoretical location into a profile loss coefficient, which is finally integrated into a scalar aerodynamic efficiency reference value. The blade angle of attack gradient parameters of the initial blade fluid domain grid model are updated using a lookup table method. An initial gradient parameter table is constructed with radial position as the independent variable and angle of attack as the dependent variable. Linear corrections are applied to the data in the table based on aerodynamic efficiency reference values: when the aerodynamic efficiency reference value is lower than the theoretical value, the angle of attack is increased by a gradient of 0.05 degrees / meter; when it is higher than the theoretical value, the angle of attack is decreased by a gradient of 0.03 degrees / meter. The mesh boundary conditions are then re-discretized after updating the angle of attack parameters.

[0063] The construction of multiphase flow coupled boundary conditions begins with the extraction of cavitation effect characteristic frequencies. Pressure fluctuation data for 30 consecutive rotational cycles is recorded during transient flow field solver operation. The Welch power spectrum estimation method is used to process the pressure data, identifying spectral peak positions with a 10Hz frequency resolution. The dominant frequency is defined as the frequency corresponding to the maximum amplitude, and the secondary frequency is the frequency component with the second largest amplitude and an interval of more than 10Hz from the dominant frequency. The natural frequencies of the blade structure are obtained through finite element modal analysis, calculating the first six modal frequency values. Spectral superposition analysis is implemented using a frequency response function: an input spectral matrix containing cavitation characteristic frequencies and a transfer function matrix of the structure's natural frequencies are constructed, and a coupled spectrum is generated through convolution. The resonance threshold range is set to a bandwidth of ±15% of the natural frequency; any coupled spectral component exceeding 200% of the reference amplitude is considered out of range. The eddy current field suppression control signal is output as a binary command, with "1" indicating that the suppression mechanism needs to be activated. Turbulent viscosity coefficient adjustment is performed based on the Reynolds stress transport equation. The transient flow field solver incorporates a k-ωSST turbulence model. When it receives a control signal, it calculates the turbulent viscosity coefficient. The time-average value is multiplied by a gain factor of 1.25, while its maximum value is limited to no more than 1.8 times the initial value. After the coefficient is updated, the eddy viscosity sub-model needs to be reinitialized, and a new round of iterative calculations needs to be started until the flow field stabilizes. This process fully realizes the dynamic response mechanism of boundary conditions, forming a closed-loop control loop.

[0064] The hardware configuration for particle image velocimetry includes two 20-megapixel high-speed cameras, positioned at a 75° solid angle outside the observation window of the experimental water tank. The laser sheet thickness is set to 1 mm, and the pulse interval is adjusted between 50-200 μs based on the flow velocity. Image analysis employs an adaptive multi-channel cross-correlation algorithm, with an initial window size of 32×32 pixels, refined to 16×16 pixels through three iterations. The separation point criterion is defined as the location where the negative value of the dot product of the wall friction stress vector and the main flow vector first exceeds a critical value. Spatial registration of the tilt angle-separation point mapping uses the Iterative Closest Point (ICP) algorithm, with an iteration convergence tolerance of 0.1 mm. The aerodynamic efficiency reference value is converted into a correction coefficient through dimensional analysis. The quantitative relationship between it and efficiency is as follows:

[0065]

[0066] in: The difference between the actual and theoretical lift coefficients The drag coefficient is used. Angle-of-attack gradient correction employs a piecewise linear function, with control points set at 10% chord length intervals in the radial position. The confidence interval test for the cavitation characteristic spectrum uses the Bootstrap resampling method, repeated 1000 times to ensure the peak frequency detection error is less than 1Hz. The turbulent viscosity coefficient adjustment process incorporates a PID controller; when the adjusted eddy viscosity leads to a higher turbulent viscosity ratio... When the value exceeds 500, a 0.9-fold attenuation feedback coefficient is triggered to maintain numerical stability. The entire implementation process is encapsulated through modular functions, and data from each step is transmitted in JSON format, forming a standardized data processing chain.

[0067] Example 2: See Figure 3 The blade structure is divided using a geometric topological decomposition algorithm, discretizing the blade spanwise and chordwise directions into rectangular monitoring units. Spanwise division is based on the radial distance from the hub to the blade tip, with partition lines set at 5% chord length steps; chordwise divisions are generated along the direction from the leading edge to the trailing edge at 10% chord length intervals. All partition lines intersect orthogonally after projection into three-dimensional space, forming a mesh-like monitoring sub-region. Each sub-region is assigned a unique code identifier, and its spatial coordinate range is stored in a deformation database. Real-time load fluctuation amplitude acquisition is achieved through a distributed sensor array, with miniature strain gauges placed at the center point of each sub-region, recording the normal and tangential stress components at a sampling frequency of 500Hz. Load data processing employs a sliding time window algorithm: a 0.2-second time window is set to calculate the peak-to-valley difference of the stress waveform, and the average of 10 consecutive cycles is taken as the fluctuation amplitude output value. The comparison analysis between this amplitude and the material yield strength is performed in an independent calculation thread, with the material yield strength value calling the standard value of the corresponding grade from the material property library.

[0068] The local deformation risk coefficient is calculated based on a relative deviation model. For each sub-region, a fluctuation amplitude is established. With yield limit Ratio function: when When the output coefficient is 0; when When linear interpolation is used, it is mapped to the interval (0,1); when A 1.2x over-limit flag is directly output. The risk coefficients of all sub-regions are updated in real time to a two-dimensional matrix, with the matrix row and column indices corresponding to the blade spanwise and chordwise coordinates. A weighted cumulative calculation is initiated using the material fatigue database interface to extract the historical cumulative damage amount for each sub-region. The damage amount is converted into weighting coefficients. Areas that have not reached the fatigue threshold Mildly damaged area Severely damaged areas Final deformation eigenvalues The generation formula is a weighted average: the numerator is the risk coefficient of each sub-region. With weight The sum of the products, with the sum of the weighting coefficients in the denominator. When the eigenvalue exceeds a preset threshold of 0.85, the system activates a red alarm signal.

[0069] The high-risk sub-region localization employs a spatial gradient optimization algorithm. First, the risk coefficient matrix is ​​denoised using Gaussian filtering. Then, the partial derivative matrices of the risk values ​​in the spanwise and chordwise directions are calculated. Points where the product of the two partial derivatives exceeds a critical value are designated as core risk sources. Three-dimensional coordinate interpolation is performed by combining the risk values ​​of eight neighboring points, achieving a localization accuracy of 0.1 mm. After the target region is determined, the flow field data extraction module is activated: the latest time-step velocity tensor field is obtained from the transient solver buffer, and the velocity gradient within a 10-grid area surrounding the target sub-region is extracted. Pressure pulsation data are collected using a boundary layer probe array. The probes are arranged in a star-shaped radial pattern centered on the target area, with the spacing scaled proportionally to the grid size, and record the peak-to-peak value and dominant frequency components of the pulsations.

[0070] The adaptive mesh optimization algorithm includes a fourth-order precision solver. Input parameters are used to construct a 3D data package: velocity gradient data is converted into a vorticity field. The curl modulus distribution and pressure fluctuation data are converted into an energy spectral density function. Initial mesh quality evaluation metrics are set as follows: element aspect ratio threshold of 2.5, skew angle threshold of 20 degrees, and Jacobian determinant lower limit of 0.6. The algorithm first constructs an objective function in the parameter space: vorticity gradient change rate accounts for 40% weight, pressure fluctuation gradient accounts for 35%, and historical deformation records account for 25%. The optimal refinement scheme is found using a quasi-Newton iterative method, with eight iterations. Each iteration updates three control parameters: the mesh node multiplication factor is adjusted between 1.2 and 3.0, the boundary layer element height reduction ratio is set to 20%-50%, and the transition zone growth rate is limited to 1.05-1.15. Finally, a local mesh refinement level parameter matrix is ​​output, where matrix elements correspond to the refinement strength coefficient at different spatial locations.

[0071] The dynamic topology engine is invoked during the mesh model reconstruction phase. The encryption level parameter matrix is ​​parsed: regions with coefficients between 1.0 and 1.5 maintain their original resolution; regions between 1.5 and 2.0 undergo Level 1 encryption (each element is split into four parts); regions above 2.0 undergo Level 2 encryption (each element is split into sixteen parts). The reconstruction process employs octree spatial decomposition technology: the target region is encapsulated as a cubic bounding box, and a spatial index tree is built according to the encryption level. Node insertion follows the Delaunay criterion, and each new node connects topological relationships using a tetrahedral mesh generation algorithm. After local updates, global smoothing is performed: the Laplacian smoothing operator is used to optimize node positions, iterating three times to ensure the mesh quality indicators meet the standards. The reconstruction time is controlled within 1 / 3 of the fluid computation time step interval to ensure the continuity of the simulation time series. The node distribution density of the new mesh model reaches three times the initial density in the risk core region, changes with an exponential decay gradient in the transition region, and remains unchanged in the edge region. The deformation characteristic value calculation thread synchronously updates the mesh topology data and enters the next monitoring cycle.

[0072] The flow field data coupling mechanism employs dual-buffered memory management. During pressure pulsation peak acquisition, a ring-shaped data buffer stores pulsation waveforms for 200 time steps. Waveform processing utilizes a moving-window Fourier transform: a window length of 60 milliseconds, a sliding step size of 10 milliseconds, and a spectral resolution of 1 Hz. The dominant frequency is determined by harmonic components whose amplitude exceeds 60% of the fundamental amplitude; if multiple harmonics exist, the three frequencies with the highest amplitudes are recorded. Adaptive optimization algorithm input interface configuration data preprocessing: eigenvalue decomposition is performed on the velocity gradient tensor to extract the direction and intensity of the maximum shear rate; the pressure pulsation spectrum is converted to a 1 / 3 octave band spectrum, and the energy proportion weight of each frequency band is calculated. Historical deformation records are retrieved from the elastic deformation database: the residual deformation vectors from the five most recent reconstructions at that location are read, and deformation trend parameters are predicted using an autoregressive model. The entire data processing chain is managed by a real-time operating system, with timing errors controlled within microseconds. Data transmission between modules uses shared memory mapping to avoid inter-process communication latency.

[0073] Example 3: See Figure 4 The initialization of the vortex field tracking module is based on an updated fluid domain grid model, which includes dynamically reconstructed locally refined regions. The node density of the grid model reaches its highest value in the blade tip region and decreases exponentially along the radial direction towards the hub. The vortex coreline extraction algorithm uses the λ² criterion, identifying the vortex structure by calculating the second invariant of the velocity gradient tensor. In the three-dimensional flow field, spatial points with negative λ² values ​​constitute the vortex core candidate point set, and a continuous spatial curve is generated by fitting using the moving least squares method. Each coreline is timestamped, recording its spatial evolution trajectory over five consecutive rotation cycles. The evolution trajectory data is stored as a parametric spline curve, including differential geometric features such as curvature and torsion. The preset vortex decay model adopts an empirical parametric form, defining the vortex intensity decay law along the axial direction. The model input parameters include Reynolds number, tip speed ratio, and angle of attack distribution, and the output is the theoretical vortex decay curve.

[0074] The similarity matching process is implemented using a dynamic time warping algorithm. The actual measured vortex trajectory and the theoretical decay curve are nonlinearly aligned along the time axis, and the spatial distance integral at corresponding points is calculated. The matching index is defined as the normalized reciprocal of the distance error; a correction mechanism is triggered when this value falls below a preset threshold of 0.85. The generation of the tip vortex correction factor employs an adaptive weighting strategy: higher weights are assigned to trajectory segments with larger matching errors, and lower weights to segments with smaller errors. The correction factor is ultimately expressed as the harmonic mean of the weighted errors of each segment, with its value controlled between 0.5 and 2.0. The coefficients of the vortex viscosity equation are adjusted proportionally based on the correction factor, and the scaled coefficients are immediately updated to the turbulence model parameter library of the transient solver.

[0075] The monitoring of turbulent kinetic energy decay rate was accomplished using a virtual sensor array deployed in the wake field. Three rows of measuring points were arranged along the flow direction, spaced at intervals of 0.2, 0.5, and 1.0 times the blade diameter. Time-series data of turbulent kinetic energy were recorded at each measuring point, and the decay characteristics were analyzed using an autocorrelation function. The decay rate was calculated as the logarithmic decrease in peak turbulent kinetic energy between adjacent measuring points divided by the flow-direction distance. The vortex core dissipation compensation coefficient was obtained by comparing the measured decay rate with the theoretical value: when the measured value was faster than the theoretical prediction, the coefficient was negative to enhance eddy viscosity; when the measured value was slower than the theoretical prediction, the coefficient was positive to suppress dissipation. This coefficient, in conjunction with the correction factor, forms a dual regulation mechanism.

[0076] The adjustment of the iteration step size for the coefficients of the eddy viscosity equation follows stability constraints. The adjusted step size factor is defined. for:

[0077]

[0078] in: This is the stability constant (default value 0.3). This is the vortex core dissipation compensation coefficient. This formula ensures that the step size adjustment and the physical properties of the flow field remain in dynamic balance. When the absolute value of the compensation coefficient is large, the step size is automatically reduced to ensure numerical convergence; when the compensation coefficient approaches zero, the step size is restored to the default value to improve computational efficiency. The residual change rate is monitored in real time during the adjustment process. If the residual decrease rate is less than 5% for three consecutive iterations, a secondary step size optimization procedure is triggered.

[0079] The spatial evolution analysis of the blade wake vortex structure employs Lagrange particle tracking technology. A swarm of labeled particles is released into the flow field, with their initial positions distributed within a 0.1-diameter radius behind the blade tip. Particle trajectories are calculated using a fourth-order Runge-Kutta integral, recording position and velocity data at each time step. The instantaneous position of the vortex core is extracted from the trajectory data, constructing a time-varying vortex core envelope. The geometric parameters of the envelope include the centerline radius of curvature, cross-sectional eccentricity, and axial elongation, which are used to assess the stability of the vortex flow field. When a sudden decrease in the radius of curvature or a sudden increase in eccentricity is detected, the local mesh refinement level is automatically increased to improve the analytical accuracy of the vortex core region.

[0080] The data interaction between the eddy current field tracking module and the dynamic mesh system adopts a two-way coupling method. The vortex kernel characteristic parameters output by the tracking module are transmitted to the mesh controller in real time, triggering anisotropic refinement in the local region. The refinement intensity is determined according to the eigenvalue distribution of the vorticity gradient tensor: the mesh size in the direction of the maximum eigenvalue is compressed to 1 / 3 of the original, while the orthogonal direction maintains the original resolution. At the same time, the mesh system feeds back the node coordinate information of the refined region to correct the spatial reference frame for vortex identification. This coupling mechanism effectively solves the problem of capturing large deformation eddy current fields, controlling the total number of meshes while maintaining computational accuracy.

[0081] The timing control of the module's operation is embedded in the iterative loop of the main solver. The eddy current field analysis subroutine is launched immediately after each physical time step. The analysis process employs a multi-level parallel computing strategy: eddy core identification is accelerated by the GPU, decay rate calculation is performed by multiple CPU threads, and data synchronization is achieved through shared memory. Under typical operating conditions, the entire tracking process takes less than 20% of the physical time step, ensuring real-time performance requirements. Historical data storage uses a sliding window management system, retaining complete field information from the most recent ten time steps for trend analysis and anomaly detection.

[0082] The calculation of the material fatigue cumulative damage value is based on the improved Miner's linear damage theory. The stress time history of each deformation monitoring sub-region is decomposed into multiple stress cycles, and the cycle counting is realized by the rain-flow method. For each identified stress cycle, the corresponding fatigue life value is obtained by querying the material S-N curve. The damage increment is calculated as the ratio of the number of cycles to the number of failure cycles, and the damage values of all cycles are accumulated to obtain the current total damage amount. The sub-region damage weight coefficient is set according to the damage amount classification: 1.0 for the undamaged region, 1.2 for the slightly damaged region (0 < D ≤ 0.3), 1.5 for the moderately damaged region (0.3 < D ≤ 0.7), and 2.0 for the severely damaged region (D > 0.7). The weight coefficient is updated every five rotation cycles, and the time-varying characteristics of the damage rate are considered during the update.

[0083] The calculation process of the weighted deformation risk value introduces spatial correlation correction. The risk coefficients of adjacent sub-regions are smoothed by the Gaussian kernel function, and the bandwidth of the kernel function is adaptively adjusted according to the grid size. The processed risk value matrix is transformed to the frequency domain by the discrete cosine transform, and after filtering out the high-frequency noise components, it is transformed back to the spatial domain. The overall blade deformation eigenvalue finally takes the 95th percentile of the weighted risk values of all sub-regions, and this statistic can effectively reflect the extreme risk situation. When the eigenvalue exceeds the threshold, the system automatically generates a risk distribution cloud map, marking the exact coordinates and damage levels of the high-risk regions.

[0084] The hardware acceleration of the eddy current field tracking module adopts a hybrid programming architecture. The core algorithm is implemented by CUDA and runs on the NVIDIA Tesla computing card, and the data preprocessing and postprocessing are executed in parallel by the multi-core CPU. The data transfer between the CPU and the GPU is optimized through the PCIe3.0 channel, and the asynchronous copy and memory pinned technology are used to reduce the latency. The dynamic scheduling of the computing tasks is managed by OpenMP, and the workload of each thread is balanced according to the real-time load. The module has a built-in self-diagnosis function, regularly checking the memory usage rate, calculation deviation, and timing performance, and an automatic recovery program is triggered in case of abnormal conditions. The entire system maintains a stable resource occupancy rate during continuous operation, and the peak video memory consumption is controlled within 8GB.

[0085] The standardized data interface established during the implementation process supports multi-physics field coupling. The flow field data is stored in the HDF5 format, including complete metadata descriptions. The grid topology information is organized by the CGNS standard to ensure compatibility with other CAE software. The damage database is built on the SQLite relational architecture, supporting complex queries and transaction processing. All data exchange protocols follow the IEEE15939 standard, meeting the requirements of engineering data traceability while ensuring processing efficiency. The corrected parameters and monitoring indicators output by the module are integrated into the total control system through the OPCUA protocol and participate in the closed-loop control of the entire thruster digital twin.

[0086] Example 4: Adaptive mesh optimization begins with the correlation analysis of fluid separation point location data and pressure pulsation peak data on the blade surface. Separation point data is acquired using particle tracking technology, including three-dimensional coordinates and a scalar value of separation intensity; pressure pulsation data comes from a dynamic pressure sensor array, recording the maximum value of the pressure time history curve at each monitoring point. A 256×256 correlation matrix is ​​constructed, with rows corresponding to chordal coordinate partitions and columns corresponding to spanwise coordinate partitions. Each matrix element is the covariance coefficient between the number of separation points and the pressure peak value within that grid cell. The covariance is calculated using a synchronous time window overlay method, with the time window width covering three blade rotation cycles. The correlation matrix is ​​processed by singular value decomposition to obtain principal components, and the eigenvalues ​​of the first ten principal components are used for spectral clustering analysis. The clustering results divide the blade surface into four characteristic regions: a leading-edge strongly coupled region, a mid-chord weakly coupled region, a trailing-edge transition region, and a blade tip separation region.

[0087] The priority ranking of grid-reinforced regions is based on the energy concentration of feature regions. A priority index is defined, see Table 1.

[0088] Table 1: Mesh Encryption Priority Parameter Table.

[0089] Region Type Boundary layer thickness coefficient Leaf tip gap coefficient Priority Index Leading edge strong coupling region 0.82 0.15 0.93 Leaf tip separation zone 0.65 0.38 0.79 Trailing edge transition zone 0.43 0.22 0.61 Mid-sine weak coupling region 0.29 0.07 0.34

[0090] The priority index is calculated through a nonlinear combination of the boundary layer thickness parameter and the blade tip clearance value. The boundary layer thickness parameter is extracted from the flow field boundary layer probe group, and the blade tip clearance value is measured in real time by a laser displacement sensor. A coupled iterative calculation is initiated using a constrained optimizer: the objective function is set as the weighted sum of the pressure pulsation gradient change rate and the vorticity gradient, with constraints including a 25% upper limit on mesh distortion and a 300% upper limit on the node number growth rate. A sequential quadratic programming method is used for five iterations, updating the local mesh refinement level parameters in each iteration.

[0091] The eddy current field tracking module synchronously acquires turbulent kinetic energy attenuation rate data during operation. Turbulent kinetic energy sampling planes are arranged at positions 0.2D, 0.5D, and 1.0D (D is the diameter) behind the blade, with an 8×8 measurement point matrix set on each plane. Three-dimensional velocity components are acquired at each measurement point using a hot-wire anemometer, and continuous 50-millisecond data segments are recorded at a sampling rate of 100kHz. The turbulent kinetic energy attenuation rate is calculated using the energy spectrum integration method: first, the velocity signals at each measurement point are subjected to FFT transformation to generate an energy spectrum, and the integral value of the 1-10kHz frequency band is taken as the turbulent kinetic energy metric. The attenuation rate is defined as the logarithmic attenuation coefficient of the peak turbulent kinetic energy between adjacent sampling planes, and the attenuation per unit distance is taken as the average of three rotation cycles.

[0092] The vortex core dissipation compensation coefficient is generated using a proportional-derivative (PD) control strategy. A standard attenuation rate baseline curve is established, with the standard reference value obtained from a table based on the Reynolds number. The compensation coefficient calculation process is as follows: the measured attenuation rate is subtracted from the standard attenuation rate to obtain the deviation; the deviation is passed through a first-order low-pass filter to eliminate high-frequency noise; the filtered result is multiplied by a proportional gain coefficient of 0.8 and then superimposed with the differential component. The differential component comes from the attenuation rate gradient change values ​​of the most recent five samples. The final compensation coefficient is limited to the interval [-0.5, 0.5] and discretized with a step size of 0.1.

[0093] The iteration step size adjustment of the coefficients in the eddy viscosity equation of the transient flow field solver is achieved through command conversion. The compensation coefficient is input into the command parser, and the parser's built-in lookup table outputs the step size adjustment: when the compensation coefficient is in the ±0.1 interval, the basic step size is maintained; in the ±(0.1-0.3) interval, the step size is scaled by a factor of 0.8; and in the ±(0.3-0.5) interval, the step size is scaled by a factor of 0.5. The adjusted iteration step size is automatically injected into the solver's control parameter stream, and a residual monitoring mechanism is triggered simultaneously. The residual monitoring window width is set to 10 iteration steps. When the residual decrease rate is lower than the set value for three consecutive iteration steps, the step size relaxation mode is automatically triggered: the increment is increased by 20% based on the current step size value, and the iteration counter is reset.

[0094] Monte Carlo resampling was used to verify the confidence of cavitation characteristic spectra. 1000 subsamples were randomly selected from the original pressure pulsation data, each containing 80% of the original data. The Welch spectrum estimation process was repeated for each subsample: a Hamming window function was used with a window length of 2048 sampling points and a 50% overlap. The peak frequency of each resampled sample was recorded in a frequency matrix with dimensions of 1000×8 (corresponding to the maximum harmonic order of the cavitation characteristic frequency). The standard deviation and confidence interval for each frequency column were calculated: the frequency variation range at a 95% confidence level was used, and the frequency drift threshold was set to 0.8 Hz. When the frequency of a resampled sample deviated from the mean by more than the threshold, that subsample was marked as invalid data.

[0095] The PID controller regulator is designed as a three-channel parallel structure. The proportional channel gain is set to 1.2, the integral time constant is 0.05 seconds, and the derivative gain coefficient is 0.3. The three channel outputs are linearly superimposed after amplitude limiting, and the superposition result is used as the incremental correction term for the turbulent viscosity coefficient. The controller input signal undergoes ramp preprocessing: linear interpolation is performed on abrupt signals within a 0.01-second time window to avoid numerical oscillations caused by step changes. The viscosity adjustment history curve is automatically recorded after each control cycle. The historical data is used for the self-tuning algorithm: when three consecutive overshoots exceeding 5% are detected, a gain coefficient reduction mechanism is triggered, with each reduction being 10% of the current value.

[0096] Data encapsulation and transmission employ a tree-structured organization. Seven main data packets are established: mesh topology data containing node coordinates and element connectivity; flow field parameter packet storing velocity vector field and pressure scalar field; structural response packet recording stress distribution data; control command packet transmitting encryption level parameters; spectral feature packet managing frequency component data; damage assessment packet containing fatigue damage values; and equipment status packet storing sensor calibration values. Each data packet is serialized using JSON-LD format, and metadata is defined using Schema.org lexical definitions. The data transmission cycle is synchronized with the main iteration step size of the fluid solver, and the amount of data transmitted each time is compressed to less than 65% of the original data.

[0097] The vortex core dissipation correction process is equipped with a three-level anomaly handling mechanism. The primary anomaly detection scans for sensor failure states, such as current overruns or signal saturation. The intermediate detection analyzes the rationality of physical quantities, including verification of turbulent kinetic energy nonnegativity and vorticity divergence. The advanced detection performs numerical stability diagnostics, monitoring residual growth rate and iterative convergence speed. When an alarm is triggered at any level, the system automatically switches to safety mode: freezing the current vortex viscosity coefficient, maintaining constant boundary layer parameters, and extending the calculation step size to 1.5 times the baseline value. Simultaneously, a data rollback procedure is initiated: reloading the complete flow field data from the three steps prior to the anomaly time point and reconstructing the computational environment in the isolation zone. The status log generated during the repair process is archived in binary format, containing information such as timestamps, anomaly codes, and recovery trajectories.

[0098] Spatial interpolation of the tip vortex correction factor employs the radial basis function method. A scalar field of the correction factor is created using three turbulent kinetic energy sampling planes as control points. Multiple quadratic surface functions are selected as the basis functions, with the support radius set to 1.5 times the sampling plane spacing. The interpolation weights are determined by solving a system of linear equations: the left-hand matrix contains the basis function values ​​for the control point spacing, and the right-hand matrix contains the measured correction factor vector. The interpolation results are reconstructed as a continuous field on the 3D vortex kernel model and then transferred to unstructured mesh nodes via barycentric coordinate mapping. The mapped data is smoothed using cubic spline interpolation to eliminate local sharp points. Finally, the correction factor is fused with the vortex viscosity equation coefficients using multiplicative coupling, manifesting as an equivalent viscosity gain coefficient in the source term of the governing equations.

[0099] Example 5: The acquisition of sub-region damage weighting coefficients begins with the collection of stress amplitude spectral density distribution data for the deformation monitoring sub-regions. A triaxial strain rosette sensor array is deployed in each monitoring region, employing a Wheatstone bridge structure and temperature compensation circuitry. A signal acquisition card synchronously records the normal and shear stress components at a sampling frequency of 200 kHz, with each recording session covering ten propeller rotation cycles. The raw stress signal is processed by an eighth-order Butterworth low-pass filter, with the cutoff frequency set to fifty times the blade passage frequency. The filtered data is then segmented and windowed using a Hanning window function, with the window length equal to the rotation cycle duration. Each data segment undergoes a fast Fourier transform, with the frequency resolution adjusted to 1 Hz. The spectral density calculation results are used to construct an amplitude spectrum matrix within the 5 Hz to 1 kHz frequency band, with rows corresponding to frequency components and columns corresponding to each sampling segment.

[0100] The calculation of the equivalent alternating stress cycle number uses a modified rainflow counting method. The continuous stress time history is decomposed into a peak-valley sequence, retaining cycle points whose amplitude exceeds 10% of the maximum stress value. The cycle counting process employs a five-point judgment rule: a convexity-concavity test is performed on five adjacent peak-valley points; four points that meet the closed-loop condition constitute a stress cycle. Each cycle records four characteristic values: peak value, valley value, average stress, and cycle number. Cycle classification is stored in a three-dimensional histogram: the first dimension is the average stress level, divided into ten levels at 5 MPa intervals; the second dimension is the stress amplitude level, divided into twenty levels at 2 MPa intervals; the third dimension is the cycle frequency, divided into five levels at 10 Hz intervals. Each cell of the histogram stores a cycle event counter.

[0101] The matching process for the material SN curves utilizes a temperature-region material database. This database contains a base SN curve and six correction coefficients: a surface finish coefficient ranging from 0.8 to 1.2, a mean stress correction coefficient using the Goodman model, a size correction coefficient ranging from 0.7 to 0.95, a temperature compensation factor for the influence of ambient temperature, a corrosion damage factor set to 1.0 to 1.5, and a residual stress correction coefficient based on X-ray diffraction measurements. The matching algorithm executes in three steps: first, it retrieves the material grade and heat treatment status of the current monitoring sub-region; second, it automatically interpolates the baseline SN curve based on the monitoring point temperature; and finally, it applies the six correction coefficients to synthesize a usable SN relationship curve. The curve data point set is converted into a continuous function using cubic spline interpolation.

[0102] The generation of damage weight coefficient update instructions is based on a linear damage accumulation model. For each stress cycle category, its damage contribution is calculated: the number of cycles divided by the fatigue life value obtained from the SN curve. The cumulative damage value is obtained by summing all category damage values. Weight coefficient assignment rules: when the cumulative damage value is below 0.05, a weight of 1.0 is assigned; in the 0.05-0.3 range, it is linearly interpolated to 1.3; in the 0.3-0.7 range, it is linearly interpolated to 1.7; and when it exceeds 0.7, the weight is fixed at 2.0. This coefficient is refreshed every five rotation cycles, taking into account the damage growth rate factor: if the damage increment in an adjacent cycle reaches 15% of the total in the previous cycle, an additional weight coefficient of 0.1 is added.

[0103] The material performance degradation compensation mechanism implements spectral load effect correction. When the average stress fluctuation exceeds the baseline value by 20%, the overload compensation program is activated: a reduction factor of 0.9 is applied to the slope term of the SN curve. In the high-temperature region (monitoring points above 150°C), the creep-fatigue interaction module is simultaneously activated: the hold-time correction for each cycle is calculated, and an equivalent damage value is added in the tensile holding section. Potential compensation is applied to monitoring points under corrosive environments: an additional damage proportionality coefficient is added based on the relationship between corrosion current density and stress intensity factor. Non-proportional correction is introduced for multiaxial stress states: a non-proportional additional damage value is calculated using phase difference measurement results as a supplement to the linear damage model.

[0104] The data storage architecture employs a hierarchical time-series database. Raw stress waveform data is stored on a high-speed solid-state drive array, preserving the time-domain waveform and timestamps. Spectral data is stored in an in-memory database, maintaining the spectral matrix of the most recent twenty samples. Cyclic counting results are managed using a columnar storage engine, with indexes created by time partitioning. The SN curve library is configured with a dual storage mechanism: hot data is stored in GPU memory, and cold data is placed in a distributed file system. Weight coefficient update logs are recorded using a Write-Ahead Log (WAL) mechanism, supporting breakpoint resumption. The entire system is configured with a three-level caching mechanism: 5 seconds of data are stored in the L1 processor cache, 30 seconds of data are stored in the L2 memory pool, and historical data is archived to the L3 hard disk drive.

[0105] The sensor network is managed using a self-organizing protocol. Sensor nodes in each deformation monitoring sub-region form an independent subnet, transmitting data via the ZigBee protocol. Relay nodes are configured at the hub, using a star topology to aggregate data from each sub-region. The acquisition and control module runs an adaptive scheduling algorithm: the sampling frequency is increased to 400kHz in high-damage areas and reduced to 50kHz in low-damage areas. The data transmission cycle is dynamically adjusted based on network load: a complete dataset is transmitted every 0.5 seconds under normal conditions; when increased risk is detected, the cycle switches to 0.1 seconds; in dangerous situations, a burst mode is activated, transmitting simplified data packets at 0.02-second intervals. Node fault detection employs a dual mechanism of heartbeat packets and data verification; abnormal nodes are automatically isolated, and neighboring nodes take over the monitoring task.

[0106] The long-term damage prediction model employs time series analysis. Cumulative damage value sequences for one hundred consecutive periods are collected for each monitoring sub-region, and an autoregressive integral moving average model is established. The model order is automatically determined according to the AIC criterion, typically using an ARIMA(2,1,1) structure. The prediction process is executed twice per period: during the rotation, parameters are updated based on current data; at the end of the rotation, predicted values ​​for the next ten periods are output. The prediction results are compared with real-time monitoring values, and model recalibration is triggered when the residual exceeds 15%. An artificial neural network predictor is overlaid in high-risk areas: a three-layer fully connected network is constructed, with the input layer containing thirty feature parameters and the output layer providing the probability distribution of damage trend categories.

[0107] Acoustic emission (AE) sensing technology is used to detect phase transitions in the microstructure of materials. Broadband AE sensors with a frequency range of 50 kHz to 1 MHz are added to key monitoring areas. The AE signals are analyzed for waveform characteristics to identify the microscopic deformation mechanisms of the material: continuous emission corresponds to dislocation slip, while sudden emission corresponds to twinning or phase transitions. A mapping model between AE energy and fatigue damage is established: when a phase transition characteristic signal is detected, the weighting coefficient is automatically increased by 0.2. A microstrain amplification algorithm is used in the active phase transition region: based on the measured austenite content, the amplification factor of the local strain readings is adjusted, with a coefficient variation range of 1.0–1.8.

[0108] The system is configured with an automatic diagnostic mode for operation and maintenance. Weekly full system verification is performed: a standard sine wave test signal is input to check the frequency response characteristics of each channel; a standard resistor network is injected to calibrate the Wheatstone bridge; and communication error rate testing and transmission delay measurement are performed. Monthly SN curve database consistency verification is performed: a version comparison is conducted with the central material library, and automatic updates are initiated when a difference exceeds 5%. Sensor base checks are performed during each propeller overhaul: a laser rangefinder verifies the sensor installation position accuracy, and repositioning is performed when the error exceeds 0.1mm. All maintenance records are written to blockchain storage to prevent historical data tampering.

[0109] It should be noted that, in this document, relational terms such as "first" and "second" are used only to distinguish one entity or operation from another, and do not necessarily require or imply any such actual relationship or order between these entities or operations. Furthermore, the terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such process, method, article, or apparatus.

[0110] Although embodiments of the invention have been shown and described, it will be understood by those skilled in the art that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the appended claims and their equivalents.

Claims

1. A simulation-based modeling method for annular propellers of underwater thrusters, characterized in that, Includes the following steps: Based on the three-dimensional geometric parameters and material property data of the annular propeller of the underwater thruster, an initial blade fluid domain mesh model was constructed using the discrete element simulation algorithm. After applying multiphase flow coupled boundary conditions to the initial blade fluid domain mesh model, the hydrodynamic load distribution data of the annular propeller during operation is obtained through a transient flow field solver. Based on the hydrodynamic load distribution data, blade deformation characteristic values ​​are generated. If the blade deformation characteristic values ​​exceed the preset deformation tolerance range, a dynamic mesh reconstruction command is activated. Based on the dynamic mesh reconstruction command, pressure pulsation data of the fluid-structure interaction interface is extracted, and the stress concentration factor of the blade structure is calculated in combination with the material property data. If the stress concentration factor exceeds the preset stress concentration factor safety threshold set based on the material property data, an abnormal stress state signal of the blade is generated. When the abnormal stress state signal of the blade is triggered, the eddy current field tracking module is invoked to update the fluid domain mesh model.

2. The underwater propeller modeling method based on simulation according to claim 1, characterized in that, The three-dimensional geometric parameters include the blade tilt angle distribution curve, hub diameter, and blade tip clearance of the annular propeller. The location data of fluid separation points on the blade surface are collected by particle image velocimetry technology. The blade tilt angle distribution curve is discretized and the rate of change of tilt angle with radial position is calculated to obtain the blade angle of attack gradient parameter; By mapping the blade tilt angle distribution curve with the fluid separation point location data in spatial coordinates, a reference value for the blade aerodynamic efficiency is obtained. The blade angle of attack gradient parameters of the initial blade fluid domain grid model are corrected based on the aerodynamic efficiency reference value.

3. The underwater propeller modeling method based on simulation according to claim 2, characterized in that, The process of constructing the multiphase flow coupled boundary conditions includes: The characteristic frequencies of cavitation effect within the rotation cycle of the annular propeller are extracted, and the characteristic frequencies of cavitation effect are superimposed with the natural frequencies of the propeller structure for spectral analysis. If the result of the spectral superposition analysis exceeds the resonance threshold range, then the eddy current field suppression control signal is output. The turbulent viscosity coefficient of the transient flow field solver is adjusted based on the eddy field suppression control signal.

4. The underwater propeller modeling method based on simulation according to claim 3, characterized in that, The method for calculating the characteristic value of the blade deformation is as follows: The blade structure is divided into several deformation monitoring sub-regions, and the real-time load fluctuation amplitude of each deformation monitoring sub-region is obtained; The deviation between the real-time load fluctuation amplitude and the material yield limit value of the corresponding sub-region is calculated to obtain the local deformation risk coefficient. The blade deformation characteristic value is output by performing a weighted summation operation on all local deformation risk coefficients.

5. The underwater propeller modeling method based on simulation according to claim 4, characterized in that, The execution process of the dynamic mesh reconstruction instruction includes: The coordinates of high-risk deformation monitoring sub-regions are located based on abnormal blade stress state signals. Acquire fluid velocity gradient data and pressure pulsation peak data in the high-risk deformation monitoring sub-region; Fluid velocity gradient data and pressure pulsation peak data are input into an adaptive mesh optimization algorithm to generate local mesh refinement level parameters. The node distribution density of the fluid domain mesh model is updated based on the local mesh refinement level parameters.

6. The underwater propeller modeling method based on simulation according to claim 5, characterized in that, The operation process of the eddy current field tracking module includes: The vortex coreline was extracted based on the updated fluid domain grid model to obtain the spatial evolution trajectory of the blade wake vortex structure. The spatial evolution trajectory is matched with the preset eddy current decay model. If the similarity matching result is lower than the eddy current stability threshold, the tip eddy correction factor is output. Adjust the eddy viscosity equation coefficients of the transient flow field solver based on the tip eddy correction factor.

7. The underwater propeller modeling method based on simulation according to claim 6, characterized in that, The weighted summation method for the local deformation risk coefficient is as follows: Obtain the cumulative material fatigue damage value for each deformation monitoring sub-region and label it as the sub-region damage weight coefficient; The local deformation risk coefficient is multiplied by the corresponding sub-region damage weight coefficient to obtain the weighted deformation risk value. The arithmetic mean of all weighted deformation risk values ​​is processed to output the overall deformation characteristic value of the blade.

8. The underwater propeller modeling method based on simulation according to claim 7, characterized in that, The execution steps of the adaptive grid optimization algorithm include: Construct a correlation matrix between the location data of fluid separation points on the blade surface and the peak pressure pulsation data; The priority ranking of grid encryption regions is determined based on the distribution of eigenvalues ​​in the correlation matrix. The boundary layer thickness parameters of the fluid domain mesh model are coupled and iteratively calculated with the blade tip clearance value according to priority to generate local mesh refinement level parameters.

9. The underwater propeller modeling method based on simulation according to claim 6, characterized in that, When the eddy field tracking module generates the blade tip vortex correction factor during execution, it simultaneously collects the turbulent kinetic energy attenuation rate of the blade tail vortex structure; the difference between the turbulent kinetic energy attenuation rate and the eddy stability threshold is calculated to obtain the vortex core dissipation compensation coefficient. The iteration step size parameter of the vortex viscosity equation coefficients of the transient flow field solver is corrected based on the vortex core dissipation compensation coefficient.

10. The underwater propeller modeling method based on simulation according to claim 7, characterized in that, The method for obtaining the sub-region damage weight coefficient includes: Extract the stress amplitude spectral density distribution data for each deformation monitoring sub-region, and calculate the equivalent number of alternating stress cycles based on the stress amplitude spectral density distribution data; The equivalent alternating stress cycle number is matched with the material SN curve, and the sub-region damage weight coefficient update command is output.

Images