Warehouse management model determination method and system and storage medium

Abstract

The invention relates to the technical field of finite element analysis, in particular to a warehouse management model determination method and system and a storage medium. The method comprises the following steps: constructing grid model data of a cargo stack; determining material attribute data of the cargo; collecting warehouse temperature and humidity data; loading the warehouse temperature and humidity data to the grid model data, and performing temperature and humidity stress field calculation according to the material attribute data to obtain initial stress field data; carrying out transportation vibration and strain analysis on the cargo stack according to the initial stress field data to obtain contact force and creep strain data and a stress-strain tensor field; performing frequency vibration mode analysis on the grid model data and the material attribute data to obtain frequency vibration mode data; and calculating the energy consumption damping ratio according to the frequency vibration mode data and the stress-strain tensor field to obtain damping ratio data. According to the invention, through the refined analysis and management of the dynamic stability of the goods stack in the high-humidity and variable-temperature environment, the scientificity and reliability of warehouse management are improved.

Description

Technical Field

[0001] The present invention relates to the technical field of finite element analysis, and particularly to a method, a system and a storage medium for determining a warehousing management model. Background Art

[0002] Under special environmental conditions, such as high-temperature and high-humidity areas, cold-chain warehousing, and warehouses near transportation hubs with frequent vibrations, the stability of goods is a major problem in warehousing management.

[0003] Most traditional stack stability analyses are based on static mechanics principles, regarding the stack as a rigid body and ignoring dynamic factors such as vibrations, impacts during transportation, and creep of materials. Static analysis cannot predict the resonance phenomenon of the stack under vibrations at specific frequencies, and resonance can lead to the rapid destruction of the stack. Creep refers to the phenomenon that the strain of materials slowly increases with time under a constant stress. Many packaging materials (such as cardboard boxes, plastic films) have obvious creep characteristics, and long-term stacking can cause the stack to deform or even collapse. Static analysis cannot consider this time-dependence.

[0004] Existing methods use empirical formulas or simple linear models for the influence of temperature and humidity, and cannot accurately reflect the complex influence of temperature and humidity changes on stack stability. The temperature and humidity in a warehouse are usually unevenly distributed, and the goods at different positions are affected differently by temperature and humidity. Existing methods often use an average value or a single value to represent the temperature and humidity of the entire warehouse, ignoring this spatial difference. Temperature and humidity changes not only cause thermal expansion and contraction and hygroscopic expansion and contraction, but also affect the mechanical properties of materials (such as elastic modulus, yield strength, etc.). Existing methods rarely consider this coupling relationship between temperature and humidity and mechanical properties.

[0005] In summary, there are problems of the limitations of static analysis and the incomplete consideration of environmental factors in the existing technology that need to be solved urgently. Summary of the Invention

[0006] Based on this, it is necessary to provide a method, a system and a storage medium for determining a warehousing management model to solve at least one of the above technical problems.

[0007] To achieve the above object, a method for determining a warehousing management model includes the following steps:

[0008] Step S1: Construct grid model data of the stack; determine material property data of the goods; collect warehouse temperature and humidity data; load the warehouse temperature and humidity data into the grid model data, calculate the temperature and humidity stress field according to the material property data to obtain initial stress field data; perform transportation vibration and strain analysis on the stack according to the initial stress field data to obtain contact force and creep strain data and stress-strain tensor field;

[0009] Step S2: Conduct frequency vibration mode analysis on the grid model data and material property data to obtain frequency vibration mode data; calculate the energy consumption damping ratio based on the frequency vibration mode data and the stress-strain tensor field to obtain damping ratio data; integrate the damping ratio data and the frequency vibration mode data into a modal parameter set;

[0010] Step S3: Conduct instability critical domain analysis on the pallet according to the modal parameter set to obtain instability critical domain data; predict the long-term stability of the creep influence factor based on the instability critical domain data and the contact force and creep strain data to obtain a stability margin factor;

[0011] Step S4: Set storage adjustment rules for different risk levels according to the stability margin factor, and optimize the real-time state storage strategy according to the storage adjustment rules to obtain a storage optimization plan.

[0012] The present invention constructs a refined finite element model of the pallet stack, and comprehensively considers the temperature and humidity environment in the warehouse, material properties, transportation vibration, and contact and creep effects, obtaining comprehensive data including the initial stress field, contact force and creep strain, and stress-strain tensor field. These data provide accurate initial conditions and boundary conditions for subsequent analysis, overcome the limitations of simplified models, ignored environmental factors and dynamic effects in traditional methods, significantly improve the accuracy and reliability of the analysis, and lay a solid foundation for the stability assessment of the storage pallet stack. Based on the mesh model, material properties, and stress-strain tensor field data obtained in step S1, frequency vibration mode analysis and energy consumption damping ratio calculation are carried out to obtain a set of modal parameters including natural frequencies, vibration modes, and damping ratios. These parameters reveal the inherent dynamic characteristics of the pallet stack and the energy dissipation mechanism during vibration. Compared with the empirical damping model used in traditional methods, the set of modal parameters obtained by the refined analysis and calculation in this step is more accurate and reliable, providing a more scientific basis for subsequent stability assessment and prediction. Based on the set of modal parameters, instability critical domain analysis is carried out to obtain instability critical domain data such as dynamic stability factors, critical inclination angles, and critical slip distances; and further considering the creep effect, long-term stability prediction is carried out to obtain a stability margin factor. This step combines the dynamic characteristics of the pallet stack, various instability modes such as tipping and sliding, and the creep effect of the material, realizing a comprehensive, dynamic, and long-term assessment of the stability of the pallet stack, overcoming the limitations of traditional static analysis methods, and providing a more reliable decision-making basis for warehouse management. According to the stability margin factor, risk levels are divided, and corresponding warehouse adjustment rules are formulated. By real-time monitoring temperature and humidity data and dynamically adjusting the warehouse strategy, a warehouse optimization plan is finally obtained. This step combines the theoretical analysis results with actual warehouse management, realizes the quantitative assessment and dynamic regulation of risks, can effectively prevent the occurrence of pallet stack instability accidents, improves the safety and management efficiency of the warehouse, reduces losses caused by damaged goods and safety accidents, and has significant economic and social benefits. Therefore, the present invention provides a method for determining a warehouse management model, which fundamentally solves the limitations of traditional static analysis and the problem of incomplete consideration of environmental factors through refined analysis and management of the dynamic stability of pallet stacks in high-humidity and variable-temperature environments, providing a more scientific and reliable basis for warehouse management.

[0013] Preferably, the temperature and humidity stress field calculation in step S1 includes:

[0014] Construct a continuous function of the temperature and humidity field according to the temperature and humidity data in the warehouse to obtain a temperature and humidity field function;

[0015] Set boundary conditions for the temperature and humidity field function to obtain a temperature and humidity field; map the temperature and humidity field to the mesh model data, and carry out temperature and humidity gradient field calculation to obtain mesh temperature and humidity data;

[0016] Calculate the thermal expansion and hygroscopic expansion strains based on the grid temperature and humidity data and material property data to obtain the thermo-hygroscopic strain field;

[0017] Calculate the initial stress field based on the thermo-hygroscopic strain field to obtain the initial stress field data.

[0018] In the present invention, by deploying a sensor network inside the warehouse to collect temperature and humidity data, and using the radial basis function interpolation method and cubic spline interpolation method to construct a temperature and humidity field function that is continuous in space and time, the limitations brought about by the traditional method of only using discrete temperature and humidity empirical values or average values are overcome, realizing a refined and dynamic description of the temperature and humidity distribution inside the warehouse, improving the accuracy and reliability of the subsequent calculation of the thermo-hygroscopic stress field, and providing a key basis for accurately analyzing the stability of the pallet stack. By setting reasonable boundary conditions (such as adiabatic, convective heat transfer, etc.) for the boundaries of the warehouse walls, floor, roof, doors, vents, etc., and accurately mapping the constructed temperature and humidity field function onto the finite element mesh model of the pallet stack, while calculating the temperature and humidity gradients of each unit, fine grid temperature and humidity data are obtained. This not only considers the influence of the warehouse environment on the temperature and humidity of the pallet stack, but also considers the non-uniform distribution of temperature and humidity inside the pallet stack, providing accurate input for the subsequent calculation of thermo-hygroscopic strain and initial stress, and further improving the analysis accuracy. Based on the fine grid temperature and humidity data and accurate material thermal expansion coefficient and hygroscopic expansion coefficient, calculate the thermal expansion strain and hygroscopic expansion strain of each unit, and superimpose them to obtain the thermo-hygroscopic strain field. This step directly correlates the temperature and humidity changes with the deformation of the pallet stack, fully considering the expansion / shrinkage characteristics of the packaging material under different temperature and humidity conditions, and providing a key basis for accurately calculating the initial stress caused by temperature and humidity. Through the finite element method, based on the thermo-hygroscopic strain field and the elastic mechanical properties of the material, calculate the initial stress field generated due to the inability of the thermo-hygroscopic strain to develop freely. This step converts the deformation caused by temperature and humidity changes into stress, obtaining the initial stress distribution of the pallet stack in the stacked state, providing real initial conditions for subsequent analyses (such as transportation vibration analysis, creep analysis), significantly improving the accuracy and reliability of the overall analysis, and avoiding the errors caused by ignoring the initial stress in the traditional method.

[0019] Preferably, the transportation vibration and strain analysis in step S1 includes:

[0020] Determine the vibration load time history curve of the transportation process, apply the vibration load time history curve to the grid model data, and perform transient dynamics analysis based on the initial stress field data to obtain the node displacement time history data;

[0021] Define the contact pairs between the pallets in the grid model data, and perform creep calculation on the contact pairs according to the material property data to obtain the contact force and creep strain data;

[0022] Based on the node displacement time history data, material property data, and contact force and creep strain data, the stress-strain tensor field is calculated to obtain the stress-strain tensor field.

[0023] In the embodiment of the present invention, vibration acceleration data is collected by installing sensors on an actual transport vehicle and converted into a vibration load time history curve, which is applied to the bottom of the pallet finite element model. Combining with the initial stress field data calculated previously, transient dynamic analysis is carried out to obtain the node displacement time history data of the pallet during transportation vibration. This method truly simulates the dynamic excitation received by the pallet during transportation and takes into account the influence of the initial stress, overcoming the limitation of traditional static analysis that ignores the vibration effect, and providing key data for accurately evaluating the dynamic response and stability of the pallet. By defining contact pairs in the pallet finite element model and using the penalty function method to simulate the contact behavior, and at the same time combining the creep model of the material (such as the Burgers model) to calculate the creep strain increment on the contact surface, the contact force and creep strain data are obtained. This step fully considers the mechanical behavior of the contact surface between pallets and the influence of the creep characteristics of the packaging material on the contact force and contact deformation, providing an important basis for accurately analyzing the stability of the pallet during long-term stacking and vibration. Based on the finite element theory, the node displacement time history data, material property data (elastic modulus, Poisson's ratio, etc.), and contact force and creep strain data are combined to calculate the strain tensor and stress tensor of each element at each time step, and the stress-strain tensor field data is formed by combination. This step comprehensively considers the influence of elastic deformation, thermo-hygroscopic strain, creep strain, and contact force of the pallet during transportation vibration, obtaining the detailed stress and strain states inside the pallet changing with time, providing a complete and accurate data basis for subsequent modal analysis, damping ratio calculation, and stability evaluation.

[0024] Preferably, the frequency and mode analysis in step S2 includes:

[0025] Extract the stiffness matrix and mass matrix from the mesh model data and material property data to obtain the stiffness and mass matrix;

[0026] Construct the natural frequency characteristic equation based on the stiffness and mass matrix to obtain the natural frequency characteristic equation;

[0027] Calculate the natural frequency and mode based on the natural frequency characteristic equation to obtain the frequency and mode data.

[0028] Based on the finite element mesh model data and material property data of the pallet stack, the present invention calculates the stiffness matrix and mass matrix of each element, and assembles them into the overall stiffness matrix K and the overall mass matrix M. This step links the geometric shape and material properties of the pallet stack with its mechanical behavior, providing the key system matrices for subsequent modal analysis and serving as the basis for frequency vibration mode analysis. Accurate stiffness and mass matrices are the prerequisite for obtaining reliable natural frequencies and vibration modes. Based on the obtained overall stiffness matrix K and overall mass matrix M, the free vibration characteristic equation of the pallet stack |K - ω 2 M| = 0 is constructed. This step transforms the vibration characteristics of the pallet stack into a mathematical equation, laying the foundation for solving the natural frequencies and vibration modes. The characteristic equation accurately describes the mechanical behavior of the pallet stack under the condition of undamped free vibration. Numerical methods such as the subspace iteration method are used to solve the characteristic equation to obtain the natural frequencies of the pallet stack and the corresponding normalized vibration mode vectors. This step reveals the inherent vibration characteristics of the pallet stack, namely at which frequencies resonance is likely to occur and the vibration patterns of the pallet stack at these frequencies. These frequency vibration mode data provide key reference bases for subsequent dynamic stability analysis, damping ratio calculation, and warehousing strategy optimization.

[0029] Preferably, the calculation of the energy consumption damping ratio in step S2 includes:

[0030] Select a typical time period for the stress-strain tensor field to obtain a modal time series;

[0031] Calculate the strain energy time history of the stress-strain tensor field according to the modal time series to obtain a strain energy curve;

[0032] Calculate the material damping dissipation energy according to the material property data and the stress-strain tensor field to obtain the material dissipation energy;

[0033] Calculate the contact friction dissipation energy according to the modal time series and the contact force and creep strain data to obtain the friction dissipation energy;

[0034] Calculate the creep dissipation energy according to the stress-strain tensor field, the modal time series, and the contact force and creep strain data to obtain the creep dissipation energy;

[0035] Conduct a hysteresis loop area analysis according to the modal time series, the material dissipation energy, the friction dissipation energy, and the creep dissipation energy to obtain the total dissipation energy;

[0036] Conduct modal energy distribution and modal dissipation calculation according to the total dissipation energy and the frequency vibration mode data to obtain modal dissipation distribution data;

[0037] Calculate the damping ratio according to the modal dissipation distribution data, the strain energy curve, and the frequency vibration mode data under the principle of energy balance to obtain the damping ratio data.

[0038] By focusing on these time periods, the analysis concentrates on the inherent damping characteristics of the system, avoiding the influence of external excitations such as initial transportation vibrations. This makes the estimation of the damping ratio more accurate and reliable, which is crucial for predicting long-term stability and the response to future disturbances. The selection criteria (based on natural frequency and minimum number of cycles) ensure that the selected data can represent the modal behavior. By focusing on these time periods, the analysis concentrates on the inherent damping characteristics of the system, avoiding the influence of external excitations such as initial transportation vibrations. This makes the estimation of the damping ratio more accurate and reliable, which is crucial for predicting long-term stability and the response to future disturbances. The selection criteria (based on natural frequency and minimum number of cycles) ensure that the selected data can represent the modal behavior. The energy dissipated due to the viscoelastic properties inherent in the material (specifically corrugated cardboard in this example) is quantified. By using a viscoelastic model (calibrated with DMA test data) and integrating the damping power over time, this calculation accurately captures the energy loss caused by internal friction within the material itself. This is a key component of the total damping and is usually an important contributor, especially in materials like cardboard. By considering the normal and tangential contact forces and the relative sliding velocity, this calculation accurately determines the energy lost through frictional sliding. This is another key component of the total damping and is particularly important in systems with many contact interfaces such as a stack of boxes. It accounts for dissipation mechanisms overlooked in the material damping calculation. By incorporating creep models (the Burgers model for cardboard, the power-law model for plastic films) and integrating the creep power over time, this calculation captures the energy loss associated with slow, continuous deformation. This is particularly important for long-term stability assessment as creep can significantly alter the stress distribution and stability of the pallet stack. It adds a time-dependent aspect to the energy dissipation. By plotting the stress-strain hysteresis curve and calculating the area enclosed by it, this analysis directly measures the energy loss per cycle. This is a crucial check on the sum of the separately calculated dissipation components (material, friction, creep). The consistency between the hysteresis loop area and the sum of the components enhances confidence in the accuracy of the damping calculation. This is a way to close the energy balance loop. By using the modal participation factor, the calculation correctly accounts for the relative contribution of each mode to the overall vibration. This is crucial because different modes may have different damping characteristics. This step provides the modal energy dissipation (ΔWi), which is necessary for calculating the modal damping ratio. The damping ratio (ζi) of each vibration mode is calculated using the modal energy dissipation (ΔWi) and the maximum strain energy (Ui) of that mode. The damping ratio is a dimensionless parameter used to quantify the rate of energy dissipation in each mode. These damping ratios are important input parameters for dynamic simulations and stability analyses and can accurately predict the response of the pallet stack to various disturbances. It is the final outcome of the entire process and provides the key output for subsequent stability calculations.

[0039] Preferably, the instability critical region analysis in step S3 includes:

[0040] Construct a dynamic stability criterion based on the modal parameter set to obtain a dynamic stability factor;

[0041] Obtain the geometric information of the pallet stack; establish a tipping instability model based on the geometric information of the pallet stack to obtain a tipping instability model;

[0042] Calculate the critical inclination angle of the tipping instability model based on the contact force and creep strain data to obtain the critical inclination angle;

[0043] Establish a slip instability model based on the contact force, creep strain data, and material property data to obtain a slip instability model;

[0044] Calculate the critical slip distance of the slip instability model to obtain the critical slip distance.

[0045] In the present invention, the vibration characteristics (natural frequency and damping ratio) of the pallet stack are integrated into a single, quantitative index - the dynamic stability factor (D). This factor provides a method for quickly evaluating the ability of the pallet stack to resist dynamic disturbances (such as vibrations during transportation). By comparing with a preset threshold, it can be determined whether the pallet stack is in a dynamically stable state, avoiding complex dynamic analysis. This criterion takes into account two unfavorable situations: too low frequency (prone to resonance) and too small damping (slow vibration attenuation). By obtaining key information such as the geometric dimensions, stacking method, and center of gravity position of the pallet stack and simplifying it into a rigid body model, a moment balance equation can be established to analyze the force situation of the pallet stack in an inclined state. This model simplifies complex physical phenomena into computable mechanical problems, making it possible to quantitatively evaluate the tipping risk. By considering the friction coefficient, normal pressure distribution, and stress redistribution caused by creep on the contact surface, this calculation can more accurately reflect the actual situation. The critical inclination angle is an intuitive index that can be used to determine whether the pallet stack will tip under given conditions, providing a direct decision-making basis for warehouse management. By considering the normal force, tangential force (friction force), and static friction coefficient at the bottom of the pallet stack, a force balance equation can be established to analyze the force situation of the pallet stack in the horizontal direction. This model takes into account the relationship between the driving force (such as horizontal inertial force) and the anti-slip force, making it possible to quantitatively evaluate the slip risk. By considering the normal force distribution, friction coefficient distribution, and the influence of creep on the contact surface, this calculation can more accurately reflect the actual situation. The critical slip distance is an intuitive index that can be used to determine whether the pallet stack will slip under a given horizontal driving force, providing a direct decision-making basis for warehouse management. It takes into account the non-uniformity of the contact surface and the creep effect, making the calculation result closer to the actual situation.

[0046] Preferably, the long-term stability prediction of the creep influence factor in step S3 includes:

[0047] Obtain the data of the instability critical domain, where the data of the instability critical domain includes the dynamic stability factor, the critical inclination angle, and the critical slip distance; establish a creep constitutive model and calibrate the parameters according to the material property data to obtain a creep model;

[0048] Determine the time discretization strategy according to the creep model to obtain a time discretization scheme;

[0049] Calculate the creep strain increment according to the time discretization scheme, the creep model, and the contact force and creep strain data to obtain creep evolution data;

[0050] Update the geometry and contact characteristics of the stack according to the creep evolution data to obtain the deformed state data;

[0051] Conduct a rebalancing analysis based on the deformed state data and perform stress redistribution calculation to obtain a rebalanced stress field;

[0052] Recalculate the dynamic stability factor according to the rebalanced stress field, the deformed state data, and the dynamic stability factor to obtain the dynamic stability time history;

[0053] Recalculate the critical inclination angle according to the deformed state data, the rebalanced stress field, and the critical inclination angle to obtain the inclination angle time history;

[0054] Recalculate the critical slip distance according to the deformed state data, the rebalanced stress field, and the critical slip distance to obtain the slip threshold time history;

[0055] Integrate the time-varying stability parameters of the dynamic stability time history, the inclination angle time history, and the slip threshold time history to obtain the time-varying stability parameters; calculate the stability margin factor according to the time-varying stability parameters to obtain the stability margin factor.

[0056] By establishing creep constitutive models (such as the Burgers model for corrugated cardboard and the power-law model for plastic film) and calibrating the parameters, the present invention can accurately describe the deformation behavior of materials under long-term loads, which is the key to predicting the creep effect. Since the creep process is time-dependent, it is necessary to discretize the time domain. Adopting non-uniform time steps (small at the beginning and gradually increasing later) can reduce the computational amount while ensuring accuracy. The implicit time integration method improves numerical stability and is suitable for long-term creep analysis. By calculating the creep strain increment according to the current stress state and creep model at each time step and accumulating it to the total creep strain, the creep deformation process of the pallet can be traced. The creep evolution data provides the basis for subsequent geometry update and rebalancing analysis. By updating the node coordinates according to the creep strain, the deformation of the pallet can be simulated. At the same time, updating the normal force and tangential force on the contact surface can reflect the stress redistribution caused by creep. The deformation state data provides the basis for subsequent rebalancing analysis and stability assessment. Since creep will cause changes in the stress state inside the pallet, it is necessary to perform a static rebalancing analysis to obtain the stress field under the new equilibrium state. The rebalanced stress field is the key input for subsequent re-evaluating stability. By re-performing the modal analysis using the updated geometry and rebalanced stress field, the new natural frequency and damping ratio can be calculated, and then the new dynamic stability factor can be obtained. The dynamic stability time history D(t) reflects the change trend of the dynamic stability of the pallet over time and provides an important basis for long-term stability prediction. By re-calculating the critical inclination angle using the updated geometry and rebalanced stress field, the change of the critical inclination angle over time can be obtained. The inclination angle time history θc(t) reflects the change trend of the anti-overturning ability of the pallet over time and provides an important basis for long-term stability prediction. By re-calculating the critical slip distance using the updated geometry and rebalanced stress field, the change of the critical slip distance over time can be obtained. The slip threshold time history (t) reflects the change trend of the anti-slip ability of the pallet over time and provides an important basis for long-term stability prediction. The dynamic stability time history, inclination angle time history, and slip threshold time history are integrated into a single, comprehensive index - the stability margin factor S(t). This factor comprehensively reflects the stability change of the pallet under creep action and provides the stability margin compared with the initial state (or safety threshold). The stability margin factor S(t) is an intuitive and easy-to-understand index that can be used to judge whether there is a risk of instability of the pallet at any time and provides the final decision-making basis for warehouse management.

[0057] Preferably, step S4 includes the following steps:

[0058] Step S41: Divide the risk levels according to the stability margin factor to obtain a risk level threshold table;

[0059] Step S42: Perform risk identification based on the risk level threshold table and the stability margin factor to obtain a risk distribution map;

[0060] Step S43: Formulate warehousing strategy adjustment rules according to the risk level threshold table to obtain warehousing adjustment rules;

[0061] Step S44: Obtain real-time warehouse temperature and humidity data; perform dynamic adjustment and optimization iteration based on the real-time warehouse temperature and humidity data, the risk distribution map, the stability margin factor, and the warehousing adjustment rules to obtain an optimized risk map;

[0062] Step S45: Generate a warehousing optimization plan according to the optimized risk map to obtain a warehousing optimization plan.

[0063] In the present invention, the value of the stability margin factor S(t) is converted into an easily understandable risk level (safe, low risk, medium risk, high risk), and a risk level threshold table is established. This enables warehouse managers to quickly and intuitively evaluate the stability state of the pallet stack without in-depth understanding of complex mechanical calculations. The risk level threshold table provides a unified standard and basis for subsequent risk assessment and warehousing strategy adjustment. Clear, specific, and operable adjustment rules are formulated for pallet stacks of different risk levels. These rules cover different levels of intervention measures from strengthening monitoring to immediate unloading and quantify them into operable parameters (such as maximum stacking height, temperature and humidity control range, etc.). The warehousing adjustment rule library provides an action guide for warehouse managers to deal with different risk situations, ensuring the effectiveness and consistency of adjustment measures. By real-time monitoring the warehouse temperature and humidity data and feeding it back into the model, the stability margin factor and risk level of the pallet stack can be continuously updated. According to the latest risk distribution map and warehousing adjustment rules, iterative adjustment is performed until the risk levels of all pallet stacks are within an acceptable range. This ensures that the warehousing strategy can adapt to environmental changes and the long-term creep effect of the pallet stack, maintaining the long-term stability of the warehouse. By providing a visual interface (optimized risk map), managers can clearly understand the optimized warehousing layout and management requirements. The final warehousing optimization plan aims to improve the safety and management efficiency of the warehouse and reduce the risk of pallet stack instability.

[0064] Preferably, the present invention also provides a warehousing management model determination system for executing the above-mentioned warehousing management model determination method. The warehousing management model determination system includes:

[0065] The multi-field coupling characterization module is used to construct the grid model data of the pallet stack; determine the material property data of the goods; collect the temperature and humidity data of the warehouse; load the temperature and humidity data of the warehouse into the grid model data, calculate the temperature and humidity stress field according to the material property data, and obtain the initial stress field data; perform transportation vibration and strain analysis on the pallet stack according to the initial stress field data, and obtain the contact force, creep strain data, and stress-strain tensor field;

[0066] The modal damping identification module is used to perform frequency and vibration mode analysis on the grid model data and material property data to obtain frequency and vibration mode data; calculate the energy consumption damping ratio according to the frequency and vibration mode data and the stress-strain tensor field to obtain damping ratio data; integrate the damping ratio data and the frequency and vibration mode data into a modal parameter set;

[0067] The stability margin evaluation module is used to perform instability critical domain analysis on the pallet stack according to the modal parameter set to obtain instability critical domain data; perform long-term stability prediction of the creep influence factor according to the instability critical domain data and the contact force and creep strain data to obtain the stability margin factor;

[0068] The warehousing strategy optimization module is used to set warehousing adjustment rules under different risk levels according to the stability margin factor, and perform real-time state warehousing strategy optimization according to the warehousing adjustment rules to obtain a warehousing optimization plan.

[0069] Preferably, a computer-readable storage medium stores a computer program, and when the computer program is executed, the above-mentioned method for determining a warehousing management model is implemented. Brief Description of the Drawings

[0070] Figure 1 It is a schematic flow chart of the steps of a method for determining a warehousing management model;

[0071] Figure 2 It is a schematic detailed implementation step flow chart of step S4 in the present invention.

[0072] The realization of the purpose, functional features, and advantages of the present invention will be further described in conjunction with the embodiments with reference to the drawings. Detailed Embodiments

[0073] The technical method of the present invention will be clearly and completely described below with reference to the drawings. Obviously, the described embodiments are some, but not all, of the embodiments of the present invention. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present invention without creative efforts shall fall within the protection scope of the present invention.

[0074] In addition, the drawings are only schematic illustrations of the present invention and are not necessarily drawn to scale. The same reference numerals in the drawings denote the same or similar parts, and thus repeated descriptions thereof will be omitted. Some of the block diagrams shown in the drawings are functional entities and do not necessarily correspond to physically or logically independent entities. The functional entities may be implemented in software form, or implemented in one or more hardware modules or integrated circuits, or implemented in different networks and / or processor methods and / or microcontroller methods.

[0075] It should be understood that although terms such as "first" and "second" may be used herein to describe various units, these units should not be limited by these terms. These terms are only used to distinguish one unit from another. For example, without departing from the scope of the exemplary embodiments, the first unit may be referred to as the second unit, and similarly the second unit may be referred to as the first unit. The term "and / or" used herein includes any and all combinations of one or more of the listed related items.

[0076] In the embodiments of the present invention, with reference to Figure 1 as shown, it is a schematic flow chart of the steps of the method for determining the warehousing management model of the present invention. In this example, the method for determining the warehousing management model includes the following steps:

[0077] Step S1: Construct the grid model data of the pallet; determine the material property data of the goods; collect the temperature and humidity data of the warehouse; load the temperature and humidity data of the warehouse into the grid model data, calculate the thermo-hygro stress field according to the material property data, and obtain the initial stress field data; perform transportation vibration and strain analysis on the pallet according to the initial stress field data, and obtain the contact force and creep strain data and the stress-strain tensor field;

[0078] In the embodiments of the present invention, first, use CAD software to construct a three-dimensional geometric model of the pallet and perform mesh division to generate finite element mesh model data; then, query or experimentally determine the physical properties (density, elastic modulus, Poisson's ratio, coefficient of thermal expansion, coefficient of hygroscopic expansion, creep parameters, etc.) of the goods and packaging materials, and associate them with the corresponding units in the mesh model to form a material property database; then, collect the temperature and humidity data in the warehouse through sensors, use the radial basis function interpolation method and the cubic spline interpolation method to construct a continuous temperature and humidity field function, map it onto the mesh model, calculate the temperature and humidity gradient, and combine the material property data to calculate the thermo-hygro strain, and then solve the equilibrium equation by the finite element method to obtain the initial stress field data; finally, apply the measured transportation vibration acceleration data to the bottom of the mesh model, perform transient dynamics analysis, and at the same time consider the viscoelastic characteristics of the material and the contact friction between the pallets to calculate the contact force and creep strain data and the stress-strain tensor field.

[0079] Step S2: Perform frequency vibration mode analysis on the mesh model data and material property data to obtain frequency vibration mode data; calculate the energy consumption damping ratio based on the frequency vibration mode data and the stress-strain tensor field to obtain damping ratio data; integrate the damping ratio data and the frequency vibration mode data into a modal parameter set;

[0080] In the embodiment of the present invention, based on the finite element theory, the element stiffness matrix and mass matrix are extracted from the mesh model data and the material property database, and assembled into the overall stiffness matrix K and the overall mass matrix M to construct the characteristic equation |K - ω 2 M| = 0; solve the characteristic equation by the subspace iteration method to obtain the natural frequency vector ω and the normalized vibration mode matrix Φ of the pallet; select the typical time period of vibration attenuation in the stress-strain tensor field data, calculate the change in strain energy of each element in each time period, and calculate the material damping dissipation energy, contact friction dissipation energy, and creep dissipation energy based on the material viscoelastic model, contact friction model, and creep model respectively. Verify the total energy dissipation through the analysis of the hysteresis loop area, and distribute the total energy dissipation to each mode according to the modal participation factor. Finally, calculate the damping ratio of each mode based on the energy balance principle to form damping ratio data; combine the natural frequency, vibration mode, and damping ratio into a modal parameter set.

[0081] Step S3: Perform instability critical domain analysis on the pallet according to the modal parameter set to obtain instability critical domain data; perform long-term stability prediction of the creep influence factor based on the instability critical domain data and the contact force and creep strain data to obtain a stability margin factor;

[0082] In the embodiment of the present invention, based on the natural frequency and damping ratio in the modal parameter set, construct a dynamic stability criterion and calculate the dynamic stability factor D; establish a mechanical model for the tipping instability of the pallet, and calculate the critical inclination angle θc according to the contact force and creep strain data; establish a mechanical model for the sliding instability of the pallet, and calculate the critical sliding distance according to the contact force and creep strain data and the material property data; establish a creep constitutive model (such as the Burgers model, power-law model) based on the material property data and perform parameter calibration; through the time discretization strategy, calculate the creep strain increment in each time step and accumulate to obtain the creep evolution data; update the pallet geometry and contact characteristics according to the creep evolution data, perform rebalancing analysis to obtain the rebalancing stress field; based on the updated state, recalculate the dynamic stability factor D(t), critical inclination angle θc(t), and critical sliding distance (t), and integrate them into time-varying stability parameters; finally, calculate the stability margin factor S(t) according to the time-varying stability parameters.

[0083] Step S4: Set the warehousing adjustment rules under different risk levels according to the stability margin factor, and optimize the real-time state warehousing strategy according to the warehousing adjustment rules to obtain a warehousing optimization plan.

[0084] In the embodiments of the present invention, according to the numerical range of the stability margin factor S(t) and practical experience, the stability risks of the stacks are divided into four levels: safe, low risk, medium risk, and high risk, and the threshold intervals of the stability margin factor corresponding to each level are set to form a risk level threshold table; according to the existing warehousing scheme and the model calculation results of steps S1 - S3, the risk level of each stack is evaluated, and a risk distribution map is generated; for different risk levels, corresponding warehousing strategy adjustment rules are formulated and quantified into operable parameters to form a warehousing adjustment rule library; according to the real-time temperature and humidity data, the risk distribution map is dynamically updated, and adjustments and optimization iterations are carried out according to the warehousing adjustment rule library until the risk levels of all stacks are within an acceptable range; finally, the final warehousing strategy is sorted into a document or a data table, and a visual interface is provided to form a warehousing optimization scheme.

[0085] Preferably, the temperature and humidity stress field calculation in step S1 includes:

[0086] Construct a continuous function of the temperature and humidity field based on the warehouse temperature and humidity data to obtain the temperature and humidity field function;

[0087] Set the boundary conditions for the temperature and humidity field function to obtain the temperature and humidity field; map the temperature and humidity field to the grid model data, calculate the temperature and humidity gradient field, and obtain the grid temperature and humidity data;

[0088] Calculate the thermal expansion and hygroscopic expansion strains based on the grid temperature and humidity data and the material property data to obtain the thermo-hygroscopic strain field;

[0089] Calculate the initial stress field based on the thermo-hygroscopic strain field to obtain the initial stress field data.

[0090] In the embodiments of the present invention, along the three directions of length, width, and height in the warehouse, a temperature and humidity sensor is arranged every 5 meters, with a total of 120 sensors arranged. The sensors collect data once per minute and transmit the data to the central server through a wireless network. The original data received by the server is discrete data points, including the sensor number (ID), spatial coordinates (x, y, z), timestamp (t), temperature value (T), and relative humidity value (RH). Preprocess these original data, including: eliminating outliers caused by sensor failures or communication errors (for example, data points with temperature values outside the range of -20°C to 60°C, or relative humidity values outside the range of 0% to 100%); interpolating missing data using the average value of adjacent time point data. The preprocessed data is used to construct a continuous function of the temperature and humidity field. The radial basis function (RBF) interpolation method is used to construct the temperature and humidity field in the spatial dimension. The multiquadric function is selected as the radial basis function, and its expression is: φ(r) = (r 2 +c2 )^(1 / 2), where r is the Euclidean distance between two points in space, and c is the shape parameter, which is determined to be 0.5 by the cross-validation method. For each moment t, the temperature field function T(x, y, z, t) is expressed as: T(x, y, z, t) = Σλ i (t) × φ(||P - P i ||), where P = (x, y, z) is the spatial coordinate of the point to be determined, and P i is the spatial coordinate of the i-th sensor, and λ i (t) is the weight coefficient corresponding to this sensor, which is obtained by solving the linear equation system Aλ = b, where A ij = φ(||P i - P j ||), and b is the vector of temperature measurement values of each sensor at this moment. The construction method of the relative humidity field function RH(x, y, z, t) is the same as that of the temperature field function. In the time dimension, the cubic spline interpolation method is used to fit the curves of temperature and relative humidity changing with time at each spatial point (x, y, z) with piecewise cubic polynomials, ensuring that the function values and the first and second derivatives are continuous at the connection points of adjacent time periods. The finally obtained temperature and humidity field functions T(x, y, z, t) and RH(x, y, z, t) can accurately describe the temperature and humidity values at any position and any moment in the warehouse.

[0091] The walls, floor, and roof of the warehouse are set as adiabatic boundaries, that is, on these boundaries, the normal gradient of the temperature and humidity field is zero. For the warehouse door, it is set as a convective heat transfer boundary, and its boundary condition equation is: where k is the thermal conductivity of the wall material (taking 0.5 W / (m·K)), is the gradient of the temperature along the boundary normal direction, h is the convective heat transfer coefficient (taking 10 W / (m 2 ·K)), T is the boundary temperature, and T extis the external environmental temperature (measured by sensors installed outside the warehouse). For the ventilation openings inside the warehouse, the air flow rate through the ventilation openings is calculated based on the size and wind speed of the ventilation openings, and is applied as a boundary condition to the temperature and humidity field function. The constructed temperature and humidity field function is mapped onto the finite element mesh model of the pallet. For each node of the mesh model, according to its spatial coordinates (x, y, z) and the current time t, the temperature value and relative humidity value of the node are directly calculated by substituting them into the temperature and humidity field functions T(x, y, z, t) and RH(x, y, z, t). For each hexahedron element, the average value of the temperature values of its 8 nodes is calculated as the average temperature of the element, and the average value of the relative humidity values of the 8 nodes is calculated as the average relative humidity of the element. For tetrahedron elements, the average value of 4 nodes is calculated. Calculate the temperature and humidity gradient fields. For each element, the central difference method is used to calculate the temperature gradient and relative humidity gradient. Taking the hexahedron element as an example, the formula for calculating the component of the temperature gradient in the x direction is: where T1 to T8 are the temperature values of the 8 nodes of the element, and Δx is the length of the element in the x direction. The temperature gradients in the y and z directions, and the components of the relative humidity gradient in the three directions are calculated using a similar method. Finally, the temperature gradient vector and the relative humidity gradient vector of each element are obtained. Store these data as grid temperature and humidity data.

[0092] The packaging material of the goods is corrugated cardboard, and its coefficient of thermal expansion α is 1.2×10 -5 / °C, and the coefficient of hygroscopic expansion β is 2.5×10 -4 / %RH. According to the grid temperature and humidity data, calculate the temperature change ΔT and relative humidity change ΔRH of each element. The temperature change is the difference between the current temperature and the reference temperature (20°C), and the relative humidity change is the difference between the current relative humidity and the reference relative humidity (50%RH). Calculate the thermal expansion strain. The corrugated cardboard is regarded as an isotropic material, and its thermal expansion strain ε th is the three normal strain components, and their values are equal: ε th xx＝ ε th yy = ε th zz = α×ΔT. The shear strain components are zero: ε th xy = ε th yz = ε th zx = 0. Calculate the hygroscopic expansion strain. The hygroscopic expansion strain ε hy of the corrugated cardboard is also the three normal strain components, and their values are equal: ε hy xx = ε hy yy = ε hyzz = β × ΔRH. The shear strain components are zero: ε th xy = ε hy yz = ε hy zx = 0. Superimpose the thermal expansion strain and the hygroscopic expansion strain of each element to obtain the thermo - hygroscopic strain ε th+hy .

[0093] The finite element method is used to calculate the initial stress field. The elastic modulus E of the corrugated board is 3 GPa, and the Poisson's ratio ν is 0.2. According to the theory of elasticity, the stress - strain relationship is: σ = D × (ε - ε th+hy ), where σ is the stress tensor, ε is the total strain tensor, ε th+hy is the thermo - hygroscopic strain tensor, and D is the elastic matrix.

[0094] Since the stack is subject to its own gravity and constraints in the stacking state, the thermo - hygroscopic strain cannot develop freely, thus generating an initial stress. By solving the finite element equilibrium equation: K × u = F, the nodal displacement vector u is obtained, where K is the global stiffness matrix and F is the equivalent nodal force vector. The equivalent nodal force vector F is caused by the thermo - hygroscopic strain, and its calculation formula is: F = ∫B T × D × ε th+ty dV, where B is the strain - displacement matrix, and the integration is carried out over the entire volume of the stack.

[0095] After obtaining the nodal displacement vector u by solving, according to the geometric equation ε = B × u, calculate the total strain ε of each element. Then, according to the stress - strain relationship σ = D × (ε - ε th+ty ), calculate the stress tensor σ of each element. The stress tensor σ contains six components: normal stresses σ xx , σyy, σzz, and shear stresses σxy, σyz, σzx. Store the stress tensor of each element to form the initial stress field data. These data record the initial stress distribution caused by temperature and humidity changes for each element. For example, the initial stress tensor of a certain element is:

[0096] σ = |0.5 MPa 0.1 MPa 0.05 MPa|

[0097] |0.1 MPa 0.6 MPa 0.08 MPa|

[0098] |0.05 MPa 0.08 MPa 0.4 MPa|

[0099] This means that the element is subject to a tensile stress of 0.5 MPa in the x - direction, a tensile stress of 0.6 MPa in the y - direction, a tensile stress of 0.4 MPa in the z - direction, and there are also small shear stresses. These initial stresses will be used as the initial conditions for subsequent analyses (such as transportation vibration analysis).

[0100] Preferably, the transportation vibration and strain analysis in step S1 includes:

[0101] Determine the vibration load time history curve of the transportation process, apply the vibration load time history curve to the mesh model data, and perform transient dynamic analysis based on the initial stress field data to obtain the node displacement time history data;

[0102] Define the contact pairs between the stacks in the mesh model data, and perform creep calculation on the contact pairs according to the material property data to obtain the contact force and creep strain data;

[0103] Perform stress-strain tensor field calculation based on the node displacement time history data, material property data, and contact force and creep strain data to obtain the stress-strain tensor field.

[0104] In the embodiment of the present invention, by installing a triaxial acceleration sensor on the chassis of an actual transport vehicle, the vibration acceleration data of the vehicle during driving on typical transport sections (including highways, urban roads, and rural roads) is recorded. The sampling frequency is set to 100 Hz, and the recording duration is 3 hours. Analyze the collected acceleration data, extract the acceleration signal in the vertical direction (Z-axis), and perform filtering processing (using a low-pass filter with a cut-off frequency of 20 Hz) to remove high-frequency noise. Use the filtered acceleration signal as the vibration load time history curve. Since the stack is mainly subjected to vertical vibration excitation during transportation, the acceleration time history curve is applied to the bottom nodes of the stack finite element mesh model. The specific method is as follows: Release the Z-direction displacement constraints of all bottom nodes and apply the displacement constraints corresponding to the acceleration time history curve. The value of the displacement constraint is obtained by integrating the acceleration signal twice with respect to time. Perform transient dynamic analysis. Use the Newmark-β method for time integration, and set the time step to 0.001 seconds. The total analysis duration is 10 seconds (select a relatively intense segment of the vibration acceleration signal). During the analysis, consider the initial stress field data calculated previously. This means that when calculating the node displacement at each time step, the initial stress field is added to the dynamic equation as an initial condition. The form of the dynamic equation is: M×ü + C×˙u + K×u = F(t) + F0, where M is the mass matrix, C is the damping matrix (using Rayleigh damping, and the damping coefficient is determined according to the subsequent modal analysis results), K is the stiffness matrix, u is the node displacement vector, ü and ˙u are the node acceleration vector and velocity vector respectively, F(t) is the external load vector (converted from the vibration acceleration time history curve), and F0 is the equivalent node force vector caused by the initial stress field. By solving the above dynamic equation, the node displacement vector u(t) at each time step is obtained. Combine the node displacement data at all time steps to form the node displacement time history data. This data records the displacement change of each node at each moment of the stack under vibration excitation.

[0105] In the finite element mesh model of the pallet stack, contact pairs are defined between the packaging boxes of adjacent goods. The face-to-face contact elements are used to simulate the contact behavior. The penalty function method is adopted for the contact algorithm, and the normal contact stiffness is set to 1×10 7 N / m, and the tangential contact stiffness is set to 5×10 6 N / m. The friction coefficient is set to 0.3. Creep calculations are performed on the contact pairs. The packaging box material is corrugated cardboard, and its creep behavior is described by the Burgers model. The Burgers model is composed of a Maxwell model and a Kelvin model in series. The Maxwell model is composed of a spring (elastic modulus E1) and a damper (viscosity coefficient η1) in series, and the Kelvin model is composed of a spring (elastic modulus E2) and a damper (viscosity coefficient η2) in parallel. Through creep tests, the parameters of the Burgers model of corrugated cardboard are measured as: E1 = 2 GPa, η1 = 50 GPa·s, E2 = 1 GPa, η2 = 10 GPa·s. At each time step, the creep strain increment on the contact surface is calculated according to the Burgers model. The calculation formula for the creep strain increment is: Δε_c = (σ / E1)×Δt + (σ / η1)×Δt + (σ / E2)×(1 - exp(-E2Δt / η2)), where σ is the normal stress on the contact surface and Δt is the time step. During the calculation process, the normal stress on the contact surface needs to be iteratively updated because creep will cause stress redistribution. The creep strain increments calculated at each time step are accumulated to obtain the total creep strain on the contact surface. At the same time, the normal pressure and tangential force (frictional force) on the contact surface at each time step are recorded to form the contact force and creep strain data. These data reflect the changes in the mechanical behavior of the contact surface under the combined action of vibration and creep of the pallet stack.

[0106] Based on the finite element theory, according to the nodal displacement time history data, material property data (elastic modulus, Poisson's ratio) of each element, and the contact force and creep strain data, the strain tensor and stress tensor of each element at each time step are calculated. Calculation of the strain tensor: First, according to the nodal displacement time history data, the shape function matrix B of each element is calculated. Then, according to the geometric equation ε = B×u, the total strain ε of each element at each time step is calculated. The total strain ε includes elastic strain, thermo-hygroscopic strain (from the calculation of the thermo-hygroscopic stress field in step S1), and creep strain (from the transportation vibration and strain analysis in step S1). Calculation of the stress tensor: According to the constitutive relationship (linear elastic model) of the material, the stress tensor σ of each element at each time step is calculated. The stress-strain relationship is: σ = D×(ε - ε th+hy - ε_c), where D is the elastic matrix, ε is the total strain, ε th+hy is the thermo-hygroscopic strain, and ε_c is the creep strain.

[0107] Combine the stress tensor and strain tensor of each time step and each element to form stress-strain tensor field data. This data completely records the changes in the stress and strain states of each point (element) inside the pallet over time during the entire transportation vibration process. For example, at t = 5 seconds, the stress tensor of a certain element is:

[0108] σ = |1.2MPa 0.3MPa 0.1MPa|

[0109] |0.3MPa 1.5MPa 0.2MPa|

[0110] |0.1MPa 0.2MPa 0.9MPa|

[0111] The strain tensor is:

[0112] ε = |0.0005 0.0001 0.00005|

[0113] |0.0001 0.0006 0.00008|

[0114] |0.00005 0.00008 0.0004|

[0115] This represents the stress and strain states of the element at this moment.

[0116] Preferably, the frequency vibration mode analysis in step S2 includes:

[0117] Extract the stiffness matrix and mass matrix from the mesh model data and material property data to obtain the stiffness mass matrix;

[0118] Construct the natural frequency characteristic equation based on the stiffness mass matrix to obtain the natural frequency characteristic equation;

[0119] Calculate the natural frequency and vibration mode according to the natural frequency characteristic equation to obtain the frequency vibration mode data.

[0120] Based on the finite element mesh model data of the pallet stack (including node coordinates, element types, element connection relationships, etc.) and material property data (elastic modulus, Poisson's ratio, density), calculate the stiffness matrix and mass matrix of each element. Taking the commonly used eight-node hexahedron element as an example, the calculation formula for the element stiffness matrix is: Ke = ∫∫∫BT × D × B × det(J)dξdηdζ, where B is the strain-displacement matrix, D is the elastic matrix, J is the Jacobian matrix, and ξ, η, ζ are local coordinates. The integration is calculated using the Gaussian quadrature method, and the number of integration points is 2×2×2. The calculation of the element mass matrix uses the lumped mass matrix method. For each element, its total mass is evenly distributed to 8 nodes. The calculation formula for the element mass is: Me = ρ × V, where ρ is the material density and V is the element volume. Assemble the stiffness matrices of all elements into the global stiffness matrix K. During the assembly process, according to the correspondence between the element node numbers and the global node numbers, each element in the element stiffness matrix is accumulated to the corresponding position in the global stiffness matrix. The global stiffness matrix K is a large sparse symmetric matrix. Assemble the mass matrices of all elements into the global mass matrix M. Since the lumped mass matrix method is used, the global mass matrix M is a diagonal matrix, and the elements on the diagonal are the masses of each node.

[0121] Based on the obtained global stiffness matrix K and global mass matrix M, construct the free vibration characteristic equation of the pallet stack. The form of the characteristic equation is: |K - ω 2 M| = 0, where ω is the natural frequency (circular frequency). This equation indicates that when the pallet stack vibrates at the natural frequency ω, the inertial force (related to ω 2 M) and the elastic restoring force (related to K) reach equilibrium. The characteristic equation can also be written in a more common form: (K - ω 2 M) × φ = 0, where φ is the mode shape vector corresponding to the natural frequency ω. The mode shape vector φ describes the relative displacement mode of each node when the pallet stack vibrates at the natural frequency ω.

[0122] Use the subspace iteration method to solve the characteristic equation (K - ω 2 M) × φ = 0. The subspace iteration method is a numerical iteration method suitable for solving partial eigenvalues and eigenvectors of large sparse matrices. Select an initial iteration vector group Φ0, and the number of its columns is greater than the number of eigenvalue-eigenvector pairs to be solved (for example, if the first 10 modes need to be solved, then select 20 initial iteration vectors).

[0123] Preferably, the calculation of the energy dissipation damping ratio in step S2 includes:

[0124] Select a typical time period for the stress-strain tensor field to obtain the modal time series;

[0125] Calculate the strain energy time history of the stress-strain tensor field according to the modal time series to obtain the strain energy curve;

[0126] Calculate the material damping dissipation energy based on the material property data and the stress-strain tensor field to obtain the material dissipation energy;

[0127] Calculate the contact friction dissipation energy according to the modal time series and the contact force and creep strain data to obtain the friction dissipation energy;

[0128] Calculate the creep dissipation energy based on the stress-strain tensor field, the modal time series, and the contact force and creep strain data to obtain the creep dissipation energy;

[0129] Conduct a hysteresis loop area analysis based on the modal time series, the material dissipation energy, the friction dissipation energy, and the creep dissipation energy to obtain the total dissipation energy;

[0130] Conduct modal energy distribution and modal dissipation calculation based on the total dissipation energy and the frequency mode shape data to obtain the modal dissipation distribution data;

[0131] Conduct damping ratio calculation based on the modal dissipation distribution data, the strain energy curve, and the frequency mode shape data under the principle of energy balance to obtain the damping ratio data.

[0132] In the embodiment of the present invention, from the stress-strain tensor field data obtained in step S1, select the time period during which the pallet shows obvious vibration attenuation characteristics during transportation vibration. The specific method is as follows: First, for each natural frequency ωi obtained from the frequency mode shape analysis in step S2, calculate its corresponding vibration period Ti = 2π / ωi. Then, find the interval in the stress-strain tensor field data where the node displacement amplitude gradually decreases with time, and this interval contains at least 5 complete vibration periods (based on the period T1 corresponding to the first-order natural frequency). For example, if the first-order natural frequency ω1 = 15 Hz, then the corresponding vibration period T1 ≈ 0.42 seconds. Select the time period from t = 2 seconds to t = 4 seconds (including approximately 4.8 vibration periods) in the stress-strain tensor field data as the typical time period of the first-order mode. Similarly, select the corresponding typical time periods for other modes. Extract the stress-strain data within the selected typical time period of each mode to form a modal time series. This series records the changes in stress and strain with time during the vibration attenuation process of the pallet in each mode.

[0133] For the modal time series of each mode, calculate the strain energy at each time step. The formula for calculating the strain energy is: U(t) = (1 / 2) × ∫Vσ(t):ε(t)dV, where σ(t) is the stress tensor, ε(t) is the strain tensor, and V is the volume of the pallet. The integration is performed over the entire volume of the pallet. Specifically, for each element of the pallet, calculate its strain energy: Ue(t) = (1 / 2) × σe(t):εe(t) × Ve, where σe(t) and εe(t) are the stress tensor and strain tensor of the element at time t (extracted from the modal time series), respectively, and Ve is the volume of the element. Then, sum up the strain energies of all elements to obtain the total strain energy U(t) of the pallet at this moment. For each time step in the modal time series, perform the above calculations to obtain the curve of the strain energy varying with time U(t). This curve reflects the fluctuations and dissipation of the strain energy during the vibration decay process of the pallet.

[0134] According to the material property data, the damping characteristics of corrugated board are described by a viscoelastic model. The constitutive relationship of the viscoelastic model is: σ(t) = E×ε(t) + η×˙ε(t), where E is the elastic modulus, η is the viscosity coefficient, and ˙ε(t) is the strain rate. Through dynamic mechanical analysis (DMA) tests, the storage modulus E' = 2.5 GPa and the loss modulus E” = 0.2 GPa of the corrugated board are measured. The viscosity coefficient η can be estimated by η = E” / ω, where ω is the angular frequency of vibration. Since the damping dissipation is mainly related to the lower-order modes, here the first natural frequency ω1 = 15 Hz is taken to calculate η, and η ≈ 2.12 MPa·s is obtained. For each time step in the modal time series, calculate the material damping dissipation power of each element: Pme(t) = σve(t):˙εve(t), where σve(t) is the viscous stress tensor and ˙εve(t) is the viscous strain rate tensor. The viscous stress tensor σve(t) = η×˙ε(t), and the viscous strain rate tensor ˙εve(t) is obtained by performing viscoelastic decomposition on the total strain rate ˙ε(t). The total strain rate ˙ε(t) is obtained by performing numerical differentiation on the strain tensor ε(t). Integrate the material damping dissipation power of each element over time to obtain the material damping dissipation energy of the element during this time period: where t1 and t2 are the start time and end time of the modal time series. Sum up the material damping dissipation energies of all elements to obtain the overall material damping dissipation energy Wm of the pallet during this time period.

[0135] According to the contact force and creep strain data obtained in step S1, extract the normal force Fn(t) and tangential force Ft(t) on the contact surface within the corresponding time period of the modal time series. At the same time, calculate the relative slip velocity vslip(t) on the contact surface. The relative slip velocity vslip(t) is obtained by numerically differentiating the relative displacement of the contact node pair. Calculate the frictional dissipation power on each contact surface: Pfe(t) = Ft(t)·vslip(t). Integrate the frictional dissipation power on each contact surface over time to obtain the frictional dissipation energy of this contact surface within this time period: Accumulate the frictional dissipation energies of all contact surfaces to obtain the overall contact frictional dissipation energy Wf of the pallet within this time period.

[0136] According to the stress-strain tensor field and the contact force and creep strain data obtained in step S1, extract the stress tensor σ(t) and creep strain εc(t) of each element within the corresponding time period of the modal time series. Calculate the creep dissipation power of each element: Pcre(t) = σ(t):˙εc(t), where ˙εc(t) is the creep strain rate tensor, which is obtained by numerically differentiating the creep strain εc(t). Integrate the creep dissipation power of each element over time to obtain the creep dissipation energy of this element within this time period: Accumulate the creep dissipation energies of all elements to obtain the overall creep dissipation energy Wcr of the pallet within this time period.

[0137] For the modal time series of each mode, plot the stress-strain hysteresis loop curve. The specific method is as follows: Select several elements in the pallet with significant stress and strain. For each element, use the normal stress in a certain direction (such as σ xx ) as the abscissa and the normal strain in this direction (such as ε xx ) as the ordinate to plot the stress-strain relationship curve within one vibration period. This curve usually presents as a closed loop, called the hysteresis loop. Calculate the area enclosed by each hysteresis loop. The size of the hysteresis loop area reflects the energy dissipation of this element within one vibration period. Use a numerical integration method (such as the trapezoidal rule) to calculate the hysteresis loop area. Average the hysteresis loop areas of all selected elements to obtain the average hysteresis loop area Wavg. This average hysteresis loop area approximately represents the total energy dissipation of the pallet in this mode. To verify the accuracy of the calculation results, compare the average hysteresis loop area Wavg with the previously calculated material damping dissipation energy Wm, contact frictional dissipation energy Wf, and creep dissipation energy Wcr. Theoretically, the total energy dissipation should be equal to the sum of these three: Wavg≈Wm + Wf + Wcr. If the deviation is large, it is necessary to check whether there are errors in the calculation process or whether there are other energy dissipation mechanisms not considered.

[0138] Based on the modal superposition principle, the total energy dissipation \(W_{total}\) (i.e., the average hysteresis loop area \(W_{avg}\) obtained in the previous step) is distributed to each mode. Assume that the vibration of the pallet can be expressed as a linear superposition of the vibrations of each order mode: \(u(t)=\sum_{i}q_{i}(t)\times\varphi_{i}\), where \(u(t)\) is the nodal displacement vector, \(q_{i}(t)\) is the modal coordinate of the \(i\)-th order mode, and \(\varphi_{i}\) is the normalized mode shape vector of the \(i\)-th order mode. Calculate the participation factor \(P_{i}\) of each order mode. The participation factor \(P_{i}\) reflects the weight of this mode in the total vibration. The calculation formula for the participation factor is: \(P_{i}=(\varphi_{i}\ T \times M\times r)\ 2 / (\varphi_{i}\ T \times M\times\varphi_{i})\), where \(\varphi_{i}\) is the normalized mode shape vector of the \(i\)-th order mode, \(M\) is the global mass matrix, and \(r\) is the excitation vector. Since the free vibration decay process is considered here, the excitation vector \(r\) can be taken as the unit vector. Distribute the total energy dissipation \(W_{total}\) according to the modal participation factor. The calculation formula for the energy dissipation \(\Delta W_{i}\) of the \(i\)-th order mode is: \(\Delta W_{i}=P_{i}\times W_{total}\). Combine all the calculated \(\Delta W_{i}\) to obtain the modal dissipation distribution data.

[0139] Based on the energy balance principle of the vibration decay of a single-degree-of-freedom system, derive the calculation formula for the damping ratio. For a single-degree-of-freedom system, the relationship between the energy dissipation \(\Delta W\) and the damping ratio \(\zeta\) in one vibration cycle is: \(\Delta W = 2\pi\zeta\omega U\), where \(\omega\) is the natural frequency and \(U\) is the maximum strain energy. Generalize the above relationship to a multi-degree-of-freedom system. For the \(i\)-th order mode, the calculation formula for its damping ratio \(\zeta_{i}\) is: \(\zeta_{i}=\Delta W_{i} / (4\pi U_{i})\), where \(\Delta W_{i}\) is the energy dissipation of the \(i\)-th order mode (obtained from the modal dissipation distribution data in the previous step), and \(U_{i}\) is the maximum strain energy of the \(i\)-th order mode. Calculation of the maximum strain energy \(U_{i}\): From the strain energy curve \(U(t)\) obtained in step 2, select the maximum value of the strain energy within the typical time period corresponding to the \(i\)-th order mode as \(U_{i}\). For each order mode, calculate its damping ratio \(\zeta_{i}\) according to the above formula. Combine all the calculated damping ratios \(\zeta_{i}\) to form the damping ratio data. For example, the damping ratios of the first three order modes are: \(\zeta_{1}=0.02\), \(\zeta_{2}=0.03\), \(\zeta_{3}=0.04\). These data reflect the energy dissipation ability of the pallet in different vibration modes.

[0140] Preferably, the instability critical domain analysis in step S3 includes:

[0141] Construct a dynamic stability criterion based on the modal parameter set to obtain a dynamic stability factor;

[0142] Obtain the geometric information of the pallet; establish a tipping instability model based on the geometric information of the pallet to obtain a tipping instability model;

[0143] Calculate the critical inclination angle of the tipping instability model based on the contact force and creep strain data to obtain the critical inclination angle;

[0144] Establish a slip instability model based on the contact force, creep strain data and material property data to obtain the slip instability model;

[0145] Calculate the critical slip distance of the slip instability model to obtain the critical slip distance.

[0146] In the embodiment of the present invention, based on the modal parameter set (including natural frequency and damping ratio) obtained in step S2, a dynamic stability criterion is constructed. Two main factors are considered: natural frequency and damping ratio. Natural frequency criterion: If the lowest natural frequency ω1 of the stack is lower than a certain threshold ωmin, it is considered that the stack is prone to resonance and the dynamic stability is insufficient. The setting of the threshold ωmin is related to the vibration frequency range during transportation. For example, if the vibration frequency of the transport vehicle mainly concentrates in the range of 5 - 10 Hz, then ωmin = 5 Hz can be set. Damping ratio criterion: If the lowest modal damping ratio ζ1 of the stack is lower than a certain threshold ζmin, it is considered that the vibration attenuation ability of the stack is weak and the dynamic stability is insufficient. The setting of the threshold ζmin is related to the damping characteristics of the stack material and actual engineering experience. For example, ζmin = 0.01 can be set. Combining the above two criteria, a dynamic stability factor D is formed. The calculation formula of D is: D = w1×(ω1 / ωmin)+w2×(ζ1 / ζmin), where w1 and w2 are weight coefficients and w1 + w2 = 1. The values of the weight coefficients reflect the relative importance of the natural frequency and damping ratio to the dynamic stability. For example, w1 = 0.6 and w2 = 0.4 can be taken. If D≥1, it is considered that the dynamic stability of the stack meets the requirements; if D<1, it is considered that the dynamic stability of the stack is insufficient.

[0147] Obtain the geometric information of the pallet stack, including the overall dimensions of the pallet stack (length L, width W, height H), the dimensions of each cargo box, the stacking pattern of the cargo boxes (such as row pattern, staggered pattern, etc.), and the centroid position of the pallet stack. This information can be obtained from the warehouse management system or through on-site measurement. Establish a mechanical model of tipping instability. Simplify the pallet stack into a rigid body and consider the case of tipping around the bottom edge. Establish a moment balance equation: M_overturning = M_restoring, where M_overturning is the overturning moment and M_restoring is the anti-overturning moment. Calculation of the overturning moment: M_overturning = m × g × h × sinθ, where m is the total mass of the pallet stack, g is the acceleration due to gravity, h is the height of the centroid of the pallet stack, and θ is the inclination angle of the pallet stack. Calculation of the anti-overturning moment: M_restoring = m × g × (W / 2) × cosθ, where W is the width of the pallet stack (assuming tipping occurs in the width direction).

[0148] Based on the contact force and creep strain data obtained in step S1, determine the friction coefficient μ and the normal pressure distribution N(x, y) on the bottom contact surface of the pallet stack. Since creep causes stress redistribution on the contact surface, the friction coefficient and the normal pressure distribution vary with time. Gradually increase the inclination angle θ of the pallet stack and calculate the overturning moment M_overturning and the anti-overturning moment M_restoring at each inclination angle. When calculating the anti-overturning moment, it is necessary to consider the normal pressure distribution N(x, y) on the contact surface and the change in normal pressure due to inclination. When the overturning moment M_overturning is greater than the anti-overturning moment M_restoring, the corresponding inclination angle is the critical inclination angle θc. Since the friction coefficient and the normal pressure distribution vary with time, the critical inclination angle θc also varies with time.

[0149] Based on the contact force and creep strain data obtained in step S1, determine the normal force Fn and the tangential force Ft (frictional force) on the bottom contact surface of the pallet stack. At the same time, based on the material property data, determine the static friction coefficient μs of the bottom contact surface of the pallet stack. Establish a mechanical model of sliding instability. Consider the case of the pallet stack sliding in the horizontal direction. Establish a force balance equation: F_driving = F_resisting, where F_driving is the driving force (such as the horizontal inertial force during transportation) and F_resisting is the anti-sliding force. Calculation of the driving force: F_driving = m × a, where m is the total mass of the pallet stack and a is the acceleration in the horizontal direction (which can be obtained from the transportation vibration data). Calculation of the anti-sliding force: F_resisting = μs × Fn, where μs is the static friction coefficient and Fn is the normal force.

[0150] Gradually increase the driving force \(F_{driving}\) in the horizontal direction and calculate the anti-slip force \(F_{resisting}\) under each driving force. When calculating the anti-slip force, the normal force distribution and friction coefficient distribution on the contact surface need to be considered. When the driving force \(F_{driving}\) is greater than the anti-slip force \(F_{resisting}\), the stack begins to slip. At this time, the corresponding horizontal displacement is the critical slip distance. Since the normal force and friction coefficient on the contact surface change with time (affected by creep), the critical slip distance also changes with time. Define the critical slip distance as the relative displacement of the bottom of the stack with respect to the ground when the stack begins to slip. Due to the non-uniformity of the contact surface and the creep effect, the critical slip distance is not a fixed value but a quantity that changes with time. By gradually increasing the horizontal displacement applied to the bottom of the stack and monitoring the horizontal reaction force, when the horizontal reaction force reaches the maximum value and then begins to decrease, the corresponding horizontal displacement is the critical slip distance.

[0151] Preferably, the long-term stability prediction of the creep influence factor in step S3 includes:

[0152] Obtain the data of the instability critical domain, where the data of the instability critical domain includes the dynamic stability factor, the critical inclination angle, and the critical slip distance; establish a creep constitutive model and calibrate the parameters according to the material property data to obtain a creep model;

[0153] Determine the time discretization strategy according to the creep model to obtain a time discretization scheme;

[0154] Calculate the creep strain increment according to the time discretization scheme, the creep model, and the contact force and creep strain data to obtain creep evolution data;

[0155] Update the geometric shape and contact characteristics of the stack according to the creep evolution data to obtain the deformation state data;

[0156] Conduct a rebalancing analysis according to the deformation state data and perform a stress redistribution calculation to obtain a rebalanced stress field;

[0157] Recalculate the dynamic stability factor according to the rebalanced stress field, the deformation state data, and the dynamic stability factor to obtain the dynamic stability time history;

[0158] Recalculate the critical inclination angle according to the deformation state data, the rebalanced stress field, and the critical inclination angle to obtain the inclination angle time history;

[0159] Recalculate the critical slip distance according to the deformation state data, the rebalanced stress field, and the critical slip distance to obtain the slip threshold time history;

[0160] Integrate the time-varying stability parameters for the dynamic stability history, inclination angle history, and slip threshold history to obtain the time-varying stability parameters; calculate the stability margin factor based on the time-varying stability parameters to obtain the stability margin factor.

[0161] In the embodiments of the present invention, obtain the instability critical domain data from the previous steps of step S3, including: the dynamic stability factor D (initial value), the critical inclination angle θc (initial value), and the critical slip distance (initial value). These data are the starting conditions for long-term stability prediction. According to the material property data, select a suitable creep constitutive model for the packaging materials of the goods (such as corrugated cardboard, plastic film, etc.). For corrugated cardboard, the Burgers model is used to describe its creep behavior. The Burgers model is composed of a Maxwell model and a Kelvin model in series. The Maxwell model is composed of a spring (elastic modulus E1) and a damper (viscosity coefficient η1) in series, and the Kelvin model is composed of a spring (elastic modulus E2) and a damper (viscosity coefficient η2) in parallel. Through creep tests, determine the Burgers model parameters of corrugated cardboard. The tests are carried out under constant temperature and humidity conditions (temperature 25°C, relative humidity 60% RH). Apply a constant tensile stress to the corrugated cardboard specimen and record the curve of its strain changing with time. By fitting the test data, obtain the Burgers model parameters: E1 = 2 GPa, η1 = 50 GPa·s, E2 = 1 GPa, η2 = 10 GPa·s. For plastic film, the power-law creep model is used to describe its creep behavior. The expression of the power-law creep model is: εc = A×σ n ×t m , where εc is the creep strain, σ is the stress, t is the time, and A, n, m are material parameters. Through creep tests, determine the power-law creep model parameters of plastic film: A = 1×10 -10 MPa -n ·s -m , n = 1.5, m = 0.3.

[0162] According to the characteristics of the selected creep models (Burgers model and power-law model), determine the time discretization strategy. Since the creep process changes rapidly in the initial stage and then gradually slows down, a non-uniform time step is adopted. The initial time step is set to 1 hour. The subsequent time steps gradually increase, and the growth coefficient is 1.1. That is, Δt i+1 = 1.1×Δt i . The total prediction time is set to 30 days (720 hours). According to the above time step setting, calculate the total number of time steps N. In each time step, use the implicit time integration method to calculate the creep strain increment. The implicit method has good numerical stability and is suitable for long-term creep analysis.

[0163] Based on the contact force and creep strain data obtained in step S1, and the time discretization scheme determined in the previous step, calculate the creep strain increment of each element within each time step. For corrugated cardboard elements, the Burgers model is used to calculate the creep strain increment. Within each time step, according to the current stress state and Burgers model parameters, calculate the creep strain increment Δεc. The calculation formula is: Δεc = (σ / E1)×Δt + (σ / η1)×Δt + (σ / E2)×(1 - exp(-E2Δt / η2)), where σ is the current stress and Δt is the current time step. For plastic film elements, the power-law model is used to calculate the creep strain increment. The calculation formula is: Δεc = A×σ n ×[(t + Δt) m -t m , where σ is the current stress, t is the current time, and Δt is the current time step. At the end of each time step, accumulate the calculated creep strain increment to the total creep strain of the element. That is, εc(t + Δt) = εc(t) + Δεc. Record the creep strain of each element at each time step to form creep evolution data. This data reflects the change in deformation over time caused by the creep effect during the long-term stacking of the pallet stack.

[0164] Update the geometry of the pallet stack according to the creep evolution data calculated in the previous step. The specific method is: for each node of the pallet stack, calculate the displacement increment of the node according to the creep strain of its adjacent elements. Then, add the initial coordinates of the node to the displacement increment to obtain the new coordinates of the node. Update the contact characteristics of the pallet stack. Since creep will cause stress redistribution on the contact surface, it is necessary to recalculate the normal force and tangential force on the contact surface. The penalty function method is used to calculate the contact force. The normal contact stiffness remains unchanged, and the tangential contact stiffness is adjusted according to the creep strain. Record the updated node coordinates, normal force and tangential force on the contact surface to form deformation state data. This data reflects the changes in the geometry and contact state of the pallet stack under the action of creep.

[0165] Based on the updated pallet geometry and contact characteristics, a static rebalancing analysis is carried out. Since creep will cause stress redistribution, it is necessary to recalculate the stress state inside the pallet. The static equilibrium equation is established: K×u = F, where K is the global stiffness matrix, u is the nodal displacement vector, and F is the nodal force vector. The nodal force vector F includes gravity, contact forces, and equivalent nodal forces caused by creep. The Newton - Raphson iteration method is used to solve the above - mentioned nonlinear equation. In each iteration, according to the current nodal displacements, the global stiffness matrix K and the nodal force vector F are updated. After the iteration converges, the new nodal displacement vector u is obtained. According to the new nodal displacement vector u, the strain of each element is calculated. Then, according to the constitutive relationship of the material (linear elastic model), the stress of each element is calculated. The calculated stress is used as the rebalanced stress field.

[0166] Based on the updated pallet geometry (from the deformed state data) and the rebalanced stress field, the modal analysis is carried out again. According to the frequency - mode analysis method described in step S2, the natural frequencies and damping ratios of the pallet in the current state are calculated. According to the new natural frequencies and damping ratios, the dynamic stability factor D is recalculated. The calculation formula is the same as the dynamic stability factor calculation formula in step S3: D = w1×(ω1 / ωmin)+w2×(ζ1 / ζmin). The dynamic stability factor D at each time step is recorded to form the dynamic stability time history D(t). This time history reflects the change trend of the dynamic stability of the pallet with time under the action of creep.

[0167] Based on the updated pallet geometry (from the deformed state data) and the rebalanced stress field, the critical inclination angle of the pallet is recalculated. According to the critical inclination angle calculation method described in step S3, the inclination angle θ of the pallet is gradually increased, and the overturning moment M_overturning and the restoring moment M_restoring at each inclination angle are calculated. During the calculation process, the updated position of the pallet's center of gravity, the normal pressure distribution on the contact surface, and the friction coefficient need to be considered. When the overturning moment M_overturning is greater than the restoring moment M_restoring, the corresponding inclination angle is the critical inclination angle θc(t) at the current moment. The critical inclination angle θc(t) at each time step is recorded to form the inclination angle time history. This time history reflects the change trend of the anti - overturning ability of the pallet with time under the action of creep.

[0168] Recalculate the critical slip distance of the stack based on the updated stack geometry (from the deformed state data) and the rebalanced stress field. According to the critical slip distance calculation method described in step S3, gradually increase the driving force F_driving in the horizontal direction and calculate the anti-slip force F_resisting at each driving force. During the calculation, the normal force distribution and friction coefficient on the updated contact surface need to be considered. When the driving force F_driving is greater than the anti-slip force F_resisting, the stack begins to slip. At this time, the corresponding horizontal displacement is the critical slip distance (t) at the current moment. Record the critical slip distance (t) at each time step to form a slip threshold time history. This time history reflects the changing trend of the stack's anti-slip ability over time under creep action.

[0169] Integrate the previously calculated dynamic stability time history D(t), inclination angle time history θc(t), and slip threshold time history (t) to form time-varying stability parameters. These three parameters reflect the stability changes of the stack under creep action from different aspects. To comprehensively evaluate the overall stability of the stack, define a stability margin factor S(t). The calculation formula for the stability margin factor S(t) is: S(t) = min{D(t) / D0, θc(t) / θc0, (t) / 0}, where D0, θc0, and 0 are the initial values (or set safety thresholds) of the dynamic stability factor, critical inclination angle, and critical slip distance, respectively. The stability margin factor S(t) represents the stability margin of the stack at the current moment relative to the initial state (or safe state). If S(t) ≥ 1, it is considered that the stack is stable at the current moment; if S(t) < 1, it is considered that the stack has a risk of instability at the current moment. Calculate the stability margin factor S(t) at each time step and use it as the final output result. This result reflects the changing trend of the overall stability of the stack over time under long-term creep action.

[0170] As an example of the present invention, refer to Figure 2 As shown, in this example, step S4 includes:

[0171] Step S41: Divide the risk levels according to the stability margin factor to obtain a risk level threshold table;

[0172] Step S42: Identify risks according to the risk level threshold table and the stability margin factor to obtain a risk distribution map;

[0173] Step S43: Formulate storage strategy adjustment rules according to the risk level threshold table to obtain storage adjustment rules;

[0174] Step S44: Obtain real-time warehouse temperature and humidity data; perform dynamic adjustment and optimization iteration according to the real-time warehouse temperature and humidity data, risk distribution map, stability margin factor, and storage adjustment rules to obtain an optimized risk map;

[0175] Step S45: Generate a warehousing optimization plan based on the optimized risk graph to obtain a warehousing optimization plan.

[0176] In the embodiment of the present invention, according to the numerical range of the stability margin factor S(t) obtained in step S3 and the actual warehouse management experience, the stability risk of the cargo stack is divided into four levels: safe, low risk, medium risk, and high risk. The threshold interval of the stability margin factor is set for each risk level. Specifically as follows:

[0177] Safety: S(t)≥1.5. This means that the stability of the cargo stack is much higher than the safety requirement and no additional measures are required.

[0178] Low risk: 1.2≤S(t)<1.5. This means that the stability of the cargo stack is slightly higher than the safety requirement, but there is a certain risk and regular monitoring is required.

[0179] Medium risk: 1.0≤S(t)<1.2. This means that the stability of the cargo stack is close to the critical state of safety requirements and certain adjustment measures need to be taken.

[0180] High risk: S(t)<1.0. This means that the stability of the cargo stack is lower than the safety requirement and there is a high risk of instability. Immediate measures must be taken.

[0181] The above risk levels and their corresponding stability margin factor threshold intervals are organized into a table to form a risk level threshold table, which will serve as the basis for subsequent risk assessment and warehousing strategy adjustment.

[0182] According to the existing storage scheme (including the stacking method, height, location, etc. of goods), combined with the model calculation of steps S1-S3, the stability margin factor S(t) of each cargo stack is obtained. The stability margin factor S(t) of each cargo stack is compared with the risk level threshold table to determine the risk level to which it belongs. For example, if S(t) = 1.3 of a certain cargo stack, then according to the risk level threshold table, the cargo stack belongs to the "low risk" level. Generate a risk distribution map. The risk distribution map uses the warehouse plan as the base map and uses different colors to represent cargo stacks of different risk levels. For example, green represents the "safe" level, blue represents the "low risk" level, yellow represents the "medium risk" level, and red represents the "high risk" level. On the risk distribution map, you can clearly see the risk status of each cargo stack in the warehouse and the distribution of high-risk areas.

[0183] According to the risk level threshold table, formulate corresponding storage strategy adjustment rules for cargo stacks with different risk levels.

[0184] Security level: Keep the existing storage plan unchanged.

[0185] Low risk level:

[0186] 1. Strengthen temperature and humidity monitoring to ensure that the temperature and humidity fluctuations in the warehouse are within an acceptable range.

[0187] 2. Regularly check the stability of the pallet stack. If there is slight deformation, adjust the position of the goods in a timely manner.

[0188] 3. Avoid placing heavy objects above the pallet stack or performing operations that may cause vibration.

[0189] Medium risk level:

[0190] 4. Reduce the height of the pallet stack. For example, reduce the pallet stack that was originally stacked in 5 layers to 4 layers.

[0191] 5. Adjust the stacking method of the goods. For example, place heavier goods at the bottom of the pallet stack and lighter goods at the top.

[0192] 6. Increase the spacing between pallet stacks.

[0193] 7. If conditions permit, take reinforcement measures for the pallet stack.

[0194] High risk level:

[0195] 8. Immediately unload some goods to reduce the height and weight of the pallet stack.

[0196] 9. Reinforce the pallet stack, such as using straps, support frames, etc.

[0197] 10. Adjust the temperature and humidity in the warehouse to reduce the creep rate.

[0198] 11. Re-evaluate the stability of the pallet stack and replace the goods or adjust the stacking position if necessary.

[0199] Quantify the above adjustment rules into specific operation parameters, such as the maximum stacking height, the weight limit of the goods, the temperature and humidity control range, etc., and organize them into a document to form a warehouse storage adjustment rule library.

[0200] Through the temperature and humidity sensors installed in the warehouse, the temperature and humidity data in the warehouse are collected in real time. The data collection frequency is once per minute. According to the real-time temperature and humidity data, update the model calculations in steps S1 - S3 to obtain the latest stability margin factor S(t) for each pallet stack. According to the latest stability margin factor S(t), update the risk distribution map. If the risk level of a certain pallet stack changes (such as from medium risk to high risk), then adjust its color accordingly on the risk distribution map.

[0201] Analyze the new risk distribution map according to the warehousing adjustment rule base to determine whether there are stacks that need to be adjusted. For the stacks that need to be adjusted, adjust them according to their risk levels and the corresponding rules in the warehousing adjustment rule base. For example, if the risk level of a certain stack changes from medium risk to high risk, it will be processed according to the adjustment rules for high risk levels (such as immediately unloading part of the goods, strengthening, adjusting temperature and humidity, etc.). After the adjustment is completed, re-run the model in steps S1 - S3, calculate the stability margin factor S(t) of each stack after the adjustment, and evaluate the adjustment effect. If the stability margin factor after the adjustment still does not meet the requirements (such as still being at a high risk level), further adjust the warehousing strategy and repeat the above steps until the risk levels of all stacks are within the acceptable range. Output the final risk distribution map as the optimized risk map. This map reflects the latest risk status of the stacks in the warehouse after dynamic adjustment and optimization iteration.

[0202] Generate the final warehousing optimization plan based on the optimized risk map and in combination with the warehousing adjustment rule base. The warehousing optimization plan includes the following contents:

[0203] Stacking method of goods: Clearly stipulate the best stacking method for each type of goods (such as row stacking, staggered stacking, etc.), as well as the maximum height, width, and length of the stack.

[0204] Location of goods: According to the optimized risk map, specify the best location of each stack in the warehouse to avoid placing high-risk stacks near key passages or exits.

[0205] Temperature and humidity control requirements: Specify the control range of temperature and humidity in the warehouse, as well as the adjustment strategies under different seasons or weather conditions.

[0206] Reinforcement measures: For the stacks that need to be reinforced, clearly stipulate the reinforcement methods (such as using straps, support frames, etc.) and the specifications of the reinforcement materials.

[0207] Monitoring and inspection: Specify the frequency and methods of regularly monitoring temperature and humidity, inspecting the stability of the stacks, as well as the handling process after problems are found.

[0208] Risk distribution map: Provide a visual interface.

[0209] Organize the above content into a detailed document or data table and provide a visual interface (such as the optimized risk map) to form the final warehousing optimization plan. This plan provides clear operation guidelines for warehouse managers and helps improve the safety and management efficiency of the warehouse.

[0210] Therefore, in any aspect, the embodiments should be regarded as exemplary and non-limiting. The scope of the present invention is defined by the appended claims rather than the above description. Therefore, all changes falling within the meaning and scope of the equivalent elements of the application documents are intended to be encompassed within the present invention.

[0211] The above description is only a specific implementation manner of the present invention, enabling those skilled in the art to understand or implement the present invention. Various modifications to these embodiments will be obvious to those skilled in the art. The general principles defined herein can be implemented in other embodiments without departing from the spirit or scope of the present invention. Therefore, the present invention will not be limited to these embodiments shown herein, but rather to the broadest scope consistent with the principles and novel features invented herein.

Claims

1. A method for determining a warehouse management model, characterized in that: The following steps are involved: Step S1: constructing grid model data of cargo stacks; determining material property data of cargo; collecting warehouse temperature and humidity data; loading warehouse temperature and humidity data into grid model data, calculating temperature and humidity stress field according to material property data, and obtaining initial stress field data; According to the initial stress field data, the transport vibration and strain analysis of the cargo stack is carried out to obtain the contact force and creep strain data and the stress-strain tensor field; Step S2: Perform frequency mode analysis on the grid model data and material property data to obtain frequency mode data; Calculate the energy dissipation damping ratio according to the frequency vibration mode data and the stress-strain tensor field to obtain the damping ratio data; integrate the damping ratio data and the frequency vibration mode data into a modal parameter set; Step S3: performing an instability critical domain analysis on the cargo stack according to the modal parameter set to obtain instability critical domain data; performing a long-term stability prediction of the creep influence factor according to the instability critical domain data and the contact force and creep strain data to obtain a stability margin factor; Step S4: Set storage adjustment rules under different risk levels according to the stability margin factor, optimize the real-time state storage strategy according to the storage adjustment rules, and obtain a storage optimization plan.

2. The method for determining a warehouse management model according to claim 1, characterized in that: The temperature and humidity stress field calculation in step S1 includes: According to the warehouse temperature and humidity data, a continuous function of the temperature and humidity field is constructed to obtain the temperature and humidity field function; The boundary conditions of the temperature and humidity field function are set to obtain the temperature and humidity field; the temperature and humidity field is mapped to the grid model data, and the temperature and humidity gradient field is calculated to obtain the grid temperature and humidity data; The thermal expansion and hygroscopic expansion strains are calculated based on the grid temperature and humidity data and material property data to obtain the thermal and humidity strain field; The initial stress field is calculated based on the thermal and moist strain field to obtain the initial stress field data.

3. The method for determining a warehouse management model according to claim 1, characterized in that: The transport vibration and strain analysis in step S1 includes: Determine the vibration load time history curve of the transportation process, apply the vibration load time history curve to the grid model data, perform transient dynamic analysis based on the initial stress field data, and obtain the node displacement time history data; The contact pairs between the cargo stacks in the mesh model data are defined. The contact pairs use surface-to-surface contact elements to simulate the contact behavior. The surfaces between adjacent cargo packaging boxes are used as contact elements. The normal contact stiffness of the contact element is set to 1×10 7 N / m, and the tangential contact stiffness is set to 5×10 6 N / m, the friction coefficient is set to 0.3, and the creep calculation of the contact pair is performed based on the material property data to obtain the contact force and creep strain data; The stress-strain tensor field is calculated based on the node displacement time history data, material property data, contact force and creep strain data to obtain the stress-strain tensor field.

4. The method for determining a warehouse management model according to claim 1, characterized in that: The frequency mode analysis in step S2 includes: Extract the stiffness matrix and mass matrix from the mesh model data and material property data to obtain the stiffness mass matrix; The inherent rate characteristic equation is constructed according to the stiffness mass matrix to obtain the inherent rate characteristic equation; The natural frequency and vibration mode are calculated according to the natural rate characteristic equation to obtain the frequency vibration mode data.

5. The method for determining a warehouse management model according to claim 1, characterized in that: The energy consumption damping ratio calculation in step S2 includes: Select typical time periods for the stress-strain tensor field to obtain modal time series; The strain energy time history of the stress-strain tensor field is calculated according to the modal time series to obtain the strain energy curve; Calculate the material damping dissipation energy according to the material property data and the stress-strain tensor field to obtain the material dissipation energy; The contact friction dissipation energy is calculated based on the modal time series and the contact force and creep strain data to obtain the friction dissipation energy; The creep dissipated energy is calculated based on the stress-strain tensor field, modal time series, contact force and creep strain data to obtain the creep dissipated energy; The total dissipated energy is obtained by performing hysteresis loop area analysis based on the modal time series, material dissipated energy, friction dissipated energy and creep dissipated energy; Perform modal energy distribution and modal dissipation calculations based on total dissipated energy and frequency vibration mode data to obtain modal dissipation distribution data; The damping ratio is calculated based on the energy balance principle according to the modal dissipation distribution data, strain energy curve and frequency vibration shape data to obtain the damping ratio data.

6. The method for determining a warehouse management model according to claim 1, characterized in that: The instability critical region analysis in step S3 includes: The dynamic stability criterion is constructed according to the modal parameter set to obtain the dynamic stability factor; Acquire geometric information of the cargo stack; establish a dumping instability model based on the geometric information of the cargo stack to obtain a dumping instability model; The critical inclination angle of the toppling instability model is calculated based on the contact force and creep strain data to obtain the critical inclination angle; A sliding instability model is established based on contact force, creep strain data and material property data to obtain a sliding instability model; The critical slip distance of the slip instability model is calculated and the critical slip distance is obtained.

7. The method for determining a warehouse management model according to claim 1, characterized in that: The long-term stability prediction of the creep influencing factor in step S3 includes: Acquire the instability critical domain data, where the instability critical domain data includes dynamic stability factor, critical inclination angle and critical slip distance; establish the creep constitutive model and calibrate the parameters according to the material property data to obtain the creep model; Determine the time discretization strategy according to the creep model and obtain the time discretization scheme; The creep strain increment is calculated according to the time discretization scheme, creep model, contact force and creep strain data to obtain the creep evolution data; Update the geometric shape and contact characteristics of the cargo stack according to the creep evolution data to obtain the deformation state data; Rebalance analysis is performed based on deformation state data, and stress redistribution calculation is performed to obtain a rebalanced stress field; Recalculate the dynamic stability factor according to the rebalanced stress field, deformation state data and dynamic stability factor to obtain the dynamic stability history; According to the deformation state data, the rebalanced stress field and the critical inclination angle, the critical inclination angle is recalculated to obtain the inclination angle time history; According to the deformation state data, the rebalanced stress field and the critical slip distance, the critical slip distance is recalculated to obtain the slip threshold time history; The time-varying stability parameters of the dynamic stability time history, the tilt angle time history and the slip threshold time history are integrated to obtain the time-varying stability parameters; the stability margin factor is calculated according to the time-varying stability parameters to obtain the stability margin factor.

8. The method for determining a warehouse management model according to claim 1, characterized in that: Step S4 includes the following steps: Step S41: classify the risk levels according to the stability margin factors to obtain a risk level threshold table; Step S42: performing risk identification according to the risk level threshold table and the stability margin factor to obtain a risk distribution map; Step S43: Formulate storage strategy adjustment rules according to the risk level threshold table to obtain storage adjustment rules; Step S44: acquiring real-time warehouse temperature and humidity data; performing dynamic adjustment and optimization iteration according to the real-time warehouse temperature and humidity data, risk distribution map, stability margin factor and storage adjustment rules to obtain an optimized risk map; Step S45: Generate a warehousing optimization plan based on the optimized risk graph to obtain a warehousing optimization plan.

9. A warehouse management model determination system, characterized in that: Used to execute the warehouse management model determination method according to claim 1, the warehouse management model determination system comprises: The multi-field coupling characterization module is used to construct the grid model data of the cargo stack; determine the material property data of the cargo; collect the temperature and humidity data of the warehouse; load the temperature and humidity data of the warehouse into the grid model data, calculate the temperature and humidity stress field according to the material property data, and obtain the initial stress field data; perform transportation vibration and strain analysis on the cargo stack according to the initial stress field data, and obtain the contact force and creep strain data and stress-strain tensor field; The modal damping identification module is used to perform frequency vibration mode analysis on the mesh model data and material property data to obtain frequency vibration mode data; calculate the energy loss damping ratio according to the frequency vibration mode data and the stress-strain tensor field to obtain the damping ratio data; integrate the damping ratio data and the frequency vibration mode data into a modal parameter set; The stability margin assessment module is used to analyze the instability critical domain of the cargo stack according to the modal parameter set to obtain the instability critical domain data; the long-term stability prediction of the creep influence factor is performed according to the instability critical domain data and the contact force and creep strain data to obtain the stability margin factor; The warehousing strategy optimization module is used to set warehousing adjustment rules under different risk levels according to the stability margin factor, and optimize the real-time warehousing strategy according to the warehousing adjustment rules to obtain the warehousing optimization plan.

10. A computer-readable storage medium storing a computer program, characterized in that: When the computer program is executed, the warehouse management model determination method according to any one of claims 1 to 8 is implemented.

Images