CN112784495B - Mechanical structure real-time fatigue life prediction method based on data driving - Google Patents

Mechanical structure real-time fatigue life prediction method based on data driving Download PDF

Info

Publication number
CN112784495B
CN112784495B CN202110116397.0A CN202110116397A CN112784495B CN 112784495 B CN112784495 B CN 112784495B CN 202110116397 A CN202110116397 A CN 202110116397A CN 112784495 B CN112784495 B CN 112784495B
Authority
CN
China
Prior art keywords
crack
boundary
representing
displacement
component
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202110116397.0A
Other languages
Chinese (zh)
Other versions
CN112784495A (en
Inventor
谢贵重
钟玉东
李�浩
杜文辽
冯世哲
邬昌军
李客
巩晓赟
王良文
刘林
张世欣
王滔
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Zhengzhou University of Light Industry
Original Assignee
Zhengzhou University of Light Industry
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Zhengzhou University of Light Industry filed Critical Zhengzhou University of Light Industry
Priority to CN202110116397.0A priority Critical patent/CN112784495B/en
Publication of CN112784495A publication Critical patent/CN112784495A/en
Application granted granted Critical
Publication of CN112784495B publication Critical patent/CN112784495B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/27—Design optimisation, verification or simulation using machine learning, e.g. artificial intelligence, neural networks, support vector machines [SVM] or training a model
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/02—Neural networks
    • G06N3/04—Architecture, e.g. interconnection topology
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/02—Neural networks
    • G06N3/08—Learning methods
    • G06N3/084—Backpropagation, e.g. using gradient descent
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/04—Ageing analysis or optimisation against ageing

Landscapes

  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Evolutionary Computation (AREA)
  • General Engineering & Computer Science (AREA)
  • Artificial Intelligence (AREA)
  • Software Systems (AREA)
  • General Physics & Mathematics (AREA)
  • Health & Medical Sciences (AREA)
  • General Health & Medical Sciences (AREA)
  • Mathematical Physics (AREA)
  • Computing Systems (AREA)
  • Molecular Biology (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Biomedical Technology (AREA)
  • Biophysics (AREA)
  • Computational Linguistics (AREA)
  • Data Mining & Analysis (AREA)
  • Geometry (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Medical Informatics (AREA)
  • Computer Hardware Design (AREA)
  • Investigating Strength Of Materials By Application Of Mechanical Stress (AREA)
  • Testing Of Devices, Machine Parts, Or Other Structures Thereof (AREA)

Abstract

本发明提出了一种基于数据驱动的机械结构实时疲劳寿命预测方法，其步骤为：首先，获得机械结构中对应的Paris模型的指数m和系数c的取值范围，并生成一系列指数和系数；其次，在机械结构上随机获取观测点，针对一组（mq，cq）进行循环，使用对偶互易边界元法对观测点进行疲劳裂纹扩展分析，获得观测点的位移和实时疲劳寿命信息，组成数据信息对；直至遍历所有（m，c），得到数据集；再将数据集输入BP神经网络中进行训练，得到BP神经网络模型；最后，采集机械结构中观测点的位移，并将观测点的位移输入BP神经网络模型中，得到观测点的实时疲劳寿命信息。本发明仅通过观测点位移就可预测机械结构的疲劳寿命，节省大量的时间和成本。

Figure 202110116397

The present invention proposes a data-driven real-time fatigue life prediction method for mechanical structures, the steps of which are: first, obtain the value ranges of the exponent m and coefficient c of the Paris model corresponding to the mechanical structure, and generate a series of exponents and coefficients ;Secondly, randomly obtain observation points on the mechanical structure, carry out a cycle for a set of (m q , c q ), use the dual reciprocal boundary element method to analyze the fatigue crack growth of the observation points, and obtain the displacement and real-time fatigue life of the observation points. information to form data information pairs; until all (m, c) are traversed, the data set is obtained; then the data set is input into the BP neural network for training, and the BP neural network model is obtained; finally, the displacement of the observation point in the mechanical structure is collected, and the The displacement of the observation point is input into the BP neural network model, and the real-time fatigue life information of the observation point is obtained. The invention can predict the fatigue life of the mechanical structure only by observing the displacement of the point, and saves a lot of time and cost.

Figure 202110116397

Description

一种基于数据驱动的机械结构实时疲劳寿命预测方法A data-driven real-time fatigue life prediction method for mechanical structures

技术领域technical field

本发明涉及机械产品的疲劳裂纹的检测与预测技术领域，特别是指一种基于数据驱动的机械结构实时疲劳寿命预测方法。The invention relates to the technical field of detection and prediction of fatigue cracks of mechanical products, in particular to a data-driven real-time fatigue life prediction method of mechanical structures.

背景技术Background technique

在机械材料加工过程，不可避免的会出现空洞、夹杂及裂纹等缺陷，在复杂工况及循环载荷作用情况下，产生应力过度集中，引发疲劳断裂，甚至会造成灾难性后果。据统计，疲劳断裂造成的破坏约占到全部力学破坏的50％-90％，这都与工程结构中存在的缺陷相关。为了预防安全事故，减少经济损失，预测机械结构的剩余疲劳寿命成为人们关注的重点。In the process of mechanical material processing, defects such as voids, inclusions and cracks will inevitably appear. Under complex working conditions and cyclic loads, excessive stress concentration will occur, causing fatigue fractures, and even catastrophic consequences. According to statistics, the damage caused by fatigue fracture accounts for about 50%-90% of the total mechanical damage, which is related to the defects existing in the engineering structure. In order to prevent safety accidents and reduce economic losses, predicting the remaining fatigue life of mechanical structures has become the focus of attention.

随着计算力学及计算机技术的发展，越来越多的数值方法被用来模拟机械结构裂纹的扩展，进而预测机械结构的疲劳寿命。使用数值方法可以精确模拟裂纹的扩展过程，掌握裂纹扩展规律，从而进行机械结构的强度分析和寿命预测。在这些数值方法中，最为有效的两种方法分别为有限元法与边界元法。With the development of computational mechanics and computer technology, more and more numerical methods are used to simulate the crack propagation of mechanical structures, and then predict the fatigue life of mechanical structures. Using the numerical method, the crack propagation process can be accurately simulated, and the crack propagation law can be mastered, so as to carry out the strength analysis and life prediction of the mechanical structure. Among these numerical methods, the two most effective methods are the finite element method and the boundary element method.

有限元法进行裂纹的扩展模拟时应力、应变精度比位移低一阶，裂纹尖端高应力区需要布置高密度网格，并且裂纹只能限制在网格线上，将裂纹面设置为单元的边界、裂尖设置为单元的节点，模拟裂纹扩展时，随着裂纹的扩展，网格需要重新划分。在有限元框架中发展起来的扩展有限元，在常规单元形函数中加入能反映不连续性的阶跃函数及裂尖渐进位移场函数。网格与裂纹相互独立，裂纹扩展不需重构网格，能方便地分析含裂纹体的强不连续问题。然而，为了获得裂纹扩展高精度模拟结果，依然需要大量的网格。边界元法作为一种半解析的数值方法，利用的面力和位移具有同阶精度的优势，可提高机械结构应变、应力分析的精度，在裂纹扩展模拟中更具优势。边界元法基本解的奇异性质，更适合用于机械结构实时疲劳寿命预测这类应力奇异性问题。When the finite element method is used for crack propagation simulation, the stress and strain accuracy is one order lower than the displacement. The high-stress area at the crack tip needs to be arranged with high-density mesh, and the crack can only be limited to the mesh line, and the crack surface is set as the boundary of the element. , and the crack tip is set as the node of the element. When simulating crack propagation, the mesh needs to be re-divided as the crack expands. The extended finite element developed in the finite element framework adds a step function that can reflect the discontinuity and the asymptotic displacement field function of the crack tip to the conventional element shape function. The mesh and the crack are independent of each other, and the crack propagation does not need to reconstruct the mesh, which can easily analyze the strong discontinuity of the cracked body. However, in order to obtain high-accuracy simulation results of crack propagation, a large number of meshes are still required. As a semi-analytical numerical method, the boundary element method utilizes the advantage of the same-order accuracy of surface force and displacement, which can improve the accuracy of mechanical structural strain and stress analysis, and is more advantageous in crack propagation simulation. The singularity of the basic solution of the boundary element method is more suitable for stress singularity problems such as real-time fatigue life prediction of mechanical structures.

由于使用Paris模型，机械结构参数的不确定性，Paris模型中的系数需要通过试验测定大量的数据点进行拟合获得，带有随机性。如果使用摄动法来处理的话，往往需要应用到先验知识。由于裂纹扩展过程参数的变换，摄动技术难以与数值方法结合，应用于裂纹扩展分析。这就给对偶互易边界元法预测机械结构实时疲劳寿命是带来了巨大挑战。Due to the uncertainty of mechanical structure parameters using the Paris model, the coefficients in the Paris model need to be obtained by fitting a large number of data points through experiments, with randomness. If the perturbation method is used to deal with it, it often needs to apply prior knowledge. Due to the transformation of the parameters in the crack growth process, the perturbation technique is difficult to combine with the numerical method and applied to the crack growth analysis. This brings a huge challenge to the real-time fatigue life prediction of mechanical structures by dual reciprocal boundary element method.

近年来，机器学习技术引起了人们的关注。反向传播神经网络(BPNN)作为人工神经网络中的一种，目前已经在基于大量物理实验或数值模拟数据集的工程预测问题中获得了成功应用。计算方法与神经网络技术相结合，在处理裂纹分析中显示出巨大的潜力。然而，关于机械结构的实时疲劳寿命预测的文献很少。In recent years, machine learning techniques have attracted attention. Back-propagation neural network (BPNN), as one of artificial neural networks, has been successfully applied in engineering prediction problems based on a large number of physical experiments or numerical simulation data sets. Computational methods combined with neural network techniques show great potential in dealing with crack analysis. However, there is little literature on real-time fatigue life prediction of mechanical structures.

发明内容SUMMARY OF THE INVENTION

针对上述背景技术中存在的不足，本发明提出了一种基于数据驱动的机械结构实时疲劳寿命预测方法，解决了现有技术中无法预测机械结构的疲劳寿命，造成机械结构的检测精度降低的技术问题。In view of the deficiencies in the above-mentioned background technologies, the present invention proposes a data-driven real-time fatigue life prediction method for mechanical structures, which solves the problem that the fatigue life of mechanical structures cannot be predicted in the prior art, which reduces the detection accuracy of mechanical structures. question.

本发明的技术方案是这样实现的：The technical scheme of the present invention is realized as follows:

一种基于数据驱动的机械结构实时疲劳寿命预测方法，其步骤如下：A data-driven real-time fatigue life prediction method for mechanical structures, the steps are as follows:

步骤一：根据试验手段获得机械结构中对应的Paris模型的指数m和系数c的取值范围，并分别在指数m和系数c的取值范围内采用优先数系生成一系列指数和系数，得到Q组Paris模型参数对(mq，cq)，其中，q＝1,2,…,Q；Step 1: Obtain the value range of the exponent m and coefficient c of the Paris model corresponding to the mechanical structure according to the test method, and use the priority number system to generate a series of exponents and coefficients within the value range of the exponent m and coefficient c respectively, and obtain Q groups of Paris model parameter pairs (m q , c q ), where q=1, 2, . . . , Q;

步骤二：在机械结构上设置观测点，针对第q组Paris模型参数对(mq，cq)进行循环，使用对偶互易边界元法对机械结构进行疲劳裂纹扩展分析，获得观测点的位移和机械结构实时疲劳寿命信息，并将观测点的位移和机械结构实时疲劳寿命信息组成数据信息对；Step 2: Set observation points on the mechanical structure, cycle (m q , c q ) for the q-th group of Paris model parameters, use the dual reciprocal boundary element method to analyze the fatigue crack growth of the mechanical structure, and obtain the displacement of the observation point and the real-time fatigue life information of the mechanical structure, and combine the displacement of the observation point and the real-time fatigue life information of the mechanical structure into a data information pair;

步骤三：循环执行步骤二，直至遍历Q组Paris模型参数对，得到数据集，并将数据集分为训练集和测试集，其中，数据集包括一系列数据信息对；Step 3: Execute step 2 in a loop until traversing Q groups of Paris model parameter pairs to obtain a data set, and divide the data set into a training set and a test set, wherein the data set includes a series of data information pairs;

步骤四：将训练集输入BP神经网络中进行训练，得到BP神经网络模型，并利用测试集对BP神经网络模型进行验证；Step 4: Input the training set into the BP neural network for training, obtain the BP neural network model, and use the test set to verify the BP neural network model;

步骤五：采集机械结构中观测点的位移，并将观测点的位移输入BP神经网络模型中，得到机械结构的实时疲劳寿命信息。Step 5: Collect the displacement of the observation point in the mechanical structure, and input the displacement of the observation point into the BP neural network model to obtain the real-time fatigue life information of the mechanical structure.

所述使用对偶互易边界元法对机械结构进行疲劳裂纹扩展分析，获得观测点的位移和实时疲劳寿命信息的方法为：The method of using the dual reciprocal boundary element method to analyze the fatigue crack growth of the mechanical structure, and obtaining the displacement of the observation point and the real-time fatigue life information is:

S2.1、构建对偶互易边界积分方程，包括位移边界积分方程和面力边界积分方程；S2.1. Construct dual reciprocal boundary integral equations, including displacement boundary integral equations and surface force boundary integral equations;

S2.2、分别对位移边界积分方程和面力边界积分方程进行离散处理，得到：S2.2. Discrete the displacement boundary integral equation and the surface force boundary integral equation respectively, and obtain:

Figure BDA0002920814570000021
Figure BDA0002920814570000021

其中，H*表示位移基本解的积分，G*表示面力基本解的积分，S表示非裂纹边界，

Figure BDA0002920814570000022
表示裂纹上边界，
Figure BDA0002920814570000023
表示裂纹下边界，uS表示非裂纹边界上位移，
Figure BDA0002920814570000024
表示裂纹上边界上的位移，
Figure BDA0002920814570000025
表示裂纹下边界上的位移，tS表示非裂纹边界上的面力，
Figure BDA0002920814570000026
表示裂纹上边界上的面力，
Figure BDA0002920814570000027
表示裂纹下边界上的面力；where H * represents the integral of the fundamental solution of displacement, G * represents the integral of the fundamental solution of surface force, S represents the non-crack boundary,
Figure BDA0002920814570000022
represents the upper boundary of the crack,
Figure BDA0002920814570000023
represents the lower boundary of the crack, u S represents the upper displacement of the non-crack boundary,
Figure BDA0002920814570000024
represents the displacement on the upper boundary of the crack,
Figure BDA0002920814570000025
represents the displacement on the lower boundary of the crack, t S represents the surface force on the non-crack boundary,
Figure BDA0002920814570000026
represents the surface force on the upper boundary of the crack,
Figure BDA0002920814570000027
represents the surface force on the lower boundary of the crack;

S2.3、对公式(5)施加边界条件，并整理得到：S2.3. Apply boundary conditions to formula (5), and get:

Figure BDA0002920814570000031
Figure BDA0002920814570000031

其中，A*表示位移或面力已知量，xS表示非裂纹边界上未知的位移或面力，

Figure BDA0002920814570000032
表示裂纹上边界上未知的位移或面力，
Figure BDA0002920814570000033
表示裂纹下边界上未知的位移或面力，bS表示非裂纹边界上已知的位移或面力合并得到的分量，
Figure BDA0002920814570000034
表示裂纹上边界上已知的位移或面力合并得到的分量，
Figure BDA0002920814570000035
表示裂纹下边界上已知的位移或面力合并得到的分量；where A * represents the known displacement or surface force, x S represents the unknown displacement or surface force on the non-crack boundary,
Figure BDA0002920814570000032
represents the unknown displacement or surface force on the upper boundary of the crack,
Figure BDA0002920814570000033
represents the unknown displacement or surface force on the lower boundary of the crack, b S represents the combined component of the known displacement or surface force on the non-crack boundary,
Figure BDA0002920814570000034
represents the component obtained by combining the known displacements or surface forces on the upper boundary of the crack,
Figure BDA0002920814570000035
Represents the component obtained by combining the known displacements or surface forces on the lower boundary of the crack;

利用高斯消去法对公式(6)求解，获得边界的位移和面力，并将边界的位移和面力代入公式(1)中获得观测点的位移；Use the Gaussian elimination method to solve the formula (6), obtain the displacement and surface force of the boundary, and substitute the displacement and surface force of the boundary into the formula (1) to obtain the displacement of the observation point;

S2.4、根据边界的位移利用相互作用积分算法计算混合模态应力强度因子：S2.4. According to the displacement of the boundary, use the interaction integral algorithm to calculate the mixed mode stress intensity factor:

Figure BDA0002920814570000036
Figure BDA0002920814570000036

其中，M(1,2)表示相互作用积分，

Figure BDA0002920814570000037
为互作用应变能密度，
Figure BDA0002920814570000038
Figure BDA0002920814570000039
表示状态二下的i方向位移分量，
Figure BDA00029208145700000310
表示状态一下的i方向位移分量，x1表示内点x的分量，Γ表示相互作用积分的路径，Aε是以裂尖为圆心、半径为R的圆，
Figure BDA00029208145700000311
表示状态一下的应变，
Figure BDA00029208145700000312
表示状态二下的应变，
Figure BDA00029208145700000313
表示状态一下的应力，
Figure BDA00029208145700000314
表示状态二下的应力；状态一表示机械结构的真实状态，状态二表示辅助状态；where M (1,2) represents the interaction integral,
Figure BDA0002920814570000037
is the interaction strain energy density,
Figure BDA0002920814570000038
Figure BDA0002920814570000039
represents the i-direction displacement component in state 2,
Figure BDA00029208145700000310
represents the displacement component in the i-direction of the state below, x 1 represents the component of the inner point x, Γ represents the path of the interaction integral, A ε is the circle with the crack tip as the center and the radius as R,
Figure BDA00029208145700000311
represents the state of the strain,
Figure BDA00029208145700000312
represents the strain in state two,
Figure BDA00029208145700000313
represents the stress at one state,
Figure BDA00029208145700000314
Represents the stress in state 2; state 1 represents the real state of the mechanical structure, and state 2 represents the auxiliary state;

S2.5、将公式(7)转化为：S2.5. Convert formula (7) into:

Figure BDA00029208145700000315
Figure BDA00029208145700000315

其中，E为杨氏模型，

Figure BDA00029208145700000316
表示状态一下的I型应力强度因子，
Figure BDA00029208145700000317
表示表示状态二下的I型应力强度因子，
Figure BDA00029208145700000318
表示状态一下的II型应力强度因子，
Figure BDA00029208145700000319
表示状态二下的II型应力强度因子，
Figure BDA00029208145700000320
表示状态一下的III型应力强度因子，
Figure BDA00029208145700000321
表示状态二下的III型应力强度因子；Among them, E is Young's model,
Figure BDA00029208145700000316
represents the I-type stress intensity factor for the first state,
Figure BDA00029208145700000317
represents the type I stress intensity factor in state two,
Figure BDA00029208145700000318
represents the type II stress intensity factor for the first state,
Figure BDA00029208145700000319
represents the type II stress intensity factor in state two,
Figure BDA00029208145700000320
represents the type III stress intensity factor for the first state,
Figure BDA00029208145700000321
represents the type III stress intensity factor in state two;

S2.6、计算混合模态应力强度因子△Keq：S2.6. Calculate the mixed mode stress intensity factor △K eq :

Figure BDA00029208145700000322
Figure BDA00029208145700000322

其中，KI表示I型应力强度因子，KII表示II型应力强度因子；Among them, K I represents the type I stress intensity factor, and K II represents the type II stress intensity factor;

S2.7、采用最大周向应力准则确定裂纹方向：S2.7. Use the maximum circumferential stress criterion to determine the crack direction:

Figure BDA00029208145700000323
Figure BDA00029208145700000323

其中，θ为裂纹扩展角，sign表示符号函数；Among them, θ is the crack propagation angle, and sign represents the sign function;

S2.8、根据混合模态应力强度因子△Keq确定裂纹扩展速率，对于Paris模型，裂纹扩展速率为：S2.8. Determine the crack growth rate according to the mixed mode stress intensity factor △K eq . For the Paris model, the crack growth rate is:

Figure BDA0002920814570000041
Figure BDA0002920814570000041

其中，a为裂纹长度，N为应力循环次数；where a is the crack length and N is the number of stress cycles;

S2.9、根据裂纹的初始长度a0，裂纹的临界长度L，在裂纹扩展步设定裂纹扩展增量△a，可知扩展总步数为

Figure BDA0002920814570000042
在裂纹扩展过程中，针对第t步中△Keq，计算相应的载荷循环次数△Nt：S2.9. According to the initial length a 0 of the crack and the critical length L of the crack, set the crack growth increment Δa in the crack growth step, it can be known that the total number of growth steps is
Figure BDA0002920814570000042
During the crack propagation process, for the ΔK eq in the t-th step, calculate the corresponding number of load cycles ΔN t :

Figure BDA0002920814570000043
Figure BDA0002920814570000043

其中，1≤t≤T；Among them, 1≤t≤T;

S2.10、将载荷循环次数△NT作为机械结构的实时疲劳寿命信息。S2.10. The number of load cycles ΔNT is used as the real-time fatigue life information of the mechanical structure.

所述位移边界积分方程包括非裂纹边界上对应的位移边界积分方程和裂纹上边界对应的位移边界积分方程；所述面力边界积分方程是指裂纹下边界对应的面力边界积分方程；The displacement boundary integral equation includes the displacement boundary integral equation corresponding to the non-crack boundary and the displacement boundary integral equation corresponding to the crack upper boundary; the surface force boundary integral equation refers to the surface force boundary integral equation corresponding to the crack lower boundary;

S2.1.1、构建源点落在非裂纹边界上对应的位移边界积分方程：S2.1.1. Construct the displacement boundary integral equation corresponding to the source point falling on the non-crack boundary:

Figure BDA0002920814570000044
Figure BDA0002920814570000044

其中，

Figure BDA0002920814570000045
表示非裂纹边界上的源点，x表示非裂纹边界上的场点，x+表示裂纹上边界的场点，x-表示裂纹下边界的场点，
Figure BDA0002920814570000046
表示非裂纹边界形状系数，
Figure BDA0002920814570000047
Figure BDA0002920814570000048
表示场点落在非裂纹边界的位移基本解，
Figure BDA0002920814570000049
表示场点落在非裂纹边界的面力基本解，
Figure BDA00029208145700000410
表示场点落在裂纹上边界的位移基本解，
Figure BDA00029208145700000411
表示场点落在裂纹下边界的位移基本解，
Figure BDA00029208145700000412
表示场点落在裂纹上边界的面力基本解，
Figure BDA00029208145700000413
表示场点落在裂纹下边界的面力基本解，
Figure BDA00029208145700000414
S表示非裂纹边界，
Figure BDA00029208145700000415
表示裂纹上边界，
Figure BDA00029208145700000416
表示裂纹下边界，
Figure BDA00029208145700000417
表示非裂纹边界上的源点处的j方向位移分量，uj(x)表示非裂纹边界上场点的j方向位移分量，tj(x)表示非裂纹边界上的场点处的j方向面力分量，uj(x+)表示裂纹上边界的场点处的j方向位移分量，uj(x-)表示裂纹下边界的场点处的j方向位移分量，tj(x+)表示裂纹上边界处的j方向面力分量，tj(x-)表示裂纹下边界处的j方向面力分量，
Figure BDA0002920814570000051
表示位移基本解，
Figure BDA0002920814570000052
表示面力基本解，
Figure BDA0002920814570000053
场点y'为x、x+或x-，
Figure BDA0002920814570000054
表示delta函数，v表示泊松比，r表示源点和场点的距离，n表示法向量，ni表示法向量i方向的分量，nj表示法向量j方向的分量，
Figure BDA0002920814570000055
表示源点
Figure BDA00029208145700000523
的i方向分量，y'i表示场点y'的i方向分量，
Figure BDA0002920814570000056
表示源点
Figure BDA00029208145700000524
的j方向分量，y'j表示场点y'的j方向分量；in,
Figure BDA0002920814570000045
represents the source point on the non-crack boundary, x represents the field point on the non-crack boundary, x + represents the field point on the upper boundary of the crack, x - represents the field point on the lower boundary of the crack,
Figure BDA0002920814570000046
represents the non-crack boundary shape factor,
Figure BDA0002920814570000047
Figure BDA0002920814570000048
represents the basic solution of the displacement where the field point falls on the non-crack boundary,
Figure BDA0002920814570000049
represents the fundamental solution of the surface force where the field point falls on the non-crack boundary,
Figure BDA00029208145700000410
represents the basic solution of the displacement where the field point falls on the upper boundary of the crack,
Figure BDA00029208145700000411
represents the basic solution of the displacement where the field point falls on the lower boundary of the crack,
Figure BDA00029208145700000412
represents the fundamental solution of the surface force where the field point falls on the upper boundary of the crack,
Figure BDA00029208145700000413
represents the fundamental solution of the surface force where the field point falls on the lower boundary of the crack,
Figure BDA00029208145700000414
S represents the non-crack boundary,
Figure BDA00029208145700000415
represents the upper boundary of the crack,
Figure BDA00029208145700000416
represents the lower boundary of the crack,
Figure BDA00029208145700000417
represents the j-direction displacement component at the source point on the non-crack boundary, u j (x) represents the j-direction displacement component at the field point on the non-crack boundary, and t j (x) represents the j-direction surface at the field point on the non-crack boundary Force component, u j (x + ) represents the j-direction displacement component at the field point on the upper crack boundary, u j (x - ) represents the j-direction displacement component at the field point on the crack lower boundary, t j (x + ) represents The j-direction force component at the upper boundary of the crack, t j (x - ) represents the j-direction force component at the lower boundary of the crack,
Figure BDA0002920814570000051
represents the displacement fundamental solution,
Figure BDA0002920814570000052
represents the basic solution of surface force,
Figure BDA0002920814570000053
Field point y' is x, x + or x - ,
Figure BDA0002920814570000054
represents the delta function, v represents the Poisson’s ratio, r represents the distance between the source point and the field point, n represents the normal vector, n i represents the component of the normal vector i in the direction, n j represents the component of the normal vector in the j direction,
Figure BDA0002920814570000055
Indicates the source point
Figure BDA00029208145700000523
The i-direction component of y' i represents the i-direction component of the field point y',
Figure BDA0002920814570000056
Indicates the source point
Figure BDA00029208145700000524
The j-direction component of y' j represents the j-direction component of the field point y';

对公式(1)进行整理合并，可得：Combining formula (1), we can get:

Figure BDA0002920814570000057
Figure BDA0002920814570000057

S2.1.2、构建源点落在裂纹上边界对应的位移边界积分方程：S2.1.2. Construct the displacement boundary integral equation corresponding to the source point falling on the upper boundary of the crack:

Figure BDA0002920814570000058
Figure BDA0002920814570000058

其中，

Figure BDA0002920814570000059
表示落在裂纹上边界的源点，
Figure BDA00029208145700000510
表示落在裂纹下边界的源点，且
Figure BDA00029208145700000511
表示两个源点的坐标相同，
Figure BDA00029208145700000512
表示裂纹上边界形状系数，
Figure BDA00029208145700000513
表示裂纹下边界形状系数，
Figure BDA00029208145700000514
表示裂纹上边界的源点处的j方向位移分量，
Figure BDA00029208145700000515
表示表示裂纹上边界的源点处的j方向位移分量，
Figure BDA00029208145700000516
表示场点落在非裂纹边界的位移基本解，
Figure BDA00029208145700000517
表示场点落在非裂纹边界的面力基本解，
Figure BDA00029208145700000518
表示场点落在裂纹下边界的位移基本解，
Figure BDA00029208145700000519
表示场点落在裂纹上边界的面力基本解；in,
Figure BDA0002920814570000059
represents the source point that falls on the upper boundary of the crack,
Figure BDA00029208145700000510
represents the source point that falls on the lower boundary of the crack, and
Figure BDA00029208145700000511
means that the coordinates of the two source points are the same,
Figure BDA00029208145700000512
is the shape factor of the upper boundary of the crack,
Figure BDA00029208145700000513
represents the shape factor of the lower boundary of the crack,
Figure BDA00029208145700000514
represents the j-direction displacement component at the source point of the upper boundary of the crack,
Figure BDA00029208145700000515
represents the j-direction displacement component at the source point representing the upper boundary of the crack,
Figure BDA00029208145700000516
represents the basic solution of the displacement where the field point falls on the non-crack boundary,
Figure BDA00029208145700000517
represents the fundamental solution of the surface force where the field point falls on the non-crack boundary,
Figure BDA00029208145700000518
represents the basic solution of the displacement where the field point falls on the lower boundary of the crack,
Figure BDA00029208145700000519
Represents the fundamental solution of the surface force where the field point falls on the upper boundary of the crack;

S2.1.3、构建源点落在裂纹下边界对应的面力边界积分方程：S2.1.3. Construct the surface force boundary integral equation corresponding to the source point falling on the lower boundary of the crack:

Figure BDA00029208145700000520
Figure BDA00029208145700000520

其中，

Figure BDA00029208145700000521
表示裂纹下边界的源点处的j方向面力分量，
Figure BDA00029208145700000522
裂纹上边界的源点处的j方向面力分量，
Figure BDA0002920814570000061
表示裂纹下边界的源点处的i方向分量，tk(x)表示非裂纹边界上的场点处的k方向面力分量，
Figure BDA0002920814570000062
表示场点落在非裂纹边界的高阶位移基本解，
Figure BDA0002920814570000063
表示场点落在非裂纹边界的高阶面力基本解，uk(x)表示非裂纹边界上场点的k方向位移分量，tk(x+)表示裂纹上边界场点处的k方向面力分量，tk(x-)表示裂纹下边界处的k方向面力分量，uk(x+)表示裂纹上边界场点处的k方向位移分量，uk(x-)表示裂纹下边界场点处的k方向位移分量，
Figure BDA0002920814570000064
G'表示剪切模量，
Figure BDA0002920814570000065
表示源点
Figure BDA0002920814570000066
的k方向分量，y'k表示场点y'的k方向分量，nk法向量k方向的分量，
Figure BDA0002920814570000067
Figure BDA0002920814570000068
in,
Figure BDA00029208145700000521
represents the j-direction force component at the source point of the lower boundary of the crack,
Figure BDA00029208145700000522
the j-direction force component at the source point of the upper boundary of the crack,
Figure BDA0002920814570000061
is the i-direction component at the source point of the lower crack boundary, t k (x) is the k-direction force component at the field point on the non-crack boundary,
Figure BDA0002920814570000062
represents the fundamental solution for higher-order displacements where the field point falls on the non-crack boundary,
Figure BDA0002920814570000063
Represents the fundamental solution of the higher-order surface force when the field point falls on the non-crack boundary, u k (x) represents the k-direction displacement component of the field point on the non-crack boundary, and t k (x + ) represents the k-direction surface at the field point on the crack upper boundary Force component, t k (x - ) is the force component in the k direction at the lower boundary of the crack, uk (x + ) is the displacement component in the k direction at the field point on the upper boundary of the crack, and u k (x - ) is the lower boundary of the crack the k-direction displacement component at the field point,
Figure BDA0002920814570000064
G' is the shear modulus,
Figure BDA0002920814570000065
Indicates the source point
Figure BDA0002920814570000066
The k-direction component of y' k represents the k-direction component of the field point y', n k is the component of the k-direction normal vector,
Figure BDA0002920814570000067
Figure BDA0002920814570000068

在裂纹的扩展过程中，当前裂纹长度为a0，设置临界长度为L，每一步扩展量为△a，共计扩展

Figure BDA0002920814570000069
步，对公式
Figure BDA00029208145700000610
进行变形后进行两端积分，得到第一步扩展的载荷循环次数△N1：In the process of crack propagation, the current crack length is a 0 , the critical length is set to L, the expansion amount of each step is △a, and the total expansion
Figure BDA0002920814570000069
step, to the formula
Figure BDA00029208145700000610
After the deformation is carried out, the two ends are integrated to obtain the number of load cycles △N 1 for the first expansion:

Figure BDA00029208145700000611
Figure BDA00029208145700000611

进行第二步扩展，得到第二步扩展的载荷循环次数△N2：Carry out the second-step expansion to obtain the load cycle times ΔN 2 of the second-step expansion:

Figure BDA00029208145700000612
Figure BDA00029208145700000612

进行第三步扩展，得到第三步扩展的载荷循环次数△N3：Carry out the third-step expansion, and obtain the load cycle times △N 3 of the third-step expansion:

Figure BDA00029208145700000613
Figure BDA00029208145700000613

同样的第四步扩展一直到第T-1步扩展，相应的载荷循环次数为△N4...△NT-1；The same fourth step is extended until the T-1 step is extended, and the corresponding load cycle times are ΔN 4 ... ΔN T-1 ;

进行第T步扩展，得到得到第T步扩展的载荷循环次数△NT：Carry out the T-th step expansion, and obtain the number of load cycles △N T for the T-th step expansion:

Figure BDA00029208145700000614
Figure BDA00029208145700000614

因此，机械结构的总寿命为△N＝△N1+△N2+...+△NT-1+△NT。Therefore, the total life of the mechanical structure is ΔN=ΔN 1 +ΔN 2 +...+ΔN T -1 +ΔNT .

所述将训练集输入BP神经网络中进行训练，得到BP神经网络模型的方法为：The described training set is input into the BP neural network for training, and the method for obtaining the BP neural network model is:

S4.1、网络初始化：确定输入层、隐含层和输出层的节点数，初始化连接权值Wi'j'和Wj'k'以及隐含层和输出层的阈值，i'＝1,2,…,n'，n'为输入层神经元数，j'＝1,2,…,l，l为隐含层神经元数，k'＝1,2,…,m'，m'为输出层神经元数；S4.1. Network initialization: determine the number of nodes in the input layer, the hidden layer and the output layer, initialize the connection weights W i'j' and W j'k' and the thresholds of the hidden layer and the output layer, i'=1 ,2,...,n',n' is the number of neurons in the input layer, j'=1,2,...,l,l is the number of neurons in the hidden layer, k'=1,2,...,m',m ' is the number of neurons in the output layer;

S4.2、计算隐含层节点的输出：S4.2. Calculate the output of hidden layer nodes:

Figure BDA0002920814570000071
Figure BDA0002920814570000071

其中，hj'表示神经元j'的输出，xi'表示输入分量，ωi'j'∈Wi'j'表示加权权值，θj'为神经元j'的阈值，f(·)为激活函数；Among them, h j' represents the output of neuron j', xi' represents the input component, ω i'j' ∈W i'j' represents the weighted weight, θ j' is the threshold of neuron j', f(· ) is the activation function;

S4.3、计算输出层节点的输出：S4.3. Calculate the output of the output layer node:

Figure BDA0002920814570000072
Figure BDA0002920814570000072

其中，yk'表示神经元k'的输出，ωj'k'∈Wj'k'表示加权权值，θk'为神经元k'的阈值；Among them, y k' represents the output of neuron k', ω j'k' ∈ W j'k' represents the weighted weight, and θ k' is the threshold of neuron k';

S4.4、计算第z个训练数据的预期输出与预测输出之间的误差ez：S4.4. Calculate the error ez between the expected output of the zth training data and the predicted output:

Figure BDA0002920814570000073
Figure BDA0002920814570000073

其中，yd,k'为神经元k'在输出层的期望输出值，yk'为神经元k'在输出层的预测输出值，则Z个输入信号的总误差为E'：Among them, y d, k' is the expected output value of the neuron k' in the output layer, y k' is the predicted output value of the neuron k' in the output layer, then the total error of the Z input signals is E':

Figure BDA0002920814570000074
Figure BDA0002920814570000074

S4.5、更新神经网络的连接权值：S4.5, update the connection weights of the neural network:

ωi'j'(p+1)＝ωi'j'(p)+△ωi'j'(p) (21)；ω i'j' (p+1)=ω i'j' (p)+△ω i'j' (p) (21);

其中，p为迭代次数，△ωi'j'(p)＝η×yi'(p)×λj'(p)为p+1次迭代中连接权值的调整部分，η为学习率，yi'(p)为神经元i'的输出信号，

Figure BDA0002920814570000075
为神经元j'的误差梯度，Xj'(p)为神经元j'的加权输入，δj'(p)为神经元j'的误差；Among them, p is the number of iterations, △ω i'j' (p)=η×y i' (p)×λ j' (p) is the adjustment part of the connection weight in p+1 iterations, and η is the learning rate , y i' (p) is the output signal of neuron i',
Figure BDA0002920814570000075
is the error gradient of neuron j', X j' (p) is the weighted input of neuron j', δ j' (p) is the error of neuron j';

S4.6、判断迭代次数是否达到最大迭代次数或者Z个输入信号的总误差E'小于误差阈值，若是，输出BP神经网络模型，否则，返回步骤S4.2进行下一次的迭代。S4.6, determine whether the number of iterations reaches the maximum number of iterations or the total error E' of the Z input signals is less than the error threshold, if so, output the BP neural network model, otherwise, return to step S4.2 for the next iteration.

本技术方案能产生的有益效果：本发明利用对偶互易边界元法生成一系列结构响应和相应的疲劳寿命数据集，并利用这些数据集对BP神经网络预测模型进行训练和测试，建立关于观测点位置的结构响应与疲劳寿命关联的实时疲劳寿命预测模型，在实际工程中通过对机械结构的观测点位置的响应试验测试来确定实时疲劳寿命信息，避免了检测精度降低的问题，节省大量的时间和成本。The beneficial effects that the technical solution can produce: the present invention uses the dual reciprocal boundary element method to generate a series of structural response and corresponding fatigue life data sets, and uses these data sets to train and test the BP neural network prediction model, and establishes a data set about the observation The real-time fatigue life prediction model that correlates the structural response of the point position with the fatigue life. In actual engineering, the real-time fatigue life information is determined by the response test of the observation point position of the mechanical structure, which avoids the problem of reducing the detection accuracy and saves a lot of time. time and cost.

附图说明Description of drawings

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。In order to explain the embodiments of the present invention or the technical solutions in the prior art more clearly, the following briefly introduces the accompanying drawings that need to be used in the description of the embodiments or the prior art. Obviously, the accompanying drawings in the following description are only These are some embodiments of the present invention. For those of ordinary skill in the art, other drawings can also be obtained according to these drawings without creative efforts.

图1为本发明的BP神经网络拓扑结构。Fig. 1 is the topology structure of the BP neural network of the present invention.

图2为本发明的流程图。Figure 2 is a flow chart of the present invention.

具体实施方式Detailed ways

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有付出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, but not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative efforts shall fall within the protection scope of the present invention.

如图2所示，本发明实施例提供了一种基于数据驱动的机械结构实时疲劳寿命预测方法，由于机械结构中对应的Paris模型指数m和系数c不确定，通过试验手段得出取值范围，在此取值范围内，按照优先数系生产关于m和c的一系列数据，并进行配对(m，c)。针对该机械结构，在容易测试部位取观测点，分别取每一组(m，c)进行循环，使用对偶互易边界元法进行疲劳裂纹扩展分析，每一步扩展都会获得观测点的位移值以及相应的实时疲劳寿命信息，组成数据信息对，存储起来，保存到一个单元中。直至所有的(m，c)循环完毕，将单元分为训练单元和测试单元。然后利用训练单元对BP神经网络模型进行训练。利用测试单元对训练好的BP神经网络模型进行验证。这样就完全建立起来观测点位移与实时疲劳寿命信息的关系，因此，可以仅通过观测点位移，利用建立好的BP神经网络模型获得相应的实时疲劳寿命信息。具体步骤如下：As shown in FIG. 2 , an embodiment of the present invention provides a data-driven real-time fatigue life prediction method for mechanical structures. Since the corresponding Paris model index m and coefficient c in the mechanical structure are uncertain, the value range is obtained by means of experiments. , within this value range, a series of data about m and c are produced according to the priority number system, and paired (m, c). For this mechanical structure, the observation points are taken at the easy-to-test parts, and each group (m, c) is taken for the cycle, and the dual reciprocal boundary element method is used to analyze the fatigue crack growth. The displacement value of the observation point and the The corresponding real-time fatigue life information is composed of data information pairs, which are stored and saved in a unit. Until all (m, c) cycles are completed, the units are divided into training units and test units. Then use the training unit to train the BP neural network model. Use the test unit to verify the trained BP neural network model. In this way, the relationship between the displacement of the observation point and the real-time fatigue life information is completely established. Therefore, the corresponding real-time fatigue life information can be obtained by using the established BP neural network model only through the displacement of the observation point. Specific steps are as follows:

步骤一：根据试验手段获得机械结构中对应的Paris模型的指数m和系数c的取值范围，并分别在指数m和系数c的取值范围内采用优先数系生成一系列指数和系数，得到Q组Paris模型参数对(mq，cq)，其中，q＝1,2,…,Q；在实际机械结构中，Paris模型指数m和系数c不确定，很难直接模拟疲劳裂纹扩展。由于摄动法需要先验知识，模型参数的分布函数是未知的。对于特定的机械结构，需要进行试验测量，并对数据进行拟合，根据拟合获得的Paris模型指数m和系数c进行疲劳裂纹扩展分析，预测疲劳寿命。然而对待大量的机械结构，往往缺乏时间和成本，难以进行批次试验拟合。因此，本发明实施例开发了一种基于数据驱动的机械结构实时疲劳寿命预测边界元法。Step 1: Obtain the value range of the exponent m and coefficient c of the Paris model corresponding to the mechanical structure according to the test method, and use the priority number system to generate a series of exponents and coefficients within the value range of the exponent m and coefficient c respectively, and obtain Q group of Paris model parameter pairs (m q , c q ), where q=1, 2, . Since the perturbation method requires prior knowledge, the distribution functions of the model parameters are unknown. For a specific mechanical structure, it is necessary to carry out experimental measurements and fit the data. According to the Paris model exponent m and coefficient c obtained by fitting, the fatigue crack growth analysis is carried out to predict the fatigue life. However, when dealing with a large number of mechanical structures, it is often lack of time and cost, and it is difficult to perform batch test fitting. Therefore, an embodiment of the present invention develops a data-driven boundary element method for real-time fatigue life prediction of mechanical structures.

步骤二：在机械结构上设定观测点，针对第q组Paris模型参数对(mq，cq)进行循环，使用对偶互易边界元法对机械结构进行疲劳裂纹扩展分析，获得观测点的位移和机械结构实时疲劳寿命信息，并将观测点的位移和机械结构实时疲劳寿命信息组成数据信息对。Step 2: Set the observation point on the mechanical structure, cycle (m q , c q ) for the qth group of Paris model parameters, and use the dual reciprocal boundary element method to analyze the fatigue crack growth of the mechanical structure to obtain the observation point. Displacement and real-time fatigue life information of mechanical structure, and the displacement of the observation point and real-time fatigue life information of mechanical structure form a data information pair.

S2.1、构建对偶互易边界积分方程，包括位移边界积分方程和面力边界积分方程；所述位移边界积分方程包括非裂纹边界上对应的位移边界积分方程和裂纹上边界对应的位移边界积分方程；所述面力边界积分方程是指裂纹下边界对应的面力边界积分方程。S2.1. Construct the dual reciprocal boundary integral equation, including the displacement boundary integral equation and the surface force boundary integral equation; the displacement boundary integral equation includes the displacement boundary integral equation corresponding to the non-crack boundary and the displacement boundary integral corresponding to the crack upper boundary Equation; the surface force boundary integral equation refers to the surface force boundary integral equation corresponding to the lower boundary of the crack.

S2.1.1、构建源点落在非裂纹边界上对应的位移边界积分方程：S2.1.1. Construct the displacement boundary integral equation corresponding to the source point falling on the non-crack boundary:

Figure BDA0002920814570000091
Figure BDA0002920814570000091

其中，

Figure BDA0002920814570000092
表示非裂纹边界上的源点，x表示非裂纹边界上的场点，x+表示裂纹上边界的场点，x-表示裂纹下边界的场点，
Figure BDA0002920814570000093
表示非裂纹边界形状系数，
Figure BDA0002920814570000094
Figure BDA0002920814570000095
表示场点落在非裂纹边界的位移基本解，
Figure BDA0002920814570000096
表示场点落在非裂纹边界的面力基本解，
Figure BDA0002920814570000097
表示场点落在裂纹上边界的位移基本解，
Figure BDA0002920814570000098
表示场点落在裂纹下边界的位移基本解，
Figure BDA0002920814570000099
表示场点落在裂纹上边界的面力基本解，
Figure BDA00029208145700000910
表示场点落在裂纹下边界的面力基本解，
Figure BDA00029208145700000911
S表示非裂纹边界，
Figure BDA00029208145700000912
表示裂纹上边界，
Figure BDA00029208145700000913
表示裂纹下边界，
Figure BDA00029208145700000914
表示非裂纹边界上的源点处的j方向位移分量，uj(x)表示非裂纹边界上场点的j方向位移分量，tj(x)表示非裂纹边界上的场点处的j方向面力分量，uj(x+)表示裂纹上边界的场点处的j方向位移分量，uj(x-)表示裂纹下边界的场点处的j方向位移分量，tj(x+)表示裂纹上边界处的j方向面力分量，tj(x-)表示裂纹下边界处的j方向面力分量，
Figure BDA00029208145700000915
表示位移基本解，
Figure BDA00029208145700000916
表示面力基本解，
Figure BDA00029208145700000917
场点y'为x、x+或x-，
Figure BDA00029208145700000918
表示delta函数，v表示泊松比，r表示源点和场点的距离，n表示法向量，ni表示法向量i方向的分量，nj表示法向量j方向的分量，
Figure BDA0002920814570000101
表示源点
Figure BDA0002920814570000102
的i方向分量，y'i表示场点y'的i方向分量，场点y'为x、x+或x-，
Figure BDA0002920814570000103
表示源点
Figure BDA00029208145700001023
的j方向分量，y'j表示场点y'的j方向分量。in,
Figure BDA0002920814570000092
represents the source point on the non-crack boundary, x represents the field point on the non-crack boundary, x + represents the field point on the upper boundary of the crack, x - represents the field point on the lower boundary of the crack,
Figure BDA0002920814570000093
represents the non-crack boundary shape factor,
Figure BDA0002920814570000094
Figure BDA0002920814570000095
represents the basic solution of the displacement where the field point falls on the non-crack boundary,
Figure BDA0002920814570000096
represents the fundamental solution of the surface force where the field point falls on the non-crack boundary,
Figure BDA0002920814570000097
represents the basic solution of the displacement where the field point falls on the upper boundary of the crack,
Figure BDA0002920814570000098
represents the basic solution of the displacement where the field point falls on the lower boundary of the crack,
Figure BDA0002920814570000099
represents the fundamental solution of the surface force where the field point falls on the upper boundary of the crack,
Figure BDA00029208145700000910
represents the fundamental solution of the surface force where the field point falls on the lower boundary of the crack,
Figure BDA00029208145700000911
S represents the non-crack boundary,
Figure BDA00029208145700000912
represents the upper boundary of the crack,
Figure BDA00029208145700000913
represents the lower boundary of the crack,
Figure BDA00029208145700000914
represents the j-direction displacement component at the source point on the non-crack boundary, u j (x) represents the j-direction displacement component at the field point on the non-crack boundary, and t j (x) represents the j-direction surface at the field point on the non-crack boundary Force component, u j (x + ) represents the j-direction displacement component at the field point on the upper crack boundary, u j (x - ) represents the j-direction displacement component at the field point on the crack lower boundary, t j (x + ) represents The j-direction force component at the upper boundary of the crack, t j (x - ) represents the j-direction force component at the lower boundary of the crack,
Figure BDA00029208145700000915
represents the displacement fundamental solution,
Figure BDA00029208145700000916
represents the basic solution of surface force,
Figure BDA00029208145700000917
Field point y' is x, x + or x - ,
Figure BDA00029208145700000918
represents the delta function, v represents the Poisson’s ratio, r represents the distance between the source point and the field point, n represents the normal vector, n i represents the component of the normal vector i in the direction, n j represents the component of the normal vector in the j direction,
Figure BDA0002920814570000101
Indicates the source point
Figure BDA0002920814570000102
The i direction component of y' i represents the i direction component of the field point y', and the field point y' is x, x + or x - ,
Figure BDA0002920814570000103
Indicates the source point
Figure BDA00029208145700001023
The j-direction component of y' j represents the j-direction component of the field point y'.

对公式(1)进行整理合并，可得：Combining formula (1), we can get:

Figure BDA0002920814570000104
Figure BDA0002920814570000104

S2.1.2、构建源点落在裂纹上边界对应的位移边界积分方程：S2.1.2. Construct the displacement boundary integral equation corresponding to the source point falling on the upper boundary of the crack:

Figure BDA0002920814570000105
Figure BDA0002920814570000105

其中，

Figure BDA0002920814570000106
表示落在裂纹上边界的源点，
Figure BDA0002920814570000107
表示落在裂纹下边界的源点，且
Figure BDA0002920814570000108
(由于裂纹上下面重合，所以两个源点的坐标相同)，
Figure BDA0002920814570000109
表示裂纹上边界形状系数，
Figure BDA00029208145700001010
表示裂纹下边界形状系数，
Figure BDA00029208145700001011
表示裂纹上边界的源点处的j方向位移分量，
Figure BDA00029208145700001012
表示表示裂纹上边界的源点处的j方向位移分量，
Figure BDA00029208145700001013
表示场点落在非裂纹边界的位移基本解，
Figure BDA00029208145700001014
表示场点落在非裂纹边界的面力基本解，
Figure BDA00029208145700001015
表示场点落在裂纹下边界的位移基本解，
Figure BDA00029208145700001016
表示场点落在裂纹上边界的面力基本解。in,
Figure BDA0002920814570000106
represents the source point that falls on the upper boundary of the crack,
Figure BDA0002920814570000107
represents the source point that falls on the lower boundary of the crack, and
Figure BDA0002920814570000108
(Because the top and bottom of the crack overlap, the coordinates of the two source points are the same),
Figure BDA0002920814570000109
is the shape factor of the upper boundary of the crack,
Figure BDA00029208145700001010
represents the shape factor of the lower boundary of the crack,
Figure BDA00029208145700001011
represents the j-direction displacement component at the source point of the upper boundary of the crack,
Figure BDA00029208145700001012
represents the j-direction displacement component at the source point representing the upper boundary of the crack,
Figure BDA00029208145700001013
represents the basic solution of the displacement where the field point falls on the non-crack boundary,
Figure BDA00029208145700001014
represents the fundamental solution of the surface force where the field point falls on the non-crack boundary,
Figure BDA00029208145700001015
represents the basic solution of the displacement where the field point falls on the lower boundary of the crack,
Figure BDA00029208145700001016
Represents the fundamental solution of the surface force where the field point falls on the upper boundary of the crack.

S2.1.3、构建源点落在裂纹下边界对应的面力边界积分方程：S2.1.3. Construct the surface force boundary integral equation corresponding to the source point falling on the lower boundary of the crack:

Figure BDA00029208145700001017
Figure BDA00029208145700001017

其中，

Figure BDA00029208145700001018
表示裂纹下边界的源点处的j方向面力分量，
Figure BDA00029208145700001019
裂纹上边界的源点处的j方向面力分量，
Figure BDA00029208145700001020
表示裂纹下边界的源点处的i方向分量，tk(x)表示非裂纹边界上的场点处的k方向面力分量，
Figure BDA00029208145700001021
表示场点落在非裂纹边界的高阶位移基本解，
Figure BDA00029208145700001022
表示场点落在非裂纹边界的高阶面力基本解，uk(x)表示非裂纹边界上场点的k方向位移分量，tk(x+)表示裂纹上边界场点处的k方向面力分量，tk(x-)表示裂纹下边界处的k方向面力分量，uk(x+)表示裂纹上边界场点处的k方向位移分量，uk(x-)表示裂纹下边界场点处的k方向位移分量，
Figure BDA0002920814570000111
G'表示剪切模量，
Figure BDA0002920814570000112
表示源点
Figure BDA0002920814570000113
的k方向分量，y'k表示场点y'的k方向分量，nk法向量k方向的分量，
Figure BDA0002920814570000114
Figure BDA0002920814570000115
in,
Figure BDA00029208145700001018
represents the j-direction force component at the source point of the lower boundary of the crack,
Figure BDA00029208145700001019
the j-direction force component at the source point of the upper boundary of the crack,
Figure BDA00029208145700001020
is the i-direction component at the source point of the lower crack boundary, t k (x) is the k-direction force component at the field point on the non-crack boundary,
Figure BDA00029208145700001021
represents the fundamental solution for higher-order displacements where the field point falls on the non-crack boundary,
Figure BDA00029208145700001022
Represents the fundamental solution of the higher-order surface force when the field point falls on the non-crack boundary, u k (x) represents the k-direction displacement component of the field point on the non-crack boundary, and t k (x + ) represents the k-direction surface at the field point on the crack upper boundary Force component, t k (x - ) is the force component in the k direction at the lower boundary of the crack, uk (x + ) is the displacement component in the k direction at the field point on the upper boundary of the crack, and u k (x - ) is the lower boundary of the crack the k-direction displacement component at the field point,
Figure BDA0002920814570000111
G' is the shear modulus,
Figure BDA0002920814570000112
Indicates the source point
Figure BDA0002920814570000113
The k-direction component of y' k represents the k-direction component of the field point y', n k is the component of the k-direction normal vector,
Figure BDA0002920814570000114
Figure BDA0002920814570000115

S2.2、分别对位移边界积分方程和面力边界积分方程进行离散处理，得到：S2.2. Discrete the displacement boundary integral equation and the surface force boundary integral equation respectively, and obtain:

Figure BDA0002920814570000116
Figure BDA0002920814570000116

其中，H*表示位移基本解积分，G*表示面力基本解积分，uS表示非裂纹边界上位移，

Figure BDA0002920814570000117
表示裂纹上边界上的位移，
Figure BDA0002920814570000118
表示裂纹下边界上的位移，tS表示非裂纹边界上的面力，
Figure BDA0002920814570000119
表示裂纹上边界上的面力，
Figure BDA00029208145700001110
表示裂纹下边界上的面力。where H * represents the basic solution integral of displacement, G * represents the basic solution integral of surface force, u S represents the displacement on the non-crack boundary,
Figure BDA0002920814570000117
represents the displacement on the upper boundary of the crack,
Figure BDA0002920814570000118
represents the displacement on the lower boundary of the crack, t S represents the surface force on the non-crack boundary,
Figure BDA0002920814570000119
represents the surface force on the upper boundary of the crack,
Figure BDA00029208145700001110
Represents the surface force on the lower boundary of the crack.

S2.3、对公式(5)施加边界条件，将未知量调到方程左边，已知量调到方程的右端，并整理得到：S2.3. Apply boundary conditions to formula (5), adjust the unknown quantity to the left side of the equation, and adjust the known quantity to the right side of the equation, and get:

Figure BDA00029208145700001111
Figure BDA00029208145700001111

其中，A*表示位移或面力的已知量，xS表示非裂纹边界上未知的位移或面力，

Figure BDA00029208145700001112
表示裂纹上边界上未知的位移或面力，
Figure BDA00029208145700001113
表示裂纹下边界上未知的位移或面力，bS表示非裂纹边界上已知的位移或面力合并得到的分量，
Figure BDA00029208145700001114
表示裂纹上边界上已知的位移或面力合并得到的分量，
Figure BDA00029208145700001115
表示裂纹下边界上已知的位移或面力合并得到的分量。where A * represents the known quantity of displacement or surface force, x S represents the unknown displacement or surface force on the non-crack boundary,
Figure BDA00029208145700001112
represents the unknown displacement or surface force on the upper boundary of the crack,
Figure BDA00029208145700001113
represents the unknown displacement or surface force on the lower boundary of the crack, b S represents the combined component of the known displacement or surface force on the non-crack boundary,
Figure BDA00029208145700001114
represents the component obtained by combining the known displacements or surface forces on the upper boundary of the crack,
Figure BDA00029208145700001115
Represents the combined component of known displacements or surface forces on the lower boundary of the crack.

利用高斯消去法对公式(6)求解，获得未知量，并进行未知量配置，获得边界的位移和面力，并将边界的位移和面力代入公式(1)中获得观测点的位移。Using the Gaussian elimination method to solve the formula (6), the unknowns are obtained, and the unknowns are configured to obtain the displacement and surface force of the boundary, and the displacement of the boundary and the surface force are substituted into the formula (1) to obtain the displacement of the observation point.

S2.4、为了评价疲劳裂纹的扩展，需要利用对偶互易边界元法计算应力强度因子，应力强度因子是判断裂纹扩展与否的重要参数。此外，还根据应力强度因子计算了裂纹扩展方向和扩展量。利用相互作用积分算法计算了混合模态应力强度因子，其相互作用积分M表达式为：S2.4. In order to evaluate the propagation of fatigue cracks, the dual reciprocal boundary element method needs to be used to calculate the stress intensity factor. The stress intensity factor is an important parameter for judging whether the crack propagates or not. In addition, the direction and amount of crack propagation were calculated based on the stress intensity factor. The mixed mode stress intensity factor is calculated by the interaction integral algorithm, and the expression of the interaction integral M is:

Figure BDA0002920814570000121
Figure BDA0002920814570000121

其中，M(1,2)表示相互作用积分，

Figure BDA0002920814570000122
为互作用应变能密度，
Figure BDA0002920814570000123
Figure BDA0002920814570000124
表示状态二下的i方向位移分量，
Figure BDA0002920814570000125
表示状态一下的i方向位移分量，x1表示内点x的分量，Γ表示相互作用积分的路径，Aε是以裂尖为圆心、半径为R的圆，
Figure BDA0002920814570000126
表示状态一下的应变，
Figure BDA0002920814570000127
表示状态二下的应变，
Figure BDA0002920814570000128
表示状态一下的应力，
Figure BDA0002920814570000129
表示状态二下的应力；状态一表示机械结构的真实状态，状态二表示辅助状态(状态二表示只有I型裂纹或者II型断裂裂纹情况)；状态二下的位移、应力、应变可以通过只有I型或者II型断裂情况下裂纹尖端的渐近位移场和应力场获得，状态一下的位移、应力、应变可以通过对偶互易边界元法的后处理过程获得。where M (1,2) represents the interaction integral,
Figure BDA0002920814570000122
is the interaction strain energy density,
Figure BDA0002920814570000123
Figure BDA0002920814570000124
represents the i-direction displacement component in state 2,
Figure BDA0002920814570000125
represents the displacement component in the i-direction of the state below, x 1 represents the component of the inner point x, Γ represents the path of the interaction integral, A ε is the circle with the crack tip as the center and the radius as R,
Figure BDA0002920814570000126
represents the state of the strain,
Figure BDA0002920814570000127
represents the strain in state two,
Figure BDA0002920814570000128
represents the stress at one state,
Figure BDA0002920814570000129
Represents the stress in state 2; state 1 represents the real state of the mechanical structure, state 2 represents the auxiliary state (state 2 represents only type I cracks or type II fracture cracks); the displacement, stress and strain in state 2 can pass through only I The asymptotic displacement field and stress field of the crack tip in the case of type or type II fracture can be obtained, and the displacement, stress and strain of the state can be obtained by the post-processing process of the dual reciprocal boundary element method.

S2.5、将公式(7)转化为：S2.5. Convert formula (7) into:

Figure BDA00029208145700001210
Figure BDA00029208145700001210

其中，E为杨氏模型，

Figure BDA00029208145700001211
表示状态一下的I型应力强度因子，
Figure BDA00029208145700001212
表示表示状态二下的I型应力强度因子，
Figure BDA00029208145700001213
表示状态一下的II型应力强度因子，
Figure BDA00029208145700001214
表示状态二下的II型应力强度因子，
Figure BDA00029208145700001215
表示状态一下的III型应力强度因子，
Figure BDA00029208145700001216
表示状态二下的III型应力强度因子。Among them, E is Young's model,
Figure BDA00029208145700001211
represents the I-type stress intensity factor for the first state,
Figure BDA00029208145700001212
represents the type I stress intensity factor in state two,
Figure BDA00029208145700001213
represents the type II stress intensity factor for the first state,
Figure BDA00029208145700001214
represents the type II stress intensity factor in state two,
Figure BDA00029208145700001215
represents the type III stress intensity factor for the first state,
Figure BDA00029208145700001216
represents the type III stress intensity factor in state two.

S2.6、根据机械结构材料的断裂韧性KIC与△Keq关系判断裂纹是否扩展，并计算混合模态应力强度因子△Keq：S2.6. According to the relationship between the fracture toughness K IC and △K eq of mechanical structural materials, determine whether the crack propagates, and calculate the mixed mode stress intensity factor △K eq :

Figure BDA00029208145700001217
Figure BDA00029208145700001217

其中，KI表示I型应力强度因子，KII表示II型应力强度因子。Among them, K I represents the type I stress intensity factor, and K II represents the type II stress intensity factor.

S2.7、在裂纹扩展的情况下，确定裂纹扩展方向。采用最大周向应力准则确定裂纹方向：S2.7. In the case of crack propagation, determine the direction of crack propagation. Use the maximum circumferential stress criterion to determine the crack direction:

Figure BDA00029208145700001218
Figure BDA00029208145700001218

其中，θ为裂纹扩展角，sign表示符号函数。Among them, θ is the crack propagation angle, and sign represents the sign function.

S2.8、用对偶边界元法进行疲劳寿命预测时，需要建立一个描述裂纹扩展过程的疲劳裂纹扩展模型。根据混合模态应力强度因子△Keq确定裂纹扩展速率，对于广义Paris模型，裂纹扩展速率为：S2.8. When using the dual boundary element method for fatigue life prediction, it is necessary to establish a fatigue crack growth model that describes the crack growth process. The crack growth rate is determined according to the mixed-modal stress intensity factor ΔK eq . For the generalized Paris model, the crack growth rate is:

Figure BDA00029208145700001219
Figure BDA00029208145700001219

Figure BDA0002920814570000131
Figure BDA0002920814570000131

其中，a为裂纹长度，N为应力循环次数。m为Paris模型指数，c为Paris模型常数，△Keq是混合模态应力强度因子，会受加载条件和材料性能的影响，m和c需要通过疲劳裂纹扩展率试验数据拟合来获得，往往只能得出m、c的范围，在本发明中，会根据m和c的范围会采用优先数系生成一组特定的参数值生产一系列的随机的数据对(m，c)。where a is the crack length and N is the number of stress cycles. m is the Paris model exponent, c is the Paris model constant, △K eq is the mixed mode stress intensity factor, which will be affected by the loading conditions and material properties, m and c need to be obtained by fitting the fatigue crack growth rate test data, often Only the ranges of m and c can be obtained. In the present invention, a series of random data pairs (m, c) are generated by using a priority number system to generate a set of specific parameter values according to the ranges of m and c.

S2.9、根据裂纹的初始长度a0，裂纹的临界长度L，在裂纹扩展步设定裂纹扩展增量△a，可知扩展总步数为

Figure BDA0002920814570000132
在裂纹扩展过程中，针对第t步中△Keq，计算相应的载荷循环次数△Nt：S2.9. According to the initial length a 0 of the crack and the critical length L of the crack, set the crack growth increment Δa in the crack growth step, it can be known that the total number of growth steps is
Figure BDA0002920814570000132
During the crack propagation process, for the ΔK eq in the t-th step, calculate the corresponding number of load cycles ΔN t :

Figure BDA0002920814570000133
Figure BDA0002920814570000133

其中，1≤t≤T。Among them, 1≤t≤T.

S2.10、将载荷循环次数△NT作为观测点的实时疲劳寿命信息。S2.10. Take the number of load cycles ΔNT as the real-time fatigue life information of the observation point.

在裂纹的扩展过程中，当前裂纹长度为a0，设置临界长度为L(当裂纹扩展长度达到L时失效)，每一步扩展量为△a，共计扩展

Figure BDA0002920814570000134
步，对公式(11b)两端积分，可以得到第一步的载荷循环次数△N1：In the process of crack propagation, the current crack length is a 0 , the critical length is set to L (when the crack propagation length reaches L), and the amount of each step is △a, the total expansion
Figure BDA0002920814570000134
step, and integrating both ends of formula (11b), the number of load cycles ΔN 1 in the first step can be obtained:

Figure BDA0002920814570000135
Figure BDA0002920814570000135

同样的第二步扩展，可以得到：The same second step expansion, you can get:

Figure BDA0002920814570000136
Figure BDA0002920814570000136

第三步扩展，可以得到：The third step is to expand, you can get:

Figure BDA0002920814570000137
Figure BDA0002920814570000137

同样的第四步一直到第T-1步扩展，相应的载荷循环次数△N4...△NT-1；The same fourth step has been extended until step T-1, and the corresponding number of load cycles ΔN 4 ... ΔN T-1 ;

第T步扩展：Step T expansion:

Figure BDA0002920814570000138
Figure BDA0002920814570000138

总的寿命N＝△N1+△N2+...+△NT-1。The total lifetime N=ΔN 1 +ΔN 2 +...+ΔN T-1 .

在非裂纹边界上布置观测点(这些观测点容易用实验仪器测量)，在扩展的每一步，都使用公式(1)计算出观测点的位移，将观测点的位移与此阶段机械结构的疲劳寿命组成数据对存储起来。Arrange observation points on the non-crack boundary (these observation points are easy to measure with experimental instruments), and at each step of expansion, use formula (1) to calculate the displacement of the observation point, and compare the displacement of the observation point with the fatigue of the mechanical structure at this stage. Lifetime composition data pairs are stored.

步骤三：循环执行步骤二，直至遍历Q组Paris模型参数对，得到数据集，并将数据集分为训练集和测试集，其中，数据集包括一系列数据信息对。Step 3: Execute step 2 in a loop until the Q groups of Paris model parameter pairs are traversed to obtain a data set, and the data set is divided into a training set and a test set, wherein the data set includes a series of data information pairs.

步骤四：将训练集输入BP神经网络中进行训练，得到BP神经网络模型，并利用测试集对BP神经网络模型进行验证。为了建立数据驱动的预测模型，首先根据试验得来的数据范围，根据优先数系生成pairs模型指数m和系数c的一系列数据，针对每一个系数对(m,c)，利用对偶互易边界元法计算出来观测点位置的结构响应与疲劳寿命配对的大量数据集，构造一个反向传播神经网络来学习这些数据集。反向传播神经网络是一种基于误差反向传播技术的多层前向神经网络由一个输入层、一个或多个隐藏层和一个输出层组成。Step 4: Input the training set into the BP neural network for training, obtain the BP neural network model, and use the test set to verify the BP neural network model. In order to establish a data-driven prediction model, first, according to the data range obtained from the experiment, a series of data of pairs model index m and coefficient c are generated according to the priority number system, and for each coefficient pair (m, c), the dual reciprocity boundary is used. The meta-method computes a large dataset of structural responses paired with fatigue life at observation point locations, and constructs a back-propagation neural network to learn from these datasets. Back-propagation neural network is a multi-layer forward neural network based on error back-propagation technology, which consists of an input layer, one or more hidden layers and an output layer.

图1中的圆圈称为偏置节点。利用数据集对BP神经网络进行训练，具体步骤如下：The circles in Figure 1 are called bias nodes. Using the dataset to train the BP neural network, the specific steps are as follows:

S4.1、网络初始化：确定输入层、隐含层和输出层的节点数，初始化连接权值Wi'j'和Wj'k'，以及隐含层和输出层的阈值，i'＝1,2,…,n'，n'为输入层神经元数，j'＝1,2,…,l，l为隐含层神经元数，k'＝1,2,…,m'，m'为输出层神经元数；S4.1, network initialization: determine the number of nodes in the input layer, the hidden layer and the output layer, initialize the connection weights W i'j' and W j'k' , and the thresholds of the hidden layer and the output layer, i'= 1,2,...,n', n' is the number of neurons in the input layer, j'=1,2,...,l, l is the number of neurons in the hidden layer, k'=1,2,...,m', m' is the number of neurons in the output layer;

S4.2、计算隐含层节点的输出：S4.2. Calculate the output of the hidden layer node:

Figure BDA0002920814570000141
Figure BDA0002920814570000141

其中，hj'表示神经元j'的输出，xi'表示输入分量，ωi'j'∈Wi'j'表示加权权值，θj'为神经元j'的阈值，f(·)为激活函数；Among them, h j' represents the output of neuron j', xi' represents the input component, ω i'j' ∈W i'j' represents the weighted weight, θ j' is the threshold of neuron j', f(· ) is the activation function;

S4.3、计算输出层节点的输出：S4.3. Calculate the output of the output layer node:

Figure BDA0002920814570000142
Figure BDA0002920814570000142

其中，yk'表示k'神经元的输出，ωj'k'∈Wj'k'表示加权权值，θk'神经元k'的阈值；Among them, y k' represents the output of k' neuron, ω j'k' ∈ W j'k' represents the weighted weight, and the threshold of θ k' neuron k';

S4.4、计算第z个训练数据的预期输出与预测输出之间的误差ez：S4.4. Calculate the error ez between the expected output of the zth training data and the predicted output:

Figure BDA0002920814570000143
Figure BDA0002920814570000143

其中，yd,k'为神经元k'在输出层的期望输出值，yk'为神经元k'在输出层的预测输出值，则Z个输入信号的总误差为E'：Among them, y d, k' is the expected output value of the neuron k' in the output layer, y k' is the predicted output value of the neuron k' in the output layer, then the total error of the Z input signals is E':

Figure BDA0002920814570000144
Figure BDA0002920814570000144

S4.5、更新神经网络的连接权值：S4.5, update the connection weights of the neural network:

ωi'j'(p+1)＝ωi'j'(p)+△ωi'j'(p) (21)ω i'j' (p+1)=ω i'j' (p)+△ω i'j' (p) (21)

其中，p为迭代次数，△ωi'j'(p)＝η×yi'(p)×λj'(p)为p+1次迭代中连接权值的调整部分，η为学习率，yi'(p)为神经元i'的输出信号，

Figure BDA0002920814570000151
为神经元j'的误差梯度，Xj'(p)为神经元j'的加权输入，δj'(p)为神经元j'的误差；Among them, p is the number of iterations, △ω i'j' (p)=η×y i' (p)×λ j' (p) is the adjustment part of the connection weight in p+1 iterations, and η is the learning rate , y i' (p) is the output signal of neuron i',
Figure BDA0002920814570000151
is the error gradient of neuron j', X j' (p) is the weighted input of neuron j', δ j' (p) is the error of neuron j';

S4.6、判断迭代次数是否达到最大迭代次数或者Z个输入信号的总误差E'小于误差阈值，若是，输出BP神经网络模型，否则，返回步骤S4.2进行下一次的迭代。S4.6, determine whether the number of iterations reaches the maximum number of iterations or the total error E' of the Z input signals is less than the error threshold, if so, output the BP neural network model, otherwise, return to step S4.2 for the next iteration.

步骤五：采集机械结构中观测点的位移，并将观测点的位移输入BP神经网络模型中，得到机械结构的实时疲劳寿命信息。Step 5: Collect the displacement of the observation point in the mechanical structure, and input the displacement of the observation point into the BP neural network model to obtain the real-time fatigue life information of the mechanical structure.

本发明是将对偶互易边界元法与BP神经网络相结合。首先，利用对偶互易边界元法模拟一系列机械结构的疲劳裂纹扩展寿命分析。在每次模拟中，利用从机械结构中采取的数据，采用优先数系生成一系列参数数据集。所获得的数据集包含固定观测点的结构响应和疲劳寿命的信息。此外，每次模拟获得的数据集被保存在一个单元中。然后，利用训练单元对BP神经网络模型进行训练。利用测试单元对训练好的BP神经网络模型进行验证。最后，可以利用训练好的BP神经网络模型在未知模型参数时实现结构疲劳寿命的实时预测。The invention combines the dual reciprocal boundary element method with the BP neural network. First, the dual reciprocal boundary element method is used to simulate the fatigue crack growth life analysis of a series of mechanical structures. In each simulation, a series of parametric data sets are generated using a preferential number system using data taken from the mechanical structure. The obtained dataset contains information on structural response and fatigue life at fixed observation points. In addition, the dataset obtained for each simulation is kept in a unit. Then, use the training unit to train the BP neural network model. Use the test unit to verify the trained BP neural network model. Finally, the trained BP neural network model can be used to realize the real-time prediction of structural fatigue life when the model parameters are unknown.

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。The above descriptions are only preferred embodiments of the present invention, and are not intended to limit the present invention. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present invention shall be included in the scope of the present invention. within the scope of protection.

Claims (3)

1. A real-time fatigue life prediction method of a mechanical structure based on data driving is characterized by comprising the following steps:
the method comprises the following steps: obtaining the value ranges of an index m and a coefficient c of a corresponding Paris model in the mechanical structure according to a test means, and generating a series of indexes and coefficients by adopting a priority coefficient in the value ranges of the index m and the coefficient c respectively to obtain Q groups of Paris model parameter pairs (m is a parameter pairq，cq) Wherein Q is 1,2, …, Q;
step two: setting observation points on the mechanical structure, aiming at the q group of Paris model parameter pairs (m)q，cq) Circulating, performing fatigue crack propagation analysis on the mechanical structure by using a dual reciprocity boundary element method to obtain the displacement of the observation point and the real-time fatigue life information of the mechanical structure, and positioning the observation pointMoving and combining the real-time fatigue life information of the mechanical structure to form a data information pair;
the method for carrying out fatigue crack propagation analysis on the mechanical structure by using the dual reciprocity boundary element method to obtain the displacement of an observation point and the real-time fatigue life information comprises the following steps:
s2.1, constructing dual reciprocal boundary integral equations, including a displacement boundary integral equation and a surface force boundary integral equation;
the displacement boundary integral equation comprises a displacement boundary integral equation corresponding to a non-crack boundary and a displacement boundary integral equation corresponding to a crack upper boundary; the surface force boundary integral equation refers to a surface force boundary integral equation corresponding to the lower boundary of the crack;
s2.1.1, constructing a displacement boundary integral equation corresponding to the source point falling on the non-crack boundary:
Figure FDA0003214002660000011
wherein,
Figure FDA0003214002660000012
representing a source point on a non-crack boundary, x representing a field point on a non-crack boundary, x+Field points, x, representing the upper boundary of the crack-The field point representing the lower boundary of the crack,
Figure FDA0003214002660000013
the non-crack boundary shape factor is represented,
Figure FDA0003214002660000014
Figure FDA0003214002660000015
a fundamental solution of displacement representing a field point falling on a non-crack boundary,
Figure FDA0003214002660000016
surface force fundamental solution representing field point falling on non-crack boundary，
Figure FDA0003214002660000017
A basic solution of the displacement representing the boundary where the field point falls on the crack,
Figure FDA0003214002660000018
a basic solution of the displacement representing the field point falling on the lower boundary of the crack,
Figure FDA0003214002660000019
the base solution of the face force representing the field point falling on the upper boundary of the crack,
Figure FDA00032140026600000110
the base solution of the face force representing the field point falling at the lower boundary of the crack,
Figure FDA00032140026600000111
Figure FDA00032140026600000112
s represents the non-crack boundary and,
Figure FDA00032140026600000113
the upper boundary of the crack is represented,
Figure FDA00032140026600000114
the lower boundary of the crack is represented,
Figure FDA00032140026600000115
representing the j-direction displacement component, u, at the origin point on the non-crack boundaryj(x) Representing the j-direction displacement component, t, of field points at non-crack boundariesj(x) Representing the force component of the j-direction surface at the field point on the non-crack boundary, uj(x+) A j-direction displacement component, u, at a field point representing an upper boundary of the crackj(x-) A j-direction displacement component, t, at a field point representing the lower boundary of the crackj(x+) Indicating the surface force in the j-direction at the upper boundary of the crackComponent, tj(x-) Representing the j-direction surface force component at the lower boundary of the crack,
Figure FDA00032140026600000116
the basic solution of the displacement is represented,
Figure FDA00032140026600000117
the basic solution of the surface force is shown,
Figure FDA0003214002660000021
the field point y' is x, x+Or x-，
Figure FDA0003214002660000022
Figure FDA0003214002660000023
Representing a delta function, v representing the Poisson's ratio, r representing the distance of the source point and the field point, n representing a normal vector, niComponent representing the direction of normal vector i, njRepresenting the component in the direction of the normal vector j,
Figure FDA0003214002660000024
Figure FDA0003214002660000025
representing a source point
Figure FDA0003214002660000026
Of i-directional component, y'iThe i-direction component of the field point y' is represented,
Figure FDA0003214002660000027
Figure FDA0003214002660000028
representing a source point
Figure FDA0003214002660000029
Of j-directional component, y'jRepresents the j-direction component of field point y';
the formula (1) is sorted and combined to obtain:
Figure FDA00032140026600000210
s2.1.2, constructing a displacement boundary integral equation corresponding to the boundary of the source point falling on the crack:
Figure FDA00032140026600000211
wherein,
Figure FDA00032140026600000212
indicating the origin of the boundary that falls on the crack,
Figure FDA00032140026600000213
indicates the source point falling on the lower boundary of the crack, and
Figure FDA00032140026600000214
the coordinates representing the two source points are the same,
Figure FDA00032140026600000215
the upper boundary shape factor of the crack is represented,
Figure FDA00032140026600000216
the crack lower boundary shape factor is expressed,
Figure FDA00032140026600000217
a j-direction displacement component at the origin of the boundary on the crack,
Figure FDA00032140026600000218
representing a j-direction displacement component at a source point representing a boundary on a crack，
Figure FDA00032140026600000219
A fundamental solution of displacement representing a field point falling on a non-crack boundary,
Figure FDA00032140026600000220
representing the base solution of the face force with the field point falling on the non-crack boundary,
Figure FDA00032140026600000221
a basic solution of the displacement representing the field point falling on the lower boundary of the crack,
Figure FDA00032140026600000222
a base solution of the face force representing the field point falling on the upper boundary of the crack;
s2.1.3, constructing a surface force boundary integral equation corresponding to the boundary of the source point falling on the lower boundary of the crack:
Figure FDA00032140026600000223
wherein,
Figure FDA00032140026600000224
the j-direction surface force component at the origin of the crack lower boundary,
Figure FDA00032140026600000225
the j-direction surface force component at the origin of the boundary on the crack,
Figure FDA00032140026600000226
i-directional component, t, at the origin of the lower boundary of the crackk(x) Representing the k-direction surface force component at the field point on the non-crack boundary,
Figure FDA00032140026600000227
a high order fundamental solution of displacement representing a field point falling on a non-crack boundary,
Figure FDA00032140026600000228
high order surface force fundamental solution, u, representing a field point falling on a non-crack boundaryk(x) Representing the k-direction displacement component, t, of the field points at the non-crack boundariesk(x+) Representing the surface force component in the k-direction, t, at the boundary field point above the crackk(x-) Represents the surface force component in the k direction at the lower boundary of the crack, uk(x+) Representing the displacement component in the k-direction, u, at the boundary field point above the crackk(x-) Representing the k-direction displacement component at the boundary field point below the crack,
Figure FDA0003214002660000031
g' represents a shear modulus and a shear modulus,
Figure FDA0003214002660000032
Figure FDA0003214002660000033
Figure FDA0003214002660000034
representing a source point
Figure FDA0003214002660000035
K-direction component of (2), y'kRepresenting the k-direction component, n, of the field point ykThe component in the direction of the normal vector k,
Figure FDA0003214002660000036
Figure FDA0003214002660000037
s2.2, respectively carrying out discrete processing on the displacement boundary integral equation and the surface force boundary integral equation to obtain:
Figure FDA0003214002660000038
wherein H*Integral representing basic solution of displacement, G*Represents the integral of the surface force base solution, S represents the non-crack boundary,
Figure FDA0003214002660000039
the upper boundary of the crack is represented,
Figure FDA00032140026600000310
denotes the lower boundary of the crack, uSIndicating the displacement on the non-crack boundary,
Figure FDA00032140026600000311
indicating the displacement on the upper boundary of the crack,
Figure FDA00032140026600000312
indicates the displacement at the lower boundary of the crack, tSIndicating the face force at the non-crack boundary,
Figure FDA00032140026600000313
indicating the face force at the upper boundary of the crack,
Figure FDA00032140026600000314
representing the face force at the lower boundary of the crack;
s2.3, applying a boundary condition to the formula (5), and finishing to obtain:
Figure FDA00032140026600000315
wherein A is*Representing a known amount of displacement or surface force, xSRepresenting an unknown displacement or face force at a non-crack boundary,
Figure FDA00032140026600000316
representing an unknown displacement or surface force on the upper boundary of the crack,
Figure FDA00032140026600000317
representing an unknown displacement or surface force at the lower boundary of the crack, bSRepresenting the combined components of known displacements or face forces at non-crack boundaries,
Figure FDA00032140026600000318
representing the combined resulting components of known displacements or face forces at the upper boundary of the crack,
Figure FDA00032140026600000319
representing the combined components of known displacements or surface forces at the lower boundary of the crack;
solving the formula (6) by a Gaussian elimination method to obtain the displacement and the surface force of the boundary, and substituting the displacement and the surface force of the boundary into the formula (1) to obtain the displacement of the observation point;
s2.4, calculating a mixed modal stress intensity factor by utilizing an interaction integral algorithm according to the displacement of the boundary:
Figure FDA0003214002660000041
wherein M is(1,2)The integral of the interaction is represented by,
Figure FDA0003214002660000042
in order to interact with the strain energy density,
Figure FDA0003214002660000043
Figure FDA0003214002660000044
representing the i-direction displacement component in state two,
Figure FDA0003214002660000045
denotes the i-direction displacement component, x, at state one1Representing the component of the inner point x, and Γ representing the interactionPath of integration, AεIs a circle with the center of a crack tip and the radius of R,
Figure FDA0003214002660000046
which represents the strain in the first state of the state,
Figure FDA0003214002660000047
which represents the strain in the second state,
Figure FDA0003214002660000048
the stress in the first state is shown,
Figure FDA0003214002660000049
representing the stress in state two; the first state represents the real state of the mechanical structure, and the second state represents the auxiliary state;
s2.5, converting the formula (7) into:
Figure FDA00032140026600000410
wherein E is a Young's model,
Figure FDA00032140026600000411
representing the type I stress intensity factor in state one,
Figure FDA00032140026600000412
represents the type I stress intensity factor in the second state,
Figure FDA00032140026600000413
representing the type II stress intensity factor in state one,
Figure FDA00032140026600000414
representing the type II stress intensity factor in state two,
Figure FDA00032140026600000415
representing the type III stress intensity factor in state one,
Figure FDA00032140026600000416
representing the type III stress intensity factor in the second state;
s2.6, calculating a mixed modal stress intensity factor delta Keq：
Figure FDA00032140026600000417
Wherein, KIDenotes the type I stress intensity factor, KIIRepresenting a type II stress intensity factor;
s2.7, determining the crack direction by adopting a maximum circumferential stress criterion:
Figure FDA00032140026600000418
wherein theta is a crack propagation angle, and sign represents a sign function;
s2.8, stress intensity factor delta K according to mixed modeeqDetermining the crack propagation rate, wherein for the Paris model, the crack propagation rate is as follows:
Figure FDA00032140026600000419
wherein a is the crack length and N is the stress cycle number;
s2.9, depending on the initial length a of the crack0The critical length L of the crack, the crack growth increment Δ a in the crack growth step, and the total number of growth steps is known to be
Figure FDA00032140026600000420
During crack propagation, for Δ K in the t stepeqCalculating the corresponding load cycle number DeltaNt：
Figure FDA0003214002660000051
Wherein T is more than or equal to 1 and less than or equal to T;
s2.10, load cycle times delta NTAs real-time fatigue life information for the mechanical structure;
step three: circularly executing the step two until the Q groups of Paris model parameter pairs are traversed to obtain a data set, and dividing the data set into a training set and a testing set, wherein the data set comprises a series of data information pairs;
step four: inputting the training set into a BP neural network for training to obtain a BP neural network model, and verifying the BP neural network model by using a test set;
step five: and acquiring the displacement of an observation point in the mechanical structure, and inputting the displacement of the observation point into the BP neural network model to obtain the real-time fatigue life information of the mechanical structure.
2. The method for predicting the fatigue life of a mechanical structure in real time based on data driving of claim 1, wherein during the crack propagation process, the current crack length is a0Setting the critical length as L and the expansion amount of each step as delta a, totaling the expansion
Figure FDA0003214002660000052
Step, pair formula
Figure FDA0003214002660000053
Performing two-end integration after deformation to obtain the load cycle number delta N of the first step of expansion1：
Figure FDA0003214002660000054
Performing the second expansion to obtain the load cycle times delta N of the second expansion2：
Figure FDA0003214002660000055
Carrying out third step expansion to obtain load cycle times delta N of the third step expansion3：
Figure FDA0003214002660000056
The same fourth step is extended until the T-1 step is extended, and the corresponding load cycle number is delta N4...ΔNT-1；
Expanding the T step to obtain the load cycle times delta N of the T step expansionT：
Figure FDA0003214002660000057
Therefore, the total life of the mechanical structure is Δ N ═ Δ N1+ΔN2+...+ΔNT-1+ΔNT。
3. The method for predicting the real-time fatigue life of the mechanical structure based on the data driving as claimed in claim 1, wherein the method for inputting the training set into the BP neural network for training to obtain the BP neural network model comprises the following steps:
s4.1, network initialization: determining the node numbers of the input layer, the hidden layer and the output layer, and initializing the connection weight Wi'j'And Wj'k'And thresholds for the hidden layer and the output layer, i ' ═ 1,2, …, n ', the number of neurons in the input layer, j ' ═ 1,2, …, l, the number of neurons in the hidden layer, k ' ═ 1,2, …, m ', the number of neurons in the output layer;
s4.2, calculating the output of the hidden layer node:
Figure FDA0003214002660000061
wherein h isj'Representing the output, x, of neuron ji'Representing the input component, ωi'j'∈Wi'j'Represents a weighted weight value, θj'Is the threshold value for the neuron j',
Figure FDA0003214002660000062
is an activation function;
s4.3, calculating the output of the output layer node:
Figure FDA0003214002660000063
wherein, yk'Represents the output of the neuron k', ωj'k'∈Wj'k'Represents a weighted weight value, θk'A threshold for neuron k';
s4.4, calculating the error e between the expected output and the predicted output of the z-th training dataz：
Figure FDA0003214002660000064
Wherein, yd,k'For the expected output value, y, of the neuron k' at the output layerk'For the predicted output value of neuron k 'at the output layer, the total error of the Z input signals is E':
Figure FDA0003214002660000065
s4.5, updating the connection weight of the neural network:
ωi'j'(p+1)＝ωi'j'(p)+Δωi'j'(p) (21)；
where p is the number of iterations, Δ ωi'j'(p)＝η×yi'(p)×λj'(p) is the adjustment part of the connection weight in p +1 iterations, η is the learning rate, yi'(p) is the output signal of neuron i',
Figure FDA0003214002660000066
error gradient, X, for neuron jj'(p) is the weighted input to neuron j', δj'(p) is the error for neuron j';
and S4.6, judging whether the iteration number reaches the maximum iteration number or whether the total error E' of the Z input signals is smaller than an error threshold, if so, outputting the BP neural network model, and if not, returning to the step S4.2 to perform the next iteration.
CN202110116397.0A 2021-01-28 2021-01-28 Mechanical structure real-time fatigue life prediction method based on data driving Active CN112784495B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110116397.0A CN112784495B (en) 2021-01-28 2021-01-28 Mechanical structure real-time fatigue life prediction method based on data driving

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110116397.0A CN112784495B (en) 2021-01-28 2021-01-28 Mechanical structure real-time fatigue life prediction method based on data driving

Publications (2)

Publication Number Publication Date
CN112784495A CN112784495A (en) 2021-05-11
CN112784495B true CN112784495B (en) 2021-09-24

Family

ID=75759305

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110116397.0A Active CN112784495B (en) 2021-01-28 2021-01-28 Mechanical structure real-time fatigue life prediction method based on data driving

Country Status (1)

Country Link
CN (1) CN112784495B (en)

Families Citing this family (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115169694B (en) * 2022-07-06 2023-05-30 天津大学 Method for realizing subcritical crack dynamic expansion and life prediction
CN115169240B (en) * 2022-07-26 2024-11-19 南京理工大学 Machine learning-based TiAl alloy fatigue crack growth life prediction method
CN117649901B (en) * 2023-11-29 2024-07-09 哈尔滨工业大学 Interaction integration method for solving stress intensity factors of cracks of revolving body
CN117669390A (en) * 2024-02-01 2024-03-08 中国石油大学(华东) Metal full-stage fatigue crack growth prediction method and system based on neural network
CN118095017B (en) * 2024-04-24 2024-07-30 中国海洋大学 A fatigue life prediction method

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104062196A (en) * 2014-01-08 2014-09-24 中国石油大学(华东) Corrosion fatigue life prediction method based on damage evolution
CN105488328A (en) * 2015-11-19 2016-04-13 北京航空航天大学 Fatigue crack growth rate prediction method based on artificial neuron network
CN105956315A (en) * 2016-05-17 2016-09-21 北京航空航天大学 Method capable of carrying out fatigue crack propagation rate estimation and life prediction
CN108763841A (en) * 2018-07-24 2018-11-06 北京航空航天大学青岛研究院 A kind of elastic failure emulation mode based on Dual boundary element and strain energy optimization analysis
CN112052615A (en) * 2020-09-07 2020-12-08 郑州航空工业管理学院 Micro-motion fatigue performance prediction method based on artificial neural network
CN112257197A (en) * 2020-10-19 2021-01-22 郑州轻工业大学 A method for evaluating working stress of micro-defects in large castings and forgings

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104062196A (en) * 2014-01-08 2014-09-24 中国石油大学(华东) Corrosion fatigue life prediction method based on damage evolution
CN105488328A (en) * 2015-11-19 2016-04-13 北京航空航天大学 Fatigue crack growth rate prediction method based on artificial neuron network
CN105956315A (en) * 2016-05-17 2016-09-21 北京航空航天大学 Method capable of carrying out fatigue crack propagation rate estimation and life prediction
CN108763841A (en) * 2018-07-24 2018-11-06 北京航空航天大学青岛研究院 A kind of elastic failure emulation mode based on Dual boundary element and strain energy optimization analysis
CN112052615A (en) * 2020-09-07 2020-12-08 郑州航空工业管理学院 Micro-motion fatigue performance prediction method based on artificial neural network
CN112257197A (en) * 2020-10-19 2021-01-22 郑州轻工业大学 A method for evaluating working stress of micro-defects in large castings and forgings

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
Singularity cancellation method for time-domain boundary element formulation of elastodynamics: A direct approach;Guizhong Xie et al.;《Applied Mathematical Modelling》;20191206;全文 *

Also Published As

Publication number Publication date
CN112784495A (en) 2021-05-11

Similar Documents

Publication Publication Date Title
CN112784495B (en) Mechanical structure real-time fatigue life prediction method based on data driving
Teng et al. Structural damage detection based on transfer learning strategy using digital twins of bridges
Chang et al. Learning to simulate and design for structural engineering
Novák et al. ANN inverse analysis based on stochastic small-sample training set simulation
CN114117840B (en) A hybrid-driven structural performance prediction method based on simulation and experimental data
Fei et al. Vectorial surrogate modeling method for multi-objective reliability design
Feng et al. Data-driven algorithm for real-time fatigue life prediction of structures with stochastic parameters
Wu et al. Prediction of stress intensity factors in pavement cracking with neural networks based on semi-analytical FEA
Liu FEA-AI and AI-AI: Two-way deepnets for real-time computations for both forward and inverse mechanics problems
Fathnejat et al. An efficient two-stage approach for structural damage detection using meta-heuristic algorithms and group method of data handling surrogate model
Wang et al. Updating multiscale model of a long-span cable-stayed bridge
Maity et al. Damage assessment in structure from changes in static parameter using neural networks
Bui-Tien et al. Damage detection in structural health monitoring using hybrid convolution neural network and recurrent neural network
CN104200005A (en) Bridge damage identification method based on neural network
CN112347670B (en) A method for predicting creep parameters of rockfill material based on neural network response surface
Aydin et al. Damage diagnosis in beam-like structures by artificial neural networks
Tohidi et al. A new predictive model for restrained distortional buckling strength of half-through bridge girders using artificial neural network
Hasançebi et al. Detailed load rating analyses of bridge populations using nonlinear finite element models and artificial neural networks
Aydin et al. Damage detection in Timoshenko beam structures by multilayer perceptron and radial basis function networks
CN117594164A (en) Metal structure residual fatigue life calculation and evaluation method and system based on digital twin
CN111783212B (en) A typical damage identification method for cable-stayed bridges
Pishro et al. Advancing ultimate bond stress–slip model of UHPC structures through a novel hybrid machine learning approach
CN114048646B (en) Asphalt surface layer damage state inversion method based on finite element correction and artificial intelligence
Bilgehan et al. Buckling load estimation of cracked columns using artificial neural network modeling technique
Li et al. Model Updating Method for Jacket Platform Considering Different Component Degradation Based on Deep Learning

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant