Body digital twin fracture risk analysis method based on dynamic bayesian network

Abstract

The application belongs to the technical field of aircraft structure risk analysis, and discloses a body digital twin fracture risk analysis method based on a dynamic Bayesian network, step one: constructing a dynamic Bayesian network; step two: extracting samples from the dynamic Bayesian network at time step t=0, generating n initial particles and respective corresponding weights; step three: increasing the time step t by 1, predicting the particle position at the time step t based on the particle position at the last time step t-1, and calculating the failure probability of each particle; calculating the expected value of the failure probability at the time step t; if the structure does not fail within the time step t, updating the weight of each particle and normalizing; step four: judging whether the time step t reaches the design service life, if yes, ending; otherwise, entering step five; step five: judging whether the structure is inspected at the time step t, if no, entering step three, if yes, performing proportional replacement and weight updating on the current particle, and entering step three after replacement and updating.

Description

Technical Field

[0001] This invention belongs to the field of aircraft structural risk analysis technology, specifically relating to a method for airframe digital twin fracture risk analysis based on dynamic Bayesian networks. Background Technology

[0002] The overall goal of aircraft digital twins is to improve the accuracy of structural health diagnosis and prediction in order to make better maintenance decisions. This can be achieved through structural fatigue fracture risk analysis, which is based on probabilistic damage tolerance analysis (or probabilistic fracture mechanics model). It incorporates various uncertainties / factors into the prediction of structural fatigue crack propagation, calculates the probability of aircraft structure failure due to fatigue fracture throughout its life cycle, and integrates aircraft usage and inspection data to diagnose and update uncertainties in the predictions.

[0003] Fracture failure probability is usually characterized by single flight failure probability (SFPOF), which refers to the probability that a structure will fail in the current flight if it has not failed in the previous flight, or the probability that a structure will fail in the next flight if it has not failed in the current flight.

[0004] The probability of fracture failure needs to be calculated using probabilistic reliability analysis methods. Commonly used methods include the first-order second-moment approximation analytical method (SORM), Monte Carlo numerical simulation, and response surface surrogate model (RSM). Each method has its own advantages and disadvantages. The advantage of the SORM method is its simplicity of calculation, but its accuracy is limited. The Monte Carlo method is the most direct and accurate, but the conventional Monte Carlo method has high computational cost and slow convergence speed, especially when there are many random variables. Summary of the Invention

[0005] The purpose of this invention is to propose a digital twin fracture risk analysis method for mechanical structures based on dynamic Bayesian networks (DBNs). By constructing a dynamic Bayesian network to track fatigue crack propagation and employing a particle filter Bayesian inference algorithm, this method can not only calculate the probability of fatigue fracture failure of the structure throughout its entire life cycle, but also integrate multi-source heterogeneous information from actual observation (or measurement) data to gradually reduce the uncertainty of state variables. This enables crack propagation tracking, failure probability prediction, and Bayesian inference (or diagnosis) to support proactive cluster structure maintenance decisions.

[0006] The technical solution of this invention:

[0007] A method for analyzing fracture risk in an organism based on a dynamic Bayesian network includes the following steps:

[0008] Step 1: Identify the various sources of uncertainty and probability distribution types involved in structural fracture risk analysis, select static and dynamic nodes, and construct a dynamic Bayesian network to track probabilistic crack propagation;

[0009] Step 2: Extract samples from the joint probability distribution of the dynamic Bayesian network at time step t=0 to generate n initial particles and their corresponding weights;

[0010] Step 3: Increment time step t by 1, use the particle filter algorithm for forward propagation, predict the particle position of time step t based on the particle position of the previous time step t-1, and keep the particle weight unchanged; calculate the failure probability of each particle according to the crack size and fracture failure criterion.

[0011] The expected value of the failure probability at time step t is obtained by weighting the failure probabilities based on the particle weights.

[0012] If the structure does not fail within time step t, then the reverse reasoning of the particle filter algorithm is used to update the weight of each particle and normalize it.

[0013] Step 4: Determine if time step t has reached the designed service life. If yes, end the process; otherwise, proceed to step 5.

[0014] Step 5: Determine whether the structure was checked at time step t. If not, proceed to step 3. If yes, determine whether the check result is known.

[0015] If the result is unknown, predict the crack inspection result based on the particle position and weight in step three, and perform proportional replacement and weight update on the current particle. After replacement and update, proceed to step three. If the result is known, update the weight of the current particle based on the crack inspection result, and proceed to step three.

[0016] Furthermore, in step one, the sources of uncertainty include: initial crack size, material fracture toughness, material crack propagation rate parameters, geometric dimensions, flight load spectrum, maximum stress per flight, crack detection probability, and crack size after repair.

[0017] Furthermore, in step one, the static node includes: fracture toughness K. c The crack propagation rate parameter θ;

[0018] Dynamic nodes include: the initial crack size a at each time step t. t 0 Flight payload history F t Stress intensity factor history ΔK t Crack propagation increment Δa t Crack size a t .

[0019] Furthermore, in step one, the dynamic Bayesian network has the following joint probability distribution:

[0020]

[0021] in, p(Δa t |ΔK t ,θ)=1.0,

[0022] Furthermore, in step two, the particle is a sample drawn from the joint probability distribution of the dynamic Bayesian network in one time step, denoted by x. t (i) The weight is represented by ω. t (i) The superscript (i) indicates the i-th particle, where i = 1 to n;

[0023] For the initial time step t=0, each particle contains fracture toughness, crack propagation rate parameters, and initial crack size.

[0024] Furthermore, in step three, the formula for calculating the failure probability of each particle is as follows:

[0025]

[0026] Among them, POF t (i) σ represents the failure probability of the i-th particle at time step t. max This represents the maximum stress encountered during each flight, and typically follows a Gumbel distribution. and It represents the crack size and fracture toughness of the i-th particle at time step t. H(·) is the critical stress that causes failure of the i-th particle under a given crack size and fracture toughness, and σ is the critical stress that causes failure of the i-th particle. max The cumulative distribution function, 1-H(σ) c ) represents the maximum stress σ during each flight. max Exceeding the critical stress σ c The probability of.

[0027] Furthermore, in step three, the reverse reasoning of the particle filter algorithm updates the weight of each particle using the following formula:

[0028]

[0029] in, and Let represent the particle weights at time step t before and after the update based on the assumption that the structure has not failed. Since the particle weights remain constant during forward propagation, ∝ indicates that they are directly proportional.

[0030] Furthermore, in step five, the crack inspection results are predicted, and the current particles are proportionally replaced and their weights updated. The specific process is as follows:

[0031] (1) Calculate the expected value of crack detection probability (PCD) based on the crack size and weight of each particle before inspection, as shown in the following formula:

[0032]

[0033] in and It represents the crack size and weight of the i-th particle at time step t before inspection, and POD(·) is the crack detection probability function.

[0034] (2) If a crack is found, repair it directly, and randomly select n particles from the current particles. repair Each particle is replaced with the initial state after repair, and the weights of the replaced particles are normalized to PCD, where n repair =n*PCD;

[0035] (3) Transfer the remaining n miss Each particle is weighted according to whether a crack is detected, and the weights are normalized to 1-PCD; where n miss =nn repair The weight update formula is as follows:

[0036]

[0037] in, and These represent the particle weights before and after the update based on the crack inspection results, respectively. After normalization,

[0038] Furthermore, in step five, the weight of the current particle is updated based on the crack inspection results, as follows:

[0039]

[0040] in, Indicates the state The following observations were made of y. t The likelihood probability;

[0041] If the crack inspection result is "crack present" or "crack absent," corresponding to "hit" or "miss" respectively, then the likelihood probability is:

[0042]

[0043] If the crack inspection result is the crack size measurement value, then... If we express that the likelihood probability is... This probability needs to be obtained by deriving the statistical model of the POD curve, and it follows a normal distribution, as shown in the following formula:

[0044]

[0045] The specific values ​​are for illustrative purposes only, and the values ​​will differ for different nondestructive testing methods.

[0046] The beneficial effects of this invention are:

[0047] This invention proposes a method for analyzing the fracture risk of a machine body based on a dynamic Bayesian network. By constructing a dynamic Bayesian network that tracks the propagation of fatigue cracks, it performs forward propagation and backward inference of various uncertain variables. This method is more versatile, balances computational efficiency and accuracy, and can integrate actual observation data to perform Bayesian inference and diagnosis of probability distributions, thereby supporting proactive cluster structure maintenance decisions. Attached Figure Description

[0048] Figure 1 A schematic diagram of a dynamic Bayesian network for tracking probabilistic crack propagation;

[0049] Figure 2 The flowchart shows the calculation process for fracture failure probability based on particle filtering.

[0050] Figure 3 A schematic diagram showing the comparison of crack size distribution before and after inspection using Bayesian inference;

[0051] Figure 4 This is a schematic diagram of the predicted single flight failure probability (SFPOF) over the entire lifespan. Detailed Implementation

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

[0053] The solution of the present invention includes the following steps:

[0054] S1: Identify the various sources of uncertainty and probability distribution types involved in structural fracture risk analysis, select static and dynamic nodes, and construct a dynamic Bayesian network (DBN) to track probabilistic crack propagation.

[0055] In one possible embodiment, in step S1, sources of uncertainty include: initial crack size, material parameters (including fracture toughness, crack propagation rate parameters, etc.), geometric dimensional information (thickness, pore size, etc.), flight load spectrum, maximum stress per flight, crack detection probability (POD curve), and crack size after repair. All sources can be discrete variables or continuous random variables.

[0056] In one possible embodiment, in step S1, the static node includes material parameters, denoted by K. c The fracture toughness is represented by θ, which represents the crack propagation rate parameter; the dynamic nodes include the initial crack size a at each time step t. t 0 Flight payload history F t Stress intensity factor history ΔK t Crack propagation increment Δa t Crack size a t ;

[0057] In one possible embodiment, in step S1, the dynamic Bayesian network has the following joint probability distribution:

[0058]

[0059] Where p(x) i |parent(x i )) represents each child node x given that the parent node is known. i The conditional probability distribution.

[0060] Since the stress intensity factor has a deterministic functional relationship with the load and crack size, therefore Since the crack propagation increment is also a deterministic function of the stress intensity factor and crack propagation rate parameter, therefore p(Δa) t |ΔK t ,θ)=1.0; Since the crack size is equal to the sum of the initial size and the crack propagation increment, therefore

[0061] Under normal circumstances, the flight load history for each takeoff and landing can be considered a fixed value (calculated based on flight parameter data collected by airborne equipment and the "flight parameter-load" model). At this time, p(F t ) = 1.0; it can also be assumed that there is an unavoidable certain error between the collected data and the "flight parameter-load" model, then p(F) = 1.0; t The fracture toughness K is a probability distribution, typically a normal distribution. c The crack propagation rate parameter θ follows a normal distribution or a log-normal distribution, and the initial crack size at time zero is given. It follows a log-normal distribution or a Weibull distribution.

[0062] S2: Extract samples from the joint probability distribution of the dynamic Bayesian network at time step t=0 to generate n initial particles and their corresponding weights;

[0063] In one possible embodiment, in step S2, the particle is a sample drawn from the joint probability distribution in one time step, denoted by x. t (i) The weight is represented by ω. t (i) The superscript (i) indicates the i-th particle, where i = 1 to n. For the initial time step t = 0, each particle contains fracture toughness, crack propagation rate parameters, and initial crack size.

[0064] S3: Increment time step t by 1, use the forward propagation of the particle filter algorithm to predict the particle position at time step t based on the particle position at the previous time step t-1, and keep the particle weights unchanged; calculate the failure probability of each particle according to the crack size and fracture failure criteria; perform a weighted average of the failure probabilities according to the particle weights to obtain the expected value of the failure probability at time step t; if the structure does not fail within time step t, use the backward reasoning of the particle filter algorithm to update and normalize the weights of each particle.

[0065] In one possible embodiment, in step S3, the failure probability of each particle is calculated. The most commonly used fracture failure criterion is the fracture toughness criterion, which is that the maximum stress intensity factor during flight exceeds the fracture toughness of the material or the maximum load during flight exceeds the critical strength of the structure. The calculation formula is as follows:

[0066]

[0067] Among them, POF t (i) σ represents the failure probability of the i-th particle at time step t. max This represents the maximum stress encountered during each flight, typically assumed to follow a Gumbel distribution. and It represents the crack size and fracture toughness of the i-th particle at time step t. H(·) is the critical stress that leads to failure given the crack size and fracture toughness. max The cumulative distribution function, 1-H(σ) c ) represents the maximum stress σ during each flight. max Exceeding the critical stress σ c The probability of.

[0068] In one possible embodiment, the formula for calculating the expected value of the failure probability in step S3 is:

[0069]

[0070] Among them, SPOF t This represents the expected failure probability at time step t. This represents the weight of the i-th particle.

[0071] In one possible embodiment, in step S3, the structure not failing at time step t means that the structure has not experienced fatigue fracture at time step t. This assumption is a premise of SFPOF because it is defined as the probability that the structure will fail in the next flight if it has not failed in the current flight.

[0072] In one possible embodiment, in step S3, the backward inference of the particle filtering algorithm updates the weight of each particle using the following formula:

[0073]

[0074] in, and Let represent the particle weights at time step t before and after the update based on the assumption that the structure has not failed. Since the particle weights remain constant during forward propagation, ∝ indicates that they are directly proportional.

[0075] S4: Determine whether time step t has reached the designed service life. If yes, end; otherwise, proceed to step S5.

[0076] S5: Determine whether the structure is checked at time step t. If not, proceed to step S3. If yes, determine whether the check result is known. If the result is unknown, predict the crack check result based on the particle position and weight in step S3 and perform proportional replacement and weight update on the current particle. After replacement and update, proceed to step S3. If the result is known, update the weight of the current particle based on the crack check result. After update, proceed to step S3.

[0077] In one possible embodiment, in step S5, the crack inspection result is predicted and the current particles are proportionally replaced and their weights are updated. The specific steps are as follows:

[0078] 1) Calculate the expected probability of crack detection (PCD) based on the crack size and weight of each particle before inspection, using the following formula:

[0079]

[0080] in and t represents the crack size and weight of the i-th particle before inspection at time step t, and POD(·) is the crack detection probability function, which is a characterization of non-destructive testing capability. It is usually represented by the cumulative distribution function of log-normal distribution or exponential distribution.

[0081] 2) If a crack is detected, repair it directly, and randomly select n particles from the current particles. repair Each particle is replaced with the initial state after repair, and the weights of the replaced particles are normalized to PCD, where n repair =n*PCD;

[0082] 3) Transfer the remaining n miss Each particle is weighted according to whether a crack was detected (formula below), and the weights are normalized to 1-PCD. Where n miss =nn repair .

[0083]

[0084] in, and These represent the particle weights before and after the update based on the crack inspection results, respectively. After normalization,

[0085] In one possible embodiment, in step S5, the weight of the current particle is updated based on the crack inspection results, specifically using the following formula:

[0086]

[0087] in, Indicates the state The following observations were made of y. t The likelihood probability.

[0088] If the crack inspection result is "crack present" or "crack absent," corresponding to "hit" or "miss" respectively, then the likelihood probability is:

[0089]

[0090] If the crack inspection result is the crack size measurement value, then... If we express that the likelihood probability is... This probability needs to be obtained by deriving the statistical model of the POD curve. It generally follows a normal distribution, as shown in the following formula:

[0091]

[0092] The specific values ​​are for illustrative purposes only, and the values ​​will differ for different nondestructive testing methods.

[0093] Example 1

[0094] Taking a critical component of the fuselage frame of an aircraft as an example, with a designed service life of 5000 takeoffs and landings, and structural inspections scheduled as follows: the first inspection at 2600 takeoffs and landings, and subsequent inspections at 1100 takeoffs and landings, the probability of fatigue fracture failure of this component over its entire service life is predicted using the technical method of this invention.

[0095] 1) Following step S1, sort out the various sources of uncertainty and probability distribution types involved in the fracture risk analysis of this part, select static nodes and dynamic nodes, and construct a dynamic Bayesian network to track the probabilistic crack propagation.

[0096] Uncertainty sources include: initial crack size, material fracture toughness, crack propagation rate parameters, maximum stress per flight, crack detection probability (POD curve), and crack size after repair. The probability distribution and parameters are listed in Table 1. Other parameters are taken as fixed values, such as geometric dimensional information (thickness, pore size, etc.).

[0097] Table 1. Probability distributions and parameters of various uncertainty sources

[0098]

[0099] Static nodes are time-independent variables, including material parameters and crack propagation rate parameters. Dynamic nodes are variables whose state changes over time, including the initial crack size, flight load history, stress intensity factor history, crack propagation increment, and failure probability at each time step.

[0100] The constructed dynamic Bayesian network framework, such as Figure 1 K c Fracture toughness is represented by C and m, crack propagation rate parameters are represented by a. t 0 This represents the initial crack size at time step t (= crack size a at time step t-1). t-1 ), use F t The flight payload history at time step t is represented by ΔK. t The stress intensity factor history at time step t is represented by Δa. t The increment of crack propagation at time step t is represented by a. t The crack size at time step t is represented by POF. t This represents the failure probability at time step t.

[0101] Figure 1In this diagram, elliptical nodes represent random nodes, meaning the variable is random given the parent node; an arrow pointing to an elliptical node indicates a conditional probability distribution. Triangular nodes represent function nodes, meaning the variable is deterministic given the parent node; an arrow pointing to a triangular node indicates a deterministic function. Rectangular nodes represent variables with observed values ​​(including load history data and crack detection data). Solid elliptical nodes represent continuous variables, and dashed elliptical nodes represent discrete variables. Solid arrows connect nodes within the same time step, and dashed arrows connect nodes across different time steps.

[0102] The joint probability distribution of the dynamic Bayesian network is as follows:

[0103]

[0104] Since the relationship between nodes is a deterministic function, therefore p(Δa t |ΔK t ,θ)=1.0, The flight load history for each takeoff and landing is assumed to be a fixed value (calculated based on flight parameter data collected by onboard equipment and the "flight parameter-load" model), therefore p(F t = 1.0.

[0105] 2) The calculation process for fracture failure probability based on particle filtering is as follows: Figure 2 First, following step S2, 10,000 initial particles and weights are extracted from the joint probability distribution of the initial crack size, fracture toughness, and crack propagation rate parameters.

[0106] 3) Following steps S3 to S5, assuming the crack inspection results at the predetermined inspection times (2600°C, 3700°C, 4800°C) are unknown, it is necessary to predict the crack inspection results and perform proportional replacement and weight update on the current particles. The specific method is as follows:

[0107] Taking the inspection after 2600 as an example, based on the size and weight of each particle before inspection, the expected value of crack detection probability PCD = 0.15 is calculated. Assuming that the crack is repaired as soon as it is detected, 1500 particles are randomly selected from the current particles to replace the initial state after repair (based on the probability distribution of crack size after repair), and the weight of the replaced particles is normalized to 0.15; the remaining 8500 particles are updated with Bayesian weights as if the crack was not detected, and the weights are normalized to 0.85. Figure 3 This is a comparison of crack size distribution before and after inspection using Bayesian inference. The solid line represents the crack size distribution before inspection, the dashed line represents the crack size distribution after weight update, and the dotted dashed line represents the crack size distribution after weight normalization.

[0108] Repeat steps S3 to S5 until time step t reaches the design service life. The SFPOF prediction result is as follows: Figure 4 As can be seen, the SFPof increases with flight time until the scheduled inspection time. Since a certain percentage of the crack size is reset after inspection, the SFPof decreases after inspection.

[0109] The above description is merely a specific embodiment of the present invention, providing a detailed description of the invention. Parts not covered herein are conventional techniques. However, the scope of protection of the present invention is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in the present invention should be included within the scope of protection of the present invention. The scope of protection of the present invention should be determined by the scope of the claims.

Claims

1. A method for analyzing fracture risk in an organism based on a dynamic Bayesian network, characterized in that: The method includes the following steps: Step 1: Identify the various sources of uncertainty and probability distribution types involved in structural fracture risk analysis, select static and dynamic nodes, and construct a dynamic Bayesian network to track probabilistic crack propagation; Step 2: Extract samples from the dynamic Bayesian network at time step t=0 to generate n initial particles and their corresponding weights; Step 3: Increment time step t by 1, and use the particle filter algorithm for forward propagation. Predict the particle position at time step t based on the particle position at the previous time step t-1, keeping the particle weights unchanged. Calculate the failure probability of each particle according to the crack size and fracture failure criterion. The calculation formula is as follows: in, Let represent the failure probability of the i-th particle at time step t. σ max This represents the maximum stress encountered during each flight, and typically follows a Gumbel distribution. and It represents the crack size and fracture toughness of the i-th particle at time step t. It is the critical stress that causes the i-th particle to fail under a given crack size and fracture toughness. H ( · )yes σ max The cumulative distribution function, 1- H ( σ c () indicates the maximum stress during each flight. σ max Exceeding the critical stress σ c The probability of; The expected value of the failure probability at time step t is obtained by weighting the failure probabilities based on the particle weights. If the structure does not fail within time step t, then the reverse reasoning of the particle filter algorithm is used to update the weight of each particle and normalize it. Step 4: Determine if time step t has reached the designed service life. If yes, end the process; otherwise, proceed to step 5. Step 5: Determine whether the structure was checked at time step t. If not, proceed to step 3. If yes, determine whether the check result is known. If the result is unknown, predict the crack inspection result based on the particle position and weight in step three, and perform proportional replacement and weight update on the current particle. After replacement and update, proceed to step three. If the result is known, update the weight of the current particle based on the crack inspection result, and proceed to step three.

2. The method according to claim 1, characterized in that: In step one, the sources of uncertainty include: initial crack size, material fracture toughness, material crack propagation rate parameters, geometric dimensions, flight load spectrum, maximum stress per flight, crack detection probability, and crack size after repair.

3. The method according to claim 2, characterized in that: In step one, the static node includes: fracture toughness. K c Crack propagation rate parameter θ ; Dynamic nodes include: the initial crack size at each time step t. 0 Flight payload history F t Stress intensity factor history ΔK t Crack propagation increment Δa t Crack size a t .

4. The method according to claim 3, characterized in that: In step one, the dynamic Bayesian network has the following joint probability distribution: in, , , .

5. The method according to claim 4, characterized in that: In step two, the particle is a sample drawn from the joint probability distribution of the dynamic Bayesian network in one time step. x t (i) The weight is represented by ω. t (i) The superscript (i) indicates the i-th particle, where i = 1 to n; For the initial time step t=0, each particle contains fracture toughness, crack propagation rate parameters, and initial crack size.

6. The method according to claim 5, characterized in that: In step three, based on the reverse reasoning of the particle filter algorithm, the weight of each particle is updated using the following formula: in, and Let represent the particle weights at time step t before and after the update based on the assumption that the structure has not failed. Since the particle weights remain constant during forward propagation, then ;∝ indicates a direct proportion.

7. The method according to claim 6, characterized in that: In step five, the crack inspection results are predicted, and the current particles are proportionally replaced and their weights are updated. The specific process is as follows: (1) Calculate the expected value of crack detection probability based on the crack size and weight of each particle before inspection. PCD, The formula is as follows: in and It represents the crack size and weight of the i-th particle at time step t before inspection. POD ( · ) is the crack detection probability function; (2) If a crack is found, repair it directly by randomly selecting from the current particles. n repair Each particle is replaced with the initial state after repair, and the weights of the replaced particles are normalized to... PCD ,in n repair = n*PCD ; (3) Put the remaining n miss Each particle is weighted according to whether a crack was detected, and the weights are normalized to 1- PCD ;in, n miss = n - n repair The weight update formula is as follows: in, and These represent the particle weights before and after the update based on the crack inspection results, respectively. After normalization, .

8. The method according to claim 7, characterized in that: In step five, the process of updating the weight of the current particle based on the crack inspection results is as follows: in, Indicates the state The following observations y t The likelihood probability; If the crack inspection result is "crack present" or "crack absent," it corresponds to hit / m If iss, then the likelihood probability is: If the crack inspection result is the crack size measurement value, then... If we express that the likelihood probability is... This probability needs to be obtained using a statistical model derived from the POD curve. It follows a normal distribution, as shown in the following formula: 。

Images