Real-time high-precision structural finite element analysis method based on reduced-order model

Abstract

The invention is suitable for the field of predictive life calculation of aeronautical equipment, and provides a real-time high-precision structural finite element analysis method based on a reduced-order model, which is used for carrying out life calculation under structural load monitoring based on life parts of an airframe. According to structural load monitoring, airborne real-time data flow, an after-flight parameter data file and test bed data collected by a sensor serve as input, import and'decoding ', and calculation of structural load monitoring is carried out through a maintenance twin modeling and maintenance decision calculation software platform driving operation model based on a reduced-order model. Carrying out finite element load monitoring data and residual life calculation on parts with life on the machine body; the method sequentially comprises the steps of data processing, stress-strain calculation, rain flow counting, material attribute calculation and stress / strain fatigue calculation.

Description

TECHNICAL FIELD

[0001] The application is suitable for the field of predictive life calculation of aviation equipment, and provides a real-time high-precision structure finite element analysis method based on a reduced-order model. BACKGROUND

[0002] The safe life is established on the basis of no initial defects of the structure, and the life only considers the crack formation life, and once the macroscopic detectable crack appears in the structure, the structure is considered to be damaged. The length of the safe life of the helicopter body structure is mainly determined according to the life of the key dangerous parts of the structure, and once the dangerous parts are damaged, a disastrous accident will be caused.

[0003] Unlike fixed-wing aircraft, the helicopter structure not only has structural damage caused by state switching (i.e. load state change) during use, but also has structural damage caused by state (no load state change). Therefore, the industry needs to use a monitoring (SUMS) module, which can calculate the total time of each flight state by using its flight state recognition calculation function, to solve the life prediction problem based on the flight state. When performing residual safe life prediction of the helicopter, the commonly used reference fatigue safe life is the fatigue safe life of the key parts under the test load spectrum, so when performing life management of the helicopter structure, it is necessary to convert the actual load-time history of the helicopter into damage under the design spectrum through equivalent conversion.

[0004] The predictive maintenance digital twin calculation model of the helicopter body structure is obtained by driving the digital twin model of the monitored structure to generate a prediction algorithm by using the collected helicopter flight state and structure health feature data, and performing state evaluation or fault diagnosis according to the obtained signals to provide the perfect state information or fault warning of the related structure. The predictive maintenance digital twin calculation model of the helicopter body structure can perform real-time and online global monitoring on the helicopter structure, reduce the downtime caused by periodic detection and maintenance, reduce the maintenance and support cost, prevent disastrous structural failure, and improve the reliability and safety of the helicopter. SUMMARY

[0005] In view of the above defects, the purpose of the present application is to provide a real-time high-precision structure finite element analysis method based on a reduced-order model, which aims to solve the problems proposed in the background art, and perform life calculation under structural load monitoring based on the body life parts; wherein the structural load monitoring takes the real-time data stream collected by the sensor, the post-flight flight parameter data file and the test bench data as input, imports, decodes, and drives the operation model by the maintenance twin modeling and maintenance decision calculation software platform based on the reduced-order model to perform structural load monitoring calculation, so as to perform finite element load monitoring data and residual life calculation on the body life parts. The method sequentially includes data processing, stress and strain calculation, rainflow counting, material property calculation and stress / strain fatigue calculation process.

[0006] Further, the method comprises the following specific steps:

[0007] S1, the database reads the current residual damage life of the service life part;

[0008] S2, read the recorded parameters from the flight parameter data, input the stress and strain digital twin model of the key part, calculate the stress and strain history of the surface node; process the stress and strain signals, interpolate the null values, and filter the stress / strain signals;

[0009] S3, rainflow counting; use rainflow counting method to decompose the stress history into several stress cycles, and then use linear cumulative damage method to calculate the contribution of each stress cycle to fatigue damage;

[0010] S4, rainflow counting stress cycle calculation, get load cycle data;

[0011] S5, use linear cumulative damage method to calculate the fatigue damage of flight task;

[0012] S6, return the residual damage life of the service life part, and then calculate the cumulative damage according to the corresponding key component load spectrum;

[0013] S7, stress fatigue life calculation.

[0014] Further, the step S2 specifically comprises the following steps:

[0015] S2.1, use flight parameters and key component load data to build stress and strain digital twin surrogate model; store the stress and strain of the key component surface node calculated by the stress and strain digital twin surrogate model in the array until all flight parameters are read, and finally get the stress and strain history of the key component; for unknown flight state parameters, use load identification twin model to get the real-time load of the part;

[0016] The load identification twin model establishes the finite element model of each key part, uses the full-order model calculation result to establish a set of orthogonal bases, calculates a set of standard orthogonal basis vectors that can best describe the physical characteristics of the full-order system, and truncates a certain number of characteristic vectors to realize the reduction of the full-order model, and get the reduced model of the strain and stress of the tail structure;

[0017] S2.2, process the stress and strain signals and convert the stress and strain;

[0018] S2.3, interpolation processing; complete the load data by interpolation method, and interpolate once between the adjacent two data points in time;

[0019] S2.4, after obtaining the measured stress data, limit and filter the stress data.

[0020] Further, the load identification twin model in step S2.1 is constructed in the following manner:

[0021] For the free surface, the stress components perpendicular to the surface direction are all zero; a coordinate system is selected, and the z-axis is perpendicular to the surface, so that the independent unknown stress components are reduced to three, as shown in the following formula:

[0022]

[0023] First, the finite element model of each key part is established. If it is a solid structure, a layer of shell element with a thickness of 0.00001 is covered on the surface of the body mesh to extract the principal stress strain and shear stress strain of the structure surface;

[0024] First, the full-order result is constructed into a snapshot matrix A i Numerical simulation results from a theoretical model, or experimental results from an actual physical process; wherein the value represents the numerical result of node b in the a-th snapshot;

[0025]

[0026] Where L is the number of snapshots selected, and N is the number of grid nodes; find a set of standard orthogonal basis vectors:

[0027] Φ={Φ1Φ2…Φ L}

[0028] The sampling data T i The error norm of its projection is equivalent to a constrained maximum value:

[0029]

[0030] Φ T Φ=I

[0031] The vector Φ is the POD canonical orthogonal basis, and the matrix Φ is the POD canonical orthogonal basis set, which is an N-row L-column matrix, N>>L, and satisfies the orthogonality;

[0032] The constrained maximum value problem is solved by the Lagrange multiplier method, and the eigenvalue problem of the following formula is used to replace the original optimization problem:

[0033] SS T x=λx

[0034] To calculate the above eigenvalue problem, the singular value decomposition method is used to obtain the optimal orthogonal basis set for the snapshot matrix S; therefore, the matrix S∈R N×L is singular value decomposed:

[0035] A = BSV T

[0036] where A ∈ R N×N and V ∈ R L×L is an orthogonal matrix, S = R L×L = diag{σ1, σ2,..., σ L} is a diagonal matrix, each element in S represents a singular value, arranged in descending order;

[0037] Since the number of grid nodes N of the numerical calculation model is much larger than the number of selected snapshot samples L, i.e. N » L, the dimension of matrix AA T is much larger than that of matrix A T A, but their non-zero eigenvalues are exactly the same, so the eigenvalues λ i and eigenvectors of the smaller matrix A T A are solved first.

[0038]

[0039] Then the eigenvectors φ T of the larger matrix SS i can be converted.

[0040]

[0041] Finally, a set of standard orthogonal basis vectors Φ that can best describe the physical characteristics of the full-order system is solved, which can express the numerical response of any node on the model through linear combination, i.e.

[0042]

[0043] α i is the modal coefficient of each basis φ i .

[0044] By truncating a certain number of eigenvectors to realize the reduction of the full-order model, according to the energy ratio principle, a sufficient number of the first r basis is selected for model reduction, then the total energy of the reduced subspace should be close to 1 compared with the energy of the full-order model, expressed as

[0045]

[0046] Taking γ = 99.99, the truncated optimal orthogonal basis can be expressed as

[0047]

[0048] Then the unknown quantity to be solved is

[0049]

[0050] First, according to the range of load DOE test, full factorial design, Latin hypercube sampling can be used to generate design points; assuming 25 groups of test, using full order results by row to construct snapshot matrix u i The numerical simulation results of 25 groups of theoretical models;

[0051] Among them, the first line is the displacement of working condition 1 The second line is the displacement of working condition 2 So on and so forth;

[0052]

[0053] Among them, 25 is the number of snapshots selected, m is the number of grid nodes, assuming the model has 58877 nodes, and u is the comprehensive displacement;

[0054] The matrix U is decomposed into three matrices by SVD:

[0055] [B, S, V T ] = SVD (A)

[0056] Solve the modal coefficients corresponding to 25 samples:

[0057]

[0058] Calculate the eigenvalue λ i , that is, the modal energy S i 2 , calculate the relationship I (r) between the number of selected orthonormal basis vectors and modal energy, get the result, when selecting 4 groups of orthonormal basis vectors, I (r) is greater than 99.99%;

[0059] Take the first 4 order eigenvalue and the corresponding modal coefficient α r for fitting the displacement field of the structure, estimate the displacement field calculation after completing once reduction needs 0.00245 seconds, that is, 400Hz;

[0060]

[0061] At the same time, through Kriging algorithm, 20 groups of load input and corresponding modal 20 groups of modal coefficients are learned, 5 groups are used as verification group, get the proxy model of any load input F and modal coefficient α r ;

[0062] Then the corresponding structure displacement under i input is;

[0063]

[0064] Similarly, the snapshot matrix of the Vomises strain of the first layer composite material and the snapshot matrix of the Vomises strain of the first layer composite material are established, and through the above data processing, the reduced order model of the strain and stress of the flat tail structure can be obtained.

[0065] Therefore, the application introduces health management and digital twin methods into the field of predictive maintenance of aviation equipment by using algorithm design, data structure design, interface design, testing and verification, optimization and improvement, etc., and combines with the predictive maintenance of aviation equipment, so as to realize the breakthrough of the preventive maintenance technology of aviation equipment. BRIEF DESCRIPTION OF DRAWINGS

[0066] Figure 1 The figure is a business architecture diagram of the digital twin calculation model of the body structure predictive maintenance;

[0067] Figure 2 The figure is a workflow diagram of the digital twin calculation model of the body structure predictive maintenance;

[0068] Figure 3 The figure is a business logic flowchart of the structure load monitoring;

[0069] Figure 4 The figure is a rainflow counting load cycle diagram;

[0070] Figure 5 The figure is a rainflow counting load-time history example diagram;

[0071] Figure 6 The figure is a modal truncation number and energy proportion relationship diagram. DETAILED DESCRIPTION

[0072] The technical solutions in the embodiments of the application will be described clearly and completely below with reference to the drawings in the embodiments of the application. Obviously, the described embodiments are only some of the embodiments of the application, but not all the embodiments of the application.

[0073] Meanwhile, in the description of the application, unless otherwise explicitly specified and limited, the terms "connected", "connected" should be understood in a broad sense, for example, can be fixedly connected, or can be detachably connected, or integrally connected; can be mechanically connected, or can be electrically connected; can be directly connected, or can be indirectly connected through an intermediate medium. For those skilled in the art, the specific meaning of the above terms in the application can be understood according to the specific circumstances.

[0074] In order to make the purpose, technical scheme and advantages of the application more clear and apparent, the application will be further described in detail below with reference to the drawings. It should be understood that the specific embodiments described herein are only used to explain the application, and are not used to limit the application.

[0075] The body structure predictive maintenance digital twin computing model (hereinafter referred to as the body structure algorithm model) is to introduce health management and digital twin means into the field of predictive maintenance of aviation equipment, and combine it with the predictive maintenance of aviation equipment, so as to realize the key breakthrough of the preventive maintenance technology of aviation equipment.

[0076] Digital twin is to make full use of physical model, sensor update, operation history and other data, integrate multi-disciplinary, multi-physical quantity, multi-scale and multi-probability simulation process, and complete mapping in virtual space, so as to reflect the whole life cycle process of the corresponding entity equipment. It can collect data through sensing devices, test and dynamically describe each element of the physical entity, and then realize the interactive mapping of physical space and virtual space.

[0077] Predictive maintenance is a maintenance based on state, which analyzes periodic or continuous state monitoring data through a prediction algorithm, determines the running condition of the equipment and estimates the remaining useful life (RUL), predicts the future development trend of the equipment state, and makes maintenance plan in advance according to the prediction.

[0078] Predictive maintenance digital twin is to integrate the concept of predictive maintenance into digital twin on the basis of digital twin. Specifically, the virtual model constructed by digital twin technology can monitor the state of the entity equipment in real time, and predict its future performance change by using historical data and simulation analysis. Once potential failure or performance decline is predicted, maintenance can be carried out in advance to avoid failure and improve the reliability and service life of the equipment. Predictive maintenance digital twin can realize real-time monitoring and fault prediction of the aircraft, ensure flight safety, replace or repair in advance, and improve the reliability and service life of the aircraft.

[0079] The helicopter body structure predictive maintenance digital twin computing model drives the digital twin model of the monitored structure to generate a prediction algorithm by collecting the flight state and structure health feature data of the helicopter, and provides the perfect state information or fault warning of the related structure according to the acquired signals. The helicopter body structure predictive maintenance digital twin computing model can realize real-time and online global monitoring of the helicopter structure, reduce the downtime caused by periodic detection and maintenance, reduce the maintenance and support cost, prevent catastrophic structural failure, and improve the reliability and safety of the helicopter.

[0080] In this process, the data acquisition, processing and transmission of the physical entity are responsible by the on-board health management system, and the collected data is different for different types. It can be divided into three categories: aircraft usage data, structure load data and structure health data.

[0081] The virtual entity can respectively perform structural use monitoring, structural load monitoring and life calculation under structural health monitoring model according to data collected by the entity and in combination with different calculation models.

[0082] The application provides a real-time high-precision structural finite element analysis method based on a reduced-order model, which is used for processing life calculation under structural load monitoring.

[0083] The airborne health management system comprises a body structure predictive maintenance digital twin computing module, sensors for monitoring flight parameters and a maintenance twin data management system for storing data, and the body structure predictive maintenance digital twin computing module has functions of structural safety life monitoring, structural use monitoring (SUMS), structural load monitoring (LMIE) and structural health monitoring (SHEM). Specifically,

[0084] The body structure predictive maintenance digital twin computing is used for performing structural safety life monitoring, structural use monitoring (SUMS), structural load monitoring (LME) and structural health monitoring (SHM) calculation by using collected helicopter flight state and structural health characteristic data.

[0085] For structural load monitoring (LME), based on the above systematic principle, the structural algorithm model comprises a structural load monitoring (LME) module and an operation model, and the operation model comprises a maintenance twin modeling based on a reduced-order model and a maintenance decision calculation software platform. The collected helicopter flight state and structural health characteristic data are used as input to drive the operation model, so as to realize structural load monitoring; the operation model provides two fatigue life calculation methods, can output structural stress strain cloud map results, generates a single-machine load history of a key component according to the structural stress strain cloud map results output by load monitoring, and is used for calculating structural fatigue life, thereby solving the problems of structural finite element fatigue component load monitoring and life calculation based on helicopter flight state and structural health characteristic data.

[0086] The structural algorithm model takes airborne real-time data flow, post-flight flight parameter data files and test bench data as input, is imported into a "decoder" through the maintenance twin data management system, is preprocessed and then connected to the operation model, the operation model is driven by the maintenance twin modeling based on a reduced-order model and the maintenance decision calculation software platform to perform structural load monitoring calculation, and the operation model outputs aviation equipment single-machine dynamic and static analysis results, so that load monitoring data and residual life of structural fatigue components are obtained and stored in the maintenance twin data management system.

[0087] The business architecture of the body structure predictive maintenance digital twin computing model is shown in the accompanyingFigure 1 as shown.

[0088] The structural algorithm model is mainly programmed in Python, uses graphic elements and visualization technology, realizes the free combination and arrangement of the predictive maintenance digital twin calculation model of the aircraft structure, and maps the aircraft type and parameters in the model to the algorithm model to complete initialization by reading the data in the maintenance twin data management system.

[0089] Since the predictive maintenance digital twin calculation model of the aircraft structure (structural algorithm model) is mainly based on multi-node calculation of a large data cluster, it has high requirements for the number and performance of servers. Therefore, the structural algorithm model is deployed on the maintenance twin modeling and maintenance decision calculation software platform based on the reduced-order model, which can share the server computing power on the cloud.

[0090] The structural algorithm model maps the aircraft type and parameters in the model to the operation model to complete initialization by reading the data in the maintenance twin data management system. That is, the structural algorithm model is driven by the maintenance twin modeling and maintenance decision calculation software platform based on the reduced-order model to perform structural load monitoring calculation, and the results are returned to the maintenance twin data management system, while the data are provided for the maintenance twin system micro-service application state monitoring, fault (life) prediction, and maintenance decision calling.

[0091] In summary, the workflow of the system from the start to the end of the service mainly includes the following steps: listening to the request, analyzing the request, creating the task, calculating the inference, returning the result, etc. The specific workflow is shown in the accompanying Figure 2 as shown.

[0092] The software environment equipped by the present application is as follows:

[0093] Operating system: openEuler 22.03LTS SP3;

[0094] Database: Dream 8, MongoDB;

[0095] File library: Hadoop (HDFS);

[0096] Running environment: Python 3.9.12 and necessary third-party algorithm packages (such as Pandas), Jdk1.8;

[0097] Micro-service: SpringCloud and nacos (service configuration and control components, or replaced by corresponding components of aviation cloud);

[0098] Trigger: xxl-job-admin (or replaced by corresponding components of aviation cloud).

[0099] The structural load monitoring module algorithm generally includes five processes of data processing, stress and strain calculation, rainflow counting, material attribute and stress fatigue (strain fatigue) calculation.

[0100] The data processing step is to process the flight parameter data and sensor data input to the algorithm, and to solve the problems of input data jump points and missing data.

[0101] The stress and strain calculation step is to calculate the stress and strain history of each component during flight using the established flight parameter-stress and strain digital twin model.

[0102] The rainflow counting step is to convert the stress and strain history into strain cycle pairs (first point value, second point value) and stress cycles (cycle number, average value and amplitude of each cycle).

[0103] The material attribute step is to configure the corresponding material data according to the selected calculation method.

[0104] The stress and strain fatigue calculation method is to calculate the fatigue life of key components using stress cycles and material attributes, and to provide the damage value and remaining fatigue life of each component.

[0105] The structural load monitoring module provides two fatigue life calculation methods, and the operator needs to determine which damage calculation method to select, and then reads the remaining safe life of each component according to the model, serial number and parts.

[0106] The specific structural load monitoring business logic is shown in the accompanying Figure 3 .

[0107] As shown in the structural load monitoring logic flowchart Figure 3 , it specifically includes the following steps:

[0108] S1, the database reads the current remaining damage life of the life-limited part. According to the model, serial number and key structural component number, the current remaining damage life of the life-limited part (the analysis object includes each life-limited part of each model of the helicopter, and part of the comparison list of each life-limited part is shown in Table 1 below) is queried.

[0109] Table 1 is a comparison table of part of the life-limited parts of the helicopter;

[0110]

[0111]

[0112] S2, read the recorded parameters from the flight parameter data, input the stress and strain digital twin model of the key part, and calculate the surface node stress and strain history. The stress and strain signals are processed, the null values are interpolated, and the stress / strain filtering processing work is performed.

[0113] Specifically includes the following steps S2.1~S2.4:

[0114] S2.1, use flight parameters and load data of each key component to construct a surrogate model that can approximately express the structural load response, i.e. stress-strain digital twin surrogate model;

[0115] Then store the stress-strain of the key components calculated by the stress-strain digital twin surrogate model in the array until all flight parameters are read, and finally obtain the stress-strain history of the key components.

[0116] For key components, according to the measured load (including corresponding flight parameters) data provided by the load department, perform intrinsic orthogonal analysis on all data to obtain eigenvectors with contribution greater than 99%, i.e. measured load of these key components; and for unknown flight state parameters, real-time load of these components can be obtained through the load identification twin model (flight parameter-load surrogate model).

[0117] The system can quickly calculate the load of key components without sacrificing accuracy through the load identification twin model (flight parameter-load surrogate model) described above; the process of real-time load calculation through the load identification twin model (flight parameter-load surrogate model) specifically includes the following steps:

[0118] S2.1.0.1, read a flight parameter data;

[0119] S2.1.0.2, load identification twin model calculates the load of each key component;

[0120] S2.1.0.3, the stress-strain digital twin surrogate model of the key components calculates the stress-strain of the surface nodes of the key components and stores them in the array;

[0121] S2.1.0.4, repeat steps S2.1.0.1~S2.1.0.3 until all flight parameters are read, and finally obtain the stress-strain history of the key components.

[0122] How to construct the load identification twin model (flight parameter-load surrogate model)

[0123] The construction process of the load identification twin model (flight parameter-load surrogate model) specifically includes the following steps:

[0124] S2.1.1, obtain data such as flight altitude, speed, attitude angle, etc. from the flight parameter system as input variables, and corresponding load data (measured or aerodynamic analysis results) as output variables, and write the input and output variables into a CSV file;

[0125] S2.1.2, divide the data into training set and test set, perform scaling features, standardization or normalization preprocessing steps.

[0126] S2.1.3, use machine learning algorithm to train the surrogate model (according to the complexity of the problem, polynomial response surface (RSM), Kriging model, neural network (NN), radial basis function, etc. can be selected as the surrogate model). The model can approximate the true input and output relationship through learning on sample data, ensuring that the model has good generalization ability, that is, the model can make accurate predictions on unseen data;

[0127] S2.1.4, use test set to evaluate the performance of the model, if the performance does not meet the requirements, repeat step (3) until the requirements are met.

[0128] Another part of the key components according to the calculation load data provided by the aerodynamic department, each discrete point flight parameter value as a feature, using BP neural network learning model for training, finally get load twin identification model.

[0129] How to build load identification twin model (flight parameter-load surrogate model)

[0130] Specifically, the specific construction method of the load identification twin model is as follows:

[0131] Since fatigue cracks are often generated and expanded from the surface of the structure, we are more concerned about the stress state of the surface. For a free surface (the surface has no surface load), it can be considered that the stress components perpendicular to the surface are all zero. Select a suitable coordinate system, let the z axis be perpendicular to the surface, then the independent unknown stress components are reduced to 3, see the following formula.

[0132]

[0133] First, establish the finite element model of each key component. If it is a solid structure, cover a layer of shell element with a thickness of 0.00001 on the surface of the body mesh to extract the principal stress strain and shear stress strain of the structure surface;

[0134] The proper orthogonal decomposition method is to use the full-order model calculation results to establish a set of orthogonal bases, and intercept the optimal base combination to make it sufficient to describe the physical characteristics of the original system. First, construct the snapshot matrix A i It can come from the numerical simulation results of theoretical models, or from experimental results of actual physical processes. Where a value represents the numerical result of node b in snapshot a, so the snapshot matrix can reflect the spatial distribution value of the physical field in different samples.

[0135]

[0136] Where L is the number of snapshots selected, and N is the number of grid nodes.

[0137] The essence of the proper orthogonal decomposition method is to find a set of standard orthogonal basis vectors:

[0138] Φ = {Φ1 Φ2 … Φ L}

[0139] The best condition is to require the sampling data T i The error norm of its projection is equivalent to a constrained maximum problem:

[0140]

[0141] Φ T Φ = I

[0142] The vector Φ is the POD canonical orthogonal basis, and the matrix Φ is the POD canonical orthogonal basis set, which is an N by L (N » L) matrix, satisfying the orthogonality as shown in equation (4).

[0143] The constrained maximum problem can be solved by the Lagrange multiplier method, and the original optimization problem is replaced by calculating the eigenvalue problem as follows:

[0144] SS T x = λx

[0145] To calculate the above eigenvalue problem, the singular value decomposition method can be used to obtain the optimal orthogonal basis set for the snapshot matrix S. Therefore, the matrix S ∈ R N×L is singular value decomposed:

[0146] A = BSV T

[0147] where A ∈ R N×N and V ∈ R N×L is an orthogonal matrix, S = R L×L = diag{σ1,σ2,…,σ L} is a diagonal matrix, and each element in S represents a singular value arranged in descending order.

[0148] Since the number of numerical calculation model grid nodes N is much larger than the number of selected snapshot samples L, i.e. N » L, the dimension of the matrix AA T is much larger than that of the matrix A T A, but their non-zero eigenvalues are exactly the same, so the eigenvalue λ T and the eigenvector i of the smaller dimension matrix A

[0149]

[0150] and the dimension of the larger matrix SST φ i

[0151]

[0152] Finally, a set of orthonormal basis vectors Φ that best describe the physical characteristics of the full-order system are solved, which can express the numerical response of any node on the model through linear combination, i.e.

[0153]

[0154] α i are the modal coefficients of each basis φ i

[0155] The order reduction of the full-order model is achieved by truncating a certain number of eigenvectors, because the contribution of each POD canonical vector in expressing the full-order model matrix can be measured by the corresponding eigenvalue, so the truncation operation is based on the principle of energy proportion, that is, select enough first r bases for model order reduction, then the total energy of the reduced subspace should be close to 1 compared with the energy of the full-order model, expressed as

[0156]

[0157] Take γ = 99.99, then the truncated optimal orthonormal basis can be expressed as

[0158]

[0159] Then the unknown quantity to be solved is

[0160]

[0161] First, according to the range of load, DOE test can be carried out, which can use full factorial design, Latin hypercube sampling (LHS) or other DOE methods to generate design points. Assuming that 25 groups of tests are planned, the snapshot matrix is constructed according to the row using the full-order results u i is the numerical simulation result of 25 groups of theoretical models. Among them, the first row is the displacement of working condition 1 The second row is the displacement of working condition 2 and so on.

[0162]

[0163] Where 25 is the number of snapshots selected, m is the number of grid nodes (assuming the model has 58877 nodes), and u is the comprehensive displacement.

[0164] The matrix U is decomposed into three matrices by SVD:

[0165] [B, S, V​T ] = SVD(A)

[0166] Solve the modal coefficients corresponding to 25 samples

[0167]

[0168] Calculate the eigenvalue λ i , that is, the modal energy S of each order i 2 , calculate the relationship I(r) between the number of selected orthonormal basis vectors and modal energy, the results can be seen in the following figure Figure 6 (modal truncation number and energy ratio relationship) When 4 orthonormal basis vectors are selected, I(r) is greater than 99.99%.

[0169] Take the first 4 order eigenvalue vectors with modal energy greater than 99% and the corresponding modal coefficients α r for fitting the displacement field of the structure, and the displacement field calculation needs 0.00245 seconds, that is, 400Hz, after completing the first order reduction.

[0170]

[0171] At the same time, through the Kriging algorithm, 20 groups of load inputs and corresponding modal 20 groups of modal coefficients are learned (5 groups as verification groups), and the proxy model of any load input F and modal coefficient α r is obtained.

[0172] Then the corresponding structure displacement under certain input is

[0173]

[0174] Similarly, the snapshot matrix of the first layer composite Vomises strain and the snapshot matrix of the first layer composite Vomises strain are established, and through the above data processing, the strain and stress reduction model of the flat tail structure can be obtained.

[0175] Load identification twin model (flight parameter-load proxy model) solution to multi-axial fatigue problem

[0176] When dealing with multi-axial fatigue problems, the load identification twin model (flight parameter-load proxy model) often converts the stress tensor into a scalar according to the strength theory, and the von Mises equivalent stress is a formula for calculating the stress state of materials in a complex stress field, as follows.

[0177]

[0178] According to the yield criterion (von Mises criterion), the stress state of the material in the complex stress field can be represented by an equivalent stress, which is called the von Mises equivalent stress, and the node stress-strain tensor field is equivalent to the signed finite element analysis (Signed von Mises) stress (the sign is consistent with the maximum principal stress / strain).

[0179] S2.2, ① Stress and strain signals are processed: when multiple direction strain / stress occurs, there should be a clear principal stress / strain direction (the scheme does not consider the biaxial strain / stress state when the principal stress direction is unknown), and it is assumed that the volume of the component is conserved.

[0180] The following is the conversion equation between true stress and strain:

[0181] σ = ε * E

[0182] In the formula: σ is the stress, unit Mpa; ε is the strain, unit mm; E is the elastic modulus of the component.

[0183] S2.3, ② Interpolation processing

[0184] The load data is completed by interpolation method, and one interpolation is performed on the adjacent two data points at time:

[0185]

[0186] In the formula: x i is the time when the data point is empty; y i is the linear interpolation data point at this time; x i-1 is the previous time; y i-1 is the data point value at the previous time; x i+1 is the next time; y i+1 is the data point value at the next time.

[0187] S2.4, ③ After obtaining the measured stress data, it should be limited and filtered.

[0188] The amplitude limiting is filtered according to the amplitude limiting requirements defined by the user, and this method is filtered by directly setting the threshold value.

[0189] For the burr signal in the strain, the Butterworth high-pass filtering method is adopted, the preset frequency threshold value is set, and the drift problem of the strain load signal is eliminated. The order and cutoff frequency required by Butterworth high-pass filtering should be reasonably set.

[0190]

[0191] In the formula: n is the filter order; ωc This is the cutoff frequency of the filter;

[0192] Specifically, the measured stress data obtained in step S2 is subjected to outlier detection, outlier completion, and smoothing filtering steps.

[0193] For outlier identification, a criterion is established that the point is identified as an outlier when the difference between the measured value and the expected value is greater than the Globus threshold, using a seven-point second-order polynomial forward difference formula.

[0194]

[0195] In the formula: i = 7, 8, 9, ..., N; y and Let be the measured value and the expected value at point i, respectively. Calculate the expected value point by point in data order. and new information With a confidence level of α = 99%, we know from the Grubbs critical value table that K(7,0.99) = 2.097. Points that satisfy the following formula can be regarded as outliers.

[0196]

[0197] The data may contain closely spaced consecutive jump points. When point k is an outlier, the following equation is satisfied: The points are also outliers.

[0198] To avoid misidentifying and eliminating step signals as outliers, we take m = 3. When k + i > k + 3, we can consider y... k y k+i , ..., y k+m All are within the normal range.

[0199] For outlier imputation: After identifying and removing points that have been determined to be outliers, the method is to use the first three points (y, y) of each outlier point as the basis for imputation. k-3 y k-2 y k-1 and the last three points y k+m+1 y k+m+2 y k+m+3 Correct outlier points for the six-point Lagrange interpolation of known values.

[0200]

[0201] In the formula: i, j ≠ k, k+1, ..., k+m; i ≠ j; l = k, k+1, ..., k+m. Correction values ​​are provided for each outlier.

[0202] For data smoothing filtering, follow these steps:

[0203] A2.1, initialize the window size to 7 (pending), prepare an empty list or array to store the filtered data;

[0204] A2.2, starting from the first element of the original data, traverse each element to the last element.

[0205] A2.3, the first 3 and last 3 data points of the data point (if present) will form a window with the current point, add all the data in the window, divide the sum by the number of valid data in the window to get the mean value;

[0206] A2.4, store the calculated mean value in the list prepared earlier as the filtered value of the current data point;

[0207] A2.5, move to the next data point and repeat steps C and D until all points of the original data are traversed.

[0208] A2.6, after traversal, return the list or array storing the filtered values.

[0209] S3, rainflow counting; use rainflow counting method to decompose the stress history into several stress cycles, and then use linear cumulative damage method to calculate the contribution of each stress cycle to fatigue damage.

[0210] Rainflow Counting Method is an algorithm for fatigue life analysis, mainly used for processing non-periodic stress-time history data, which can identify complete stress cycles in stress history, record each cycle to provide stress range and average stress. The input items of rainflow counting are: stress (or strain) history, material fatigue limit, hysteresis filtering threshold (default is 70% of fatigue strength limit, set to 0 to skip);

[0211] The default value of hysteresis filtering threshold refers to "Statistical loads data for MD-82 / 83 aircraft in commercial opertaion", which is set to 70% of fatigue strength limit, and if the hysteresis filtering threshold is 0, it is skipped.

[0212] Hysteresis filtering is used to remove small amplitude loads that do not cause damage to simplify the load spectrum. According to the input (or default) hysteresis filtering threshold, traverse the signal data, remove the small amplitude signal part less than the threshold, and save the remaining history after filtering.

[0213] The specific steps are as follows:

[0214] S3.1, peak-valley value extraction

[0215] Traverse the stress history, find all the local maximum and minimum points, these points will be used as turning points (i.e. points where the stress direction changes), the steps are as follows:

[0216] S3.1.1, from the second element of the sequence, compare the current point with its previous and next points in size.

[0217] S3.1.2, if the value of the current point is greater than the values of its previous and next points, the current point is considered as a peak. If the value of the current point is less than the values of its previous and next points, the current point is considered as a valley.

[0218] S3.1.3, record the identified peak and valley and their position information.

[0219] S3.2, rainflow counting

[0220] Let X be the current range, Y be the range before X; the counting rule of rainflow counting is as follows:

[0221] S3.2.1, read the next peak or valley, if it exceeds the total number of peak and valley points, go to step S3.2.6;

[0222] S3.2.2, if there are less than 3 points, go to step S3.2.1; use the latest 3 peaks and valleys that have not been discarded to determine the ranges X and Y;

[0223] S3.2.3, compare the absolute values of ranges X and Y, ① if X < Y, go to step S3.2.1, ② if X ≥ Y, go to step S3.2.4;

[0224] S3.2.4, count the range Y as one cycle, discard the peak and valley values of the range Y, go to step S3.2.2;

[0225] S3.2.5, the remaining cycles will count each range that has not been counted before as half a cycle;

[0226] S3.2.6, if there are remaining cycles, count them in reverse from the end of the sequence. If there is only one range left (i.e. one incomplete cycle), it can be counted as half or one complete cycle.

[0227] After identifying all possible complete cycles, there may be remaining data points that do not form complete cycles, these remaining points may form an incomplete half cycle, ignore this half cycle. Finally, get the stress cycle pairs and strain cycle pairs, store them in the form of first point value, second point value.

[0228] An example is used to explain the above counting rules, as shown in the attached Figure 5 The rainflow counting example diagram shows a load-time history.

[0229] The counting rules are as follows:

[0230] V1, Y = |A-B|, X = |B-C|; because X > Y, count |A-B| as 1 cycle, and discard A and B points, see Fig. (b). Note that this 1 cycle is formed by pairing part of the range A-B and the range B-C;

[0231] V2, Y = |C-D|, X = |D-E|; because X < Y, continue;

[0232] V3, Y = |D-E|, X = |E-F|; because X < Y, continue;

[0233] V4, Y = |E-F|, X = |F-G|; because X > Y, count |E-F| as 1 cycle, and discard E and F points, see Fig. (c);

[0234] V5, Y = |C-D|, X = |D-G|; because X > Y, count |C-D| as 1 cycle, and discard C and D points, see Fig. (d);

[0235] V6, Y = |G-H|, X = |H-I|; because X < Y, continue to count down to the end;

[0236] V7, Y = |H-I|, X = |G-H|; because X > Y, count |H-I| as 1 cycle, and discard E and F points, see Fig. (e);

[0237] V8, end.

[0238] S3.3, cycle counting

[0239] The cycle pairs recorded in the last step (i.e. peak-valley points) can be used to calculate the mean and amplitude, both of which are only integers, and the formula is as follows:

[0240]

[0241] Finally, the cycle count, cycle mean, and cycle amplitude arrays are obtained.

[0242] S3.4 Rainflow matrix

[0243] An 8X8 two-dimensional array (or similar data structure) is established, where the rows represent different amplitude ranges, the columns represent different mean value ranges, and the values in the array represent the number of cycles that occur in that amplitude and mean value range.

[0244] S4, stress cycle calculation of rainflow counting, to obtain load cycle data; In this scheme, the two-parameter cycle counting method is used, and the two parameters include σ max , σmin , only peak and valley are converted according to stress time history, and rearranged to start from the highest peak or the lowest valley (the larger one in absolute value). The peak, valley recording the specific path of "raindrop" flow in the rainflow counting process is converted into a complete load cycle. See the rainflow counting load cycle diagram in the attached Figure 4 .

[0245] The corresponding stress range from the load cycle: Δσ = σ max -σ min

[0246] Average stress amplitude:

[0247] Finally, the cycle number n of each stress level (stress range, average stress amplitude) is counted i , and the average strain is calculated by goodman, Soderberg or Gerbe transformation.

[0248] S5, linear cumulative damage method to calculate the fatigue damage of flight mission;

[0249] According to the linear cumulative damage method, the results of rainflow counting are accumulated to obtain:

[0250]

[0251] In the formula: N if is the fatigue life (failure cycle number) of the S-N curve of the component under the combination of stress amplitude and average stress level i; n i is the stress cycle number of the load level iD i is the cumulative damage under a flight load cycle n i .

[0252] S6, the remaining damage life of the life component is returned to obtain D i According to the corresponding key component load spectrum for cumulative calculation.

[0253] D total =D ht -D i

[0254] In the formula: D total The remaining damage life of the key component, cycle number; D ht The remaining damage life of the component history, cycle number; D i The cumulative damage under a flight load cycle n i .

[0255] S7, for stress fatigue life calculation;

[0256] The system first reads the material properties of the corresponding components and algorithms in the material database, and the SN curve fitting database is stored in the material database to provide efficient data retrieval, update and management functions.

[0257] For data acquisition of SN curve, the system input items include component number and damage calculation method (stress method, stress-strain method), and the output is SN curve, fatigue strength coefficient σ` f , fatigue plasticity coefficient ε` f , fatigue strength index b, and fatigue plasticity index c.

[0258] The system input items include historical cumulative damage, cumulative flight hours, stress cycles (cycle number, average value and amplitude of each cycle), average stress calculation method (Goodman, Gerber or Soderberg), SN curve, and dispersion coefficient. Stress fatigue life calculation is performed according to the above data; the calculation includes the following steps:

[0259] S7.1, stress correction

[0260] If the average stress effect has been considered in the S / N curve (i.e. the fatigue design uses a series of S / N curves for different stress ratios, and the unassigned stress ratio is obtained by interpolation according to the adjacent stress ratio), this step can be skipped. Otherwise, the average stress correction must be considered.

[0261] The Goodman method is used to connect the two points (σ m , σ a ) and (σ b , 0) in the coordinate system with a straight line, and the intersection point with the y-axis is the σ -1 . It can also be calculated by formula:

[0262]

[0263] Where σ b is the tensile strength of the material.

[0264] Replace σ b with σ s (yield strength of the material), which is the Soderberg method, which is more conservative.

[0265] And the Gerber calculation is as follows:

[0266]

[0267] Select one average stress correction method, traverse the entire stress cycle (average and amplitude), the average value and amplitude of each loop are corrected. That is, the stress cycle (cycle, the average value and amplitude of each cycle) is converted into equivalent stress cycle (cycle, equivalent stress).

[0268] S7.2, damage and life calculation

[0269] The fatigue cumulative damage criterion which is practically useful and widely used in engineering is linear cumulative damage theory (often referred to as Palmgram-Miner theory, or simply referred to as Miner theory). Its formula is

[0270]

[0271] Where n i Stress cycle, N i The cycle number corresponding to the damage time of equivalent stress.

[0272] The total damage of a component is equal to the historical cumulative damage of the component plus the cumulative damage of the current flight.

[0273] If the equation of the SN curve is Nσ m =C, then the damage value d is calculated as:

[0274]

[0275] If the SN curve is given in the form of multiple coordinate points (stress, cycle), the cycle number corresponding to a certain equivalent stress can be calculated by interpolation.

[0276] The current total cumulative damage Where The historical cumulative damage.

[0277]

[0278] S7.3, strain fatigue life calculation

[0279] The strain fatigue life calculation method refers to the local stress-strain method. Practice shows that fatigue failure mainly depends on the plastic deformation of the stress concentration part of the component. The strain fatigue analysis method focuses on the local stress-strain behavior of the stress concentration part. Therefore, it is theoretically more in line with the fatigue failure mechanism. At the same time, this method can conveniently consider the influence of loading sequence, and the fatigue properties data of the material are less, so it is valued by engineering, and has become one of the important methods for estimating the fatigue crack formation life of structure.

[0280] The local stress-strain method is based on the low cycle fatigue test data to estimate the life, and can get satisfactory results for the structure details under high stress, whose local stress exceeds the yield limit of the material.

[0281] The local stress-strain method is based on the low cycle fatigue test data to estimate the life, and can get satisfactory results for the structure details under high stress, whose local stress exceeds the yield limit of the material.

[0282] The strain life curve is expressed by the following formula:

[0283]

[0284] In the formula, σ` f is the fatigue strength coefficient, whose value is approximately equal to the static tensile fracture strength σ f , b is the fatigue strength index, ε` f is the fatigue plasticity coefficient, and c is the fatigue plasticity index.

[0285] The cyclic stress-strain curve can be expressed by the following formula:

[0286]

[0287] K′ is the cyclic strength coefficient, equal to n′ is the cyclic strain hardening index, equal to

[0288] The hysteresis loop curve is expressed by the following formula:

[0289]

[0290] The specific calculation steps are as follows:

[0291] S7.3.1, Calculate the steady-state loop of stress-strain response

[0292] Read the strain rainflow counting cycle pair, from 0 to the first point, consider that the strain is loaded from 0 to ε1, which can be expressed by the cyclic stress-strain curve , whose ε a = ε1-0, then σ1 can be calculated.

[0293]

[0294] The stress-strain path from the first to the second point is expressed by the hysteresis loop curve If the first point value ε1 of the cycle pair is greater than the second point value ε2, it is unloading, and the local stress σ2 is calculated as follows:

[0295]

[0296] σ2 = σ1 - Δσ

[0297] If the cycle is less than the second point value ε2, then the loading, the value of the local stress is calculated as follows

[0298] σ2 = σ1 - Δσ

[0299] According to the above results, σ m and ε a

[0300]

[0301] S7.3.2, damage and life calculation

[0302] According to the Manson-Coffin formula, the fatigue life of each load level and the corresponding fatigue damage are calculated, and the fatigue life of the structure detail is estimated by using the Miner fatigue cumulative damage theory. The Manson-Coffin formula using the Morrow elastic stress linear correction is as follows:

[0303]

[0304] The Newton-Raphson method iteration method can be used to calculate the cycle number under ε a and σ m

[0305]

[0306] The total damage is calculated as follows

[0307]

[0308] The current total cumulative damage Wherein is the historical cumulative damage.

[0309]

[0310] Of course, the present application can have other various embodiments, without departing from the spirit and essence of the present application, those skilled in the art can make various corresponding changes and modifications according to the present application, but these corresponding changes and modifications should all belong to the protection scope of the claims attached to the present application.​

Claims

1. A real-time high-precision structural finite element analysis method based on a reduced-order model, characterized in that ； Based on the body has life piece, under the structure load monitoring, carry out life calculation; wherein, the structure load monitoring takes the airborne real-time data stream collected by the sensor, the flight parameter data file, the test bench data as the input, imports, "decodes", by the maintenance twin modeling based on the reduced order model and the maintenance decision calculation software platform drives the operation model to carry out the calculation of the structure load monitoring, thereby to the body on the life piece carries out the finite element load monitoring data and the residual life calculation; The method sequentially comprises data processing, stress and strain calculation, rainflow counting, material attribute calculation, and stress / strain fatigue calculation processes.

2. The reduced order model based real-time high-accuracy structural finite element analysis method of claim 1, wherein, The method comprises the following specific steps: S1, the database reads the current residual damage life of the life piece; S2, read the recorded parameters from the flight parameter data, input the stress and strain digital twin model of the key component, calculate the surface node stress and strain history; process the stress and strain signals, interpolate the null values, and filter the stress / strain; S3, rainflow counting; use rainflow counting method to decompose the stress history into several stress cycles, and then use linear cumulative damage method to calculate the contribution of each stress cycle to fatigue damage; S4, stress cycle calculation of rainflow counting, get load cycle data; S5, calculate the fatigue damage of flight task by using linear cumulative damage method; S6, return the residual damage life of the life piece, and then calculate the cumulative damage according to the corresponding key component load spectrum; S7, stress fatigue life calculation.

3. The reduced order model based real-time high-accuracy structural finite element analysis method of claim 2, wherein, The step S2 specifically comprises the following steps: S2.1, construct a stress and strain digital twin surrogate model using flight parameters and key component load data; store the stress and strain of the key component surface node calculated by the stress and strain digital twin surrogate model in an array until all flight parameters are read, and finally obtain the stress and strain history of the key component; for unknown flight state parameters, obtain the real-time load of the component through the load identification twin model; The load identification twin model establishes finite element models of each key component, calculates a set of orthogonal bases using the full-order model calculation results, obtains a set of standard orthogonal basis vectors that can best describe the physical characteristics of the full-order system, truncates a certain number of characteristic vectors to realize the reduction of the full-order model, and obtains the reduced order model of the tail structure strain and stress; S2.2, process the stress and strain signals and convert the stress and strain; S2.3, interpolation processing; complete the load data by interpolation method, and interpolate between the adjacent two data points; S2.4, limit and filter the measured stress data.

4. The reduced order model based real-time high-accuracy structural finite element analysis method of claim 3, wherein, The construction method of the load identification twin model in step S2.1 is as follows: For free surface, the stress components perpendicular to the surface are all zero; select the coordinate system, let the z axis be perpendicular to the surface, then the independent unknown stress components are reduced to 3, as follows: First, establish finite element models of each key component, if it is a solid structure, cover a layer of shell element with thickness of 0.00001 on the surface of the body mesh for extracting the principal stress and strain and shear stress and strain of the structure surface; First, construct a snapshot matrix from the full-order results A i numerical simulation results from a theoretical model, or experimental results from an actual physical process; wherein value a numerical result representing a node numbered b in the a-th snapshot; Wherein L is the selected snapshot number, N is the grid node number; find a set of standard orthogonal basis vectors: Φ = {Φ1 Φ2 Φ L} Requesting sampling data T i The error norm of its projection is equivalent to a constrained maximum: Φ T Φ = I The vector Φ is a POD orthonormal basis, and the matrix Φ is a set of POD orthonormal bases, which is an N-row L-column matrix, N > > L, and satisfies the orthogonality; The maximum value problem of the constraint is solved by using the Lagrange multiplier method, and the original optimization problem is replaced by calculating the eigenvalue problem of the following formula: SS T x = λx To compute the above eigenvalue problem, singular value decomposition is applied to the snapshot matrix S to obtain the optimal orthonormal basis set; thus the matrix S ∈ R N×L Singular value decomposition is performed: A = BSV T where A ∈ R N×N and V ∈ R L×L is an orthogonal matrix, S = R L×L = diag{σ1, σ2,..., σ L} is a diagonal matrix, each element of S represents a singular value, arranged in descending order. Since the number of grid nodes N in the numerical computation model is much larger than the number of selected snapshot samples L, i.e., N >> L, therefore, matrix AA T The dimension is much larger than that of matrix A T The dimension of A is known, but the two matrices have the same non-zero eigenvalues, so we first solve for the smaller dimension matrix A. T The eigenvalues ​​of A are λ i With feature vectors Further, the characteristic vector φ T of the large dimension matrix SS i can be calculated. Finally, a set of standard orthonormal basis vectors Φ that can best describe the physical characteristics of the full-order system is solved, which can express the numerical response of any node on the model through linear combination, that is, α i are modal coefficients of each base φ i ​ By truncating a certain number of eigenvectors to realize the order reduction of the full-order model, according to the energy proportion principle, the first r basis is selected to realize the order reduction of the model, then the total energy of the reduced subspace should be close to 1 compared with the energy of the full-order model, expressed as Taking γ = 99.99, the truncated optimal orthonormal basis can be expressed as Then the unknown to be solved is First, according to the range of load, DOE test can be carried out, and full factorial design, Latin hypercube sampling can be used to generate design points; assuming 25 groups of tests are planned, the snapshot matrix is constructed according to the row by using full order results u i The numerical simulation results of 25 groups of theoretical models Wherein, the 1st row is the displacement of working condition 1 The 2nd row is the displacement of working condition 2 And so on; Where 25 is the selected snapshot number, m is the number of grid nodes, assuming that the model has 58877 nodes, and u is the comprehensive displacement; The matrix U is decomposed into three matrices by SVD: [B, S, V T ] = SVD(A) Solve the modal coefficients corresponding to the 25 samples: Computing eigenvalue λ i , i.e. each order modal energy S i 2 , computing the relationship I(r) between the selected standard orthogonal basis vector number and modal energy, obtaining the result that I(r) is greater than 99.99% when 4 groups of orthogonal base vectors are selected; The first four eigenvectors with more than 99% modal energy and the corresponding modal coefficients α r The displacement field for fitting the structure, the displacement field calculation after one order reduction is completed in 0.00245 seconds, that is, 400Hz; Meanwhile, the Kriging algorithm is used to learn the load input of 20 groups of samples and the corresponding modal 20 groups of modal coefficients, 5 groups are used as the verification group, and the proxy model of the arbitrary load input F and the modal coefficient a is obtained r The corresponding structure displacement under i input is Similarly, the snapshot matrix of the Vomises strain of the first layer composite material and the snapshot matrix of the Vomises strain of the first layer composite material are established, and through the above data processing, the strain and stress order reduction model of the flat tail structure can be obtained.

5. The reduced order model based real-time high accuracy structural finite element analysis method of claim 2, wherein, The measured stress data obtained in step S2 are subjected to an outlier judgment, an outlier filling and a smoothing filtering step.

6. The reduced order model based real-time high accuracy structural finite element analysis method of claim 2, wherein, The step S3 specifically comprises the following steps: S3.1, peak and valley value extraction; traverse the stress history to find all local maximum and minimum points, which will be used as turning points; S3.2, rainflow counting; finally obtain stress cycle pairs and strain cycles; S3.3, cycle statistics; obtain the cycle number, cycle mean, and cycle amplitude array; S3.4 rainflow matrix; establish an 8X8 two-dimensional array.

7. The reduced order model based real-time high accuracy structural finite element analysis method of claim 2, wherein, The stress fatigue life calculation specifically comprises the following steps: S7.1, stress correction; use the formula: where σ b is the tensile strength of the material; replace σbwith σs; traverse the entire stress cycle and convert the stress cycle to an equivalent stress cycle; S7.2, use the formula to calculate damage and life; where n i Number of stress cycles, N i Number of cycles at failure for equivalent stress; The historical cumulative damage of the component plus the cumulative damage of the current flight cycle equals the total damage of the component; If the equation of the SN curve is Nσ m = C, then the damage value d is calculated by the formula: If the SN curve is given in the form of multiple coordinate points, the cycle number corresponding to a certain equivalent stress can be calculated by interpolation; the current total cumulative damage wherein is the historical cumulative damage; S7.3, strain fatigue life calculation; Strain life curve calculation: where σ` f is the fatigue strength coefficient, which has a value approximately equal to the static tensile breaking strength σ f , b is the fatigue strength exponent, ε` f is the fatigue plasticity coefficient, and c is the fatigue plasticity exponent; The cyclic stress-strain curve can be expressed by the following formula: K ′ is the cyclic strength coefficient, equal to n ′ is the cyclic strain hardening exponent, equal to The hysteresis loop is represented by the following equation:

Images