Contact system temperature rise and short circuit dynamic response simulation method and system

Abstract

The invention discloses a contact system temperature rise and short circuit dynamic response simulation method and system, and relates to the technical field of high-voltage electric appliance simulation. Comprising the steps that a finite element model and material attribute parameters of a contact system are obtained, steady-state electric-thermal coupling analysis is conducted, an initial temperature distribution field is obtained, transient multi-physics coupling analysis under the short-circuit working condition is conducted on the basis of the finite element model with the initial temperature distribution field as the initial condition, and in each time step of transient multi-physics coupling analysis, a transient multi-physics coupling analysis result is obtained. According to the contact force at the current moment and the evolution state of the contact interface, dynamic contact resistance is calculated, calculation of the dynamic contact resistance is based on the dynamic response characteristic of the contact interface under the force-thermal load, and short-circuit dynamic response data of the contact system is obtained based on the result of transient multi-physics coupling analysis. Closed-loop coupling of electricity, heat, force and microcosmic contact physics is achieved, and the dynamic response prediction precision and the numerical value stability of the contact system under the extreme short-circuit working condition are improved.

Description

Technical Field

[0001] This invention relates to the field of high-voltage electrical appliance simulation technology, and in particular to a simulation method and system for temperature rise and short-circuit dynamic response of a contact system. Background Technology

[0002] Vacuum circuit breakers are important protective devices in power systems, and the dynamic behavior of their contact system under short-circuit faults affects the breaking success rate. Under the impact of tens of thousands of amperes of short-circuit current, the contact interface will experience electromagnetic repulsion impact, rapid accumulation of Joule heat, and rapid evolution of micro-contact spots. These physical processes are complex and interconnected.

[0003] Existing simulation methods for contact systems are based on classical Holm contact theory, treating contact resistance as a static function dependent on contact pressure, i.e., a static algebraic model. During simulation, a quasi-static coupling strategy is employed, meaning changes in contact area and contact pressure are instantaneously synchronized, neglecting the time-response characteristics of the contact interface. This static model has limitations in handling high-frequency transient short-circuit conditions, failing to accurately describe the hysteresis effect of the contact response and numerical distortion under extreme conditions. Summary of the Invention

[0004] Objective of the invention: In view of the above-mentioned problems in the prior art, this application provides a simulation method and system for temperature rise and short-circuit dynamic response of contact system.

[0005] Technical solution: A simulation method for temperature rise and short-circuit dynamic response of a contact system, comprising:

[0006] Obtain the finite element model and material property parameters of the contact system;

[0007] Steady-state electro-thermal coupling analysis was performed based on the finite element model and material property parameters to obtain the initial temperature distribution field.

[0008] Using the initial temperature distribution field as the initial condition, transient multiphysics coupling analysis under short-circuit conditions is carried out based on the finite element model.

[0009] In each time step of the transient multiphysics coupling analysis, the dynamic contact resistance is calculated based on the contact force and the evolution state of the contact interface at the current moment. The calculation of the dynamic contact resistance is based on the dynamic response characteristics of the contact interface under force-thermal load.

[0010] Based on the results of transient multiphysics coupling analysis, the short-circuit dynamic response data of the contact system are obtained.

[0011] A simulation system is provided for implementing a simulation method for temperature rise and short-circuit dynamic response of a contact system. The system includes:

[0012] The model acquisition module is used to acquire the finite element model and material property parameters of the contact system;

[0013] The steady-state analysis module is used to perform steady-state electro-thermal coupling analysis based on the finite element model and material property parameters to obtain the initial temperature distribution field.

[0014] The transient coupling module is used to perform transient multiphysics coupling analysis under short-circuit conditions based on the finite element model, using the initial temperature distribution field as the initial condition.

[0015] The contact evolution calculation module is used to calculate the dynamic contact resistance at each time step of the transient multiphysics coupling analysis based on the contact force and the evolution state of the contact interface at the current moment. The calculation of the dynamic contact resistance is related to the dynamic response characteristics of the contact interface under force-thermal load.

[0016] The response analysis module is used to obtain the short-circuit dynamic response data of the contact system based on the results of the transient multiphysics coupling analysis.

[0017] Beneficial effects: This invention achieves closed-loop coupling of electricity, heat, force and microscopic contact physics, improving the accuracy and numerical stability of dynamic response prediction of contact system under extreme short-circuit conditions. Attached Figure Description

[0018] Figure 1 This is a flowchart of the simulation method for temperature rise and short-circuit dynamic response of the contact system in this application.

[0019] Figure 2 This is a flowchart illustrating the calculation of dynamic contact resistance in this application.

[0020] Figure 3 This is a flowchart illustrating how the initial temperature distribution field is obtained in this application.

[0021] Figure 4 This is a flowchart of the transient multiphysics coupling analysis of this application.

[0022] Figure 5 This is a flowchart illustrating the data processing after obtaining the short-circuit dynamic response data of the contact system in this application. Detailed Implementation

[0023] To enable those skilled in the art to better understand the present invention, the technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings of the embodiments of the present invention. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort should fall within the scope of protection of the present invention.

[0024] It should be noted that the terms "first," "second," etc., in the specification and accompanying drawings of this invention are used to distinguish similar objects and are not necessarily used to describe a specific order or sequence. It should be understood that such data can be interchanged where appropriate so that embodiments of the invention described herein are implemented in a sequence other than those illustrated or described herein. Furthermore, the terms "including" and "having," and any variations thereof, are intended to cover non-exclusive inclusion; for example, a process, method, system, product, or apparatus that includes a series of steps or units is not necessarily limited to those explicitly listed, but may include other steps or units not explicitly listed or inherent to the process, method, product, or apparatus.

[0025] To address the aforementioned issues, the applicant conducted in-depth searches and analyses, and discovered:

[0026] While some existing studies have considered multi-physics coupling, they often assume that material properties are constant or use simplified heat source models when dealing with short-circuit high-current conditions, without considering the influence of high temperature on the microstructure of contact materials or the quantum tunneling effect at the moment of contact separation.

[0027] Correspondingly, the microscopic deformation of the actual contact interface involves elastoplastic rheology, and the change in contact area lags behind the change in contact force. Existing static formulas cannot characterize the asymmetric hysteresis characteristics, leading to biases in the prediction of contact separation time and dynamic contact resistance. Furthermore, existing models do not consider the thermal damage softening effect caused by high temperature accumulation, nor do they correct for the singularity that the contact resistance tends to infinity at the moment of separation. This results in voltage spikes or non-convergence in the calculation results when simulating the zero-crossing point of the short-circuit current or the instant of contact repulsion, making it difficult to reflect the withstand capability of the contact system.

[0028] To solve these problems, combined with Figures 1 to 5 The present invention will be specifically described through the following embodiments.

[0029] In some embodiments, a scheme for simulating the temperature rise and short-circuit dynamic response of a contact system is provided. This scheme considers the interaction of electromagnetic, thermal, structural, and microscopic contact physics mechanisms to predict the dynamic behavior of the contact system under extreme short-circuit conditions. This application performs steady-state electro-thermal analysis based on the finite element model of the contact system; in other words, it obtains the initial temperature field and performs transient multiphysics coupling analysis under short-circuit conditions.

[0030] Furthermore, a contact state evolution model is introduced, using a differential evolution equation containing an asymmetric time constant to describe the dynamic hysteresis process of the microscopic contact area changing with contact force. Based on this, a thermal damage degradation factor and a tunneling effect correction term are introduced to physically correct the resistance calculations at the moment of high-temperature softening and contact separation, specifically including:

[0031] Step 101: Obtain the finite element model and material property parameters of the contact system.

[0032] Specifically, the process of constructing the finite element model is completed in commercial finite element simulation software, such as ANSYS finite element analysis software. Based on the actual geometric dimensions of the contact system, a model is established including the moving contact, stationary contact, conductive rod, contact fingers, and insulating support structure. For structures far from the contact area, a sparser mesh is used; for the contact interface and its vicinity, local mesh refinement is applied. For example, the mesh size for the contact area is set between 0.1 mm and 0.5 mm to capture the physical field gradient changes at microscopic contact points. Regarding the element type selection, elements supporting multiphysics coupling are used. For example, in ANSYS finite element analysis software, SOLID226 or SOLID227 elements are selected. These elements have voltage, temperature, and displacement degrees of freedom and are used to solve the electro-thermal-structural coupling equations.

[0033] Furthermore, regarding material property parameters, the physical characteristics of each component material are defined, particularly its nonlinear characteristics as a function of temperature. The contact material is a copper-chromium alloy, such as CuCr25 or CuCr50. Taking CuCr25 as an example, its density is approximately 8.9 g / cm³. 3 The Poisson's ratio is approximately 0.3. Thermoelectric parameters are defined as functions of temperature, such as resistivity and thermal conductivity. For example, at 20°C, the resistivity is set to approximately 3.5 μΩ·cm, and the thermal conductivity to approximately 350 W / (m·K); while at 800°C, the resistivity increases, and the thermal conductivity decreases. Next, parameters such as the material's Young's modulus, coefficient of thermal expansion, and specific heat capacity are defined for the calculation of structural and thermal fields. These parameters can be input into the simulation system in tabular form, and the system will automatically interpolate the corresponding material properties based on the temperature of the current node during the calculation.

[0034] Step 102: Based on the finite element model and material property parameters, perform steady-state electro-thermal coupling analysis to obtain the initial temperature distribution field.

[0035] Specifically, boundary conditions are applied to both ends of the circuit in the finite element model, such as applying a rated current to one end and grounding the other. Convective heat transfer boundary conditions are applied to the outer surface of the model to simulate natural air cooling or forced air cooling. For example, the ambient temperature is set to 20°C, and the natural convection heat transfer coefficient is set to 10 W / (m²·K).

[0036] Next, the solver uses direct coupling or iterative coupling to solve the steady-state electrical conduction equation and heat conduction equation. The current flowing through the conductor generates Joule heat, which is added to the heat conduction equation as a heat source term. As the temperature increases, the material resistivity increases, altering the current distribution and the magnitude of the Joule heat. The system iterates until the temperature field distribution reaches convergence. The converged result is the initial temperature distribution field, representing the thermal state of the contact system just before a short-circuit fault occurs.

[0037] Step 103: Using the initial temperature distribution field as the initial condition, perform transient multiphysics coupling analysis under short-circuit conditions based on the finite element model.

[0038] After completing the steady-state analysis, the simulation enters the transient analysis phase. The node temperature data obtained from the steady-state analysis is mapped and loaded into the initial time step of the transient analysis. Next, the circuit boundary conditions are modified, replacing the input current load from the rated current with a preset short-circuit current waveform. The short-circuit current waveform is a sine wave superimposed with a decaying DC component, and its peak value can reach tens of thousands of amperes.

[0039] Furthermore, the time step for transient analysis is set to the microsecond level to capture rapidly changing physical quantities during short circuits, such as 10 microseconds or less. For structural dynamics solutions, large deformation is enabled, and contact pair properties, such as the friction coefficient, are set. The solver at this stage needs to simultaneously process the electromagnetic field equations, heat conduction equations, and structural dynamics equations, calculating the Lorentz force, Joule heating, and dynamic displacement of the contacts generated by the short-circuit current.

[0040] Step 104, wherein in each time step of the transient multiphysics coupling analysis, the dynamic contact resistance is calculated based on the contact force and the evolution state of the contact interface at the current moment. The calculation of the dynamic contact resistance is based on the dynamic response characteristics of the contact interface under force-thermal load.

[0041] Within each time step, the finite element solver calculates the Lorentz force based on the current current distribution and applies it to the moving contact, causing movement of the moving contact and changes in contact pressure. The system reads the current contact force values ​​at the contact interface and updates the contact resistance by combining them with the unique evolutionary state model of the contact interface.

[0042] The evolutionary state depends not only on the current contact pressure but also on the historical loading path and thermal damage state of the contact interface. The system introduces an internal state variable to characterize the proportion of the microscopic contact area; this state variable evolves over time and exhibits a memory effect. Calculating the dynamic contact resistance using this method reflects the hysteresis characteristics of the contact interface during rapid loading and unloading, as well as the thermal softening effect caused by high temperatures. The calculated dynamic contact resistance is assigned to the contact unit, affecting the current density distribution and contact heat generation at that location, achieving a closed-loop coupling of electrical, magnetic, mechanical, thermal, and contact microphysics.

[0043] Step 105: Based on the results of transient multiphysics coupling analysis, obtain the short-circuit dynamic response data of the contact system.

[0044] Specifically, after a series of iterative calculations over a series of time steps, the simulation ends and outputs complete dynamic response data. This data includes the temperature rise curves of each measuring point in the contact system over time; the displacement curve of the moving contact over time, which is used to determine the contact repulsion phenomenon; and the contact resistance curve over time.

[0045] Next, by analyzing this data, designers assess the contact system's resistance to welding and its dynamic and thermal stability under specified short-circuit currents. For example, they examine the displacement data at the moment of contact separation to determine the contact separation time; and they check the peak value of the temperature rise curve to determine whether it exceeds the material's melting point or softening temperature.

[0046] According to one aspect of this application, a contact resistance calculation model based on differential evolution equations is also proposed, comprising:

[0047] Step 201: Calculate the dynamic contact resistance based on the current contact force and the evolution state of the contact interface. Specifically, this includes introducing a contact state variable ξ that characterizes the normalized effective contact area of ​​the contact interface.

[0048] Here, a physical quantity ξ is defined to quantify the microscopic contact state, i.e., the contact state variable. From a microscopic perspective, a smooth contact surface is actually composed of countless tiny peaks. When two contacts come into contact, the actual current flow area is limited to the patchy region where these tiny peaks contact each other, i.e., the effective contact area A. _eff The contact state variable ξ is defined as the effective contact area A. _eff With nominal contact area A _apparent The ratio of (or the maximum possible contact area when closed) is a dimensionless scalar between 0 and 1.

[0049] When ξ approaches 1, it indicates that the contact interface is in a tight fit and the conductive channel is open; when ξ approaches 0, it indicates that the contact interface is in a separated or weak contact state and the conductive channel is close to closed. This is because the traditional contact resistance formula is based on the contact force F. _c With resistance R _c In the meantime, the change in contact area is ignored. However, in reality, the change in contact area lags behind the change in contact force due to the viscoelasticity and plastic deformation of the material. Therefore, ξ is introduced as an intermediate state variable to decouple the effect of contact force on resistance into two processes: the force changes the contact area, and the contact area determines the resistance.

[0050] In another embodiment, the contact state variable ξ can also be ξ = A. _eff / A _ref ;where A _eff A is the sum of the areas of all microscopic contact spots participating in conductivity at the current moment; _ref For reference contact force F _ref The effective contact area after reaching steady state under the action.

[0051] Step 202: Determine the equilibrium contact state ξ of the contact interface based on the contact force at the current moment. _eq .

[0052] Equilibrium contact state ξ _eq This indicates that under a given constant contact force F _c Under the influence of the action, the contact interface reaches a stable contact state after a required time. This is the target value or attractor in the contact state evolution process. This state is determined based on classical contact mechanics theories, such as the Greenwood-Williamson (GW) model or Hertzian contact theory.

[0053] Accordingly, the equilibrium contact state ξ of the contact interface is determined. _eq The equilibrium contact state is calculated using a power-law relation, and its expression is:

[0054] ξ _eq =(F _c / F _ref ) n ;

[0055] Among them, F _c F is the contact force at the current moment. _ref The reference contact force is denoted by η, and the contact deformation index is denoted by η. The value of η is determined based on the hardness and surface microstructure of the contact material: when the contact is mainly elastic deformation, η is 2 / 3; when the contact is mainly plastic deformation, η is 1; and when the contact is a mixed elastic-plastic deformation, η is between 2 / 3 and 1. _c The macroscopic contact force at the current moment is calculated by transient analysis using the finite element model; F_ref The characteristic pressure required to bring the contact interface to a closed state is taken.

[0056] Based on this, the value of the contact deformation index *n* is determined by the hardness and surface microstructure of the contact material. For example, when the deformation of the contact micro-protrusion is elastic, according to Hertz's contact theory, the contact area is proportional to the 2 / 3 power of the pressure, and *n* is set to 2 / 3, approximately 0.667. When the contact pressure is high and the micro-protrusion undergoes plastic yielding, the contact area and pressure are linearly related, and *n* is set to 1. In practical engineering, the contact state is between elastic and plastic, i.e., a mixed elastic-plastic deformation mode, and *n* is set to between 2 / 3 and 1. For another example, for copper contacts that have undergone surface hardening through multiple separation and joining operations, *n* is set to 0.8. Through this segmented or interpolated parameter setting, the model can adapt to contact conditions with different material hardness and surface roughness. Among these parameters, the contact force *F*... _ref Determined using experimental calibration or theoretical estimation methods.

[0057] Step 203: Invoke the differential evolution equation of the contact state variable ξ with respect to time, and update the contact state variable ξ for the next time step based on the differential evolution equation;

[0058] Traditional static models assume that ξ is instantaneously equal to ξ _eq The changes in contact states require time, constituting a dynamic evolution process. Therefore, a first-order differential evolution equation is constructed:

[0059] dξ / dt=(ξ _eq -ξ) / τ;

[0060] Where dξ / dt is the rate of change of the contact state variable with time; t is time; ξ _eq ξ represents the equilibrium contact state calculated based on the current contact force, i.e., the target value of the evolution; ξ represents the current contact state variable; τ represents the time constant, characterizing how quickly the contact state tends to the equilibrium state.

[0061] The evolution equation states that the rate of change of the contact state is proportional to the deviation between the current state and the target state. When the contact force suddenly increases, ξ... _eq During the leap, the actual ξ does not change immediately, but rather, with τ as the characteristic time constant, it moves exponentially towards the new ξ. _eq Asymptotic. This equation conforms to the characteristics of the Voigt-Kelvin model in materials rheology and can capture the delayed effect of contact response.

[0062] In numerical simulation, the differential equation is discretized and solved. Based on this equation, within each simulation time step Δt, the solution is obtained by considering ξ(t) from the previous time step and the currently calculated ξ. _eqThen, calculate ξ(t+Δt) for the next moment. Here, ξ(t) represents the contact state determined in the previous moment, and ξ... _eq ξ(t+Δt) represents the current equilibrium contact state, and ξ(t+Δt) represents the contact state at the next moment.

[0063] Step 204: Calculate the dynamic contact resistance based on the updated contact state variable ξ and the current temperature of the contact interface.

[0064] Specifically, the updated microscopic contact state variable ξ is obtained and mapped to the macroscopic dynamic contact resistance R. _c Contact resistance is physically composed of contraction resistance and film resistance, and the mapping relationship is as follows:

[0065] R _c =R _c0 ×ξ -α ;

[0066] Among them, R _c Indicates dynamic contact resistance; R _c0 This represents the reference contact resistance. It is the contact resistance value when ξ=1, which is determined by the resistivity of the material and the contact geometry; α is the contact resistance exponent, which ranges from 0.5 to 1.0 and depends on the distribution pattern of the contact spots.

[0067] Furthermore, temperature affects resistivity, therefore a temperature term is introduced, and the relationship is as follows:

[0068] R _c =R _c0 ×(1+γ×(T _c -T _0 ))×ξ -α +R _tunnel ;

[0069] Among them, T _c T represents the current temperature of the contact interface. _0 Reference temperature; R _tunnel ξ represents the tunneling effect resistance term; γ is the temperature coefficient of resistance of the material. Through the relationship including the temperature term, the contact state variable ξ and temperature T are... _c Commonly affecting resistance R _c .

[0070] In this application, the dynamic contact resistance can also be:

[0071] R _c =R _c，_ref ×ξ -0.5 ×sqrt(1+α _ρ ×(T _c -T _ref ));

[0072] Among them, R _c，_ref For ξ=1 and T _c =T _ref The reference contact resistance value at that time; ξ is the contact state variable; α _ρ T is the temperature coefficient of resistivity of the contact material; _c The contact interface temperature; T _ref The reference temperature is sqrt, which stands for square root.

[0073] Step 205, wherein the differential evolution equation characterizes the contact state variable ξ as it evolves with time toward the equilibrium contact state ξ. _eq A dynamic process of gradual change.

[0074] Specifically, in the quasi-static method, the time term dt is ignored, and the waveform of the change in contact resistance is synchronized with the waveform of the contact force (distinguishing only by a proportionality coefficient or exponential relationship). However, due to the existence of differential equations, even if the contact force F... _c It is an ideal step signal, and the calculated contact resistance R _c A transition process will also occur; when the contact force F _c When high-frequency oscillation occurs, the calculated R _c It exhibits amplitude decay and phase lag. This dynamic characteristic is consistent with physical experimental observations, such as the high-frequency micro-vibration of the contacts when the short-circuit current generates an electrodynamic repulsion force. This model can avoid non-physical, drastic resistance jumps in the simulation results, thus improving the stability of numerical calculations.

[0075] Correspondingly, the differential evolution equation also possesses an intrinsic natural regularization effect. When the contact force F _c When the short-circuit repulsion force suddenly drops to zero, ξ _eq It also becomes zero. If an algebraic model is used, the contact resistance will instantaneously jump to infinity, causing numerical collapse. In this evolution equation, the change of ξ is affected by the decay time constant τ. _dec The constraint, whose decay behavior approximately satisfies:

[0076] ξ(t)=ξ _0 ×exp(-t / τ _dec );

[0077] Where ξ(t) is the contact state variable at time t; _0 ξ is the value of ξ at the instant the contact force drops to zero; τ _dec is the decay time constant; t is the time from that instant; exp(...) is the natural exponential function.

[0078] Based on this, time t is a finite value, ξ(t) is a positive number, and R _c Since the values ​​are finite, singularities will not occur. This evolutionary equation provides the first layer of numerical protection.

[0079] Since static algebraic formulas cannot reflect the dynamic hysteresis effect under short-circuit impacts when calculating contact resistance, this embodiment introduces normalized contact state variables and their differential evolution equations to establish a dynamic model of the gradual change of the microscopic physical state of the contact interface over time.

[0080] According to another aspect of this application, an implementation method for asymmetric hysteresis characteristics and numerical solutions is also proposed. In actual physical processes, during loading, the crushing of the microscopic protrusions of the contact involves plastic deformation, which is relatively fast; during unloading, elastic recovery relies on elastic rebound, which is relatively slow. This embodiment introduces an asymmetric time constant to describe this phenomenon and proposes an exponential integral algorithm for numerical stability.

[0081] Step 301: The differential evolution equation contains growth and decay terms. That is, the time constant in the differential evolution equation exhibits asymmetric hysteresis characteristics, including the growth time constant τ. _inc and decay time constant τ _dec .

[0082] The time constant τ in the evolution equation is not a single constant, but rather the growth time constant τ _inc and decay time constant τ _dec The two time constants are defined based on the changing trend of the contact state.

[0083] In other words, as the contact pressure increases, the micro-contact area increases, i.e., the contact state variable ξ increases, and the micro-protrusions undergo yielding and rheological changes. This process is determined by the plastic flow velocity of the material and responds rapidly. The corresponding time constant is defined as the growth time constant τ. _inc When the contact pressure decreases, the microscopic contact area decreases, meaning the contact state variable ξ decreases. The separation of the contact spots is based on the release of elastic potential energy within the material. This rebound process is hindered by surface adhesion and the cold welding effect, resulting in a slow response. The corresponding time constant at this time is defined as the decay time constant τ. _dec .

[0084] By distinguishing between two constants, the model simulates an asymmetric hysteresis loop that loads quickly and unloads slowly, which can predict the microscopic dynamic behavior of the contacts near the zero-crossing point of the short-circuit current.

[0085] Step 302: The growth time constant corresponding to the growth term is less than the decay time constant corresponding to the decay term, characterizing the hysteresis characteristics of the contact interface during loading and unloading. Growth time constant τ _inc Less than the decay time constant τ _dec This is used to simulate a faster loading process on the touch interface than a faster unloading process.

[0086] Specifically, τ _inc <τ _decThis aligns with actual operating conditions. For example, in a copper contact system, the growth time constant τ... _inc The time constant τ is set between 1 microsecond and 10 microseconds. _dec The time is set between 50 and 100 microseconds. When the electrodynamic repulsion force generated by the short-circuit current increases rapidly, the contact resistance decreases rapidly because τ... _inc Small resistance and fast response; however, when the current crosses zero and the repulsive force disappears rapidly, the contact resistance does not immediately return to a high resistance state, but remains in a low resistance state for a period of time, because τ _dec Larger size, slower recovery.

[0087] This effect manifests macroscopically as the dynamic engagement or suspension of the contacts, which can explain why, under certain rapid transient conditions, the contacts separate under force but do not arc. This dynamic process cannot be captured if a single time constant is used.

[0088] Furthermore, the value of the time constant is based on the spatial resolution of the model. For a microscale model with a single micro-convexity as the object of analysis, its growth time constant τ _inc The value is between 1 microsecond and 10 microseconds, and the decay time constant τ _dec The value ranges from 50 microseconds to 200 microseconds. For a macroscopic model where the entire contact surface is the object of analysis, the equivalent time constant increases due to the statistical dispersion effect of the micro-protrusion response time. In this embodiment, τ _inc Set to 5 microseconds, τ _dec Set to 80 microseconds. In specific simulations, the parameter values ​​at the corresponding scale are selected based on the spatial resolution, and τ must be satisfied. _inc <τ _dec .

[0089] Step 303, the differential evolution equation is expressed as:

[0090] dξ / dt=H(ξ _eq -ξ)×((ξ _eq -ξ) / τ _inc )+H(ξ-ξ _eq )×((ξ _eq -ξ) / τ _dec );

[0091] Where dξ / dt is the rate of change of the contact state variable with time, t is time; H(...) is the unit step function; ξ is the current contact state variable; _eq The equilibrium contact state is calculated based on the current contact force; τ _inc τ is the growth time constant, representing the rate at which the contact area increases; _dec τ is the decay time constant, characterizing the rate at which the contact area decreases, and τ _inc <τ _dec .

[0092] The differential evolution equation is expressed as a piecewise function, when ξ < ξ _eq dξ / dt=(ξ _eq -ξ) / τ _inc When ξ>ξ _eq dξ / dt=(ξ _eq -ξ) / τ _dec ;

[0093] Based on the current contact state quantity ξ and the target equilibrium contact state ξ _eq The relative magnitudes of ξ and ξ determine the time constant of the evolution equation. That is, when ξ < ξ _eq This indicates that the current effective contact area is less than the target equilibrium area, and the system is in the loading or growth phase. A smaller growth time constant τ is selected. _inc This makes dξ / dt larger, allowing ξ to catch up with ξ quickly. _eq When ξ>ξ _eq When the current effective contact area is greater than the target equilibrium area, the system is in the unloading or decay phase, and a larger decay time constant τ is selected. _dec This makes the absolute value of dξ / dt smaller, and ξ slowly decreases towards ξ. _eq .

[0094] The piecewise linear differential equations construct a direction-dependent dynamic system. The introduction of nonlinearity ensures that the system's output depends not only on the current input but also on the input's historical path.

[0095] Step 304: Update the contact state variable ξ for the next time step based on the differential evolution equation, specifically using an exponential time integral scheme:

[0096] ξ(t+Δt)=ξ _eq +(ξ(t)-ξ _eq )×exp(-Δt / τ);

[0097] Where t is the current time, Δt is the time step, and τ is the growth time constant or decay time constant selected based on the current evolution trend.

[0098] To simulate absolute stability at any time step, this embodiment employs an exponential time integral scheme.

[0099] Assume that within the current time step Δt, the equilibrium state ξ _eq The time constant τ remains constant. The differential equation is dξ / dt = (ξ... _eq -ξ) / τ has a solution, which can be expressed in discrete form as:

[0100] ξ _new =ξ _eq +(ξ_old -ξ _eq )×exp(-Δt / τ _curr );

[0101] Where, ξ _new Let ξ be the contact state variable value at the next time step t+Δt; _old τ is the value at the current time t; Δt is the simulation time step; τ _curr According to the criterion of the differential evolution equation, the value is τ. _inc or τ _dec .

[0102] For example, at a certain moment, the contact system is in a steady state, ξ _old =0.5. Due to the impact of the short-circuit current, the contact force suddenly increases. The new equilibrium contact state value is calculated using the power-law formula, i.e., ξ. _eq =0.8.

[0103] Set the growth time constant τ _inc =10 microseconds, simulation time step Δt=2 microseconds.

[0104] Because ξ _old <ξ _eq (i.e., 0.5 < 0.8), indicating that it is in the loading stage, τ = 10 microseconds is selected.

[0105] Substitute the values ​​into the exponential integral formula to calculate the value at the next time step:

[0106] ξ _new =0.8+(0.5-0.8)×exp(-2 / 10)=0.5544.

[0107] Therefore, after one time step, ξ increases from 0.5 to 0.5544, instead of jumping directly to 0.8. With subsequent time steps, ξ will approach 0.8 exponentially. If a traditional algebraic model were used, the value would jump directly to 0.8, losing the dynamic transition from 0.3.

[0108] As an alternative, if the simulation time step Δt is much smaller than the time constant τ, for example, Δt = 0.01 × τ, the first-order explicit Euler method can also be used for approximate calculation:

[0109] ξ _new =ξ _old +(ξ _eq -ξ _old )×(Δt / τ);

[0110] When Δt > 2 × τ, the Euler method will cause the calculation results to diverge, while the exponential integral scheme converges for any Δt. When thermal softening correction terms are involved, the equilibrium contact state is replaced with the effective equilibrium state after thermal softening correction. Therefore, the complete exponential integral discretization formula includes mechanical evolution terms and thermal softening correction terms. The corresponding calculation process includes:

[0111] Calculate the effective equilibrium contact state including thermal softening correction:

[0112] ξ _eq，_eff k =(F _c k / F _ref ) n ×Θ _soft (T _max k );

[0113] Where, ξ _eq，_eff k The effective equilibrium contact state at time step k; F _c k F is the contact force at the k-th time step. _ref The reference contact force is Θ; n is the contact deformation index; Θ _soft (T _max k The value of the thermal softening factor is calculated based on the historical highest temperature at the k-th time step. This factor characterizes the effect of high temperature on the material's yield strength, leading to a decrease and an increase in the effective contact area under the same contact force. Its value is greater than or equal to 1; the superscript k indicates the time step number. The effective time constant τ is determined based on the current evolution trend. _eff When ξ k <ξ _eq，_eff k , τ _eff =τ _inc When ξ k ≥ξ _eq，_eff k , τ _eff =τ _dec .

[0114] Based on this, exponential time integration is performed to obtain the contact state variables for the next time step:

[0115] ξ k+1 =ξ _eq，_eff k +(ξ k -ξ _eq，_eff k )×exp(-Δt / τ _eff );

[0116] Where, ξ k+1ξ represents the contact state variable value at the (k+1)th time step; k ξ represents the contact state variable value at the k-th time step; _eq，_eff k The effective equilibrium contact state includes thermal softening correction; Δt is the time step in seconds; τ _eff is the effective time constant, in seconds; exp(...) is the natural exponential function.

[0117] The mechanically driven equilibrium evolution, the equilibrium upsizing caused by thermal softening, and the asymmetric time response are unified in the exponential integral.

[0118] Thermal softening factor Θ _soft The calculation formula is:

[0119] Θ _soft (T _max )=1+β×exp((T _max -T _soft ) / T _scale );

[0120] Among them, T _soft T is the softening initiation temperature of the material. _scale β is a temperature-scale constant used to control the rate at which the softening effect increases with temperature; β is the softening strength coefficient. Taking copper as an example, let T be... _soft =400℃, T _scale =100℃, β=0.5. When T _max Less than T _soft When the exponent term is very small, Θ _soft The value is approximately 1, meaning no thermal softening effect occurs; when T _max More than T _soft At that time, Θ _soft The rapid increase leads to the effective equilibrium contact state ξ. _eq，_eff The increase reflects the physical phenomenon that high-temperature softening leads to a decrease in the material's yield strength and an increase in the contact area under the same contact force. We can take ξ as... _eq，_eff =min(ξ _eq，_eff ,1), so that the contact state variable does not exceed 1.

[0121] Accordingly, the system monitors the temperature T of the contact interface at each time step of the simulation. _curr Update the historical highest temperature variable T _max Its formula is:

[0122] T _max (t)=max(T _max (t-Δt), T _curr (t));

[0123] Among them, T _max(t) represents the updated historical high temperature; T _max (t-Δt) represents the historical highest temperature at the previous time step; T _curr (t) represents the current contact interface temperature; max(...) is a function to take the larger value. _max These are local historical variables maintained independently for each contact node. The corresponding T values ​​vary depending on the current density and heat dissipation conditions at different contact node locations. _max The softening levels differ. In the finite element method, the integration point of each contact element needs to be stored separately for T. _max Scalar value.

[0124] Thermal softening factor Θ _soft The independent variable is the historical highest temperature T. _max When the temperature exceeds the material's recrystallization temperature (300-450°C for copper-chromium alloys), an annealing effect occurs, resulting in the elimination of dislocation networks, grain growth, and a decrease in material hardness and yield strength. This metallurgical process is thermodynamically irreversible; even if the temperature subsequently decreases, the microstructural changes will not be reversed. Therefore, the degree of softening depends on the highest temperature experienced in the material's history.

[0125] The revised formula for calculating the effective equilibrium contact state can explain why, after a short-circuit impact, the contact resistance of the contact becomes lower or the contact becomes tighter. This is because high-temperature softening causes the material to undergo greater plastic deformation under the same contact force, resulting in an irreversible increase in the effective contact area.

[0126] In this embodiment, the effective equilibrium contact state with thermal softening correction can also be calculated using the following formula:

[0127] ξ _eq，_eff k =(F _c k / F _ref n ×[1+β×exp((T _max k -T _soft ) / T _scale )];

[0128] Among them, T _soft T is the softening initiation temperature of the material. _scale F is the temperature scale constant. _c k F is the contact force at the k-th time step. _ref Reference contact force; superscript n represents the contact deformation index; β represents the softening strength coefficient; T _max k This represents the highest historical temperature at the k-th time step; the superscript k is the time step number.

[0129] In another alternative implementation, a contact degradation model based on a thermal damage attenuation term is also proposed. This implementation is related to the thermal softening factor Θ. _soft The methods for correcting the equilibrium contact state differ: the thermal softening factor is increased by ξ. _eq The first term reflects the effect of high temperature increasing the contact area; the second term, on the other hand, directly acts on the contact state variable ξ itself, causing it to decay with damage accumulation, reflecting how irreversible structural damage caused by high temperature can reduce contact effectiveness. These two approaches are applicable to different temperature ranges and damage mechanisms, specifically including:

[0130] Step 401, in this optional embodiment, the differential evolution equation further includes a thermal damage attenuation term, which is used to characterize the attenuation effect of irreversible structural damage caused by extreme high temperature on contact effectiveness.

[0131] The differential evolution equation is expressed as:

[0132] dξ / dt=f _mech (ξ,ξ _eq )-(ξ / τ _deg )×Θ _deg (T _max );

[0133] Among them, f _mech (ξ,ξ _eq ) represents the sum of the growth and decay terms based on the mechanical mechanism; τ _deg Θ is the thermal degradation time constant; _deg (T _max ) is the thermal damage degradation function, whose value is related to the highest historical temperature T at the contact interface. _max Attenuation term -(ξ / τ) _deg )×Θ _deg (T _max ) drives ξ to decay downwards, Θ _deg The larger the value of ξ, the faster the decay rate of ξ, indicating more severe structural damage.

[0134] Specifically, during a high-current short circuit, when the temperature of the contact spot rises sharply, exceeding the material's softening point and approaching or reaching its melting point (e.g., above 800°C for copper-chromium alloys), localized melting, surface erosion, and microstructural damage occur at the contact spot. This thermal damage effect differs from the thermal softening effect: thermal softening softens the material and increases the contact area, while thermal damage causes irreversible structural degradation of the contact spot, leading to a decrease in contact effectiveness. Therefore, in the attenuation term, Θ... _deg The increase in ξ drives a decrease, indicating that the effective integrity of the contact interface is irreversibly reduced due to overheating damage.

[0135] Step 402, in this optional embodiment, Θ _deg (T _max) is the thermal damage degradation function, whose value is related to the highest historical temperature T at the contact interface. _max When the historical highest temperature T _max When the material's damage temperature threshold is exceeded, the thermal damage degradation function Θ _deg (T _max The value of ) increases exponentially with temperature.

[0136] Thermal damage degradation function Θ _deg The highest temperature ever recorded (T) _max The monotonically increasing function is used to characterize the degree of structural degradation of the contact spots caused by extreme high temperatures. In actual operating conditions, even if the temperature has decreased at the current moment, the structural damage to the contact spots will be preserved if extreme high temperatures have been reached in the past.

[0137] Specifically, the system monitors the temperature T of the contact interface at each time step of the simulation. _curr Update the historical highest temperature variable T _max Its formula is:

[0138] T _max (t)=max(T _max (t-Δt), T _curr (t));

[0139] Among them, T _max (t) represents the updated historical high temperature; T _max (t-Δt) represents the historical highest temperature at the previous time step; T _curr (t) represents the current contact interface temperature; max(...) is the function to take the larger value.

[0140] Based on T _max Construct a thermal damage degradation function to calculate Θ _deg :

[0141] Θ _deg (T _max )=β _deg ×exp((T _max -T _dam ) / T _scale，_deg );

[0142] Among them, T _dam This is the damage initiation temperature of the material, and its value is higher than the softening initiation temperature T in the thermal softening factor. _soft This corresponds to the extreme high temperature range; T _scale，_deg β is the damage temperature scale constant; _deg Let T be the damage strength coefficient. Taking copper-chromium alloy as an example, let T be... _dam =700℃, T _scale，_deg =150℃, β _deg =0.3. When T_max Much smaller than T _dam At that time, Θ _deg Approximately 0, the thermal damage attenuation term has no effect; when T _max More than T _dam At that time, Θ _deg The rapid increase drives ξ to accelerate its decay, reflecting the physical process of irreversible failure of the contact spot due to local melting and microstructure damage.

[0143] Accordingly, T _max This refers to a local historical variable maintained independently for each contact node. The T values ​​at different contact node locations vary due to differences in current density and heat dissipation conditions. _max The damage levels also differ depending on the contact element. In the finite element method implementation, the integration point T of each contact element needs to be stored separately. _max Scalar value.

[0144] According to another aspect of this application, a modification to the tunneling effect is also proposed, including:

[0145] Determine whether the contact state variable ξ is less than the tunneling effect activation threshold;

[0146] If so, then a tunneling effect resistance term is introduced to correct the contact resistance, and the dynamic contact resistance R _c The calculation is as follows:

[0147] R _c =R _c0 ×(1+γ×(T _c -T _0 ))×ξ -α +R _tunnel ;

[0148] R _tunnel =A×exp(B×(1-ξ / ξ _tunnel ));

[0149] Among them, R _c0 γ is the reference contact resistance; T is the temperature coefficient of resistance of the material; _c T represents the current temperature of the contact interface. _0 Reference temperature; α is the contact resistance index; R _tunnel For tunneling effect resistance; ξ _tunnel The tunneling effect activation threshold is defined as the threshold value. When the contact state variable ξ is less than this threshold, it indicates that the contact interface is in a microscopic separation transition stage. In this embodiment, ξ... _tunnelSetting it to 0.05 indicates that when the normalized effective contact area is less than 5%, the contact interface is considered to have entered the microscopic separation transition stage, at which point the tunneling effect resistance correction term is activated; A and B are tunneling effect constants. is the tunneling effect resistance reference value, characterizing the order of tunneling resistance of the contact interface in a fully microscopic separation state; B is the tunneling effect attenuation coefficient, which is highly correlated with the contact interface barrier. A has dimensions, while B is a dimensionless parameter.

[0150] At the moment of contact separation, a tunneling effect resistance R is introduced. _tunnel Regularization correction is applied to the dynamic contact resistance. This is done by adjusting the shrinkage resistance term R. _c0 ×ξ -α Multiplied by a temperature correction factor (1+γ×(T)) _c -T _0 The temperature dependence that keeps the contraction resistance term consistent between the tunneling effect branch and the non-tunneling effect branch allows the dynamic contact resistance formula to be in the form ξ=ξ _tunnel Continuity at that point.

[0151] If not, ξ≧ξ _tunnel The dynamic contact resistance is calculated as follows:

[0152] R _c =R _c0 ×ξ -α ×(1+γ×(T _c -T _0 )),

[0153] The tunneling effect resistance term is not introduced.

[0154] In numerical simulations, as the contacts tend to separate, the contact state variable ξ approaches 0. According to formula R... _c =R _c0 ×ξ -α When ξ→0, R _c →∞. In numerical calculations, infinite resistance indicates a singularity in the conductance matrix, which can cause the solver to crash or produce non-physical voltage spikes.

[0155] When two metal surfaces are extremely close but not fully in contact, for example, with a gap of less than 5 nanometers, electrons tunnel through the potential barrier via the quantum tunneling effect, forming a weak tunneling current. This phenomenon shows that although the resistance is large, it is not infinite. To address the numerical stability issue, a parallel tunneling resistance term R is introduced. _tunnel Alternatively, you can add regularization terms.

[0156] In this embodiment, ξ _tunnel Setting it to 0.05 indicates that when the normalized effective contact area is less than 5%, it is considered that the contact interface has entered the microscopic separation transition stage, at which point the tunneling effect resistance correction term is activated.

[0157] Step 403: The revised dynamic contact resistance calculation formula satisfies the following: when ξ approaches 0, the resistance value tends to a finitely large tunnel resistance value, thus avoiding divergence in numerical calculation.

[0158] Specifically, the tunnel resistance can be considered as being in parallel with the contraction resistance, or a small amount can be added to the denominator.

[0159] The corrected formula is:

[0160] R _c，_total =R _c0 ×(ξ+ξ _reg ) -α ;

[0161] Where, ξ _reg It is a regularization constant, characterizing the ratio of the equivalent contact area contributed by the tunneling effect, and its value is set at 10. -6 Up to 10 -4 between.

[0162] In some scenarios, let R be... _c0 =10μΩ, α=1.

[0163] Under the uncorrected operating condition, when ξ decreases to 0, R _c =10×0 -1 =∞, the calculation overflows.

[0164] Under the revised operating conditions: Set ξ _reg =10 -5 .

[0165] When ξ=0, i.e., separation, R _c，_total =10×(0+10 -5 ) -1 =10×10 5 =1000000μΩ=1Ω.

[0166] When ξ=0.5, i.e., normal contact, R _c，_total =10×(0.5+0.00001) -1 ≈20μΩ.

[0167] Based on this, within the normal contact range, ξ is greater than ξ. _reg Correction term ξ _reg The effect can be ignored; for example, the difference between 20μΩ and 19.9996μΩ is negligible. At the moment of separation, for example, when ξ→0, limiting the resistance to 1Ω makes the high-voltage circuit simulation approximate as an open circuit, which can make the matrix invertible and allow the simulation to pass through the moment of contact separation.

[0168] According to one aspect of this application, a thermal softening factor can also be introduced, which characterizes the amplification effect of the decrease in material yield strength caused by high temperature on the contact area, and its value varies with the highest historical temperature T. _max The thermal softening factor increases with the increase of ξ, and its value ranges from 1 to 1. When no thermal softening occurs, the thermal softening factor is approximately equal to 1. When the historical highest temperature exceeds the softening temperature of the material, the thermal softening factor is greater than 1. This factor acts on the equilibrium contact state ξ. _eq This increases the effective equilibrium contact state after correction, reflecting the physical effect of thermal softening causing an irreversible increase in the contact area.

[0169] According to another aspect of this application, a multiphysics coupling implementation process is also proposed, specifically, realizing the transfer and iterative solution of multi-field data (electric-thermal-mechanical-structural) within the finite element framework, including:

[0170] Step 501: Based on the finite element model and material property parameters, perform steady-state electro-thermal coupling analysis to obtain the initial temperature distribution field, including:

[0171] Based on the rated current operating condition, the current continuity equation is solved to obtain the current density distribution and Joule heat source density of the contact system.

[0172] By applying the Joule heat source density to the heat conduction equation, the temperature distribution can be obtained by solving the equation.

[0173] Update the resistivity and thermal conductivity properties of the material based on the temperature distribution;

[0174] The temperature distribution solution and update process is repeated in a bidirectional coupled iteration until the temperature distribution converges. The converged temperature distribution is then used as the initial temperature distribution field. In other words, the electrical conduction equation and the heat conduction equation are iteratively solved until the temperature field distribution converges, and the converged temperature field is used as the initial temperature distribution field.

[0175] Specifically, in the steady-state analysis phase, the system needs to satisfy current conservation and energy conservation, and its governing equations are:

[0176] Electric conduction equation (current conservation): div(σ(T)×grad(V))=0;

[0177] Where div is the divergence operator, grad is the gradient operator, V is the electric potential, and σ(T) is the conductivity as a function of temperature.

[0178] Heat conduction equation (energy conservation): div(k(T)×grad(T))+Q _joule =0;

[0179] Where k(T) is the thermal conductivity as a function of temperature, T is the temperature, and Q is the thermal conductivity. _jouleThis is the Joule heat source term. The Joule heat source term is calculated from the electric field, and its calculation formula is:

[0180] Q _joule =σ(T)×|grad(V)| 2 ;

[0181] The equation is passed through conductivity σ(T) and Joule heat Q. _joule Tightly coupled.

[0182] Accordingly, the iterative solution process includes:

[0183] Assuming the entire model is at room temperature, for example 20°C, and based on the material properties σ(T) at room temperature... _0 Solving the electrical conduction equation yields the preliminary potential distribution V. _0 Joule heat distribution Q _0 .

[0184] Next, Q _0 The load is applied as a heat source to the heat conduction equation, and a new temperature distribution T is obtained by solving for it. _1 .

[0185] Based on this, according to the new temperature T _1 The electrical conductivity σ(T) of the updated material _1 ), and by resolving the electrical conduction equation, V is obtained. _1 and Q _1 .

[0186] This process is repeated until the norm of the temperature difference between two adjacent iterations is less than a preset tolerance, such as 0.1℃. The output is the convergent temperature field T. _final This is the initial state for subsequent transient analysis.

[0187] Step 502, perform transient multiphysics coupling analysis under short-circuit conditions based on the finite element model, including:

[0188] At each time step, the Lorentz force on the contact system is calculated based on the short-circuit current curve;

[0189] By applying the Lorentz force to the structural dynamics equations, the nodal displacements of the contact system and the contact forces at the contact interfaces can be obtained.

[0190] The process involves determining the contact separation state of the contact system based on nodal displacement and calculating the evolution of the contact force transmitted to the contact interface. In other words, at each time step, the Lorentz force obtained from electromagnetic field analysis is mapped to the moving contact, and the structural dynamics equations are solved to obtain the contact displacement and contact force.

[0191] Specifically, upon entering transient analysis, the time dimension and inertial effects must be introduced. The structural dynamics equations take the following form:

[0192] M×(d) 2 u) / (dt) 2 )+C×(du) / (dt)+K×u=F _ext ;

[0193] Where M is the mass matrix, C is the damping matrix, K is the stiffness matrix, u is the displacement vector, and F... _ext Let be the external force vector, (d) 2 u) / (dt) 2 (du) / (dt) is the second derivative of displacement with respect to time, i.e., acceleration; (du) / (dt) is the first derivative of displacement with respect to time, i.e., velocity.

[0194] Correspondingly, the external force F _ext The Lorentz force F generated by the short-circuit current _Lorentz The calculation of the Lorentz force is based on Ampere's law:

[0195] F _Lorentz =∫(J×B)dV;

[0196] Where J is the current density vector, B is the magnetic induction intensity vector, ∫ represents the integral, and dV is the volume element.

[0197] Specifically, a sequential coupling strategy is adopted for coupling implementation, that is, at each time step t _i Initially, the current distribution is calculated based on the contact state at the previous moment, and the electromagnetic force is obtained by solving Maxwell's equations. Further, this electromagnetic force is applied as a nodal load to the structural mesh. The structural dynamics equations are then solved using time integration algorithms such as the Newmark-β method or the Hilber-Hughes-Taylor (HHT) method, updating the displacement u(t) of the moving contact. _i ) and contact pressure F at the contact interface _c (t _i ).

[0198] Among them, the magnetic induction intensity B can be obtained by solving Ampere's law in Maxwell's equations from the short-circuit current distribution. This solution process can be automatically completed using the electromagnetic field solver in the finite element model, which is a conventional numerical calculation method for electromagnetic fields in this field.

[0199] Step 503, the transient multiphysics coupling analysis also includes a thermal feedback update step: using the calculated dynamic contact resistance and the short-circuit current flowing through the contact system, the Joule heat source power at the contact interface is updated.

[0200] Joule thermal power P at the contact interface _contact The calculation formula is:

[0201] P _contact =I _short2 ×R _c ;

[0202] Among them, I _short R is the instantaneous short-circuit current flowing through this contact point. _c This represents the dynamic contact resistance. In the finite element model, this power is treated as a surface heat source applied to the surface of the contact element, or as a volume heat source applied inside the contact element. For example, the contact element volume is V. ol Then the rate of body heat generation is q=P _contact / V ol Among them, the dynamic contact resistance R_c directly determines the concentrated heat generation power at the contact interface.

[0203] Step 504: Apply the Joule heat source power to the transient heat conduction equation and solve for the contact interface temperature in the next time step.

[0204] Specifically, the transient heat conduction equation differs from the steady-state equation in that it includes a heat capacity term and can describe the change in temperature over time:

[0205] ρ×C _p ×(dT / dt)=div(k×grad(T))+Q _joule +Q _contact ;

[0206] Where ρ is density, C _p For specific heat capacity, Q _contact The calculated heat source at the contact surface is dT / dt, which is the first derivative of temperature T with respect to time t.

[0207] Solving this equation yields the next time step t. _(i+1) Temperature field distribution T(t) _(i+1) This process includes the heat generated by volume resistance and the contact resistance R. _c The resulting localized high temperature.

[0208] Step 505: Feedback the contact interface temperature of the next time step to the calculation process of dynamic contact resistance, forming a two-way closed-loop coupling of multiple fields of electricity, heat and force.

[0209] Specifically, the calculated new temperature T(t) _(i+1) It will be through the thermal damage function Θ _deg (T) and the temperature coefficient of resistance γ affect the contact resistance R at the next moment. _c The calculation involves considering that temperature changes alter the Young's modulus E(T) of the material. As temperature increases, the Young's modulus decreases, the material softens, and this leads to increased contact deformation and contact area under the same electromagnetic repulsion. Therefore, the stiffness matrix K(T) is updated in the structural analysis of the next time step.

[0210] Based on this, a strongly coupled loop is constructed, namely the resistor R. _c → Current I & Heat source Q → Temperature T → Material properties (E, σ) & Damage Θ → Contact state quantity ξ → Resistance R _c This fully closed-loop mechanism can capture the electro-thermal-mechanical oscillations during a short circuit.

[0211] For example, if computational resources are limited, weak coupling or unidirectional coupling may be used, such as ignoring the effect of temperature on Young's modulus. This implementation sacrifices computational accuracy, especially under extreme conditions close to the melting point.

[0212] As an example, a design application of a bump contact is also proposed, specifically including:

[0213] Step 601: The moving contact surface in the contact system is provided with several protrusions to form a multi-point contact structure;

[0214] Obtaining the finite element model of the contact system includes: calling a refined geometric model containing convex structures, and performing local mesh refinement on the convex structures and their contact areas to capture the evolution of the contact state of each convex structure during short-circuit vibration. In other words, the contact interface of the contact system has several convex structures, and the area where the convex structures are located is subjected to local mesh refinement.

[0215] Specifically, in high-performance vacuum circuit breakers or high-current contactors, a specific array of bumps is designed on the surface of the moving contact to control the trajectory of the electric arc and prevent disordered adhesion of the contact surfaces. These bumps are typically cylindrical, frustum-shaped, or spherical, with a height on the order of millimeters. For example, a single moving contact may have three or more bumps evenly distributed on its surface.

[0216] Accordingly, special meshing is required for the convex region when establishing the finite element model. Because the short-circuit current concentrates through the convex, the current density is high (e.g., thousands of amperes per square millimeter) and the contact stress gradient is large. Using a conventional mesh size will fail to capture the current contraction effect and local Joule heat concentration. Therefore, a multi-level sub-model or local partitioning technique is employed to divide the cylindrical region containing the convex and the underlying matrix region into independent, refined mesh regions.

[0217] Next, if the overall grid size of the contact is 2 mm, the grid size of the convex region is refined to 0.1 mm or smaller, so that each convex contact surface has at least 5×5 unit nodes. This refinement process allows for sufficiently dense numerical integration points to calculate the shunting ratio and contact resistance evolution at each convex.

[0218] Step 602, after obtaining the short-circuit dynamic response data of the contact system, also includes:

[0219] Extract the contact separation time and peak temperature rise data for each protrusion;

[0220] The vibration stability of the convex structure is evaluated based on the contact separation time, and its short-circuit withstand capability is evaluated based on the peak temperature rise. In other words, the separation time t of the contact system is calculated based on short-circuit dynamic response data. _sep and the highest temperature rise T _peak .

[0221] Specifically, after completing the simulation calculation, the system automatically extracts two performance evaluation indicators, with the separation time t being the first indicator. _sep and the highest temperature rise T _peak Separation time t _sep This refers to the time required from the onset of the short-circuit current, i.e., t=0, until the displacement of the moving contact first exceeds a preset separation threshold, indicating physical contact disconnection. This index represents the contact system's ability to resist electrodynamic repulsion, where t... _sep The longer the contact, the tighter the closure and the better the resistance to repulsion.

[0222] Maximum temperature rise T _peak This refers to the highest temperature value ever recorded on all contact bumps throughout the entire simulation time window. This indicator is related to the contact's resistance to welding, i.e., Tcontact. _peak The lower the value, the less likely the contacts are to weld together.

[0223] Step 603: Based on the evaluation results, set the number, distribution location, or geometry of the protrusions until the vibration resistance and temperature rise performance indicators are met. Adjust the number, location, or geometry of the protrusion structure and repeat the simulation until the optimal separation time t is found. _sep Maximum and highest temperature rise T _peak Minimal bump design scheme.

[0224] Accordingly, multiple candidate solutions are generated by changing the parameters of the convex points.

[0225] For example, a reference scheme is generated, in which the moving contact surface is flat and has no protrusions.

[0226] Furthermore, the separation time t _sep =2.5ms, highest temperature rise T _peak =600℃. This design has a low temperature rise, but the contact is unstable under certain vibration conditions.

[0227] For example, a 3-bump scheme can be generated by evenly distributing three circular bumps with a diameter of 5 mm along the circumference.

[0228] Specifically, due to the reduction in contact points, the contact pressure increases, the contact becomes tighter, and the contact resistance decreases. Separation time t _sep The extension increases the repulsive force, but due to the concentrated current, the maximum temperature rises by T._peak rise.

[0229] For example, six small bumps with a diameter of 3 millimeters are distributed along the circumference, and the height of the bumps is increased.

[0230] Next, the separation time t _sep Maintaining a duration of 3.1ms, it exhibits good repulsion. However, due to the increased number of current shunt paths, the maximum temperature rises by T. _peak The temperature is reduced to 720°C. Therefore, this setting maintains the anti-repulsion performance and controls the temperature rise.

[0231] As another alternative implementation method, besides manually adjusting the parameters, the simulation program can be integrated into an automation platform, such as Easet, with the objective function set as follows:

[0232] J=w1×(1 / t _sep )+w2×T _peak ;

[0233] Where w1 is the weighting coefficient for the separation time term, and w2 is the weighting coefficient for the maximum temperature rise term; w1 and w2 contain corresponding dimension transformation operators used to map the time and temperature terms to dimensionless penalty scores. The globally optimal convex point layout parameters are automatically found using a genetic algorithm or particle swarm optimization algorithm.

[0234] In a contact system containing multiple bumps, each bump maintains its own independent contact state variable ξ. _i The subscript i represents the convex point number, based on the local contact force F it bears. _c，_i and local temperature T _c，_i The evolution equations are solved. The coupling between the convex points arises naturally from the global field equations of the finite element model: the resistance R of each convex point. _c，_i Each bump receives its own current, which in turn affects the Joule heat and temperature of each bump; correspondingly, each bump shares the total contact force, and the bumps soften due to heat, leading to ξ. _i As the force increases, the force shared by the convex point increases, while the force on the other convex points decreases accordingly. The indirect coupling achieved by this global field equation allows the simulation to automatically capture the load transfer and heat concentration effects between multiple convex points.

[0235] According to another aspect of this application, a contact system simulation system, device, and storage medium are proposed, comprising:

[0236] A simulation system for temperature rise and short-circuit dynamic response of a contact system based on a contact state evolution model includes: a model acquisition module for acquiring the finite element model and material property parameters of the contact system; and a steady-state analysis module for performing steady-state electro-thermal coupling analysis based on the finite element model and material property parameters to obtain the initial temperature distribution field.

[0237] Specifically, the model acquisition module is a graphical user interface (GUI) front-end that supports users importing CAD geometric files in formats such as STEP and IGES, and provides a material database interface for users to select material parameters such as copper, silver, or copper-chromium alloys. This module also integrates a mesh generation engine that can automatically generate suitable finite element meshes.

[0238] The steady-state analysis module encapsulates numerical algorithms for solving elliptic partial differential equations, such as the preprocessed conjugate gradient method (PCG), for quickly calculating the potential and temperature distribution under steady-state conditions.

[0239] The transient coupling module is used to perform transient multiphysics coupling analysis under short-circuit conditions based on the finite element model, using the initial temperature distribution field as the initial condition.

[0240] The transient coupling module includes a time integral controller for managing the simulation time step. This module can read the temperature field data output by the steady-state analysis module and map it to step 0 of the transient analysis. This module has a multi-physics interface, which can coordinate the use of the electromagnetic field solver, thermal field solver, and structural dynamics solver to achieve real-time exchange of data between fields.

[0241] The contact evolution calculation module is used to calculate the dynamic contact resistance at each time step of the transient multiphysics coupling analysis, based on the current contact force and the evolution state of the contact interface. The calculation of dynamic contact resistance is related to the dynamic response characteristics of the contact interface under force-thermal load.

[0242] The contact evolution calculation module incorporates a contact evolution model algorithm, including a differential equation solver, asymmetric hysteresis logic, and a thermal damage correction function. This module runs as a user subroutine or a standalone plugin, intercepting the solver's contact force data at each time step. After internal calculations, it returns the corrected contact resistance matrix or heat source vector to the main solver, achieving a seamless integration of microscopic contact physics and macroscopic finite element calculation.

[0243] The response analysis module is used to obtain short-circuit dynamic response data of the contact system based on the results of transient multiphysics coupling analysis.

[0244] Specifically, the response analysis module is responsible for data post-processing and visualization, converting massive amounts of node data (such as temperature, displacement, and voltage) into intuitive graphs, contour maps, or animations. This module has preset dedicated report templates for the contact system, automatically calculating and displaying indicators such as separation time, peak temperature, and contact bounce count.

[0245] This application implements the simulation method in the form of hardware carrier and software product, and is an assembled physical device based on system architecture, hardware configuration and computer program product.

[0246] For example, an electronic device is provided, including a memory and a processor. The memory stores a computer program, and the processor executes the computer program to implement a simulation method for dynamic response of contact system temperature rise and short circuit.

[0247] This invention also relates to a physical electronic device for running simulation software. This electronic device is typically a high-performance workstation or server. Considering the computational resource requirements of multiphysics coupling calculations, the device is equipped with a multi-core processor with a clock speed of 3.0 GHz or higher, 16 or more cores, and supports parallel computing.

[0248] The storage system includes not only non-volatile storage media, such as solid-state drives (SSDs), with a capacity of 2TB or more for storing the operating system and simulation software, but also high-speed random access memory (RAM) with a capacity of 128GB or more for loading stiffness matrices and solving temporary data. Furthermore, to accelerate graphics rendering and some matrix operations, the device is also equipped with a professional-grade graphics card.

[0249] For example, a computer-readable storage medium is provided, on which a computer program is stored, which, when executed by a processor, implements a simulation method for dynamic response of contact system temperature rise and short circuit.

[0250] This application is delivered in the form of a software product and stored on a computer-readable storage medium. The storage medium may be a physical CD-ROM / DVD, USB flash drive, or external hard drive, or it may be network-based cloud storage space or a download link.

[0251] The computer program contains all the source code or compiled binary code that implements the simulation methods described above. For example, it may contain a core algorithm library written in C++ or Fortran for solving differential evolution equations, and user interface logic written in Python. When the program is installed and loaded into computer memory, it controls the computer's processor to execute a sequence of instructions, including model building, mesh generation, equation assembly, iterative solving, and result post-processing, thus reproducing the physical behavior of the contact system in the digital world.

[0252] It should be noted that the various specific technical features described in the above embodiments can be combined in any suitable manner without contradiction. To avoid unnecessary repetition, the present invention will not describe the various possible combinations separately.

Claims

1. A simulation method for temperature rise and short-circuit dynamic response of a contact system, characterized in that, include: Obtain the finite element model and material property parameters of the contact system; Steady-state electro-thermal coupling analysis was performed based on the finite element model and material property parameters to obtain the initial temperature distribution field. Using the initial temperature distribution field as the initial condition, transient multiphysics coupling analysis under short-circuit conditions is carried out based on the finite element model. In each time step of the transient multiphysics coupling analysis, the dynamic contact resistance is calculated based on the contact force and the evolution state of the contact interface at the current moment. The calculation of the dynamic contact resistance is based on the dynamic response characteristics of the contact interface under force-thermal load. Based on the results of transient multiphysics coupling analysis, the short-circuit dynamic response data of the contact system are obtained.

2. The method according to claim 1, characterized in that, The dynamic contact resistance is calculated based on the current contact force and the evolution of the contact interface, specifically including: Introduce a contact state variable to characterize the normalized effective contact area of ​​the contact interface; Determine the equilibrium contact state of the contact interface based on the contact force at the current moment; Invoke the differential evolution equation of the contact state variables with respect to time, and update the contact state variables for the next time step based on the differential evolution equation; Calculate the dynamic contact resistance based on the updated contact state variables and the current temperature of the contact interface; The differential evolution equation characterizes the dynamic process by which the contact state variables gradually change towards the equilibrium contact state over time.

3. The method according to claim 2, characterized in that, The differential evolution equation contains growth and decay terms. The growth time constant corresponding to the growth term is smaller than the decay time constant corresponding to the decay term, which characterizes the hysteresis characteristics of the contact interface during loading and unloading. The differential evolution equation is expressed as: dξ / dt=H(ξ _eq -ξ)×((ξ _eq -ξ) / τ _inc )+H(ξ-ξ _eq )×((ξ _eq -ξ) / τ _dec ); Where dξ / dt is the rate of change of the contact state variable with time, t is time; H(...) is the unit step function; ξ is the current contact state variable; _eq The equilibrium contact state is calculated based on the current contact force; τ _inc τ is the growth time constant, representing the rate at which the contact area increases; _dec τ is the decay time constant, characterizing the rate at which the contact area decreases, and τ _inc <τ _dec .

4. The method according to claim 2, characterized in that, The dynamic contact resistance is calculated based on the updated contact state variable ξ and the current temperature of the contact interface, specifically including: Determine whether the contact state variable ξ is less than the tunneling effect activation threshold; If so, then a tunneling effect resistance term is introduced to correct the contact resistance, and the dynamic contact resistance R _c The calculation is as follows: R _c =R _c0 ×(1+γ×(T _c -T _0 ))×ξ -α +R _tunnel ; R _tunnel =A×exp(B×(1-ξ / ξ _tunnel )); Among them, R _c0 The reference contact resistance is α; the contact resistance index is R. _tunnel For tunneling effect resistance; ξ _tunnel The tunneling effect activation threshold is defined as the threshold at which the contact state variable ξ is less than, indicating that the contact interface is in the microscopic separation transition stage; A and B are tunneling effect constants; γ is the temperature coefficient of resistance of the material; T _c T represents the current temperature of the contact interface. _0 This is a reference temperature.

5. The method according to claim 1, characterized in that, Steady-state electro-thermal coupling analysis was performed based on the finite element model and material property parameters to obtain the initial temperature distribution field, including: Based on the rated current operating condition, the current continuity equation is solved to obtain the current density distribution and Joule heat source density of the contact system. By applying the Joule heat source density to the heat conduction equation, the temperature distribution can be obtained by solving the equation. Update the resistivity and thermal conductivity properties of the material based on the temperature distribution; Repeat the steps to perform bidirectional coupling iterations until the temperature distribution converges, and use the converged temperature distribution as the initial temperature distribution field.

6. The method according to claim 1, characterized in that, Transient multiphysics coupling analysis under short-circuit conditions based on the finite element model includes: At each time step, the Lorentz force on the contact system is calculated based on the short-circuit current curve; By applying the Lorentz force to the structural dynamics equations, the nodal displacements of the contact system and the contact forces at the contact interfaces can be obtained. The process of calculating the evolution of the contact force transmitted to the contact interface is based on determining the contact separation state of the contact system according to the nodal displacement.

7. The method according to claim 2, characterized in that, The contact state variable ξ for the next time step is updated based on the differential evolution equation, specifically solved using an exponential time integral scheme: ξ(t+Δt)=ξ _eq +(ξ(t)-ξ _eq )×exp(-Δt / τ); Where t is the current time, Δt is the time step, and τ is the growth time constant or decay time constant selected based on the current evolution trend.

8. The method according to claim 1, characterized in that, The moving contact surface in the contact system is provided with several protrusions to form a multi-point contact structure; Obtaining the finite element model of the contact system includes: calling a refined geometric model containing the convex structure, performing local mesh refinement on the convex and its contact area, and capturing the evolution behavior of the contact state of each convex during short-circuit vibration.

9. The method according to claim 8, characterized in that, After obtaining the short-circuit dynamic response data of the contact system, the following is also included: Extract the contact separation time and peak temperature rise data for each protrusion; The vibration resistance of the bump structure is evaluated based on the contact separation time, and the short-circuit withstand capability of the bump structure is evaluated based on the peak temperature rise. Based on the evaluation results, the number, distribution location, or geometric dimensions of the protrusions are set until the vibration resistance and temperature rise performance indicators are met.

10. A simulation system, said system being used to implement the simulation method for temperature rise and short-circuit dynamic response of a contact system as described in any one of claims 1 to 9, characterized in that, The system includes: The model acquisition module is used to acquire the finite element model and material property parameters of the contact system; The steady-state analysis module is used to perform steady-state electro-thermal coupling analysis based on the finite element model and material property parameters to obtain the initial temperature distribution field. The transient coupling module is used to perform transient multiphysics coupling analysis under short-circuit conditions based on the finite element model, using the initial temperature distribution field as the initial condition. The contact evolution calculation module is used to calculate the dynamic contact resistance at each time step of the transient multiphysics coupling analysis based on the current contact force and the evolution state of the contact interface. The calculation of the dynamic contact resistance is related to the dynamic response characteristics of the contact interface under force-thermal load. The response analysis module is used to obtain the short-circuit dynamic response data of the contact system based on the results of the transient multiphysics coupling analysis.

Images