Wooden packaging box product optimization design simulation method based on digital twinning

Abstract

This invention discloses a simulation method for optimizing the design of wooden packaging boxes based on digital twins, relating to the fields of wooden packaging box design and computer-aided engineering simulation technology. The method includes: collecting the geometric and physical properties of the product to be packaged and the initial design parameters of the wooden packaging box; constructing a multi-physics digital twin of the product and a digital twin of the packaging box structure; inputting historical transportation environmental load spectra into the product twin; obtaining the product's acceleration and stress spatiotemporal distribution through transient response simulation; and drawing a map of the product's vulnerable and sensitive areas. Mapping the map to the contact boundary of the packaging box twin; and conducting multi-condition simulations of drop and vibration coupling between the product and packaging to obtain the cumulative evolution sequence of damage to key components of the packaging box. Using this sequence as the constraint objective function, iteratively optimizing the initial parameters of the packaging box using a parametric optimization engine to obtain an optimized parameter combination. Updating the packaging box digital twin based on the new parameters; cyclically simulating and verifying until convergence requirements are met; and outputting the final wooden packaging box design scheme.

Description

Technical Field

[0001] This invention relates to the field of wooden packaging box design and computer-aided engineering simulation technology, specifically a simulation method for optimizing the design of wooden packaging box products based on digital twins. Background Technology

[0002] For product protection design of wooden packaging boxes, existing methods typically rely on engineers' experience for structural selection and dimensional estimation, followed by independent static or dynamic verification of the packaging box using finite element software. In this process, the packaging object is often oversimplified as a rigid mass block or only its outer geometry is considered, ignoring the differences in dynamic response exhibited by the product itself under transportation stress and the weaknesses in its internal structure. This simplification results in packaging design failing to accurately respond to the protection needs of specific parts of the product, resorting only to overall thickening of the board or adding redundant lining to meet safety margins, leading to material waste and increased manufacturing costs. Furthermore, existing design processes are mostly open-loop verification, meaning that after selecting design parameters, a simulation verification is performed, and parameters are manually adjusted based on the verification results for further verification. This lacks a closed-loop optimization mechanism that directly links product damage status with packaging structural parameters. Under this model, designers struggle to efficiently obtain the optimal solution that balances lightweighting and structural integrity from numerous parameter combinations. Each iteration requires manual modeling and setting of simulation conditions, resulting in long design cycles and difficulty in exhaustively exploring the optimization space.

[0003] The problems that this solution needs to address are: how to accurately identify the vulnerable parts of a product under transportation loads within a coupled simulation framework, and thereby drive the differentiated parameter configuration of the packaging structure; and how to quantify the material damage accumulation process obtained in the simulation into constraints that can drive the optimization engine to automatically seek optimization. Summary of the Invention

[0004] This paper proposes a simulation method for optimizing the design of wooden packaging boxes based on digital twins. The method realistically reproduces the dynamic stress and acceleration spatiotemporal distribution of the product under transportation loads in the digital twin environment, automatically identifies the product's vulnerable areas, and constructs a fully coupled simulation model of the product and packaging using the vulnerable area map as a mapping benchmark. Then, the damage accumulation evolution sequence of the key components of the wooden packaging box is extracted and transformed into a constraint-driven parametric optimization engine to automatically find the optimal solution in the variable space. Finally, the design scheme of the wooden packaging box is output after closed-loop convergence verification, realizing the synergy between the packaging structure's precise protection of specific parts of the product and lightweight design.

[0005] To achieve the above objectives, the present invention provides the following technical solution: The present invention provides a simulation method for optimized design of wooden packaging boxes based on digital twins, the method comprising the following:

[0006] The process involves acquiring the geometric model and physical properties of the product to be packaged, as well as the initial design parameters of the wooden crate, to construct a multiphysics digital twin of the product and a structural digital twin of the wooden crate. Historical transportation environmental load spectra are applied to the multiphysics digital twin of the product, and the spatiotemporal distribution of the product's acceleration and stress response is extracted through transient response simulation, generating a vulnerability sensitive area map. This vulnerability sensitive area map is mapped to the contact boundaries of the wooden crate structural digital twin, and a multi-condition simulation of drop and vibration coupled to the product is performed to obtain the damage accumulation evolution sequence of key components of the wooden crate. Using the damage accumulation evolution sequence as the objective function constraint, a parametric optimization engine is invoked to iteratively optimize the initial design parameters of the wooden crate, generating an optimized design parameter combination. The structural digital twin of the wooden crate is updated based on the optimized design parameter combination, and simulation verification is repeated until the convergence condition is met, outputting the final wooden crate design scheme.

[0007] As a preferred technical solution of the present invention, in the process of constructing a digital twin, mesh nodes and material density, elastic modulus, and Poisson's ratio are extracted from the product's 3D model to generate a finite element representation of the product. Simultaneously, based on the initial design parameters of the wooden packaging box, the box's geometric boundaries, board thickness distribution, and connection node stiffness matrix are generated to construct the initial finite element mesh of the wooden packaging box. The finite element representation of the product and the initial finite element mesh of the wooden packaging box are initially assembled through contact pair definitions to form the assembly's ground state. Material properties are parameterized and written into the digital twin data layer, thereby achieving a high-fidelity mapping between the product and packaging in digital space, providing a unified model foundation for subsequent coupled simulations.

[0008] When generating a vulnerability and sensitive area map for a product, this technical solution extracts multi-axis acceleration time histories from historical transportation load spectra and obtains typical load condition segments by frequency-amplitude statistical clustering. Each typical load condition segment is then applied sequentially to the product's multiphysics digital twin, and the transient dynamic equations are solved, outputting the acceleration response time series and equivalent stress contour map for each node. Nodes in the equivalent stress contour map that exceed the material's yield stress threshold are clustered as stress hotspot regions, and acceleration-sensitive parts are marked based on the peak distribution of the acceleration response time series. These stress hotspot regions and acceleration-sensitive parts are then merged to form the product vulnerability and sensitive area map. This map accurately identifies vulnerable parts of the product during transportation, enabling packaging protection design to focus on critical areas.

[0009] Furthermore, during the multi-condition simulation of the product-packaging coupling, guided by the vulnerability sensitive area map of the product, the stress hotspots on the surface of the product's multiphysics digital twin are nodally matched with the contact surfaces of the inner lining of the wooden packaging box's structural digital twin, establishing contact pair constraints. Based on these contact pair constraints, two conditions are set: an initial velocity impact condition based on drop height and a random vibration condition based on power spectral density. The explicit dynamic response of the product-packaging system is jointly solved under both conditions, recording the stress-strain state of each board element of the wooden packaging box at each load step. This approach directly transforms product vulnerability information into local boundary conditions of the packaging, significantly improving the simulation's accuracy in capturing the real damage evolution process.

[0010] Preferably, in obtaining the damage accumulation evolution sequence, the maximum principal stress and stress cycle number corresponding to each load step are extracted from the stress-strain state of each board element to construct a stress spectrum matrix. Based on the wood material cumulative damage criterion, the stress spectrum matrix is ​​transformed into a damage increment sequence according to the loading time sequence. Spatial interpolation is performed on the damage increment sequence along the board thickness direction and the plane direction to obtain a three-dimensional distribution of damage accumulation evolution per board element. The top few boards are selected as key components in descending order of damage value, and the damage accumulation evolution sequence of the key components is output. This sequence accurately reflects the progressive damage process of various parts of the packaging box under coupled loads, providing clear guidance on weak points for optimized design.

[0011] As another preferred technical solution of the present invention, when calling the parametric optimization engine for iterative optimization, the board thickness, stiffener layout parameters, and connection node stiffness of the wooden packaging box are set as optimization variables to construct a parametric model; the maximum damage value in the cumulative evolution sequence of damage of key components is used as the constraint violation degree, and the total mass of the packaging box is used as the optimization objective to establish a constraint optimization mathematical model. A sequential approximation optimization algorithm based on a radial basis function surrogate model is adopted to generate initial sample points in the variable space and calculate the true response, iteratively updating the surrogate model until a design parameter combination that satisfies the damage constraint and minimizes the mass is obtained. This method effectively achieves lightweight design while ensuring the damage resistance performance of the packaging.

[0012] During the digital twin update process, the plate thickness and stiffener layout parameters in the optimized design parameter combination are mapped to the corresponding geometric features of the wooden packaging box structure digital twin, and the finite element mesh is re-divided. The optimized connection node stiffness matrix replaces the corresponding stiffness parameters in the original structural digital twin, and the material property data layer is updated. The updated wooden packaging box structure digital twin and the product multiphysics digital twin are reassembled to generate a new round of simulation-ready coupled digital twin. This closed-loop update mechanism ensures that each optimization iteration is carried out on a model that accurately reflects the current design state.

[0013] To generate the historical transportation environmental load spectrum, this invention acquires the original triaxial acceleration signals from the actual transportation route spectrum, and extracts the characteristic parameters of the impact event and steady-state vibration segment through time-frequency analysis and envelope detection. The peak acceleration and pulse width of the impact event are proportionally mixed and reconstructed with the power spectral density of the steady-state vibration segment to generate a statistically consistent time-domain sample of the transportation environmental load, which serves as the historical transportation environmental load spectrum. This load spectrum fully reproduces the mechanical characteristics of the actual transportation conditions, providing reliable input conditions for simulation.

[0014] When setting up drop and vibration scenarios, the initial velocity at the moment of contact was calculated based on the expected drop height of the wooden packaging box and applied to the bottom nodes of the digital twin of the wooden packaging box structure. A gravity field was also applied to the entire model. The acceleration power spectral density curve of random transportation vibrations was extracted from the specifications and decomposed into load amplitudes at various frequency points. These loads were then applied to the constraint boundary points of the digital twin of the wooden packaging box structure to simulate random vibration inputs. This loading method comprehensively covers typical dynamic hazards during transportation.

[0015] In the convergence judgment phase, the optimized design parameter combination of the current round is used as input to re-execute the product-packaging coupled multi-condition simulation, obtaining the damage accumulation evolution sequence of key components in this round. The relative change rate of the maximum damage value in this round is compared with that in the previous round. If the change rate is lower than a set threshold for two consecutive times and the damage values ​​of all key components are within the safe range, the convergence condition is met; otherwise, iteration continues. This convergence criterion based on the damage change rate ensures that the output scheme achieves a stable optimal balance between protective performance and structural quality.

[0016] The technical effects and advantages provided by the present invention in the above technical solution are as follows:

[0017] By constructing a multiphysics digital twin of the product and applying historical transportation environmental load spectra for transient response simulation, a vulnerability and sensitive area map of the product is generated. This allows for the precise identification of the stress characteristics and acceleration-sensitive parts of the product under transportation vibration and impact. This map reflects the spatial heterogeneity of the product's internal structural dynamic response, replacing the conventional simplistic treatment of treating the product as a rigid body or homogeneous mass block. Based on this map, node matching and contact pair constraints are established between the product's stress hotspot areas and the contact surfaces of the wooden packaging box lining. This enables the packaging's structural digital twin to accurately bear the load transfer from the product's vulnerable parts at the contact boundaries, concentrating packaging protection resources on areas where the product truly needs cushioning and support, and avoiding over-design of non-sensitive areas.

[0018] Based on the stress-strain state extracted from multi-condition coupled simulation, a stress spectrum matrix is ​​constructed and transformed into a damage accumulation evolution sequence of key components along the time history and spatial dimension according to the wood material cumulative damage criterion. This sequence quantitatively characterizes the progressive damage process of each panel under the combined action of drop and vibration loads, and its spatiotemporal distribution information directly reflects the strength of the protective performance of each component of the packaging under the initial design parameters. The maximum damage value in this damage accumulation evolution sequence is used as the constraint violation degree, and the total mass of the packaging box is used as the optimization objective. A constraint optimization mathematical model is constructed and a parametric optimization engine is driven to perform iterative optimization. During the optimization process, design variables such as panel thickness, stiffener layout, and connection node stiffness are continuously adjusted under the feedback of the damage accumulation sequence. The surrogate model approximation algorithm guides the search direction to always be toward the parameter combination that satisfies the damage safety constraint and minimizes the mass. This closed-loop mechanism eliminates the lag and blindness of manual trial and adjustment, and realizes a quantitative trade-off and automatic optimization between structural damage resistance and lightweight objectives. Attached Figure Description

[0019] To more clearly illustrate the technical solutions in the embodiments of this application or the prior art, the drawings used in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments recorded in this invention. For those skilled in the art, other drawings can be obtained based on these drawings.

[0020] Figure 1 This is a flowchart of a simulation method for optimizing the design of wooden packaging boxes based on digital twins;

[0021] Figure 2 This is a schematic diagram illustrating the process of constructing a digital twin of the product and its wooden packaging box;

[0022] Figure 3 This is a flowchart for generating a vulnerability and sensitive area map of a product;

[0023] Figure 4 This is a flowchart of the update and optimization iteration convergence determination process for the digital twin structure of wooden packaging boxes;

[0024] Figure 5 This is a statistical chart showing the distribution of grid dimensions of wooden packaging box shell units and cross-sectional dimensions of reinforcing beam units. Detailed Implementation

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

[0026] See Figure 1 This invention provides a simulation method for optimizing the design of wooden packaging boxes based on digital twins. The method first obtains the geometric model and physical properties of the product to be packaged, as well as the initial design parameters of the wooden packaging box. Based on these, a multi-physics digital twin of the product and a structural digital twin of the packaging box are constructed. Then, a historical transportation environmental load spectrum generated using statistical consistency is applied to the product's multi-physics digital twin. The spatiotemporal distributions of the product's acceleration and stress responses are extracted through transient response simulation, generating a vulnerability sensitive area map. This vulnerability sensitive area map is mapped to the contact boundary of the wooden packaging box's structural digital twin, establishing a coupling relationship between the product and the packaging box. A multi-condition simulation of drop and vibration is then performed to obtain the damage accumulation evolution sequence of key components of the wooden packaging box. Using the damage accumulation evolution sequence as an objective function constraint, a parametric optimization engine is constructed and invoked to iteratively optimize the initial design parameters of the wooden packaging box, generating an optimized combination of design parameters. The digital twin of the wooden packaging box structure is updated based on the optimized design parameter combination. The updated digital twin of the wooden packaging box structure is then reassembled with the product multiphysics digital twin. Coupled simulation and verification are repeated until the preset convergence conditions are met, and the final wooden packaging box design scheme is output.

[0027] Example 1:

[0028] In specific implementation, please refer to Figure 2 The process of constructing a multiphysics digital twin of the product and a structural digital twin of the wooden packaging box involves three steps.

[0029] The first step involves extracting mesh nodes and material density, elastic modulus, and Poisson's ratio from the product's 3D model to generate a finite element representation of the product. The 3D model originates from the CAD model file of the product to be packaged; the file format can be STEP or IGES. Geometric cleanup is performed on the 3D model, removing non-load-bearing features such as fillets, chamfers, and threaded holes with diameters smaller than a preset threshold, resulting in a simplified product geometry. Second-order tetrahedral elements are used to mesh the simplified product geometry. The element size is set based on the overall product dimensions and the stress wave propagation wavelength, with the mesh size ranging from one-tenth to one-eighth of the minimum wavelength. After meshing, the node coordinates and element topology relationships of all elements are extracted. The node coordinates are sequentially numbered to obtain a node list, and the element topology relationships are represented as a node index list. Simultaneously, the density, elastic modulus, and Poisson's ratio of the corresponding material are read from the product material database. The node list, element topology relationships, density, elastic modulus, and Poisson's ratio are combined into a structured data volume, which constitutes the finite element representation of the product. Density, elastic modulus, and Poisson's ratio are written as material property scalars into the attribute card corresponding to each unit, and can be called and modified in parameterized form in the digital twin data layer later.

[0030] The second step involves generating the box's geometric boundary, board thickness distribution, and connection node stiffness matrix based on the initial design parameters of the wooden packaging box, thus constructing the initial finite element mesh for the box. The initial design parameters include the box's external length, width, and height; the initial thickness of each board panel; the cross-sectional shape, placement, and spacing of the reinforcing ribs; and the initial connection stiffness values ​​at nail or tenon joints. The box's geometric boundary is determined by its external dimensions, and a cuboid envelope is generated as the reference surface for the box's external boundary. The board thickness distribution is composed of the thickness values ​​of each side panel, top panel, bottom panel, end panel, and inner lining panel, with each panel assigned an independent thickness variable in the parametric definition. For each board panel, shell elements are used for discretization, with the shell's neutral surface taken from the board's geometric mid-surface. The shell element's thickness direction defines the corresponding board thickness variable. The mesh size of the shell elements is taken as one-twentieth to one-tenth of the box's minimum external dimensions. The reinforcing ribs are simulated using beam elements, with the beam element cross-sectional properties calculated based on the reinforcing rib's cross-sectional shape and size. The reinforcing rib nodes are connected to the shell element nodes through shared nodes or rigid constraints. The connection node stiffness matrix is ​​used to characterize the semi-rigid properties of the corners of the box and the connections between panels. Each connection node defines stiffness coefficients with six degrees of freedom, including three translational stiffnesses and three rotational stiffnesses. Initially, based on the connection type (e.g., nail spacing, nail diameter, and wood dowel bearing strength), the initial translational and rotational stiffness values ​​are calculated using empirical formulas and assembled into a diagonal connection node stiffness matrix. The meshes of shell and beam elements are combined, and each element is assigned corresponding material orientation, thickness, and section properties. Spring elements defined by the stiffness matrix are applied at the connections, thus constructing the initial finite element mesh of the wooden packaging box.

[0031] The third step involves initially assembling the finite element representation of the product and the initial finite element mesh of the wooden packaging box through contact pair definitions, forming the ground state of the assembly, and parameterizing the material properties into the digital twin data layer. During assembly, the node coordinate system of the product's finite element representation and the node coordinate system of the wooden packaging box's finite element mesh are unified to the same global coordinate system. Based on the relative positions of the product with the lining and blocks in the actual packaging structure, the product mesh is positioned at a preset position within the wooden packaging box's internal space. Then, contact pairs are established. The contact pair type is a functionable face-to-face contact, and the friction coefficient is taken as the measured static friction coefficient between the wood and the product surface material. One side of the contact pair is the set of potential contact area nodes on the surface of the product's finite element representation, and the other side is the set of relative nodes on the surface of the wooden packaging box's lining or block shell unit. For each contact pair, the normal contact stiffness and tangential friction coefficient are defined. The normal contact stiffness is taken as the product factor of the material's equivalent modulus under the contact surface and the element's characteristic length. After assembly, the ground state of the assembly, containing the product mesh, the wooden packaging box mesh, and the contact pair definitions, is generated. Simultaneously, the product material properties (density, modulus of elasticity, Poisson's ratio) and the wooden packaging box material properties (modulus of elasticity parallel to the grain, modulus of elasticity across the grain, shear modulus, density, yield strength, etc.) are written into the digital twin data layer in the form of parameterized name-value pairs. The digital twin data layer uses a key-value pair structure for storage, where the key is a string containing the component identifier and attribute name, and the value is the corresponding physical quantity value. Subsequently, during the iterative optimization of the initial design parameters of the wooden packaging box, modifying the corresponding parameter values ​​in the digital twin data layer can directly drive the attribute updates in the finite element model.

[0032] See Figure 5In the figure, the horizontal axis represents the element size in millimeters, and the vertical axis represents the frequency of elements within the corresponding size range. The legend distinguishes between shell element mesh sizes and beam element cross-sectional sizes. The shell element mesh sizes are concentrated in the range of 10 mm to 60 mm, exhibiting a multi-peak fluctuation trend, with the higher frequency size range mainly distributed between 30 mm and 60 mm, indicating that the shell element meshing follows the setting specification of one-twentieth to one-tenth of the minimum external dimensions of the box. The beam element cross-sectional size data are concentrated in the range of approximately 30 mm to 40 mm, with the frequency peaks clearly concentrated in this range, indicating that the stiffener cross-sectional sizes are concentrated in this range, consistent with the description in Example 1 regarding the beam element cross-sectional properties being set based on the shape and size of the stiffener cross-section. The two element size ranges partially overlap, with the shell element size having a lower frequency in the 30 mm to 40 mm range and coinciding with the beam element cross-sectional size, reflecting the design requirement that shell elements and beam elements share nodes at intersections during meshing. Overall, the figure specifically illustrates the distribution characteristics of mesh size and cross-sectional size of shell and beam elements in the digital twin of the wooden packaging box structure in Example 1, verifying the rationality of the parametric design and mesh generation strategy.

[0033] Example 2:

[0034] In specific implementation, please refer to Figure 3 The process of generating a vulnerability and sensitive area map of a product involves constructing a historical transportation environment load spectrum, extracting typical load condition segments, solving the transient dynamics of the product's multiphysics digital twin, and clustering and labeling vulnerable areas.

[0035] The historical transportation environmental load spectrum is generated as follows: Raw triaxial acceleration signals from the actual transportation route spectrum are collected. During the acquisition process, a triaxial accelerometer is fixedly mounted on a rigid bracket on the floor of the transport vehicle near the actual placement position of the wooden packing crates. The acceleration amplitude in three orthogonal directions throughout the transportation process is continuously recorded at a sampling frequency of 512 Hz. These three directions are denoted as longitudinal, lateral, and vertical. After acquisition, time-frequency analysis and envelope detection are performed on the raw triaxial acceleration signals. Time-frequency analysis uses a short-time Fourier transform with a Hanning window as the window function. The window length is set to 1024 sampling points, and the overlap rate is set to 50%, yielding the spectral distribution at each time point. Envelope detection uses the Hilbert transform method to calculate the magnitude of the analytic signal from each of the three raw triaxial acceleration signals, obtaining the acceleration envelope. Impact events are detected on the envelope using a preset amplitude threshold, which is three times the root mean square value of the acceleration throughout the transportation process. Whenever the envelope amplitude exceeds a preset amplitude threshold, the time interval between the threshold threshold start point and the fall back below the threshold endpoint is recorded as an impact event, and the peak acceleration and pulse width parameters of each impact event are extracted. The signal segment between adjacent impact events is divided into steady-state vibration segments, and the power spectral density of each steady-state vibration segment is calculated using the Welch method. The segment length is set to 2048 sampling points, the overlap rate is set to 50%, and the Hanning window function is used. The power spectral densities of multiple steady-state vibration segments are arithmetically averaged to obtain the average power spectral density of the steady-state vibration segment. Thus, the peak acceleration and pulse width of the impact event, and the power spectral density of the steady-state vibration segment constitute the characteristic parameters of the historical transportation environment load spectrum.

[0036] The peak acceleration, pulse width, and power spectral density of the steady-state vibration segment of the impact event are proportionally mixed and reconstructed to generate a statistically consistent time-domain sample of the transportation environmental load, which serves as the historical transportation environmental load spectrum. During the mixed reconstruction, the power spectral density of the steady-state vibration segment is transformed into a steady-state random vibration time-domain sequence using inverse Fourier transform. A random phase angle is introduced during the transformation, and this random phase angle follows a uniform distribution in the interval [0, 2π]. The impact events are inserted into the steady-state random vibration time-domain sequence according to the time intervals of their original occurrence sequence. The waveform of each impact event is approximated using a half-sine pulse, with the amplitude of the half-sine pulse taken as the peak acceleration of the corresponding impact event and the duration taken as the pulse width of the corresponding impact event. During insertion, a weighted smooth transition is applied between the start and end times of the half-sine pulse and the corresponding times in the steady-state random vibration time-domain sequence. The width of the transition superposition region is one-tenth of the half-sine pulse width, and the weighted smoothing method is a linear cross-gradual transition, resulting in the mixed-reconstructed transportation environmental load time-domain sample. This time-domain sample of the transportation environmental load maintains complete consistency with the original actual transport route spectrum in three statistical dimensions: impact amplitude distribution, impact interval distribution, and steady-state frequency domain energy distribution.

[0037] After generating the historical transportation environmental load spectrum, multi-axis acceleration time histories are extracted from it. Since the hybrid reconstructed transportation environmental load time-domain sample itself is a continuous multi-axis acceleration time history covering several hours, several multi-axis acceleration time history segments are extracted from this sample according to a fixed time window. Each segment is 60 seconds long, and the time interval between adjacent segments is 30 seconds. Then, typical load condition segments are obtained by frequency-amplitude statistical clustering. During clustering, for each multi-axis acceleration time history segment, the root mean square vector of the acceleration in the three axes within the segment is calculated, and a one-dimensional energy index is synthesized. The formula for calculating the one-dimensional energy index is the square root of the sum of the squares of the three-axis root mean square values. The one-dimensional energy indices of all segments are used as clustering features, and K-means clustering is performed using a K-value of three, with Euclidean distance as the distance metric. After clustering, the multi-axis acceleration time history segment corresponding to the central feature value of each class is selected as the typical load condition segment of that class, resulting in three typical load condition segments, corresponding to low-intensity vibration condition, medium-intensity vibration condition and high-intensity impact condition, respectively.

[0038] Typical load case segments are sequentially applied to the product's multiphysics digital twin. During loading, the triaxial acceleration time series extracted from three multiaxial acceleration time history segments are used as acceleration boundary conditions, uniformly applied as body loads to the three translational degrees of freedom of all nodes in the product's multiphysics digital twin. When solving the transient dynamic equations, the Newmark implicit time integration method is used, with the time step being one-fifth of the sampling interval of the triaxial acceleration time series in the typical load case segment. Rayleigh damping is used as the damping model, and the mass coefficient and stiffness coefficient of Rayleigh damping are calculated from the product's first two modal frequencies and a set modal damping ratio, which is set to 500%. At each time step, the triaxial acceleration response component values ​​of all nodes in the product's multiphysics digital twin are output and stored as an acceleration response time series according to node number. Simultaneously, the equivalent stress component values ​​at all element integration points of each time step are output and stored as an equivalent stress time series according to element number. The maximum value of the equivalent stress time series is taken in the time domain to obtain the maximum equivalent stress value of each unit in the entire loading process. The maximum equivalent stress values ​​of all units are then mapped to the surface of the product's multiphysics digital twin to form an equivalent stress cloud map.

[0039] Nodes exceeding the material yield stress threshold in the equivalent stress cloud map are clustered into stress hotspot regions. The material yield stress threshold is taken as the tensile yield strength of the product material under static load conditions at room temperature. All surface nodes in the equivalent stress cloud map are traversed, and nodes with the maximum equivalent stress value exceeding the material yield stress threshold are selected as stress over-threshold nodes. Density-based spatial clustering is performed on discrete stress over-threshold nodes. The spatial clustering algorithm used is DBSCAN, with a neighborhood radius three times the average size of the surface units of the product's multiphysics digital twin, and a minimum number of neighboring nodes of three. Each cluster of stress over-threshold nodes constitutes a stress hotspot region. Simultaneously, acceleration-sensitive parts are marked based on the peak distribution of the acceleration response time series. For the triaxial acceleration response time series of each node, the absolute values ​​of the three axial acceleration peaks are extracted, and the maximum absolute value of the three axial acceleration peaks is taken as the acceleration peak response of that node. All nodes are sorted from largest to smallest acceleration peak responses, and the nodes in the top 15% of acceleration peak responses are selected as acceleration-sensitive nodes. Density-based spatial clustering is also performed on discrete acceleration-sensitive nodes. The neighborhood radius is three times the average size of the surface unit of the product's multiphysics digital twin, and the minimum number of neighborhood nodes is three. Each cluster of acceleration-sensitive nodes constitutes an acceleration-sensitive part.

[0040] After obtaining stress hotspot regions and acceleration-sensitive areas, these regions are merged to form a product vulnerability sensitive area map. During merging, for cases where the same node may belong to both a stress hotspot region and an acceleration-sensitive area, overlapping nodes are retained only once in the final product vulnerability sensitive area map, and the overlapping area is simultaneously marked with both stress hotspot type and acceleration-sensitive type. The product vulnerability sensitive area map is represented on the 3D surface of the product's multiphysics digital twin using zoned color annotations, with different colors corresponding to pure stress hotspot regions, pure acceleration-sensitive areas, and stress-acceleration coupled sensitive regions, respectively.

[0041] Example 3:

[0042] In practice, the process of performing multi-condition simulation of drop and vibration of product-packaging coupling is divided into three parts: contact pair constraint establishment, condition setting, and joint solution.

[0043] When establishing contact pair constraints, the vulnerability-sensitive area map of the product is used as a guide to perform node matching between the stress hotspot regions on the surface of the product's multiphysics digital twin and the inner lining contact surface of the wooden packaging box structure digital twin. The stress hotspot region on the surface of the product's multiphysics digital twin is a continuously distributed set of nodes, each with an independent identifier detached from the product's multiphysics digital twin. The inner lining contact surface of the wooden packaging box structure digital twin is composed of the surface node set of the inner lining shell unit. The matching process uses a minimum distance projection algorithm. For each node in the stress hotspot region, the shortest distance from that node to all surface projection points of the shell units on the inner lining contact surface is calculated. If the shortest distance is less than a preset gap tolerance, a contact pair is established between the node in the stress hotspot region and the surface node of the shell unit where the shortest projection point is located on the inner lining contact surface. The preset gap tolerance is one-tenth of the average side length of the inner lining shell unit. Each contact pair comprises a master surface node and a slave surface node. The master surface node is taken from the product surface node in the stress hotspot region, and the slave surface node is taken from the corresponding shell element surface node on the contact surface of the wooden crate lining. After establishing the contact pairs, contact attributes are set for each pair. Contact attributes include normal contact behavior and tangential contact behavior. The normal contact behavior adopts the penalty function method, and the normal contact stiffness is taken as the product of the harmonic mean of the elastic modulus of the bottom material in the stress hotspot region and the cross-grain elastic modulus of the lining wood, multiplied by the characteristic side length of the master surface element. The tangential contact behavior adopts the Coulomb friction model, and the friction coefficient is taken as the static friction coefficient measured by the inclined plane sliding experiment between the product surface material and the wooden crate lining wood. The experimentally measured static friction coefficient value is selected in the range of 0.35 to 0.55. After the contact pair constraints are set, a data block containing all contact pair definitions is generated, and the data block is written into the corresponding digital twin data layer of the assembly.

[0044] Based on contact pair constraints, two conditions are set: an initial velocity impact condition based on drop height and a random vibration condition based on power spectral density.

[0045] When setting the initial velocity impact condition based on the drop height, first convert the expected drop height of the wooden packaging box into the initial velocity at the moment of contact. The conversion relationship is expressed by the following formula:

[0046]

[0047] in, The initial velocity applied at the moment of contact to the bottom node of the digital twin of the wooden packaging box structure is expressed in millimeters per second. This represents the acceleration due to gravity, taken as a constant value of 9800 mm / s². This indicates the expected drop height of the wooden crate, in millimeters. (Expected drop height of wooden crate) The value is determined based on the total mass rating of the package, and the specific rules are as follows: When the total mass of the package is less than or equal to 20 kg, the expected drop height is... Take 800 mm; when the total weight of the package is greater than 20 kg and less than or equal to 50 kg, the expected drop height is... Take 600 mm; when the total weight of the package exceeds 50 kg, the expected drop height is... Take 400 mm. Obtain the initial velocity at the moment of contact. Subsequently, this initial velocity is applied as a concentrated initial velocity load to the bottom nodes of the digital twin of the wooden packaging box structure. The velocity vector direction is vertically downwards, i.e., along the component in the global coordinate system aligned with the direction of gravity. Simultaneously, a gravitational field is applied as a uniform volume load to all elements of both the product's multiphysics digital twin and the wooden packaging box structure's digital twin. The magnitude and direction of the gravitational acceleration are the same as described above. Consistent. The simulation time for the drop condition is taken from the moment of contact with the product upon impact until the damping of repeated collisions between the product and the wooden packaging box decreases to the point where the system's kinetic energy is less than one-thousandth of the initial kinetic energy. The maximum simulation time is preset to 0.1 seconds.

[0048] When setting up random vibration conditions based on power spectral density, the acceleration power spectral density curve of random vibration in transportation is extracted from the transportation test specifications. The transportation test specifications adopt the random vibration test spectrum of highway transportation specified in the national standard. This spectrum is expressed as a piecewise linear form of frequency-power spectral density values, covering a frequency range of 1 Hz to 200 Hz. The acceleration power spectral density curve is discretized into 256 frequency points in a non-uniform manner. During discretization, the frequency points are denser in sections with rapid changes in power spectral density and sparser in sections with gradual changes. Each frequency point corresponds to a frequency value and a power spectral density value. For each frequency point, the load amplitude at the corresponding frequency point is calculated using a set frequency resolution. The load amplitude is taken as the square root of the product of the power spectral density value and the frequency resolution. The load amplitudes at the 256 frequency points are synthesized into a time-domain sample of random vibration acceleration through inverse Fourier transform. A random phase angle is introduced into the inverse Fourier transform, and the random phase angle follows a uniform distribution in the interval of 0 to 2π. The generated random vibration acceleration time-domain sample lasts for 60 seconds with a time step of 0.001 seconds. The random vibration acceleration time-domain sample is applied to the constraint boundary points of the digital twin of the wooden packaging box structure in the form of acceleration basic excitation. The constraint boundary points are selected from the node set in the contact area between the bottom surface of the wooden packaging box and the transport platform. The acceleration time history in the corresponding direction is applied to the three translational directions at the constraint boundary points, while the three rotational directions are constrained and fixed.

[0049] The explicit dynamic response of the product-packaging system is solved jointly under two operating conditions. During the joint solution, the initial velocity impact load of the drop condition and the acceleration base excitation load of the random vibration condition are sequentially concatenated into a complete combined load time history. The first segment of the combined load time history is the drop impact segment, the duration of which is determined by the simulation time of the drop condition; the second segment is the random vibration segment, with a duration of 60 seconds. A 0.2-second zero-load transition segment is set between the two load segments to allow sufficient attenuation of the residual vibration after the drop impact. The central difference explicit time integration method is used to solve the motion equations of the product-packaging system for the combined load time history. The time step is half the ratio of the characteristic length of the smallest element of the product's multiphysics digital twin to the material's sound velocity, and does not exceed one-fifth of the sampling interval of the combined load time history to ensure numerical stability. After each load step is solved, all plate shell elements in the digital twin of the wooden packaging box structure are traversed. The stress and strain components at the center integration points of the upper and lower surfaces of each shell element are obtained. The six stress components and six strain components are stored according to the load step number, element number, and integration point position. The stress-strain state of each plate element of the wooden packaging box under each load step is recorded. The storage format is a structured data table consisting of load step number, element number, integration point identifier, six stress tensor components, and six strain tensor components. This data table will serve as the basic input for subsequent damage accumulation calculations.

[0050] Example 4:

[0051] In practical implementation, the process of extracting the maximum principal stress and stress cycle number corresponding to each load step from the stress-strain state of each plate element and constructing the stress spectrum matrix is ​​as follows: For each plate shell element in the digital twin of the wooden packaging box structure, the six stress components at the center integration points of the upper and lower surfaces of the plate shell element are sequentially read from the recorded stress-strain state data table for each load step. The six stress components are the normal stress components. , , and shear stress components , , Using the six stress components mentioned above, the three principal stress values ​​at the integration point of each load step are calculated by solving a cubic equation. The value with the largest value is selected as the maximum principal stress, denoted as . For each load step, the maximum principal stress of the plate shell element is taken as the larger of the maximum principal stress at the center integration point of the upper surface and the maximum principal stress at the center integration point of the lower surface. This yields the time series of the maximum principal stress of the plate shell element over the entire combined load time history.

[0052] Rainflow counting was performed on the maximum principal stress time series to extract all closed stress cycles. During rainflow counting, the maximum principal stress time series was first filtered for extreme points to remove minute fluctuations with amplitudes less than 0.5% of the global maximum peak value in the maximum principal stress time history. Then, a four-point rainflow counting rule was applied to the filtered extreme point series to extract each complete stress cycle one by one. For each extracted stress cycle, the stress amplitude of the stress cycle was recorded. mean stress The actual number of stress cycles, n, is also considered. After the rainflow counting is complete, all stress cycle information is combined into a stress spectrum matrix. Each row of the stress spectrum matrix corresponds to a stress cycle, and the stress amplitude is stored sequentially within each row. mean stress And the number of iterations, n.

[0053] Based on the cumulative damage criterion for wood materials, the stress spectrum matrix is ​​transformed into a damage increment sequence according to the loading time sequence. The loading time sequence refers to the order in which each stress cycle is extracted from the time series of the maximum principal stress. The cumulative damage criterion for wood materials adopts the linear cumulative damage theory. The damage increment generated by each stress cycle is calculated by the following formula, and the damage increments are arranged in chronological order:

[0054]

[0055] in, This represents the damage increment generated by the p-th stress cycle arranged in the loading sequence, and is dimensionless. The actual number of occurrences of the p-th stress cycle is obtained by the rainflow counting method; This represents the allowable number of cycles corresponding to the p-th stress cycle on the fatigue life curve of the wood material, and is dimensionless. Allowable number of cycles The method for determining it is: based on the stress amplitude of the p-th stress cycle. and stress ratio Find the SN curve family for wood materials. The SN curve family for wood materials is based on stress ratio. The parameters are obtained by fitting multiple sets of constant amplitude fatigue test data, and the SN curves of each set are expressed in power function form. .in Indicates the stress ratio as The fatigue strength coefficient at that time, expressed in megapascals (MPa). Indicates the stress ratio as The fatigue strength index at stress p is dimensionless. It is based on the stress amplitude during the p-th stress cycle. Reverse calculation ,like If it is zero, then the corresponding We directly set the value to zero. After calculating the damage increment for all stress cycles, we obtain the damage increment sequence.

[0056] Spatial interpolation is performed along the thickness and planar directions of the damage increment sequence to obtain a three-dimensional distribution of damage accumulation evolution per plate. Spatial interpolation along the thickness direction uses the damage increments at multiple integration points of the shell element along the thickness direction. For each plate shell element, five Gaussian integration points are set in the thickness direction of the standard shell element, numbered from top to bottom as integration point one to integration point five. Using the damage increments at these five Gaussian integration points, the damage increment at any thickness coordinate within the shell element is obtained through Lagrange interpolation. Spatial interpolation along the planar direction uses the damage increment values ​​at the nodes of the plate shell element, and the damage increment at any point within the plate surface is obtained through two-dimensional shape function interpolation. The interpolated damage increments are successively accumulated in the time domain to obtain a three-dimensional spatial distribution of the cumulative damage value per plate, with each coordinate point corresponding to a total damage value. Total damage value This is the sum of the damage increments for all load steps at that point.

[0057] The method of selecting the top few plates as key components in descending order of damage value is as follows: For each independent plate, calculate the total damage value at all spatial interpolation points within all the shell elements contained in that plate. The maximum value is taken as the representative damage value of the plate. All plates are sorted from largest to smallest representative damage values, and the top M plates with the largest representative damage values ​​are selected as critical components, with M set to 5. The reason for choosing M as 5 is that in common failure modes of wooden crates, damage is concentrated in a few bottom load-bearing plates, side support plates, or corner connection plates; choosing 5 ensures coverage of high-risk areas along all major load-bearing paths. For each selected critical component, the damage increment sequence under each load step or stress cycle, as well as the final accumulated three-dimensional damage distribution data, are output, constituting the damage accumulation evolution sequence of the critical component.

[0058] When using the parametric optimization engine to iteratively optimize the initial design parameters of the wooden packaging box, the thickness of the board material, the layout parameters of the stiffeners, and the stiffness of the connection nodes are set as optimization variables to construct a parametric model. The board thickness optimization variable includes the thickness value of each board piece, with a continuous value range of 6 mm to 25 mm. The stiffener layout parameters include the stiffener spacing and the stiffener cross-sectional dimensions. The stiffener spacing is a discrete variable with a step size of 50 mm, ranging from 100 mm to 400 mm. The stiffener cross-sectional dimensions include the stiffener width and stiffener height, with the stiffener width ranging from 20 mm to 60 mm and the stiffener height ranging from 10 mm to 40 mm, both being continuous variables. The connection node stiffness optimization variables include the translational stiffness and rotational stiffness proportionality coefficients of each connection node. The proportionality coefficients are multiples of the initial connection node stiffness, with a continuous value range of 0.5 to 2.0. These optimization variables are denoted as the design variable vector x, and a parameterized model is established so that it can automatically update the board thickness distribution, stiffener geometry and position, and connection node stiffness matrix of the digital twin of the wooden packaging box structure based on any set of values ​​of x.

[0059] The maximum damage value in the cumulative damage evolution sequence of key components is used as the constraint violation degree, and the total mass of the packaging box is used as the optimization objective to establish a constrained optimization mathematical model. The maximum damage value is taken as the maximum value among the representative damage values ​​of all key components, denoted as . Total weight of packaging boxes The wood density is calculated by multiplying the volume of each panel by the wood density and summing the results. The wood density is the measured air-dry density of the tree species used in the wooden packaging boxes. The constrained optimization mathematical model is expressed as: minimizing the objective function... (x), satisfying the constraints. (x)-1.0≤0, and the design variable vector x is within the range of values ​​of each variable.

[0060] A sequential approximation optimization algorithm based on a radial basis function surrogate model is used for the solution. The radial basis function surrogate model uses a Gaussian kernel function as the basis function, with the kernel function taking the form φ(r) = exp(-(ε·r)²), where r is the Euclidean distance between points, and the shape parameter ε is determined by minimizing the prediction error using leave-one-out cross-validation. First, an initial set of sample points is generated in the variable space using Latin hypercube sampling, with the number of initial sample points being ten times the dimension of the design variable vector. Then, a product-packaging coupled multi-condition simulation is used to calculate the total mass of the actual packaging box corresponding to each initial sample point. and maximum damage value This forms the initial training dataset. Using the initial training dataset, the total mass of the packing boxes is constructed. Radial basis function surrogate model and maximum damage value The radial basis function surrogate model is then used. The sequential approximation iteration proceeds: on the two current surrogate models, the constraint expectation improvement criterion is used as an addition strategy to find the point with the maximum expected improvement that simultaneously satisfies the expected reduction objective function and the constraints. The constraint expectation improvement value of the candidate points is calculated, and the point with the largest constraint expectation improvement value is selected as the next expensive simulation point. The true value of this point is obtained by re-executing the product-packaging coupled multi-condition simulation. and The value is calculated and this new sample point is added to the training dataset. The two radial basis function surrogate models are retrained using the updated training dataset, completing one iteration. This iterative process is repeated until the maximum damage value is obtained from two consecutive iterations. All values ​​are less than or equal to 1.0 and the total weight of the packaging box is less than or equal to 1.0. If the relative decrease is less than two per thousand, the optimization is considered to have converged, and the current combination of design parameters that satisfies the damage constraint and has the minimum quality is output.

[0061] Example 5:

[0062] In specific implementation, please refer to Figure 4 The process of updating the digital twin of the wooden packaging box structure based on the optimized design parameter combination includes geometric parameter mapping, mesh re-division, stiffness parameter replacement, data layer update, and reassembly.

[0063] The optimized design parameter combination maps the plate thickness parameters to the corresponding geometric features of the digital twin of the wooden packaging box structure. The optimized design parameter combination includes optimized plate thickness values ​​for each plate component, which directly replace the corresponding plate thickness variables in the initial design parameters of the wooden packaging box. For each plate component, the thickness attribute value of the corresponding shell element in the digital twin of the wooden packaging box structure is modified to the optimized plate thickness value. The stiffener layout parameters are then mapped to the corresponding geometric features of the digital twin of the wooden packaging box structure. The optimized stiffener layout parameters include stiffener spacing, stiffener cross-sectional width, and stiffener cross-sectional height. The stiffener spacing value determines the positioning of the stiffener in the plate plane, while the stiffener cross-sectional width and height together determine the cross-sectional shape and size of the stiffener. During mapping, based on the optimized stiffener spacing, the positioning lines of the stiffeners are regenerated in the geometric model of the corresponding plate component. The spacing of the positioning lines is equal to the optimized stiffener spacing value, and the cross-sectional moment of inertia and cross-sectional area properties of the beam element are recalculated based on the optimized stiffener cross-sectional width and height.

[0064] After completing the geometric parameter mapping, the digital twin of the wooden packaging box structure is re-meshed. During re-meshing, the global mesh size reference remains unchanged, still taking one-twentieth to one-tenth of the box's minimum external dimensions. For cases where changes in stiffener spacing lead to variations in stiffener positioning lines, shell element nodes and stiffener beam element nodes are forced to share nodes at intersections during the re-meshing process to ensure displacement coordination. After re-meshing, a new list of node coordinates, element topology relationships, and boundary condition sets are generated, replacing the original mesh data in the wooden packaging box structure's digital twin.

[0065] The optimized connection node stiffness matrix replaces the corresponding stiffness parameters in the digital twin of the log packaging box structure. The optimized connection node stiffness matrix is ​​obtained by multiplying the optimized translational stiffness scaling factor and rotational stiffness scaling factor by the corresponding translational and rotational stiffness values ​​in the initial connection node stiffness matrix, respectively. For each connecting spring element in the digital twin of the log packaging box structure, the six degrees of freedom stiffness coefficients corresponding to that connecting spring element are read. The optimized translational stiffness scaling factor is multiplied by the original translational stiffness coefficient, and the optimized rotational stiffness scaling factor is multiplied by the original rotational stiffness coefficient to obtain the updated connection node stiffness coefficient, which is then written into the properties of the connecting spring element.

[0066] Subsequently, the material property data layer is updated. In the digital twin data layer, the values ​​of the plate thickness parameter are modified to the optimized plate thickness values, the values ​​of the stiffener cross-section width and stiffener cross-section height are modified to the optimized corresponding values, and the translational stiffness ratio coefficient and rotational stiffness ratio coefficient of the connection node stiffness, as well as the latest values ​​of each calculated stiffness component, are written into the corresponding key-value pair storage entries.

[0067] The updated digital twin of the wooden packaging crate structure and the product multiphysics digital twin are reassembled to generate a new simulation-ready coupled digital twin. During reassembly, the same positioning method as the initial assembly is used to position the product multiphysics digital twin to a preset position inside the updated wooden packaging crate structure digital twin. Then, contact pair constraints are re-established, with the same rules as the initial assembly. The friction coefficient is the static friction coefficient between the product surface material and the inner lining wood of the wooden packaging crate. The normal contact stiffness is the product of the harmonic mean of the elastic modulus of the bottom material of the stress hotspot region and the transverse elastic modulus of the inner lining wood multiplied by the characteristic side length of the principal surface element. After assembly, the coupled digital twin data is written to the simulation-ready state storage area.

[0068] The specific steps for repeating simulation verification until the convergence condition is met are as follows: Using the optimized design parameter combination of the current round as input, re-execute the product-packaging coupled drop and vibration multi-condition simulation on the updated coupled digital twin. The simulation process includes setting an initial velocity impact condition based on drop height and a random vibration condition based on power spectral density, jointly solving the explicit dynamic response of the product-packaging system, and recording the stress-strain state of each board element of the wooden packaging box under each load step. Extract the maximum principal stress and stress cycle number corresponding to each load step from the stress-strain state of each board element, construct a stress spectrum matrix, and convert the stress spectrum matrix into a damage increment sequence according to the loading time sequence based on the wood material cumulative damage criterion. Spatial interpolation is performed on the damage increment sequence along the thickness direction and the plane direction to obtain the three-dimensional distribution of damage accumulation evolution per board element in this round. The top 5 boards in descending order of damage value are selected as the key components of this round, the damage accumulation evolution sequence of the key components in this round is output, and the maximum damage value of this round is extracted from this sequence.

[0069] The relative rate of change between the maximum damage value in this round and the maximum damage value in the previous round is calculated using the following formula:

[0070]

[0071] in, Indicates the iteration round number. It is an integer greater than or equal to 2; Indicates the first Maximum damage value of wheel relative to the first The relative rate of change of the maximum damage value of the wheel, dimensionless; Indicates the first The maximum damage value obtained from wheel simulation is dimensionless. Indicates the first The maximum damage value obtained from the simulation is dimensionless. If the first... Wheel and the first The relative rate of change of the wheel Less than the set threshold And the first Wheel and the first The relative rate of change of the wheel It is also less than the set threshold. If the relative rate of change is lower than a set threshold for two consecutive occurrences, then it is considered that the relative rate of change is lower than a set threshold for two consecutive occurrences. (Set threshold) The value of 0.02 is based on the fact that when the relative change of the maximum damage value is less than two percent, the structural response of the wooden packaging box has reached a stable state, and the changes in design parameters caused by continued iteration are minimal, which meets the engineering convergence requirements.

[0072] In addition to determining that two consecutive relative rates of change are below a set threshold, it is also necessary to verify that the damage values ​​of all critical components are within a safe range. The safe range is defined as the cumulative total damage value of all critical components. All are less than or equal to 1.0. Cumulative total damage value The maximum total damage value of all spatial interpolation points within the corresponding plate is used. If this maximum value does not exceed 1.0, the critical component is considered to be within the safe range. If both conditions are met simultaneously—that the relative change rate is less than 0.02 for two consecutive iterations and that the damage values ​​of all critical components are within the safe range—then the convergence condition is met, the iteration stops, and the currently optimized design parameter combination is output as the final wooden packaging box design scheme. If either condition is not met, the optimized design parameter combination of the current round is used as the initial design parameters for the next round of iteration. The parametric optimization engine's optimization process is re-executed, the digital twin of the wooden packaging box structure is updated again, and coupled multi-condition simulation is performed. This process is repeated until the convergence condition is met.

[0073] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any changes or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application.

Claims

1. A simulation method for optimizing the design of wooden packaging boxes based on digital twins, characterized in that, include: Obtain the geometric model and physical properties of the product to be packaged, as well as the initial design parameters of the wooden packaging box, and construct a multiphysics digital twin of the product and a structural digital twin of the wooden packaging box; Apply historical transportation environment load spectrum to the product's multiphysics digital twin, extract the spatiotemporal distribution of the product's acceleration and stress response through transient response simulation, and generate a vulnerability sensitive area map of the product; The vulnerability-sensitive area map of the product is mapped to the contact boundary of the digital twin of the wooden packaging box structure. The drop and vibration multi-condition simulation of the product-packaging coupling is performed to obtain the damage accumulation evolution sequence of the key components of the wooden packaging box. Using the damage accumulation evolution sequence as the objective function constraint, the parametric optimization engine is invoked to iteratively optimize the initial design parameters of the wooden packaging box, generating an optimized combination of design parameters; The digital twin of the wooden packaging box structure is updated based on the optimized design parameter combination, and the simulation verification is repeated until the convergence condition is met, and the final wooden packaging box design scheme is output.

2. The simulation method for optimized design of wooden packaging boxes based on digital twins according to claim 1, characterized in that, The steps for constructing a multiphysics digital twin of the product and a structural digital twin of the wooden packaging box include: Extract mesh nodes and material density, elastic modulus and Poisson's ratio from the 3D model of the product to generate the finite element representation of the product. Based on the initial design parameters of the wooden packaging box, the geometric boundary of the box, the thickness distribution of the board and the stiffness matrix of the connection nodes are generated, and the initial finite element mesh of the wooden packaging box is constructed. The finite element representation of the product and the initial finite element mesh of the wooden packaging box are initially assembled through contact pair definition to form the ground state of the assembly, and the material properties are parameterized and written into the digital twin data layer.

3. The simulation method for optimized design of wooden packaging boxes based on digital twins according to claim 1, characterized in that, The steps for generating the product vulnerability and sensitive area map include: Multi-axis acceleration time history is extracted from historical transportation environment load spectrum, and typical load condition segments are obtained by frequency-amplitude statistical clustering. Typical load conditions are applied sequentially to the multiphysics digital twin of the product, the transient dynamic equations are solved, and the acceleration response time series and equivalent stress contour map of each node are output. Nodes exceeding the material yield stress threshold in the equivalent stress cloud map are clustered into stress hotspot regions. Acceleration-sensitive parts are marked based on the peak distribution of the acceleration response time series. The stress hotspot regions and acceleration-sensitive parts are then merged to form a vulnerability-sensitive area map of the product.

4. The simulation method for optimized design of wooden packaging boxes based on digital twins according to claim 3, characterized in that, The steps for performing multi-condition simulation of drop and vibration coupling of product-packaging include: Guided by the vulnerability-sensitive area map of the product, the stress hotspots on the surface of the product's multiphysics digital twin are matched with the contact surfaces of the inner lining of the wooden packaging box structure digital twin to establish contact pair constraints. Based on contact pair constraints, two conditions are set: an initial velocity impact condition based on drop height and a random vibration condition based on power spectral density. The explicit dynamic response of the product-packaging system is solved jointly under two working conditions, and the stress-strain state of each board element of the wooden packaging box is recorded under each load step.

5. The simulation method for optimized design of wooden packaging boxes based on digital twins according to claim 4, characterized in that, The steps for obtaining the damage accumulation evolution sequence of key components of the wooden packaging box include: The maximum principal stress and stress cycle number corresponding to each load step are extracted from the stress-strain state of each plate element to construct a stress spectrum matrix; Based on the cumulative damage criterion for wood materials, the stress spectrum matrix is ​​transformed into a damage increment sequence according to the loading time sequence; Spatial interpolation is performed on the damage increment sequence along the plate thickness direction and the plane direction to obtain the three-dimensional distribution of damage accumulation evolution in units of plates. The top few plates are selected as key components in descending order of damage value, and the damage accumulation evolution sequence of the key components is output.

6. The simulation method for optimized design of wooden packaging boxes based on digital twins according to claim 5, characterized in that, The steps of iteratively optimizing the initial design parameters of the wooden packaging box by calling the parametric optimization engine include: The thickness of the wooden packaging box, the layout parameters of the reinforcing ribs, and the stiffness of the connection nodes are set as optimization variables to construct a parametric model; The maximum damage value in the cumulative evolution sequence of damage of key components is used as the constraint violation degree, and the total mass of the packaging box is used as the optimization objective to establish a constraint optimization mathematical model. A sequential approximation optimization algorithm based on a radial basis function surrogate model is adopted to generate initial sample points in the variable space and calculate the true response. The surrogate model is iteratively updated until the design parameter combination that satisfies the damage constraint and minimizes the quality is obtained.

7. The simulation method for optimized design of wooden packaging boxes based on digital twins according to claim 1, characterized in that, The steps of updating the digital twin of the wooden packaging box structure based on the optimized design parameters include: The plate thickness and stiffener layout parameters in the optimized design parameter combination are mapped to the corresponding geometric features of the digital twin of the wooden packaging box structure, and the finite element mesh is redefined. Replace the corresponding stiffness parameters in the digital twin of the log packaging box structure with the optimized connection node stiffness matrix, and update the material property data layer; The updated digital twin of the wooden packaging box structure is reassembled with the product's multiphysics digital twin to generate a new simulation-ready coupled digital twin.

8. The simulation method for optimized design of wooden packaging boxes based on digital twins according to claim 1, characterized in that, The historical transportation environment load spectrum is generated as follows: The original triaxial acceleration signals in the actual operation spectrum were collected, and the characteristic parameters of the impact event and steady-state vibration segment were extracted by time-frequency analysis and envelope detection. The peak acceleration, pulse width, and power spectral density of the steady-state vibration segment of the impact event are proportionally mixed and reconstructed to generate a statistically consistent time-domain sample of the transportation environmental load, which serves as the historical transportation environmental load spectrum.

9. The simulation method for optimized design of wooden packaging boxes based on digital twins according to claim 4, characterized in that, The steps of setting the initial velocity impact condition based on drop height and the random vibration condition based on power spectral density include: The initial velocity at the moment of contact is calculated based on the expected drop height of the wooden packaging box and applied to the bottom nodes of the digital twin of the wooden packaging box structure. The gravity field is then applied to the entire model. The acceleration power spectral density curve of random vibration during transportation is extracted from the specifications, decomposed into the load amplitude at each frequency point, and applied to the constraint boundary points of the digital twin of the wooden packing box structure to simulate random vibration input.

10. The simulation method for optimized design of wooden packaging boxes based on digital twins according to claim 1, characterized in that, The steps of repeating simulation verification until the convergence condition is met include: Using the optimized design parameter combination of the current round as input, the product-packaging coupled multi-condition simulation is re-executed to obtain the damage accumulation evolution sequence of the key components in this round; Compare the relative rate of change of the maximum damage value in this round with the maximum damage value in the previous round. If the rate of change is lower than the set threshold for two consecutive times and the damage values ​​of all key components are within the safe range, the convergence condition is met; otherwise, continue iterating.

Images