A method for calculating wave propagation velocity of rigid contact network based on finite element simulation

Abstract

The present invention discloses a method for calculating the wave propagation velocity of a rigid contact network based on finite element simulation, specifically: according to the design parameters of the rigid contact network, a finite element model of the rigid contact network is established by using ANCF beam units and three-dimensional spring units; concentrated loads and global gravity loads are applied to the rigid contact network model, and its equilibrium state under the joint action of the gravity load and the concentrated load is calculated; the concentrated load is removed instantaneously, and the free vibration process of the rigid contact network is calculated by using the Newmark method; a spatiotemporal distribution diagram of wave propagation of the rigid contact network is constructed, and the wave propagation velocity is calculated by using a method of tracking the transfer trajectory of the wave peak value. The present invention provides a method for calculating the wave propagation velocity on a rigid contact network based on finite element simulation, which can avoid the field condition restrictions and high test costs faced by physical tests, and at the same time avoid the problem of difficulty in accurately tracking the vibration start time of each position in traditional methods.

Description

Technical Field

[0001] The present invention relates to the technical field of rigid contact network wave behavior simulation technology, in particular to a method for calculating the wave propagation speed of a rigid contact network based on finite element simulation. Background Art

[0002] The rigid catenary system is a crucial component for obtaining power for electrified trains in tunnels. Stable contact between the rigid catenary system and the pantograph is crucial for safe and stable train operation. With the continuous development of urban rail transit and mountain tunnel railways, the requirements for the operating speed of rigid catenary systems are constantly increasing. As operating speeds increase, the wave propagation behavior on the rigid catenary increasingly impacts the current-collecting quality of the catenary-pantograph system. When train speeds approach the wave propagation speed of the rigid catenary, wave resonance can even occur, seriously threatening the safe operation of the system. The wave propagation speed is an important indicator for measuring the wave propagation behavior of rigid catenary systems and a key parameter to consider when designing high-speed rigid catenary systems. However, due to high testing costs and limited on-site tunnel conditions, direct measurement of the wave propagation speed on actual lines is difficult. Therefore, it is necessary to propose a numerical method to calculate the wave propagation speed of rigid catenary systems to provide a technical reference for high-speed rigid catenary design. Summary of the Invention

[0003] To accurately extract the wave propagation velocity of a rigid catenary, this paper proposes a method for calculating the wave propagation velocity of a rigid catenary based on finite element simulation. The method first establishes a finite element model of the rigid catenary using the absolute nodal coordinate method (ANCF) beam element. A concentrated load is then applied to the model and its deformation state is calculated. The concentrated load is then instantaneously released, and the free vibration process of the rigid catenary is calculated. Finally, the spatiotemporal distribution of the rigid catenary wave is constructed and the wave propagation velocity is extracted.

[0004] The present invention provides a method for calculating the wave propagation velocity of a rigid contact network based on finite element simulation, which specifically includes the following steps:

[0005] Step 1: Based on the rigid catenary design parameters, a rigid catenary finite element model is established using ANCF beam elements and three-dimensional spring elements;

[0006] Step 2: Apply concentrated load and global gravity load to the rigid catenary finite element model, and calculate the equilibrium state of the rigid catenary under the combined action of gravity load and concentrated load;

[0007] Step 3: Based on the equilibrium state, the concentrated load is removed instantaneously, and the free vibration process of the rigid contact network is calculated using the Newmark method;

[0008] Step 4: Construct a spatiotemporal distribution diagram of the rigid contact network fluctuations and calculate the fluctuation propagation velocity by tracing the peak transfer trajectory of the fluctuations.

[0009] Furthermore, step 1 includes:

[0010] 1.1: Use ANCF beam elements to discretize the rigid catenary busbar and contact wire, and construct the ANCF beam element coordinates, stiffness matrix and mass matrix; the ANCF beam element coordinate format is:

[0011] e A =[e1 e2 e3 e4 e5 e6 e7 e8 e9 e 10 e 11 e 12 ] T

[0012] in,

[0013] e1=x| χ＝0 e2=y| χ＝0 e3=z| χ＝0

[0014] e7=x| χ＝L e8=y| χ＝L e9=z| χ＝L

[0015] Where x, y, and z are the coordinates of any point on the ANCF beam element in three-dimensional space, L is the original length of the element, and χ is a local variable on the element that varies from 0 to L. The coordinates of any point on the ANCF beam element can be calculated using the following formula:

[0016] r = [xyz] T =Se A

[0017] S is the shape function of the ANCF beam element, written as:

[0018] S=[s1I3 s2I3 s3I3 s4I3]

[0019] Where I3 is the third-order unit matrix, s i for:

[0020]

[0021] Where ξ = χ / L.

[0022] Stiffness matrix K of ANCF beam element A and the mass matrix M A They are:

[0023]

[0024] K A =K L +K T

[0025] in,

[0026]

[0027]

[0028] Where m A is the mass of the ANCF beam element, E is Young's elastic modulus, A is the cross-sectional area of ​​the ANCF beam element, I is the section moment of inertia of the ANCF beam element, and ε is the axial strain of the ANCF beam element;

[0029] 1.2: Use three-dimensional linear spring elements to discretize the suspension structure and construct the mass matrix M of the three-dimensional spring element S and the stiffness matrix K S :

[0030]

[0031]

[0032] Among them, m s is the equivalent mass of the three-dimensional spring element, k x , k y , k z are the stiffness coefficients of the three-dimensional spring unit in the directions of the X, Y, and Z coordinate axes respectively; since the base end of the three-dimensional spring unit is fixed and the movable end is connected to the ANCF beam unit, the coordinates of the three-dimensional spring unit are:

[0033] e S =[xyz] T

[0034] 1.3: Select the unit type for each unit of the rigid contact network according to its actual structure, assemble the coordinate vectors and mass matrices of all units, and construct the system coordinate vector e and system mass matrix M of the rigid contact network. The unit types include ANCF beam units and three-dimensional spring units.

[0035] Furthermore, step 2 includes the following steps:

[0036] 2.1: Calculate the element concentrated load vector F caused by the concentrated load on the ANCF beam element C :

[0037] F C =S' T [0 F C0] T

[0038] Among them, F C is the concentrated load size, S' is the shape function matrix of the ANCF beam unit subjected to the concentrated load; then, according to the corresponding relationship between the degrees of freedom of each unit of the rigid contact network, the unit concentrated load vector is assembled into the concentrated load vector P of the entire system. C ;

[0039] 2.2: Calculate the element gravity load vector F caused by gravity acting on the ANCF beam element G :

[0040]

[0041] Among them, g is the gravity constant and ρ is the material density. Then, according to the corresponding relationship between the degrees of freedom of each unit of the rigid contact network, the unit gravity load vector is assembled into the overall gravity load vector P of the system. G ;

[0042] 2.3: Calculate the elastic internal force vector of the rigid catenary system. The elastic internal force vector of the rigid catenary system is composed of the elastic internal force of the ANCF beam unit and the elastic internal force vector of the three-dimensional spring unit;

[0043] The calculation formula of the elastic internal force of the ANCF beam element is:

[0044] Q A =K A e A

[0045] The calculation formula of the elastic internal force of the three-dimensional spring element is:

[0046] Q S =K S (e S -e0)

[0047] Among them, e0 is the initial coordinate vector of the spring unit; according to the corresponding relationship between the degrees of freedom of each unit of the rigid contact network, the elastic internal force vector Q of the rigid contact network system is assembled;

[0048] 2.4: Calculate the tangent stiffness matrix of the rigid catenary system. The tangent stiffness matrix of the rigid catenary system is composed of the tangent stiffness matrix of the ANCF beam unit and the tangent stiffness matrix of the three-dimensional spring unit. The tangent stiffness matrix of the ANCF beam unit is

[0049]

[0050] The tangent stiffness matrix of the three-dimensional spring element and the stiffness matrix K of the three-dimensional spring element SThe same; according to the corresponding relationship of the degree of freedom of each unit of the rigid contact network, the unit tangent stiffness matrix is ​​assembled to obtain the rigid contact network system tangent stiffness matrix J;

[0051] 2.5: Calculate the displacement increment ΔU of the rigid contact network under the unbalanced action of the system elastic internal force and external load. The calculation formula is:

[0052] ΔU=J -1 (P C +P G -Q)

[0053] 2.6: If max(ΔU)≤the first convergence tolerance, output the system coordinate vector e of the rigid contact network in equilibrium and proceed to the next step; otherwise, update the rigid contact network system coordinate vector: e=e+ΔU, and go to step 2.3.

[0054] Furthermore, step 3 includes the following steps:

[0055] 3.1: Calculate the intermediate matrices a1, a2, and a3 required for Newmark iteration. The calculation formula is:

[0056]

[0057] Where Δt is the Newmark iteration time step, β and γ are Newmark parameters; C is the system damping matrix, which can be calculated by Rayleigh damping;

[0058] 3.2: Initialize the Newmark iteration step index i = 1 and initialize the system coordinate vector e i , velocity vector and the acceleration vector Among them, since the rigid contact network is still in a static state, the initial velocity vector is a zero vector, and the calculation formula of the initial acceleration vector is:

[0059]

[0060] 3.3: Preset e i+1 =e i .

[0061] 3.4: Calculate the system tangent stiffness matrix J i+1 and elastic internal force vector Q i+1 ;

[0062] 3.5: Calculate the system coordinate increment, the formula is

[0063]

[0064] 3.6: Determine if max(Δe i+1)≤the second convergence tolerance, then the output system coordinates e i+1 , and proceed to step 3.7; otherwise, update e i+1 =e i+1 +Δe i+1 , and return to step 3.4;

[0065] 3.7: If i+1 < the total number of Newmark iterations, then according to the Newmark assumption, update and The formula is

[0066]

[0067] Then set i=i+1 and return to step 3.3; otherwise, complete the calculation and exit the loop.

[0068] Furthermore, step 4 includes the following steps:

[0069] 4.1: Assemble the vertical displacement data of the rigid contact network during the free vibration process obtained in step 3 into a wave propagation space-time distribution diagram according to spatial distribution and time sequence; the vertical displacement data is the system coordinate vector e output in step 3 i+1 The vertical coordinate difference between the system coordinate vector e output in step 2;

[0070] 4.2: Extract the displacement peak point in the fluctuation spatiotemporal distribution diagram;

[0071] 4.3: Use a linear function to fit the trajectory of the peak points, calculate the time-space slope of each trajectory line and take the average value, which is the wave propagation speed.

[0072] Beneficial technical effects of the present invention

[0073] 1. The method for calculating the wave propagation velocity of the rigid contact network proposed in the present invention can accurately calculate the wave propagation velocity of the rigid contact network based on finite element simulation, which can avoid the field condition restrictions and high test costs faced by physical tests.

[0074] 2. The present invention uses ANCF beam elements and three-dimensional spring elements to model the contact network, busbar, and suspension structure respectively. After modeling, they are assembled according to actual conditions, making the constructed finite element model more accurate;

[0075] 3. The present invention simulates the equilibrium state under a concentrated load and the free vibration process under instantaneous removal of the concentrated load, recording the system coordinates in the equilibrium state and the system coordinates at the position of maximum vibration during free vibration to obtain the peak displacement. The wave propagation velocity is then determined by tracing the displacement peak point's transfer trajectory. The method of the present invention achieves high computational accuracy, avoiding the difficulties encountered in traditional methods in accurately tracking the onset of vibration at each location and the susceptibility of the results to dispersion interference. This invention can provide methodological support for the study of rigid contact network dynamics. BRIEF DESCRIPTION OF THE DRAWINGS

[0076] Figure 1 This is a flow chart of a method for calculating the wave propagation velocity of a rigid contact network based on finite element simulation according to an embodiment of the present invention;

[0077] Figure 2 Schematic diagram of the ANCF beam unit used in an embodiment of the present invention;

[0078] Figure 3 A schematic diagram of a three-dimensional spring unit used in an embodiment of the present invention;

[0079] Figure 4 This is the equilibrium state of the rigid contact network under the action of gravity load and concentrated load in the embodiment of the present invention;

[0080] Figure 5 This is the spatiotemporal distribution diagram of the wave propagation and the wave peak trajectory of the rigid contact network in an embodiment of the present invention. DETAILED DESCRIPTION

[0081] The present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments. This embodiment is implemented based on the technical solution of the present invention and provides a detailed implementation process, but the protection scope of the present invention is not limited to the embodiment.

[0082] This embodiment provides a method for calculating the wave propagation velocity of a rigid contact network based on finite element simulation. The specific implementation process is as follows: Figure 1 As shown:

[0083] Step 1: Based on the rigid catenary design parameters, a rigid catenary finite element model is established using ANCF beam elements and three-dimensional spring elements;

[0084] 1.1: Use ANCF beam elements to discretize the rigid catenary busbar and contact wire, and construct the element coordinates, stiffness matrix and mass matrix of the ANCF beam element;

[0085] The coordinate format of ANCF beam element is:

[0086] e A =[e1 e2 e3 e4 e5 e6 e7 e8 e9 e10 e 11 e 12 ] T

[0087] in,

[0088] e1=x| χ＝0 e2=y| χ＝0 e3=z χ＝0

[0089] e7=x| χ＝L e8=y| χ＝L e9=z χ＝L

[0090] Where x, y, and z are the coordinates of any point on the ANCF beam element in three-dimensional space, L is the original length of the ANCF beam element, and χ is a local variable on the ANCF beam element that varies from 0 to L. The coordinates of any point on the ANCF beam element can be calculated using the following formula:

[0091] r = [xyz] T =Se A

[0092] S is the shape function of the ANCF beam element, written as:

[0093] S=[s1I3 s2I3 s3I3 s4I3]

[0094] Where I3 is the third-order unit matrix, s i for:

[0095]

[0096] Where ξ=χ / L;

[0097] Stiffness matrix K of ANCF beam element A and the mass matrix M A They are:

[0098]

[0099] K A =K L +K T

[0100] in,

[0101]

[0102]

[0103] Where m Ais the mass of the ANCF beam element, E is Young's elastic modulus, A is the cross-sectional area of ​​the ANCF beam element, I is the section moment of inertia of the ANCF beam element, and ε is the axial strain of the ANCF beam element;

[0104] 1.2: Use three-dimensional spring elements to discretize the suspension structure and construct the mass matrix M of the three-dimensional spring element S and the stiffness matrix K S :

[0105]

[0106]

[0107] Among them, m s is the equivalent mass of the three-dimensional spring unit, that is, the equivalent mass of the suspension structure unit, k x , k y , k z are the stiffness coefficients of the three-dimensional spring unit in the XYZ coordinate axis directions respectively; since the base end of the three-dimensional spring unit is fixed and the active end is connected to the ANCF beam unit, the coordinates of the three-dimensional spring unit can be directly written as

[0108] e S =[xyz] T

[0109] 1.3: Select the unit type for each unit of the rigid contact network according to its actual structure, assemble the coordinate vectors and mass matrices of all units, and construct the system coordinate vector e and system mass matrix M of the rigid contact network. The unit types include ANCF beam units and three-dimensional spring units.

[0110] It is equivalent to dividing the rigid contact network into several units, selecting the matching ANCF beam unit or three-dimensional spring unit according to the actual structure of the rigid contact network, and then assembling the entire rigid contact network using the ANCF beam unit or three-dimensional spring unit according to the connection relationship of each unit.

[0111] Since it is a prior art in this field to use a built unit model for assembly and calculation based on the actual structure of the rigid contact network, it will not be described in detail here.

[0112] In this embodiment, the contact network, busbar, and suspension structure are modeled respectively by ANCF beam units and three-dimensional spring units. After modeling, they are assembled according to actual conditions, making the constructed finite element model more accurate.

[0113] Step 2: Apply concentrated load and global gravity load to the rigid catenary model, and calculate its equilibrium state under the combined action of gravity load and concentrated load.

[0114] In this embodiment, concentrated loads are applied to the contact network and busbars, and the self-gravity of the contact network and busbars is added to calculate the equilibrium state under the combined action of the gravity load and the concentrated load; then (step 3), the concentrated load is instantly removed to calculate the free vibration of the rigid contact network.

[0115] 2.1: Calculate the element concentrated load vector F caused by the concentrated load on the ANCF beam element C , which can be expressed as

[0116] F C =S' T [0 F C 0] T

[0117] Among them, F C is the concentrated load size, S' is the shape function matrix of the ANCF beam unit subjected to the concentrated load, and the shape function matrix here is the same as the shape function matrix in step 1.1; then, according to the corresponding relationship between the degrees of freedom of each unit of the rigid contact network, the unit concentrated load vector is assembled into the concentrated load vector P of the entire system C ; This step is to assemble the concentrated load applied to an ANCF beam element into the entire system.

[0118] 2.2: Calculate the element gravity load vector F caused by gravity acting on the ANCF beam element G , and its calculation formula is

[0119]

[0120] Among them, g is the gravity constant and ρ is the material density. Then, according to the corresponding relationship between the degrees of freedom of each unit of the rigid contact network, the unit gravity load vector is assembled into the overall gravity load vector P of the system. G ; This step is to calculate the overall gravity of the contact network and busbar in the rigid contact network system.

[0121] 2.3: Calculate the elastic internal force vector of the rigid catenary system. The elastic internal force vector of the rigid catenary system is composed of the elastic internal force vector of the ANCF beam unit and the elastic internal force vector of the three-dimensional spring unit; the elastic internal force Q of the ANCF beam unit A The calculation formula is

[0122] Q A =K A e A

[0123] Elastic internal force Q of three-dimensional spring element S The calculation formula is

[0124] Q S =K S (e S-e0)

[0125] Among them, e0 is the initial coordinate vector of the spring unit; according to the corresponding relationship between the degrees of freedom of each unit of the rigid contact network, the elastic internal force vector Q of the rigid contact network system is assembled;

[0126] 2.4: Calculate the tangent stiffness matrix of the rigid catenary system. The tangent stiffness matrix of the rigid catenary system is composed of the tangent stiffness matrix of the ANCF beam unit and the tangent stiffness matrix of the three-dimensional spring unit. The tangent stiffness matrix of the ANCF beam unit is

[0127]

[0128] The tangent stiffness matrix of the three-dimensional spring element and the stiffness matrix K of the three-dimensional spring element S The same; according to the corresponding relationship of the degrees of freedom of each unit of the rigid contact network, the unit tangent stiffness matrix is ​​assembled to obtain the tangent stiffness matrix J of the rigid contact network system;

[0129] 2.5: Calculate the displacement increment ΔU of the rigid contact network under the unbalanced action of elastic internal force and external load. The calculation formula is:

[0130] ΔU=J -1 (P C +P G -Q) (18)

[0131] 2.6: If max(ΔU) ≤ the first convergence tolerance, output the system coordinate vector e of the rigid catenary in equilibrium and proceed to the next step. Otherwise, update the rigid catenary system coordinate vector: e = e + ΔU, and go to step 2.3. That is, in equilibrium, record the rigid catenary system coordinate vector e.

[0132] Step 3: Based on the equilibrium shape of the rigid contact network calculated in step 2, the concentrated load is removed instantaneously, and the free vibration process of the rigid contact network is calculated using the Newmark method.

[0133] 3.1: Calculate the intermediate matrices a1, a2, and a3 required for Newmark iteration. The calculation formula is:

[0134]

[0135] Where Δt is the Newmark iteration time step, β and γ are Newmark parameters. Preferably, β and γ can be set to 0.25 and 0.5, respectively. Those skilled in the art can also choose other values ​​based on actual conditions. M is the system mass matrix, and C is the system damping matrix, which can be calculated using Rayleigh damping.

[0136] 3.2: Initialize the Newmark iteration step index i = 1 and initialize the system coordinate vector e i , velocity vector and the acceleration vector Among them, since the rigid contact network is still in a static state, the initial velocity vector is a zero vector, and the calculation formula of the initial acceleration vector is:

[0137]

[0138] 3.3: Preset e i+1 =e i .

[0139] 3.4: Calculate the system tangent stiffness matrix J i+1 and elastic internal force vector Q i+1 .

[0140] 3.5: Calculate the system coordinate increment, the formula is

[0141]

[0142] 3.6: Determine if max(Δe i+1 )≤the second convergence tolerance, then the output system coordinates e i+1 , and proceed to step 3.7; otherwise, update e i+1 =e i+1 +Δe i+1 , and return to step 3.4.

[0143] It is worth noting that the first convergence tolerance and the second convergence tolerance may be the same or different, and those skilled in the art may select one according to actual needs.

[0144] 3.7: If i+1 < the total number of Newmark iterations, then according to the Newmark assumption, update and The formula is

[0145]

[0146] Then set i=i+1 and return to step 3.3; otherwise, complete the calculation and exit the loop.

[0147] Step 4: Construct a spatiotemporal distribution diagram of the rigid contact network fluctuations and calculate the fluctuation propagation velocity by tracing the peak transfer trajectory of the fluctuations.

[0148] 4.1: Assemble the vertical displacement data of the rigid contact network during the free vibration process obtained in step 3 into a wave propagation space-time distribution diagram according to spatial distribution and time sequence; the vertical displacement data is the system coordinate vector e output in step 3 i+1The vertical coordinate difference between the system coordinate vector e output in step 2;

[0149] 4.2: Extract the displacement peak point of the fluctuation space-time distribution;

[0150] 4.3: Use a linear function to fit the trajectory of the peak points, calculate the time-space slope of each trajectory line and take the average value, which is the wave propagation speed.

[0151] In another embodiment of the present invention, a finite element simulation is performed on a 61-span straight-segment rigid contact network to calculate the wave propagation velocity of the rigid contact network.

[0152] Specifically, a rigid contact network finite element model is established according to the method of the aforementioned embodiment, and the specific parameters are shown in Table 1.

[0153] A concentrated downward load of 200N is applied to the mid-span of the 31st span. The equilibrium state under the combined action of gravity load and concentrated load is as follows: Figure 4 After that, the concentrated load is removed instantly, the Newmark iteration time step is 0.001s, and the method of the above embodiment is used to construct the spatial and temporal distribution diagram of the rigid contact network, as shown in FIG. Figure 5 As shown (for the convenience of display, the attached Figure 5 The displacements were normalized).

[0154] Table 1 Rigid contact network parameters in the embodiment

[0155]

[0156] Extraction of attached Figure 5 The displacement peak point is shown as the "△" point in the figure. The figure diverges to the left and right obliquely upward directions with the white bright spot at the concentrated load position as the center. The divergent direction is the propagation direction of the vibration wave. The line connecting the displacement peak points in the same propagation direction is the peak point connection trajectory.

[0157] A linear function is used to fit the peak point connection trajectory, and the time-space slope of each trajectory line is calculated and the average value is taken. Figure 5 The time-space slope is calculated to be 81.66 based on the data, that is, the wave propagation speed of the rigid contact network in this embodiment is 81.66 m / s.

[0158] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions and improvements made within the spirit and principles of the present invention should be included in the scope of protection of the present invention.

Claims

1. A method for calculating the wave propagation velocity of a rigid contact network based on finite element simulation, characterized in that: The following steps are involved: Step 1: According to the rigid contact network design parameters, the rigid contact network finite element model is established using ANCF beam elements and three-dimensional spring elements; Step 2: Apply concentrated load and global gravity load to the rigid contact network finite element model, and calculate the equilibrium state of the rigid contact network under the combined action of gravity load and concentrated load; Step 3: Based on the equilibrium form, the concentrated load is removed instantaneously, and the free vibration process of the rigid contact network is calculated using the Newmark method; Step 4: Construct the spatiotemporal distribution diagram of the rigid contact network fluctuations, and calculate the fluctuation propagation speed by tracking the peak transfer trajectory of the fluctuations; The step 3 is as follows: 3.1: Calculate the intermediate matrix required for Newmark iteration , and , the calculation formula is: ， in, is the Newmark iteration time step, and is the Newmark parameter; is the system damping matrix, which can be calculated by Rayleigh damping; M is the system mass matrix; 3.2: Initialize Newmark iteration step index , initialize the system coordinate vector , velocity vector and the acceleration vector ; Among them, since the rigid contact network is still in a static state, the initial velocity vector is a zero vector, and the calculation formula for the initial acceleration vector is: ， 3.3: Presets ; 3.4: Calculate the system tangent stiffness matrix and the elastic internal force vector ; 3.5: Calculate the system coordinate increment, the formula is ， 3.6: Determine if The second convergence tolerance, the output system coordinate vector , and proceed to step 3.7; otherwise, update , and return to step 3.4; 3.7: Determine if The total number of Newmark iterations is then updated according to the Newmark assumption. and , the formula is ， Then make And return to step 3.3; otherwise, complete the calculation and exit the loop; The step 4 is specifically as follows: 4.1: The vertical displacement data of the rigid contact network during the free vibration process obtained in step 3 are assembled into a wave propagation space-time distribution diagram according to spatial distribution and time sequence; the vertical displacement data is the system coordinate vector output in step 3 The system coordinate vector output from step 2 The vertical coordinate difference between ; 4.2: Extract the displacement peak point in the fluctuation space-time distribution diagram; 4.3: Use a linear function to fit the trajectory of the peak points, calculate the time-space slope of each trajectory line and take its average value, which is the wave propagation speed.

2. The method for calculating the wave propagation velocity of a rigid contact network based on finite element simulation according to claim 1 is characterized in that: Step 1 includes: 1.1: Use ANCF beam elements to discretize the rigid contact network busbar and contact wire, and construct the unit coordinates, stiffness matrix and mass matrix of the ANCF beam element; The coordinate format of ANCF beam element is: ， in, ， in, , , is the coordinate of any point on the ANCF beam element in three-dimensional space, is the original length of the ANCF beam element, The ANCF beam element is from 0 to The coordinates of any point on the ANCF beam element can be calculated by the following formula: ， is the shape function of the ANCF beam element, written as: ， in is the third-order identity matrix, for: ， in ; Stiffness matrix of ANCF beam element and the mass matrix They are: ， ， in, ， ， In the formula, is the mass of the ANCF beam element, E is Young's elastic modulus, A is the cross-sectional area of ​​the ANCF beam element, I is the section moment of inertia of the ANCF beam element, is the axial strain of the ANCF beam element; 1.2: Use three-dimensional spring units to discretize the suspension structure and construct the mass matrix of the three-dimensional spring unit and the stiffness matrix : ， ； in, is the equivalent mass of the three-dimensional spring element, , , are the stiffness coefficients of the three-dimensional spring unit in the directions of the three coordinate axes XYZ; since the base end of the three-dimensional spring unit is fixed and the active end is connected to the ANCF beam unit, the coordinates of the three-dimensional spring unit are ； 1.3: Select the unit type of each unit of the rigid contact network according to its actual structure, assemble the coordinate vectors and mass matrices of all units, and construct the system coordinate vector of the rigid contact network and the system mass matrix , the unit types include ANCF beam unit and three-dimensional spring unit.

3. The method for calculating the wave propagation velocity of a rigid contact network based on finite element simulation according to claim 1 or 2, characterized in that: The step 2 is specifically as follows: 2.1: Calculate the element concentrated load vector caused by the concentrated load on the ANCF beam element : ， in, is the concentrated load size, S' is the shape function matrix of the ANCF beam unit subjected to the concentrated load; then, according to the corresponding relationship between the degrees of freedom of each unit of the rigid contact network, the unit concentrated load vector is assembled into the concentrated load vector of the whole system ; 2.2: Calculate the element gravity load vector caused by gravity acting on the ANCF beam element : ， in, is the gravitational constant, is the material density; then, according to the corresponding relationship between the degrees of freedom of each unit of the rigid contact network, the unit gravity load vector is assembled into the overall gravity load vector of the system ; 2.3: Calculate the elastic internal force vector of the rigid catenary system. The elastic internal force vector of the rigid catenary system is composed of the elastic internal force vector of the ANCF beam unit and the elastic internal force vector of the three-dimensional spring unit; the elastic internal force Q of the ANCF beam unit A The calculation formula is: ， Elastic internal force Q of three-dimensional spring element S The calculation formula is: ， in, is the initial coordinate vector of the spring unit; according to the corresponding relationship between the degrees of freedom of each unit of the rigid contact network, the elastic internal force vector of the rigid contact network system is assembled ; 2.4: Calculate the tangent stiffness matrix of the rigid catenary system. The tangent stiffness matrix of the rigid catenary system is composed of the tangent stiffness matrix of the ANCF beam unit and the tangent stiffness matrix of the three-dimensional spring unit; the tangent stiffness matrix of the ANCF beam unit is ， Tangent stiffness matrix of three-dimensional spring element and stiffness matrix of three-dimensional spring element Same; according to the corresponding relationship of the unit degrees of freedom of the rigid contact network, the unit tangent stiffness matrix is ​​assembled to obtain the tangent stiffness matrix of the rigid contact network system ; 2.5: Calculation of displacement increment of rigid contact network under the unbalanced action of elastic internal force and external load of the system , the calculation formula is: ， 2.6: Determine if , then the system coordinate vector of the rigid catenary under equilibrium state is output And proceed to the next step; otherwise, update the rigid contact network system coordinate vector: , and go to step 2.3.

Images