Method and system for optimizing degradable tableware material formulation

Abstract

This invention relates to the field of computer-aided process design technology, specifically to a method and system for optimizing the formulation of biodegradable tableware materials. The method includes the following steps: generating equilibrium molecular conformations based on atomic coordinates and calculating initial formulation component parameters; converting modulus parameters to generate theoretical elastic recovery strain parameters; mapping models to calculate deformation and marking diffusion characteristic length parameters; calculating diffusion coefficients and flux parameters to update local hydrogen ion concentration indices; updating reaction rate constants and integrating them to generate finalized formulation structural parameters. In this invention, a microscopic molecular dynamics model is constructed to analyze interface voids to quantitatively assess compatibility; relaxation time spectra are used to invert elastic recovery strain to predict deformation risk; autocatalytic reaction and diffusion theories are combined to track hydrogen ion concentration evolution; a dynamic mapping between degradation and failure is established; and full life-cycle verification is completed in a virtual environment, significantly shortening the R&D cycle and reducing trial-and-error costs.

Description

TECHNICAL FIELD

[0001] The application relates to the technical field of computer-aided process design, in particular to a degradable tableware material formula optimization method and system. BACKGROUND

[0002] The technical field of computer-aided process design mainly involves assisting process designers in completing the conversion process from product design data to manufacturing instruction data by using computer software and hardware technology. The core matters of this field include the formulation of product manufacturing process routes, the detailed design of process contents, the selection of process equipment, and the calculation of process parameters. Through the construction of process knowledge base and logical reasoning mechanism, digital modeling and simulation of raw material processing methods, processing sequence and resource allocation are realized. The traditional degradable tableware material formula optimization method refers to the adjustment of the proportion of raw materials such as polylactic acid or starch-based biodegradable materials in tableware manufacturing. It is usually completed by using the laboratory manual trial and error method. The researchers select polylactic acid particles, starch, and plasticizers according to experience, weigh and mix them, melt and blend the mixture using a twin-screw extruder, and then extrude and granulate. Then, the particles are prepared into standard samples by a hot press forming machine. The mechanical strength and degradation period of the samples are measured in a tensile testing machine and a degradation test box. Finally, the technicians manually correct the weight of each component according to the physical values obtained from the experiment and repeat the weighing, mixing, extruding, and testing steps to determine the production formula.

[0003] The existing technology relies on manual experience to select raw materials and physically mix and try. It needs to go through a tedious extrusion granulation, hot press forming, and destructive physical testing process, resulting in a long formula development cycle and huge material loss. The experimental data can only reflect the macro mechanical results and cannot reveal the failure mechanism at the microscopic molecular level. It is difficult to predict the long-term degradation behavior and structural stability of the material in different environments at the early stage of forming. This causes the formula adjustment to lack theoretical guidance and makes it difficult to achieve precise life control, which cannot meet the strict requirements of industrial production for degradable controllability. SUMMARY

[0004] In order to achieve the above purpose, the application adopts the following technical scheme: a degradable tableware material formula optimization method, comprising the following steps:

[0005] S1: performing energy minimization calculation based on atomic coordinate data and potential energy function parameters to generate an equilibrium state molecular conformation model, discretizing and aggregating to generate interface gap cluster data, calculating the volume of the interface gap cluster data and the chain segment transition energy level potential barrier parameters, and generating initial formula component parameters;

[0006] S2: based on the initial formula component parameters, calculate the storage modulus and loss modulus parameters, convert to generate discrete relaxation time spectrum data, accumulate the component weight of the relaxation period greater than the cooling and curing characteristic time length parameter, generate the theoretical elastic recovery strain parameter;

[0007] S3: map the theoretical elastic recovery strain parameter to the deformation amount calculated by the three-dimensional geometric model of the dinner plate, discretize it into voxel grid nodes, calculate the node-to-boundary distance and mark it as the diffusion characteristic length parameter;

[0008] S4: call the initial formula component parameters, extract the interface void cluster data to calculate the substrate diffusion coefficient parameter, combine the diffusion characteristic length parameter to calculate the physical outflow and chemical generation flux parameter, compare the difference and update the local hydrogen ion concentration index of the voxel grid node;

[0009] S5: update the local hydrogen ion concentration index to the reaction rate constant, calculate the mass loss cumulative time length parameter and the time change gradient parameter, mark the structure failure position coordinates, and integrate to generate the shaped formula structure parameter.

[0010] As a further scheme of the present application, the initial formula component parameters include starch particle volume fraction, polylactic acid matrix mass ratio, interface modifier grafting rate, the theoretical elastic recovery strain parameter includes in-mold shrinkage rate value, size rebound after demolding, warping deformation index, the diffusion characteristic length parameter includes node-to-boundary Euclidean distance, effective diffusion path length, cross-sectional geometric shape factor, the local hydrogen ion concentration index includes cumulative acid product molar amount, local protonation degree index, hydrolysis reaction activity concentration, and the shaped formula structure parameter includes optimal component addition ratio, product macroscopic geometric wall thickness, and reinforcement rib distribution topology data.

[0011] As a further scheme of the present application, the initial formula component parameter acquisition step is specifically:

[0012] S101: based on the atomic coordinate data and potential energy function parameters of starch and polylactic acid molecules, perform energy minimization iteration calculation, monitor the total potential energy convergence state, stop iteration when the total potential energy derivative is less than the preset convergence standard, lock the bond length and bond angle data between atoms, establish a three-dimensional space structure with atomic interaction force in force balance state, and generate a microcosmic stable configuration of the composite system;

[0013] S102: discretize the microcosmic stable configuration of the composite system into microcosmic grid units, set a virtual probe radius parameter, traverse the microcosmic grid units to perform space occupation detection, mark the blank grid units not covered by the atomic van der Waals radius, aggregate adjacent connected blank grid units to form independent irregular cavity structures, extract the cavity structure geometric boundary coordinates, and establish a microcosmic free volume topology set;

[0014] S103: Call the microscopic free volume topology set, calculate the cumulative free volume fraction and average pore size of the cavity structure and mark them as volume parameters, calculate the work done by polymer chain segments to overcome intermolecular forces when crossing the cavity structure and mark them as chain segment transition energy level barrier parameters, and match the component concentration configuration based on the volume parameters and chain segment transition energy level barrier parameters to generate initial formulation component parameters.

[0015] As a further aspect of the present invention, the process for obtaining the preset convergence criterion is as follows:

[0016] The potential energy function parameters are called to resolve the bond stretching constant and filter out the stiffness coefficient with the largest value. The atomic coordinate data of starch and polylactic acid molecules are traversed to calculate the Euclidean distance between covalently connected atomic pairs. The minimum bond length feature value is extracted, and a fixed displacement tolerance ratio constant is obtained. The product of the minimum bond length feature value and the displacement tolerance ratio constant is calculated and marked as the maximum allowable residual displacement. The multiplication operation of the maximum stiffness coefficient and the maximum allowable residual displacement is performed to calculate the critical force amplitude for maintaining atomic dynamic equilibrium. The floating-point machine precision of the current computing environment is detected and converted into the minimum energy gradient limit. The values ​​of the critical force amplitude and the minimum energy gradient limit are compared, and the maximum value of the two is selected and set as the preset convergence criterion.

[0017] As a further aspect of the present invention, the step of obtaining the theoretical elastic recovery strain parameter specifically includes:

[0018] S201: Based on the initial formula component parameters, construct a virtual melt model, apply broadband sinusoidal shear boundary conditions, calculate stress response waveform data, decompose it into in-phase elastic component data and out-of-phase viscous component data through Fourier transform, calculate energy storage modulus parameters and loss modulus parameters, and generate dynamic rheological response modulus data.

[0019] S202: Call the dynamic rheological response modulus data, establish a discretized integral kernel function matrix that correlates the frequency modulus and relaxation intensity, perform iterative error minimization operation to invert the relaxation intensity distribution, determine the modulus contribution weight at different time scales, and generate discrete relaxation time spectrum data;

[0020] S203: Obtain the characteristic cooling and curing time parameters of the injection molding process, traverse the discrete relaxation time spectrum data, identify slow relaxation units whose time constant is greater than the cooling and curing time value, accumulate the modulus contribution weights corresponding to the slow relaxation units, calculate the residual elastic deformation potential energy, and generate theoretical elastic recovery strain parameters.

[0021] As a further aspect of the present invention, the step of obtaining the diffusion feature length parameter specifically includes:

[0022] S301: Call the theoretical elastic recovery strain parameters and map them point by point to the finite element mesh nodes of the three-dimensional geometric model of the plate. Construct a set of mechanical equilibrium equations containing geometric nonlinearity, perform numerical iteration to solve the three-dimensional deformation displacement vector of each node under residual stress, superimpose the displacement vector to the original coordinates to reconstruct the mesh shape, and generate macroscopic deformation geometric topology data.

[0023] S302: Based on the macroscopic deformation geometric topology data, define the closed entity space, set the isotropic discretization step size, perform a global voxelization scan operation, convert the continuous geometric entity into a discrete voxel mesh array, remove the external background mesh and extract the three-dimensional Cartesian coordinates of the geometric center of the internal voxel unit, and establish a voxelized node space coordinate set.

[0024] S303: Call the voxelized node spatial coordinate set, construct a spatial neighborhood index structure to traverse the outer surface mesh of the macroscopic deformation geometric topology data, retrieve the nearest neighbor boundary projection point corresponding to each voxel node, calculate the minimum Euclidean straight-line distance from the voxel center to the projection point, mark the distance value as a scalar field attribute, and generate a diffusion feature length parameter.

[0025] As a further aspect of the present invention, the step of obtaining the local hydrogen ion concentration index specifically includes:

[0026] S401: Call the initial formulation component parameters, extract the void volume ratio data based on the interface void cluster data, construct an effective diffusion transport model of solute in porous media, calculate the migration rate of substances in polymer matrix and mark it as substrate diffusion coefficient parameter, combine with diffusion characteristic length parameter, calculate the physical flux value of acidic substances migrating from the interior to the surface per unit time, and generate physical efflux flux parameter of acidic products.

[0027] S402: Call the initial formulation component parameters, extract the molar concentration value of hydrolyzable ester bonds in the polylactic acid molecular chain, obtain the preset hydrolysis reaction rate constant, perform multiplication operation to calculate the chemical reaction rate of ester bond breaking to generate carboxyl terminus, quantify the amount of acidic terminus group material generated per unit volume per unit time, and generate acidic terminus group chemical generation flux parameter.

[0028] S403: Calculate the numerical difference between the chemical generation flux parameter of the acidic terminal group and the physical efflux flux parameter of the acidic product, determine the net accumulation rate of acidic substances in the local area, add the difference data to the local acidic substance concentration parameter at the current time step, convert the molar concentration into a hydrogen ion activity value using the acid dissociation equilibrium constant, update the acidity state of the voxel grid node, and generate a local hydrogen ion concentration index.

[0029] As a further aspect of the present invention, the process for obtaining the preset hydrolysis reaction rate constant is specifically as follows:

[0030] The initial formulation component parameters are called to analyze the ester bond chemical structure type of polylactic acid molecules. The standard ester bond hydrolysis activation energy parameter and pre-exponential factor physical quantity corresponding to the target chemical structure type are matched. The absolute temperature value set in the current simulation scenario is obtained. The standard molar gas constant is called to calculate the product of the absolute temperature value and the standard molar gas constant to generate a thermal noise energy benchmark. The division operation between the standard ester bond hydrolysis activation energy parameter and the thermal noise energy benchmark is performed. The result of the operation is inverted and the natural exponential function value is calculated to obtain the effective collision probability of molecules. The multiplication operation between the pre-exponential factor physical quantity and the effective collision probability of molecules is performed to generate the intrinsic rate value of chemical bond breaking. The intrinsic rate value is established as the preset hydrolysis reaction rate constant.

[0031] As a further aspect of the present invention, the step of obtaining the structural parameters of the finalized formula specifically includes:

[0032] S501: Call the local hydrogen ion concentration index to simulate the evolution of the autocatalytic reaction rate. Use numerical iteration to solve and dynamically update the hydrolysis reaction rate constant at each time step to simulate the mass loss during the material degradation process. When the remaining mass fraction is lower than the preset functional threshold, record the corresponding time span and mark it as the cumulative mass loss duration parameter to generate hydrolysis reaction kinetic evolution data.

[0033] S502: Call the hydrolysis reaction kinetic evolution data, extract the time evolution curve of the hydrolysis reaction rate constant, calculate the slope change of adjacent time steps as the time change gradient parameter, detect the gradient value of each voxel grid node, identify the location where the gradient undergoes a step change as a potential fracture point, extract the corresponding three-dimensional spatial geometric coordinates, and generate a material structure failure topology set.

[0034] S503: By comprehensively analyzing the hydrolysis reaction kinetic evolution data and the material structure failure topology set, a feature dataset integrating the degradation lifetime in the time dimension and the failure distribution in the spatial dimension is constructed. The dataset is then mapped back to the original formulation component space as a constraint condition for formulation screening. Combined with the degradation cycle and structural integrity requirements of the product design, the component ratio and geometric shape are screened to generate the finalized formulation structure parameters.

[0035] A biodegradable plate material formulation optimization system, including:

[0036] The molecular conformation optimization module performs energy minimization calculations based on atomic coordinate data and potential energy function parameters, generates an equilibrium molecular conformation model, discretizes and aggregates interface void cluster data, calculates the volume of the interface void cluster data and the potential barrier parameters of the chain segment transition energy level, and generates the initial formulation component parameters.

[0037] The elastic strain calculation module calculates the energy storage and loss modulus parameters based on the initial formula component parameters, converts them into discrete relaxation time spectrum data, accumulates the component weights of relaxation periods greater than the cooling and solidification characteristic duration parameters, and generates theoretical elastic recovery strain parameters.

[0038] The geometric discretization module maps the theoretical elastic recovery strain parameters to the three-dimensional geometric model of the plate to calculate the deformation, discretizes it into voxel mesh nodes, calculates the distance from the node to the boundary and marks it as the diffusion feature length parameter.

[0039] The concentration dynamic update module calls the initial formulation component parameters, extracts the interface void cluster data to calculate the substrate diffusion coefficient parameters, combines the diffusion characteristic length parameters to calculate the physical expulsion and chemical generation flux parameters, compares the difference and updates the local hydrogen ion concentration index of the voxel grid nodes.

[0040] The structural failure integration module updates the reaction rate constant by the local hydrogen ion concentration index, calculates the cumulative mass loss duration parameter and time change gradient parameter, marks the coordinates of the structural failure location, and integrates to generate the finalized formula structural parameters.

[0041] Compared with the prior art, the advantages and positive effects of the present invention are as follows:

[0042] In this invention, a microscopic molecular dynamics model is constructed to analyze the interfacial void clusters and chain segment transition energy barriers, enabling quantitative assessment of the material's microscopic compatibility and precise initial formulation design. The elastic recovery strain during the molding process is inverted using relaxation time spectra to predict geometric deformation risks. Combining autocatalytic reaction kinetics and porous media diffusion theory, the evolution of local hydrogen ion concentration and the migration path of acidic products are tracked in real time, establishing a dynamic mapping relationship between mass loss and structural mechanical failure. Full life cycle performance verification is completed in a virtual environment, shortening the R&D cycle and reducing trial and error costs, ensuring that the product has excellent structural stability while meeting degradation requirements. Attached Figure Description

[0043] To more clearly illustrate the technical solutions in the embodiments of the present invention, the accompanying drawings used in the description of the embodiments will be briefly introduced below. Obviously, the accompanying drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0044] Figure 1 This is a schematic diagram of the steps of the present invention;

[0045] Figure 2 This is a detailed schematic diagram of S1 of the present invention;

[0046] Figure 3This is a detailed schematic diagram of S2 of the present invention;

[0047] Figure 4 This is a detailed schematic diagram of S3 of the present invention;

[0048] Figure 5 This is a detailed schematic diagram of S4 of the present invention;

[0049] Figure 6 This is a detailed schematic diagram of S5 of the present invention;

[0050] Figure 7 This is a system module diagram of the present invention. Detailed Implementation

[0051] The technical solution of the present invention will now be described with reference to the accompanying drawings.

[0052] To make the technical problems, technical solutions and advantages of the present invention clearer, a detailed description will be given below in conjunction with the accompanying drawings and specific embodiments.

[0053] Please see Figure 1 This invention provides a method for optimizing the formulation of biodegradable plate materials, including the following steps:

[0054] S1: Perform energy minimization calculations based on atomic coordinate data and potential energy function parameters to generate an equilibrium molecular conformation model. Discretize and aggregate the interface void cluster data to generate interface void cluster data and calculate the volume and chain segment transition energy level barrier parameters of the interface void cluster data to generate initial formulation component parameters.

[0055] S2: Based on the initial formulation component parameters, calculate the energy storage and loss modulus parameters, convert them into discrete relaxation time spectrum data, accumulate the component weights of relaxation periods greater than the cooling and solidification characteristic duration parameters, and generate theoretical elastic recovery strain parameters.

[0056] S3: Map the theoretical elastic recovery strain parameters to the three-dimensional geometric model of the plate to calculate the deformation, discretize it into voxel mesh nodes, calculate the distance from the node to the boundary and mark it as the diffusion feature length parameter;

[0057] S4: Call the initial formulation component parameters, extract the interface void cluster data to calculate the substrate diffusion coefficient parameters, combine the diffusion characteristic length parameters to calculate the physical expulsion and chemical generation flux parameters, compare the difference and update the local hydrogen ion concentration index of the voxel grid nodes.

[0058] S5: Update the reaction rate constant with the local hydrogen ion concentration index, calculate the cumulative mass loss duration parameter and time change gradient parameter, mark the coordinates of the structural failure location, and integrate to generate the finalized formula structural parameters.

[0059] The initial formulation component parameters include starch granule volume fraction, polylactic acid matrix mass ratio, and interface modifier grafting rate. The theoretical elastic recovery strain parameters include in-mold shrinkage rate, post-demolding dimensional rebound, and warpage index. The diffusion characteristic length parameters include node-to-boundary Euclidean distance, effective diffusion path length, and cross-sectional geometry factor. The local hydrogen ion concentration index includes cumulative acidic product molar amount, local protonation index, and hydrolysis reaction activity concentration. The final formulation structural parameters include optimal component addition ratio, product macroscopic geometric wall thickness, and reinforcing rib distribution topology data.

[0060] Please see Figure 2 The specific steps for obtaining the initial formulation component parameters are as follows:

[0061] S101: Based on the atomic coordinate data and potential energy function parameters of starch and polylactic acid molecules, perform energy minimization iterative calculation, monitor the convergence state of the total potential energy, stop iteration when the derivative of the total potential energy is less than the preset convergence standard, lock the bond length and bond angle data between atoms, establish a three-dimensional spatial structure in which the atomic interaction force is in a state of force equilibrium, and generate the microscopic steady-state configuration of the composite system.

[0062] The system invokes a pre-built molecular dynamics simulation engine, loads the initial atomic coordinate file containing starch polysaccharide segments and polylactic acid (PLA) segments, and matches it with general polymer potential function parameters such as COMPASSII or PCFF. It then analyzes the defined bond stretching, angular bending, dihedral torsion, and non-bonded interaction force field parameters. During this process, the convergence logic for the energy minimization iterative calculation follows a specific pre-defined convergence criterion acquisition procedure. Specifically, it calls the potential function parameters, analyzes the bond stretching constant matrix, iterates through all defined chemical bond types, and selects the stiffness coefficient with the largest value. For example, the stretching vibration stiffness coefficient of the carbon-oxygen double bond in the PLA backbone is the maximum value in the defined system, and its value is set to 5000 kcal / mol / square Å. Subsequently, it iterates through the atomic coordinate data of starch and PLA molecules, calculates the Euclidean distances between all covalently bonded atom pairs, and extracts the minimum bond length eigenvalue. For example, the minimum bond length measured at the hydrogen-carbon atom connection is 1.09 Å. A preset, fixed displacement tolerance proportionality constant is obtained, typically set to 1 x 10⁻⁵, to limit the perturbation level of atoms near their equilibrium positions. Next, the product of the minimum bond length eigenvalue and the displacement tolerance proportionality constant is calculated, i.e., 1.09 multiplied by 0.00001, resulting in 0.0000109 Å. This result is marked as the maximum permissible residual displacement. The process continues by multiplying the maximum stiffness coefficient by the maximum permissible residual displacement, i.e., 5000 multiplied by 0.0000109, to calculate the critical force amplitude required to maintain atomic dynamic equilibrium, which is 0.0545 kcal / mol / Å. Simultaneously, the floating-point machine precision of the current computing environment is checked; for example, 1 x 10⁻¹⁵ in a double-precision floating-point environment, it is converted to the physical dimensions of the energy gradient as the minimum energy gradient limit. Finally, the critical force amplitude of 0.0545 is compared with the minimum energy gradient threshold, and the maximum value of 0.0545 is selected as the preset convergence criterion. During the iteration process, the total potential energy derivative is monitored in real time. When the calculated root mean square force gradient is less than the preset convergence criterion of 0.0545, the system is considered to have reached thermodynamic equilibrium, the iteration stops, the interatomic bond lengths and bond angles are locked, a three-dimensional spatial structure with atomic interaction forces in equilibrium is established, and the final microscopic steady-state configuration of the composite system is generated.

[0063] S102: Discretize the microscopic steady-state configuration of the composite system into microscopic grid units, set the virtual probe radius parameter, traverse the microscopic grid units to perform space occupancy detection, mark the blank grid units not covered by the atomic van der Waals radius, aggregate adjacent connected blank grid units to form independent irregular cavity structures, extract the geometric boundary coordinates of the cavity structures, and establish a microscopic free volume topology set.

[0064] The microscopic steady-state configuration of the composite system is discretized into microscopic grid cells. The discretization operation first defines a three-dimensional grid matrix covering the entire molecular configuration's three-dimensional bounding box, setting the side length resolution of the microscopic grid cells, for example, 0.5 Å, thus dividing the continuous spatial domain into discrete sets of voxels. Subsequently, a virtual probe radius parameter is set, based on the van der Waals radius of a helium atom or a specific gas molecule, for example, a value of 1.2 Å. The scanning program traverses each microscopic grid cell to perform space occupancy detection. For any grid center point in the i-th row, j-th column, k-th layer of the matrix, the Euclidean distance between it and the nuclear coordinates of its nearest neighbor atoms is calculated. If this distance is less than the sum of the van der Waals radius of the corresponding atom and the virtual probe radius, the grid is determined to be occupied by an atomic entity; otherwise, if the distance is greater than the sum of the aforementioned radii, the position is marked as a blank grid cell not covered by the atomic van der Waals radius. After completing the full-domain scan, a connected component analysis algorithm is executed to aggregate adjacent connected blank grid cells, merging all blank voxels sharing faces or edges to form independent irregular cavity structures. The geometric boundary coordinates of these cavity structures are extracted, and the set of vertices on the closed surface of each cavity is identified, establishing a microscopic free volume topology set. For example, in a system containing 1 million grid cells, 150,000 blank grid cells are detected and labeled, forming 350 independent void clusters through aggregation. Each cluster consists of hundreds to thousands of consecutive blank voxels. The spatial distribution data of these clusters constitutes the microscopic free volume topology set, providing the geometric basis for subsequent diffusion performance prediction.

[0065] S103: Call the microscopic free volume topology set, calculate the cumulative free volume fraction and average pore size of the cavity structure and mark them as volume parameters, calculate the work done by polymer chain segments to overcome intermolecular forces when crossing the cavity structure and mark them as chain segment transition energy level barrier parameters, match the component concentration configuration based on the volume parameters and chain segment transition energy level barrier parameters, and generate the initial formulation component parameters.

[0066] Using a microscopic free volume topology set, the number of blank micro-grid units contained in each independent irregular cavity structure is first counted. This number is then multiplied by the volume of a single grid unit to calculate the cumulative free volume fraction of the cavity structure. For example, if the volume of a single grid unit is 0.125 cubic angstroms, and a cavity contains 800 grid units, then the cavity volume is 100 cubic angstroms. The total volume of all cavities is summed and divided by the total volume of the simulation system to obtain the free volume ratio, such as 10.125. Simultaneously, the average pore size is calculated based on the cavity geometry and, together with the volume fraction, is labeled as a volume parameter. Next, the work done by polymer chain segments to overcome intermolecular forces as they traverse the cavity structure is calculated. The energy increment required for a chain segment to transition from a high-density region to an adjacent cavity is calculated using the Lennard-Jones potential energy formula, and this increment is labeled as the chain segment transition energy level barrier parameter. For example, the calculated average barrier for polylactic acid chain segment transitions is 45 kilojoules per mole. Finally, based on the volume parameters and chain segment transition energy level barrier parameters, the component concentration configuration is matched, and the starch particle volume fraction and polylactic acid matrix mass ratio are adjusted according to a preset structure-activity relationship database to generate initial formulation component parameters. For example, when the free volume ratio is detected to be less than 8%, the grafting rate of the interface modifier is automatically increased to 2.5 to optimize the interface bonding, as shown in Table 1. Table 1 lists examples of initial formulation component parameters derived from microstructure parameters.

[0067] Table 1 Example of initial formulation component parameters

[0068]

[0069] As shown in Table 1, this step determines the specific ratio of starch, polylactic acid and modifier based on the calculated microstructure characteristics, which serves as the input benchmark for subsequent macroscopic performance simulation.

[0070] Please see Figure 3 The specific steps for obtaining the theoretical elastic recovery strain parameters are as follows:

[0071] S201: Based on the initial formulation component parameters, a virtual melt model is constructed, broadband sinusoidal shear boundary conditions are applied, stress response waveform data is calculated, and Fourier transform is used to decompose it into in-phase elastic component data and out-of-phase viscous component data. The energy storage modulus parameter and loss modulus parameter are calculated, and dynamic rheological response modulus data is generated.

[0072] Based on the initial formulation component parameters, a virtual melt model was constructed, and a non-equilibrium molecular dynamics simulation box containing polydisperse polymer chains was established, with periodic boundary conditions set. Subsequently, a broadband sinusoidal shear boundary condition was applied, imposing a shear deformation of the form strain equal to the maximum amplitude multiplied by a sine function (angular frequency multiplied by time) on the upper surface of the simulation box, covering a wide frequency range from 0.1 radians per second to 100 radians per second. The stress response waveform data generated during the simulation was recorded in real time and decomposed into in-phase elastic component data and out-of-phase viscous component data through Fourier transform. Specifically, the amplitude of the component in the stress waveform that is in phase with the strain waveform was extracted, divided by the strain amplitude, and the storage modulus parameter (G') was calculated; the amplitude of the component with a 90-degree phase difference from the strain waveform was extracted, and the loss modulus parameter (G'') was calculated, thereby generating dynamic rheological response modulus data. For example, at a shear frequency of 10 radians per second and an applied strain amplitude of 0.05, the measured stress response amplitude was 5 Pascals, with a phase difference of 0.6 radians. Decomposition calculations yielded a storage modulus of 82.53 Pascals and a loss modulus of 56.46 Pascals. By iterating through all test frequency points, a series of discrete data points corresponding to frequencies and moduli were generated, providing fundamental data for constructing a viscoelastic property model of the material.

[0073] S202: Call the dynamic rheological response modulus data, establish the discretized integral kernel function matrix of the correlation frequency modulus and relaxation intensity, perform iterative error minimization operation to invert the relaxation intensity distribution, determine the modulus contribution weight at different time scales, and generate discrete relaxation time spectrum data;

[0074] A discretized integral kernel function matrix relating the frequency modulus and relaxation intensity is established. This process is based on the generalized Maxwell model, expressing the storage modulus as the integral of the square of the product of relaxation intensity and frequency plus relaxation time, divided by the square of the product of frequency and relaxation time. An iterative error minimization operation is performed to invert the relaxation intensity distribution. Using nonlinear least squares or Tikhonov regularization algorithms, the relaxation modulus amplitude at preset relaxation time nodes is solved, and the modulus contribution weights at different time scales are determined, generating discrete relaxation time spectrum data. For example, setting a logarithmically equally spaced relaxation time series from 0.01 seconds to 1000 seconds, inversion calculations reveal a significant peak at a relaxation time of 1 second, corresponding to a relaxation modulus intensity of 1500 Pascals, while the intensity at a relaxation time of 100 seconds is 200 Pascals. These data points constitute the time-response characteristic spectrum of different molecular motion modes within the material, accurately describing the dynamic process of the material's transition from a solid-like to a liquid-like state.

[0075] S203: Obtain the characteristic duration parameters of the cooling and solidification of the injection molding process, traverse the discrete relaxation time spectrum data, identify the slow relaxation units whose time constant is greater than the value of the cooling and solidification duration, accumulate the modulus contribution weights corresponding to the slow relaxation units, calculate the residual elastic deformation potential energy, and generate theoretical elastic recovery strain parameters.

[0076] The cooling and solidification characteristic duration parameter of the injection molding process is obtained. This parameter is determined by the actual injection molding cycle, for example, set to 20 seconds, representing the time it takes for the material to cool from the molten state to below the glass transition temperature within the mold. The discrete relaxation time spectrum data is traversed to identify slow relaxation units with time constants greater than the cooling and solidification duration, i.e., all spectral components with relaxation time t greater than 20 seconds are selected. Subsequently, the modulus contribution weights corresponding to the slow relaxation units are accumulated, and the modulus intensity values ​​corresponding to these long relaxation times are summed to calculate the residual elastic deformation potential energy corresponding to the stress that was not completely released at the moment of cooling completion. Assuming that there are three components in the relaxation spectrum with time constants of 50 seconds, 100 seconds, and 500 seconds, their corresponding modulus contributions are 100 Pascals, 50 Pascals, and 10 Pascals, respectively, an accumulation calculation is performed to obtain a total effective modulus contribution of 160 Pascals. Combined with the applied rheological strain level, the theoretical elastic recovery strain parameter is generated. For example, the calculated theoretical elastic recovery strain is 0.035, indicating that the material has a potential dimensional springback tendency of 3.5% after demolding. This parameter directly quantifies the influence of microscopic viscoelastic properties on macroscopic dimensional stability.

[0077] Please see Figure 4 The specific steps for obtaining the diffusion feature length parameter are as follows:

[0078] S301: Call the theoretical elastic recovery strain parameters, map them point by point to the finite element mesh nodes of the three-dimensional geometric model of the plate, construct a set of mechanical equilibrium equations containing geometric nonlinearity, perform numerical iteration to solve the three-dimensional deformation displacement vector of each node under residual stress, superimpose the displacement vector to the original coordinates to reconstruct the mesh shape, and generate macroscopic deformation geometric topology data.

[0079] The theoretical elastic recovery strain parameters are called and mapped point by point to the finite element mesh nodes of the 3D geometric model of the plate. The CAD model of the plate is imported and tetrahedral meshed to generate a finite element model containing tens of thousands of nodes. The theoretical elastic recovery strain of 0.035 calculated above is used as a pre-strain tensor and loaded onto each mesh node. A set of mechanical equilibrium equations with geometric nonlinearity is constructed, considering the large deformation characteristics of the material, and numerical iteration is performed to solve the 3D deformation displacement vector of each node under residual stress. The Newton-Raphson iterative method is used to solve the equilibrium equations until the residual force at the node is less than the preset tolerance. After the solution is completed, the displacement vector is superimposed to the original coordinates to reconstruct the mesh shape. For example, if the original coordinates of an edge node are (100, 0, 5) and the calculated displacement vector is (-0.5, 0, 0.2), then the reconstructed coordinates are (99.5, 0, 5.2). All nodes are traversed to complete the coordinate update and generate macroscopic deformation geometric topology data, which accurately describes the actual warping and shrinkage of the plate after cooling due to the release of internal stress.

[0080] S302: Based on macroscopic deformation geometric topology data, define the closed entity space, set the isotropic discretization step size, perform global voxelization scanning operation, convert the continuous geometric entity into a discrete voxel mesh array, remove the external background mesh and extract the three-dimensional Cartesian coordinates of the geometric center of the internal voxel unit, and establish a voxelized node space coordinate set.

[0081] Based on macroscopic deformation geometric topology data, a closed solid space is defined, and the external boundary surface of the deformed mesh is identified to confirm that it constitutes a watertight closed volume. An isotropic discretization step size is set, for example, a voxel edge length of 1 mm. A global voxelization scan operation is performed, using ray casting or parity checking rules to convert the continuous geometric solid into a discrete voxel mesh array. During this process, the external background mesh is removed, retaining only the voxel elements located inside the plate solid, and the three-dimensional Cartesian coordinates of the geometric centers of the internal voxel elements are extracted to establish a voxelized node spatial coordinate set. For example, for a plate model with a diameter of 200 mm and a thickness of 2 mm, after voxelization, approximately 60,000 effective voxel nodes are generated. Each node records its absolute position coordinates (x, y, z) in space, and these coordinate data constitute the spatial carrier for subsequent chemical degradation simulations.

[0082] S303: Call the voxelized node spatial coordinate set, construct a spatial neighborhood index structure to traverse the outer surface mesh of the macroscopic deformation geometric topology data, retrieve the nearest neighbor boundary projection point corresponding to each voxel node, calculate the minimum Euclidean straight distance from the voxel center to the projection point, mark the distance value as a scalar field attribute, and generate the diffusion feature length parameter.

[0083] The voxelized node spatial coordinate set is invoked to construct a spatial neighborhood index structure, such as a KD-Tree or Octree structure, to accelerate nearest neighbor search. The external surface mesh of the macroscopic deformation geometry topology data is traversed, and for each voxel node in the coordinate set, its nearest neighbor boundary projection point on the surface mesh is retrieved. The minimum Euclidean distance from the voxel center to the projection point is calculated. For example, for a voxel node with center coordinates (50, 50, 1), the nearest surface point coordinates are (50, 50, 2), resulting in a calculated distance of 1 mm. This distance value is labeled as a scalar field property and assigned to the voxel node. After completing the global calculation, a diffusion characteristic length parameter is generated. This parameter intuitively reflects the depth of each material element from the external environment; a smaller distance indicates closer proximity to the surface, making it easier for acidic products to be expelled, while a larger distance indicates deeper layers, making it easier for acidic substances to accumulate.

[0084] Please see Figure 5 The specific steps for obtaining the local hydrogen ion concentration index are as follows:

[0085] S401: Call the initial formulation component parameters, extract the void volume ratio data based on the interface void cluster data, construct an effective diffusion transport model of solute in porous media, calculate the migration rate of substances in polymer matrix and mark it as the substrate diffusion coefficient parameter, combine the diffusion characteristic length parameter, calculate the physical flux value of acidic substances migrating from the interior to the surface per unit time, and generate the physical efflux flux parameter of acidic products.

[0086] The initial formulation component parameters are called, and the void volume ratio data, for example, 5%, is extracted based on the interface void cluster data. An effective diffusion transport model of the solute in the porous medium is constructed using the Maxwell-Eucken model or the Bruggeman approximation. The diffusion coefficient of the pure substrate (e.g., 1 x 10⁻¹⁰ m² / s) and porosity are input, and the migration rate of the substance in the polymer matrix is ​​calculated and labeled as the substrate diffusion coefficient parameter. It is assumed that the effective diffusion coefficient is reduced to 0.8 x 10⁻¹⁰ m² / s based on the model calculation. Combining the diffusion characteristic length parameter, for example, a node with a characteristic length of 1 mm, the physical flux of acidic substances migrating from the interior to the surface per unit time is calculated. According to Fick's first law, the flux equals the diffusion coefficient multiplied by the concentration gradient. If the local concentration gradient is 10 mol / m³ / m, the calculated flux is 0.8 x 10⁻⁹ mol / m² / s, generating the physical efflux flux parameter for acidic products.

[0087] S402: Call the initial formula component parameters, extract the molar concentration of hydrolyzable ester bonds in the polylactic acid molecular chain, obtain the preset hydrolysis reaction rate constant, perform multiplication operation to calculate the chemical reaction rate of ester bond breaking to generate carboxyl terminus, quantify the amount of acidic terminal group material generated per unit volume per unit time, and generate acidic terminal group chemical generation flux parameters.

[0088] The initial formulation component parameters are used to extract the molar concentration of hydrolyzable ester bonds in the polylactic acid (PLA) molecular chain, for example, an initial concentration of 15 mol / L. A preset hydrolysis rate constant is obtained, the process of which involves rigorous physicochemical derivation: First, the initial formulation component parameters are used to analyze the ester bond chemical structure type of the PLA molecule (e.g., α-hydroxy ester), and the corresponding standard ester bond hydrolysis activation energy parameter (e.g., 80,000 joules per mole) and pre-exponential factor physical quantity (e.g., 1 x 10⁸ per second). The absolute temperature value set in the current simulation scenario is obtained, for example, 333.15 Kelvin (60 degrees Celsius), and the standard molar gas constant 8.314 joules per mole per Kelvin is used. The product of the absolute temperature value and the standard molar gas constant, i.e., 333.15 multiplied by 8.314, is calculated to generate a thermal noise energy benchmark, approximately 2769.8 joules per mole. The standard ester bond hydrolysis activation energy parameter of 80000 was divided by the thermal noise energy reference of 2769.8, yielding a result of approximately 28.88. This result was then inverted (-28.88), and the natural exponential function value (e to the power of -28.88) was calculated to obtain the effective collision probability of the molecule, which is approximately 2.86 x 10^-13. Finally, the pre-exponential factor physical quantity of 1 x 10^8 was multiplied by the effective collision probability of the molecule to generate the intrinsic rate of chemical bond breaking, which is approximately 2.86 x 10^-5 per second. This value was established as the preset hydrolysis reaction rate constant. Based on this constant, the chemical reaction rate for the formation of carboxyl-terminal groups from ester bond breaking was calculated using multiplication operations. The amount of acidic terminal groups produced per unit volume per unit time was quantified, generating the chemical flux parameters for the formation of acidic terminal groups, as shown in Table 2.

[0089] Table 2 Example of Calculation of Kinetic Parameters for Hydrolysis Reaction

[0090]

[0091] As shown in Table 2, the reaction rate constant at a specific temperature was determined through a rigorous thermodynamic calculation process, and then the theoretical flux of acidic substances was obtained by combining the concentration data.

[0092] S403: Calculate the numerical difference between the chemical generation flux parameter of acidic terminal groups and the physical efflux flux parameter of acidic products, determine the net accumulation rate of acidic substances in the local area, accumulate the difference data to the local acidic substance concentration parameter at the current time step, convert the molar concentration into hydrogen ion activity value using the acid dissociation equilibrium constant, update the acidity state of the voxel grid nodes, and generate a local hydrogen ion concentration index.

[0093] First, the data dimensions are standardized, converting the physical efflux flux into a volumetric degradation rate. If the previously calculated efflux flux is 0.8 x 10⁻⁹ mol / m² / s, and considering the voxel unit characteristic length of 1 mm, the concentration decrease rate caused by efflux is 0.8 x 10⁻⁶ mol / L / s. This is then subtracted from the chemical generation flux in Table 2 (4.29 x 10⁻⁴ mol / L / s), resulting in a difference of approximately 4.28 x 10⁻⁴ mol / L / s. This indicates that the net accumulation rate of acidic substances in the local area is positive and relatively large. This difference is then added to the local acidic substance concentration parameter at the current time step. If the concentration at the previous time step was 0.1 mol / L and the time step length was 10 seconds, the updated concentration increases by 0.00428 mol / L. Subsequently, the molar concentration was converted into a hydrogen ion activity value using the acid dissociation equilibrium constant (Ka, e.g., 1.3 x 10⁻⁴). According to the weak acid ionization formula, the hydrogen ion concentration is approximately equal to the square root of Ka multiplied by the acid concentration. Substituting this into the calculation, the hydrogen ion concentration is approximately 0.0037 mol / L. The acidity state of the voxel grid nodes is then updated to generate a local hydrogen ion concentration index. This process dynamically simulates the "internal acidic autocatalysis" phenomenon caused by diffusion hysteresis within the material.

[0094] Please see Figure 6 The specific steps for obtaining the structural parameters of the finalized formula are as follows:

[0095] S501: Call the local hydrogen ion concentration index to simulate the evolution of the autocatalytic reaction rate. Use numerical iteration to solve and dynamically update the hydrolysis reaction rate constant at each time step to simulate the mass loss during the material degradation process. When the remaining mass fraction is lower than the preset functional threshold, record the corresponding time span and mark it as the cumulative mass loss duration parameter to generate hydrolysis reaction kinetic evolution data.

[0096] By substituting the local hydrogen ion concentration index into the autocatalytic reaction rate evolution process, and employing the autocatalytic kinetic equation, the new reaction rate constant is equal to the base rate constant plus the product of the catalytic coefficient and the hydrogen ion concentration. For example, if the base constant is 2.86 x 10⁻⁵, the catalytic coefficient is 0.5, and the hydrogen ion concentration is 0.0037, then the updated rate constant is 1.88 x 10⁻³, significantly improving the rate. The hydrolysis reaction rate constant is dynamically updated at each time step using numerical iteration to simulate the mass loss during material degradation. The decrease in polymer molecular weight within each time step is calculated and converted into mass loss. When the remaining mass fraction falls below a preset functional threshold (e.g., 80% of the initial mass), the total time elapsed from the start of the simulation to this point is recorded, for example, 1200 hours, and marked as the cumulative mass loss duration parameter, generating hydrolysis reaction kinetic evolution data. This data chain comprehensively records the time history of the material from intact to functional failure.

[0097] S502: Call the hydrolysis reaction kinetic evolution data, extract the time evolution curve of the hydrolysis reaction rate constant, calculate the slope change of adjacent time steps as the time change gradient parameter, detect the gradient value of each voxel grid node, identify the location where the gradient abruptly changes as the potential fracture point, extract the corresponding three-dimensional spatial geometric coordinates, and generate the material structure failure topology set.

[0098] Using hydrolysis reaction kinetics data, the time evolution curve of the hydrolysis reaction rate constant is extracted; this curve typically exhibits an exponential upward trend. The slope change of adjacent time steps is calculated as the time-varying gradient parameter, i.e., the derivative of the rate constant with respect to time is calculated. The gradient values ​​of each voxel grid node are checked individually, with a gradient threshold set, for example, 0.001 per square second. When the gradient value of a node exceeds this threshold, a step abrupt change occurs, indicating a violent self-accelerating degradation reaction at that location, identifying it as a potential fracture point. The corresponding three-dimensional spatial geometric coordinates are extracted, for example (35, 40, 2), generating a material structural failure topology set. This set spatially depicts the distribution of weak regions within the material where structural collapse first occurs.

[0099] S503: By comprehensively analyzing the hydrolysis reaction kinetic evolution data and the material structure failure topology set, a feature dataset is constructed that integrates the degradation lifetime in the time dimension and the failure distribution in the spatial dimension. The dataset is then mapped back to the original formulation component space as a constraint condition for formulation screening. Combined with the degradation cycle and structural integrity requirements of the product design, the component ratio and geometric shape are screened to generate the finalized formulation structural parameters.

[0100] The evolution data and failure topology set are analyzed to construct a feature dataset, which is then mapped back to the original formulation component space as a constraint for formulation screening. Combining the degradation cycle requirements of the product design (e.g., greater than 90 days and less than 180 days) and structural integrity requirements (e.g., no early fracture points in critical stress areas), for example, if the current formulation results in a degradation lifetime of 50 days (not meeting the requirements), the formulation is automatically adjusted by reducing the starch content or increasing the interface modifier, and the above process is iterated again. Finally, when the calculation results meet all constraints, the component ratios (e.g., 20% starch, 78% PLA, 2% modifier) ​​and geometric wall thickness design (e.g., 2.5 mm) are locked, generating the finalized formulation structural parameters, as shown in Table 3.

[0101] Table 3. Results of Screening Structural Parameters for Finalized Formulation

[0102]

[0103] As shown in Table 3, through multiphysics field coupling simulation and iterative optimization, the final determined structural parameters of the finalized formula not only meet the requirements of environmental degradation, but also significantly improve the structural durability of the product, realizing the digital finalization of the formula design.

[0104] Please see Figure 7 A biodegradable tableware material formulation optimization system, including:

[0105] The molecular conformation optimization module performs energy minimization calculations based on atomic coordinate data and potential energy function parameters, generates an equilibrium molecular conformation model, discretizes and aggregates interface void cluster data, calculates the volume of the interface void cluster data and the potential barrier parameters of the chain segment transition energy level, and generates the initial formulation component parameters.

[0106] The elastic strain calculation module calculates energy storage and loss modulus parameters based on the initial formula component parameters, converts them into discrete relaxation time spectrum data, accumulates the component weights of relaxation periods greater than the cooling and solidification characteristic duration parameters, and generates theoretical elastic recovery strain parameters.

[0107] The geometric discretization module maps the theoretical elastic recovery strain parameters to the three-dimensional geometric model of the plate to calculate the deformation, discretizes it into voxel mesh nodes, calculates the distance from the node to the boundary and marks it as the diffusion feature length parameter.

[0108] The concentration dynamic update module calls the initial formulation component parameters, extracts interface void cluster data to calculate the substrate diffusion coefficient parameters, combines the diffusion characteristic length parameters to calculate the physical expulsion and chemical generation flux parameters, compares the difference and updates the local hydrogen ion concentration index of the voxel grid nodes.

[0109] The structural failure integration module updates the reaction rate constant with the local hydrogen ion concentration index, calculates the cumulative mass loss duration parameter and time change gradient parameter, marks the coordinates of the structural failure location, and integrates to generate the finalized formula structural parameters.

[0110] The above description is merely a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the technical scope disclosed in the present invention should be included within the scope of protection of the present invention. Therefore, the scope of protection of the present invention should be determined by the scope of the claims.

Claims

1. A method for optimizing the formulation of biodegradable tableware materials, characterized in that, Includes the following steps: S1: Perform energy minimization calculations based on atomic coordinate data and potential energy function parameters to generate an equilibrium molecular conformation model. Discretize and aggregate the interface void cluster data to generate interface void cluster data and calculate the volume and chain segment transition energy level barrier parameters of the interface void cluster data to generate initial formulation component parameters. S2: Based on the initial formulation component parameters, calculate the energy storage and loss modulus parameters, convert them into discrete relaxation time spectrum data, accumulate the component weights of relaxation periods greater than the cooling and solidification characteristic duration parameters, and generate theoretical elastic recovery strain parameters. S3: Map the theoretical elastic recovery strain parameters to the three-dimensional geometric model of the plate to calculate the deformation, discretize it into voxel mesh nodes, calculate the distance from the node to the boundary and mark it as the diffusion feature length parameter; S4: Call the initial formulation component parameters, extract the interface void cluster data to calculate the substrate diffusion coefficient parameters, combine the diffusion characteristic length parameters to calculate the physical exhaust and chemical generation flux parameters, compare the difference and update the local hydrogen ion concentration index of the voxel grid nodes; S5: Update the reaction rate constant with the local hydrogen ion concentration index, calculate the cumulative mass loss duration parameter and time change gradient parameter, mark the coordinates of the structural failure location, and integrate to generate the finalized formula structural parameters.

2. The method for optimizing the formulation of biodegradable plate materials according to claim 1, characterized in that, The initial formulation component parameters include starch particle volume fraction, polylactic acid matrix mass ratio, and interface modifier grafting rate. The theoretical elastic recovery strain parameters include in-mold shrinkage rate, post-demolding dimensional rebound, and warpage index. The diffusion characteristic length parameters include node-to-boundary Euclidean distance, effective diffusion path length, and cross-sectional geometry factor. The local hydrogen ion concentration index includes cumulative acidic product molar amount, local protonation degree index, and hydrolysis reaction activity concentration. The final formulation structural parameters include optimal component addition ratio, product macroscopic geometric wall thickness, and reinforcing rib distribution topology data.

3. The method for optimizing the formulation of biodegradable plate materials according to claim 1, characterized in that, The specific steps for obtaining the initial formulation component parameters are as follows: S101: Based on the atomic coordinate data and potential energy function parameters of starch and polylactic acid molecules, perform energy minimization iterative calculation, monitor the convergence state of the total potential energy, stop iteration when the derivative of the total potential energy is less than the preset convergence standard, lock the bond length and bond angle data between atoms, establish a three-dimensional spatial structure in which the atomic interaction force is in a state of force equilibrium, and generate the microscopic steady-state configuration of the composite system. S102: Discretize the microscopic steady-state configuration of the composite system into microscopic grid units, set the virtual probe radius parameter, traverse the microscopic grid units to perform space occupancy detection, mark the blank grid units not covered by the atomic van der Waals radius, aggregate adjacent connected blank grid units to form independent irregular cavity structures, extract the geometric boundary coordinates of the cavity structures, and establish a microscopic free volume topology set. S103: Call the microscopic free volume topology set, calculate the cumulative free volume fraction and average pore size of the cavity structure and mark them as volume parameters, calculate the work done by polymer chain segments to overcome intermolecular forces when crossing the cavity structure and mark them as chain segment transition energy level barrier parameters, and match the component concentration configuration based on the volume parameters and chain segment transition energy level barrier parameters to generate initial formulation component parameters.

4. The method for optimizing the formulation of biodegradable plate material according to claim 3, characterized in that, The specific process for obtaining the preset convergence criterion is as follows: The potential energy function parameters are called to resolve the bond stretching constant and filter out the stiffness coefficient with the largest value. The atomic coordinate data of starch and polylactic acid molecules are traversed to calculate the Euclidean distance between covalently connected atomic pairs. The minimum bond length feature value is extracted, and a fixed displacement tolerance ratio constant is obtained. The product of the minimum bond length feature value and the displacement tolerance ratio constant is calculated and marked as the maximum allowable residual displacement. The multiplication operation of the maximum stiffness coefficient and the maximum allowable residual displacement is performed to calculate the critical force amplitude for maintaining atomic dynamic equilibrium. The floating-point machine precision of the current computing environment is detected and converted into the minimum energy gradient limit. The values ​​of the critical force amplitude and the minimum energy gradient limit are compared, and the maximum value of the two is selected and set as the preset convergence criterion.

5. The method for optimizing the formulation of biodegradable plate material according to claim 3, characterized in that, The specific steps for obtaining the theoretical elastic recovery strain parameters are as follows: S201: Based on the initial formula component parameters, construct a virtual melt model, apply broadband sinusoidal shear boundary conditions, calculate stress response waveform data, decompose it into in-phase elastic component data and out-of-phase viscous component data through Fourier transform, calculate energy storage modulus parameters and loss modulus parameters, and generate dynamic rheological response modulus data. S202: Call the dynamic rheological response modulus data, establish a discretized integral kernel function matrix that correlates the frequency modulus and relaxation intensity, perform iterative error minimization operation to invert the relaxation intensity distribution, determine the modulus contribution weight at different time scales, and generate discrete relaxation time spectrum data; S203: Obtain the characteristic cooling and curing time parameters of the injection molding process, traverse the discrete relaxation time spectrum data, identify slow relaxation units whose time constant is greater than the cooling and curing time value, accumulate the modulus contribution weights corresponding to the slow relaxation units, calculate the residual elastic deformation potential energy, and generate theoretical elastic recovery strain parameters.

6. The method for optimizing the formulation of biodegradable plate material according to claim 5, characterized in that, The specific steps for obtaining the diffusion feature length parameter are as follows: S301: Call the theoretical elastic recovery strain parameters and map them point by point to the finite element mesh nodes of the three-dimensional geometric model of the plate. Construct a set of mechanical equilibrium equations containing geometric nonlinearity, perform numerical iteration to solve the three-dimensional deformation displacement vector of each node under residual stress, superimpose the displacement vector to the original coordinates to reconstruct the mesh shape, and generate macroscopic deformation geometric topology data. S302: Based on the macroscopic deformation geometric topology data, define the closed entity space, set the isotropic discretization step size, perform a global voxelization scan operation, convert the continuous geometric entity into a discrete voxel mesh array, remove the external background mesh and extract the three-dimensional Cartesian coordinates of the geometric center of the internal voxel unit, and establish a voxelized node space coordinate set. S303: Call the voxelized node spatial coordinate set, construct a spatial neighborhood index structure to traverse the outer surface mesh of the macroscopic deformation geometric topology data, retrieve the nearest neighbor boundary projection point corresponding to each voxel node, calculate the minimum Euclidean straight-line distance from the voxel center to the projection point, mark the distance value as a scalar field attribute, and generate a diffusion feature length parameter.

7. The method for optimizing the formulation of biodegradable plate material according to claim 6, characterized in that, The specific steps for obtaining the local hydrogen ion concentration index are as follows: S401: Call the initial formulation component parameters, extract the void volume ratio data based on the interface void cluster data, construct an effective diffusion transport model of solute in porous media, calculate the migration rate of substances in polymer matrix and mark it as substrate diffusion coefficient parameter, combine with diffusion characteristic length parameter, calculate the physical flux value of acidic substances migrating from the interior to the surface per unit time, and generate physical efflux flux parameter of acidic products. S402: Call the initial formulation component parameters, extract the molar concentration value of hydrolyzable ester bonds in the polylactic acid molecular chain, obtain the preset hydrolysis reaction rate constant, perform multiplication operation to calculate the chemical reaction rate of ester bond breaking to generate carboxyl terminus, quantify the amount of acidic terminus group material generated per unit volume per unit time, and generate acidic terminus group chemical generation flux parameter. S403: Calculate the numerical difference between the chemical generation flux parameter of the acidic terminal group and the physical efflux flux parameter of the acidic product, determine the net accumulation rate of acidic substances in the local area, add the difference data to the local acidic substance concentration parameter at the current time step, convert the molar concentration into a hydrogen ion activity value using the acid dissociation equilibrium constant, update the acidity state of the voxel grid node, and generate a local hydrogen ion concentration index.

8. The method for optimizing the formulation of biodegradable plate material according to claim 7, characterized in that, The specific process for obtaining the preset hydrolysis reaction rate constant is as follows: The initial formulation component parameters are called to analyze the ester bond chemical structure type of polylactic acid molecules. The standard ester bond hydrolysis activation energy parameter and pre-exponential factor physical quantity corresponding to the target chemical structure type are matched. The absolute temperature value set in the current simulation scenario is obtained. The standard molar gas constant is called to calculate the product of the absolute temperature value and the standard molar gas constant to generate a thermal noise energy benchmark. The division operation between the standard ester bond hydrolysis activation energy parameter and the thermal noise energy benchmark is performed. The result of the operation is inverted and the natural exponential function value is calculated to obtain the effective collision probability of molecules. The multiplication operation between the pre-exponential factor physical quantity and the effective collision probability of molecules is performed to generate the intrinsic rate value of chemical bond breaking. The intrinsic rate value is established as the preset hydrolysis reaction rate constant.

9. The method for optimizing the formulation of biodegradable plate material according to claim 7, characterized in that, The specific steps for obtaining the structural parameters of the finalized formula are as follows: S501: Call the local hydrogen ion concentration index to simulate the evolution of the autocatalytic reaction rate. Use numerical iteration to solve and dynamically update the hydrolysis reaction rate constant at each time step to simulate the mass loss during the material degradation process. When the remaining mass fraction is lower than the preset functional threshold, record the corresponding time span and mark it as the cumulative mass loss duration parameter to generate hydrolysis reaction kinetic evolution data. S502: Call the hydrolysis reaction kinetic evolution data, extract the time evolution curve of the hydrolysis reaction rate constant, calculate the slope change of adjacent time steps as the time change gradient parameter, detect the gradient value of each voxel grid node, identify the location where the gradient undergoes a step change as a potential fracture point, extract the corresponding three-dimensional spatial geometric coordinates, and generate a material structure failure topology set. S503: By comprehensively analyzing the hydrolysis reaction kinetic evolution data and the material structure failure topology set, a feature dataset integrating the degradation lifetime in the time dimension and the failure distribution in the spatial dimension is constructed. The dataset is then mapped back to the original formulation component space as a constraint condition for formulation screening. Combined with the degradation cycle and structural integrity requirements of the product design, the component ratio and geometric shape are screened to generate the finalized formulation structure parameters.

10. A biodegradable plate material formulation optimization system, characterized in that, The system is used to implement the method for optimizing the formulation of biodegradable plate materials according to any one of claims 1-9, the system comprising: The molecular conformation optimization module performs energy minimization calculations based on atomic coordinate data and potential energy function parameters, generates an equilibrium molecular conformation model, discretizes and aggregates interface void cluster data, calculates the volume of the interface void cluster data and the potential barrier parameters of the chain segment transition energy level, and generates the initial formulation component parameters. The elastic strain calculation module calculates the energy storage and loss modulus parameters based on the initial formula component parameters, converts them into discrete relaxation time spectrum data, accumulates the component weights of relaxation periods greater than the cooling and solidification characteristic duration parameters, and generates theoretical elastic recovery strain parameters. The geometric discretization module maps the theoretical elastic recovery strain parameters to the three-dimensional geometric model of the plate to calculate the deformation, discretizes it into voxel mesh nodes, calculates the distance from the node to the boundary and marks it as the diffusion feature length parameter. The concentration dynamic update module calls the initial formulation component parameters, extracts the interface void cluster data to calculate the substrate diffusion coefficient parameters, combines the diffusion characteristic length parameters to calculate the physical expulsion and chemical generation flux parameters, compares the difference and updates the local hydrogen ion concentration index of the voxel grid nodes. The structural failure integration module updates the reaction rate constant by the local hydrogen ion concentration index, calculates the cumulative mass loss duration parameter and time change gradient parameter, marks the coordinates of the structural failure location, and integrates to generate the finalized formula structural parameters.

Images