Fan tower multi-source load coupling simulation test method

Abstract

This invention discloses a multi-source load coupling simulation test method for wind turbine towers, comprising: establishing a numerical model of the entire wind turbine system; dividing the numerical model of the entire wind turbine system into a physical substructure to be tested and a numerical substructure to be simulated, wherein the numerical substructure includes a dynamic model of a servo control system; generating a reduced-order numerical substructure model that is dynamically equivalent and can be calculated in real time based on the numerical substructure; in a real-time closed loop, based on external multi-source load input, repeatedly performing the following operations: measuring the actual interface response of the physical substructure at the physical-numerical interface; inputting this response as a boundary condition into the reduced-order numerical substructure model, and calculating the interface target command for the next time step in real time; and applying the interface target command to the physical substructure via a loading system. This invention does not require modification of the simulation software source code, has strong compatibility, and can reproduce the damage evolution and limit state of the structure under multi-source coupled loads.

Description

Technical Field

[0001] This invention relates to the field of structural engineering testing technology, and in particular to a multi-source load coupling simulation test method for wind turbine towers. Background Technology

[0002] Large offshore wind turbine structures operate under complex environmental loads from multiple sources, including wind, waves, currents, and earthquakes, and there is a coupling effect between their servo control system and structural dynamics. Reproducing and evaluating their dynamic response and damage evolution mechanisms under extreme conditions in the laboratory can help improve the safety and economy of structural design.

[0003] Existing hybrid simulation technologies for wind turbine structures typically treat components such as the tower as physical substructures, while establishing the remaining parts as numerical models. To ensure real-time computation, simplified linear models can be used for the numerical substructures. Regarding model order reduction, mainstream methods require invasive modifications to the simulation software's source code to extract the system matrix. At the loading and execution level, some experiments still employ open-loop loading or quasi-static loading methods, failing to establish a high-speed, bidirectional dynamic information loop between the physical and numerical substructures.

[0004] Existing technologies mainly face problems with insufficient fidelity in numerical models and inadequate realism in the coupling of physical processes. Therefore, further research and innovation are needed to address these issues. Summary of the Invention

[0005] The purpose of this invention is to provide a multi-source load coupling simulation test method for wind turbine towers, in order to solve at least one of the aforementioned problems.

[0006] According to one aspect of this application, a multi-source load coupling simulation test method for wind turbine towers includes:

[0007] A numerical model of the wind turbine system is established, and the numerical model of the wind turbine system is divided into a physical substructure to be tested and a numerical substructure to be simulated. The numerical substructure includes the dynamic model of the servo control system.

[0008] Based on numerical substructure, a reduced-order numerical substructure model that is dynamically equivalent and can be computed in real time is generated.

[0009] In the real-time closed loop, based on external multi-source load input, the following operations are repeated: the actual interface response of the physical substructure at the physical-numerical interface is measured; the actual interface response is used as a boundary condition input to the reduced-order numerical substructure model, and the interface target instruction for the next time step is calculated in real time; the interface target instruction is applied to the physical substructure via the loading system.

[0010] Beneficial effects: This invention requires no modification to the simulation software source code, has strong compatibility, and can reproduce the damage evolution and ultimate state of structures under multi-source coupled loads. The related technical effects will be described in detail below with reference to specific embodiments. Attached Figure Description

[0011] Figure 1 A flowchart of a multi-source load coupling simulation test method for wind turbine towers provided in this application embodiment.

[0012] Figure 2 The execution flowchart of the non-intrusive projection framework provided in the embodiments of this application within each computation time step.

[0013] Figure 3 This is a flowchart illustrating the generation of a reduced-order numerical substructure model, as provided in an embodiment of this application.

[0014] Figure 4 A flowchart of a closed-loop real-time simulation including physical feedback provided for embodiments of this application. Detailed Implementation

[0015] To enable those skilled in the art to better understand the present invention, the technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings of the embodiments of the present invention. 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 should fall within the scope of protection of the present invention.

[0016] It should be noted that the terms "first," "second," etc., in the specification and accompanying drawings of this invention are used to distinguish similar objects and are not necessarily used to describe a specific order or sequence. It should be understood that such data can be interchanged where appropriate so that embodiments of the invention described herein can be implemented in sequences other than those illustrated or described herein. Furthermore, the terms "including" and "having," and any variations thereof, are intended to cover non-exclusive inclusion; for example, a process, method, system, product, or apparatus that includes a series of steps or units is not necessarily limited to those steps or units explicitly listed, but may include other steps or units not explicitly listed or inherent to such processes, methods, products, or apparatus.

[0017] To address the aforementioned issues, the applicant conducted in-depth searches and analyses, and discovered:

[0018] First, numerical models are usually static and cannot track the performance evolution of physical specimens, such as stiffness degradation, due to cumulative damage, leading to distortion in the later stages of the simulation. At the same time, open-loop or weakly coupled loading methods cut off the real-time feedback of the structure's actual response to the servo control system, making it impossible to reproduce the transient strong coupling effect between the servo and the structure under conditions such as emergency shutdown.

[0019] Secondly, achieving model order reduction and real-time computation relies on invasive modifications to the simulation software source code, resulting in poor compatibility and high technical barriers. Furthermore, the online calculation of the Jacobian matrix in the implicit integration algorithm will constrain real-time performance.

[0020] In the following text, the py spring is used to characterize the nonlinear relationship between the horizontal resistance of the soil and the horizontal displacement of the pile, the tz spring is used to characterize the nonlinear relationship between the pile side friction and the axial displacement of the pile, and the Qz spring is used to characterize the nonlinear relationship between the pile end resistance and the pile end settlement.

[0021] The Newmark-β method is used in conjunction with the Newmark-β method.

[0022] The primitive numerical substructure solver can be simplified to the primitive solver.

[0023] System modeling and partitioning correspond to establishing a numerical model of the entire wind turbine system and partitioning it.

[0024] The execution process of real-time hybrid loading corresponds to the execution process of real-time closed loop.

[0025] The snapshot matrix corresponds to the full-order response snapshot matrix.

[0026] Real-time hybrid simulation is a simulation experiment.

[0027] Degrees of freedom condensation means retaining only a few key locations while eliminating other internal variables through mathematical relationships. This results in a smaller model and faster computation. Alternatively, it refers to reducing the dimensionality of high-dimensional internal degrees of freedom in a structural or soil numerical model through matrix transformations or modal truncation, retaining only the important degrees of freedom related to external load inputs and boundary interactions. This reduces the model size and improves real-time computational efficiency without compromising computational accuracy.

[0028] New interest ε _k The corresponding sequence is the new information sequence.

[0029] To solve the above problems, combined with Figures 1 to 4 The present invention will be specifically described through the following embodiments.

[0030] On the one hand, the multi-source load coupling simulation test method for wind turbine towers can include main stages such as system modeling and partitioning, reduced-order model generation, real-time hybrid loading, and response analysis. Systems implementing the method of this invention include:

[0031] The physical substructure to be tested, such as a segmental specimen of a wind turbine tower;

[0032] Multi-degree-of-freedom electro-hydraulic servo loading system for applying force and displacement;

[0033] Sensor acquisition systems are used to measure physical quantities such as force, displacement, acceleration, and strain.

[0034] A deterministic real-time computing platform, which can employ a high-speed control platform based on an industrial real-time processor, such as dSPACE, Speedgoat, or NIPXI systems, is used to run reduced-order models and perform closed-loop control.

[0035] Accordingly, the method of the present invention can be implemented by performing the following steps:

[0036] Step 101: Based on the pre-acquired wind turbine structural parameters, servo control parameters, and environmental load definition parameters, establish a numerical model of the wind turbine system. Divide the numerical model of the wind turbine system into a physical substructure to be tested and a numerical substructure to be simulated. The numerical substructure includes the dynamic model of the servo control system.

[0037] Before establishing a numerical model of the entire wind turbine system, it is necessary to collect and organize multi-source parameters in advance, including at least:

[0038] Structural parameters, such as the geometry, material properties, and mass distribution of the tower, blades, and nacelle;

[0039] Servo control parameters, such as the control law, gain, and limiter threshold of the pitch and torque controllers;

[0040] Environmental load definition parameters, such as wind field characteristics, wave spectrum and ground motion records.

[0041] Furthermore, the numerical model of the wind turbine system is a multiphysics coupling model, used to reproduce the dynamic behavior of wind turbines in real marine environments in a computer.

[0042] Optionally, the model can be specifically composed of the following coupled sub-models:

[0043] The aerodynamic model uses the blade element momentum theory (BEM) to calculate the aerodynamic forces of the wind turbine.

[0044] The servo control model uses logic including a proportional-integral (PI) controller to actively control the pitch and generator torque.

[0045] The elastodynamic model uses the finite element method based on Timoshenko beam theory to model the flexible deformation of the tower and blades.

[0046] A hydrodynamic model was used to calculate the wave forces on the tower foundation using the Morrison equations.

[0047] Soil-structure interaction model.

[0048] Alternatively, the soil-structure interaction model can employ py, tz, and Qz nonlinear springs to simulate the nonlinear restoring force characteristics of the foundation.

[0049] In some implementations, to simplify calculations, the soil-structure interaction model can also be represented by an equivalent linearized impedance matrix.

[0050] After the model is built, a virtual-to-real substructure is defined, that is, determining which parts of the model are handled by the physical specimen and which parts are retained in the numerical simulation. This invention retains the part containing the servo control system in the numerical substructure, providing a foundation for realizing servo-structure coupled simulation.

[0051] Step 102: Based on the numerical substructure, generate a reduced-order numerical substructure model that is dynamically equivalent and can be computed in real time.

[0052] Because the number of degrees of freedom of numerical substructures is enormous, typically reaching thousands or even tens of thousands, the time required to solve one step is much greater than the step size required for real-time experiments, such as 1 to 2 milliseconds. Therefore, it cannot be directly used for real-time closed-loop loading.

[0053] In other words, the computation time required for the multiple Newton iterations and the solution of large linear equations in its implicit time integration algorithm is much greater than the step size required by real-time experiments.

[0054] Therefore, it is necessary to reduce the order of the numerical substructure model, that is, to use mathematical transformations to project the high-dimensional dynamic equations into a low-dimensional subspace, which can characterize the main dynamic characteristics of the system. While ensuring the equivalence of dynamic characteristics, the number of computational degrees of freedom is reduced, so that the single-step solution time of the model can meet the requirements of real-time computing.

[0055] After generating the reduced-order numerical substructure model, the model and its offline pre-computed parameterized Jacobian matrix model are deployed together to a deterministic real-time computing platform for no-load closed-loop joint debugging tests. Only after verifying that the real-time computing capability and the response characteristics of the loading system meet the requirements can the formal real-time hybrid loading stage begin.

[0056] The specific implementation methods, such as the non-intrusive order reduction framework and the multi-criteria adaptive modal truncation criterion, will be described in subsequent embodiments.

[0057] Step 103: In the real-time closed loop, based on predefined external multi-source load inputs, repeat the following operations:

[0058] Measure the actual interface response of the physical substructure at the physical-numerical interface;

[0059] The actual interface response is used as a boundary condition input into the reduced-order numerical substructure model. Combined with the external multi-source load input at the current moment, the interface target command for the next time step is calculated in real time.

[0060] The interface target instructions are applied to the physical substructure via the loading system.

[0061] Alternatively, the above steps are executed cyclically on a real-time computing platform with a fixed time step. The external multi-source load input, i.e., the load time series generated based on a predefined target load case matrix, is used to evaluate the wind turbine's service performance throughout its entire life cycle. Specifically, it may include:

[0062] Normal operating conditions, such as normal power generation and normal shutdown; extreme operating conditions, such as encountering extreme turbulent winds and extreme ocean waves;

[0063] Transient operating conditions, such as emergency shutdown and power grid failure;

[0064] Multiple disaster coupled conditions, such as the combined action of wind, waves and earthquakes.

[0065] Specifically, at the beginning of each time step, force sensors and displacement sensors deployed on the physical substructure collect the actual interface response at the physical-numerical interface and transmit it to the real-time computing platform via the data acquisition system. Upon receiving this interface response, the real-time computing platform uses it as boundary conditions, combined with the current external multi-source load input, to drive the reduced-order numerical substructure model to perform a one-step time integration, thereby calculating the interface target command to be applied to the interface at the next time step.

[0066] Furthermore, the target instructions of this interface drive the electro-hydraulic servo loading system to load the physical substructure via a digital-to-analog converter and a controller.

[0067] Based on this, the measurement-solution-loading cycle is repeated continuously at a speed of milliseconds, forming a physical-numerical closed loop, which realizes the real-time simulation of the dynamic response of the wind turbine structure under the coupling effect of complex external loads and internal servo control.

[0068] On the other hand, it describes alternative implementation schemes for the physical-numerical substructure partitioning used to achieve high-fidelity simulation, including:

[0069] Step 201: In system modeling and partitioning, the scope of the physical substructure is defined as including the wind turbine tower body and the foundation transition section connected to the foundation.

[0070] This division incorporates the local nonlinear mechanical behavior caused by flange or grouting connections within the basic transition section into the testing of the physical substructure, thereby improving the simulation fidelity of the damaged area.

[0071] Accordingly, hybrid simulation methods typically cut off the connection between the tower bottom and the foundation, including the entire foundation portion in the numerical model. This approach assumes that the boundary conditions at the tower bottom are linear or modelable; however, the connection area between the wind turbine tower and the foundation, i.e., the foundation transition section, is a region with strong structural nonlinearity, stress concentration, and susceptibility to damage.

[0072] In one alternative partitioning scheme, the lower boundary of the physical substructure is extended downwards to encompass the entire basic transition section. This scheme can capture the complex and difficult-to-accurately numerically modeled local nonlinear mechanical behavior within the basic transition section through physical experiments.

[0073] For example, on a large-diameter monopile foundation with flange connections, the nonlinear behavior within the foundation transition section is mainly manifested in the contact and separation of the flanges, changes in bolt preload, and potential bolt loosening or fatigue fracture. Treating this section as a physical specimen allows for direct measurement of bolt strain and flange opening, thus obtaining the boundary condition response.

[0074] Alternatively, on the basis of a grouted jacket, the nonlinear behavior in the foundation transition section mainly stems from the nonlinear constitutive relationship of the grouting material, the bond-slip between the grouting layer and the sleeve, and the propagation of microcracks under high stress.

[0075] Based on this, the grouting connection section can be used as a physical specimen to reproduce its stiffness degradation and energy dissipation characteristics under cyclic loading, and its fidelity is difficult to achieve with simplified numerical models.

[0076] In summary, by incorporating the basic transition section into the physical substructure, the complex physical processes in the damaged area can be directly placed in the physical testing environment, thereby improving the accuracy of predicting the structural limit state and fatigue life.

[0077] Step 202: In addition to the interface between the top of the tower and the bottom of the foundation transition section, the physical-numerical interface also includes at least one intermediate loading section set along the height direction of the physical substructure.

[0078] The intermediate loading section defines at least one horizontal translational degree of freedom, which is used to obtain the equivalent load command at each intermediate loading section by the reduced-order numerical substructure model in real-time hybrid loading, so as to apply distributed environmental loads such as wind pressure to the physical substructure tube more realistically.

[0079] Since wind loads are not concentrated at the top of the tower, but are distributed non-uniformly along the height of the tower, applying equivalent force and moment only at the top of the tower can simulate the overall bending moment and shear force response of the structure, but cannot accurately reproduce the local stress distribution and higher-order vibration mode response along the tower.

[0080] To address this issue, distributed loading points can be optionally set on the physical substructure. Specifically, in addition to the upper interface defined at the top of the tower connecting to the numerical control module, and the lower interface defined at the bottom of the foundation transition section connecting to the numerical foundation, several locations along the height of the tower are selected as intermediate loading sections. For example, for a 30-meter-high tower physical specimen, an intermediate loading section can be set at heights of 10 meters and 20 meters.

[0081] Each intermediate loading section needs to define at least two horizontal translational degrees of freedom, used to apply horizontal forces in the downwind and crosswind directions, respectively.

[0082] During real-time hybrid loading, the reduced-order numerical substructure model not only calculates the target command at the interface but also simultaneously calculates the equivalent nodal loads obtained from distributed wind pressure integration at each intermediate loading section. The load commands are synchronously sent to the loading actuators located at the corresponding sections for application.

[0083] In some alternative implementations, the intermediate loading section can also increase the torsional degree of freedom about the central axis of the tower.

[0084] In this case, introducing an intermediate loading section allows the distributed wind loads that originally acted on the tower in the numerical model to be reproduced on the physical specimen, making the deformation and stress distribution of the physical specimen closer to the actual situation. This can be used to capture the structural dynamic response under wind shear, tower shadow effect, and higher-order modes.

[0085] Accordingly, in this embodiment, the optional set of degrees of freedom topology D of the physical-numerical interface is... _I This can be formally expressed as:

[0086] D _I ={DOF _top (6),DOF _bottom (6),DOF _mid1 (2),DOF _mid2 (2)};

[0087] Among them, DOF _top (6) Represents the six degrees of freedom of the interface at the top of the tower, including three translational and three rotational degrees of freedom; DOF _bottom (6) represents the six degrees of freedom of the lower intersection interface at the bottom of the basic transition section; while DOF _mid1(2) and DOF _mid2 (2) represents the two horizontal translational degrees of freedom of each of the two intermediate loading sections.

[0088] On the other hand, the technical framework used to explain the generation of the reduced-order numerical substructure model is the non-intrusive projection framework. Specifically, this includes:

[0089] Step 301, the process of generating the reduced-order numerical substructure model, is specifically implemented through a non-intrusive projection framework. The non-intrusive projection framework performs the following operations within each computation time step:

[0090] Map the state vector in the reduced-order degree-of-freedom space at the current moment to the full-order degree-of-freedom space;

[0091] The original numerical substructure solver is invoked to solve the nonlinear terms, including aerodynamic forces, hydrodynamic forces, and servo control forces, in one step in the full-order degree-of-freedom space to obtain the full-order residual forces.

[0092] The full-order residual force is projected back into the reduced-order degree-of-freedom space to solve for the state vector at the next moment.

[0093] The complete internal nonlinear computation mechanism is preserved without modifying the original solver source code upon which the numerical substructure depends.

[0094] The initial state vector is determined based on preset initial conditions, and the state vectors at subsequent times are the solution results of the previous time step.

[0095] This step addresses the issue of traditional model reduction methods, such as the standard Galerkin projection method, which require intrusive modifications to the source code of large simulation software to extract its internal mass, stiffness, and damping matrices.

[0096] The simulation software can be engineering-grade wind turbine aerodynamics and load simulation software, such as OpenFAST and Bladed.

[0097] For example, the non-intrusive projection framework calls the original numerical substructure solver as a black box, and its implementation can be further decomposed into the following iteratively executed sub-steps:

[0098] Accordingly, a mapping extension of the state vector is performed. At any computation time step t _k This yields the state vector, i.e., the coordinate vector q, in the reduced-order degree-of-freedom space. _k By using the pre-calculated projection basis matrix Φ, the low-dimensional state vector is mapped to the full-order degree-of-freedom space to obtain the full-order physical displacement vector u. _k Their relationship can be expressed as:

[0099] u _k=Φ×q _k .

[0100] In the above, the displacement patterns of all nodes in the numerical substructure in physical space can be reconstructed from a few modal participation coefficients.

[0101] Furthermore, the original solver is invoked to solve for the full-order residual force. The reconstructed full-order physical displacement vector u _k and the velocity vector u obtained by the time integration algorithm _k’ As input parameters, the original, unmodified numerical substructure solver is invoked. This solver performs mechanical calculations; for example, the aerodynamic module calculates unsteady aerodynamic forces based on the current blade position and wind speed, the servo module adjusts the pitch angle based on rotational speed feedback, and the hydrodynamic module calculates wave forces, etc.

[0102] After performing one step of calculation, it outputs the full-order residual force vector r at the current moment. _k This characterizes the degree to which the system does not satisfy dynamic equilibrium in the current state, and can be expressed as:

[0103] r _k =F _ext (u _k ,u _k’ ,t _k )-F _int (u _k ,u _k’ );

[0104] In the formula, F _ext Let F be the external load vector. _int This is the force vector within the system, which includes inertial force, damping force, and elastic restoring force.

[0105] Furthermore, the full-order residual force is projected back into the reduced-order space. The full-order residual force vector is then transformed by the transpose Φ of the projection basis matrix. T Perform a left multiplication, i.e., execute the Galerkin projection operation, to obtain the reduced residual force vector.

[0106] According to the principles of dynamics, the reduced-order residual force vector must be zero, i.e., Φ T ×r _k =0. This equation is the reduced-dimensional dynamic equation that needs to be solved on a real-time computing platform.

[0107] The coordinate vector at the next time step can be obtained by using a time integration algorithm, such as the implicit Newmark-β method, on this low-dimensional equation.

[0108] By employing an expansion-solve-compression loop, the coupled dynamics characteristics can be preserved in the reduced-order model by utilizing its nonlinear computational capabilities without affecting the internal implementation of the original solver. Furthermore, the non-intrusive approach makes the method of this invention both engineering-practical and software-compatible, allowing integration with various commercial or open-source professional simulation software.

[0109] As an example, when constructing the projected basis matrix, it is necessary to determine which vibration modes with internal degrees of freedom should be retained, i.e., the multi-criteria adaptive modal truncation criterion. This can be achieved using the following method:

[0110] Before applying this criterion, a data foundation for analysis, namely a full-order response snapshot matrix, needs to be generated. Accordingly, offline real-time simulation is performed using a numerical model of the entire wind turbine system.

[0111] For example, the load cases used in the simulation should cover the target load case matrix, especially extreme load cases and multi-hazard coupled load cases containing rich frequency components. During the simulation, the displacement response time series of all internal degrees of freedom of the numerical substructure are recorded at high frequency. The time series data are then compiled into a large matrix, which is the full-order response snapshot matrix. This matrix reflects the dynamic behavior characteristics of the numerical substructure under various excitations.

[0112] In some embodiments, the projected basis matrix may be constructed using a Craig-Bampton-like method to preserve information about the degrees of freedom at the physical-numerical interface and to perform modal condensation of the internal degrees of freedom. In this case, the full-order displacement vector u is divided into internal degree-of-freedom displacements ui. _i and interface degree of freedom displacement u _b Then its reduced-order coordinates q _k The relationship is:

[0113] u _i =Φ _k ×q _k +Φ _c ×u _b ;

[0114] In the formula, Φ _k It is a matrix composed of the fixed interface mode vectors to be screened, Φ _c It is a static constraint mode matrix.

[0115] Step 401, generating a reduced-order numerical substructure model, also includes:

[0116] The mode order for preserving the degrees of freedom within the numerical substructure is determined by using a multi-criteria adaptive criterion.

[0117] Based on the determined preserved mode order, modal truncation and degree of freedom condensation are performed on the numerical substructure to obtain a reduced-order numerical substructure model.

[0118] The multi-criteria adaptive criterion comprehensively evaluates the spectral characteristics of multi-source load inputs and the energy contribution of each mode under load excitation, enabling the reduced-order model to have sufficient accuracy in characterizing multi-source coupling effects such as wind, waves, earthquakes, and servo control.

[0119] In this step, the wind, wave, seismic, and servo control loads faced by the wind turbine structure exhibit vastly different energy distribution characteristics in the frequency domain. For example, wind load energy is mainly concentrated in the low-frequency band, wave load in the mid-to-low-frequency band, while seismic load and some servo control transient responses may excite structural vibrations over a considerable frequency band.

[0120] In some embodiments, a single-mode cutoff criterion is employed, such as simply selecting a cutoff frequency. This may result in the omission of higher-order modes that are crucial to the response to a certain load, or the retention of redundant modes that contribute little to the overall response. Therefore, this invention employs a multi-criteria comprehensive evaluation method to improve the rationality and completeness of the selected mode set from multiple dimensions.

[0121] Step 402, the multi-criteria adaptive criterion specifically includes at least one of the following criteria:

[0122] Based on the power spectral density analysis of multi-source load input, modes whose natural frequencies cover the main frequency distribution bands of load energy were initially screened.

[0123] Based on the coordinate response energy of each mode under different load conditions, the modes whose cumulative energy is greater than the linear superposition value of the response under each single-source independent excitation are further screened out and forcibly retained under multi-source joint excitation, thus capturing the nonlinear coupling amplification effect between multi-source loads.

[0124] Alternatively, based on the modal responses of each order obtained by simulating the logarithmic substructure under multi-source joint excitation and each single-source independent excitation, the mode whose response under multi-source joint excitation is greater than the product of the linear superposition value of the response under each single-source independent excitation and the preset coupling amplification threshold is identified and forcibly retained, thereby capturing the nonlinear coupling amplification effect between multi-source loads.

[0125] Specifically, the multi-criteria system is implemented according to a hierarchical logic from coarse to fine. Accordingly, a multi-source load spectral coverage criterion is applied for initial screening. This criterion ensures that the frequency range of the selected modes covers the main distribution intervals of all external load energy.

[0126] In practice, Fourier transforms are performed on the time series of loads such as wind, wave, and earthquake in the target load case matrix to calculate their power spectral density functions. The power spectral densities of each load are then superimposed to obtain the overall energy-frequency distribution map. The upper limit frequency f is determined based on this map. _maxThis frequency covers, for example, 99% of the total energy of multi-source loads.

[0127] In some scenarios, assuming that the wind load energy is mainly distributed in the 0-0.5 Hz range, the wave load energy is distributed in the 0.05-0.3 Hz range, and the seismic load energy is distributed in the 0.1-25 Hz range, then the f under this load condition... _max The frequency is 25 Hz. The initial screening rule is to retain all inherent frequencies f. _n Satisfy f _n ≤α×f _max The modes are initially retained, where α is the safety factor, for example, 1.5, which means that all modes with natural frequencies below 37.5 Hz are initially retained.

[0128] Based on the initially selected set of modes, a modal energy participation criterion is applied for further screening. This criterion is used to quantify the energy contribution of each mode in the actual structural response, eliminating modes that, although within the frequency range, are hardly excited.

[0129] In practice, the generated full-order response snapshot matrix is ​​projected onto the initially selected modal vectors to obtain the coordinate response time series q for each mode. _j (t). Next, the response energy index for each mode is calculated, such as the mean square value E of its coordinates. _j =mean(q _j (t) 2 );

[0130] In the formula, mean() represents taking the average value.

[0131] Based on this, all candidate modes are sorted according to their energy contribution, and the cumulative energy contribution percentage is calculated. Modes that cause the cumulative energy contribution to reach a preset threshold, such as 99.9%, are retained.

[0132] For example, if calculations show that the cumulative energy contribution of the first 48 modes has reached 99.9%, then under this criterion, only the 48 modes are retained.

[0133] As an alternative, a cross-coupling sensitivity criterion can be applied for supplementary screening. This criterion is used to identify and force the retention of modes that are important for reflecting the nonlinear coupling effect of multi-source loads. Modes may not be active under a single load and are therefore easily ignored by the aforementioned energy criterion.

[0134] In some embodiments, comparative simulations are required. For example, three full-order simulations are performed: one under wind load only, one under wave load only, and one under combined wind and wave action. For each mode, the response energy E under combined excitation is calculated. _j_coupled And the sum of the response energies under each single-source excitation (E_j_wind +E _j_wave ).

[0135] In some scenarios, if the ratio C of a certain mode is found... _j =E _j_coupled / (E _j_wind +E _j_wave If the value is greater than 1, for example, exceeding the threshold of 1.2, it indicates that the mode has been amplified by 20% under coupling effect and should be identified as a cross-coupling sensitive mode and forcibly retained, regardless of how small its energy contribution is under a single load.

[0136] Based on this, the offline-determined set of retained modes is the union of modes selected by the energy participation criterion and modes identified by the cross-coupling sensitivity criterion. It should be understood that this set serves as the initial mode basis for real-time hybrid loading. During the loading process, this mode set can be dynamically expanded based on online error monitoring results to achieve adaptive behavior.

[0137] In some embodiments, it is illustrated how a reduced-order numerical substructure model possesses online self-calibration and adaptive evolution capabilities to address the performance evolution of the physical substructure during long-term experiments. Accordingly, this embodiment includes:

[0138] Step 501: The reduced-order numerical substructure model has online self-correction capability, which is used to continuously track the changes in physical substructure performance caused by accumulated damage during the execution of real-time hybrid loading.

[0139] Among these, performance changes include at least the degradation of the equivalent stiffness of the physical substructure.

[0140] In this step, the initially generated reduced-order numerical substructure model is established based on the initial state of the physical substructure being intact, and belongs to the static model.

[0141] In some scenarios, after prolonged fatigue loading or ultimate load application, the physical substructure will experience cumulative damage, such as the plastic development of the material, loosening of bolts in connection areas, and the initiation and propagation of microcracks in welded areas. This damage will directly lead to changes in the mechanical properties of the physical substructure, such as degradation of its equivalent stiffness and alteration of its damping characteristics.

[0142] In this situation, if the numerical model cannot be updated accordingly, the mismatch between the physical and numerical systems will gradually increase, potentially leading to distorted experimental results or even system instability. Therefore, endowing the reduced-order model with online self-calibration capabilities, enabling it to continuously track the state changes of the physical specimen like a learning system, is beneficial for achieving long-term hybrid simulations.

[0143] In some embodiments, the physical parameters identified by the Extended Kalman Filter (EKF) are primarily used for the following purposes:

[0144] Update the stiffness or damping terms in the reduced-order dynamic equations that are related to the physical-numerical interface coupling to make the reduced-order model more accurate in predicting the interface response.

[0145] This serves as the driving information for model mismatch assessment and adaptive modal order adjustment.

[0146] Correspondingly, the non-intrusive framework calls the original solver and only simulates the physical processes inside the numerical substructure. The state changes of the physical substructure are directly obtained through sensor measurements and then applied as boundary conditions. Therefore, the original solver does not need to be modified based on the EKF identification results.

[0147] Step 502, Execute online self-correction capability, which can be achieved through a state estimation algorithm;

[0148] The state estimation algorithm combines the coordinates of the reduced-order numerical substructure model and the physical parameters to be identified to form an augmented state vector. In each computation time step, the actual interface response of the physical substructure is used as the observation value to predict and update the augmented state vector, thereby realizing online identification and correction of the physical parameters.

[0149] In this embodiment, the Extended Kalman Filter (EKF) algorithm can be used to achieve online self-calibration. The specific implementation process is as follows:

[0150] Accordingly, an augmented state vector X is defined, which not only contains the coordinates q representing the dynamic response of the model, but also the important physical parameters p to be identified, which can characterize the performance changes of the physical substructures. This augmented state vector can be expressed as:

[0151] X _k =[q T _k ,p T _k ] T ;

[0152] In the formula, the physical parameter p can be one or more scalars. For example, to track the stiffness degradation of the tower foundation, p can be set as the equivalent stiffness coefficient describing the moment-rotation relationship at the tower base connection. Here, the subscript k represents the discrete time step k, and the superscript T indicates the vector transpose operation.

[0153] Within each computation time step, the EKF algorithm performs a prediction step and an update step.

[0154] For example, in the prediction step, the algorithm predicts the prior state vector at the next moment based on the dynamic equations and parameter evolution model of the reduced-order model. The evolution of parameter p is typically modeled as a random walk process, meaning it is assumed to remain constant over a short period, but allows for small perturbations.

[0155] In another example, the update step feeds back real-world information to the model. In this step, the actual interface response z of the physical substructure at the physics-numerical interface is... _k For example, the interface reaction force vector measured by a force sensor is used as an observation in the EKF algorithm.

[0156] Simultaneously, using the reduced-order model and the current prior state estimate, the interface reaction vector h(X) predicted by the model is calculated. _k The difference between the model predictions and the actual observed values, i.e., ε. _k =z _k -h(X _k This is known as information or prediction error.

[0157] In summary, the EKF algorithm uses this information and combines it with the system's covariance information to calculate a better Kalman gain, which is then used to correct the prior state vector to obtain the posterior state estimate.

[0158] In the update step of the EKF algorithm, it is necessary to calculate the Jacobian matrix of the observation function h with respect to the augmented state vector X. Since the observation is an interface reaction force, it can be expressed by coordinates q and physical parameters p through the stiffness relations and dynamic equations of the reduced-order model. Therefore, the Jacobian matrix can be decomposed as follows:

[0159] The partial derivatives with respect to coordinates can be obtained analytically using the chain rule of the reduced stiffness matrix and the time integral scheme;

[0160] The partial derivatives of the physical parameters can be analytically derived from the way parameter p participates in the restoring force term or approximated numerically through finite difference.

[0161] In one implementation, the Unscented Kalman Filter (UKF) can be used instead of the EKF to avoid explicitly calculating the Jacobian matrix.

[0162] Through this process, not only is the estimated value of the coordinate q corrected to be closer to the actual dynamics, but the physical parameter p to be identified is also continuously and incrementally updated in the direction that can minimize the prediction error, thus realizing online tracking and identification of physical stiffness degradation.

[0163] Step 503: Quantify the model mismatch between the reduced-order numerical substructure model and the physical substructure using the prediction error sequence between the model predictions and the observed values ​​output by the state estimation algorithm.

[0164] Furthermore, when the model mismatch continues to exceed a preset threshold, it is determined that the current retained modal order is insufficient to characterize the state changes of the physical substructure, triggering an adaptive increase in the retained modal order of the reduced-order numerical substructure model.

[0165] When a physical substructure undergoes significant nonlinear behavior or local damage, its dynamic characteristics may undergo a qualitative change. Updating only a few parameters is no longer sufficient to describe its new behavioral patterns. In this case, it is necessary to increase the degrees of freedom of the model, i.e., to evolve the model structure itself.

[0166] In this case, a real-time quantitative index of model mismatch is constructed using the EKF algorithm byproduct, namely the innovation sequence. Because the innovation ε... _k The dimensions and covariance of the expression change over time, and directly using its amplitude is not stable. Therefore, the normalized squared innovation can be used as the mismatch index d. 2 _k The calculation formula is as follows:

[0167] d 2 _k =ε T _k ×S -1 _k ×ε _k ;

[0168] In the formula, S _k This is the information covariance matrix calculated by the EKF algorithm. This index is a dimensionless scalar; under ideal model-fitting conditions, its sequence follows a chi-square distribution with degrees of freedom equal to the dimension of the observations. S -1 _k The superscript -1 in the matrix S indicates that the superscript -1 indicates the superscript -1 in the matrix S. _k The inverse operation of ε. T _k The superscript T in ε indicates that the superscript T is relative to ε. _k The transpose of d. 2 _k The superscript 2 in the text indicates that the value of d is 2. _k Perform the squaring operation.

[0169] Based on this, a preset threshold can be set. For example, the threshold can be set according to the 95% confidence limit of the chi-square distribution. During real-time loading, the mismatch index is continuously monitored. If the index continuously and significantly exceeds the preset threshold within a certain time window, the system determines that the current model mismatch is no longer due to parameter drift, but rather to insufficient model structure.

[0170] Once this judgment is made, the system will automatically trigger a model update event, which will activate the alternate modes that were identified in the offline phase but not activated in the initial model, add them to the projection basis matrix, and adaptively increase the degrees of freedom of the reduced-order model so that it can capture and characterize the new and more complex dynamic behaviors of the physical substructure.

[0171] In other words, based on parameter self-calibration, the adaptive evolution of the model structure was further realized.

[0172] In other scenarios, when performing adaptive modal order increases, the initial conditions for the newly activated modes can be determined in the following way:

[0173] Accordingly, the full-order displacement vector u reconstructed at the current moment using the existing reduced-order basis is... current Projected onto the newly added modal vector φ new The initial coordinate value q of the newly activated mode is obtained. new =φ T new ×M×u current , where M is the mass matrix; similarly, the initial velocity values ​​of the new activated modes are calculated using the full-order velocity reconstruction vector at the current moment.

[0174] For the augmented state vector in the state estimation algorithm, the covariance matrix of its extended dimension is initialized as a large diagonal matrix, reflecting the high uncertainty of the new activated mode state, so that the filter can quickly converge to the correct state through the observation data.

[0175] By using an initialization method, the continuity of time integrals during modal expansion can be maintained, avoiding closed-loop transient shocks caused by unreasonable initial conditions.

[0176] Step 504 further includes at least one of the following measures to enhance the robustness of online self-calibration:

[0177] Before inputting the actual interface response as an observation into the state estimation algorithm, sensor signal preprocessing is performed to remove outliers.

[0178] Pre-defined physical boundary constraints are applied to the important physical parameters identified by the state estimation algorithm to prevent non-physical results from occurring during parameter updates.

[0179] Optionally, preprocessing can be performed on the observed signals acquired by the sensors. In actual operating conditions, factors such as environmental electromagnetic interference and transient sensor failures may cause the acquired interface response signals to contain outliers or anomalous jumps. If outliers are directly input into the EKF algorithm, they will contaminate the parameter estimation results. Therefore, before inputting the signal into the algorithm, a sliding window filtering strategy, such as median filtering, is applied to smooth the signal and remove individual impulse interference.

[0180] Optionally, physical boundary constraints can be imposed on the parameter identification results. For example, the identified equivalent stiffness coefficients cannot be physically negative. After each update step of the EKF, the newly obtained parameter estimates are checked.

[0181] If it is less than a preset minimum positive value, such as 0.01 times its initial value, it is forcibly corrected to that minimum positive value. This constraint prevents the identification process from diverging into non-physical regions due to numerical issues, ensuring the convergence and stability of the algorithm.

[0182] For systems with high nonlinearity, the state estimation algorithm can also employ the Unscented Kalman Filter (UKF). The UKF transmits the mean and covariance of the state distribution through deterministically sampled sigma points, avoiding the linearization and differentiation of the nonlinear function required by the EKF, thus improving estimation accuracy and adaptability to strong nonlinearities.

[0183] As an optional implementation method to ensure that the non-intrusive order reduction framework meets stringent real-time requirements, the following are some commonly used methods:

[0184] Step 601: Based on the full-order response snapshot obtained through pre-simulation, offline pre-compute the parameterized model of the reduced-order Jacobian matrix with respect to the state variables;

[0185] In real-time hybrid loading, the parameterized model is queried or interpolated to obtain the current reduced-order Jacobian matrix, avoiding online calculation and projection of the full-order Jacobian matrix and meeting real-time requirements.

[0186] In a closed-loop system with real-time hybrid loading, each computation step needs to be completed within a short time budget, such as 1-2 milliseconds. When using implicit time integration algorithms with good numerical stability, such as the Newmark-Beta method, it is necessary to calculate the tangent stiffness matrix of the system, i.e., the Jacobian matrix, and solve the linear equations at each time step or iteration.

[0187] For a non-intrusive framework, obtaining the reduced-order Jacobian matrix J _reduced The Jacobian matrix J needs to be calculated or assembled in the full-order space first. _full Then through projection operation J _reduced =Φ T ×J _full ×Φ reduces its dimensionality.

[0188] The above operations, especially J _full The computation of reduced-order Jacobian matrices is quite time-consuming. To address this issue, this invention proposes a solution that shifts the computational burden from online to offline methods, namely, constructing a computationally inexpensive parameterized model or surrogate model for the reduced-order Jacobian matrix. This process consists of two stages: offline construction and online invocation.

[0189] During the offline build phase, full-order response snapshot data is utilized. For each time point in the snapshot data, corresponding to the instantaneous state of the system, the full-order Jacobian matrix J can be calculated for that state. _full The projection yields its corresponding reduced-order Jacobian matrix J._reduced .

[0190] Furthermore, by traversing all snapshot points, the corresponding dataset can be obtained, where each data point contains a state vector and its corresponding reduced-order Jacobian matrix.

[0191] Here, state variables are physical quantities used to parameterize the Jacobian matrix, which may include coordinates, velocities, and physical-numerical interface displacements that can affect the nonlinear behavior of the system in the reduced-order model.

[0192] Next, a parameterized model is constructed based on the dataset. Data clustering methods can be used. For example, each reduced-order Jacobian matrix can be treated as a data point in a high-dimensional space, and the K-means clustering algorithm can be used to divide the matrix into N categories. Furthermore, the center matrix of each category is calculated. The resulting parameterized model consists of N center matrices and their corresponding state variable space regions.

[0193] The parameterized model is not limited to the above methods and can also take other forms, such as a response surface model based on multinomial regression, or a small, trained neural network model whose input is the state variable and whose output is an approximation of the reduced-order Jacobian matrix.

[0194] During the online invocation phase, i.e., in real-time hybrid loading, the current state variable is retrieved whenever a Jacobian matrix is ​​needed at each time step. Further, the offline-built parameterized model is queried or interpolated. If a clustering model is used, this invocation process is a relatively fast query process; the system determines which cluster the current state variable belongs to and directly uses the center matrix of that region as an approximation of the reduced-order Jacobian matrix for the current time step. The computational cost of this query operation is lower than performing a complete matrix calculation and projection online, and can be completed in microseconds.

[0195] It should be understood that by constructing offline and querying online, this invention preprocesses time-consuming computational tasks, ensuring that the online computation step size meets real-time requirements and providing a guarantee for the operation of the entire hybrid simulation system.

[0196] In other scenarios, the implicit Newmark-β method can be used as the time integration algorithm in real-time closed loop because it has unconditional numerical stability and is adaptable to numerical rigidity problems that may occur after model order reduction.

[0197] In some other scenarios, if the condition number of the reduced-order model is good and the real-time step size is small enough, the computationally more efficient fourth-order Runge-Kutta method can be used as an alternative.

[0198] According to one aspect of this application, in a real-time hybrid loading closed loop, a multi-mechanism fusion adaptive delay compensation scheme can be adopted to address the response delay problem of the loading system. Specifically:

[0199] Prior to this, offline dynamic characteristic calibration of the multi-degree-of-freedom loading system is required. This process is used to obtain the dynamic model of each loading actuator. A command signal with a known waveform, such as a swept-frequency sine wave or a broadband white noise signal, can be applied to an individual actuator, and the actual displacement response of the actuator can be measured synchronously using a high-precision displacement sensor, such as a linear variable differential transformer (LVDT).

[0200] Accordingly, by performing system identification and analysis on the input command signal and the output response signal, a transfer function model can be established for the actuator. This model describes the dynamic process of the actuator from receiving a command to generating actual displacement, including its inherent delay, bandwidth, and gain characteristics. This transfer function model can provide a foundation for the subsequent implementation of an inverse transfer function feedforward compensation mechanism.

[0201] Step 701, in real-time hybrid loading, before the loading system applies the interface target instruction, it also includes preprocessing the interface target instruction through an adaptive delay compensator;

[0202] Among them, the adaptive delay compensator is used to identify and quantify the equivalent delay caused by physical inertia of the loading system online based on the comparison between the target command of the interface and the measured response of the loading system. It dynamically adjusts its internal compensation parameters according to the real-time identified equivalent delay, actively compensates for the phase lag introduced by the loading system, and suppresses the risk of instability of the closed-loop system.

[0203] Correspondingly, in any physically loaded system, due to the physical inertia of its actuators, the response cannot instantly follow the command, and there will inevitably be a delay. Actuators include servo valves and hydraulic cylinders.

[0204] In a closed-loop system of real-time hybrid simulation, the phase lag introduced by the response delay is equivalent to introducing negative damping into the coupled system. When its magnitude exceeds the physical damping of the system, it can easily lead to closed-loop instability, experimental divergence, and affect experimental safety.

[0205] Therefore, proactive compensation is implemented for the loading system's latency. The adaptive latency compensator is a dynamic compensation unit with online learning and adjustment capabilities. During the experiment, it continuously preprocesses the target instructions on the interface, generating advanced and corrected instructions. This ensures that the loading system's physical response synchronously tracks the original, unprocessed target instructions after receiving the preprocessed instructions.

[0206] Step 702, the adaptive delay compensator incorporates at least one of the following compensation mechanisms:

[0207] By using the interface target instruction sequence from multiple past time steps, an extrapolation polynomial is fitted to predict instruction values ​​that are ahead of the current time as compensated instructions.

[0208] The system compares the target command on the interface with the measured response of the actuator of the loading system obtained by the sensor in real time, calculates the equivalent delay and gain attenuation at the current loading frequency and amplitude online, and corrects the extrapolation duration and amplitude of the command based on the equivalent delay and gain attenuation.

[0209] The interface target command is filtered by the inverse function of the actuator transfer function obtained through offline calibration, and the dynamic characteristics of the actuator are pre-canceled in the frequency domain.

[0210] At any given computation time step, the mechanism utilizes the target instruction values ​​from the most recent historical moments stored in memory, for example. Based on these historical data points, the system fits a low-order polynomial. Further, it extrapolates the time length τ forward over this polynomial. _est Calculate the preprocessed instruction value u _comp (t _k )=p(t _k +τ _est ). τ here _est This is the estimated delay time of the system, and p() corresponds to the fitting polynomial.

[0211] It should be understood that this mechanism is computationally simple and resource-efficient, but its compensation effect depends on the delay τ. _est The accuracy of the estimate.

[0212] Furthermore, the online delay and gain estimation compensation mechanism can be used to provide dynamically updated delay estimates τ. _est In this mechanism, the equivalent delay of the actuator varies with the frequency and amplitude of the command signal. Therefore, during the loading process, the mechanism compares in real time the original target command sequence sent to the actuator with the actual actuator response sequence measured by sensors.

[0213] Optionally, this comparison can be implemented using online sliding window cross-correlation analysis. That is, within a sliding window of a set length, the time offset that maximizes the cross-correlation function is found, thus obtaining the equivalent delay τ within the current window. _eq Based on the comparison of the energy or amplitude of the two signals within the window, the equivalent gain attenuation γ is obtained. _eq The equivalent delay and equivalent gain attenuation obtained online are used in real time to update the compensator parameters, dynamically adjust the extrapolation duration, and perform gain correction on the compensation command.

[0214] For example, the extrapolation duration is dynamically set to τ. _eq The calculated compensation command is multiplied by the gain correction factor 1 / γ._eq .

[0215] Building upon this, the inverse transfer function feedforward compensation mechanism can be considered as an alternative, working in conjunction with the aforementioned mechanisms. It utilizes the actuator transfer function model H(s). When processing the target command, this mechanism calculates the inverse model H of the transfer function. _inv The formula (s) = 1 / H(s) is implemented as a digital filter. After time-domain extrapolation compensation, the target command is further filtered by this inverse model filter. This is used to pre-cancele the dynamic characteristics of the actuator itself.

[0216] It should be understood that when a command, after being filtered by the inverse model, is input to the actuator, its output will reproduce the original command. In practical applications, a regularized approximate inverse model that is effective within the target frequency band is used, along with a low-pass filter, to avoid excessive amplification of high-frequency noise.

[0217] According to another aspect of this application, an advanced loading control strategy may be optionally employed during real-time hybrid loading. This includes:

[0218] Step 801, during real-time hybrid loading, a frequency-dependent force-displacement hybrid control strategy is adopted. This strategy uses a preset crossover frequency as the boundary, i.e.:

[0219] In the low-frequency band below the crossover frequency, a displacement control mode is adopted, in which the displacement command calculated by the reduced-order numerical substructure model is applied to the physical substructure, and the measured interface force is used as feedback input to the reduced-order numerical substructure model.

[0220] In the high-frequency band above the crossover frequency, a force control mode is adopted, in which the force command calculated by the reduced-order numerical substructure model is applied to the physical substructure, and the measured interface displacement is used as feedback input to the reduced-order numerical substructure model.

[0221] At the crossover frequency, a digital complementary filter is used to achieve a smooth transition between the two control modes.

[0222] For wind turbine towers, a single control mode cannot meet the requirements of wideband hybrid simulation. Therefore, the proposed frequency-dependent force-displacement hybrid control strategy combines the advantages of both modes, automatically switching the dominant control mode across different frequency bands to achieve globally superior control performance.

[0223] The strategy uses a preset crossover frequency to delineate the dominant frequency bands of the two control modes. Optionally, this crossover frequency setting is related to the dynamic performance of the loaded actuator.

[0224] One exemplary setting method is to set it to one-third of the actuator's closed-loop bandwidth in force control mode. For example, if the force control bandwidth of a large hydraulic servo actuator is 15 Hz, the crossover frequency can be set to 5 Hz.

[0225] In the low-frequency band below the crossover frequency, such as in the 0 to 5 Hz range, the system primarily employs a displacement control mode. Specifically, at each computation time step, the reduced-order numerical substructure model calculates the target interface displacement command for the next time step. This displacement command is sent to the controller of the loading system, driving the actuators to move the physical substructure to the target displacement. Force sensors positioned at the connection between the actuators and the physical substructure measure the actual interface force between them in real time. This measured interface force, as the actual boundary condition, is fed back to the reduced-order numerical substructure model for solving the problem in the next time step.

[0226] In this frequency band, direct control of displacement can ensure the matching of boundary conditions between the numerical model and the physical specimen, which helps to achieve the accuracy and stability of quasi-static loading.

[0227] In other scenarios, at higher frequency bands above the crossover frequency, such as above 5 Hz, the system seamlessly switches to force control mode as the dominant mode. Specifically, at each computation time step, the reduced-order numerical substructure model calculates the interface target force command for the next time step; this force command is sent to the controller of the loading system, driving the actuators to apply the target force to the physical substructure; displacement sensors placed at the interface measure the actual displacement response of the physical substructure under this force in real time; this measured interface displacement, as the actual boundary condition, is fed back to the reduced-order numerical substructure model for solving the next time step.

[0228] In this frequency band, direct control of force can avoid the conflict between the actuator and the structural inertial force, realize the reproduction of high-frequency dynamic force, and ensure the stability and fidelity of dynamic loading.

[0229] Alternatively, a digital complementary filter technique can be employed. The reduced-order model simultaneously calculates the target displacement command and the target force command at each time step. A low-pass filter is applied to the target displacement command, while a complementary high-pass filter is applied to the target force command. The crossover frequency of the two filters is set to a preset crossover frequency, such as 5 Hz. The command sent to the loading system controller is the superposition of the two filtered command signals.

[0230] In this way, at low frequencies, the weight of displacement commands is close to 1, while the weight of force commands is close to 0; as the frequency increases and crosses the cross frequency, the weight of displacement commands smoothly decays to 0, while the weight of force commands smoothly increases to 1.

[0231] The above implementation scheme based on complementary filters ensures that the transfer of control between the two modes is continuous and smooth, avoiding transient oscillations in the system that may be caused by hard switching, and guaranteeing the stability and reliability of the entire hybrid loading process.

[0232] On the other hand, this paper describes how to reproduce the transient coupling effect between the servo control system and structural dynamics through physical-numerical coupling. See below for details:

[0233] Step 901, the real-time simulation of the transient response of the servo control system is achieved through a closed loop that includes physical feedback, in which:

[0234] The state changes of the physical substructure are fed back to the reduced-order numerical substructure model via the actual interface response.

[0235] The change in state causes a change in the controller input signal in the dynamic model of the servo control system;

[0236] The dynamic model of the servo control system generates a response action based on the changed input signal, and the response action causes changes in aerodynamic or mechanical loads.

[0237] The changing loads are solved by the reduced-order numerical substructure model and then applied back to the physical substructure via the interface target command, forming a complete servo-structure coupling effect.

[0238] In other embodiments, by using real-time closed-loop, the real-time dynamic response of the physical substructure is fed back to the reduced-order numerical substructure model. The reduced-order numerical substructure model adjusts the control commands of its internal servo control system in real time according to the feedback, and applies the servo load changes caused by the adjustment of the control commands back to the physical substructure. This can reproduce the real transient coupling behavior between the servo and the structure under transient conditions such as emergency shutdown.

[0239] Taking a typical emergency shutdown scenario as an example, suppose that during the test, extreme gust loads act on the system, causing the bottom connection area of ​​the tower column, which is a physical substructure, to enter the plastic deformation stage. At this time, the transient coupling information flow between the servo and the structure will be transmitted according to the following chain:

[0240] Correspondingly, changes in the physical state and feedback occur. These changes in the physical state are immediately measured by sensors placed at the physical-numerical interface. For example, the increase in interfacial reaction force measured by a force sensor begins to lag behind the increase in displacement, or the interfacial displacement measured by a displacement sensor becomes larger under the same load increment.

[0241] Furthermore, the model state is solved and triggered. The actual interface response measured by the sensors is used as a boundary condition and input into the reduced-order numerical substructure model in the next millisecond time step. Because the model senses a softer physical boundary, the dynamic state of its solved numerical substructure will change. For example, the acceleration at the top of the tower or the rotor speed may experience unexpected and severe overshoot, exceeding the servo system safety threshold.

[0242] Based on this, the actions are transmitted to the servo control logic. State variables such as overshoot tower top acceleration or rotor speed are input in real-time into the servo control system dynamics model contained within the reduced-order model. This servo model, for example, the control logic module implemented in the Simulink modeling algorithm development environment, will immediately determine that the system has entered an emergency state and trigger its preset emergency shutdown protection procedure.

[0243] For example, one emergency stop action is that the servo system issues a command requiring the pitch angle of all three blades to be rapidly adjusted from the current operating angle to a 90-degree feathering position at maximum speed.

[0244] Furthermore, abrupt changes and applications of servo loads occur. Pitch control is executed within the numerical model, causing a sudden and drastic reduction in the blade angle of attack, resulting in a sharp drop in the aerodynamic thrust and torque calculated by the aerodynamic model acting on the entire rotor. This servo load is immediately reflected in the solution results of the interface target command at the next time step. For example, the thrust command applied to the tower top will decrease. The new, abrupt command is then applied back to the physical tower, which has already entered plastic deformation, via the loading system.

[0245] Through the above closed-loop feedback chain, the highly nonlinear and instantaneous strong coupling effect between the structural dynamics characteristics and the servo control system can be captured in transient events such as emergency shutdown.

[0246] Furthermore, after completing the loading of all preset target operating conditions, the system records data from the physical substructure sensors and the numerical substructure calculation process. Subsequent data processing and analysis may optionally include the following steps:

[0247] Accordingly, interface compatibility verification is performed. This is used to verify the quality of the hybrid loading test and quantitatively evaluate the tracking accuracy between the target command and the actual response at the physical-numerical interface throughout the entire test. Specifically, the time series of the interface target command recorded by the real-time computing platform and sent to the loading system is compared with the time series of the actual interface response acquired by sensors. For example, the target displacement command at the top of the tower is compared with the measured displacement, and the target reaction force command at the lower interface is compared with the measured reaction force. Evaluation metrics may include the root mean square error, peak error, and cross-correlation coefficient between the two.

[0248] For example, in a mixed-load test, the cross-correlation coefficient of displacement tracking should be higher than 0.99, while the error of force tracking should be within an acceptable range, thus verifying the effectiveness of the delay compensator and control strategy.

[0249] Furthermore, the overall dynamic response of the entire machine is reconstructed. After real-time hybrid loading, the measured response of the physical substructure and the reduced-order coordinate response of the numerical substructure are obtained, and the two sets of data are fused into the overall dynamic response corresponding to the full-order model degrees of freedom.

[0250] In practice, the response of the physical substructure, such as strain and acceleration at any location within it, is directly obtained from its sensor measurements. For the numerical substructure, such as blades and nacelles, the full-field response is obtained by reconstructing the response by left-multiplying the reduced-order coordinate time series by the projection basis matrix.

[0251] Based on this, by splicing the measured physical substructure response with the reconstructed numerical substructure response, a dataset that can be used for overall structural strength verification, fatigue life assessment, and limit state analysis can be obtained.

[0252] Furthermore, a quantitative analysis of the multi-source load coupling effect is conducted. This invention can capture the nonlinear coupling effect between multi-source loads. To quantify this effect, a coupling amplification factor index can be optionally introduced. The calculation of this index requires comparative experiments.

[0253] For example, to analyze the wind-wave coupling effect, three separate experiments are required:

[0254] One case involved wind load alone, another involved wave load alone, and the third involved a combination of wind and wave loads.

[0255] Furthermore, for a specific response quantity, such as the bending moment M at the bottom of the tower, the formula for calculating the coupling amplification factor is:

[0256] Coupling factor _Factor =M _peak_coupled / (M _peak_wind_only +M _peak_wave_only );

[0257] Among them, M _peak_coupled It is the peak bending moment at the base of the tower measured under combined loading conditions, while M _peak_wind_only and M _peak_wave_only These are the peak values ​​measured under two single-source loading conditions, respectively.

[0258] If the coefficient value is greater than 1.0, for example, 1.25, it indicates that the wind-wave coupling effect increases the peak bending moment at the tower base by 25%, a risk that would be ignored by the linear superposition method. This indicator can provide a basis for the fatigue resistance and ultimate bearing capacity design of wind turbine structures.

[0259] Furthermore, a three-level comparative verification can be optionally performed.

[0260] The comparison and verification with pure numerical simulation were conducted. The response results obtained by the hybrid simulation method of this invention were compared with the results obtained by pure software simulation using a numerical model of the wind turbine system without order reduction. The high degree of agreement between the two results proves that the order reduction framework, delay compensation, and load control strategy adopted in this invention do not introduce large system errors.

[0261] The method of this invention is compared with traditional experimental methods. The results are compared with those of traditional pseudo-dynamic tests or shaking table tests using open-loop loading based on predefined load spectra. The comparison reveals that the method of this invention can capture physical phenomena that traditional methods might miss, such as the strong servo-structure coupling effect during emergency shutdown.

[0262] The method was compared with actual field measurement data. When conditions permitted, the results obtained by the method under specific sea conditions were compared with data obtained from on-site monitoring of real offshore wind turbines under similar sea conditions. The good agreement with the on-site data demonstrates the reproducibility of the method.

[0263] In some embodiments, calling the original solver means, under a given full-order displacement and velocity state, invoking the force calculation functions of each physics module within the solver to evaluate the internal and external forces under the current state, rather than driving the solver to perform its complete implicit time-step solution, including Newton iterations.

[0264] For the aerodynamic module, this call performs forward kinematics calculation of blade element momentum based on the given blade position and wind speed.

[0265] For the servo module, control logic decisions are made based on given sensor inputs.

[0266] For the hydrodynamic module, the Morison equation evaluation can be performed based on a given location.

[0267] The computational cost of a single force assessment in each module is less than that of multiple Newton iterations and solving large linear equation systems in implicit time steps. The time consumption of a single call is usually in the sub-millisecond range, which can meet the step size requirements of real-time calculation.

[0268] For example, the method of the present invention may also be:

[0269] The tower section containing the basic transition section is classified as a physical substructure to improve the simulation fidelity of the damaged area;

[0270] A non-intrusive projection framework is used to generate a reduced-order numerical substructure model. This model has online self-correction capability based on state estimation to track physical damage, and real-time performance is ensured by offline pre-computing of the Jacobian matrix.

[0271] Furthermore, in closed-loop loading with adaptive delay compensation, the physical response is fed back to the numerical model in real time to reproduce the servo-structure transient coupling effect.

[0272] For example, regarding the setting of the covariance matrix in the EKF algorithm, the method for setting the process noise covariance matrix Q and the observation noise covariance matrix R in the state estimation algorithm is as follows:

[0273] The observation noise covariance matrix R is initialized based on the nominal accuracy of the sensor used. For example, if the force sensor's measurement accuracy is 0.5% of full scale, then the corresponding diagonal element in R is set to the square of the sensor's range multiplied by 0.5%. In the process noise covariance matrix Q, the diagonal element corresponding to coordinate q can be set to a small value, such as the square of 1% of the expected response amplitude; the diagonal element corresponding to the parameter p to be identified reflects the expected rate of change of the parameter, for example, set to the square of 0.1% of the initial stiffness value, allowing for small adjustments to the parameter at each step. In practical applications, the values ​​of Q and R can be optimized through offline pre-simulation.

[0274] For example, the following security protection measures are also included during the execution of the real-time closed loop:

[0275] When the real-time computing platform detects that the single-step calculation time exceeds the preset step size budget, the system starts the degradation mode, such as temporarily freezing the interface target instructions of the previous step to maintain loading continuity, and resuming normal calculation in subsequent time steps.

[0276] When an abnormal drop in hydraulic pressure of the loading system is detected or the actuator displacement exceeds the safety limit, the system automatically triggers an emergency unloading procedure, and all actuators slowly retract to the zero position to protect the physical specimen and equipment.

[0277] When communication between the real-time computing platform and the load controller is interrupted for more than a preset timeout threshold (e.g., 5 milliseconds), the load controller independently performs a safety action to maintain the current position until communication is restored or the operator intervenes manually.

[0278] This invention introduces a reduced-order model with online self-correction capabilities. By embedding a state estimation algorithm, the measured response of the physical substructure is used as an observation, and physical parameters, including equivalent stiffness, are continuously identified and corrected online. This allows the numerical model to evolve dynamically, capturing the performance degradation process of the physical specimen in real time, thus solving the problem of simulation distortion.

[0279] Furthermore, a real-time closed-loop loading system was also constructed. The dynamic response of the physical substructure is fed back to the numerical model in real time, directly driving the logical decisions of its internal servo control system. Changes in commands generated by the control system (such as emergency pitch) will instantly change the load acting on the physical substructure, forming a high-speed, bidirectional coupled information flow, reproducing the coupled physical process under extreme conditions.

[0280] Accordingly, a non-intrusive projection framework is adopted, treating the original solver as a black box, thus avoiding source code modification and improving the method's versatility. Furthermore, by pre-compiling its parameterized model with respect to state variables offline, online computation is transformed into efficient querying or interpolation, ensuring the system's hard real-time performance.

[0281] It should be noted that the various specific technical features described in the above embodiments can be combined in any suitable manner without contradiction. To avoid unnecessary repetition, the present invention will not describe the various possible combinations separately.

Claims

1. A fan tower multi-source load coupling simulation test method, characterized in that, include: A numerical model of the wind turbine system is established, and the numerical model of the wind turbine system is divided into a physical substructure to be tested and a numerical substructure to be simulated. The numerical substructure includes the dynamic model of the servo control system. Based on numerical substructure, a reduced-order numerical substructure model that is dynamically equivalent and can be computed in real time is generated. In the real-time closed loop, based on external multi-source load input, the following operations are repeated: the actual interface response of the physical substructure at the physical-numerical interface is measured; the actual interface response is used as a boundary condition input to the reduced-order numerical substructure model, and the interface target instruction for the next time step is calculated in real time; the interface target instruction is applied to the physical substructure via the loading system.

2. The method of claim 1, wherein, When establishing and dividing the numerical model of the wind turbine system, the scope of the physical substructure is defined as including the wind turbine tower body and the foundation transition section connected to the foundation. Through this division, the local nonlinear mechanical behavior caused by flange connection or grouting connection in the foundation transition section is included in the test of the physical substructure.

3. The method of claim 1, wherein, The reduced-order numerical substructure model is generated, specifically through a non-intrusive projection framework. This framework performs the following operations within each computation time step: Map the state vector in the reduced-order degree-of-freedom space at the current moment to the full-order degree-of-freedom space; The original numerical substructure solver is invoked to solve the nonlinear terms, including aerodynamic forces, hydrodynamic forces, and servo control forces, in one step in the full-order degree-of-freedom space to obtain the full-order residual forces. Project the full-order residual force back into the reduced-order degree-of-freedom space and solve for the state vector at the next moment. The internal nonlinear computation mechanism is preserved without modifying the source code of the original numerical substructure solver on which the numerical substructure depends.

4. The method of claim 1, wherein, Generating a reduced-order numerical substructure model also includes: The mode order for preserving the internal degrees of freedom of the numerical substructure is determined by a multi-criteria adaptive criterion. Based on the determined preserved mode order, modal truncation and degree of freedom condensation are performed on the numerical substructure to obtain a reduced-order numerical substructure model. The multi-criteria adaptive criterion comprehensively evaluates the spectral characteristics of multi-source load inputs and the energy contribution of each mode under load excitation.

5. The method according to claim 4, characterized in that, The multi-criteria adaptive criterion specifically includes at least one of the following criteria: Based on the power spectral density analysis of multi-source load input, modes whose natural frequencies cover the main frequency distribution bands of load energy were selected. Based on the coordinate response energy of each mode under different load conditions, a set of modes whose cumulative energy contribution reaches a threshold is selected. It identifies and forces the retention of modes whose response under multi-source joint excitation is greater than the linear superposition value of the response under each single-source independent excitation, thus capturing the nonlinear coupling amplification effect between multi-source loads.

6. The method according to claim 1, characterized in that, The reduced-order numerical substructure model has online self-correction capability, continuously tracking the performance changes of the physical substructure caused by accumulated damage during real-time closed-loop execution; among which, the performance changes include at least the degradation of the equivalent stiffness of the physical substructure.

7. The method according to claim 6, characterized in that, The online self-calibration capability is achieved through a state estimation algorithm. The state estimation algorithm combines the coordinates of the reduced-order numerical substructure model and the physical parameters to be identified to form an augmented state vector. In each computation time step, the actual interface response of the physical substructure is used as the observation value to predict and update the augmented state vector.

8. The method according to claim 1, characterized in that, In the real-time closed loop, before the loading system applies the interface target command, it also includes preprocessing the interface target command through an adaptive delay compensator. The adaptive delay compensator is used to identify and quantify the equivalent delay of the loading system due to physical inertia online, and dynamically adjust its internal compensation parameters according to the real-time identified equivalent delay to compensate for the phase lag introduced by the loading system.

9. The method according to claim 2, characterized in that, In addition to the interface at the top of the tower and the bottom of the foundation transition section, the physical-numerical interface also includes at least one intermediate loading section set along the height direction of the physical substructure. The intermediate loading section defines at least one horizontal translational degree of freedom, and in real-time hybrid loading, distributed environmental loads such as wind pressure are applied to the physical substructure cylinder.

10. The method according to claim 1, characterized in that, The real-time simulation of the transient response of the servo control system, included in the numerical substructure, is achieved through a closed loop that incorporates physical feedback, in which: The state changes of the physical substructure are fed back to the reduced-order numerical substructure model via the actual interface response. The change in state causes a change in the controller input signal in the dynamic model of the servo control system; The dynamic model of the servo control system generates a response action based on the changed input signal, and the response action causes changes in aerodynamic or mechanical loads. The changing loads are solved by the reduced-order numerical substructure model and then applied back to the physical substructure via interface target commands, forming a servo-structure coupling effect.

Images