Steel rail hidden damage crack initiation simulation prediction method considering three-dimensional geometric irregularity of rail surface

Abstract

The invention discloses a steel rail hidden damage crack initiation simulation prediction method considering three-dimensional geometric irregularity of a rail surface. The method comprises the following steps: step 1, acquiring wheel-rail dynamic response and strain rate space-time distribution characteristics; 2, calculating the local elastic-plastic contact stress of the wheel rail; 3, establishing a circulating rolling three-dimensional local finite element model, and solving a steel rail plastic stress-strain solution; and 4, solving the steel rail hidden damage crack initiation angle and the crack initiation life by adopting a plastic stress-strain solution. According to the steel rail hidden damage crack initiation simulation prediction method considering the three-dimensional geometric irregularity of the rail surface, the three-dimensional geometric irregularity of the surface of the steel rail is considered, the application range of traditional dynamics is expanded, and initiation and development of steel rail hidden damage cracks can be predicted more accurately; and a guidance basis is provided for formulating a reasonable grinding period and grinding depth on site and selecting a material of a newly-replaced steel rail, and the method has a wide application prospect.

Description

Technical Field

[0001] The invention belongs to the technical field of wheel-rail relationship and material performance in railway engineering, and relates to a method for determining service performance requirements of rails. Background Art

[0002] There are two main methods for predicting wheel-rail rolling contact fatigue damage: one is engineering application-oriented, starting from the perspective of the vehicle-rail system, using the vehicle-rail coupling dynamics model to simulate the dynamic interaction between the wheel and rail (such as normal force, tangential force, creep rate, contact spot size, etc.), and then using the damage prediction method oriented to engineering applications to predict wheel-rail rolling contact fatigue damage. Therefore, in theory, the details of wheel-rail rolling contact are omitted for ease of application. The second is mechanism-oriented, starting from the perspective of wheel-rail dynamic interaction and continuous medium elastic-plastic mechanics, combined with the micromechanical behavior of the material, to reveal the initiation mechanism of fatigue damage. The wheel-rail rolling contact pair is subjected to large cyclic high loads, and accurate and detailed cyclic stress, plastic deformation analysis and appropriate fatigue damage prediction model are the key to clarify the rolling contact fatigue mechanism. However, the existing methods still have shortcomings in predicting wheel-rail rolling contact fatigue damage: on the one hand, the three-dimensional geometric unevenness of the rail surface is not considered in the wheel-rail refined modeling, and the stress-strain field in the contact body cannot be accurately calculated; on the other hand, the numerical algorithm and material constitutive model used in the current elastic-plastic contact model have not yet considered the dynamic stress-strain relationship related to the strain rate. Therefore, there is currently a lack of a simulation prediction method for the initiation of hidden cracks in rails that can simultaneously consider the three-dimensional geometric irregularities of the rail surface and high strain rates.

[0003] At present, the rapid prediction methods for wheel-rail rolling contact fatigue include:

[0004] 1) Based on the stability diagram, the concept of saturation is introduced and expanded to obtain the extended stability diagram;

[0005] 2) Using the normal stress and tangential stress distribution on the rail surface and the potential function method, the stress-strain field distribution inside the wheel-rail body under steady-state rolling conditions is solved;

[0006] 3) Combined with the extended stability diagram in step 1, select an appropriate fatigue model, and use the stress-strain amplitude obtained in step 2 to predict fatigue life based on the critical plane method.

[0007] On the one hand, the rapid prediction method for wheel-rail rolling contact fatigue does not take into account the three-dimensional geometric unevenness of the rail surface, and errors will occur in the calculation results when calculating the elastic-plastic stress-strain field, resulting in inaccurate elastic-plastic contact loads, which will greatly affect the life prediction; on the other hand, the material constitutive model used in the prediction of this method does not take the strain rate into account, but the strain rate has a great influence on the critical damage value, thereby affecting the life prediction. Summary of the invention

[0008] In order to overcome the above-mentioned defects of the prior art, the present invention proposes a simulation prediction method for the initiation of hidden cracks in rails taking into account the three-dimensional geometric unevenness of the rail surface.

[0009] A simulation prediction method for the initiation of hidden cracks in rails considering the three-dimensional geometric irregularity of the rail surface, characterized in that the method comprises

[0010] Step 1: Obtain the wheel-rail dynamic response and strain rate temporal and spatial distribution characteristics;

[0011] Step 2: Calculate the local elastic-plastic contact stress between wheel and rail;

[0012] Step 3: Establish a cyclic rolling three-dimensional local finite element model to solve the rail plastic stress-strain solution;

[0013] Step 4 uses the plastic stress-strain solution to solve the rail hidden crack initiation angle and crack initiation life.

[0014] Wherein, the wheel-rail dynamic response and strain rate spatiotemporal distribution characteristics obtained in step 1 include:

[0015] (1.1) Numerical simulation of three-dimensional geometric irregularities of track surface measured by actual measurements;

[0016] (1.2) Establish a vehicle-track coupling dynamics model to solve the wheel-rail dynamic response under rolling contact;

[0017] (1.3) Establish a subsurface strain calculation model and solve the temporal and spatial distribution characteristics of the rail strain rate under rolling contact.

[0018] Among them, in the numerical simulation of the three-dimensional geometric irregularity of the track surface measured in (1.1),

[0019] (1.1.1) In order to obtain the three-dimensional unevenness distribution data of rail hidden damage, the rail surface is scanned by a composite three-dimensional scanner (such as HandySCAN3D), and the data is processed by post-processing software. The specific steps are as follows:

[0020] ① Post marking points near the hidden damage of the rail to determine the spatial position of the rail, and calibrate the instrument accuracy based on the marking points;

[0021] ②Control the 3D scanner to move continuously to obtain the 3D information of the position where the laser passes, thereby forming continuous 3D track surface point cloud data;

[0022] ③ Import the 3D point cloud data into the post-processing software, select the reference object to determine the overall coordinate system of the point cloud data, and finally obtain the correct coordinates of the rail surface point cloud data;

[0023] ④ Interpolate the randomly distributed point cloud data and integrate the point cloud data into a matrix expression to ensure that the coordinates along the horizontal and vertical directions remain consistent, so as to introduce the track surface data considering three-dimensional unevenness into the subsequent dynamic model;

[0024] ⑤ Through the filter, the point cloud data is processed with high frequency along the longitudinal direction of the rail surface to eliminate the problem of fluctuation of point cloud data caused by the accuracy problem of the 3D scanner.

[0025] (1.1.2) The real three-dimensional geometric point cloud data containing hidden damage on the rail is discretized into grid units along the longitudinal direction, and the wheel is discretized into rotating body grid units along the longitudinal direction. The wheel-rail spatial geometry solution method is used to obtain the wheel-rail contact points and the three-dimensional asymmetric contact geometry gap distribution.

[0026] Among them, the vehicle-track coupling dynamic model is established in (1.2) to solve the wheel-rail dynamic response under rolling contact.

[0027] (1.2.1) The vehicle model considers the car body, frame, wheelset and suspension system. The car body and bogie are considered as rigid bodies, and the suspension system is considered as spring and damping units. The track model includes rails, fastener systems and ballast. The rails are simulated using the Euler beam model, and the fastener system is considered as a spring-damper unit.

[0028] (1.2.2) The wheel-rail contact geometric clearance obtained in step (1.1) is brought into the dynamic model. In this way, the wheel-rail contact model is used to solve the accurate wheel-rail rolling contact behavior under dynamic response, and the dynamic parameters required for each integral step are obtained, including the wheel-rail contact force and creep rate.

[0029] Among them, the above (1.3) establishes a subsurface strain calculation model and solves the temporal and spatial distribution characteristics of the rail strain rate under rolling contact.

[0030] (1.3.1) While calculating the wheel-rail dynamic response in (1.2), the strain field results at each integration step are calculated in real time. This is the innovation of this patent compared with the traditional dynamic model. First, the normal contact stress matrix p is obtained at each integration step. N and the tangential contact stress matrix s N , and the influence factor matrix D under normal stress is obtained by using the influence factor method N And the influence factor matrix T under the action of tangential stress N , the subscript N represents the number of grids divided in the contact area based on the boundary element method, and then the above matrix is ​​operated as follows for numerical calculation:

[0031] First, the contact stress matrix and the influence factor matrix are topologically converted into calculation matrices with twice the number of grids in the transverse and longitudinal directions, where the stress matrix p for the topology is2N and 2N The elements are filled with 0 values, and the impact factor matrix D 2N and T 2N Elements are filled using circular convolution;

[0032] Then, the matrix is ​​transformed from the time domain to the frequency domain using discrete Fourier transform.

[0033] (1.3.2) Calculate the product of the topological normal contact stress matrix and the normal influence factor matrix in the frequency domain, and the product of the tangential contact stress matrix and the tangential influence factor matrix in the frequency domain, and then use the inverse Fourier transform to convert the calculation results in the frequency domain to the time domain:

[0034]

[0035] Where IFFT stands for inverse Fourier transform, ρ is the calculation result of topology in half space to twice the number of grids, q and r correspond to the x, y, and z directions respectively; then the matrix elements of the last N+1 items of the matrix topology are discarded to obtain the stress value σ corresponding to each rail target point qr , which contains 9 stress components, and then the strain component results can be calculated by generalized Hooke's law;

[0036] (1.3.3) The strain rate is calculated by the strain field results at each integration step. The specific steps are as follows:

[0037] ① Select the rail surface strain rate calculation point as the target point O, so that the wheel moves at a unit distance (unit distance = integration step Δt p × running speed v) moving;

[0038] ② Determine the calculation area along the longitudinal direction with the target point O as the center and divide the grid. When the wheel passes through the left boundary x of the calculation area L When the wheel passes through the right boundary x of the calculation area, the wheel-rail contact stress matrix is ​​interpolated to the target point coordinates to calculate the strain field. R When , it means the wheel has left and the strain rate is no longer calculated;

[0039] ③ Divide the time history of the wheel passing through the calculation area into T discrete time points

[0040] t 1,2,…,T ∈(x L / v,x R / v), calculate the strain rate results at each point;

[0041] ④ The strain rate reference value of the target point O is the maximum value among a series of point strain rates, so the strain rate at the target point O is solved by the following formula:

[0042]

[0043] In the formula, Indicates t i The strain result at the moment. i+1 -t i It can be expressed as the strain rate calculation step size set in the dynamics calculation.

[0044] Wherein, in the step 2: calculating the local elastic-plastic contact stress between wheel and rail,

[0045] The dynamic parameters and strain rate results calculated by the dynamic model are used as

[0046] Input quantity, the wheel-rail elastic-plastic rolling contact equation and sliding equation are proposed:

[0047]

[0048] Where P n and P τ represents normal and tangential stress; A is defined as the load-displacement influence coefficient; h is the normal gap when the two contact bodies are in rigid contact; ω is the sliding amount; g I is the tangential contact boundary; and It manifests itself as plastic deformation in the normal and tangential directions during contact; c x 、c y , are the longitudinal, transverse and spin creep rates,

[0049] The contact equations (c) and (d) are extended to consider the influence of plastic deformation of the contact interface on the contact solution; the tangential sliding equation (e) is extended to consider the deformation rate caused by stress fluctuations in dynamic contact, which was ignored in previous calculations due to the steady-state assumption. These two items are extended to reflect the mutual coupling influence of contact stress and elastic-plastic deformation during wheel-rail impact and stress

[0050] Effect of fluctuations on the strain field.

[0051] Solve the above equations and establish the wheel-rail elastic-plastic rolling contact model to solve the elastic-plastic contact stress. The specific steps are as follows:

[0052] (2.1) Establishing elastic contact mechanics model to calculate local elastic solution

[0053] Using a graded loading method, the wheel-rail macro-dynamic parameters are introduced into the elastic contact model to calculate the local contact stress under the current load conditions, and the sub-surface stress-strain calculation model established in (1.1.3) is used to calculate the internal stress field results of the rail;

[0054] (2.2) Based on the Von Mises yield criterion, the plasticity occurrence area is determined and the plastic strain is solved:

[0055] The yield judgment formula (f) is established. When the yield function f>0, it means that the material in this area exceeds the yield limit and plastic deformation occurs; otherwise, it is still in elastic deformation.

[0056]

[0057] Where δ qr is the Kronecker symbol, σ y is the yield stress, and its value increases with the strain rate.

[0058] The strain rate calculated by (1.1.3) is used as the input parameter, and the material constitutive model related to the strain rate at different times is established in combination with the Johnson-Cook method. The radial return algorithm proposed by Simon and Hughes is used to solve the plastic strain.

[0059] (2.3) Using the eigenstrain method to solve the internal residual stress and surface plastic deformation

[0060] The internal area of ​​the rail is discretely divided into rectangular grid units of equal size, and the plastic characteristics in each unit are assumed to be constant. The influence factors characterizing the relationship between plastic strain and residual stress given by Liu et al. are used to solve the residual stress of each unit. The residual stress generated by each unit is superimposed to obtain the residual stress distribution in the plastic area. Finally, the plastic deformation of the rail surface is solved with the help of residual stress. The idea is similar to that of residual stress, and both are calculated with the help of the influence factor method.

[0061] (2.4) Obtaining wheel-rail elastic-plastic rolling contact stress

[0062] The surface plastic deformation ( and ) is fed back to equations (c) to (e), and the wave equation (g) is introduced. This equation is used to solve the size of the wave displacement under unit displacement excitation, so the plastic deformation can be introduced into equation (g), and the wave displacement u of the potential contact area can be obtained by using the superposition principle. Repeat steps (a) to (c) until the contact stress calculated by the two iterations converges, and it is considered that the current load is balanced. Continue to increase the load until the wheel-rail contact force calculated by (1.1) is reached, and the calculation can be terminated to obtain the elastic-plastic contact stress at the current moment.

[0063]

[0064] Where λ and μ are Lame constants, r and z represent radial and vertical coordinates, and ρ is the material density;

[0065] Wherein, the step 3: establishing a cyclic rolling three-dimensional local finite element model to solve the rail plastic stress-strain solution,

[0066] A three-dimensional local finite element model is established, longitudinal constraints are imposed on both end surfaces of the rail, and full constraints are adopted on the bottom surface of the rail. The distribution of the plastic stress-strain field of the rail needs to consider the constitutive relationship of the material, and the rail material is adjusted by modifying the UMAT subroutine in combination with the actual situation.

[0067] In this model, the translation of normal stress and tangential stress in the rolling direction is used to simulate the cyclic rolling process. Considering the non-steady-state rolling contact at the three-dimensional geometric unevenness of the rail, the normal stress and tangential stress both vary with time, and the load history is taken from the elastic-plastic contact stress distribution at each moment in step 1. The cyclic moving load on the rail surface is implemented in this patent by modifying the inp file of the finite element software through a self-compiled program. At the initial position point A, the contact load gradually increases through 10 incremental steps, and then the contact load gradually moves longitudinally to the end point B. The contact load moves a distance of 1mm in one incremental step, and the total rolling distance is 100mm. After 100 steps of rolling, the contact spot reaches point B and finally decreases to zero through 10 incremental steps. In order to save calculation costs without affecting the calculation results, five load cycles can be used in the calculation, and the plastic stress-strain field distribution of the rail can be obtained through this model.

[0068] Wherein, the step 4 uses the plastic stress-strain solution to solve the rail hidden crack initiation angle and crack initiation life.

[0069] (4.1) When predicting the crack initiation life of rails, an appropriate rolling contact fatigue life prediction model is selected according to the local material response of the rails, such as the Jiang-Sehitoglu model suitable for the ratchet effect. This model considers that fatigue damage is caused by the plastic strain energy in the critical plane of the material, and its fatigue damage increment can be expressed as:

[0070]

[0071] Where <> is the Macauley bracket; dD is the fatigue damage increment; σmr, σ0, σf are the material memory parameter, fatigue strength limit and true fracture stress respectively; σ and τ are the normal compressive stress and shear stress respectively; dε p and dγ p are the plastic strain increments corresponding to σ and τ respectively; b and m are material constants; dY is the plastic strain energy increment in the critical material plane;

[0072] (4.2) Use the critical plane method to calculate the maximum value of fatigue damage parameters and predict the initiation angle and crack life of rail hidden damage.

[0073] The corresponding fatigue damage parameter value can be calculated for each plane at any point in the material, and the plane where the maximum value of the fatigue damage parameter is located is the critical plane. According to the plastic stress-strain results of the rail under each cyclic load calculated by the finite element, the damage value generated at a certain point can be obtained. When the damage of a point on the critical surface reaches the critical damage value D0, the point is failed. And considering the influence of strain rate on ductility, a structure consistent with the site is obtained, and a correction equation related to strain rate is introduced in the critical damage value, such as formula (j):

[0074]

[0075] In the formula, and are strain rate and reference strain rate; m1 is an undetermined coefficient, which is determined by the phenomenon observed on site and indoor rolling test;

[0076] According to the stress state at one point, the plastic stress-strain result in step 3 is transformed into a tensor coordinate with a step size of 1° to obtain the stress-strain state on any plane at this point, and the fatigue damage parameter value of the plane is solved according to equations (h), (i) and (j). During the coordinate transformation, the engineering shear strain output by the finite element method should be converted into a tensor shear strain. At the same time, the fatigue damage parameter values ​​of all target points are calculated, and the proportion of the two components, normal stress-strain and shear stress-strain, is quantitatively analyzed to find the maximum value of the fatigue damage parameter. The plane angle corresponding to the maximum value of the fatigue damage material is the crack initiation angle;

[0077] The crack initiation life Nf is defined as:

[0078]

[0079] Where ΔD is the cumulative damage value under a cyclic loading.

[0080] The final rail hidden damage cracking life is obtained by solving the definition formula of cracking life.

[0081] There are two main simulation methods for dynamic simulation of rail surface fatigue damage: multi-body dynamics and transient finite element. The multi-body dynamics method is widely used due to its better computational efficiency, while the traditional dynamic method is still based on multiple assumptions. For example, it assumes that the lateral contact position of the wheel and track profile remains unchanged at the unevenness, only considers the two-dimensional geometric unevenness in the longitudinal and vertical directions, and ignores the influence of the lateral unevenness width on the wheel-rail contact position conversion and the wheel-rail impact force.

[0082] The invention proposes a simulation prediction method for the initiation of hidden cracks in rails that takes into account the three-dimensional geometric unevenness of the rail surface. This method takes into account the three-dimensional geometric unevenness of the rail surface, expands the application scope of traditional dynamics, and can more accurately predict the initiation and development of hidden cracks in rails. It provides guidance for formulating reasonable grinding cycles, grinding depths, and material selection for new rails on site, and has broad application prospects. BRIEF DESCRIPTION OF THE DRAWINGS

[0083] Figure 1 Schematic diagram of the local finite element model of rail.

[0084] Figure 2a This is a schematic diagram of the point cloud data of the measured hidden damage on the rail surface.

[0085] Figure 2b Schematic diagram of profile deviation calculation at different positions.

[0086] Figure 2c Schematic diagram of the track surface hidden damage matrix point cloud data formed by interpolation.

[0087] Figure 3 Schematic diagram of the change curve of wheel-rail macro-dynamic parameters with hidden damage and unevenness.

[0088] Figure 4 Schematic diagram of the strain rate variation curve of the hidden injury unevenness.

[0089] Figure 5 Schematic diagram of the wheel-rail elastic-plastic contact stress distribution at the hidden damage and unevenness.

[0090] Figure 6 Schematic diagram of rail plastic stress-strain calculation results. DETAILED DESCRIPTION

[0091] A method for simulating and predicting the initiation of hidden cracks in rails taking into account the three-dimensional geometric irregularities of the rail surface provided by the present invention is described in detail below in conjunction with the accompanying drawings and specific embodiments.

[0092] Figure 1 As shown, a simulation prediction method for the initiation of hidden cracks in rails considering the three-dimensional geometric unevenness of the rail surface is characterized in that the method includes

[0093] Step 1: Obtain the wheel-rail dynamic response and strain rate temporal and spatial distribution characteristics;

[0094] Step 2: Calculate the local elastic-plastic contact stress between wheel and rail;

[0095] Step 3: Establish a cyclic rolling three-dimensional local finite element model to solve the rail plastic stress-strain solution;

[0096] Step 4 uses the plastic stress-strain solution to solve the rail hidden crack initiation angle and crack initiation life.

[0097] Wherein, the wheel-rail dynamic response and strain rate spatiotemporal distribution characteristics obtained in step 1 include:

[0098] (1.1) Numerical simulation of three-dimensional geometric irregularities of track surface measured by actual measurements;

[0099] (1.2) Establish a vehicle-track coupling dynamics model to solve the wheel-rail dynamic response under rolling contact;

[0100] (1.3) Establish a subsurface strain calculation model and solve the temporal and spatial distribution characteristics of the rail strain rate under rolling contact.

[0101] Among them, in the numerical simulation of the three-dimensional geometric irregularity of the track surface measured in (1.1),

[0102] (1.1.1) In order to obtain the three-dimensional unevenness distribution data of rail hidden damage, the rail surface is scanned by a composite three-dimensional scanner (such as HandySCAN3D), and the data is processed by post-processing software. The specific steps are as follows:

[0103] ① Post marking points near the hidden damage of the rail to determine the spatial position of the rail, and calibrate the instrument accuracy based on the marking points;

[0104] ②Control the 3D scanner to move continuously to obtain the 3D information of the position where the laser passes, thereby forming continuous 3D track surface point cloud data;

[0105] ③ Import the 3D point cloud data into the post-processing software, select the reference object to determine the overall coordinate system of the point cloud data, and finally obtain the correct coordinates of the rail surface point cloud data;

[0106] ④ Interpolate the randomly distributed point cloud data and integrate the point cloud data into a matrix expression to ensure that the coordinates along the horizontal and vertical directions remain consistent, so as to introduce the track surface data considering three-dimensional unevenness into the subsequent dynamic model;

[0107] ⑤ Through the filter, the point cloud data is processed with high frequency along the longitudinal direction of the rail surface to eliminate the problem of fluctuation of point cloud data caused by the accuracy problem of the 3D scanner.

[0108] (1.1.2) The real three-dimensional geometric point cloud data containing hidden damage on the rail is discretized into grid units along the longitudinal direction, and the wheel is discretized into rotating body grid units along the longitudinal direction. The wheel-rail spatial geometry solution method is used to obtain the wheel-rail contact points and the three-dimensional asymmetric contact geometry gap distribution.

[0109] Among them, the vehicle-track coupling dynamic model is established in (1.2) to solve the wheel-rail dynamic response under rolling contact.

[0110] (1.2.1) The vehicle model considers the car body, frame, wheelset and suspension system. The car body and bogie are considered as rigid bodies, and the suspension system is considered as spring and damping units. The track model includes rails, fastener systems and ballast. The rails are simulated using the Euler beam model, and the fastener system is considered as a spring-damper unit.

[0111] (1.2.2) The wheel-rail contact geometric clearance obtained in step (1.1) is brought into the dynamic model. In this way, the wheel-rail contact model is used to solve the accurate wheel-rail rolling contact behavior under dynamic response, and the dynamic parameters required for each integral step are obtained, including the wheel-rail contact force and creep rate.

[0112] Among them, the above (1.3) establishes a subsurface strain calculation model and solves the temporal and spatial distribution characteristics of the rail strain rate under rolling contact.

[0113] (1.3.1) While calculating the wheel-rail dynamic response in (1.2), the strain field results at each integration step are calculated in real time. This is the innovation of this patent compared with the traditional dynamic model. First, the normal contact stress matrix p is obtained at each integration step. N and the tangential contact stress matrix s N , and the influence factor matrix D under normal stress is obtained by using the influence factor method N And the influence factor matrix T under the action of tangential stress N , the subscript N represents the number of grids divided in the contact area based on the boundary element method, and then the above matrix is ​​operated as follows for numerical calculation:

[0114] First, the contact stress matrix and the influence factor matrix are topologically converted into calculation matrices with twice the number of grids in the transverse and longitudinal directions, where the stress matrix p for the topology is 2N and 2N The elements are filled with 0 values, and the impact factor matrix D 2N and T 2N Elements are filled using circular convolution;

[0115] Then, the matrix is ​​transformed from the time domain to the frequency domain using discrete Fourier transform.

[0116] (1.3.2) Calculate the product of the topological normal contact stress matrix and the normal influence factor matrix in the frequency domain, and the product of the tangential contact stress matrix and the tangential influence factor matrix in the frequency domain, and then use the inverse Fourier transform to convert the calculation results in the frequency domain to the time domain:

[0117]

[0118] Where IFFT stands for inverse Fourier transform, ρ is the calculation result of topology in half space to twice the number of grids, q and r correspond to the x, y, and z directions respectively; then the matrix elements of the last N+1 items of the matrix topology are discarded to obtain the stress value σ corresponding to each rail target point qr , which contains 9 stress components, and then the strain component results can be calculated by generalized Hooke's law;

[0119] (1.3.3) The strain rate is calculated by the strain field results at each integration step. The specific steps are as follows:

[0120] ① Select the rail surface strain rate calculation point as the target point O, so that the wheel moves at a unit distance (unit distance = integration step Δt p × running speed v) moving;

[0121] ② Determine the calculation area along the longitudinal direction with the target point O as the center and divide the grid. When the wheel passes through the left boundary x of the calculation area L When the wheel passes through the right boundary x of the calculation area, the wheel-rail contact stress matrix is ​​interpolated to the target point coordinates to calculate the strain field. R Time, table

[0122] The wheel is off and the strain rate is no longer calculated;

[0123] ③ Divide the time history of the wheel passing through the calculation area into T discrete time points

[0124] t 1,2,…,T ∈(x L / v,x R / v), calculate the strain rate results at each point;

[0125] ④ The strain rate reference value of the target point O is the maximum value among a series of point strain rates.

[0126] The strain rate at point O is solved by the following formula:

[0127]

[0128] In the formula, Indicates t i The strain result at the moment. i+1 -t i It can be expressed as the strain rate calculation step size set in the dynamics calculation.

[0129] Wherein, in the step 2: calculating the local elastic-plastic contact stress between wheel and rail,

[0130] The dynamic parameters and strain rate results calculated by the dynamic model are used as input quantities, and the wheel-rail elastic-plastic rolling contact equation and sliding equation are proposed:

[0131]

[0132] Where P n and P τ represents normal and tangential stress; A is defined as the load-displacement influence coefficient; h is the normal gap when the two contact bodies are in rigid contact; ω is the sliding amount; g I

[0133] is the tangential contact boundary; and It manifests itself as plastic deformation in the normal and tangential directions during contact; c x 、c y , are the longitudinal, transverse and spin creep rates,

[0134] The contact equations (c) and (d) are extended to consider the influence of plastic deformation of the contact interface on the contact solution; the tangential sliding equation (e) is extended to consider the deformation rate caused by stress fluctuations in dynamic contact, which was ignored in previous calculations due to the steady-state assumption. These two items are extended to reflect the mutual coupling influence of contact stress and elastic-plastic deformation during wheel-rail impact and stress

[0135] Effect of fluctuations on the strain field.

[0136] Solve the above equations and establish the wheel-rail elastic-plastic rolling contact model to solve the elastic-plastic contact stress. The specific steps are as follows:

[0137] (2.1) Establishing elastic contact mechanics model to calculate local elastic solution

[0138] Using a graded loading method, the wheel-rail macro-dynamic parameters are introduced into the elastic contact model to calculate the local contact stress under the current load conditions, and the sub-surface stress-strain calculation model established in (1.1.3) is used to calculate the internal stress field results of the rail;

[0139] (2.2) Based on the Von Mises yield criterion, the plasticity occurrence area is determined and the plastic strain is solved:

[0140] The yield judgment formula (f) is established. When the yield function f>0, it means that the material in this area exceeds the yield limit and plastic deformation occurs; otherwise, it is still in elastic deformation.

[0141]

[0142] Where δ qr is the Kronecker symbol, σ y is the yield stress, and its value increases with the strain rate.

[0143] The strain rate calculated by (1.1.3) is used as the input parameter, combined with

[0144] The Johnson-Cook method establishes the material constitutive model related to strain rate at different times, and the radial return algorithm proposed by Simon and Hughes is used to solve the plastic strain;

[0145] (2.3) Using the eigenstrain method to solve the internal residual stress and surface plastic deformation

[0146] The internal area of ​​the rail is discretely divided into rectangular grid units of equal size, and the plastic characteristics in each unit are assumed to be constant. The influence factors characterizing the relationship between plastic strain and residual stress given by Liu et al. are used to solve the residual stress of each unit. The residual stress generated by each unit is superimposed to obtain the residual stress distribution in the plastic area. Finally, the plastic deformation of the rail surface is solved with the help of residual stress. The idea is similar to that of residual stress, and both are calculated with the help of the influence factor method.

[0147] (2.4) Obtaining wheel-rail elastic-plastic rolling contact stress

[0148] The surface plastic deformation ( and ) is fed back to equations (c) to (e), and the wave equation (g) is introduced. This equation is used to solve the size of the wave displacement under unit displacement excitation, so the plastic deformation can be introduced into equation (g), and the wave displacement u of the potential contact area can be obtained by using the superposition principle. Repeat steps (a) to (c) until the contact stress calculated by the two iterations converges, and it is considered that the current load is balanced. Continue to increase the load until the wheel-rail contact force calculated by (1.1) is reached, and the calculation can be terminated to obtain the elastic-plastic contact stress at the current moment.

[0149]

[0150] Where λ and μ are Lame constants, r and z represent radial and vertical coordinates, and ρ is the material density;

[0151] Wherein, the step 3: establishing a cyclic rolling three-dimensional local finite element model to solve the rail plastic stress-strain solution,

[0152] Establish a three-dimensional local finite element model such as Figure 1 As shown in the figure, longitudinal constraints are applied to both end faces of the rail, and full constraints are applied to the bottom of the rail. The distribution of the plastic stress-strain field of the rail needs to consider the constitutive relationship of the material, and the material of the rail is adjusted by modifying the UMAT subroutine in combination with the actual situation.

[0153] In this model, the translation of normal stress and tangential stress in the rolling direction is used to simulate the cyclic rolling process. Considering the non-steady-state rolling contact at the three-dimensional geometric unevenness of the rail, the normal stress and tangential stress both vary with time, and the load history is taken from the elastic-plastic contact stress distribution at each moment in step 1. The cyclic moving load on the rail surface is implemented in this patent by modifying the inp file of the finite element software through a self-compiled program. At the initial position point A, the contact load gradually increases through 10 incremental steps, and then the contact load gradually moves longitudinally to the end point B. The contact load moves a distance of 1mm in one incremental step, and the total rolling distance is 100mm. After 100 steps of rolling, the contact spot reaches point B and finally decreases to zero through 10 incremental steps. In order to save calculation costs without affecting the calculation results, five load cycles can be used in the calculation, and the plastic stress-strain field distribution of the rail can be obtained through this model.

[0154] Wherein, the step 4 uses the plastic stress-strain solution to solve the rail hidden crack initiation angle and crack initiation life.

[0155] (4.1) When predicting the crack initiation life of rails, an appropriate rolling contact fatigue life prediction model is selected according to the local material response of the rails, such as the Jiang-Sehitoglu model suitable for the ratchet effect. This model considers that fatigue damage is caused by the plastic strain energy in the critical plane of the material, and its fatigue damage increment can be expressed as:

[0156]

[0157] Where <> is the Macauley bracket; dD is the fatigue damage increment; σmr, σ0, σf are the material memory parameter, fatigue strength limit and true fracture stress respectively; σ and τ are the normal compressive stress and shear stress respectively; dε p and dγ p are the plastic strain increments corresponding to σ and τ respectively; b and m are material constants; dY is the plastic strain energy increment in the critical material plane;

[0158] (4.2) Use the critical plane method to calculate the maximum value of fatigue damage parameters and predict the initiation angle and crack life of rail hidden damage.

[0159] The corresponding fatigue damage parameter value can be calculated for each plane at any point in the material, and the plane where the maximum value of the fatigue damage parameter is located is the critical plane. According to the plastic stress-strain results of the rail under each cyclic load calculated by the finite element, the damage value generated at a certain point can be obtained. When the damage of a point on the critical surface reaches the critical damage value D0, the point is failed. And considering the influence of strain rate on ductility, a structure consistent with the site is obtained, and a correction equation related to strain rate is introduced in the critical damage value, such as formula (j):

[0160]

[0161] In the formula, and are strain rate and reference strain rate; m1 is an undetermined coefficient, which is determined by the phenomenon observed on site and indoor rolling test;

[0162] According to the stress state at one point, the plastic stress-strain result in step 3 is transformed into a tensor coordinate with a step size of 1° to obtain the stress-strain state on any plane at this point, and the fatigue damage parameter value of the plane is solved according to equations (h), (i) and (j). During the coordinate transformation, the engineering shear strain output by the finite element method should be converted into a tensor shear strain. At the same time, the fatigue damage parameter values ​​of all target points are calculated, and the proportion of the two components, normal stress-strain and shear stress-strain, is quantitatively analyzed to find the maximum value of the fatigue damage parameter. The plane angle corresponding to the maximum value of the fatigue damage material is the crack initiation angle;

[0163] The crack initiation life Nf is defined as:

[0164]

[0165] Where ΔD is the cumulative damage value under a cyclic loading.

[0166] The final rail hidden damage cracking life is obtained by solving the definition formula of cracking life.

[0167] Example

[0168] The invention method is explained by taking the three-dimensional unevenness of the rail surface measured on a certain line as an example.

[0169] Step 1: According to step 1.1, obtain the hidden damage point cloud data through the actual measurement composite 3D scanner, such as Figure 2a As shown; select the reference coordinate system to determine the point cloud coordinates, convert the irregular data points into a matrix form through interpolation, and obtain a regular surface model, such as Figure 2b As shown; the track surface data points are filtered to eliminate high-frequency fluctuations, such as Figure 2c shown.

[0170] Step 2: Establish a multi-body dynamics model as described in step 1.2 to calculate the wheel-rail macro contact force and creep rate distribution curve, as shown in Figure 3 shown.

[0171] Step 3: Establish the subsurface strain field calculation model according to step 1.3, and use formula (d) to calculate the strain rate variation curve of the hidden damage unevenness, such as Figure 4 shown.

[0172] Step 4: According to step 2.1, the calculation step of the elastic-plastic contact mechanics calculation model is set to one tenth of the total macroscopic wheel-rail force. The elastic contact model is used to calculate the elastic contact solution under the current load increment. The subsurface stress calculation model is used to calculate the equivalent stress distribution inside the rail. The plasticity occurrence area is determined through step 2.2.

[0173] Step 5: Based on the strain rate curve calculated in step 3, the material constitutive model under different strain rates is established in combination with the Johnson-Cook method, and the plastic strain is calculated using the intrinsic strain method; the residual stress and surface plastic deformation in the plastic occurrence area are solved using the intrinsic strain method, and finally the elastic-plastic contact stress distribution at each moment is calculated, such as Figure 5 shown.

[0174] Step 6: The elastic-plastic contact stress distribution at each moment calculated in Step 5 is imported into the three-dimensional local finite element model to calculate the final plastic stress-strain result, such as Figure 6 As shown, it is found that the rail material response at this time is a ratchet effect.

[0175] Step 7: The rail material response is ratchet effect, and the prediction model uses Jiang-Sehitoglu model. The calculated crack initiation angle is 50.3°, and the crack initiation life is 1724472 times (number of train passes).

[0176] Finally, it should be noted that the above embodiments are only used to describe the technical solution of the present invention rather than to limit the technical method. The present invention can be extended to other modifications, changes, applications and embodiments in application, and therefore it is believed that all such modifications, changes, applications and embodiments are within the spirit and teaching scope of the present invention.

Claims

1. A simulation prediction method for the initiation of hidden cracks in rails considering the three-dimensional geometric irregularity of the rail surface, characterized in that: The method comprises Step 1: Obtain the wheel-rail dynamic response and strain rate temporal and spatial distribution characteristics; Step 2: Calculate the local elastic-plastic contact stress between wheel and rail; Step 3: Establish a cyclic rolling three-dimensional local finite element model to solve the rail plastic stress-strain solution; Step 4 uses the plastic stress-strain solution to solve the rail hidden crack initiation angle and crack initiation life.

2. The method for simulating and predicting the initiation of hidden cracks in rails considering the three-dimensional geometric irregularity of the rail surface according to claim 1, characterized in that: In step 1, the wheel-rail dynamic The temporal and spatial distribution characteristics of mechanical response and strain rate include (1.1) Numerical simulation of three-dimensional geometric irregularities of track surface measured by actual measurements; (1.2) Establish a vehicle-track coupling dynamics model to solve the wheel-rail dynamic response under rolling contact; (1.3) Establish a subsurface strain calculation model and solve the temporal and spatial distribution characteristics of the rail strain rate under rolling contact.

3. The method for simulating and predicting the initiation of hidden cracks in rails considering the three-dimensional geometric irregularity of the rail surface according to claim 2, characterized in that: In the numerical simulation of the three-dimensional geometric irregularity of the track surface measured in (1.1), (1.1.1) In order to obtain the three-dimensional unevenness distribution data of rail hidden damage, the rail surface is scanned by a composite three-dimensional scanner, and the data is processed by post-processing software. The specific steps are as follows: ① Post marking points near the hidden damage of the rail to determine the spatial position of the rail, and calibrate the instrument accuracy based on the marking points; ②Control the 3D scanner to move continuously to obtain the 3D information of the position where the laser passes, thereby forming continuous 3D track surface point cloud data; ③ Import the 3D point cloud data into the post-processing software, select the reference object to determine the overall coordinate system of the point cloud data, and finally obtain the correct coordinates of the rail surface point cloud data; ④ Interpolate the randomly distributed point cloud data and integrate the point cloud data into a matrix expression to ensure that the coordinates along the horizontal and vertical directions remain consistent, so as to introduce the track surface data considering three-dimensional unevenness into the subsequent dynamic model; ⑤ Through the filter, the point cloud data is processed with high frequency along the longitudinal direction of the rail surface to eliminate the problem of fluctuation of point cloud data caused by the accuracy problem of the 3D scanner. (1.1.2) The real three-dimensional geometric point cloud data containing hidden damage on the rail is discretized into grid units along the longitudinal direction, and the wheel is discretized into rotating body grid units along the longitudinal direction. The wheel-rail spatial geometry solution method is used to obtain the wheel-rail contact points and the three-dimensional asymmetric contact geometry gap distribution.

4. The method for simulating and predicting the initiation of hidden cracks in rails considering the three-dimensional geometric irregularity of the rail surface according to claim 2, characterized in that: (1.2) Build the vehicle - Track coupling dynamics model, solving the wheel-rail dynamic response under rolling contact, (1.2.1) The vehicle model includes the car body, frame, wheelset and suspension. The car body and bogie are considered as rigid bodies, and the suspension system is considered as spring and damping units. The track model includes rails, fastener systems and ballast. The rails are simulated using the Euler beam model, and the fastener system is considered as a spring-damper unit. (1.2.2) The wheel-rail contact geometric clearance obtained in step (1.1) is brought into the dynamic model. In this way, the wheel-rail contact model is used to solve the accurate wheel-rail rolling contact behavior under dynamic response, and the dynamic parameters required for each integral step are obtained, including the wheel-rail contact force and creep rate.

5. The method for simulating and predicting the initiation of hidden cracks in rails considering the three-dimensional geometric irregularity of the rail surface according to claim 2, characterized in that: In (1.3) the subsurface strain calculation model is established and the temporal and spatial distribution characteristics of the rail strain rate under rolling contact are solved. (1.3.1) While calculating the wheel-rail dynamic response in (1.2), the strain field results at each integration step are calculated in real time. This is the innovation of this patent compared with the traditional dynamic model. First, the normal contact stress matrix p is obtained at each integration step. N and the tangential contact stress matrix s N , and the influence factor matrix D under normal stress is obtained by using the influence factor method N And the influence factor matrix T under the action of tangential stress N , the subscript N represents the number of grids divided in the contact area based on the boundary element method, and then the above matrix is ​​operated as follows for numerical calculation: First, the contact stress matrix and the influence factor matrix are topologically converted into calculation matrices with twice the number of grids in the transverse and longitudinal directions, where the stress matrix p for the topology is 2N and 2N The elements are filled with 0 values, and the impact factor matrix D 2N and T 2N The elements are filled by circular convolution; then, the matrix is ​​transformed from the time domain to the frequency domain using discrete Fourier transform. (1.3.2) Calculate the product of the topological normal contact stress matrix and the normal influence factor matrix in the frequency domain, and the product of the tangential contact stress matrix and the tangential influence factor matrix in the frequency domain, and then use the inverse Fourier transform to convert the calculation results in the frequency domain to the time domain: Where IFFT stands for inverse Fourier transform, ρ is the calculation result of topology in half space to twice the number of grids, q and r correspond to the x, y, and z directions respectively; then the matrix elements of the last N+1 items of the matrix topology are discarded to obtain the stress value σ corresponding to each rail target point qr , which contains 9 stress components, and then the strain component results can be calculated by generalized Hooke's law; (1.3.3) The strain rate is calculated by the strain field results at each integration step. The specific steps are as follows: ① Select the rail surface strain rate calculation point as the target point O, so that the wheel moves at a unit distance (unit distance = integration step Δt p × running speed v) moving; ② Determine the calculation area along the longitudinal direction with the target point O as the center and divide the grid. When the wheel passes through the left boundary x of the calculation area L When the wheel passes through the right boundary x of the calculation area, the wheel-rail contact stress matrix is ​​interpolated to the target point coordinates to calculate the strain field. R When , it means the wheel has left and the strain rate is no longer calculated; ③ Divide the time history of the wheel passing through the calculation area into T discrete time points t 1,2,…,T ∈(x L / v,x R / v), calculate the strain rate results at each point; ④ The strain rate reference value of the target point O is the maximum value among a series of point strain rates, so the strain rate at the target point O is solved by the following formula: In the formula, Indicates t i The strain result at the moment t i+1 -t i It can be expressed as the strain rate calculation step size set in the dynamics calculation.

6. The method for simulating and predicting the initiation of hidden cracks in rails considering the three-dimensional geometric irregularity of the rail surface according to claim 1, characterized in that: In step 2: calculating the local elastic-plastic contact stress between wheel and rail, The dynamic parameters and strain rate results calculated by the dynamic model are used as input quantities, and the wheel-rail elastic-plastic rolling contact equation and sliding equation are proposed: Where P n and P τ represents normal and tangential stress; A is defined as the load-displacement influence coefficient; h is the normal gap when the two contact bodies are in rigid contact; ω is the sliding amount; g I is the tangential contact boundary; and It manifests itself as plastic deformation in the normal and tangential directions during contact; c x 、c y , are the longitudinal, transverse and spin creep rates, The contact equations (c) and (d) are extended to take into account the influence of the plastic deformation of the contact interface on the contact solution; the tangential sliding equation (e) is extended to take into account the deformation rate caused by stress fluctuations in dynamic contact, which was ignored in previous calculations due to the steady-state assumption. Solve the above equations and establish the wheel-rail elastic-plastic rolling contact model to solve the elastic-plastic contact stress. The specific steps are as follows: (2.1) Establish an elastic contact mechanics model to calculate the local elastic solution. Use a graded loading method to bring the wheel-rail macrodynamic parameters into the elastic contact model to calculate the local contact stress under the current load conditions. Use the subsurface stress-strain calculation model established in (1.1.3) to calculate the internal stress field results of the rail. (2.2) Based on the Von Mises yield criterion, the plasticity occurrence area is determined and the plastic strain is solved: The yield judgment formula (f) is established. When the yield function f>0, it means that the material in this area exceeds the yield limit and plastic deformation occurs; otherwise, it is still in elastic deformation. Where δ qr is the Kronecker symbol, σ y is the yield stress, which increases with the strain rate. The strain rate calculated by (1.1.3) is used as the input parameter, and the material constitutive model related to the strain rate at different times is established in combination with the Johnson-Cook method. The radial return algorithm proposed by Simon and Hughes is used to solve the plastic strain. (2.3) Using the eigenstrain method to solve the internal residual stress and surface plastic deformation The internal area of ​​the rail is discretely divided into rectangular grid units of equal size, and the plastic characteristics in each unit are assumed to be constant. The influencing factors characterizing the relationship between plastic strain and residual stress given by Liu et al. are used to solve the residual stress of each unit. The residual stress generated by each unit is superimposed to obtain the residual stress distribution in the plastic area; finally, the plastic deformation of the rail surface is solved with the help of residual stress. (2.4) Obtaining wheel-rail elastic-plastic rolling contact stress The surface plastic deformation ( and ) is fed back to equations (c) to (e), and the wave equation of equation (g) is introduced. This equation is used to solve the size of the wave displacement under unit displacement excitation, so the plastic deformation can be introduced into equation (g), and the wave displacement u of the potential contact area can be obtained by using the superposition principle; repeat steps (a) to (c) until the contact stress calculated by the two iterations converges, and it is considered that the current load has reached equilibrium. Continue to increase the load until the wheel-rail contact force calculated by (1.1) is reached, and then the calculation can be terminated to obtain the elastic-plastic contact stress at the current moment; Where λ and μ are Lame constants, r and z represent radial and vertical coordinates, and ρ is the material density.

7. The method for simulating and predicting the initiation of hidden cracks in rails considering the three-dimensional geometric irregularity of the rail surface according to claim 1, characterized in that: Step 3: Establish a cyclic rolling three-dimensional local finite element model to solve the rail plastic stress-strain solution. A three-dimensional local finite element model is established, longitudinal constraints are imposed on both longitudinal ends of the rail, and full constraints are adopted on the bottom of the rail. The constitutive relationship of the material needs to be considered to obtain the distribution of the plastic stress-strain field of the rail.

8. The method for simulating and predicting the initiation of hidden cracks in rails considering the three-dimensional geometric irregularity of the rail surface according to claim 1, characterized in that: In step 4, the plastic stress-strain solution is used to solve the initiation angle of hidden cracks in the rail and the crack initiation life. (4.1) When predicting the crack initiation life of rails, an appropriate rolling contact fatigue life prediction model is selected according to the local material response of the rails, such as the Jiang-Sehitoglu model suitable for the ratchet effect. This model considers that fatigue damage is caused by the plastic strain energy in the critical plane of the material, and its fatigue damage increment can be expressed as: Where <> is the Macauley bracket; dD is the fatigue damage increment; σmr, σ0, σf are the material memory parameter, fatigue strength limit and true fracture stress respectively; σ and τ are the normal compressive stress and shear stress respectively; dε p and dγ p are the plastic strain increments corresponding to σ and τ respectively; b and m are material constants; dY is the plastic strain energy increment in the critical material plane; (4.2) Use the critical plane method to calculate the maximum value of fatigue damage parameters and predict the initiation angle and crack life of rail hidden damage. The corresponding fatigue damage parameter value can be calculated for each plane at any point in the material, where the plane with the maximum value of the fatigue damage parameter is the critical plane; the damage value generated at a certain point can be obtained according to the plastic stress-strain results of the rail under each cyclic load calculated by the finite element method. When the damage of a point on the critical surface reaches the critical damage value D0, the point fails; and considering the influence of strain rate on ductility, a structure consistent with the site is obtained, and a correction equation related to strain rate is introduced into the critical damage value, such as formula (j): In the formula, and are strain rate and reference strain rate; m1 is an undetermined coefficient, which is determined by the phenomenon observed on site and indoor rolling test; According to the stress state of a point, the plastic stress-strain result in step 3 is transformed into a tensor coordinate with a step size of 1° to obtain the stress-strain state on any plane of this point, and the fatigue damage parameter value of the plane is solved according to equations (h), (i) and (j). During the coordinate transformation, the engineering shear strain output by the finite element method should be converted into a tensor shear strain. At the same time, the fatigue damage parameter values ​​of all target points are calculated, and the proportions of the two components of normal stress-strain and shear stress-strain are quantitatively analyzed to find the maximum value of the fatigue damage parameter; the plane angle corresponding to the maximum value of the fatigue damage material is the crack initiation angle; The crack initiation life Nf is defined as: Where ΔD is the cumulative damage value under one cyclic loading, The final rail hidden damage cracking life is obtained by solving the definition formula of cracking life.

Images