Accurate calculation method of three-dimensional contact force between pantograph and rigid catenary

Abstract

The application discloses a kind of accurate calculation methods of three-dimensional contact force between pantograph and rigid catenary, specifically: using ANCF beam unit and three-dimensional linear bar unit to establish rigid catenary finite element model, using reduction mass block method to establish pantograph model, pantograph-rigid catenary motion equation is constructed;Using central difference method to predict the next time pantograph-rigid catenary system motion state, if prediction state is contact, then using Newton iteration method to solve the accurate three-dimensional contact force of current time;If prediction state is offline, then current time contact force is zero;Judge whether to complete the calculation of all time steps, if it has been completed, then output bow net contact force sequence and end;If not completed, then update system state, update iteration time step index, return prediction state judgment.The application can obtain accurate three-dimensional contact force and displacement result of bow net system, and can effectively avoid the difficulty in selecting contact stiffness and the problem of virtual penetration in traditional penalty function method.

Description

Technical Field

[0001] This invention belongs to the field of dynamic simulation technology for electrified railways, and in particular relates to a method for accurately calculating the three-dimensional contact force between the pantograph and the rigid contact wire. Background Technology

[0002] The pantograph-rigid contact wire system is an important way for electrified trains in tunnels to obtain electrical energy. At present, in the field of pantograph-contact wire dynamics research, the finite element method is mainly used to establish the pantograph-rigid contact wire system model, and the motion equations of the model are solved by the step-by-step integration method. C. Vera et al. (C. Vera, B. Suarez, J. Paulin, et al. Simulation model for the study of overhead rail current collector systems dynamics, focused on the design of a new conductor rail[J]. Veh. Syst. Dyn. 2006, 44: 595-614) established a finite element model of the rigid contact wire using the finite element software ANSYS. A. Bautista et al. (A. Bautista, J. Montesinos, P. Pintado. Dynamic interaction between pantograph and rigid overhead lines using a coupled FEM-Multibody procedure[J]. Mech. Mach. Theory. 2016, 97: 100-111) established a finite element model of the rigid contact network using beam elements and simulated the dynamic contact between the pantograph and the catenary. Zhou Ning et al. (Zhou Ning, Zou Huan, Zou Dong, et al. Adaptability study of simulation model of urban rail transit pantograph-catenary system[J]. Journal of Southwest Jiaotong University, 2017, 52(2): 408-423) compared the variable stiffness model and the finite element model of the rigid contact network and proved that the results obtained by the finite element model are more accurate under high-speed conditions. L. Chen et al. (L. Chen, F. Duan, Y. Song, et al. Assessment of dynamic interaction performance of high-speed pantograph and overhead conductor rail system[J]. IEEE Trans. Instrum. Meas. 2022, 71: 1-14) established a finite element model of a rigid contact network using the absolute nodal coordinate method (ANCF) and solved the dynamic contact force of the pantograph-catenary system using the Newmark method. Existing modeling and simulation studies of pantograph-rigid contact network systems generally use the penalty function method to simulate the pantograph-catenary contact. However, the contact stiffness required by this method is an empirical parameter, making it difficult to determine a suitable value. Furthermore, the penalty function method can cause virtual penetration of the pantograph-catenary contact pair, leading to inaccurate system displacement.Furthermore, current simulation methods only define the contact force in the vertical direction, without considering the change in the direction of the contact force in three-dimensional space caused by pantograph-catenary vibration, or the influence of tangential friction. Summary of the Invention

[0003] To obtain accurate three-dimensional dynamic contact force between the pantograph and the rigid contact network, this invention provides a precise calculation method for the three-dimensional contact force between the pantograph and the rigid contact network.

[0004] The present invention provides a method for accurately calculating the three-dimensional contact force between a pantograph and a rigid contact wire, comprising the following steps:

[0005] Step 1: Based on the design parameters of the rigid contact network, establish a finite element model of the rigid contact network using ANCF beam elements and three-dimensional linear rod elements. Based on the pantograph parameters, establish a pantograph model using the reduced mass block method, and construct the motion equations of the pantograph-rigid contact network.

[0006] S11: The rigid contact wire busbar and contact wire are discretized using ANCF beam elements. The coordinates of a single element are:

[0007] e A =[e1 e2 e3 e4 e5 e6 e7 e8 e9 e 10 e 11 e 12 ] T (1)

[0008] in,

[0009]

[0010] Where L is the original length of the element, x is a local variable on the element that varies from 0 to L, and r1, r2, and r3 are the coordinates of any point on the element in three-dimensional space, calculated by the following formula:

[0011] r = [r1 r2 r3] T =Se A (3)

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

[0013] S=[s1I3 s2I3 s3I3 s4I3] (4)

[0014] Where I3 is a third-order identity matrix, s i for:

[0015]

[0016] Where ξ = x / L.

[0017] Stiffness matrix K of ANCF element A and mass matrix M A They are respectively:

[0018]

[0019] K A =K L +K T (7)

[0020] in

[0021]

[0022]

[0023] In the formula, m A Let E be the mass of the element, E be the Young's modulus, A be the cross-sectional area of ​​the ANCF beam element, I be the moment of inertia of the ANCF beam element, and ε be the axial strain of the beam element.

[0024] S12: Employs a three-dimensional linear spring element discrete suspension structure with a mass matrix M S and stiffness matrix K S They are respectively:

[0025]

[0026]

[0027] Where, m S k is the equivalent mass of the suspension structure. x k y k z These are the stiffness coefficients of the spring element in the X, Y, and Z coordinate axes, respectively.

[0028] S13: Based on the element connection relationship, assemble the stiffness matrix of the ANCF beam element and the three-dimensional spring element to construct the stiffness matrix and mass matrix of the rigid contact network as a whole.

[0029] S14: Establish the reduced mass model of the pantograph, with the following equations:

[0030]

[0031] In the formula, m1, m2, and m3 are the masses of the bow head, upper frame, and lower frame, respectively; c1, c2, and c3 are the damping of the bow head, upper frame, and lower frame, respectively; k1, k2, and k3 are the stiffness of the bow head, upper frame, and lower frame, respectively; y1, y2, and y3 are the displacements of the bow head, upper frame, and lower frame, respectively; F C and F LThese are the contact force between the pantograph and the catenary, and the static lifting force on the pantograph, respectively.

[0032] S15: The kinematic equations of the pantograph-rigid contact wire system are as follows:

[0033]

[0034] In the formula, Let M and C be the acceleration and velocity vectors of the system as a whole, respectively; let P be the external load vector acting on the system; and let Q be the elastic internal force vector of the system as a whole. The formula for calculating the elastic internal force of the ANCF beam element is:

[0035] Q A =K A e A (14)

[0036] The method for calculating the internal forces of a three-dimensional linear spring element is as follows:

[0037] Q S =K S U S (15)

[0038] Among them, U S Let be the displacement vector of the spring element.

[0039] Step 2: Use the central difference method to predict the motion state of the pantograph-rigid contact network system at the next moment. If the predicted state is "contact", proceed to step 3. If the predicted state is "offline", the contact force in the current step is zero vector, and proceed to step 4.

[0040] S21: Based on the central difference assumption, calculate the position of the pantograph and rigid contact network system at the next moment under the non-contact state. The calculation formula is as follows:

[0041]

[0042] Among them, e i and e i+1 Let be the coordinate vectors of the pantograph-catenary system at the current and next time steps, respectively, and let i be the iteration time step, representing the state of the variable at time i. In this explanation, let i be the current time, i-1 be the previous time, i+1 be the next time, and Δt be the iteration time step. At the initial time, i.e., when i = 0, the coordinates of the system at the previous time step can be calculated using the following formula:

[0043]

[0044] Where e0 is the initial coordinate vector of the system, which is obtained by discretizing the rigid contact wire; The initial velocity of the system can be set as the zero vector; The initial acceleration vector of the system is calculated using the following formula:

[0045]

[0046] S22: Based on the position of the pantograph-catenary system at the next moment, determine the contact state between the pantograph and the catenary at the next moment. The determination rule is as follows:

[0047]

[0048] Where T = [0 1 0], and G is the contact constraint matrix, written as:

[0049] G = [0 6n×3 S c 0 6(N-n-2)×3 S p (20)

[0050] Define the pantograph's movement direction as forward and the opposite direction as backward; n is the number of all ANCF nodes behind the ANCF units in contact on the rigid contact network; N is the total number of ANCF nodes on the rigid contact network; S c S is the shape function matrix of the ANCF element currently in contact. p This is an auxiliary matrix used to calculate the vertical position of the bow head, written as:

[0051]

[0052] If the predicted state is "offline", then the contact force F at the current moment i =[0 0 0] T The coordinates of the pantograph-catenary system at the next moment are calculated using e from equation (16). i+1 If the predicted state is "contact", proceed to step 4; if the predicted state is "contact", proceed to step 3.

[0053] Step 3: Use Newton's iterative method to solve for the precise three-dimensional contact force at the current moment.

[0054] S31: Assemble the contact constraint matrix G at the current moment. i Normal contact force direction vector n n，i tangential friction force direction vector n i，i The direction vector of the resultant contact force n i The method for calculating the direction vector of the normal contact force is as follows:

[0055]

[0056] Among them, A cw and A p Let be the tangential vectors of the contact wire and the pantograph head, respectively. In the pantograph mass model, considering only the vertical vibration of the pantograph, the tangential vector of the pantograph head is written as:

[0057] A p =[0 0 1] T (twenty three)

[0058] The formula for calculating the tangential vector of the contact wire is:

[0059]

[0060] Among them, e cw Let x be the coordinate vector of the ANCF element in contact with the pantograph. cw These are the local coordinates of the contact point on this ANCF element.

[0061] The direction of the tangential friction force is the same as the direction of the relative velocity between the bow head and the contact line at the contact point; therefore, the formula for calculating the direction vector of the tangential friction force is:

[0062]

[0063] Where v is the relative velocity between the bow head and the contact line at the contact point, the calculation formula is:

[0064]

[0065] Among them, v p =[v t 0 0] T v t The influence of the vehicle's speed on the train's speed is demonstrated through precise calculation of the direction of the tangential friction force.

[0066] Using Coulomb's law of friction to define contact force, and assuming the coefficient of sliding friction is μ, the direction vector of the resultant contact force is:

[0067]

[0068] The expression for the total contact force is:

[0069] F = f·n (28)

[0070] Where f is the amplitude of the total contact force.

[0071] S32: Assemble the pantograph-catenary system considering the tangent stiffness matrix J and the generalized unbalanced load increment in the equation of motion under contact constraints. The calculation formula is as follows:

[0072]

[0073]

[0074] in

[0075]

[0076] S33: Calculate the generalized coordinate increment vector of the system and determine whether the iteration has converged; the formula for calculating the generalized coordinate increment vector is:

[0077]

[0078] Determine if the convergence condition max(ΔU) ≤ tolerance is met. If the convergence condition is met, output the contact force F at the current moment. i =f i ·n i And update the system coordinate vector e at the next time step. i+1 If the convergence condition is not met, then update f. i With e i+1 Then return to step S31.

[0079] Step 4: Determine if the pantograph-catenary contact force calculation for all time steps has been completed. If completed, output the pantograph-catenary contact force time series and end the process. If not completed, update the pantograph position according to the vehicle speed, update i = i + 1, and return to Step 2.

[0080] The beneficial technical effects of this invention compared to the prior art are as follows:

[0081] 1. The contact force calculation method proposed in this invention does not require providing empirical parameters similar to contact stiffness, nor does it have the problem of virtual penetration, and can obtain accurate contact force and pantograph-catenary system displacement results;

[0082] 2. In the process of calculating the contact force, this invention fully considers the three-dimensional change of the normal contact force caused by pantograph-catenary vibration; it also considers the tangential friction force, which can obtain an accurate three-dimensional contact force vector. At the same time, through the accurate calculation of the tangential friction force, the influence of train speed can be represented. Attached Figure Description

[0083] Figure 1 This is a schematic diagram of the contact and three-dimensional contact force between the bow and the fireproof wire in this invention.

[0084] Figure 2 This is a comparison between the Y-direction component of the bow-catenary contact force obtained in an embodiment of the present invention and the detection data.

[0085] Figure 3 The X-direction component of the bow-catenary contact force obtained in the embodiments of the present invention is shown.

[0086] Figure 4 The Z-direction component of the bow-catenary contact force obtained in the embodiments of the present invention is shown.

[0087] Figure 5 This is a comparison of the dynamic height of the contact line at the contact point and the bow head obtained in an embodiment of the present invention. Detailed Implementation

[0088] The present invention will now be described in further detail with reference to the accompanying drawings and specific embodiments.

[0089] This embodiment is a simulation case of dynamic contact between a straight section of rigid contact wire and a pantograph in one anchor segment. The vehicle speed is 60 km / h, and the contact force between the pantograph and the contact wire is calculated using the method proposed in this invention. Specific parameters in this embodiment are shown in Tables 1 and 2, and the coefficient of friction between the pantograph head and the contact wire is 0.3. In this embodiment, the train is assumed to move along the positive X direction, the Y direction is perpendicular, and the Z direction is perpendicular to the centerline of the track.

[0090] Table 1. Pantograph parameters in the embodiments

[0091]

[0092] Table 2 Parameters of Rigid Contact Network in Examples

[0093]

[0094] The specific implementation process of this embodiment is as follows:

[0095] Step 1: Based on the rigid contact network design parameters shown in Table 2, establish a finite element model of the rigid contact network using ANCF beam elements and three-dimensional linear rod elements. Based on the pantograph parameters, establish the pantograph model using the reduced mass block method, and construct the motion equations of the pantograph-rigid contact network.

[0096] S11: The rigid contact wire busbar and contact wire are discretized using ANCF beam elements. The coordinates of a single element are:

[0097] e A =[e1 e2 e3 e4 e5 e6 e7 e8 e9 e 10 e 11 e 12 ] T (1)

[0098] in,

[0099]

[0100] Where L is the original length of the element, x is a local variable on the element that varies from 0 to L, and r1, r2, and r3 are the coordinates of any point on the element in three-dimensional space, calculated by the following formula:

[0101] r = [r1 r2 r3] T =Se A (3)

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

[0103] S=[s1I3 s2I3 s3I3 s4I3] (4)

[0104] Where I3 is a third-order identity matrix, s i for:

[0105]

[0106] Where ξ = x / L.

[0107] Stiffness matrix K of ANCF element A and mass matrix M A They are respectively:

[0108]

[0109] K A =K L +K T (7)

[0110] in

[0111]

[0112]

[0113] In the formula, m A Let E be the mass of the element, E be the Young's modulus, A be the cross-sectional area of ​​the ANCF beam element, I be the moment of inertia of the ANCF beam element, and ε be the axial strain of the beam element.

[0114] S12: Employs a three-dimensional linear spring element discrete suspension structure with a mass matrix M S and stiffness matrix K S They are respectively:

[0115]

[0116]

[0117] Where, m S k is the equivalent mass of the suspension structure. x k y k z These are the stiffness coefficients of the spring element in the X, Y, and Z coordinate axes, respectively.

[0118] S13: Based on the element connection relationship, assemble the stiffness matrix of the ANCF beam element and the three-dimensional spring element to construct the stiffness matrix and mass matrix of the rigid contact network as a whole.

[0119] S14: Establish the reduced mass model of the pantograph, with the following equations:

[0120]

[0121] In the formula, m1, m2, and m3 are the masses of the bow head, upper frame, and lower frame, respectively; c1, c2, and c3 are the damping of the bow head, upper frame, and lower frame, respectively; k1, k2, and k3 are the stiffness of the bow head, upper frame, and lower frame, respectively; y1, y2, and y3 are the displacements of the bow head, upper frame, and lower frame, respectively; F C and F L These are the contact force between the pantograph and the catenary, and the static lifting force on the pantograph, respectively.

[0122] S15: The kinematic equations of the pantograph-rigid contact wire system are as follows:

[0123]

[0124] In the formula, Let M and C be the acceleration and velocity vectors of the system as a whole, respectively; let P be the external load vector acting on the system; and let Q be the elastic internal force vector of the system as a whole. The formula for calculating the elastic internal force of the ANCF beam element is:

[0125] Q A =K A e A (14)

[0126] The method for calculating the internal forces of a three-dimensional linear spring element is as follows:

[0127] Q S =K S U S (15)

[0128] Among them, U S Let be the displacement vector of the spring element.

[0129] Step 2: Use the central difference method to predict the motion state of the pantograph-rigid contact network system at the next moment. If the predicted state is "contact", proceed to step 3. If the predicted state is "offline", the contact force in the current step is zero vector, and proceed to step 4.

[0130] S21: Based on the central difference assumption, calculate the position of the pantograph and rigid contact network system at the next moment under the non-contact state. The calculation formula is as follows:

[0131]

[0132] Among them, e i and e i+1Let be the coordinate vectors of the pantograph-catenary system at the current and next time steps, respectively, and let i be the iteration time step, representing the state of the variable at time i. In this explanation, let i be the current time, i-1 be the previous time, i+1 be the next time, and Δt be the iteration time step. At the initial time, i.e., when i = 0, the coordinates of the system at the previous time step can be calculated using the following formula:

[0133]

[0134] Where e0 is the initial coordinate vector of the system, which is obtained by discretizing the rigid contact wire; The initial velocity of the system can be set as the zero vector; The initial acceleration vector of the system is calculated using the following formula:

[0135]

[0136] S22: Based on the position of the pantograph-catenary system at the next moment, determine the contact state between the pantograph and the catenary at the next moment. The determination rule is as follows:

[0137]

[0138] Where T = [0 1 0], and G is the contact constraint matrix, written as:

[0139] G = [θ] 6n×3 S c 0 6(N-n-2)×3 S p (20)

[0140] Define the pantograph's movement direction as forward and the opposite direction as backward; n is the number of all ANCF nodes behind the ANCF units in contact on the rigid contact network; N is the total number of ANCF nodes on the rigid contact network; S c S is the shape function matrix of the ANCF element currently in contact. p This is an auxiliary matrix used to calculate the vertical position of the bow head, written as:

[0141]

[0142] If the predicted state is "offline", then the contact force F at the current moment i =[0 0 0] T The coordinates of the pantograph-catenary system at the next moment are calculated using e from equation (16). i+1 If the predicted state is "contact", proceed to step 4; if the predicted state is "contact", proceed to step 3.

[0143] Step 3: Use Newton's iterative method to solve for the precise three-dimensional contact force at the current moment.

[0144] S31: Assemble the contact constraint matrix G at the current moment.i Normal contact force direction vector n n，i tangential friction force direction vector n 1，i The direction vector of the resultant contact force n i The method for calculating the direction vector of the normal contact force is as follows:

[0145]

[0146] Among them, A cw and A p Let be the tangential vectors of the contact wire and the pantograph head, respectively. In the pantograph mass model, considering only the vertical vibration of the pantograph, the tangential vector of the pantograph head is written as:

[0147] A p =[0 0 1] T (twenty three)

[0148] The formula for calculating the tangential vector of the contact wire is:

[0149]

[0150] Among them, e cw Let x be the coordinate vector of the ANCF element in contact with the pantograph. cw These are the local coordinates of the contact point on this ANCF element.

[0151] The direction of the tangential friction force is the same as the direction of the relative velocity between the bow head and the contact line at the contact point; therefore, the formula for calculating the direction vector of the tangential friction force is:

[0152]

[0153] Where v is the relative velocity between the bow head and the contact line at the contact point, the calculation formula is:

[0154]

[0155] Among them, v p =[v t 0 0] T v t The influence of the vehicle's speed on the train's speed is demonstrated through precise calculation of the direction of the tangential friction force.

[0156] Using Coulomb's law of friction to define the contact force, and assuming the coefficient of sliding friction is μ (0.3 in this embodiment), the direction vector of the resultant contact force is:

[0157]

[0158] The expression for the total contact force is:

[0159] F = f·n (28)

[0160] Where f is the amplitude of the total contact force. Figure 1 The diagram shows the contact force between the bow and the catenary and the three-dimensional contact force.

[0161] S32: Assemble the pantograph-catenary system considering the tangent stiffness matrix J and the generalized unbalanced load increment in the equation of motion under contact constraints. The calculation formula is as follows:

[0162]

[0163]

[0164] in

[0165]

[0166] S33: Calculate the generalized coordinate increment vector of the system and determine whether the iteration has converged; the formula for calculating the generalized coordinate increment vector is:

[0167]

[0168] Determine if the convergence condition max(ΔU) ≤ tolerance is met. If the convergence condition is met, output the contact force F at the current moment. i =f i ·n i And update the system coordinate vector e at the next time step. i+1 If the convergence condition is not met, then update f. i With e i+1 Then return to step S31.

[0169] Step 4: Determine if the pantograph-catenary contact force calculation for all time steps has been completed. If completed, output the pantograph-catenary contact force time series and end the process. If not completed, update the pantograph position according to the vehicle speed, update i = i + 1, and return to Step 2.

[0170] The pantograph-catenary contact force calculated through the above steps is as follows: Figures 2-4 As shown. Since actual line testing can only detect vertical contact force, the Y-direction component (i.e., vertical component) of the contact force obtained in the embodiment is compared with the detection data obtained from a line with the same parameters, as shown. Figure 2 As shown in Table 3, the statistical values ​​of contact force are compared.

[0171] Table 3 Comparison of the Y-direction component of the contact force obtained in the examples with the statistical values ​​of the detection data.

[0172]

[0173]

[0174] It can be seen that the fluctuation range of the contact force obtained in the embodiment is highly consistent with that of the time history curve of the detection data, and the deviation of the statistical index is also below 3%, proving that the contact force calculation method provided by the present invention can accurately calculate the contact force between the pantograph and the rigid contact wire. Furthermore, Figure 3 The X-direction component of the contact force obtained in the example is... Figure 4 The Z-component of the contact force obtained in the embodiment is shown. It can be seen that the method of the present invention can effectively represent the change of contact force in three-dimensional space, and can accurately calculate the tangential friction force component. Figure 5 A comparison of the dynamic heights of the pantograph head and the contact wire at the contact point obtained in the example shows that the height of the pantograph head and the height of the contact wire at the contact point are basically the same, and the pantograph head overlift (i.e., the pantograph head height minus the contact wire height) is within 10. -9 The m-order magnitude is due to calculation error, indicating that the method of the present invention can obtain accurate displacement results of the pantograph-catenary system and can effectively avoid the problem of virtual penetration at the contact position in the traditional penalty function method.

Claims

1. A method for accurate calculation of three-dimensional contact forces between a pantograph and a rigid catenary, characterized in that, Comprising the following steps: Step 1: According to the design parameters of rigid catenary, the ANCF beam element and three-dimensional linear bar element are used to establish the finite element model of rigid catenary, and according to the pantograph parameters, the lumped mass method is used to establish the pantograph model, and the pantograph-rigid catenary motion equation is constructed; S11: The ANCF beam element is used to discretize the rigid catenary busbar and contact wire, and the single unit coordinate is: （1） Wherein, （2） where L is the element length, is a local variable on the element varying from 0 to L, , , are the coordinates of any point on the element in three-dimensional space, calculated from the following equations: （3） For the shape function of the ANCF element, written as: （4） wherein is a three-order identity matrix, is: （5） wherein ; Stiffness matrix of an ANCF element and mass matrix are respectively: （6） （7） Wherein, （8） （9） wherein is the mass of the element, is the Young's modulus of elasticity, is the cross-sectional area of the ANCF beam element, is the cross-sectional moment of inertia of the ANCF beam element, is the axial strain of the beam element; S12: Discretize the suspension structure by three-dimensional linear spring elements, whose mass matrix and stiffness matrix are respectively （10） （11） wherein, is the equivalent mass of the suspension structure, , , are the stiffness coefficients of the spring unit in the XYZ three coordinate axis directions, respectively; S13: According to the connection relationship of the unit, the stiffness matrix of the ANCF beam element and the three-dimensional spring element is assembled, and the stiffness matrix and the mass matrix of the whole rigid catenary are constructed; S14: The lumped mass model of the pantograph is established, and its equation is as follows: （12） wherein, , , M, M and M are the masses of the head, upper frame and lower frame, respectively; , , C, C and C are the dampings of the head, upper frame and lower frame, respectively; , , K, K and K are the stiffnesses of the head, upper frame and lower frame, respectively; , , X, X and X are the displacements of the head, upper frame and lower frame, respectively; and F and F are the pantograph-catenary contact force and the static lifting force on the pantograph, respectively. S15: The pantograph-rigid catenary system kinematic equation is constructed as follows: （13） where, , are the acceleration and velocity vectors of the system as a whole, respectively, , are the mass and damping matrices of the system as a whole, respectively; is the external load vector on the system; is the internal elastic force vector of the system as a whole; the ANCF beam element elastic internal force calculation formula is: （14） The three-dimensional linear spring element internal force calculation method is: （15） wherein is the spring unit displacement vector; Step 2: The central difference method is used to predict the next time pantograph-rigid catenary system motion state; S21: Based on the central difference assumption, the position of the pantograph and the rigid catenary system at the next time under the condition of no contact is calculated, and the calculation formula is: （16） wherein, and are the coordinate vectors of the pantograph-catenary system at the current time and the next time, respectively, i is the iteration time step number, denotes the state of the variable at the i time, in this explanation, it is assumed that i is the current time, is the last time, is the next time, is the iteration time step, at the initial time, i.e. , the system last time coordinate can be calculated by the following formula: （17） wherein, is the initial coordinate vector of the system, obtained from the discretization of the rigid catenary; is the initial velocity of the system, which can be set as a zero vector; is the initial acceleration vector of the system, the calculation formula of which is: （18） S22: According to the next time position of the pantograph-catenary system, the contact state between the pantograph and the catenary at the next time is judged, and the judgment rule is: （19） wherein , is the contact constraint matrix, written as: （20） The pantograph motion direction is defined as forward and the opposite direction as backward, NANCF is the number of all ANCF nodes behind the ANCF element in contact with the rigid catenary, NANCF is the number of all ANCF nodes behind the ANCF element in contact with the rigid catenary, NANCF is the number of all ANCF nodes behind the ANCF element in contact with the rigid catenary, NANCF is the number of all ANCF nodes behind the ANCF element in contact with the rigid catenary, （21） If the predicted state is "off-line", the contact force at the current time instant is calculated using the formula (16) and the process goes to step 4; if the predicted state is "contact", the process goes to step 3. If the predicted state is "off-line", the contact force at the current time instant Step 3: The Newton iteration method is used to solve the accurate three-dimensional contact force at the current time; Step 4: Judge whether the pantograph-catenary contact force calculation of all time steps is completed, if it is completed, output the pantograph-catenary contact force sequence and end; If not completed, update the pantograph position according to the vehicle running speed, update the iteration time step index, and return to step 2.

2. The method for accurate calculation of three-dimensional contact forces between a pantograph and a rigid catenary according to claim 1, characterized in that, The step 3 is specifically: S31: Assemble the contact constraint matrix of the current time , normal contact force direction vector , tangent friction force direction vector , contact force resultant direction vector The method for calculating the normal contact force direction vector is: （22） where, and are the tangential vectors of the contact line and the pantograph head, respectively. In the pantograph mass model, only the pantograph vertical vibration is considered, and the tangential vector of the pantograph head is written as: （23） The tangential vector calculation formula of the contact wire is: （24） wherein, is the coordinate vector of the ANCF element in contact with the pantograph, is the local coordinate of the contact point on this ANCF element; The tangential friction force direction is consistent with the relative velocity direction between the pantograph and the contact wire at the contact point; Therefore, the tangential friction force direction vector calculation formula is: （25） wherein is the relative speed between the arch and the contact line at the contact point, calculated as （26） wherein , is the train running speed, the influence of the vehicle running speed is shown by the accurate calculation of the direction of the tangential friction force; The contact force is defined by Coulomb's law of friction, assuming a constant coefficient of sliding friction The direction vector of the contact force resultant is （27） The total contact force expression is: （28） wherein is the magnitude of the total contact force; S32: Assemble the pantograph-catenary system considering the tangent stiffness matrix of the motion equation under contact constraints with generalized unbalanced load increments The calculation formula is: （29） （30） Wherein, （31） S33: Calculate the system generalized coordinate increment vector, and judge whether the iteration is converged; The generalized coordinate increment vector calculation formula is: （32） judging the convergence condition , if the convergence condition is satisfied, outputting the contact force at the current time , and updating the system coordinate vector at the next time ; if the convergence condition is not satisfied, updating and , and returning to step S31.

Images