Quasi-steady-state numerical simulation method for evolution characteristic of icing form of railway overhead line system
By using a quasi-steady-state numerical simulation method to study the evolution characteristics of icing morphology in railway overhead contact lines, the shortcomings of existing technologies in simulating the movement of freezing rain droplets and the icing phase transition process have been addressed, enabling accurate prediction and efficient calculation of icing morphology in overhead contact lines.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- CENT SOUTH UNIV
- Filing Date
- 2026-01-12
- Publication Date
- 2026-04-24
AI Technical Summary
Existing technologies are insufficient to accurately predict the evolution of icing patterns on railway overhead contact lines, especially under complex wind and rain coupling conditions. Traditional experimental testing and numerical simulation methods cannot fully consider the motion characteristics of freezing rain droplets and the icing phase transition process, making it difficult to balance computational efficiency and accuracy.
A quasi-steady-state numerical simulation method for the evolution characteristics of icing morphology in railway catenary is adopted. This method involves numerical simulation modeling, injecting freezing rain droplets using the Lagrange method, handling droplet collisions using the SSD standard breakage model, simulating the icing process using an ice accumulation model, and using the RBF interpolation function to achieve dynamic updates of the ice layer shape. The calculation is optimized by combining a dual-timescale strategy.
It realizes the simulation of the entire process of icing morphology of overhead contact lines, improves computational efficiency and accuracy, can simulate freezing rain processes under different meteorological conditions, reduces computational costs, and meets the needs of actual engineering.
Smart Images

Figure CN121920276A_ABST
Abstract
Description
Technical Field
[0001] This invention relates to the field of railway electrification power supply systems, and in particular to a quasi-steady-state numerical simulation method for the evolution characteristics of icing morphology in railway overhead contact lines. Background Technology
[0002] The railway overhead contact system is a crucial component of the high-speed railway electrification power supply system, and its safe and stable operation directly affects the reliability of train power supply. Under extreme weather conditions such as freezing rain, the surface of the contact wires is highly susceptible to icing, leading to increased contact resistance and a deterioration of the current collection relationship between the pantograph and the contact wire. Severe icing can cause pantograph arcing, arc interruption, and even power outages, seriously threatening train operation safety. Current research on the problem of contact wire icing mainly focuses on the following two categories:
[0003] 1) Experimental testing method: An icing experiment was conducted on a small section of the contact wire under a simulated freezing rain environment to observe the thickness and morphology of the ice. However, such experiments are limited by site and equipment conditions, making it difficult to reproduce complex wind and rain coupling conditions, and it is impossible to obtain data on the entire process of ice layer evolution over time.
[0004] 2) Numerical simulation methods: Currently, computational fluid dynamics (CFD) methods are used to simulate the impact and icing of freezing rain droplets. However, most existing simulations only consider simple unidirectional fluid-structure interaction, neglecting the breakup, adhesion, and liquid film phase transition evolution processes after droplet impact. Therefore, the prediction of ice layer shape is not accurate enough and cannot meet the needs of practical engineering. Moreover, the real-time changes in ice shape require a large amount of computational resources, making it difficult to balance computational efficiency and accuracy.
[0005] In summary, there is currently a lack of a comprehensive numerical simulation method that can simultaneously reflect the motion characteristics of freezing rain droplets, the icing phase transition process, and update the ice layer shape in real time. Therefore, it is urgent to propose new technical solutions to overcome the above shortcomings in order to accurately predict the evolution of icing patterns on railway overhead contact lines. Summary of the Invention
[0006] The purpose of this invention is to address the shortcomings of the aforementioned background technology by providing a numerical simulation method that fully considers the impact of freezing rain droplets, liquid film phase change, and the evolution of icing shape, so as to accurately predict the evolution law of icing morphology of railway catenary.
[0007] To achieve the above objectives, this invention provides a quasi-steady-state numerical simulation method for the evolution characteristics of icing morphology in railway overhead contact lines, comprising the following steps:
[0008] S1, numerical simulation modeling, establish air flow field, define temperature, humidity, and wind speed boundary conditions, and simulate the flow field around the contact wire under different meteorological conditions;
[0009] S2, In the numerical simulation model, freezing rain droplets are injected using the Lagrange method, and the parameters of the freezing rain droplets, including initial velocity, diameter distribution, temperature and specific heat capacity, are set to reproduce different freezing rain intensities and environmental characteristics;
[0010] S3, the trajectory of the discrete phase droplets in the air flow field is constantly evolving. When the droplets move to the vicinity of the contact wire surface, it is calculated whether they will collide with the contact wire surface based on the air flow field velocity and the droplet inertia.
[0011] S4. After impact, some droplets adhere to the surface of the contact network, converge to form a thin liquid film and spread on the surface. As time goes by and the surface temperature decreases, the water in the liquid film gradually freezes into ice. An ice accumulation model is used to simulate the freezing process of the liquid film. The ice accumulation model calculates the parameters of each attached droplet or liquid film unit based on the thin film energy equation, and calculates the latent heat exchange and extracts the freezing parameters to construct the growth mode of contact network icing.
[0012] S5 uses the RBF interpolation function to smooth the mesh, realizing the dynamic evolution of the ice layer shape and obtaining numerical simulation results.
[0013] Furthermore, when establishing the airflow field in S1, the mass conservation equation, the momentum equation of the continuous phase, the energy conservation equation, and the water vapor transport equation are followed.
[0014] Furthermore, in S2, a Lagrange framework is used to track the discrete phase of freezing rain droplets, following the particle motion equation and the particle thermal balance equation.
[0015] Furthermore, in S2, the bidirectional coupling source terms of the discrete phase and the continuous phase are calculated and set.
[0016] Furthermore, the collision process between the droplets and the contact wire surface in S3 is handled using the SSD standard breakage model. The SSD standard breakage model is based on dimensionless parameters, including the Weber number, Ohnesorge number, and incident angle. The conditions for breakup and the diameter distribution, velocity distribution and number of droplets after breakup are given.
[0017] Furthermore, the thin film thickness evolution in S4 is expressed as follows:
[0018]
[0019] in, For film thickness, The tangential divergence along the surface, The wall velocity is... For the mass flux of foreign attached droplets. To freeze the consumption of mass flux, This refers to the mass flux lost due to slippage.
[0020] Furthermore, the thin film energy equation is expressed as:
[0021]
[0022] in, The density of liquid water, Specific heat capacity of the liquid The convective heat transfer coefficient is... For the film temperature, For the thermal conductivity of liquid, For free flow temperature, This is the latent heat of phase transition.
[0023] Furthermore, when using the ice accumulation model, the ice-liquid interface satisfies the Stefan condition:
[0024]
[0025] in, The density of ice, For ice thickness, The thermal conductivity of ice.
[0026] Furthermore, in S5, nodes on the contact wire surface and its adjacent areas are selected. A normal displacement is assigned to each surface node based on the thickness of the newly added ice layer. Then, the RBF interpolation function is used to calculate the position adjustment of each node within the entire computational domain, achieving continuous tracking of the ice growth process. The formula for calculating the node displacement within the domain using the RBF interpolation function is:
[0027]
[0028] in, The coordinates of the target node. Let be the displacement vector of the target node. The number of all boundary nodes involved in the interpolation. These are the RBF interpolation coefficients. For the selected RBF kernel function, The coordinates of the boundary nodes, Let be the Euclidean distance from the target node to the boundary node. These are linear polynomial terms.
[0029] Furthermore, the numerical simulation process is decomposed into two time scales: small-step fluid-solid-thermal coupling propulsion and large-step ice layer shape update. Small-step propulsion is used to advance the rapid changes in air flow field, droplet motion, and liquid film icing, while larger time steps are used for macroscopic changes in ice layer shape. After accumulating the preset time interval, a mesh deformation is performed in one concentrated operation.
[0030] The above-described solution of the present invention has the following beneficial effects:
[0031] The quasi-steady-state numerical simulation method for the evolution characteristics of icing morphology of railway catenary provided by this invention fully considers the entire process from microscopic droplet impact to macroscopic ice layer formation on the catenary surface, including the entire process of droplet impact, breakage, liquid film formation, phase change and ice layer morphology growth. Compared with the prior art, it avoids the prediction deviation caused by simplification of a single link. At the same time, it realizes the updating of the ice layer shape during the simulation process by dynamically deforming the mesh through the RBF interpolation function, thus overcoming the limitations of the traditional fixed boundary method.
[0032] The quasi-steady-state numerical simulation method for the evolution characteristics of icing morphology of railway catenary provided by this invention allows users to customize boundary conditions such as temperature, humidity, and wind speed to simulate the flow field around the catenary under different meteorological conditions. When injecting freezing rain droplets using the Lagrange method, users can set parameters such as initial velocity, diameter distribution, temperature, and specific heat capacity to reproduce different freezing rain intensities and environmental characteristics.
[0033] The quasi-steady-state numerical simulation method for the evolution characteristics of icing morphology of railway catenary provided by this invention adopts a dual-time-scale strategy of quasi-steady-state evolution and grid deformation scheduling, which significantly reduces the computational cost of long-term series simulation. While ensuring the computational accuracy of each stage, it reduces unnecessary geometric update frequency, making it possible to simulate the actual freezing rain icing process that lasts for several hours. The computational efficiency is significantly improved compared with single small-step coupled calculation.
[0034] Other beneficial effects of the present invention will be described in detail in the following detailed description section. Attached Figure Description
[0035] Figure 1 This is a flowchart of the steps of the present invention. Detailed Implementation
[0036] The following specific examples illustrate the implementation of this disclosure. Those skilled in the art can easily understand other advantages and effects of this disclosure from the content disclosed in this specification. Obviously, the described embodiments are only a part of the embodiments of this disclosure, and not all of them. This disclosure can also be implemented or applied through other different specific embodiments, and the details in this specification can also be modified or changed based on different viewpoints and applications without departing from the spirit of this disclosure. It should be noted that, in the absence of conflict, the following embodiments and features in the embodiments can be combined with each other. Based on the embodiments in this disclosure, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this disclosure.
[0037] It should be noted that various aspects of embodiments within the scope of the appended claims are described below. It will be apparent that the aspects described herein can be embodied in a wide variety of forms, and any particular structure and / or function described herein is merely illustrative. Based on this disclosure, those skilled in the art will understand that one aspect described herein can be implemented independently of any other aspect, and two or more of these aspects can be combined in various ways. For example, any number of aspects set forth herein can be used to implement the device and / or practice the method. Additionally, this device and / or method can be implemented using structures and / or functionalities other than one or more of the aspects set forth herein.
[0038] It should also be noted that the illustrations provided in the following embodiments are merely schematic representations of the basic concept of this disclosure. The illustrations only show components relevant to this disclosure and are not drawn according to the actual number, shape, and size of components in implementation. In actual implementation, the type, quantity, and proportion of each component can be arbitrarily changed, and the component layout may be more complex. Furthermore, specific details are provided in the following description to facilitate a thorough understanding of the examples. However, those skilled in the art will understand that the described aspects can be practiced without these specific details.
[0039] like Figure 1 As shown, an embodiment of the present invention provides a quasi-steady-state numerical simulation method for the evolution characteristics of icing morphology in railway overhead contact lines, comprising the following steps:
[0040] S1. (In a suitable numerical simulation software) Construct a numerical simulation model, establish the air flow field, and define boundary conditions such as temperature, humidity, and wind speed to simulate the flow field around the contact wire under different meteorological conditions.
[0041] It should be noted that airflow fields typically include continuous phases such as air and water vapor. Therefore, when establishing an airflow field, the mass conservation equation must be followed:
[0042]
[0043] in, It is a velocity vector. For quality source items and , which represents the rate of gas mass generation or consumption within a local volume.
[0044] The momentum equation for a continuous phase is:
[0045]
[0046] in, air density, For air pressure, Aerodynamic viscosity, For turbulent viscosity, It is the gravitational acceleration vector. This is the momentum source term applied to the discrete relative to the continuous phase.
[0047] The energy conservation equation is:
[0048]
[0049] in, For isobaric specific heat capacity, For temperature, Thermal conductivity, The term represents the energy source, indicating the heat exchange effect of droplets releasing / absorbing latent heat.
[0050] The equation for water vapor transport is:
[0051]
[0052] in, It represents the mass fraction of water vapor. For the effective diffusion coefficient, This refers to the mass source term generated by droplet evaporation / condensation.
[0053] S2. In the numerical simulation model, freezing rain droplets are injected using the Lagrange method, and parameters such as the initial velocity, diameter distribution, temperature, and specific heat capacity of the freezing rain droplets are set to reproduce different freezing rain intensities and environmental characteristics.
[0054] In this embodiment, the freezing rain droplets are set as a discrete phase, and the discrete phase of the freezing rain droplets is tracked using a Lagrange framework. The particle motion equation is:
[0055]
[0056]
[0057] in, For the particle position, For particle velocity, For particle mass, The density of liquid water, The air resistance to the droplet. To add mass force, For other forces.
[0058] At the same time, the particle thermal balance equation must be followed:
[0059]
[0060] in, For particle temperature, Specific heat capacity of the liquid The convective heat transfer coefficient is... For particle surface area, The mass of phase transition per unit time. This is the latent heat of phase transition.
[0061] Considering the bidirectional coupling source terms of the discrete and continuous phases, further calculations are performed. , , and as follows:
[0062]
[0063]
[0064]
[0065]
[0066] in, This represents the volume of a single mesh cell in the numerical simulation model. Therefore, S1 and S2 complete the construction and initial setup of the numerical simulation model.
[0067] S3, the discrete phase droplets move in the air flow field, and their trajectory evolves continuously under the action of forces such as gravity and air resistance. When the droplets move to the vicinity of the contact wire surface, it is calculated whether they will collide with the contact wire surface based on the air flow field velocity and the droplet inertia.
[0068] In this embodiment, the collision process between the droplet and the contact wire surface is handled using the SSD standard breakage model. The SSD standard breakage model is based on dimensionless parameters, including the Weber number (…). The ratio of the inertial force acting on a droplet to its surface tension is a measure of the force exerted on the droplet. The Ohnesorge number (…) (Describing the coupling of viscosity and surface tension) and the angle of incidence. When the Weber number of droplet impacts exceeds a critical value, breakup is determined to occur. It is assumed that the broken droplets follow a Rosin-Rammler distribution to describe the complete size spectrum. Based on the law of conservation of mass, the number of generated droplets is inferred from the total mass of the parent droplet and the size distribution of the daughter droplets. In terms of velocity distribution, the daughter droplets mainly inherit the tangential momentum of the parent droplet, while the normal velocity is distributed according to a random distribution or statistical model after considering energy dissipation. This gives the conditions for breakup and the diameter distribution, velocity distribution, and number of the broken droplets. , The initial droplet diameter, The impact normal velocity, is the surface tension coefficient of the droplet.
[0069] Furthermore, custom breakup criteria can be introduced into numerical simulations, combining experimental data or field observations to improve the accuracy of secondary droplet distribution after droplet breakup. When customizing these criteria, it is possible to use... When the composite function exceeds the experimentally calibrated threshold, secondary fragmentation or splashing behavior is actively triggered. After triggering, the standard subdroplet distribution function and bounce / slip model can be overridden, and the correlation between the secondary droplet size distribution, ejection angle and velocity obtained by directly fitting experimental data or field observations can be injected, thereby achieving a high degree of consistency between simulation results and real physical phenomena.
[0070] S4. After impact, some droplets adhere to the surface of the contact network, converge to form a thin liquid film and spread on the surface. As time goes by and the surface temperature decreases, the water in the liquid film gradually freezes into ice. An ice accumulation model is used to simulate the freezing process of the liquid film.
[0071] In this embodiment, the icing model calculates key parameters such as the freezing rate of each attached droplet or liquid film unit based on the thin film energy equation, ensuring a balance between the latent heat released by the droplets and the heat absorbed by the surrounding environment. The icing model can also calculate latent heat exchange and extract icing parameters (ice thickness, density, phase change rate, etc.), thereby constructing a growth pattern for contact wire icing. When using the icing model, the thin film thickness evolution can be expressed as:
[0072]
[0073] in, For film thickness, The tangential divergence along the surface, The wall velocity is... For the mass flux of foreign attached droplets. To freeze the consumption of mass flux, This refers to the mass flux lost due to slippage.
[0074] The thin film energy equation can be expressed as:
[0075]
[0076] in, For the film temperature, For the thermal conductivity of liquid, This refers to the free-flow temperature.
[0077] In addition, when using the ice accumulation model, the ice-liquid interface also needs to satisfy the Stefan condition (one-dimensional normal approximation), which can be expressed as:
[0078]
[0079] in, The density of ice, For ice thickness, The thermal conductivity of ice, It is a spatial coordinate, specifically referring to the direction perpendicular to the ice-liquid interface. This represents the temperature gradient along the normal direction of the interface.
[0080] S5 uses the RBF interpolation function (radial basis function) to smooth the mesh, realizing the dynamic evolution of the ice layer shape and obtaining numerical simulation results.
[0081] In this embodiment, nodes on the contact wire surface and its adjacent areas are selected. A normal displacement is assigned to each surface node based on the thickness of the newly added ice layer. Then, the position adjustment of each node within the entire computational domain is calculated using the RBF interpolation function. Therefore, it can update the shape of the ice layer on the contact wire surface while maintaining mesh quality, achieving continuous tracking of the ice growth process. The formula for calculating the node displacement within the domain using the RBF interpolation function is:
[0082]
[0083] in, The coordinates of the target node. Let be the displacement vector of the target node. The number of all boundary nodes involved in the interpolation. These are the RBF interpolation coefficients. For the selected RBF kernel function, The coordinates of the boundary nodes, Let be the Euclidean distance from the target node to the boundary node. These are linear polynomial terms.
[0084] It should be noted that by moving nodes instead of adding nodes, the growth process of the ice layer can be tracked smoothly and continuously while ensuring computational accuracy and efficiency as much as possible. This avoids changes in the grid topology and thus ensures the quality of the numerical simulation.
[0085] As a preferred implementation, this embodiment decomposes the coupled evolution of "flow field-droplet-liquid film icing" into two time scales: small-step fluid-solid-thermal coupling propagation and large-step ice layer shape update. A small-step approach is adopted. This approach aims to advance rapidly changing processes such as airflow, droplet motion, and liquid film icing, ensuring the capture of instantaneous dynamic details. Simultaneously, macroscopic changes in the ice layer's shape are observed using larger time steps. To update, that is, after accumulating to a certain time interval. Afterwards, a single mesh deformation is performed to significantly save computational resources and further improve computational efficiency, demonstrating the characteristics of quasi-steady-state numerical simulation in this method.
[0086] In summary, the quasi-steady-state numerical simulation method for the evolution characteristics of icing morphology of railway catenary provided in this embodiment fully considers the entire process from microscopic droplet impact to macroscopic ice layer formation on the catenary surface, including the entire process of droplet impact, breakage, liquid film formation, phase change, and ice layer morphology growth. Compared with the prior art, it avoids prediction deviations caused by simplification of a single link. At the same time, by using the RBF interpolation function to dynamically deform the mesh, it realizes real-time updates of the ice layer shape during the simulation process, overcoming the limitations of the traditional fixed boundary method.
[0087] Meanwhile, this method also allows users to customize boundary conditions such as temperature, humidity, and wind speed to simulate the flow field around the contact network under different meteorological conditions. When injecting freezing rain droplets using the Lagrange method, users can set parameters such as initial velocity, diameter distribution, temperature, and specific heat capacity to reproduce different freezing rain intensities and environmental characteristics.
[0088] In addition, this method also adopts a dual timescale strategy of quasi-steady-state evolution and grid deformation scheduling, which significantly reduces the computational cost of long-term series simulation. While ensuring the computational accuracy of each stage, it reduces unnecessary geometric update frequency, making it possible to simulate the actual freezing rain icing process that lasts for several hours. The computational efficiency is significantly improved compared with single small-step coupled calculation.
[0089] The technical features of the above embodiments can be combined in any way. For the sake of brevity, not all possible combinations of the technical features in the above embodiments are described. However, as long as there is no contradiction in the combination of these technical features, they should be considered to be within the scope of this specification.
[0090] The above embodiments are merely illustrative of several implementation methods of this application, and their descriptions are relatively specific and detailed, but they should not be construed as limiting the scope of the application. It should be noted that those skilled in the art can make various modifications and improvements without departing from the concept of this application, and these all fall within the protection scope of this application. Therefore, the protection scope of this application should be determined by the appended claims.
Claims
1. A quasi-steady-state numerical simulation method for the evolution characteristics of icing morphology in railway overhead contact lines, characterized in that, Includes the following steps: S1, numerical simulation modeling, establish air flow field, define temperature, humidity, and wind speed boundary conditions, and simulate the flow field around the contact wire under different meteorological conditions; S2, In the numerical simulation model, freezing rain droplets are injected using the Lagrange method, and the parameters of the freezing rain droplets, including initial velocity, diameter distribution, temperature and specific heat capacity, are set to reproduce different freezing rain intensities and environmental characteristics; S3, the trajectory of the discrete phase droplets in the air flow field is constantly evolving. When the droplets move to the vicinity of the contact wire surface, it is calculated whether they will collide with the contact wire surface based on the air flow field velocity and the droplet inertia. S4. After impact, some droplets adhere to the surface of the contact network, converge to form a thin liquid film and spread on the surface. As time goes by and the surface temperature decreases, the water in the liquid film gradually freezes into ice. An ice accumulation model is used to simulate the freezing process of the liquid film. The ice accumulation model calculates the parameters of each attached droplet or liquid film unit based on the thin film energy equation, and calculates the latent heat exchange and extracts the freezing parameters to construct the growth mode of contact network icing. S5 uses the RBF interpolation function to smooth the mesh, realizing the dynamic evolution of the ice layer shape and obtaining numerical simulation results.
2. The quasi-steady-state numerical simulation method for the evolution characteristics of icing morphology of railway catenary according to claim 1, characterized in that, When establishing the airflow field in S1, the mass conservation equation, the momentum equation of the continuous phase, the energy conservation equation, and the water vapor transport equation are followed.
3. The quasi-steady-state numerical simulation method for the evolution characteristics of icing morphology of railway catenary according to claim 1, characterized in that, In S2, a Lagrange framework is used to track the discrete phase of freezing rain droplets, following the particle motion equation and the particle thermal balance equation.
4. The quasi-steady-state numerical simulation method for the evolution characteristics of icing morphology of railway catenary according to claim 1, characterized in that, Calculate and set the bidirectional coupling source terms for the discrete and continuous phases in S2.
5. The quasi-steady-state numerical simulation method for the evolution characteristics of icing morphology of railway catenary according to claim 1, characterized in that, The collision process between the droplets and the contact wire surface in S3 is handled using the SSD standard breakage model. The SSD standard breakage model is based on dimensionless parameters, including the Weber number, Ohnesorge number, and incident angle. The conditions for breakup and the diameter distribution, velocity distribution and number of droplets after breakup are given.
6. The quasi-steady-state numerical simulation method for the evolution characteristics of icing morphology of railway catenary according to claim 1, characterized in that, The thin film thickness evolution in S4 is represented as follows: in, For film thickness, The tangential divergence along the surface, For the wall flow velocity, For the mass flux of foreign attached droplets. To freeze the mass flux consumed, This refers to the mass flux lost due to slippage.
7. The quasi-steady-state numerical simulation method for the evolution characteristics of icing morphology of railway catenary according to claim 6, characterized in that, The thin film energy equation is expressed as: in, The density of liquid water, Specific heat capacity of the liquid The convective heat transfer coefficient is... For the film temperature, For the thermal conductivity of liquid, For free flow temperature, This is the latent heat of phase transition.
8. The quasi-steady-state numerical simulation method for the evolution characteristics of icing morphology of railway catenary according to claim 7, characterized in that, When using the ice accumulation model, the ice-liquid interface satisfies the Stefan condition: in, The density of ice, For ice thickness, The thermal conductivity of ice.
9. The quasi-steady-state numerical simulation method for the evolution characteristics of icing morphology of railway catenary according to claim 1, characterized in that, In S5, nodes on the contact wire surface and its adjacent areas are selected. A normal displacement is assigned to each surface node based on the thickness of the newly added ice layer. Then, the RBF interpolation function is used to calculate the position adjustment of each node within the entire computational domain, achieving continuous tracking of the ice growth process. The formula for calculating the node displacement within the domain using the RBF interpolation function is as follows: in, The coordinates of the target node. Let be the displacement vector of the target node. The number of all boundary nodes involved in the interpolation. These are the RBF interpolation coefficients. For the selected RBF kernel function, These are the coordinates of the boundary nodes. Let be the Euclidean distance from the target node to the boundary node. These are linear polynomial terms.
10. A quasi-steady-state numerical simulation method for the evolution characteristics of icing morphology in railway catenary according to any one of claims 1-9, characterized in that, The numerical simulation process is decomposed into two time scales: small-step fluid-solid-thermal coupling propulsion and large-step ice layer shape update. Small-step propulsion is used to advance the rapid changes in air flow field, droplet motion and liquid film icing, while larger time steps are used for macroscopic changes in ice layer shape. After accumulating the preset time interval, a mesh deformation is performed in one concentrated operation.