A numerical simulation method for avalanche impact protection structure based on bidirectional coupling

Abstract

The application discloses a kind of avalanche impact protection structure numerical simulation methods based on two-way coupling, it is related to mountain disaster protection technical field, including: avalanche source modeling and initialization: discrete element theory is used to simulate avalanche process, realize the whole process dynamics simulation and reproduction of avalanche disaster from start to movement to impact protection structure;Protection structure modeling and parameter setting: based on general finite element simulation platform, modeling is carried out to avalanche protection structure;Discrete element-finite element contact coupling setting: unit stress, strain data are solved;Numerical simulation control and solution;Impact response extraction and result analysis: the systematic analysis of structure response and impact characteristics is carried out, to evaluate the dynamic performance and safety margin of structure under the action of avalanche impact;The application can be used to predict the impact effect of avalanche disaster on downstream structure, and provide theoretical basis and technical support for the design and performance evaluation of protection structure.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of mountain disaster prevention, and particularly relates to a numerical simulation method of avalanche impact protection structure based on bidirectional coupling. BACKGROUND

[0002] Avalanche disaster has the characteristics of strong burst, great impact force and long propagation distance, which seriously threatens the safety of mountain traffic, engineering facilities and personnel. In recent years, with the continuous advancement of infrastructure construction in plateau, high-cold and high-mountain areas, the influence of avalanche disaster on engineering structures has been increasingly concerned. Accurate simulation of the interaction process between avalanche and structure is of great significance for the scientific design of protection structure and disaster risk assessment.

[0003] At present, the numerical simulation methods of avalanche disaster mainly include smoothed particle hydrodynamics (SPH) method, material point method (MPM), finite volume method (FVM), discrete element method (DEM) and finite element method (FEM) and the like. Among them, the discrete element method can effectively depict the contact behavior and flow process between avalanche particles, and is suitable for simulating the multi-particle dynamic characteristics of avalanche; and the finite element method has wide application in structure response analysis and mechanical behavior simulation. In order to realize high-precision coupling simulation of avalanche and structure, researchers have begun to try to build DEM-FEM coupling model, and model the avalanche particle group and structure response respectively, so as to more truly reproduce the dynamic response and energy transfer mechanism in the impact process.

[0004] Although the existing researches have revealed the particle dynamics and structure response characteristics in the avalanche impact process to some extent, there is still a lack of a general, efficient and applicable DEM-FEM coupling simulation method for complex terrain and various structure types, especially in the simulation of the impact effect of avalanche on flexible protection net, bridge pier and other engineering structures. The existing numerical simulation methods of avalanche have obvious shortcomings in depicting the complex dynamic interaction between avalanche and downstream protection structure, which are embodied in the following aspects:

[0005] (1) The quasi-multiphase model based on depth average mainly adopts one-dimensional or two-dimensional approximation simplification, which is difficult to capture the non-continuous motion characteristics and local structure effect between avalanche particles, and cannot accurately reflect the key dynamic behaviors such as particle size distribution and local impact concentration, thereby limiting its applicability in protection structure response assessment.

[0006] (2) The existing simulation research based on discrete element-finite element coupling mainly focuses on the impact process of debris flow or landslide-accumulation body, and the applicable objects are mainly viscous substances or large-size blocks, whose parameter setting and physical properties are difficult to represent the characteristics of high-speed, low-viscosity and small particle size of avalanche material, so it is not suitable for accurate modeling and protection analysis of avalanche disaster.

[0007] (3) Some studies simulate the transformation process of rock-ice avalanche to debris flow, but do not consider the whole process of interaction of avalanche impact structure, especially lack of systematic description of multi-scale dynamic process such as avalanche collision, structural deformation response and energy transmission path, which leads to limited guiding role in the design of protective structure.

[0008] In summary, the existing technical solutions have made some progress in simulating the overall motion behavior of avalanche, the transformation process of debris flow or the interaction between debris flow and structure, but still have significant limitations: first, most studies can only describe the overall flow characteristics of avalanche or debris flow, and cannot reveal the discrete dynamic behavior of avalanche material from the particle scale; second, some continuous medium models cannot accurately capture the nonlinear contact and impact process between avalanche particles and protective structures; third, the existing numerical methods generally ignore the dynamic response and cumulative damage characteristics of the structure, making it difficult to truly reflect the whole process of avalanche-structure interaction. Therefore, it is urgent to develop a numerical simulation method that can realize the two-way coupling between particle-scale avalanche dynamics and structural response, in order to improve the scientificity and safety of protective engineering design. The present invention is proposed in this background to solve the shortcomings of existing methods in precision, applicability and engineering application. SUMMARY

[0009] To solve the problems in the prior art, the purpose of the present invention is to provide a numerical simulation method for avalanche impact protective structure based on two-way coupling, which can be used to predict the impact effect of avalanche disaster on downstream structures and provide theoretical basis and technical support for the design and performance evaluation of protective structures.

[0010] To achieve the above-mentioned purpose, the technical solution adopted by the present invention is as follows: a numerical simulation method for avalanche impact protective structure based on two-way coupling, comprising the following steps:

[0011] Step 1, avalanche source modeling and initialization: using discrete element theory to simulate the avalanche process, first obtaining the basic physical parameters of avalanche particles; then initializing the particle group in the geometric boundary through random particle generation and filling algorithm; setting the initial velocity field, which can also be driven by the gravity of the slope surface to form a velocity field, realizing the whole process of dynamic simulation and reproduction of avalanche disaster from start to movement to impact on protective structure;

[0012] Step 2, modeling and parameter setting of protective structure: modeling the avalanche protective structure based on nonlinear finite element method, selecting appropriate solid elements, shell elements, beam elements and string elements to divide the grid for different types of structural components, establishing a finite element model, and setting material parameters, contact parameters and boundary constraint conditions for different components according to the structural characteristics;

[0013] Step 3, avalanche impact protection structure bidirectional coupling setting: the bidirectional coupling analysis between the avalanche and the protection structure is realized through the impact contact algorithm between the discrete element and the finite element; the real-time update of the snow particle position is realized through the capture of the particle displacement, velocity and acceleration; meanwhile, the search algorithm is combined to establish the contact relationship between the snow particles and the protection structure by updating the finite element node position and the motion information according to the dynamic equation; on this basis, the contact force between the avalanche particles and the protection structure is updated, the contact state is updated and the bidirectional coupling between the avalanche and the protection structure is realized through the physical characteristics and the contact parameters of the discrete element and the finite element contact position.

[0014] Step 4, numerical simulation control and solution: the integral time step is set, the step range meets the stability condition, the simulation total time is set according to the actual avalanche impact duration; in the simulation process, the energy output function is started, the change of the key energy items such as the impact kinetic energy, the structure internal energy, the sliding energy and the damping energy is recorded and monitored, which is used to evaluate the energy conversion process between the avalanche particles and the structure and the impact damage risk, the explicit time integral method is used for solution, the response variables such as the stress, the strain and the displacement of the structure at each time step are continuously output and are used for subsequent damage evaluation, structure performance analysis and optimization design.

[0015] Step 5, impact response extraction and result analysis: after the calculation simulation is completed, the key response indexes of the avalanche impact protection structure system are extracted, including the maximum impact force, the deformation response, the energy absorption rate, the node acceleration and the stress concentration area distribution of the structure; the velocity field, the impact frequency distribution and the contact force time history of the avalanche particles are analyzed, the main impact area of the avalanche is identified; the local and overall response time history of the protection structure is established in combination with the finite element result, and the stress characteristics and the damage mode of different structure components are analyzed.

[0016] As a further improvement of the application, in step 1, the basic physical parameters of the avalanche at least include the avalanche scale, the source distribution, the accumulation density and the friction coefficient.

[0017] As a further improvement of the application, in step 1, according to the inter-particle collision force, the gravity effect, the particle translation and rotation equations in the motion process of the snow particles, the equations are as follows:

[0018] ;

[0019] Wherein, is the mass of the i-th particle, is the particle acceleration, is the external force, is the moment of inertia, is the angular acceleration, is the torque applied to the i-th particle;

[0020] The explicit algorithm is used to solve the motion state of the particle:

[0021] (1) Time step Acceleration, angular acceleration update, the acceleration and angular acceleration of the i-th particle at time t are calculated as:

[0022]

[0023] where, , , respectively represent the particle displacement at time t, , , , , respectively represent the particle rotation angle at time t, , ,

[0024] (2) Half-step Velocity update, the velocity and angular velocity of the i-th particle at time t are calculated as:

[0025]

[0026] where, , respectively represent the velocity and angular velocity of the particle at time t; (3) Half-step

[0027] Displacement, rotation angle update, the displacement and rotation angle of the particle at time t are calculated as:

[0028]

[0029] where, , respectively represent the velocity and angular velocity of the particle at time t; (4) Explicit dynamics algorithm start condition: when

[0030] , the initial half-step velocity of the particle is represented by the following formula:

[0031] ​​​​​​​​​​​​​​​

[0032] wherein, , respectively as the speed of the particle, angular velocity start condition, , , , for moment known initial speed of the particle, initial acceleration, initial angular velocity, initial angular acceleration.

[0033] As a further improvement of the application, in step 2, for different types of structure, appropriate finite element unit is selected for modeling, including: for rigid structure, three-dimensional solid element or shell element is used to establish the finite element grid of structural members; for flexible structure, nonlinear beam element or equivalent shell element is used for modeling, and the mechanical properties of connecting members are defined.

[0034] As a further improvement of the application, in step 2, after modeling is completed, according to the mechanical properties of different structural materials, corresponding model parameters are set, including:

[0035] (1) for steel protective structure, material density, elastic modulus, Poisson's ratio, yield strength, ultimate strength, ultimate strain and damping coefficient are determined to ensure that the elastic-plastic behavior and energy dissipation characteristics of the material can be accurately reflected; (2) for concrete protective structure, material density, elastic modulus, Poisson's ratio, compressive strength, maximum aggregate size, and damping coefficient are set in combination with the damping characteristics of the structure to fully reflect the nonlinearity and damage evolution process of concrete; (3) for flexible structure and its connecting members, the equivalent stiffness, ultimate bearing capacity and slip performance parameters of nonlinear connection are set to accurately simulate the mechanical response and deformation capacity of the connecting members; in addition, the change characteristics of material parameters under the influence of environmental factors such as temperature and strain rate should also be considered to ensure the applicability and accuracy of the model.

[0036] As a further improvement of the application, in step 2, the boundary conditions of the structure model are set according to the actual working conditions, including: the nodes at the bottom of the structure are set as fixed constraints, and contact or sliding boundary is set between the structure and the foundation.

[0037] As a further improvement of the application, in step 2, the dynamics control equation of the avalanche protection structure simulation is as follows:

[0038] ;

[0039] wherein, , , respectively as the mass matrix, damping matrix, stiffness matrix of the structure, respectively as the displacement vector, velocity vector and acceleration vector of the node, is the external force vector;

[0040] To improve the computational efficiency and stability, the explicit second-order central difference scheme is adopted to discretize the time integration; the acceleration and velocity of the node are expressed as:

[0041] ;

[0042] According to the above formula, the dynamic control equation of the avalanche protection structure simulation is rewritten as the following formula:

[0043] ;

[0044] wherein, , are the node displacement vectors at the time of , ; by the form, the structure displacement is directly solved at each time step, without matrix inversion, high computational efficiency and suitable for strong nonlinear impact situation;

[0045] In order to start the explicit time integration algorithm, the starting displacement at the time of is initialized; based on Taylor expansion, the starting condition of the explicit dynamic algorithm and the initial half-step velocity of the structure unit are calculated by the following formula:

[0046] ;

[0047] wherein, , are the known initial position, initial velocity and initial acceleration of the finite element node at the time of .

[0048] As a further improvement of the present application, in step 3, the contact coupling of discrete element and finite element is realized by the contact algorithm based on penalty function; based on global search algorithm, when the discrete element particle and the finite element surface penetrate, the normal contact force is expressed as:

[0049] ;

[0050] wherein, is the penetration depth between the discrete element and the finite element, is the contact stiffness, is the contact damping coefficient, is the relative velocity of the contact node;

[0051] To ensure numerical stability and physical accuracy, by comparing the element-based stiffness between the discrete element and the finite element and the stability-based stiffness to determine the equivalent stiffness k of the contact interface as follows:

[0052] ;

[0053] where the cell-based stiffness k is determined based on the material parameters and cell characteristics as follows: e

[0054] ;

[0055] where, is a penalty function scaling factor, is the bulk modulus, is the contact area, is the solid cell volume, is the shell cell characteristic length;

[0056] Stiffness based on stability is related to the effective mass of the node and the global time step:

[0057] ;

[0058] where, is the effective mass of the node; the effective mass of the node is the smaller value of the two interacting node masses: ;

[0059] Effective angular frequency is calculated as:

[0060] ;

[0061] Considering the Coulomb friction between the discrete element and the finite element, the tangential contact force is related to the normal contact force as follows:

[0062] ;

[0063] where the friction coefficient is exponentially interpolated according to the relative velocity between the discrete element particle and the finite element node to simulate the smooth transition between static friction and dynamic friction:

[0064] ;

[0065] where, is a constant that controls the transition rate of the friction coefficient; when , the friction coefficient is equal to the static friction coefficient.

[0066] As a further improvement of the present application, step 5 is specifically as follows:​

[0067] (1) Firstly, the dynamic response time-history data of the key components of the protection structure are extracted from the finite element calculation results, including the impact position, anchoring position, impact deformation of the connected position structure component, peak internal force, stress and strain level and damage state, the local yielding, damage or anchoring failure mode of the structure in the impact process is analyzed, and the accurate evaluation of the safety performance of the structure is realized;

[0068] (2) Secondly, the maximum impact load of the avalanche on the protection structure is determined, the contact force time-history curve of the whole impact process is extracted, the peak value of the contact force is determined as the maximum impact force index, which is used for the strength checking and impact resistance performance evaluation of the structure, and the rationality and accuracy of the two-way coupled numerical model are verified through comparison with theoretical analysis and test data;

[0069] (3) Finally, the sensitivity analysis is carried out based on different working conditions, the key control parameters affecting the impact response of the structure are identified by comparing the differences of the structure responses under various working conditions; the local and overall response time-history model of the protection structure is constructed based on the avalanche impact force, deformation of the protection structure and stress concentration region distribution, and the stress characteristics and potential damage mechanism of each structure component are analyzed in depth, so as to provide a scientific basis for the structure optimization design.

[0070] The avalanche-structure numerical simulation method based on the discrete element-finite element coupling theory of the present application solves the problems of imperfect description of dynamic behavior, difficulty in reflecting particle effect and insufficient structure response precision in the existing technology in the modeling of the interaction between avalanche disaster and downstream structure. The method can accurately simulate the dynamic response characteristics of the protection structure under the impact of the avalanche by constructing a three-dimensional particle scale avalanche dynamics model, combining with the structure finite element model, and coupling to realize the numerical simulation of the whole process of avalanche particle impact, energy transmission and structure response, and can provide reliable basis for structure performance evaluation and disaster prevention design.

[0071] The beneficial effects of the present application are:

[0072] The present application proposes an avalanche-structure interaction numerical simulation method based on discrete element-finite element coupling to solve the problems of insufficient particle scale dynamics simulation, insufficient coupling between avalanche and prevention structure, and difficulty in accurately capturing the spatio-temporal distribution of impact force and structure deformation mechanism in the existing avalanche prevention technology, which significantly improves the scientificity and practicality of the design of avalanche disaster prevention structure, and achieves the following technical effects:

[0073] (1) Improve simulation accuracy and truly reflect avalanche dynamics: The present application uses the discrete element method to accurately simulate the motion of particles in an avalanche, combined with the finite element method to simulate the nonlinear dynamic response of the prevention structure, realizing the multi-physical field coupling of particle scale and structure response. The simulation results show that the velocity distribution, collision frequency and force time history of the avalanche particles can be accurately captured, and the transient and spatial variation characteristics of the avalanche impact load are truly reflected.

[0074] (2) Reveal the avalanche-structure coupling mechanism and structure energy dissipation law: Through coupling numerical simulation and model test verification, the present application reveals the development process of structure deformation, sliding and damage in the avalanche impact process, and clearly defines the stress transfer path and energy dissipation mechanism of the key components of the structure. The comparison between the test and the numerical simulation shows that the energy dissipation efficiency of the structure can be improved by 15%-25%, significantly enhancing the impact resistance of the prevention structure and prolonging the service life of the structure.

[0075] (3) Promote the progress of avalanche disaster prevention technology in high-cold mountainous areas and enhance social security: The present application provides scientific and technological support for the prevention of avalanche disasters in high-cold and high-altitude areas with complex terrain, supporting precise assessment of avalanche risk and customized design of protection structures. The engineering case designed by the present application method realizes a 30% reduction in the peak value of avalanche impact force in the avalanche danger area, effectively protecting the safety of mountain transportation, residents and infrastructure, and improving the ability of social stability and sustainable development. BRIEF DESCRIPTION OF DRAWINGS

[0076] Figure 1 is a flowchart of an embodiment of the present application;

[0077] Figure 2 is a schematic diagram of an avalanche impact protection structure in an embodiment of the present application.

[0078] REFERENCE NUMERALS

[0079] 1, three-dimensional terrain of avalanche risk area, 2, avalanche source, 3, protection structure, 4, protection object. DETAILED DESCRIPTION

[0080] The embodiments of the present application will be described in detail below with reference to the accompanying drawings.

[0081] EMBODIMENT

[0082] As shown in Figure 1 , a numerical simulation method for avalanche impact protection structure based on bidirectional coupling includes the following steps:

[0083] Step 1: Avalanche source modeling and initialization:

[0084] Discrete element method (DEM) is a numerical method for simulating the dynamics of non-continuous media with large deformation, which can capture the collision, rolling and accumulation behaviors of avalanche particles. To simulate the avalanche process using DEM, the basic physical parameters of the avalanche need to be obtained first, including the avalanche scale, source distribution, and packing density (e.g., snow volume 5000 m³, snow particle size range 0.01 m to 0.1 m, dry snow density 100-300 kg / m³, wet snow density 400-600 kg / m³, friction coefficient 0.1-0.5). In addition, the particle group is initialized within the geometric boundary through a random particle generation and filling algorithm. It is recommended that the number of particles be less than 10 9 to ensure the stability of the avalanche dynamics and computational efficiency. The initial velocity field can be set (usually estimated using fluid dynamics parameters such as mass flow rate , average velocity , volume fraction , etc.), or a velocity field can be formed by gravity driving on the slope surface to realize the whole process of avalanche disaster from start to movement to impact on protective structures. Considering the inter-particle collision force and gravity acting on the snow particles during movement, the translational and rotational equations of the particles are as follows:

[0085]

[0086] where is the mass of the i-th particle, is the particle acceleration, is the external force, is the moment of inertia, is the angular acceleration, is the torque applied to the i-th particle.

[0087] The explicit algorithm is used to solve the motion state of the particles:

[0088] (1) Time step acceleration and angular acceleration update, calculate the i-th particle acceleration and angular acceleration at time t:

[0089]

[0090] where , , represent the particle displacement at times , , , respectively; , , represent the particle velocity at times ,​ , the particle's angular displacement at time t;

[0091] (2) Half-step velocity update, the particle's velocity at time t and angular velocity are respectively expressed as:

[0092]

[0093] wherein, , respectively represent the particle's velocity and angular velocity at time t;

[0094] (3) Half-step displacement and angular displacement update, the particle's displacement and angular displacement at time t are calculated as:

[0095]

[0096] wherein, , respectively represent the particle's velocity and angular velocity at time t;

[0097] (4) Explicit dynamics algorithm start condition: when the particle's initial half-step velocity at time t is expressed as:

[0098]

[0099] wherein, , respectively as the particle's velocity and angular velocity start condition, , , , the particle's initial velocity, initial acceleration, initial angular velocity, and initial angular acceleration at time t are known.

[0100] Step 2: Protective structure modeling and parameter setting:

[0101] ​​​This step is based on the nonlinear finite element method (such as ANSYS / LS-DYNA, Abaqus, etc.) to model the protective structure. For different types of structural components, select appropriate solid elements, shell elements, beam elements, and mesh elements to divide the grid, establish the finite element model, and set different component material parameters, contact parameters, and boundary constraint conditions according to the structure characteristics: for rigid structures (such as reinforced concrete retaining wall, arched shed tunnel, etc.), three-dimensional solid elements or shell elements are usually used to establish the finite element grid of structural components; for flexible structures (such as flexible protective net, steel cable anchor rod system, etc.), nonlinear beam elements or equivalent shell elements can be used for modeling, and the mechanical properties of connecting components such as steel cable, anchor rod, and ring net are defined to truly reflect their stress and deformation behavior. After modeling is completed, set the corresponding model parameters according to the mechanical properties of different structural materials. For steel protective structures, determine the material density, elastic modulus, Poisson's ratio, yield strength, ultimate strength, ultimate strain, and damping coefficient to ensure accurate reflection of the material's elastic-plastic behavior and energy dissipation characteristics. The elastic modulus is generally set to about 2.1×10¹¹ Pa, the yield strength of Q235 steel is about 235 MPa, the Poisson's ratio is 0.3, and the damping coefficient is valued according to experience. For concrete protective structures, the material density, elastic modulus, Poisson's ratio, compressive strength, maximum aggregate size, and damping coefficient are set in combination with the structural damping characteristics to fully reflect the nonlinear and damage evolution process of concrete. The elastic modulus is usually between 2.0×10¹ 0 Pa and 3.5×10¹ 0 Pa, the compressive strength can be set to 30 MPa corresponding to the C30 grade, the Poisson's ratio is about 0.2-0.25, and the damping can be taken as 0.03-0.05. For flexible structures and their connecting components, set the equivalent stiffness, ultimate bearing capacity, and slip performance parameters of the nonlinear connection to accurately simulate the mechanical response and deformation capacity of the connecting components. In addition, the variation characteristics of material parameters under the influence of environmental factors such as temperature and strain rate should also be considered to ensure the applicability and accuracy of the model. Finally, set the boundary conditions of the structure model according to the actual working conditions, such as setting the nodes at the bottom of the structure as fixed constraints, setting the contact or sliding boundary between the structure and the foundation, etc. to ensure that the model's constraint conditions are consistent with the actual conditions during the process of being impacted by the avalanche load. The dynamics control equation of the avalanche protective structure simulation can be expressed as:

[0102]

[0103] where, , 、 are the mass matrix, damping matrix, and stiffness matrix of the structure, are the node displacement vector, velocity vector, and acceleration vector, is the external force vector.

[0104] To improve the computational efficiency and stability, the explicit second-order central difference scheme is adopted to discretize the time integration. The acceleration and velocity of a node can be expressed as:

[0105]

[0106] According to the above formula, the dynamic governing equation of avalanche protection structure simulation can be rewritten as follows:

[0107]

[0108] where, , are the displacement vectors of the nodes at time , . Through the form, the structure displacement can be directly solved explicitly at each time step, without matrix inversion, high computational efficiency and suitable for strong nonlinear impact situation. In order to start the explicit time integration algorithm, the starting displacement at time needs to be initialized. Based on Taylor expansion, the starting condition of explicit dynamic algorithm and the initial half-step velocity of structure element can be calculated by the following formula:

[0109]

[0110] where, , are the known initial position, initial velocity, initial acceleration of the finite element node at time .

[0111] Step 3: Avalanche impact protection structure bi-directional coupling setting:

[0112] The bi-directional coupling analysis between avalanche and protection structure is realized through the impact contact algorithm between discrete elements and finite elements. By capturing the displacement, velocity and acceleration of the particles, the position of the snow particles is updated in real time. At the same time, the node position and motion information of the finite element are updated through the dynamic equation, and the search algorithm is used to establish the contact relationship between the snow particles and the protection structure. On this basis, through the physical properties and contact parameters of the discrete element and finite element contact position, the contact force between the avalanche particles and the protection structure is updated, the contact state is updated, and the bi-directional coupling between the avalanche and the protection structure is realized.

[0113] The coupling of discrete element (DEM) and finite element (FEM) can be realized through the contact algorithm based on penalty function. Based on the global search algorithm, when the DEM particle and the FEM surface penetrate, the normal contact force is expressed as:

[0114]

[0115] where, is the penetration depth between discrete element and finite element, is the contact stiffness, is the contact damping coefficient, is the relative velocity of contact nodes;

[0116] To ensure numerical stability and physical accuracy, the equivalent stiffness of contact interface k is determined by comparing the element-based stiffness and the stability-based stiffness as follows:

[0117]

[0118] where, the element-based stiffness k e is determined by material parameters and element characteristics:

[0119]

[0120] where, is the penalty function scaling factor, is the bulk modulus, is the contact area, is the solid element volume, is the shell element characteristic length;

[0121] The stability-based stiffness k is related to the effective mass of nodes and the global time step:

[0122]

[0123] where, is the effective node mass; the effective node mass is taken as the smaller value of two interacting node masses: ;

[0124] The effective angular frequency of the contact interface is calculated as:

[0125] ;

[0126] Considering the Coulomb friction between discrete element and finite element, the relationship between the tangential contact force and the normal contact force can be expressed as:

[0127]

[0128] where, the friction coefficient The relative velocity between the discrete element particle and the finite element node is exponentially interpolated to simulate the smooth transition between static and dynamic friction:

[0129] ;

[0130] wherein, is a constant that controls the transition rate of the friction coefficient; when , the friction coefficient is equal to the static friction coefficient.

[0131] On the basis of the above theoretical framework, the discrete element particles and the finite element structure are coupled through a coupling interface (such as the DE_TO_SURFACE_COUPLING coupling interface based on ANSYS / LS-DYNA). The contact coupling relationship between the discrete element particles and the finite element structure is established. The contact model parameters, including the contact stiffness, the friction coefficient, and the bonding strength (such as considering the freezing or snow bonding of particles), are set to truly reflect the transmission of the avalanche impact force.

[0132] Step 4: Numerical simulation control and solution:

[0133] Under the action of the avalanche impact, the numerical solution process needs to satisfy the basic physical laws of momentum conservation, mass conservation, and energy conservation. In the process of discrete element and finite element coupling simulation, the position and velocity information of the particles and the DEM time step are continuously tracked to determine the distance relationship between the particles and the finite element structure nodes in real time, thereby determining whether contact occurs. When the contact condition is met, the contact force is calculated according to the contact algorithm, and it is applied to the structure element node, and then the position and velocity of the structure node are updated through the finite element dynamics balance equation to realize the dynamic response simulation of the structure deformation. To ensure the numerical stability of the explicit integration process, the integration time step should be reasonably set, and the step range should satisfy the stability condition. The total simulation time needs to be set according to the actual duration of the avalanche impact, which is usually 1-10 seconds. During the simulation process, the system energy output function needs to be turned on to record and monitor the changes of key energy items such as impact kinetic energy, structure internal energy, sliding energy, and damping energy, so as to evaluate the energy conversion process between the avalanche particles and the structure and the impact damage risk. The entire coupling system is solved by using the explicit time integration method, which can continuously output the response variables such as stress, strain, and displacement of the structure at each time step, and can be used for subsequent damage assessment, structure performance analysis, and optimization design.

[0134] Step 5: Impact response extraction and result analysis:

[0135] After completing the computational simulation, key response indicators of the avalanche impact protection structure system were extracted, including the maximum impact force, deformation response, energy absorption rate, nodal acceleration, and stress concentration area distribution. The velocity field, impact frequency distribution, and contact force time history of avalanche particles were analyzed to identify the main impact zones. Based on the finite element analysis results, local and overall response time histories of the protection structure were established to analyze the stress characteristics and failure modes of different structural components. First, dynamic response time history data of key components of the protection structure were extracted from the finite element calculation results, including impact deformation, peak internal force, stress-strain levels, and damage states of structural components at impact locations, anchorage locations, and connection locations. Failure modes such as local yielding, failure, or anchorage failure during the impact process were analyzed to achieve an accurate assessment of the structural safety performance. Second, the maximum impact load of the avalanche on the protection structure was determined, and the contact force time history curve of the entire impact process was extracted. The peak contact force was determined as the maximum impact force indicator for structural strength verification and impact resistance performance evaluation. Simultaneously, the rationality and accuracy of the two-way coupled numerical model were verified by comparing it with theoretical analysis and experimental data. Finally, sensitivity analysis was conducted under different working conditions. By comparing the differences in structural response under each condition, key control parameters affecting the structural impact response were identified. Combining avalanche impact force, protective structure deformation, and stress concentration area distribution, local and overall response time history models of the protective structure were constructed. This allowed for in-depth analysis of the stress characteristics and potential failure mechanisms of each structural component, providing a scientific basis for structural optimization design.

[0136] To verify the feasibility and effectiveness of the "numerical simulation method for avalanches based on discrete element-finite element coupling" proposed in this embodiment, such as... Figure 2 As shown, Figure 2 In this model, 1 represents the 3D topography of the avalanche risk area, 2 represents the avalanche source (Discrete Element Method / DEM), 3 represents the protective structure (Finite Element Method / FEM), and 4 represents the protected object (transportation route). A typical high-altitude avalanche hazard area is selected as the research object, and numerical simulation experiments are conducted. A full-process dynamic response analysis is then performed, combining real terrain, structural parameters, and avalanche conditions. The entire implementation process, following the described method, includes steps such as model preparation, geometric construction, parameter assignment, simulation control and solution, and dynamic response extraction.

[0137] Following step 1: Based on high-resolution remote sensing imagery and field survey data, the typical slope topography of the study area was reconstructed, and a digital elevation model (DEM) was used to construct the geometric boundaries of the simulated terrain. The simulation domain included the avalanche source area, the sliding path area, and the protective structure area. The source area was located at the top of the slope, with gravity driving the avalanche's descent. The initial avalanche velocity range was set to 10–25 m / s, corresponding to a typical wet avalanche. Avalanche particles were simulated using the discrete element method, with a particle radius of 10 ± 8 mm and a density of 500 kg / m³. A spherical particle group based on compressed particle contact was used, considering interparticle friction (coefficient of friction). = 0.5) and cohesion (cohesion strength = 2 kPa).

[0138] According to Step 2: The protective structure is modeled using the finite element method, and a flexible protective net is selected as the typical structural form. The equivalent membrane element is used to construct the flexible net, and the wire element is used to simulate the real steel wire rope stretching response. The net is arranged transversely along the slope bottom and is tensioned between the two side anchor points. The net height is 3 m, the total width is 10 m, and the anchor points are set as fixed constraints. The protective net and the avalanche particles are set as a contact coupling surface, allowing the particles to penetrate and generate resistance and net surface deformation. The net surface material parameters: the equivalent Young's modulus is = 50 MPa, the Poisson's ratio = 0.3, considering the geometric nonlinear large deformation response. The coupling interface uses the master-slave node contact algorithm, allowing the snow particles to interact with the structure during the impact process and transmit the impact load in real time.

[0139] According to Step 3: The snow particles use a soft ball contact model, and a spring-damper element is introduced to characterize the contact response between particles and between particles and the structure. The contact stiffness is calculated according to the particle material properties, and the contact damping is used to dissipate kinetic energy. The friction force calculation follows the Coulomb friction theory, and the cohesion force is realized using the critical shear criterion. The structure part uses a linear elastic material model, which can be extended to consider damage or plastic evolution in the later stage. The overall model uses a coupled explicit dynamics framework, and uses a unified mass-force format to calculate the mutual response of the structure and particles at each time step.

[0140] According to Step 4: After the model is established, the total simulation time is set to 5-10 seconds, and the integration time step is controlled between 1×10⁻ 6 and 5×10⁻ 5 seconds to meet the stability conditions (CFL conditions) in the contact dynamic process. The particle velocity, displacement, and contact state are updated in real time during the simulation, and the coupling equation is automatically called when the particle and the structure are in contact, completing the mechanical transmission of the particle-structure impact response. The explicit central difference integration method is used for numerical solution, which has the characteristics of high efficiency and stability, and is suitable for handling complex dynamic problems such as large deformation, multiple contacts, and multiple particle systems. The simulation process outputs the avalanche motion trajectory, particle distribution state, and structure displacement cloud diagram, which is convenient for post-processing analysis.

[0141] According to step 5: after the simulation is finished, the stress time history of the nodes at the top and middle anchoring positions of the structure, the maximum displacement, the acceleration response, and the overall deformation figure at different time steps are extracted from the model, and the impact force time history and the maximum impact force of the particles on the structure are further counted and analyzed for the spatiotemporal distribution law. Finally, the sensitivity analysis is performed by changing the particle diameter, the initial velocity, and the structure stiffness and other parameters, and it is found that the particle size is the most critical factor affecting the peak impact force, and the reduction of the structure stiffness can effectively improve the deformation capacity of the protective structure and reduce the risk of local damage.

[0142] The avalanche flow depth is generally not more than 10 m, the average speed is 5-25 m / s, and the density is 150-500 kg / m 3 According to the cohesion, water content and density, speed and terrain of snow, the avalanche can be divided into two flow modes of inertia type and gravity type. The inertia type avalanche has a higher speed (>10 m / s), can pass over the protective engineering at a lower position of the terrain, and can absorb the surrounding air during the descending process to form a density stratification. The gravity type avalanche has a lower speed (≤10 m / s), closely moves along the terrain, and presents a granular or viscous flow state. The snow temperature is a key parameter for controlling the cohesion and overall friction of snow. When the snow temperature is lower than-1°C, the inertia effect of the avalanche is dominant, and when the snow temperature is higher than-1°C, the gravity effect is dominant.

[0143] Since it is very difficult to directly measure the properties of the snow inside the avalanche, the parameters of the snow can be valued according to the reported mechanical properties of the snow as shown in Table 1: the particle radius = 10±8 mm, the density = 500 kg / m³, the average bulk density ranges from 338-379 kg / m³. The compressible particle model is adopted, the elastic modulus = 10 5 Pa, the cohesive strength , and the friction coefficient are taken as = 0.5.

[0144] Table 1 Reference values of the parameters of the heavy density avalanche material

[0145]

[0146] wherein the bulk density of the snow refers to the total mass of the snow particles (density ) in a unit volume in the stacking state.

[0147] The coupling simulation method described in the present application can truly reflect the avalanche-structure interaction process, has good universality, stability and engineering guidance value, and provides a theoretical basis and a simulation platform for the design and optimization of the disaster prevention structure in the high mountain area.

[0148] The above embodiments only express the specific implementation of the present application, which is described in more detail and specifically, but cannot be understood as the limitation of the patent scope of the present application. It should be noted that for ordinary skilled in the art, without departing from the concept of the present application, several modifications and improvements can be made, which are within the protection scope of the present application.

Claims

1. A numerical simulation method for avalanche impact protection structures based on bidirectional coupling, characterized in that, Includes the following steps: Step 1, Avalanche Source Modeling and Initialization: The discrete element method is used to simulate the avalanche process. First, the basic physical parameters of the avalanche particles are obtained; then, the particle swarm is initialized within the geometric boundary through a random particle generation and filling algorithm. An initial velocity field is set up, or a velocity field is formed by the gravity of the slope to realize the dynamic simulation and reproduction of the entire process of avalanche disaster from initiation to movement and impact on the protective structure; Step 2, Protection Structure Modeling and Parameter Setting: The avalanche protection structure is modeled based on the nonlinear finite element method. For different types of structural components, appropriate solid elements, shell elements, beam elements, and cable elements are selected to generate meshes and establish finite element models. Material parameters, contact parameters, and boundary constraints of different components are set according to structural characteristics. Step 3: Setting up the bidirectional coupling of the avalanche impact protection structure: The bidirectional coupling analysis between the avalanche and the protection structure is realized through the impact contact algorithm between discrete elements and finite elements; the snow particle position is updated in real time by capturing particle displacement, velocity, and acceleration; at the same time, the position and motion information of the finite element nodes are updated through the dynamic equations, and the contact relationship between the snow particles and the protection structure is established by combining the search algorithm; on this basis, the contact force between the avalanche particles and the protection structure is updated by the physical characteristics and contact parameters of the contact position of the discrete element and finite element, and the contact state is updated to realize the bidirectional coupling between the avalanche and the protection structure. In step 3, the contact coupling between the discrete element and the finite element is achieved through a contact algorithm based on a penalty function; based on a global search algorithm, when the discrete element particle penetrates the surface of the finite element, the normal contact force is expressed as: ； in, The penetration depth between the discrete element and the finite element method. For contact stiffness, The contact damping coefficient is... The relative velocity of the contact nodes; To ensure numerical stability and physical accuracy, the stiffness of elements in the discrete element method and the finite element method is compared. Stiffness based on stability The equivalent stiffness φ of the contact interface is determined as follows: ； Among them, the stiffness of the element is 𝑘 e Determined by material parameters and element characteristics: ； in, For the penalty stiffness factor, Bulk modulus The contact area is For the volume of a solid unit, The characteristic length of the shell element; Stiffness based on stability Related to the effective quality of nodes and the overall time step: ； in, Effective node quality; Effective node quality Take the smaller of the masses of the two interacting nodes: ; Effective angular frequency The calculation formula is: ； Considering the Coulomb friction between the discrete element method and the finite element method, the tangential contact force... Contact force with normal direction The relationship is represented as: ； Among them, the coefficient of friction Exponential interpolation is performed based on the relative velocities of discrete element particles and finite element nodes to simulate the smooth transition between static and dynamic friction: ； in, The constant is used to control the transition rate of the friction coefficient; when At this time, the coefficient of friction is always equal to the static coefficient of friction; Step 4, Numerical Simulation Control and Solution: Set the integration time step, with the step range satisfying the stability condition. The total simulation duration is set based on the actual avalanche impact duration. During the simulation, the energy output function is enabled to record and monitor the changes in key energy terms such as impact kinetic energy, structural internal energy, slip energy, and damping energy. This is used to assess the energy conversion process between avalanche particles and the structure and the risk of impact damage. An explicit time integration method is used for solution, continuously outputting the stress, strain, and displacement response variables of the structure at each time step. These are then used for subsequent damage assessment, structural performance analysis, and optimization design. Step 5, Impact Response Extraction and Result Analysis: After completing the calculation simulation, extract the key response indicators of the avalanche impact protection structure system, including the maximum impact force, deformation response, energy absorption rate, nodal acceleration, and stress concentration area distribution; analyze the velocity field, impact frequency distribution, and contact force time history of avalanche particles to identify the main impact zones of the avalanche; combine the finite element results to establish the local and overall response time histories of the protection structure, and analyze the stress characteristics and failure modes of different structural components.

2. The numerical simulation method for avalanche impact protection structures based on bidirectional coupling according to claim 1, characterized in that, In step 1, the basic physical parameters of an avalanche include at least avalanche size, source distribution, bulk density, and friction coefficient.

3. The numerical simulation method for avalanche impact protection structures based on bidirectional coupling according to claim 1, characterized in that, In step 1, based on the interparticle collision forces and gravity acting on the snow particles during their motion, the equations for particle translation and rotation are as follows: ； in, Let i be the mass of the i-th particle. For particle acceleration, As an external force, For the moment of inertia, Angular acceleration, The torque applied to the i-th particle; An explicit algorithm is used to solve for the motion state of the particles: (1) Time step Acceleration and angular acceleration updates, calculation of time t. Particle acceleration With angular acceleration : ； in, , , They represent , , The particle displacement at that moment; , , They represent , , The particle rotation angle corresponding to the given moment; (2) Half step Speed ​​updates Time of the first The speed of each particle With angular velocity They are represented as follows: ； in, , They represent The velocity and angular velocity of the particles at any given moment; (3) Half step Displacement and rotation updates, calculations Particle displacement at any time With corner for: ； in, , They represent The velocity and angular velocity of the particles at any given moment; (4) Explicit dynamics algorithm start-up condition: When At time t, the initial half-step velocity of the particle is expressed by the following formula: ； in, , These serve as the starting conditions for particle velocity and angular velocity, respectively. , , , for The initial velocity, initial acceleration, initial angular velocity, and initial angular acceleration of the particle are known at any given time.

4. The numerical simulation method for avalanche impact protection structures based on bidirectional coupling according to claim 1, characterized in that, In step 2, appropriate finite element elements are selected for modeling different types of structures. Specifically, for rigid structures, finite element meshes of structural members are established using three-dimensional solid elements or shell elements; for flexible structures, nonlinear beam elements or equivalent shell elements are used for modeling, and the mechanical properties of connecting members are defined.

5. The numerical simulation method for avalanche impact protection structures based on bidirectional coupling according to claim 4, characterized in that, In step 2, after completing the modeling, the corresponding model parameters are set according to the mechanical properties of different structural materials, specifically including: (1) For steel protective structures, determine the material density, elastic modulus, Poisson's ratio, yield strength, ultimate strength, ultimate strain and damping coefficient to ensure that the material's elastoplastic behavior and energy dissipation characteristics can be accurately reflected; (2) For concrete protective structures, determine the material density, elastic modulus, Poisson's ratio, compressive strength, maximum aggregate particle size, and set the damping coefficient in combination with the structural damping characteristics to fully reflect the nonlinearity and damage evolution process of concrete; (3) For flexible structures and their connectors, set the equivalent stiffness, ultimate bearing capacity and slip performance parameters of nonlinear connections to accurately simulate the mechanical response and deformation capacity of the connectors; in addition, the variation characteristics of material parameters under the influence of environmental factors such as temperature and strain rate should also be considered to ensure the applicability and accuracy of the model.

6. The numerical simulation method for avalanche impact protection structures based on bidirectional coupling according to claim 5, characterized in that, In step 2, setting the boundary conditions of the structural model according to the actual working conditions specifically includes: setting the bottom nodes of the structure as fixed constraints, and setting contact or sliding boundaries between the structure and the foundation.

7. The numerical simulation method for avalanche impact protection structures based on bidirectional coupling according to claim 1, characterized in that, In step 2, the dynamic control equations for the avalanche protection structure simulation are as follows: ； in, , , These are the mass matrix, damping matrix, and stiffness matrix of the structure, respectively. These are the displacement vector, velocity vector, and acceleration vector of the node, respectively. This is the vector of external forces; To improve computational efficiency and stability, an explicit second-order central difference scheme is used to discretize the time integral; the acceleration and velocity of the nodes are represented as follows: ； Based on the above equation, the dynamic governing equation for the simulation of avalanche protection structures can be rewritten as follows: ； in, , They are respectively , The nodal displacement vectors at each time step; the structural displacements are solved explicitly at each time step using the given formula. It does not require matrix inversion, has high computational efficiency, and is suitable for strongly nonlinear shock situations; To initiate the explicit time integration algorithm, the following steps are required: Initiation displacement at time Initialization is performed; based on Taylor expansion, the startup conditions of the explicit dynamic algorithm and the initial half-step velocity of the structural unit are calculated by the following formula: ； in, , They are respectively The known initial position, initial velocity, and initial acceleration of the finite element nodes at any given time.

8. The numerical simulation method for avalanche impact protection structures based on bidirectional coupling according to claim 1, characterized in that, Step 5 is described in detail below: (1) First, extract the dynamic response time history data of key components of the protective structure from the finite element calculation results, including the impact deformation, peak internal force, stress and strain level and damage state of the structural components at the impact location, anchorage location, and connection location. Analyze the local yielding, failure or anchorage failure failure modes of the structure during the impact process to achieve an accurate assessment of the structural safety performance. (2) Secondly, determine the maximum impact load of the avalanche on the protective structure, extract the contact force time history curve of the entire impact process, determine the peak contact force as the maximum impact force index, and use it for structural strength verification and impact resistance performance evaluation. At the same time, verify the rationality and accuracy of the two-way coupled numerical model by comparing it with theoretical analysis and experimental data. (3) Finally, based on different working conditions, sensitivity analysis was carried out. By comparing the differences in structural response under each working condition, the key control parameters affecting the structural impact response were identified. Combining avalanche impact force, protective structure deformation and stress concentration area distribution, local and overall response time history models of the protective structure were constructed. The stress characteristics and potential failure mechanisms of each structural component were analyzed in depth, providing a scientific basis for structural optimization design.

Images