Suspension bridge main cable nonlinear numerical simulation modeling method

Abstract

The application relates to the field of bridge engineering and computational mechanics, and discloses a nonlinear numerical simulation modeling method for a main cable of a suspension bridge. The method comprises the following steps: a finite element model of the main cable of the suspension bridge is established, and a group of internal state variables are introduced on integral points of cable elements of the finite element model, the internal state variables comprising a section topological configuration field, a friction interface state memory field, a slip activation field and a slip kinematics field; in an incremental step of finite element analysis, according to a preset coupling evolution control equation, the update values of the internal state variables driven by total strain increments are solved; then, the real stress tensor and the tangent stiffness matrix are calculated based on the updated internal state variables, and the calculation is returned to a main analysis program for balance iteration; through the calculation of each incremental step, the tracking and accumulation of the history-dependent response are realized. The application can capture the complex phenomena such as wear, slip and accumulation configuration rearrangement between steel wires, and improves the physical fidelity and calculation stability of numerical simulation.

Description

Technical Field

[0001] This invention relates to the fields of bridge engineering and computational mechanics, specifically to a nonlinear numerical simulation modeling method for the main cable of a suspension bridge. Background Technology

[0002] The main cable of a suspension bridge is the most critical load-bearing component. The accurate simulation of the mechanical behavior of the main cable is directly related to the accuracy of the assessment of the safety, applicability and durability of the entire bridge structure under loads such as wind, earthquake and vehicles. The main cable is made of thousands or even tens of thousands of high-strength steel wires tightly compressed together. The internal structure of the main cable is extremely complex, which makes the macroscopic mechanical response of the main cable exhibit significant nonlinearity, anisotropy and history dependence.

[0003] In existing numerical simulation practices, to simplify calculations, the main cable is often treated as an equivalent truss or cable element, and its constitutive relationship is assumed to be linear elastic. This approach completely ignores the complex physical mechanisms such as contact, friction, and relative slippage between the steel wires inside the main cable, and therefore fails to capture important characteristics such as energy dissipation (i.e., damping) and stiffness nonlinearity generated by these mechanisms, resulting in significant deviations in the prediction of the bridge structure's dynamic response and long-term deformation.

[0004] To improve simulation accuracy, some studies have begun to incorporate the frictional effects between steel wires into the models, typically using classical elastoplastic theory or specific friction models. However, these improved models still significantly simplify the description of physical reality. They generally fail to consider the historical evolution of the contact section under long-term cyclic loading; that is, the interface polishing and frictional performance degradation effects caused by fretting wear are generally ignored. Furthermore, the actual stacking pattern of steel wires within the main cable section is not static; local rearrangement may occur under severe vibration, and the impact of this change in topology on the mechanical properties of the main cable is often overlooked by existing models. Simultaneously, the transition from a globally locked state to a fully slipped state is a gradual process, while existing models often treat it as an abrupt event, failing to accurately reflect the gradual activation characteristics of slip.

[0005] Therefore, this invention proposes a nonlinear numerical simulation modeling method for the main cable of a suspension bridge to address the shortcomings of existing technologies. Summary of the Invention

[0006] To address the shortcomings of existing technologies, this invention provides a nonlinear numerical simulation modeling method for the main cable of a suspension bridge. This method solves the problem that numerical simulation methods for the main cable of suspension bridges are usually based on simplified linear elastic assumptions or only consider a single frictional slip effect. These methods cannot accurately and comprehensively capture the nonlinear mechanical behavior of the main cable under complex loads, which is caused by the interaction of thousands of steel wires within it, resulting in multi-physics coupling and history dependence. These complex nonlinear mechanical behaviors include wear degradation between steel wires, rearrangement of stacking configurations, and a gradual process from viscous locking to macroscopic slip. Their combined effects significantly affect the stiffness, damping, and long-term performance of the main cable.

[0007] To achieve the above objectives, the present invention provides a nonlinear numerical simulation modeling method for the main cable of a suspension bridge, comprising the following steps:

[0008] Step S1: Establish the initial finite element model of the main cable of the suspension bridge, and initialize a set of internal state variables at the integration point of each cable element in the finite element model, wherein the internal state variables include:

[0009] The cross-sectional topology field characterizes the wire stacking pattern within the cross section, the friction interface state memory field characterizes the cumulative wear effect between the wires, the slip activation field characterizes the degree of transition from the locked state to the slip state of the cross section, and the slip kinematic field quantifies the irreversible deformation caused by the relative slip between the wires.

[0010] Step S2: In any incremental step of the finite element analysis, obtain the total strain increment corresponding to the incremental step;

[0011] Step S3: Based on the preset, mutually coupled evolution control equations, solve for the updated values ​​of each internal state variable driven by the total strain increment;

[0012] Step S4: Based on the sliding kinematic field updated by one of the internal state variables, calculate the true stress tensor and the tangent stiffness matrix corresponding to the current state at the integration point at the end of the current increment step.

[0013] Step S5: Transfer the calculated true stress tensor and tangent stiffness matrix to the master finite element analysis program;

[0014] Step S6: After the calculation of one incremental step is completed, the updated internal state variables obtained in the incremental step are used as the initial conditions for the next incremental step, and steps S2 to S5 are repeated until all incremental steps are completed, thus completing the nonlinear numerical simulation modeling of the finite element model of the suspension bridge main cable.

[0015] Preferably, after step S2 and before step S3, the method further includes a step of performing elasticity prediction calculation:

[0016] Based on the assumption that the sliding kinematic field remains unchanged, a trial stress tensor is calculated according to the converged state at the end of the previous increment step and the total strain increment of the current increment step, and the trial stress tensor is used as the mechanical driving force for the evolution of the internal state variables.

[0017] Preferably, in step S3, the mutually coupled evolutionary governing equations include solving the evolution of the cross-sectional topological configuration field, specifically:

[0018] Calculate the cumulative slip norm obtained by integrating the evolution rate of the slip kinematic field over time;

[0019] The cumulative slip norm is compared with a preset critical threshold that depends on the current topology. If the critical threshold is exceeded, the transition probability from the current topology to each target topology is calculated based on the potential energy function associated with each topology.

[0020] The cross-sectional topology field is determined and updated based on the aforementioned transition probability distribution.

[0021] Preferably, in step S3, the mutually coupled evolution control equations include solving the coupling relationship between the friction interface state memory field and the slip activation field, specifically including:

[0022] Solve for the degradation evolution of the state memory field of the friction interface driven by frictional dissipation energy;

[0023] Solve for the activation evolution of the slip activation field, wherein the degenerate evolution updated friction interface state memory field, by influencing the calculation of slip driving energy, in turn affects the evolution of the slip activation field.

[0024] Preferably, the activation evolution of the slip activation field is controlled by the phase field equation, wherein:

[0025] The driving force for the activation evolution originates from the reduction of the pseudo-free energy functional;

[0026] The pseudo-free energy functional includes a slip driving energy term, the value of which depends on the stress state and the critical stress threshold.

[0027] The critical stress threshold is set based on the friction interface state memory field.

[0028] Preferably, in step S3, the mutually coupled evolution control equations include solving the nonlinear modulation coupling of the slip activation field to the slip kinematic field, specifically including:

[0029] An evolutionary law for the sliding kinematic field is established, which stipulates that the evolution rate of the sliding kinematic field depends on the sliding rate multiplier.

[0030] The value of the slip rate multiplier is calculated from the value of the slip activation field using a nonlinear modulation function.

[0031] Preferably, the evolution law of the slip kinematic field is a flow law, wherein the evolution rate of the slip kinematic field is equal to the product of the slip rate multiplier and the gradient of the slip yield function with respect to the true stress tensor.

[0032] Preferably, in step S3, solving the mutually coupled evolution control equations involves executing a return mapping algorithm; the return mapping algorithm includes discretizing the first-order ordinary differential equations describing the evolution of each internal state variable through an implicit time integration scheme to form a nonlinear algebraic equation system.

[0033] Preferably, the system of nonlinear algebraic equations is solved using the Newton-Raphson iterative method, which includes the following steps:

[0034] In each iteration, a residual vector consisting of the constraints of each internal state variable is constructed, and the Jacobian matrix is ​​calculated based on the constructed residual vector. The Jacobian matrix is ​​obtained by taking the partial derivatives of the residual vector with respect to all the internal state variables to be solved.

[0035] Solve the system of linear equations consisting of the Jacobian matrix and the residual vector to obtain the incremental correction vector of the internal state variables. Use the incremental correction vector to update the internal state variables until convergence.

[0036] In step S4, the tangent stiffness matrix corresponding to the current updated state is calculated. The tangent stiffness matrix is ​​a consistent tangent stiffness matrix. Furthermore, the calculation of the consistent tangent stiffness matrix uses the Jacobian matrix obtained in the Newton-Raphson iteration method.

[0037] This invention provides a nonlinear numerical simulation modeling method for the main cable of a suspension bridge. It has the following beneficial effects:

[0038] 1. This invention introduces four internal state variables into the finite element model of the main cable of a suspension bridge: the cross-sectional topological configuration field, the friction interface state memory field, the slip activation field, and the slip kinematic field. This allows the finite element model to move beyond the assumptions of ideal elasticity or single friction and comprehensively and precisely describe the complex mechanical behavior of the main cable, which is determined by the evolution of the steel wire stacking pattern, the wear and degradation of the contact interface, and the gradual process from local viscosity to overall slip. As a result, the simulation results are closer to physical reality and can provide a more reliable basis for the refined design and long-term performance evaluation of suspension bridges.

[0039] 2. This invention establishes evolution control equations that are mutually coupled among the internal state variables. In particular, it regulates the slip activation critical threshold through the friction interface state memory field and nonlinearly modulates the evolution rate of the slip kinematic field through the slip activation field. This makes the current response of the finite element model closely related to the entire loading history. For example, the cumulative wear caused by historical vibration loads will lower the threshold for future slip. This capture of path dependence effect enables this method to truly reflect the key performance changes of the main cable, such as stiffness degradation and damping evolution, caused by reciprocating loads during long-term service.

[0040] 3. The numerical implementation scheme adopted in this invention, in particular, solves the coupled evolution equations through the return mapping algorithm and solves the discrete nonlinear algebraic equations through Newton-Raphson iteration, ensuring the robustness of the algorithm. More importantly, by calculating the consistent tangent stiffness matrix that is completely consistent with the constitutive integral algorithm, the quadratic convergence rate of the global finite element analysis is guaranteed. This significantly reduces the number of iterations required for each incremental step when dealing with highly nonlinear problems, greatly shortening the computation time and making it computationally possible to perform long-term and complex nonlinear time history analysis on the full-bridge model. Attached Figure Description

[0041] Figure 1 This is a flowchart of the nonlinear numerical simulation modeling method for the main cable of a suspension bridge according to the present invention;

[0042] Figure 2 This is a flowchart of the coupled evolution solution within a single incremental step according to the present invention;

[0043] Figure 3 This is a schematic diagram of the system architecture of the present invention;

[0044] Figure 4 This is a schematic diagram of the structure of the method of the present invention that incorporates the global loop of finite element analysis;

[0045] Figure 5 This is a schematic diagram illustrating the probabilistic transformation of the cross-sectional topological configuration field of the present invention;

[0046] Figure 6 This is a schematic diagram illustrating the modulation effect of the slip activation field on the slip rate according to the present invention.

[0047] The attached figures are labeled as follows: 10, structural discretization module; 20, internal state field definition module; 30, coupled evolution solution module; 40, mechanical response update module. Detailed Implementation

[0048] The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. 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 are within the scope of protection of the present invention.

[0049] See attached document Figure 3 This invention provides a modeling system for implementing the nonlinear numerical simulation modeling method for the main cable of a suspension bridge. The modeling system runs on a computer device, which includes hardware such as a processor, memory, and storage. The system provided by this invention includes: a structure discretization module 10, an internal state field definition module 20, a coupled evolution solution module 30, and a mechanical response update module 40.

[0050] The function of the structure discretization module 10 is to establish a one-dimensional spatial path in the computer based on the geometric parameters of the main cable of the suspension bridge to be analyzed, and to discretize the spatial path into a series of one-dimensional cable elements connected end to end, forming a finite element mesh for subsequent calculations.

[0051] The function of the internal state field definition module 20 is to define a set of preset internal state variables at the integration point of each cable element generated by the structure discretization module 10 to describe the internal physical state of the cross section. The set of internal state variables specifically includes: the cross section topology configuration field characterizing the wire stacking mode within the cross section; the friction interface state memory field characterizing the cumulative wear effect between the wires; the slip activation field characterizing the degree of transition of the cross section from the locked state to the slip state; and the slip kinematic field quantifying the irreversible deformation caused by the relative slip between the wires.

[0052] The function of the coupled evolution solution module 30 is to solve the changes of each internal state variable defined by the internal state field definition module 20 in the current analysis step based on the preset multiphysics field coupled evolution control equation and the macroscopic mechanical driving quantity (such as the total strain increment) transmitted by the external finite element analysis program. The solution process calculates the updated values ​​of the cross-sectional topology configuration field, friction interface state memory field, slip activation field and slip kinematic field under the current mechanical driving.

[0053] The function of the mechanical response update module 40 is to calculate the true stress tensor of the cable element integration point at the end of the current analysis step based on the updated value of the slip kinematic field calculated by the coupled evolution solution module 30, and further calculate the tangent stiffness matrix corresponding to the current state. The calculated true stress tensor and tangent stiffness matrix are then returned to the external finite element analysis program for the equilibrium iteration of the global structure and the calculation of subsequent analysis steps.

[0054] See attached document Figure 1 With appendix Figure 4 This invention provides a nonlinear numerical simulation modeling method for the main cable of a suspension bridge, which is executed in a computer device and includes the following steps:

[0055] Step S101: Establish a geometrically nonlinear finite element model of the main cable, and initialize a set of internal state variables at the integration point of each cable element in the finite element model. The internal state variables include the cross-sectional topology configuration field, the friction interface state memory field, the slip activation field, and the slip kinematic field.

[0056] Step S102: At the start of any incremental step in the finite element analysis, obtain the total strain increment corresponding to that incremental step from the main analysis program as the mechanical input quantity driving the evolution of the internal state variables;

[0057] Step S103: Based on the preset, mutually coupled evolution control equations, solve for the updated values ​​of each internal state variable driven by the total strain increment; the process of solving for these updated values ​​involves coupled calculations of the discrete jump of the cross-sectional topological configuration field, the slow degradation of the friction interface state memory field, the rapid activation of the slip activation field, and the inelastic growth of the slip kinematic field.

[0058] Step S104: Based on the updated sliding kinematics obtained from the solution, calculate the true stress tensor at the integration point at the end of the current increment step, and calculate the tangent stiffness matrix corresponding to the current updated state;

[0059] Step S105: Return the calculated true stress tensor and tangent stiffness matrix to the main analysis program for global structural equilibrium iteration;

[0060] Step S106: Repeat steps S102 to S105 until all load increment steps or time increment steps are calculated, thereby obtaining the response of the main cable throughout the analysis process.

[0061] The technical details involved in the above steps will be explained in detail below.

[0062] See attached document Figure 5 To be continued Figure 6 In step S101, a geometrically nonlinear finite element model of the main cable is established, and a set of internal state variables is initialized at each integration point of the cable element in the finite element model. The specific implementation of this step is as follows:

[0063] S101-1: In a computer-aided engineering software environment, a macroscopic one-dimensional finite element model is established based on the design geometric path, cross-sectional properties, and material parameters of the main cable to be analyzed. This finite element model consists of a series of geometrically nonlinear cable elements or beam elements capable of handling large displacement and large rotation effects. Each cable element or beam element contains one or more integration points, which serve as the carrier for all subsequent physical quantity calculations.

[0064] S101-2: At each of the aforementioned integration points, define and initialize a set of internal state variables; this set of internal state variables constitutes a state vector S, which completely describes the internal physical state of the main cable section at that integration point; the state vector S includes:

[0065] S={Θ,μ s ,φ,ε s};

[0066] in:

[0067] Θ is the cross-sectional topological configuration field, a discrete integer scalar. Θ is used to characterize the microscopic packing pattern of thousands of steel wires inside the cross-section. As a functional generalized technical feature, its specific implementation may include, but is not limited to: setting Θ=1 to represent an ideal hexagonal close-packed configuration; setting Θ=2 to represent a local square packing configuration caused by construction or vibration; setting Θ=3 to represent a loose packing configuration with local voids. The introduction of this cross-sectional topological configuration field provides a physical basis for the subsequent definition of anisotropic slip laws and wear rates related to topological surveying.

[0068] μ s The friction interface state memory field is a continuous scalar; μ s Used to quantify the long-term wear state of the contact interface between steel wires, μ s It is proportional to the static friction threshold required to activate slip; the friction interface state memory field is a slow variable that records the cumulative degradation history. Its physical meaning is that as the energy dissipation caused by cyclic slip accumulates, the steel wire surface is worn and polished, making it easier to activate slip in the future.

[0069] φ is the slip activation field, a continuous scalar with a value in the range of [0, 1]. φ acts as a switch variable that controls the slip process, used to smoothly characterize the transition state of the cross section from fully locked to fully slipped. Specifically, φ=0 corresponds to the fully locked state in which the friction between the steel wires in the cross section is sufficient to lock all the steel wires and no relative slip occurs; φ=1 corresponds to the fully slipped state in which the static friction between the steel wires in the cross section has been overcome and is in a state of sufficient relative sliding; while 0<φ<1 corresponds to the transition state in which some steel wires begin to slip or the entire cross section is in a stick-slip transition.

[0070] ε sFor the sliding kinematic field, ε is a second-order tensor; s As part of the internal state variables, it is directly used to quantify the inelastic strain caused by the irreversible relative slip between the wires; the slip kinematic field is the inelastic component of the total strain tensor, and the evolution of the slip kinematic field represents the permanent deformation of the cross section due to internal slip.

[0071] S101-3: Initial conditions are set for the aforementioned internal state variables; at the initial moment of the analysis (t=0), it is assumed that the main cable is in an initial equilibrium state under no load or only its own weight, and this initial equilibrium state is considered to be without damage, wear, or slippage; therefore, the initial value S0 of the state vector is set as follows:

[0072] S0={Θ0,μ s,0 ,0,0};

[0073] Where: Θ0 is the initial topology assumed based on the main cable design or construction process; for example, a new main cable is usually set to an ideal hexagonal close-packed configuration, i.e., Θ0 = 1; μ s,0 The initial frictional parameters are those of the brand-new, unworn steel wire contact interface; the initial value of the slip activation field φ is set to 0, indicating that the main cable is initially in a fully locked state; the slip kinematic field ε... s The initial value is set to tensor 0, indicating that the main cable does not have inelastic deformation caused by slip in the initial state.

[0074] In the finite element analysis process, the specific implementation of step S102 is as follows:

[0075] S102-1: At the beginning of each new load increment step or time increment step, for example from time t n to t n+1 In this method, at each integration point of the cable element in the finite element model of the main cable, the total strain increment tensor Δε calculated by the cable element within that increment step is received from the main analysis program; simultaneously, the tensor at which the integration point was located at the end of the previous increment step (time t) is read and loaded. n The internal state variables that have converged and been stored include the total strain tensor ε. n and sliding kinematic field ε s,n .

[0076] S102-2: Perform elastic prediction calculations; the core assumption of this calculation is that within the current increment step, the material behavior is entirely elastic, i.e., no new internal slip has occurred; based on this assumption, a tentative stress state is calculated, which will serve as the basis for determining whether inelastic evolution needs to be activated; the elastic prediction calculation further includes: first, calculating the stress state at the end of the current increment step (time t)... n+1 The total strain tensor of the trial :

[0077] ;

[0078] In the formula, ε n Let time t n The total strain tensor, where Δε is the total strain increment tensor of the current increment step.

[0079] Then, based on the aforementioned assumption of elastic behavior, that is, the sliding kinematic field remains constant during this trial step (ε) s =ε s,n ), calculate the trial stress tensor at the integration point. :

[0080] ;

[0081] In the formula: For time t n+1 The calculated test stress tensor; C el ε is the fourth-order elastic stiffness tensor of the material, representing the inherent elastic properties of the material when slippage does not occur; s,n For time t n The value of the converged slip kinematic field, i.e., the accumulated inelastic strain; in the equation of the trial stress tensor, This represents the double-dot-product or double-contraction operation, a standard operation in tensor mechanics used to describe how a fourth-order tensor acts on a second-order tensor, producing another second-order tensor. If component form (i.e., tensor exponent notation) is used, the equation for this trial stress tensor can be written more specifically as follows: ;

[0082] Where i, j, k, and l are the component indices of the tensor, and their values ​​are all in the range {1, 2, 3}. The values ​​1, 2, and 3 represent the x, y, and z coordinate axes, respectively. It is a component of the test stress tensor; (C el ) ijkl It is a component of the fourth-order elastic stiffness tensor; It is a component of the elastic strain tensor.

[0083] S102-3: The calculated test stress tensor As the core mechanical driving force for the evolution of the internal state within the current incremental step, the experimental stress tensor will be used in the subsequent step S103 to evaluate the slip activation condition, calculate the slip driving force, and then solve the evolution of all internal state variables (including the cross-sectional topology field, the friction interface state memory field, the slip activation field, and the slip kinematic field). This elastic prediction step is the starting point for all subsequent inelastic calculations.

[0084] In step S103, for the multiphysics coupled evolution solution, a core component of this step is solving the evolution of the cross-sectional topological configuration field Θ; this evolution process is constructed as a probabilistic discrete state transition event driven by accumulated internal deformation, and its specific implementation is as follows:

[0085] S103-1: Calculate the cumulative slip; within the current increment step, first calculate a scalar to quantify the degree of cumulative wire rearrangement within the cross-section, namely the cumulative slip norm p. acc The cumulative slip norm is obtained by integrating the evolution rate norm of the slip kinematic field over time.

[0086] ;

[0087] In the formula: p acc It represents the accumulated glide norm from the last change in topology to the current moment; For the sliding kinematic field ε s The first derivative with respect to time, i.e., the slip strain rate tensor; The double dot product operation represents a tensor.

[0088] S103-2: Execute state transition condition judgment; calculate the cumulative slip norm p. acc With a preset topology Θ n critical threshold p crit (Θ n The critical threshold is compared to the current topological configuration; it represents the ability of the current topological configuration to resist deformation instability. For example, more stable topological configurations (such as ideal hexagonal close packing) have a higher critical threshold.

[0089] If p acc ≤p crit (Θ n If no topological transformation is triggered within the current increment step, the cross-sectional topological field remains unchanged, i.e., Θ n+1 =Θ n And continue to solve for other internal state variables;

[0090] If p acc >p crit (Θ nIf the condition is met, the triggering condition for the state transition is determined, and the next step S103-3 is executed.

[0091] S103-3: Calculate the transition probability and select a new topological configuration; this step is a probabilistic selection process with functional generalization, and its specific implementation may include:

[0092] First, for all possible target topological configurations Θ j (For example, a hexagonal close packing Θ=1 may transform into a square packing Θ=2 or a loose packing Θ=3), calculate the change from the current topological configuration Θ. n Transformation to the target topology Θ j The transition probability P nj The transition probability is calculated based on a law similar to the energy distribution law in statistical thermodynamics:

[0093] P nj =Z -1 exp((W pot (Θ n )-W pot (Θ j )) / E eff );

[0094] In the formula: P nj To obtain from topological configuration Θ n Transformation to topological configuration Θ j The transition probability of W; pot (Θ) is a predefined potential energy function associated with each topological configuration Θ, used to represent the relative stability of different packing modes; generally, denser topological configurations have lower potential energy. eff is an effective energy parameter, related to the intensity or strain rate of the current load, used to characterize the energy input intensity driving the system to overcome the energy barrier; Z is the normalization factor, also known as the partition function, which is calculated by summing over all possible target topologies to ensure that the sum of all transformation probabilities is 1.

[0095] ;

[0096] Then, based on the calculated probability distribution {P} nj The final target topology Θ is determined using a random number selection algorithm (such as the Monte Carlo method). n+1 .

[0097] S103-4: Perform state update and reset; assign the new topology selected in step S103-3 to the current topology, i.e., Θ. n+1 =Θ j Simultaneously, the cumulative slip norm p used to trigger the judgment will be... accThe value is reset to zero; the updated topology Θ n+1 It will serve as a key input parameter in the subsequent process of solving the friction interface state memory field and the sliding kinematic field evolution, in order to update its evolution law.

[0098] In step S103, following the solution of the cross-sectional topological configuration field, a coupled evolution solution is performed. The core of this process lies in solving the history-dependent frictional interface state memory field μ. s The feedback coupling relationship between the transient slip activation field φ and the field is solved through the following sub-steps:

[0099] S103-5: Solving for the friction interface state memory field μ s The slow degradation evolution is driven by the frictional dissipation energy generated when relative slippage occurs between the wires. The specific implementation method is as follows:

[0100] First, calculate the instantaneous frictional power dissipation. The instantaneous frictional dissipation power is obtained by a double dot product operation between the stress tensor at the current integration point and the evolution rate of the slip kinematic field.

[0101] ;

[0102] In the formula: σ represents the instantaneous frictional dissipation power per unit volume; σ is the stress tensor at the current integration point. For the sliding kinematic field ε s The first derivative with respect to time, i.e., the slip strain rate tensor.

[0103] Then, based on the calculated instantaneous frictional dissipation power, the friction interface state memory field μ is updated. s Its evolution follows a first-order ordinary differential equation, which can be represented as a degenerate process:

[0104] ;

[0105] In the formula: For the state memory field μ of the friction interface s First derivative with respect to time; C μ (Θ) is a degradation rate coefficient dependent on the current cross-sectional topological configuration field Θ, used to characterize the difference in wear rate under different wire packing modes; μ s,res Let μ be a preset constant. s The residual or saturation value, i.e., the steady state after sufficient wear; H(·) is the Heaviside step function, and the introduction of H(·) ensures that only when energy dissipation occurs ( The degradation process is only activated when the time is right, thus ensuring the irreversibility of the degradation process.

[0106] S103-6: Solve for the rapid activation evolution of the slip activation field φ, and reflect its relationship with the friction interface state memory field μ. s The coupling relationship; this evolution process is governed by a phase-field equation, such as an Allen-Cahn type phase-field equation, whose driving force originates from the reduction of the system's pseudo-free energy; the evolution equation of φ is:

[0107] ;

[0108] In the formula: M is the first derivative of the slip activation field φ with respect to time. φ These are the dynamic coefficients of phase field evolution, controlling the rate of the activation process; The pseudo-free energy functional of the system, its functional derivative with respect to φ This constitutes the thermodynamic driving force of evolution.

[0109] pseudo-free energy functional The density function is derived from the chemical energy term f. chem and gradient energy term f grad Composition. Among them, the chemical energy term f chem It is the memory field μ of the interface state with friction. s The key to feedback coupling; its specific form is:

[0110] f chem (φ,σ,μ s ,Θ)=W(φ)+h(φ)G drive (σ,μ s ,Θ);

[0111] In the formula: W(φ)=C W φ 2 (1-φ) 2 Let C be a double-well function, where C W φ is a constant, and this double-well function ensures that the glide activation field is energetically stable in both the φ=0 (locked) and φ=1 (glide) states; h(φ) is a monotonically increasing interpolation function, for example, h(φ)=φ 2 (3-2φ), used to smoothly apply driving force; G drive It is the slip driving energy, h(φ), which directly establishes a connection with the stress state and the friction memory field; as a functionally generalized technical feature, its specific implementation can be a slip criterion function f. slip The slip criterion function indicates the extent to which the current stress state exceeds the threshold for slip activation:

[0112] G drive (σ,μ s ,Θ)=τcrit (μ s )-τ eq (σ,Θ);

[0113] In the formula: τ eq (σ,Θ) is an equivalent stress whose specific form depends on the current topological configuration Θ, for example, it can be the equivalent shear stress acting on a potential slip surface; τ crit (μ s The critical stress threshold required to activate slip is explicitly defined as the friction interface state memory field μ. s A function, such as a simple linear relationship τ crit (μ s )=k μ μ s , where k μ It is a scaling factor; through this setting, when μ is calculated in step S103-5 s When the critical stress threshold τ decreases due to wear, the critical stress threshold τ at this point is... crit It also decreases accordingly.

[0114] In step S103, after solving for the history-dependent variables, the evolution of instantaneous behavior is solved. The core of this process lies in solving the relationship between the slip activation field φ and the slip kinematic field ε. s Nonlinear modulation coupling; this process is achieved through the following sub-steps:

[0115] S103-7: Establishing the sliding kinematics field ε s The evolutionary law, which takes the form of a flow law similar to viscoplastic mechanics, has an evolution rate... It is given by the following formula:

[0116] ;

[0117] In the formula: For the sliding kinematic field ε s The first derivative with respect to time, i.e., the slip strain rate tensor, represents the rate of increase of inelastic deformation; ξ slip (σ,Θ) is the slip yield function, whose form depends on the current stress tensor σ and the cross-sectional topology field Θ; this slip yield function defines the stress boundary conditions for the slip to begin. It is the gradient of the slip yield function with respect to the stress tensor. The direction of slip strain development is defined, and this direction is usually orthogonal to the yield surface; It is a non-negative scalar, called the slip rate multiplier or consistency parameter. The magnitude of the slip rate multiplier directly determines the speed of slip evolution.

[0118] S103-8: Establishing the slip rate multiplier The modulation coupling relationship between the sliding activation field φ and the slip rate multiplier. This step is crucial for achieving a smooth transition from locked to activated during the slip process, directly incorporating the value of the slip activation field φ into the slip rate multiplier. In the calculation; as a functional generalization of technical features, its specific implementation can be:

[0119] ;

[0120] Where: h γ (φ) is a nonlinear modulation function of the glide activation field φ; this nonlinear modulation function is explicitly designed as a monotonically increasing function of φ, for example, a power-law form h γ (φ)=A γ φ n A γ And n are material constants; through this monotonically increasing function, the value of the slip activation field φ directly controls the slip rate multiplier. The order of magnitude; when φ=0, h γ (0)=0, making Even if the slip yield condition is met at this moment, the slip rate is zero, and the section remains locked; as φ evolves from 0 to 1, h γ The value of (φ) increases accordingly, thereby allowing the slip rate The corresponding increase enables a smooth start to the sliding process; Macaulay brackets are defined as follows: Its function is to ensure that stress only occurs when the stress state is outside the yield surface (i.e., ξ). slip When >0), the slip rate multiplier Only then can it be a positive value; otherwise, it is zero. visc is a viscosity regularization parameter used to control the sensitivity of the slip rate to the degree to which the stress exceeds the yield surface and to improve the stability of numerical calculations; m is a material constant, which is the viscosity exponent.

[0121] S103-9: Integrating the evolution equation of the slip activation field φ with the slip kinematic field ε s The evolution equations are solved in a coupled manner; in the iterative solution process of each incremental step, the evolution equations of the slip activation field defined in the aforementioned steps S103-6 (i.e., Allen-Cahn type equations) and the evolution equations of the slip kinematic field defined in this step are solved simultaneously as a whole, nonlinear set of ordinary differential equations; after solving this set of ordinary differential equations, the updated slip activation field φ is obtained. n+1 and sliding kinematic field ε s,n+1 This completes the full calculation of instantaneous behavioral evolution.

[0122] See attached document Figure 2In step S103, the coupled evolution control equations of the aforementioned internal state variables are numerically solved within one incremental step. This process constitutes a return mapping algorithm, the specific implementation of which is as follows:

[0123] S103-10: The system of first-order ordinary differential equations describing the evolution of each internal state variable is discretized using a time integration scheme; as a functional generalization technique, the specific implementation of this integration scheme can be an implicit integration scheme with good numerical stability, such as the backward Euler method; using this backward Euler method, from time t... n to t n+1 (Time step is Δt=t) n+1 -t n The evolution problem of ) is transformed into a system of nonlinear algebraic equations that require a solution.

[0124] S103-11: To solve this system of nonlinear algebraic equations, the system is constructed as a residual vector. All components of this residual vector R should be zero when the solution converges; this residual vector contains all the internal state variables to be solved at time t. n+1 Constraints to be satisfied: ;

[0125] In the formula, all variables with subscript n+1 are unknowns to be solved within the current increment step, such as ε. s,n+1 μ s,n+1 φ n+1 etc.; among which:

[0126] For the residual components of the sliding kinematic field, ensure that their evolution follows the discretized flow rules;

[0127] For the residual components of the friction interface state memory field, ensure that its evolution follows the degenerate equation after discretization;

[0128] For the residual components of the slip-activated field, ensure that their evolution follows the discretized phase field equations;

[0129] , and At time t n+1 The slip rate multipliers, memory field evolution rate, and activation field evolution rate are calculated using formulas given in the preceding steps, and they are themselves time t. n+1 Functions that handle other internal state variables.

[0130] S103-12: This nonlinear residual equation system R=0 is solved using an iterative method. As a functional generalization technique, this method can be implemented using a Newton-Raphson iterative method. This method starts with an elastic trial solution and approximates the true solution through a series of linearized iterations. In the k-th iteration, the following linear equation system is solved:

[0131] J (k) δy (k+1) =-R (k) ;

[0132] In the formula: R (k) Let δy be the residual vector calculated under the known state in the k-th iteration. (k+1) This is the increment correction vector for the internal state variables to be solved in the (k+1)th iteration; this increment correction vector contains the correction amounts for all the internal state variables to be solved, for example... J (k) The Jacobian matrix, also known as the consistent tangent operator, is calculated in the k-th iteration state. This Jacobian matrix is ​​obtained by taking the partial derivatives of the residual vector R with respect to all the internal state variables y to be solved.

[0133] ;

[0134] The analytical or numerical calculation of the Jacobian matrix is ​​crucial for the implementation of the Newton-Raphson method, and its accuracy directly affects the convergence speed of the iteration.

[0135] S103-13: After solving for the incremental correction vector δy (k+1) Then, update the internal state variables to be solved:

[0136] y (k+1) =y (k) +δy (k+1) .

[0137] S103-14, after each iteration update, calculate the current internal state variable y. (k+1) The residual vector norm ||R|| (k +1) ‖, and compare it with a preset convergence tolerance ε tol Compare; if ||R (k+1) ||<ε tol If the iteration converges, the current internal state variable y is determined to be... (k+1) This is the final solution for the current increment step, ε. s,n+1 μ s,n+1 φ n+1If the condition is met, then proceed to step S103-12 and use the updated internal state variables for the next iteration until the convergence condition is met or the preset maximum number of iterations is reached.

[0138] In step S104, based on the coupled evolution solution obtained in step S103, at time t n+1 The converged internal state variables are used to update the mechanical response; the specific implementation of this step is as follows:

[0139] S104-1: Calculate the true stress tensor at the integration point at the end of the current increment step; this calculation is based on the assumption that the total strain can be decomposed into the sum of elastic strain and inelastic slip strain; after the iterative solution in step S103 converges, the current time t has been obtained. n+1 The final sliding kinematic field ε s,n+1 Based on this, the true stress tensor σ at the current moment is calculated. n+1 :

[0140] σ n+1 =C el :(ε n+1 -ε s,n+1 );

[0141] In the formula: σ n+1 For time t n+1 The final true stress tensor; C el ε is the fourth-order elastic stiffness tensor of the material; n+1 Let t be the current time. n+1 The total strain tensor, ε n+1 For ε n +Δε;ε s,n+1 The solution obtained in step S103 at time t n+1 The value of the converged sliding kinematic field.

[0142] The true stress tensor σ n+1 It reflects the true stress state of the integration point under the current total strain after all internal state evolutions have been experienced.

[0143] S104-2: Calculate the tangent stiffness matrix corresponding to the currently updated state; to ensure the quadratic convergence of the global Newton-Raphson iteration in nonlinear finite element analysis, it is necessary to calculate a tangent stiffness matrix C that is completely consistent with the complex coupled constitutive model described in this invention. alg This tangent stiffness matrix is ​​also referred to as the algorithmic tangent modulus in the literature. This tangent stiffness matrix is ​​defined as the true stress tensor σ at the current moment. n+1 With respect to the total strain tensor ε n+1 Total differential: ;

[0144] The derivation and calculation of the tangent stiffness matrix, as a functional generalization technique, relies on the application of the chain rule and implicit function theorem to the entire implicit solution process. Since in the solution process of step S103, all internal state variables y={ε s ,μ s ,φ} T It has been solved as the total strain tensor ε n+1 The implicit function, which is derived from the convergent residual equation R(y(ε)). n+1 ),ε n+1 The result is determined by ) = 0; taking the total differential of this residual equation, we get: ;

[0145] Therefore, the sensitivity of the internal state variables to the total strain tensor can be obtained: ;

[0146] Substituting this relationship into the stress tensor σ n+1 From the total differential expression, we finally obtain the formula for calculating the tangent stiffness matrix: .

[0147] in, It is by The components corresponding to the slip kinematic field are extracted; this calculation process explicitly considers how small changes in the total strain tensor, through coupled evolutionary rules, cause corresponding changes in all internal state variables, and are ultimately reflected in the rate of change of the stress tensor; the precise tangent stiffness matrix C alg will be compared with the true stress tensor σ n+1 Together, they are returned to the main analysis program in step S105.

[0148] In step S105, the modeling method of the present invention is typically run as a user-defined material subroutine within a main finite element analysis program. In this step, the updated mechanical response information calculated in step S104 is returned to the main analysis program for it to complete a global equilibrium iteration. The specific implementation of this step is as follows:

[0149] S105-1: At each integration point, after completing the iterative solution of the internal state variables and the updated calculation of the mechanical response, the method of this invention combines two core results: the true stress tensor σ at the current moment. n+1 and the consistent tangent stiffness matrix C corresponding to this state alg The process returns or writes to the memory address specified by the main analysis program through a predefined interface; this process is executed on all integration points of all cable elements or beam elements in the finite element model that use the method of this invention within the current increment step.

[0150] S105-2: The main analysis program receives the true stress tensor σ returned from all integration points. n+1 Then, by numerically integrating the volume of each cable element or beam element, the internal nodal force vectors of each cable element or beam element are calculated. Subsequently, the main analysis program assembles the internal nodal force vectors of all cable elements or beam elements into the global internal force vector F of the entire structure in the current state, following the standard finite element assembly process. int .

[0151] S105-3: Simultaneously, the main analysis program receives the consistent tangent stiffness matrix C returned by all integration points. alg Subsequently, by numerically integrating the volume of each cable element or beam element and following a standard assembly process, the global tangential stiffness matrix K of the entire structure in its current state is constructed. T Using the consistent tangent stiffness matrix calculated in step S104-2 can ensure that the global Newton-Raphson iteration has the characteristic of quadratic convergence, thus improving computational efficiency.

[0152] S105-4: The main analysis program utilizes the assembled global internal force vector F int And the global external force vector F corresponding to the known current load increment step ext Calculate the global unbalanced force vector, i.e., the residual vector R, for the current iteration step: R = F ext -F int .

[0153] S105-5: The main analysis program utilizes the residual vector R and the global tangent stiffness matrix K. T Solve a large system of linear equations to obtain the correction ΔU:K for the global nodal displacements. T ΔU=R.

[0154] S105-6: The main analysis program applies the calculated displacement correction ΔU to the nodal displacements of the finite element model to complete one global equilibrium iteration. Subsequently, the main analysis program checks whether the norm of the global residual vector or the norm of the displacement correction satisfies the preset convergence criterion. If not, the main analysis program will call the method described in this invention again based on the updated displacement field (i.e., repeat steps S102 to S105) to enter the next global equilibrium iteration until convergence.

[0155] In step S106, the calculation process constituted in steps S102 to S105 is embedded into the overall time or load increment step loop of the main analysis program to simulate the response of the main cable throughout the analysis process; the specific implementation of this step is as follows:

[0156] S106-1: In a load or time increment step (e.g., from time t) nto t n+1 In the calculation of ), the incremental step is considered complete after the main analysis program has completed all global equilibrium iterations and confirmed global convergence; at this point, at each integration point of the main cable model, the current time t n+1 All internal state variables, including the final true stress tensor σ n+1 Total strain tensor ε n+1 And the updated state vector S (at the current time t) calculated by the method of this invention. n+1 The value is denoted as S n+1 ={Θ n+1 ,μ s,n+1 ,φ n+1 ,ε s,n+1 All of these are solidified and stored.

[0157] S106-2: The main analysis program's analysis controller then determines whether all preset load increment steps have been completed or the total analysis time has been reached; if not, the analysis will proceed to the next increment step, for example, from time t. n+1 to t n+2 .

[0158] S106-3: At the start of a new increment step, the main analysis program calculates the new global displacement increment based on the applied load or time step, thereby obtaining the total strain increment at each integration point; subsequently, the main analysis program calls the method described in this invention again, using the new total strain increment as the driving force input, and loads all state vectors S stored in step S106-1 at the end of the previous increment step. n+1 This serves as the initial condition for this round of calculations.

[0159] S106-4: Repeat the complete process from step S102 to step S105: Start by obtaining the new total strain increment, perform elastic prediction, calculate the coupled evolution of each internal state variable driven by the total strain increment by returning the mapping algorithm and iterative solution, update the true stress tensor and consistent tangent stiffness matrix at the integration point, and return them to the main analysis program for global equilibrium iteration of the current increment step.

[0160] S106-5: This incremental step cycle process is the repeated execution of steps S102 to S105 until all preset load incremental steps or time incremental steps have been calculated. In this way, the method provided can track and accumulate the nonlinear and history-dependent responses inside the main cable caused by wire slippage, wear degradation and topology rearrangement step by step.

[0161] Finally, when the entire analysis was completed, the results obtained not only included the macroscopic responses of the main cable such as displacement, internal force, and stress at various moments, but also fully recorded the detailed evolution history of each physical field (section topology field, friction interface state memory field, slip activation field, and slip kinematic field) throughout the analysis process, thus realizing a refined numerical simulation of the nonlinear behavior of the main cable.

[0162] Although embodiments of the invention have been shown and described, it will be understood by those skilled in the art that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the appended claims and their equivalents.

Claims

1. A nonlinear numerical simulation modeling method for the main cable of a suspension bridge, characterized in that, Includes the following steps: Step S1: Establish the initial finite element model of the main cable of the suspension bridge, and initialize a set of internal state variables at the integration point of each cable element in the finite element model, wherein the internal state variables include: The cross-sectional topology field characterizes the wire stacking pattern within the cross section, the friction interface state memory field characterizes the cumulative wear effect between the wires, the slip activation field characterizes the degree of transition from the locked state to the slip state of the cross section, and the slip kinematic field quantifies the irreversible deformation caused by the relative slip between the wires. Step S2: In any incremental step of the finite element analysis, obtain the total strain increment corresponding to the incremental step; Step S3: Based on the preset, mutually coupled evolution control equations, solve for the updated values ​​of each internal state variable driven by the total strain increment; Step S4: Based on the sliding kinematic field updated by one of the internal state variables, calculate the true stress tensor and the tangent stiffness matrix corresponding to the current state at the integration point at the end of the current increment step. Step S5: Transfer the calculated true stress tensor and tangent stiffness matrix to the master finite element analysis program; Step S6: After the calculation of one incremental step is completed, the updated internal state variables obtained in the incremental step are used as the initial conditions for the next incremental step, and steps S2 to S5 are repeated until all incremental steps are completed, thus completing the nonlinear numerical simulation modeling of the finite element model of the suspension bridge main cable. The evolutionary governing equations describe: the degenerative evolution of the friction interface state memory field driven by frictional dissipation energy; the activation evolution of the slip activation field modulated by the friction interface state memory field and stress state; the evolution of the cross-sectional topological configuration field driven by the cumulative slip norm; and the evolution of the slip kinematic field modulated nonlinearly by the slip activation field. The cumulative slip norm is obtained by integrating the evolution rate of the slip kinematic field over time.

2. The nonlinear numerical simulation modeling method for the main cable of a suspension bridge according to claim 1, characterized in that, After step S2 and before step S3, the process also includes the step of performing elasticity prediction calculations: Based on the assumption that the sliding kinematic field remains unchanged, a trial stress tensor is calculated according to the converged state at the end of the previous increment step and the total strain increment of the current increment step, and the trial stress tensor is used as the mechanical driving force for the evolution of the internal state variables.

3. The nonlinear numerical simulation modeling method for the main cable of a suspension bridge according to claim 1, characterized in that, Step S3 involves solving for the evolution of the cross-sectional topological configuration field, specifically including: The cumulative slip norm is compared with a preset critical threshold that depends on the current topology. If the critical threshold is exceeded, the transition probability from the current topology to each target topology is calculated based on the potential energy function associated with each topology. The cross-sectional topology field is determined and updated based on the aforementioned transition probability distribution.

4. The nonlinear numerical simulation modeling method for the main cable of a suspension bridge according to claim 1, characterized in that, In step S3, the coupling relationship between the friction interface state memory field and the slip activation field is solved, specifically including: Solve for the degradation evolution of the state memory field of the friction interface driven by frictional dissipation energy; Solve for the activation evolution of the slip activation field, wherein the slip driving energy is determined based on the friction interface state memory field after degradation evolution, and the slip activation field is updated based on the determined slip driving energy. The value of the slip driving energy depends on the equivalent stress determined by the stress state and the cross-sectional topology field and the critical stress threshold determined by the friction interface state memory field.

5. The nonlinear numerical simulation modeling method for the main cable of a suspension bridge according to claim 4, characterized in that, The activation evolution of the slip activation field is governed by the phase field equation, where: The driving force for the activation evolution originates from the reduction of the pseudo-free energy functional; The density function of the pseudo-free energy functional consists of a chemical energy term and a gradient energy term. The chemical energy term includes a slip-driven energy related term, which is specifically the product of a monotonically increasing function related to the slip activation field and the slip-driven energy.

6. The nonlinear numerical simulation modeling method for the main cable of a suspension bridge according to claim 1, characterized in that, Step S3 involves solving the nonlinear modulation coupling of the slip activation field to the slip kinematic field, specifically including: An evolutionary law for the sliding kinematic field is established, which stipulates that the evolution rate of the sliding kinematic field depends on the sliding rate multiplier. The value of the slip rate multiplier is calculated from the value of the slip activation field using a nonlinear modulation function.

7. The nonlinear numerical simulation modeling method for the main cable of a suspension bridge according to claim 6, characterized in that, The evolution law of the slip kinematic field is a flow law, wherein the evolution rate of the slip kinematic field is equal to the product of the slip rate multiplier and the gradient of the slip yield function with respect to the true stress tensor, and the slip yield function depends on the true stress tensor and the cross-sectional topological configuration field.

8. The nonlinear numerical simulation modeling method for the main cable of a suspension bridge according to claim 1, characterized in that, In step S3, solving the mutually coupled evolution control equations involves executing a return mapping algorithm. The return mapping algorithm includes discretizing the first-order ordinary differential equations describing the evolution of each internal state variable through an implicit time integration scheme to form a nonlinear algebraic equation system.

9. The nonlinear numerical simulation modeling method for the main cable of a suspension bridge according to claim 8, characterized in that, Solving systems of nonlinear algebraic equations using the Newton-Raphson iterative method, which includes the following steps: In each iteration, a residual vector consisting of the constraints of each internal state variable is constructed, and the Jacobian matrix is ​​calculated based on the constructed residual vector. The Jacobian matrix is ​​obtained by taking the partial derivatives of the residual vector with respect to all the internal state variables to be solved. Solve the system of linear equations consisting of the Jacobian matrix and the residual vector to obtain the incremental correction vector of the internal state variables; The internal state variables are updated using the incremental correction vector until convergence.

10. The nonlinear numerical simulation modeling method for the main cable of a suspension bridge according to claim 9, characterized in that, In step S4, the tangent stiffness matrix corresponding to the current updated state is calculated. The tangent stiffness matrix is ​​a consistent tangent stiffness matrix. Furthermore, the calculation of the consistent tangent stiffness matrix uses the Jacobian matrix obtained in the Newton-Raphson iteration method.

Images