Parallel acceleration method for axle coupling calculation

Abstract

The invention discloses a parallel acceleration method for axle coupling calculation, and belongs to the field of vehicle bridge coupling vibration and scientific calculation acceleration, and the method comprises the steps: employing a power explicit integration method, and solving the displacement and speed of a train, a steel rail, a sleeper and a railway ballast at the current time step in a multi-thread manner at the same time; according to the position of the train on the bridge under the current time step, wheel-rail contact points of all the carriages of the train are calculated at the same time in a multi-thread mode through a trace method, and wheel-rail interaction force of all the carriages of the train is calculated according to the wheel-rail contact points of all the carriages of the train; and based on the train subsystem kinematic equation and the track subsystem kinematic equation, solving the acceleration of each carriage of the train and the acceleration of the track, and based on the bridge subsystem kinematic equation, solving the displacement, the speed and the acceleration of each contact point of the bridge. According to the invention, the characteristic of multiple calculation cores of a modern computer processor is fully utilized, the utilization rate of the computer processor is improved, and the traditional axle coupling calculation process is efficiently accelerated.

Description

Technical Field

[0001] The invention belongs to the field of vehicle-bridge coupled vibration and scientific computing acceleration, and in particular relates to a parallel acceleration method for vehicle-bridge coupled computing. Background Art

[0002] As railway bridges develop towards larger spans and higher flexibility, the bridges deform significantly under train loads, causing the bridge track alignment to change at the same time, thereby changing the wheel-rail contact state and affecting the train's driving stability. Therefore, given that the bridge, train, and track structure on the bridge form a highly coupled dynamic system, it is necessary to establish a refined dynamic model for each subsystem and carry out vehicle-bridge coupling solution.

[0003] The solution of vehicle-bridge coupled dynamic system often involves multiple subsystems, multiple time steps, large matrix iterations, etc., especially for long-span bridges. The ultra-high structural freedom makes the solution of the bridge system have higher requirements on computer performance or takes longer computing time. As modern computer processors have multi-computing core characteristics, the traditional vehicle-bridge coupling serial solution has low computing performance utilization, and cannot give full play to the advantages of multi-computing cores. The overall utilization of computer processors is generally 20%. Therefore, it is particularly important to improve the utilization of computer processors, make full use of multi-core and multi-threaded characteristics, and accelerate vehicle-bridge coupling calculations. Summary of the invention

[0004] In view of the above-mentioned deficiencies in the prior art, the present invention provides a parallel acceleration method for vehicle-bridge coupling calculation, which fully utilizes the multi-computing core characteristics of modern computer processors, improves the utilization rate of computer processors, and efficiently accelerates the traditional vehicle-bridge coupling calculation process.

[0005] In order to achieve the above-mentioned invention object, the technical solution adopted by the present invention is: a parallel acceleration method for vehicle-bridge coupling calculation, comprising the following steps:

[0006] S1. Construct the kinematic equations of the bridge subsystem, the track subsystem and the train subsystem, and assign initial values ​​to the speed, displacement and acceleration of the bridge, track and train;

[0007] S2, calculate the position of the train on the bridge at the current time step;

[0008] S3, based on the speed, displacement and acceleration of the track and train in the previous time step, the dynamic explicit integration method is used to solve the displacement and speed of the train and the displacement and speed of the track in the current time step simultaneously through multiple threads;

[0009] S4, extracting the estimated displacement value and estimated speed value of the bridge at the wheel set position of the vehicle according to the position of the train on the bridge, the displacement and speed of the train, and the displacement and speed of the track at the current time step;

[0010] S5, according to the position of the train on the bridge at the current time step and the estimated displacement value of the bridge at the position of the vehicle wheel pair, the track method is used to simultaneously calculate the wheel-rail contact points of each carriage of the train in multiple threads, and the wheel-rail interaction force of each carriage of the train is calculated according to the wheel-rail contact points of each carriage of the train;

[0011] S6. Calculate the bridge-rail interaction force based on the estimated displacement value and velocity value of the bridge at the wheel set position of the vehicle and the displacement and velocity of the track;

[0012] S7, substituting the wheel-rail interaction force of each train car into the train subsystem kinematic equation, solving the acceleration of each train car, and solving the acceleration of the track based on the track subsystem kinematic equation according to the bridge-rail interaction force and the wheel-rail interaction force, and solving the displacement, velocity and acceleration of each contact point of the bridge based on the bridge subsystem kinematic equation according to the estimated value of the bridge velocity at the vehicle wheel position by using the implicit integration method;

[0013] S8, judging whether the train displacement, bridge displacement and track displacement are all converged, if so, outputting the dynamic response of the train, the dynamic response of the bridge and the dynamic response of the track, and directly proceeding to step S9, otherwise, updating the initial value according to the current dynamic response of the train, the dynamic response of the bridge and the dynamic response of the track, and returning to step S2;

[0014] S9. Determine whether the train is within the calculation range. If so, end the process. Otherwise, return to step S2.

[0015] Furthermore, the step S2 is specifically as follows:

[0016] S201, calculating the running distance of the train according to the speed, displacement and acceleration of the train in the previous time step, and corresponding to the specific position on the bridge based on the running distance of the train, that is, the estimated position of the train;

[0017] S202, obtaining the track irregularity and bridge deformation at the estimated position of the train, and adding the track irregularity and bridge deformation to obtain the position of the train on the bridge at the current time step:

[0018] D Bi =X Bi +r i

[0019] r i =S v (Ω)+G(n 0 )

[0020]

[0021] Among them, D Bi is the position of the train on the bridge at the current time step; X Bi is the bridge deformation at the estimated location; r i The track irregularity at the estimated position; S v (Ω) is the power spectrum density function of track height irregularity; G(n 0 ) is the road surface spectrum value; A v is the roughness constant; Ω c , Ω and Ω r are cutoff frequencies.

[0022] Furthermore, the expression of the dynamic explicit integration method in step S3 is:

[0023]

[0024] Among them, X n+1 is the displacement of n+1 time step; X n is the displacement of n time steps; is the velocity of n time steps; Δt is the time step length; ψ and All are integral parameters; is the acceleration of n time steps; is the acceleration of the n-1 time step; is the velocity at time step n+1.

[0025] Furthermore, the expression of the wheel-rail contact point of each carriage of the train in step S5 is:

[0026]

[0027] l x =-cosφ w sinψ w

[0028] l y =cosφ w cosψ w

[0029] l z = sinφ w

[0030]

[0031] Among them, x c is the abscissa of the wheel-rail contact point; O is the center of the wheel rolling circle 2 The horizontal axis of x , ly and l z All are intermediate variables; R w is the radius of the wheel rolling circle; tg is the tangent function; δ w is the wheel tread contact angle; y c is the ordinate of the wheel-rail contact point; O is the center of the wheel rolling circle 2 The vertical coordinate; y w is the lateral displacement of the wheelset; z c is the vertical coordinate of the wheel-rail contact point; O is the center of the wheel rolling circle 2 The vertical coordinate of w is the rolling angle of the wheelset; ψ w is the wheel set's yaw angle; d w It is the transverse coordinate of each rolling circle of the wheel tread in the wheelset coordinate system.

[0032] Furthermore, the wheel-rail interaction force in step S5 includes the wheel-rail contact normal force and creep force.

[0033] Furthermore, the expression of the wheel-rail contact normal force is:

[0034]

[0035] δN L =δZ L / cos(δ wL +φ w )

[0036] δN R =δZ R / cos(δ wR -φ w )

[0037] Among them, P N (t) is the normal force of the wheel-rail contact at the current time step; t is time; G is the wheel-rail contact constant; R is the radius of the wheel rolling circle; δN(t) is the relative normal compression of the wheel-rail contact point; δN L is the relative normal compression between the left wheel and rail; δZ L is the relative vertical displacement of the wheel / rail at the left contact point; δ wL is the contact angle of the left wheel tread; φ w is the rolling angle of the wheelset; δN R is the relative normal compression between the right wheel and rail; δZ R is the relative vertical displacement of the wheel / rail at the right contact point; δ wR is the contact angle of the right wheel tread.

[0038] Furthermore, the creep force is expressed as:

[0039]

[0040] G(ab)=2G w G r / (G w +G r )

[0041]

[0042] Among them, F x ′ is the longitudinal creep force after correction using Shen’s theory; ε is the correction coefficient; F x is the initial longitudinal creep force; F y ′ is the lateral creep force after correction using Shen’s theory; F y is the initial lateral creep force; M z ′ is the rotation creep torque after correction using Shen’s theory; M z is the initial rotation creep torque; F R ′ is the resultant force of the lateral and longitudinal creep forces after nonlinear correction; F R is the resultant force of the lateral and longitudinal creep forces; f is the friction coefficient between the wheel and rail; N is the total normal force of the wheel-rail contact; f 11 is the longitudinal creep coefficient; ξ x is the longitudinal creep rate; f 22 is the lateral creep coefficient; ξ y is the lateral creep rate; f 23 is the rotational lateral creep coefficient; ξ sp is the rotation creep rate; f 33 is the rotation creep coefficient; G(ab) is the composite shear modulus of the wheelset and track material; a is the major semi-axis of the contact ellipse; b is the minor semi-axis of the contact ellipse; C 11 , C 22 , C 23 and C 33 All are Kalker coefficients; G w is the shear modulus of the wheelset material, G r is the shear modulus of the track material; V w1 The contact ellipse on the wheel tread is e 1 Axis movement speed; V r1 is the contact ellipse e on the rail 1 The movement speed corresponding to the axis; V w2 The contact ellipse on the wheel tread is e 2 Axis movement speed; V r2 is the contact ellipse e on the rail 2 The movement speed corresponding to the axis; Ω w3 is the contact ellipse around the wheel tread 3Axis movement speed; Ω r3 is the contact ellipse e on the rail 3 The corresponding movement speed of the axis; V is the nominal forward speed of the wheelset on the rail.

[0043] Furthermore, the step S6 is specifically as follows:

[0044] S601. Obtain the second-order ordinary differential equation of rail vibration:

[0045]

[0046] in, is the second-order derivative of the vertical rail canonical coordinate; t is time; E r is the Young's modulus of the rail; I ry is the moment of inertia of the rail section around the y-axis; ρ r is the mass density of the rail; A r is the cross-sectional area of ​​the rail; k is the cut-off modal order number; π is the circumference of a circle; l is the length of the rail; q zk (t) is the vertical rail canonical coordinate; N s N is the number of sleepers within the calculated length of the rail; w The number of axles within the calculated length of the rail; F sVi is the vertical reaction force of the ith support point of the rail; Z k () is the vibration mode function of the rail in the vertical direction; P V j is the vertical force acting on the rail by the jth wheel; N Z is the cutoff mode order of the vertical vibration mode of the rail; is the second-order derivative of the transverse rail canonical coordinate; I z is the moment of inertia of the rail section around the z axis; q yk (t) is the transverse rail canonical coordinate; F sHi is the lateral reaction force of the ith support point of the rail; Y k () is the vibration mode function of the rail in the transverse direction; P H j is the lateral force exerted by the jth wheel on the rail; N Y is the cutoff mode order in the transverse direction of the rail; is the second-order derivative of the canonical coordinates of the torsion rail; G r is the shear modulus of the rail; J rt is the torsional inertia moment of the rail section; J r0 is the polar moment of inertia of the rail section; q tk (t) is the canonical coordinate of the twisted rail; F sTi is the torsional reaction force of the ith support point of the rail; Φ k () is the vibration mode function of rail torsion; P T j is the torque of the jth wheel acting on the rail; x siis the x coordinate of the ith support point of the rail; Pj is the x-coordinate of the j-th wheelset of the rail; N T is the cutoff mode order of the rail torsional vibration mode;

[0047] S602. According to the second-order ordinary differential equation of rail vibration, the actual vibration displacement, velocity and acceleration of the rail at the current moment are obtained:

[0048]

[0049]

[0050] Among them, z r (x,t),y r (x,t),φ r (x, t) are the vertical, lateral and torsional displacements of the actual vibration of the rail at the x-coordinate at time t; are respectively the vertical, lateral and torsional velocities of the actual vibration of the rail at the x-coordinate at time t; is the vertical, lateral and torsional acceleration of the actual vibration of the rail at the x-coordinate at time t;

[0051] S603, obtaining the bridge-rail force according to the wheel-rail force, the actual vibration displacement, velocity and acceleration of the rail at the current moment.

[0052] Furthermore, the step S7 is specifically as follows:

[0053] S701, substituting the wheel-rail interaction force of each train car into the kinematic equation of the train subsystem, and solving the acceleration of each train car simultaneously by multi-threading;

[0054] S702, according to the bridge-track interaction force, use the interpolation method to form the equivalent load of the bridge node;

[0055] S703, substituting the equivalent load of the bridge node and the wheel-rail interaction force of each train carriage into the kinematic equation of the track subsystem to solve the acceleration of the track;

[0056] S704. Based on the estimated value of the bridge speed at the wheel set position of the vehicle, the equivalent load of the bridge node and the kinematic equation of the bridge subsystem, the implicit integration method is used to solve the displacement, velocity and acceleration of each contact point of the bridge simultaneously in multiple threads.

[0057] The beneficial effects of the present invention are as follows: the present invention establishes a refined dynamic model for trains, bridges and track structures and constructs a dynamic solution process method. The present invention fully considers the independence of the calculation conditions between the subsystems and the spatial independence between the vehicles during the calculation process, constructs a parallel acceleration method suitable for vehicle-bridge coupling calculation, fully utilizes the multi-core characteristics of modern computer processors, improves the utilization rate of computer processors, accelerates the vehicle-bridge coupling calculation process, and is suitable for more modern intelligent transportation scenarios with high frequency and high efficiency of data. BRIEF DESCRIPTION OF THE DRAWINGS

[0058] Figure 1 The present invention is a flow chart of the method.

[0059] Figure 2 It is a schematic diagram of the wheel-rail contact geometric relationship in a specific application scenario of the present invention.

[0060] Figure 3 It is a schematic diagram of the right wheel-rail contact point in a specific application scenario of the present invention.

[0061] Figure 4 The figure is a parallel calculation flow chart of the vehicle-bridge coupling calculation program according to an embodiment of the present invention. DETAILED DESCRIPTION

[0062] The specific implementation modes of the present invention are described below so that those skilled in the art can understand the present invention. However, it should be clear that the present invention is not limited to the scope of the specific implementation modes. For those of ordinary skill in the art, as long as various changes are within the spirit and scope of the present invention as defined and determined by the attached claims, these changes are obvious, and all inventions and creations utilizing the concept of the present invention are protected.

[0063] like Figure 1 and Figure 4 As shown, in one embodiment of the present invention, a parallel acceleration method for vehicle-bridge coupling calculation includes the following steps:

[0064] S1. Construct the kinematic equations of the bridge subsystem, the track subsystem and the train subsystem, and assign initial values ​​to the speed, displacement and acceleration of the bridge, track and train;

[0065] S2, calculate the position of the train on the bridge at the current time step;

[0066] S3, based on the speed, displacement and acceleration of the track and train in the previous time step, the dynamic explicit integration method is used to solve the displacement and speed of the train and the displacement and speed of the track in the current time step simultaneously through multiple threads;

[0067] S4, extracting the estimated displacement value and estimated speed value of the bridge at the wheel set position of the vehicle according to the position of the train on the bridge, the displacement and speed of the train, and the displacement and speed of the track at the current time step;

[0068] S5, according to the position of the train on the bridge at the current time step and the estimated displacement value of the bridge at the position of the vehicle wheel pair, the track method is used to simultaneously calculate the wheel-rail contact points of each carriage of the train in multiple threads, and the wheel-rail interaction force of each carriage of the train is calculated according to the wheel-rail contact points of each carriage of the train;

[0069] S6. Calculate the bridge-rail interaction force based on the estimated displacement value and velocity value of the bridge at the wheel set position of the vehicle and the displacement and velocity of the track;

[0070] S7, substituting the wheel-rail interaction force of each train car into the train subsystem kinematic equation, solving the acceleration of each train car, and solving the acceleration of the track based on the track subsystem kinematic equation according to the bridge-rail interaction force and the wheel-rail interaction force, and solving the displacement, velocity and acceleration of each contact point of the bridge based on the bridge subsystem kinematic equation according to the estimated value of the bridge velocity at the vehicle wheel position by using the implicit integration method;

[0071] S8, judging whether the train displacement, bridge displacement and track displacement are all converged, if so, outputting the dynamic response of the train, the dynamic response of the bridge and the dynamic response of the track, and directly proceeding to step S9, otherwise, updating the initial value according to the current dynamic response of the train, the dynamic response of the bridge and the dynamic response of the track, and returning to step S2;

[0072] S9. Determine whether the train is within the calculation range. If so, end the process. Otherwise, return to step S2.

[0073] In this embodiment, a suspension bridge with a main span of one kilometer is used as the engineering background, and the specific method of establishing the vehicle-bridge coupling serial program is as follows: the main cable and the hanger are simulated using a 6-degree-of-freedom spatial rod unit, a bridge finite element model is established, and the mass, stiffness and damping matrix of the bridge structure are obtained. The rails are modeled using an infinitely long elastic Euler beam. The stiffness of the sleeper itself is much greater than that of the ballast and fasteners. In the model, it is considered as a rigid body to consider the lateral and vertical translation of the sleeper and the rotation around the line direction. The ballast is a piled bulk body, which is loaded as a cone. The train vehicle adopts a multi-rigid body dynamics model.

[0074] Identify the mass, stiffness and damping characteristics of each subsystem of the bridge, railway train and railway line and establish the kinematic equations of each subsystem. The specific subsystems include the train subsystem, track subsystem and bridge subsystem. The kinematic equations of each subsystem are as follows:

[0075] Kinematic equations of the train subsystem:

[0076] Kinematic equations of the track subsystem:

[0077] Kinematic equations of the bridge subsystem:

[0078] Where M V , C V , K V are the mass matrix, damping matrix, and stiffness matrix of the train subsystem respectively; M T , C T , K T are the mass matrix, damping matrix, and stiffness matrix of the track subsystem respectively; M B , C B , K B are the mass matrix, damping matrix, and stiffness matrix of the bridge subsystem respectively; X V , are the displacement, velocity and acceleration vectors of the train subsystem respectively; X T , are the displacement, velocity, and acceleration vectors of the track subsystem respectively; X B , are the displacement, velocity and acceleration vectors of the bridge subsystem respectively; F V-T 、F T-V is the interaction force between the train and the track; F B-T 、F T-B is the interaction force between the track and the bridge.

[0079] The step S2 is specifically as follows:

[0080] S201, calculating the running distance of the train according to the speed, displacement and acceleration of the train in the previous time step, and corresponding to the specific position on the bridge based on the running distance of the train, that is, the estimated position of the train;

[0081] S202, obtaining the track irregularity and bridge deformation at the estimated position of the train, and adding the track irregularity and bridge deformation to obtain the position of the train on the bridge at the current time step:

[0082] D Bi =X Bi +r i

[0083] r i =S v (Ω)+G(n 0 )

[0084]

[0085] Among them, D Biis the position of the train on the bridge at the current time step; X Bi is the bridge deformation at the estimated location; r i The track irregularity at the estimated position; S v (Ω) is the power spectrum density function of track height irregularity; G(n 0 ) is the road surface spectrum value; A v is the roughness constant; Ω c , Ω and Ω r are cutoff frequencies.

[0086] In this embodiment, firstly, the running distance of the vehicle at each time step is calculated according to the running speed of the train and the specific position on the bridge corresponding to the distance is determined; the track unevenness or road unevenness at the position is read in and compared with the bridge deformation X at the position. Bi Add them together to get the actual spatial position of the point.

[0087] The expression of the dynamic explicit integration method in step S3 is:

[0088]

[0089] Among them, X n+1 is the displacement of n+1 time step; X n is the displacement of n time steps; is the velocity of n time steps; Δt is the time step length; ψ and All are integral parameters; is the acceleration of n time steps; is the acceleration of the n-1 time step; is the velocity at time step n+1.

[0090] In this embodiment, the vehicle, rails, sleepers and ballast are all solved by the dynamic explicit integration method. In each solution process, only the system displacement, velocity and acceleration at the previous moment are used as input. There is no coupling relationship between the systems. Therefore, 4 independent threads can be allocated and the OmpSection instruction in OpenMP is used to realize parallel computing of the four subsystems.

[0091] In this embodiment, the dynamic display integral of the train includes three parts: the train body, the bogie and the wheels, all of which are obtained through the display integration method, thereby extracting the displacement estimation value and speed estimation value of the bridge at the position of the vehicle wheelset. The bridge displacement estimation value directly affects and The center of the wheel rolling circle is O 2 The horizontal and vertical coordinates of the bridge are shown in Figure 5. The estimated velocity value is not used in S5 and is mainly used for the subsequent implicit integration of the bridge's own displacement, velocity, and acceleration.

[0092] The expression of the wheel-rail contact point of each carriage of the train in step S5 is:

[0093]

[0094] l x =-cosφ w sinψ w

[0095] l y =cosφ w cosψ w

[0096] l z = sinφ w

[0097]

[0098] Among them, x c is the abscissa of the wheel-rail contact point; O is the center of the wheel rolling circle 2 The horizontal axis of x , l y and l z All are intermediate variables; R w is the radius of the wheel rolling circle; tg is the tangent function; δ w is the wheel tread contact angle; y c is the ordinate of the wheel-rail contact point; O is the center of the wheel rolling circle 2 The vertical coordinate; y w is the lateral displacement of the wheelset; z c is the vertical coordinate of the wheel-rail contact point; O is the center of the wheel rolling circle 2 The vertical coordinate of w is the rolling angle of the wheelset; ψ w is the wheel set's yaw angle; d w It is the transverse coordinate of each rolling circle of the wheel tread in the wheelset coordinate system.

[0099] In this embodiment, for the train-track interaction, the actual bridge position D at the train travel position is obtained. Bi According to the contact geometric relationship between the train wheel tread and the rail head tread, the wheel-rail contact point C is found by the trace method. The geometric relationship of the wheel-rail contact point C is as follows: Figure 2 As shown, the wheel-rail contact point C is within three planes, including O 2 The rolling circle plane of the point, with O 2 The radius of the wheel rolling circle is R w is a sphere with radius O 1 -O 1'-C'-C plane; the coordinates of the wheel-rail contact point C in the absolute coordinate system can be derived through the equations of the three planes (x c ,y c , z c ).

[0100] When the wheelset lateral displacement y w 、Head shaking angle ψ w 、Roll angle φ w When a certain time is reached, by gradually changing the horizontal coordinates d of each rolling circle of the wheel tread in the wheelset coordinate system w , the wheel-rail spatial contact trace at a certain moment can be constructed; the calculation of the wheel-rail contact point is affected by the movement displacement of the rail pair and the rail and the change of the track irregularity value, and the wheel-rail contact point needs to be recalculated at every moment.

[0101] The wheel-rail interaction force in step S5 includes the wheel-rail contact normal force and the creep force.

[0102] In this embodiment, the calculation of the wheel-rail interaction force includes two parts: searching for the wheel-rail spatial contact point using the trace method and calculating the wheel-rail contact normal force and creep force. The wheel-rail contact force between the track and the train only depends on the relative spatial position of the wheelset and the track at each moment. Each car is independent of each other, so the parallel method is to allocate an independent thread for the calculation of each car; this embodiment includes 16 cars, so 16 threads are allocated.

[0103] The expression of the wheel-rail contact normal force is:

[0104]

[0105] δN L =δZ L / cos(δ wL +φ w )

[0106] δN R =δZ R / cos(δ wR -φ w )

[0107] Among them, P N (t) is the normal force of the wheel-rail contact at the current time step; t is time; G is the wheel-rail contact constant; R is the radius of the wheel rolling circle; δN(t) is the relative normal compression of the wheel-rail contact point; δN L is the relative normal compression between the left wheel and rail; δZ L is the relative vertical displacement of the wheel / rail at the left contact point; δ wL is the contact angle of the left wheel tread; φ w is the rolling angle of the wheelset; δN Ris the relative normal compression between the right wheel and rail; δZ R is the relative vertical displacement of the wheel / rail at the right contact point; δ wR is the contact angle of the right wheel tread.

[0108] The expression of creep force is:

[0109]

[0110]

[0111] G(ab)=2G w G r / (G w +G r )

[0112]

[0113] Among them, F x ′ is the longitudinal creep force after correction using Shen’s theory; ε is the correction coefficient; F x is the initial longitudinal creep force; F y ′ is the lateral creep force after correction using Shen’s theory; F y is the initial lateral creep force; M z ′ is the rotation creep torque after correction using Shen’s theory; M z is the initial rotation creep torque; F R ′ is the resultant force of the lateral and longitudinal creep forces after nonlinear correction; F R is the resultant force of the lateral and longitudinal creep forces; f is the friction coefficient between the wheel and rail; N is the total normal force of the wheel-rail contact; f 11 is the longitudinal creep coefficient; ξ x is the longitudinal creep rate; f 22 is the lateral creep coefficient; ξ y is the lateral creep rate; f 23 is the rotational lateral creep coefficient; ξ sp is the rotation creep rate; f 33 is the rotation creep coefficient; G(ab) is the composite shear modulus of the wheelset and track material; a is the major semi-axis of the contact ellipse; b is the minor semi-axis of the contact ellipse; C 11 , C 22 , C 23 and C 33 All are Kalker coefficients; G w is the shear modulus of the wheelset material, G r is the shear modulus of the track material; V w1 The contact ellipse on the wheel tread is e 1 Axis movement speed; V r1 is the contact ellipse e on the rail1 The movement speed corresponding to the axis; V w2 The contact ellipse on the wheel tread is e 2 Axis movement speed; V r2 is the contact ellipse e on the rail 2 The movement speed corresponding to the axis; Ω w3 is the contact ellipse around the wheel tread 3 Axis movement speed; Ω r3 is the contact ellipse e on the rail 3 The corresponding movement speed of the axis; V is the nominal forward speed of the wheelset on the rail.

[0114] In this embodiment, according to the Hertz nonlinear elastic contact theory, the wheel-rail contact normal force is obtained by the wheel-rail relative compression; according to the lateral displacement, yaw angle, roll angle, lateral displacement, vertical displacement, torsion angle and track irregularity of the wheelset at a certain moment, the vertical relative displacement between the left and right wheels and rails at that moment is obtained according to the calculation method of the wheel-rail contact point mentioned above, and then the vertical relative displacement between the left and right wheels and rails at zero moment is subtracted; the vertical relative displacement is projected to the normal of the wheel-rail contact point to obtain the normal compression between the left and right wheels and rails, that is:

[0115] δN L =δZ L / cos(δ wL +φ w )

[0116] δN R =δZ R / cos(δ wR -φ w )

[0117] In the formula, δ wL , δ wR ,φ w are the contact angle of the left and right wheel treads and the wheel pair roll angle respectively; δZ L , δZ R are the vertical relative displacements of the wheel and rail at the left and right contact points respectively; δN L , δN R are the relative normal compression between the left and right wheel rails, respectively;

[0118] Using Hertz nonlinear elastic contact theory, the vertical force P between the wheel and rail is determined N (t), vertical force P N (t) is the normal force of wheel-rail contact:

[0119]

[0120] Where R is the radius of the wheel rolling circle; δN(t) is the normal elastic compression of the wheel-rail contact point; G is the wheel-rail contact constant, for a conical tread: G = 4.57R -0.149 ×10 -8 (m / N 2 / 3 ), for wear-type tread: G = 3.86R -0.115 ×10 -8 (m / N 2 / 3 ).

[0121] Regarding the calculation of wheel-rail creep force, Kalker's linear creep theory is first used for calculation. Considering that the creep between the wheel and rail may reach saturation during vehicle operation, the Shen-Hedrick-Elkins theory is used for nonlinear correction.

[0122] With the center C of the wheel-rail contact ellipse as the origin, the wheel-rail contact spot coordinate system Ce is established 1 e 2 e 3 ,like Figure 3 As shown, the right wheel-rail contact point is taken as an example; e 1 The axis coincides with the X axis, which is the forward direction of the wheel; 2 The axis is located in the wheel-rail contact plane, and the angle between it and the Y axis is δ at the left contact point. L +φ, at the right wheel-rail contact point is δ R -φ;e 3 The axis is normal to the wheel-rail contact plane;

[0123] The step S6 is specifically as follows:

[0124] S601. Obtain the second-order ordinary differential equation of rail vibration:

[0125]

[0126] in, is the second-order derivative of the vertical rail canonical coordinate; t is time; E r is the Young's modulus of the rail; I ry is the moment of inertia of the rail section around the y-axis; ρ r is the mass density of the rail; A r is the cross-sectional area of ​​the rail; k is the cut-off modal order number; π is the circumference of a circle; l is the length of the rail; q zk (t) is the vertical rail canonical coordinate; N s N is the number of sleepers within the calculated length of the rail; w The number of axles within the calculated length of the rail; F sVi is the vertical reaction force of the ith support point of the rail; Z k () is the vibration mode function of the rail in the vertical direction; PVj is the vertical force exerted by the jth wheel on the rail; N Z is the cutoff mode order of the vertical vibration mode of the rail; is the second-order derivative of the transverse rail canonical coordinate; I z is the moment of inertia of the rail section around the z axis; q yk (t) is the transverse rail canonical coordinate; F sHi is the lateral reaction force of the ith support point of the rail; Y k () is the vibration mode function of the rail in the transverse direction; P H j is the lateral force exerted by the jth wheel on the rail; N Y is the cutoff mode order in the transverse direction of the rail; is the second-order derivative of the canonical coordinates of the torsion rail; G r is the shear modulus of the rail; J rt is the torsional inertia moment of the rail section; J r0 is the polar moment of inertia of the rail section; q tk (t) is the canonical coordinate of the twisted rail; F sTi is the torsional reaction force of the ith support point of the rail; Φ k () is the vibration mode function of rail torsion; P T j is the torque of the jth wheel acting on the rail; x si is the x coordinate of the ith support point of the rail; Pj is the x-coordinate of the j-th wheelset of the rail; N T is the cutoff mode order of the rail torsional vibration mode;

[0127] S602. According to the second-order ordinary differential equation of rail vibration, the actual vibration displacement, velocity and acceleration of the rail at the current moment are obtained:

[0128]

[0129] Among them, z r (x,t),y r (x,t),φ r (x, t) are respectively the vertical, lateral and torsional displacements of the actual vibration of the rail at the x-coordinate at time t; are respectively the vertical, lateral and torsional velocities of the actual vibration of the rail at the x-coordinate at time t; is the vertical, lateral and torsional acceleration of the actual vibration of the rail at the x-coordinate at time t;

[0130] S603, obtaining the bridge-rail force according to the wheel-rail force, the actual vibration displacement, velocity and acceleration of the rail at the current moment.

[0131] The step S7 is specifically as follows:

[0132] S701, substituting the wheel-rail interaction force of each train car into the kinematic equation of the train subsystem, and solving the acceleration of each train car simultaneously by multi-threading;

[0133] S702, according to the bridge-track interaction force, use the interpolation method to form the equivalent load of the bridge node;

[0134] S703, substituting the equivalent load of the bridge node and the wheel-rail interaction force of each train carriage into the kinematic equation of the track subsystem to solve the acceleration of the track;

[0135] S704. Based on the estimated value of the bridge speed at the wheel set position of the vehicle, the equivalent load of the bridge node and the kinematic equation of the bridge subsystem, the implicit integration method is used to solve the displacement, velocity and acceleration of each contact point of the bridge simultaneously in multiple threads.

[0136] In this embodiment, the interaction force and acceleration of each component in the train multi-body dynamics model are calculated by solving the vehicle motion equation established based on the D'Alembert principle. The interaction force and acceleration of each component of each car depend on the contact force transmitted from the rail to the wheelset, so each car is independent of each other, and the parallel method is to allocate an independent thread for the calculation of each car.

[0137] In this embodiment, the Newmark-β numerical integration method of the bridge structure motion equation includes the number of loops for the number of degrees of freedom of the bridge structure, and a dynamic thread allocation method is used for parallelization of multiple independent loops not involving iterations.

[0138] In this embodiment, the rails are subjected to the wheel-rail contact force transmitted by the train wheelset, and are gradually transmitted downward to the fasteners, sleepers and ballast bed, and finally to the bridge;

[0139] Since the bridge model is established using the finite element method, the position of the sleeper position point generally does not coincide with the bridge node. Therefore, the interpolation method is used to form the equivalent load of the bridge node. Considering that the bridge will deform under external loads and be constrained by the supports, the following formula is used for interpolation fitting:

[0140] Y b =N 1 (x)u yi +N 2 (x)θ zi +N 3 (x)u yj +N 4 (x)θ zj

[0141] Z b =N 1 (x)u zi -N 2 (x)θ yi+N 3 (x)u zj -N 4 (x)θ yj

[0142] φ b =N 5 (x)θ xi +N 6 (x)θ xj

[0143]

[0144] Where Y b is the lateral displacement of the centroid of the bridge unit under the sleeper; Z b is the vertical displacement of the centroid of the bridge unit under the sleeper; φ b is the torsional displacement of the bridge unit centroid under the sleeper; N 1 (x), N 2 (x), N 3 (x), N 4 (x), N 5 (x), N 6 (x) is the shape function of the spatial beam element; L is the element length; x is the longitudinal coordinate of the sleeper position in the bridge element coordinate system; u yi 、u yj 、u zi 、u zi ,θ zi ,θ zj ,θ yi ,θ yj ,θ xi ,θ xj are the lateral, vertical and torsional displacements of the left and right nodes of each unit, and i and j are the components of the two nodes of each unit;

[0145] The equivalent bridge node load F is obtained by multiplying the shape function of the spatial beam element and the sleeper load. T-B ;

[0146] F yi =N 1 (x)F gy

[0147] F zi =N 1 (x)F gz

[0148] F yj =N 3 (x)F gy

[0149] F zj =N3 (x)F gz

[0150] M xi =N 5 (x)M gx

[0151] M yi =-N 2 (x)F gz

[0152] M zi =N 2 (x)F gy

[0153] M xj =N 6 (x)M gx

[0154] M yj =-N 4 (x)F gz

[0155] M zj =N 4 (x)F gy

[0156] In the formula, F yi 、F zi 、F yj 、F zj 、M xi 、M yi 、M zi 、M xj 、M yj 、M zj The forces F transmitted to the bridge by the track are T-B Components in the three directions of lateral, vertical and torsional directions and at the two nodes i and j of each element;

[0157] Bridge equivalent node load F T-B Substitute the right-hand side of the bridge subsystem motion equation into the equation; calculate the displacement, velocity and acceleration of the bridge subsystem using the Newmark-β method;

[0158] The displacement, velocity and acceleration corresponding to the bridge subsystem are solved:

[0159] The Newmark-β method integration format for solving the motion equations for the bridge subsystem is:

[0160]

[0161] Where α and β are integration parameters, which can generally be taken as 0.25 and 0.5 respectively; Δt is the time integration step; n represents the nth integration step;

[0162] The displacement of the bridge subsystem X B ,speed Acceleration vector As X n , Substituting this into the equation, we can solve the equation of motion of the bridge subsystem;

[0163] In the above formula Utilizing X n+1 To represent, then respectively Substitute into the bridge subsystem dynamics equation at time n+1 and simplify to obtain X n+1 The solution formula is:

[0164]

[0165] In the formula, F n+1 is the bridge node load at time n+1; M, C, and K are the mass matrix, damping matrix, and stiffness matrix, respectively;

[0166] X n+1 The solution formula is:

[0167]

[0168] In the formula, a 6 =Δt(1-β), a 7 =βΔt, solving the above equation yields X n+1 ; α is the integral parameter, which is 0.25; β is the integral parameter, which is 0.5;

[0169] The solution formula is as follows:

[0170]

[0171] The obtained train wheel-rail contact force F T-V Substitute the right-hand side of the kinematic equation of the train subsystem into the train wheel-rail contact force F T-V , bridge equivalent node load F T-B Substitute the right-hand side of the kinematic equation of the track subsystem into the equation; calculate the displacement, velocity and acceleration corresponding to the track subsystem and train subsystem through the fast display integration method, and obtain the dynamic response of each subsystem at time n+1.

[0172] The wheel-rail contact force between the track and the train is spatially parallelized for each car in each time step;

[0173] The interaction forces and accelerations of the components in the train multi-body dynamics model are calculated spatially and in parallel for each car in each time step;

[0174] The interaction force and acceleration of each component in the train multi-body dynamics model are calculated by solving the vehicle motion equation based on the D'Alembert principle. The interaction force and acceleration of each component of each car depend on the contact force transmitted from the rail to the wheelset, so each car is independent of each other. The parallel method of the OmpSection directive in OpenMP is used to allocate independent threads for the calculation of each car. This embodiment includes 16 cars, so 16 threads are allocated;

[0175] The multiple loops of the Newmark-β method for the bridge structure motion equations are parallelized using a dynamic thread allocation method.

[0176] For bridge structures, there are many node units and the number of degrees of freedom is often high. When solving the Newmark-β numerical integration method, the number of loops corresponding to the order of freedom is involved. Therefore, the OmpDo instruction in OpenMP can be used to dynamically open threads for each loop, so that idle threads can efficiently execute loops.

[0177] Establish a parallel computing structure that takes into account the time-space iterative characteristics of vehicle-bridge coupling. The vehicle-bridge coupling calculation is an iterative process with multiple time steps. The integral step size is selected as 0.0001s in this embodiment, and the time for the train to cross the bridge is more than 20s, so it involves huge time steps. Opening and closing parallel domains in each time step will waste a lot of time and computing power. Therefore, in this example, before performing a multi-time step loop, a 16-thread parallel domain is constructed, and the OpenMP Single instruction is used for the serial code, and the main thread is used for execution. When encountering a wheel-rail force solution for a cycle such as multiple trains, enter the parallel space, re-enter the serial code when contact is made, and use the main thread for execution. The parallel structure of the vehicle-bridge coupling calculation is as follows: Figure 2 shown.

[0178] This method is based on the OpenMP method. Aiming at the independence of the calculation of each part in the vehicle-bridge coupling calculation, a parallel structure process is proposed, which greatly improves the computer processor utilization rate during the calculation and improves the vehicle-bridge coupling calculation speed. In this embodiment, the serial code execution time is about 120 minutes, and the parallel code execution time is 54 minutes, which greatly reduces the calculation time.

Claims

1. A parallel acceleration method for vehicle-bridge coupling calculation, characterized in that: The following steps are involved: S1. Construct the kinematic equations of the bridge subsystem, the track subsystem and the train subsystem, and assign initial values ​​to the speed, displacement and acceleration of the bridge, track and train; S2, calculate the position of the train on the bridge at the current time step; S3, based on the speed, displacement and acceleration of the track and train in the previous time step, the dynamic explicit integration method is used to solve the displacement and speed of the train and the displacement and speed of the track in the current time step simultaneously through multiple threads; S4, extracting the estimated displacement value and estimated speed value of the bridge at the wheel set position of the vehicle according to the position of the train on the bridge, the displacement and speed of the train, and the displacement and speed of the track at the current time step; S5, according to the position of the train on the bridge at the current time step and the estimated displacement value of the bridge at the position of the vehicle wheel pair, the track method is used to simultaneously calculate the wheel-rail contact points of each carriage of the train in multiple threads, and the wheel-rail interaction force of each carriage of the train is calculated according to the wheel-rail contact points of each carriage of the train; S6. Calculate the bridge-rail interaction force based on the estimated displacement value and velocity value of the bridge at the wheel set position of the vehicle and the displacement and velocity of the track; S7, substituting the wheel-rail interaction force of each train car into the train subsystem kinematic equation, solving the acceleration of each train car, and solving the acceleration of the track based on the track subsystem kinematic equation according to the bridge-rail interaction force and the wheel-rail interaction force, and solving the displacement, velocity and acceleration of each contact point of the bridge based on the bridge subsystem kinematic equation according to the estimated value of the bridge velocity at the vehicle wheel position by using the implicit integration method; S8, judging whether the train displacement, bridge displacement and track displacement are all converged, if so, outputting the dynamic response of the train, the dynamic response of the bridge and the dynamic response of the track, and directly proceeding to step S9, otherwise, updating the initial value according to the current dynamic response of the train, the dynamic response of the bridge and the dynamic response of the track, and returning to step S2; S9. Determine whether the train is within the calculation range. If so, end the process. Otherwise, return to step S2.

2. The parallel acceleration method for vehicle-bridge coupling calculation according to claim 1, characterized in that: The step S2 is specifically as follows: S201, calculating the running distance of the train according to the speed, displacement and acceleration of the train in the previous time step, and corresponding to the specific position on the bridge based on the running distance of the train, that is, the estimated position of the train; S202, obtaining the track irregularity and bridge deformation at the estimated position of the train, and adding the track irregularity and bridge deformation to obtain the position of the train on the bridge at the current time step: D Bi =X Bi +r i r i =S v (Ω)+G(n0) Among them, D Bi is the position of the train on the bridge at the current time step; X Bi is the bridge deformation at the estimated location; r i The track irregularity at the estimated position; S v (Ω) is the power spectrum density function of track height irregularity; G(n0) is the road surface spectrum value; A v is the roughness constant; Ω c , Ω and Ω r are cutoff frequencies.

3. The parallel acceleration method for vehicle-bridge coupling calculation according to claim 1, characterized in that: The expression of the dynamic explicit integration method in step S3 is: Among them, X n+1 is the displacement of n+1 time step; X n is the displacement of n time steps; is the velocity of n time steps; Δt is the time step length; ψ and All are integral parameters; is the acceleration of n time steps; is the acceleration of the n-1 time step; is the velocity at time step n+1.

4. The parallel acceleration method for vehicle-bridge coupling calculation according to claim 1, characterized in that: The expression of the wheel-rail contact point of each carriage of the train in step S5 is: l x =-cosφ w sinψ w l y =cosφ w soψ w the z =sinφ w Among them, x c is the abscissa of the wheel-rail contact point; is the horizontal coordinate of the center O2 of the wheel rolling circle; l x , l y and l z All are intermediate variables; R w is the radius of the wheel rolling circle; tg is the tangent function; δ w is the wheel tread contact angle; y c is the ordinate of the wheel-rail contact point; y is the ordinate of the center O2 of the wheel rolling circle; w is the lateral displacement of the wheelset; z c is the vertical coordinate of the wheel-rail contact point; is the vertical coordinate of the center O2 of the wheel rolling circle; φ w is the rolling angle of the wheelset; ψ w is the wheel set's yaw angle; d w It is the transverse coordinate of each rolling circle of the wheel tread in the wheelset coordinate system.

5. The parallel acceleration method for vehicle-bridge coupling calculation according to claim 1, characterized in that: The wheel-rail interaction force in step S5 includes the wheel-rail contact normal force and the creep force.

6. The parallel acceleration method for vehicle-bridge coupling calculation according to claim 5, characterized in that: The expression of the wheel-rail contact normal force is: δN L =δZ L / cos(δ wL +φ w ) δN R =δZ R / cos(δ wR -f w ) Among them, P N (t) is the normal force of the wheel-rail contact at the current time step; t is time; G is the wheel-rail contact constant; R is the radius of the wheel rolling circle; δN(t) is the relative normal compression of the wheel-rail contact point; δN L is the relative normal compression between the left wheel and rail; δZ L is the relative vertical displacement of the wheel / rail at the left contact point; δ wL is the contact angle of the left wheel tread; φ w is the rolling angle of the wheelset; δN R is the relative normal compression between the right wheel and rail; δZ R is the relative vertical displacement of the wheel / rail at the right contact point; δ wR is the contact angle of the right wheel tread.

7. The parallel acceleration method for vehicle-bridge coupling calculation according to claim 5, characterized in that: The expression of creep force is: G(ab)=2G w G r / (G w +G r ) Among them, F x ′ is the longitudinal creep force after correction using Shen’s theory; ε is the correction coefficient; F x is the initial longitudinal creep force; F y ′ is the lateral creep force after correction using Shen’s theory; F y is the initial lateral creep force; M z ′ is the rotation creep torque after correction using Shen’s theory; M z is the initial rotation creep torque; F R ′ is the resultant force of the lateral and longitudinal creep forces after nonlinear correction; F R is the resultant force of the lateral and longitudinal creep forces; f is the friction coefficient between the wheel and rail; N is the total normal force of the wheel-rail contact; f 11 is the longitudinal creep coefficient; ξ x is the longitudinal creep rate; f 22 is the lateral creep coefficient; ξ y is the lateral creep rate; f 23 is the rotational lateral creep coefficient; ξ sp is the rotation creep rate; f 33 is the rotation creep coefficient; G(ab) is the composite shear modulus of the wheelset and track material; a is the major semi-axis of the contact ellipse; b is the minor semi-axis of the contact ellipse; C 11 , C 22 , C 23 and C 33 All are Kalker coefficients; G w is the shear modulus of the wheelset material, G r is the shear modulus of the track material; V w1 V is the movement speed of the contact ellipse on the wheel tread along the e1 axis; r1 V is the motion speed corresponding to the axis of the contact ellipse e1 on the rail; w2 V is the movement speed of the contact ellipse on the wheel tread along the e2 axis; r2 is the motion speed corresponding to the contact ellipse e2 axis on the rail; Ω w3 is the speed of the contact ellipse on the wheel tread around the e3 axis; Ω r3 is the motion speed corresponding to the axis of the contact ellipse e3 on the rail; V is the nominal forward speed of the wheelset on the rail.

8. The parallel acceleration method for vehicle-bridge coupling calculation according to claim 1, characterized in that: The step S6 is specifically as follows: S601. Obtain the second-order ordinary differential equation of rail vibration: in, is the second-order derivative of the vertical rail canonical coordinate; t is time; E r is the Young's modulus of the rail; I ry is the moment of inertia of the rail section around the y-axis; ρ r is the mass density of the rail; A r is the cross-sectional area of ​​the rail; k is the cut-off modal order number; π is the circumference of a circle; l is the length of the rail; q zk (t) is the vertical rail canonical coordinate; N s N is the number of sleepers within the calculated length of the rail; w The number of axles within the calculated length of the rail; F sVi is the vertical reaction force of the ith support point of the rail; Z k () is the vibration mode function of the rail in the vertical direction; P V j is the vertical force acting on the rail by the jth wheel; N Z is the cutoff mode order of the vertical vibration mode of the rail; is the second-order derivative of the transverse rail canonical coordinate; I z is the moment of inertia of the rail section around the z axis; q yk (t) is the transverse rail canonical coordinate; F sHi is the lateral reaction force of the ith support point of the rail; Y k () is the vibration mode function of the rail in the transverse direction; P H j is the lateral force exerted by the jth wheel on the rail; N Y is the cutoff mode order of the rail in the transverse direction; is the second-order derivative of the canonical coordinates of the torsion rail; G r is the shear modulus of the rail; J rt is the torsional inertia moment of the rail section; J r0 is the polar moment of inertia of the rail section; q tk (t) is the canonical coordinate of the twisted rail; F sTi is the torsional reaction force of the ith support point of the rail; Φ k () is the vibration mode function of rail torsion; P T j is the torque of the jth wheel acting on the rail; x si is the x coordinate of the ith support point of the rail; Pj is the x-coordinate of the j-th wheelset of the rail; N T is the cutoff mode order of the rail torsional vibration mode; S602. According to the second-order ordinary differential equation of rail vibration, the actual vibration displacement, velocity and acceleration of the rail at the current moment are obtained: Among them, z r (x,t),y r (x,t),φ r (x, t) are respectively the vertical, lateral and torsional displacements of the actual vibration of the rail at the x-coordinate at time t; are respectively the vertical, lateral and torsional velocities of the actual vibration of the rail at the x-coordinate at time t; is the vertical, lateral and torsional acceleration of the actual vibration of the rail at the x-coordinate at time t; S603, obtaining the bridge-rail force according to the wheel-rail force, the actual vibration displacement, velocity and acceleration of the rail at the current moment.

9. The parallel acceleration method for vehicle-bridge coupling calculation according to claim 1, characterized in that: The step S7 is specifically as follows: S701, substituting the wheel-rail interaction force of each train car into the kinematic equation of the train subsystem, and solving the acceleration of each train car simultaneously by multi-threading; S702, according to the bridge-track interaction force, use the interpolation method to form the equivalent load of the bridge node; S703, substituting the equivalent load of the bridge node and the wheel-rail interaction force of each train carriage into the kinematic equation of the track subsystem to solve the acceleration of the track; S704. Based on the estimated value of the bridge speed at the wheel set position of the vehicle, the equivalent load of the bridge node and the kinematic equation of the bridge subsystem, the implicit integration method is used to solve the displacement, velocity and acceleration of each contact point of the bridge simultaneously in multiple threads.

Images