A simulation method for evaluating filler dispersion in high polymer materials

Abstract

This invention relates to the field of computational simulation technology for polymer materials, and particularly to a simulation method for evaluating filler dispersion in polymer materials. The method includes: acquiring at least two simulation snapshots generated during the potential energy convergence phase of a molecular dynamics simulation; dividing the simulation box space into three-dimensional grid regions; for each three-dimensional grid region, extracting polymer chain segments with a distance greater than a preset cutoff radius from the filler phase surface as a background reference frame, and calculating their first displacement vector within the inter-frame time interval; calculating the second displacement vector for the filler phase within the inter-frame time interval; calculating a local stability index based on the difference between the first and second displacement vectors; and statistically analyzing the local stability indices of all three-dimensional grid regions to generate a global stability assessment result, determining whether the filler dispersion is in a metastable state. This invention can identify pseudo-steady-state risks in advance during the potential energy convergence phase, and the assessment results are not affected by the filler shape, making it suitable for highly filled systems.

Description

Technical Field

[0001] This invention relates to the field of computational simulation technology for polymer materials, and more particularly to a simulation method for evaluating filler dispersion in polymer materials. Background Technology

[0002] In the design and development of polymer composites, adding nano- or micro-sized fillers (such as silica, carbon black, carbon nanotubes, graphene, etc.) to the polymer matrix to improve the mechanical, electrical, thermal, or barrier properties of the material has become a routine practice in industry. The uniformity and stability of filler dispersion directly determine the final service performance of the composite material. To shorten the R&D cycle and reduce trial-and-error costs, molecular dynamics simulation methods are commonly used in this field to virtually simulate the dispersion process of fillers in polymer melts, and the dispersion effect of fillers is evaluated by analyzing the simulation trajectory data. Existing evaluation methods mainly include: calculating the radial distribution function of the filler to determine the degree of filler aggregation; calculating the change of the total potential energy of the system over time to determine whether the simulation has reached thermodynamic equilibrium; and visually observing the simulation snapshot to qualitatively determine the dispersion morphology.

[0003] However, the existing assessment methods share a common technical blind spot: in highly packed systems or systems with strong interfacial interactions, the movement of polymer chains is severely hindered near the filler surface, making the simulated system prone to thermodynamic metastable states. Specifically, within a limited simulation time window (e.g., less than 100 nanoseconds), the total potential energy curve of the system tends to flatten, and the peak shape of the radial distribution function no longer shows any visually perceptible change. Existing software or analysis procedures thus determine that the filler dispersion has reached a stable state and terminate the calculation. However, in actual physical processes, if the external shear field is removed or the settling time is extended, the filler in this metastable state will undergo slow secondary aggregation due to entropy-induced cavitation or interfacial slippage, leading to a significant discrepancy between the simulation assessment conclusions and experimental observations. Current technology lacks an assessment method that can effectively predict whether the filler dispersion is truly stable and whether there is a risk of secondary aggregation based solely on simulated trajectory data in the early stages of potential energy convergence. Summary of the Invention

[0004] To overcome the above shortcomings, this invention provides a simulation method for evaluating filler dispersion in polymer materials, aiming to improve the problem that existing molecular dynamics simulation evaluation methods cannot identify that the filler dispersion is in a thermodynamic metastable state after the system potential energy converges, thus failing to provide early warning of the risk of secondary agglomeration.

[0005] This invention provides the following technical solution: a simulation method for evaluating filler dispersion in polymer materials, comprising the following steps: S1. Obtain at least two simulation snapshots generated during the system potential energy convergence phase of the molecular dynamics simulation, and divide the simulation box space into a three-dimensional grid region with a preset side length. S2. For each of the three-dimensional grid regions, extract the polymer chain segments in the grid region that are at a distance greater than a preset cutoff radius from the surface of the filler phase as the background reference system, and calculate the first displacement vector of the background reference system in the inter-frame time interval based on the at least two simulated snapshots. S3. For the filler phase within the three-dimensional grid region, calculate the second displacement vector of the filler phase within the inter-frame time interval based on the at least two simulated snapshots; S4. Based on the difference between the first displacement vector and the second displacement vector, calculate the local stability index, which characterizes the degree of decoupling between the filler phase and the bulk polymer chain segments within the three-dimensional grid region. S5. Calculate the local stability index of all three-dimensional grid regions, generate a global stability assessment result, and determine whether the packing dispersion is in a metastable state based on the global stability assessment result.

[0006] Preferably, in step S1, the step of obtaining at least two simulation snapshots generated during the system potential energy convergence phase of the molecular dynamics simulation specifically includes: Read the original trajectory file output by the molecular dynamics simulation and extract the energy curve data of the total potential energy of the system as a function of time; In the energy curve data, identify the continuous time period in which the fluctuation amplitude of the total potential energy of the system is continuously less than a preset energy threshold, and mark this continuous time period as the potential energy convergence stage of the system. From the original trajectory file, extract all simulation snapshots whose timestamps are within the potential energy convergence phase of the system to form a candidate snapshot set; From the candidate snapshot set, at least two simulated snapshots whose time interval satisfies the preset inter-frame interval condition are selected as the at least two simulated snapshots output.

[0007] Preferably, in step S1, the step of dividing the simulated box space into three-dimensional mesh regions with preset side lengths specifically includes: Read any one of the at least two simulated snapshots and obtain the box length in the three dimensions of the simulated box space in that simulated snapshot. Obtain the preset grid side length parameters, calculate the ratio of the box length in each of the three dimensions to the grid side length parameters, and round each ratio up to obtain the number of grid divisions in the three dimensions. Based on the number of mesh divisions in the three dimensions, the simulation box is divided into multiple three-dimensional mesh regions by equal intervals in the three dimensions. Assign a unique spatial index identifier to each of the three-dimensional mesh regions.

[0008] Preferably, in step S2, the step of extracting polymer chain segments within the grid region that are at a distance greater than a preset cutoff radius from the surface of the filler phase as the background reference specifically includes: For the current three-dimensional mesh region, obtain the monomer coordinates of all polymer chain segments within the three-dimensional mesh region, and obtain the surface atomic coordinates of all filler phase atoms within the three-dimensional mesh region; For each monomer of the polymer chain segment, calculate the shortest spatial distance from the monomer to the surface atoms of the filler phase; The shortest spatial distance is compared with the preset cutoff radius, and individual units whose shortest spatial distance is greater than the preset cutoff radius are selected to form a candidate individual set; Determine whether the proportion of the number of monomers in the candidate monomer set to the total number of polymer chain segment monomers in the three-dimensional grid region is greater than a preset proportion threshold. If the ratio is greater than the preset threshold, the polymer chain segment to which the candidate monomer set belongs is extracted as the background reference system of the three-dimensional grid region; If the ratio is less than or equal to the preset ratio threshold, then the parent coarse-grained grid region to which the three-dimensional grid region belongs is retrieved, and the polymer chain segments in the parent coarse-grained grid region that satisfy the condition that the shortest spatial distance is greater than the preset cutoff radius are extracted as the background reference system of the three-dimensional grid region.

[0009] Preferably, in step S2, the step of calculating the first displacement vector of the background reference frame within the inter-frame time interval based on the at least two simulated snapshots specifically includes: For the current 3D mesh region, determine the first simulation snapshot with the earlier timestamp and the second simulation snapshot with the later timestamp from the at least two simulation snapshots; Obtain the first centroid coordinates of the background reference system in the first frame simulation snapshot, and obtain the second centroid coordinates of the background reference system in the second frame simulation snapshot; Calculate the spatial vector difference between the second centroid coordinates and the first centroid coordinates, and determine the spatial vector difference as the first displacement vector of the background reference system within the inter-frame time interval.

[0010] Preferably, in S3, the step of calculating the second displacement vector of the filler phase within the inter-frame time interval based on the at least two simulated snapshots specifically includes: For the filler phase within the current three-dimensional mesh region, calculate the radius of gyration tensor of the filler phase, and determine the anisotropy of the filler phase based on the eigenvalues ​​of the radius of gyration tensor. Determine whether the anisotropy degree is greater than a preset anisotropy degree threshold; If the anisotropy is less than or equal to the preset anisotropy threshold, then the centroid displacement vector of the filler phase within the inter-frame time interval is obtained, and the centroid displacement vector is determined as the second displacement vector. If the anisotropy is greater than the preset anisotropy threshold, the overall displacement field of the filler phase during the inter-frame time interval is calculated, the centroid translation component is decomposed from the overall displacement field, and the centroid translation component is determined as the second displacement vector.

[0011] Preferably, in step S4, the step of calculating the local stability index, which characterizes the degree of decoupling between the filler phase and the bulk polymer chain segments within the three-dimensional grid region, specifically includes: For the current three-dimensional mesh region, calculate the spatial vector difference between the second displacement vector and the first displacement vector to obtain the difference displacement vector; Calculate the vector magnitude of the difference displacement vector, and calculate the vector magnitude of the first displacement vector; Calculate the ratio of the vector magnitude of the difference displacement vector to the vector magnitude of the first displacement vector, and determine this ratio as the local stability index of the three-dimensional mesh region.

[0012] Preferably, in step S5, the step of generating the global stability assessment result specifically includes: Collect all the calculated local stability indices of the three-dimensional mesh regions to form a set of local stability indices; Calculate the statistical characteristic values ​​of the set of local stability indices; The statistical characteristic values ​​are output as the global stability assessment results.

[0013] Preferably, in step S5, the step of determining whether the filler dispersion is in a metastable state based on the global stability assessment result specifically includes: Obtain the global stability assessment result and compare the global stability assessment result with a preset stability threshold; If the global stability assessment result is less than the preset stability threshold, the filler dispersion is determined to be in a true steady state, and a first assessment conclusion characterizing the dispersion stability is output. If the global stability assessment result is greater than or equal to the preset stability threshold, the filler dispersion is determined to be in a metastable state, and a second assessment conclusion characterizing the dispersion instability is output.

[0014] The present invention has the following beneficial effects: 1. In this invention, by calculating the degree of motion decoupling between the filler phase and the background reference system, rather than simply relying on the static convergence criterion of total potential energy or radial distribution function, the independent drift trend of the filler relative to the bulk polymer chain segment can be captured in the early stage when the total energy of the system has tended to stabilize. Thus, a metastable state warning is issued before secondary aggregation actually occurs, effectively avoiding the prediction deviation between the simulation evaluation conclusion and the actual physical process.

[0015] 2. In this invention, before calculating the displacement vector of the packing, the anisotropy of the packing is determined by the radius of gyration tensor. For packings with high aspect ratio or two-dimensional lamellar packing, only the translational component of the center of mass is extracted for subsequent calculation. This eliminates the noise contribution of the in-situ rotation of slender packing caused by thermal disturbance to the displacement signal, so that the local stability index can truly reflect the degree of danger of the packing migrating towards agglomeration, and significantly reduces the false positive rate.

[0016] 3. In this invention, when extracting the background reference system, a parent coarse-grained grid backtracking mechanism is introduced to address the extreme case where dense filler leads to a lack of bulk polymer chain segments in the grid area. The average chain segment displacement at a larger spatial scale is used as a substitute reference system, avoiding calculation interruptions or division-to-zero anomalies caused by the lack of a reference system. This ensures the robustness and continuous operation capability of the evaluation method under any filler content. Attached Figure Description

[0017] Figure 1 This is a schematic flowchart of a simulation method for evaluating filler dispersion in polymer materials proposed in this invention. Detailed Implementation

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

[0019] This invention provides a simulation method for evaluating filler dispersion in polymer materials, such as... Figure 1 As shown, it includes the following steps: S1. Obtain at least two simulation snapshots generated during the system potential energy convergence phase of the molecular dynamics simulation, and divide the simulation box space into a three-dimensional grid region with a preset side length.

[0020] Furthermore, in S1, the steps of obtaining at least two simulation snapshots generated during the system potential energy convergence phase of the molecular dynamics simulation specifically include: Read the original trajectory file output by the molecular dynamics simulation and extract the energy curve data of the total potential energy of the system as a function of time; In the energy curve data, identify the continuous time period in which the fluctuation amplitude of the total potential energy of the system is continuously less than the preset energy threshold, and mark this continuous time period as the potential energy convergence stage of the system. From the original trajectory file, extract all simulation snapshots whose timestamps fall within the system potential energy convergence phase to form a candidate snapshot set; From the candidate snapshot set, select at least two simulated snapshots whose time interval meets the preset inter-frame interval condition, and output them as at least two simulated snapshots.

[0021] Furthermore, in S1, the step of dividing the simulated box space into three-dimensional mesh regions with preset side lengths specifically includes: Read any one of the at least two simulated snapshots and obtain the box length in the three dimensions of the simulated box space in that simulated snapshot. Obtain the preset grid side length parameters, calculate the ratio of the box length to the grid side length parameters in each of the three dimensions, and round each ratio up to obtain the number of grid divisions in each of the three dimensions; Based on the number of mesh divisions in the three dimensions, the simulation box is divided into three equal-distance subdivisions in the three dimensions to generate multiple three-dimensional mesh regions. Assign a unique spatial index identifier to each 3D grid region.

[0022] Specifically, the raw trajectory file output by the molecular dynamics simulation software is read. This raw trajectory file contains the atomic coordinates and total potential energy of the simulated system at each time point. Energy curve data showing the change of the total potential energy over time is extracted from this raw trajectory file. This energy curve data consists of a series of data pairs consisting of a time point and the corresponding total potential energy value. The fluctuation amplitude of the total potential energy within the sliding time window is calculated from the energy curve data. This fluctuation amplitude can be characterized by the standard deviation of the total potential energy within the sliding time window, and its calculation formula is: ; in, Let be the standard deviation of the total potential energy of the system within a sliding time window centered at time t. The number of data points contained within the sliding time window. For a moment The total potential energy of the system, It is the arithmetic mean of the total potential energy of the system within the sliding time window.

[0023] When the standard deviation is consistently less than the preset energy threshold The length of the continuous time interval exceeds the preset minimum convergence time. When this continuous time period is reached, it is marked as the potential energy convergence phase of the system. A preset energy threshold is used. Based on the total number of atoms and temperature settings of the simulation system, the value is typically three to five times the root mean square fluctuation of the potential energy in the final stage of the simulation. A minimum convergence time is preset. Used to eliminate the brief plateau period that occurs during the descent of the potential energy curve, it is usually set to five to ten nanoseconds.

[0024] From the original trajectory file, extract all simulation snapshots whose timestamps fall within the potential energy convergence phase of the system, forming a candidate snapshot set. Each simulation snapshot must contain at least the three-dimensional coordinate data of all atoms within the simulation box at that timestamp and the box length data in three dimensions of the simulation box. From the candidate snapshot set, select at least two simulation snapshots whose time interval satisfies a preset inter-frame interval condition as at least two simulation snapshots to be used in subsequent steps. The preset inter-frame interval condition is that the time difference between two adjacent selected simulation snapshots satisfies the following: ; in, The timestamp of the simulated snapshot selected for the k-th frame. The timestamp of the simulated snapshot selected for the (k+1)th frame. Preset minimum frame interval. The typical value for the time interval is one to five nanoseconds to ensure that the filler phase generates an observable displacement signal with respect to the background reference frame between two frames. At least two selected simulation snapshots can be used directly for subsequent displacement vector calculations, or multiple simulation snapshots can be uniformly extracted from the candidate snapshot set to form a snapshot sequence, which can then be used for time averaging of the multi-frame displacement vectors to reduce thermal noise.

[0025] Read any one of at least two simulation snapshots, and obtain the box lengths in the three dimensions of the simulated box space from the header information or boundary coordinate information of that simulation snapshot, denoted as . Since the size of the simulated box space usually remains constant or fluctuates only within a very small range during the system's potential energy convergence phase, selecting any frame of the simulation snapshot to read the box size is sufficient to meet the accuracy requirements of subsequent mesh generation.

[0026] Obtain the preset grid side length parameters, denoted as... Mesh side length parameters The value of is determined based on the average particle size of the filler phase, and is usually set to two to three times the average particle size of the filler to ensure that each three-dimensional mesh region contains a sufficient amount of polymer chain monomers, while avoiding the local details from being submerged by the overall set average due to an excessively large mesh scale. The ratio of the box length to the mesh side length parameter in each of the three dimensions is calculated, and each ratio is rounded up to obtain the number of mesh divisions in each of the three dimensions. The calculation formula is as follows: ; ; ; in, , , These represent the number of grid divisions in the three dimensions, with symbols... This indicates the rounding up operation. The rounding up operation ensures that the mesh completely covers the entire simulation box space, leaving no boundary gaps.

[0027] Based on the number of grid divisions in three dimensions , , The simulated box space is divided into three equal-distance sections along its three dimensions. Specifically, along the x-direction... Divide the data by step size, in the y-direction with... Segment by step size, in the z direction with Segmentation is performed based on the step size. After segmentation, a total of [number] segments are generated. × × Each 3D mesh region is a cuboid subspace with side lengths in each of the three dimensions. .

[0028] Assign a unique spatial index identifier to each 3D mesh region. This unique spatial index identifier can be a triple consisting of the indices of the 3D mesh region in the three dimensions, in the form of (i,j,k), where i satisfies 1≤i≤k. The integer j is a set of integers that satisfy 1 ≤ j ≤ 1. An integer k such that 1 ≤ k ≤ The unique spatial index is an integer. This unique spatial index is used to identify the same 3D mesh region in subsequent steps and to aggregate the various types of data calculated within that 3D mesh region to the corresponding storage location. The one-to-one correspondence between the unique spatial index of all 3D mesh regions and their spatial locations is recorded in a mesh index table.

[0029] Through the above steps, an effective simulation snapshot of the potential energy convergence stage was extracted from the original trajectory file, and the simulation box space was divided into a three-dimensional grid array that facilitates local quantitative analysis. This provides a data foundation and spatial positioning basis for the independent extraction of the background reference system and calculation of the first displacement vector in each three-dimensional grid region in the subsequent step S2.

[0030] S2. For each three-dimensional grid region, extract the polymer chain segment in the grid region that is at a distance greater than the preset cutoff radius from the surface of the filler phase as the background reference system, and calculate the first displacement vector of the background reference system in the inter-frame time interval based on at least two simulated snapshots.

[0031] Furthermore, in S2, the step of extracting polymer chain segments within the grid region that are at a distance greater than a preset cutoff radius from the filler phase surface as the background reference specifically includes: For the current three-dimensional mesh region, obtain the monomer coordinates of all polymer chain segments within the three-dimensional mesh region, and obtain the surface atomic coordinates of all filler phase atoms within the three-dimensional mesh region; For each monomer in a polymer chain segment, calculate the shortest spatial distance from that monomer to the surface atoms of the filler phase. The shortest spatial distance is compared with the preset cutoff radius, and the individual units whose shortest spatial distance is greater than the preset cutoff radius are selected to form a candidate individual set; Determine whether the proportion of the number of monomers in the candidate monomer set to the total number of polymer chain segments in the three-dimensional grid region is greater than a preset proportion threshold. If the ratio exceeds a preset threshold, the polymer chain segment to which the candidate monomer set belongs will be extracted as the background reference system for that three-dimensional grid region. If the ratio is less than or equal to a preset threshold, the parent coarse-grained grid region to which the three-dimensional grid region belongs is retrieved, and the polymer chain segments in the parent coarse-grained grid region that satisfy the condition that the shortest spatial distance is greater than a preset cutoff radius are extracted as the background reference system of the three-dimensional grid region.

[0032] Furthermore, in S2, the step of calculating the first displacement vector of the background reference frame within the inter-frame time interval based on at least two simulated snapshots specifically includes: For the current 3D mesh region, determine the first simulation snapshot with the earlier timestamp and the second simulation snapshot with the later timestamp from at least two simulation snapshots; Obtain the first centroid coordinates of the background reference frame in the first frame of the simulation snapshot, and obtain the second centroid coordinates of the background reference frame in the second frame of the simulation snapshot; Calculate the spatial vector difference between the second centroid coordinates and the first centroid coordinates, and determine this spatial vector difference as the first displacement vector of the base reference frame within the inter-frame time interval.

[0033] Specifically, for the current 3D mesh region, the unique spatial index identifier assigned in the above steps is used to locate the atomic data within this 3D mesh region. The monomer coordinates of all polymer chain segments within this 3D mesh region are obtained, forming a set of polymer monomer coordinates, denoted as . Where the superscript m represents the sequence number of the polymer monomer. Simultaneously, the surface atomic coordinates of all filler phase atoms within this three-dimensional grid region are obtained, forming a set of filler surface coordinates, denoted as . The superscript n represents the serial number of the surface atom of the filler phase. The method for identifying the surface atom of the filler phase is as follows: if the distance from a certain filler phase atom to the geometric center of its corresponding filler particle is greater than the average distance of the inner layer atoms of the filler particle, or if the coordination number of the filler phase atom is less than the saturation coordination number of the bulk phase atoms of the filler phase, then the filler phase atom is marked as a surface atom of the filler phase.

[0034] For polymer monomer coordinate sets Monomer coordinates of each polymer monomer in Calculate the set of coordinates from the monomer to the packing surface. The shortest spatial distance between atoms on the surface of all filler phases. The formula for calculating this shortest spatial distance is: ; in, Let m be the shortest spatial distance from the m-th polymer monomer to the surface atoms of the filler phase, denoted by . This represents the Euclidean norm, and `min` indicates the minimum value operation. To improve computational efficiency, a KD-Tree spatial indexing structure can be used to index the coordinate set of the packing surface. Preprocessing reduces the time complexity of a single shortest distance query to logarithmic levels.

[0035] The shortest spatial distance of each polymer monomer was calculated. With preset cutoff radius Compare. Preset cutoff radius. The value is the cutoff distance of the van der Waals interaction between the polymer chain segment and the filler phase, typically set to 0.8 nm to 1.2 nm. Screening is performed to select those that meet the requirements. > The conditions for high molecular weight monomers constitute a candidate monomer set. .

[0036] Statistical candidate monomer set The number of monomers in the mixture is denoted as Simultaneously, the total number of polymer chain segments and monomers within this three-dimensional grid region is counted and denoted as . The proportion of candidate monomers to the total number of monomers is calculated using the following formula: ; Determine whether the ratio R is greater than the preset ratio threshold. Preset ratio threshold The typical value is between 10% and 20%. If R > This indicates that there are sufficient bulk polymer chain segments far from the filler phase surface within the current 3D mesh region, and the first step is executed. If R ≤ If the result is negative, it indicates that the filler phase is densely distributed within the current three-dimensional grid region, lacking sufficient bulk polymer chain segments, and the second step should be executed.

[0037] The first step is: when R> At that time, the candidate singleton set The extracted polymer chain segments are used as the background reference frame for this three-dimensional grid region. The extracted background reference frame consists of several complete polymer chains or polymer chain segments. These segments are spatially located far from the filler phase surface, and their motion characteristics approximate those of bulk chain segments in a pure polymer melt, objectively reflecting the background thermal motion level of the polymer matrix in this local region. The second step is: when R ≤ At that time, the parent coarse-grained mesh region to which the 3D mesh region belongs is retrieved. The parent coarse-grained mesh region refers to the spatial region obtained by dividing the mesh with a larger mesh side length parameter above the mesh division level in step S1. Specifically, adjacent 2×2×2 3D mesh regions can be merged into a parent coarse-grained mesh region. Within this parent coarse-grained mesh region, the above steps are repeated to satisfy... > The extracted polymer chain segments serve as the background reference frame for the given 3D mesh region. If a background reference frame meeting the required proportions cannot be extracted within the parent coarse-grained mesh region, the process can continue backtracking to higher-level coarse-grained mesh regions, up to the entire simulation box space. This backtracking mechanism ensures that a valid background reference frame can be extracted for each 3D mesh region regardless of filler content, avoiding computational interruptions due to a missing reference frame.

[0038] For the current 3D mesh region, from at least two simulation snapshots obtained in step S1, determine the first simulation snapshot with the earlier timestamp and the second simulation snapshot with the later timestamp. Record the timestamp of the first simulation snapshot as... The timestamp of the second frame simulated snapshot is recorded as The inter-frame time interval is: ; Obtain the first centroid coordinates of the background reference frame extracted through the above steps in the first frame of the simulation snapshot. And obtain the second centroid coordinates of the background reference frame in the second frame simulation snapshot. The formula for calculating the centroid coordinates of the background reference frame is: ; in, This represents the total number of polymer monomers included in the background reference system. Let be the coordinate vector of the m-th polymer monomer. Substituting the monomer coordinates from the first and second simulated snapshots into the above formula, respectively, yields the coordinates of the first centroid. Second centroid coordinates .

[0039] The spatial vector difference between the coordinates of the second centroid and the coordinates of the first centroid is calculated using the following formula: ; The spatial vector difference The baseline reference frame is determined during the inter-frame time interval. The first displacement vector within the frame. This first displacement vector characterizes the net displacement of the bulk polymer chain segments far from the filler phase surface within the current three-dimensional grid region during the inter-frame time interval, reflecting the amplitude and direction of the background thermal motion of the polymer matrix in this local region.

[0040] S3. For the filler phase within the three-dimensional grid region, calculate the second displacement vector of the filler phase within the inter-frame time interval based on at least two simulated snapshots.

[0041] Furthermore, in S3, the step of calculating the second displacement vector of the filler phase within the inter-frame time interval based on at least two simulated snapshots specifically includes: For the filler phase within the current three-dimensional mesh region, calculate the radius of gyration tensor of the filler phase, and determine the anisotropy of the filler phase based on the eigenvalues ​​of the radius of gyration tensor. Determine whether the anisotropy degree is greater than the preset anisotropy degree threshold; If the anisotropy is less than or equal to the preset anisotropy threshold, the centroid displacement vector of the filler phase in the inter-frame time interval is obtained, and the centroid displacement vector is determined as the second displacement vector. If the anisotropy is greater than the preset anisotropy threshold, the overall displacement field of the filler phase during the inter-frame time interval is calculated, the centroid translation component is decomposed from the overall displacement field, and the centroid translation component is determined as the second displacement vector.

[0042] Specifically, for the filler phase within the current 3D mesh region, the atomic coordinates of all filler phase atoms contained in that filler phase are extracted. Assume that the filler phase consists of... It is composed of atoms of a filler phase, and its atomic coordinate set is as follows: Where i is the atomic number of the filler phase, and its value ranges from 1 to i. First, calculate the centroid coordinates of the packing phase. The calculation formula is as follows: ; Based on the centroid coordinates, the radius of gyration tensor S of the packing phase is calculated. The radius of gyration tensor S is a 3×3 symmetric matrix, and its matrix elements are calculated using the following formula: ; in, and These represent the spatial coordinate directions, with values ​​of x, y, and z; For the i-th filler phase atom in Coordinate components in the direction; For the filler phase centroid in Coordinate components in the direction.

[0043] Eigenvalue decomposition of the radius of gyration tensor S yields three eigenvalues, denoted as . And satisfy These three eigenvalues ​​characterize the spatial extent of the packing phase along three mutually orthogonal principal axis directions. The anisotropy of the packing phase is determined based on these eigenvalues. The calculation formula is as follows: ; in, To prevent extremely small constants with a denominator of zero, a typical value is taken as... Anisotropy The larger the value, the more the shape of the filler phase deviates from a spherical shape, that is, the larger the aspect ratio or the thinner the laminations.

[0044] The calculated anisotropy degree Compared with the preset anisotropy threshold Comparison is performed. A preset anisotropy threshold is used. The typical value ranges from 1.5 to 2.0. If ≤ This indicates that the shape of the packing phase is approximately spherical, and its rotational motion has negligible interference with the measurement of the center of mass displacement. If > If the phase is high aspect ratio fibrous packing or high aspect ratio two-dimensional lamellar packing, the in-situ rotation caused by thermal disturbance will produce significant end displacement, and the rotation component needs to be separated from the overall displacement.

[0045] when ≤ At that time, the centroid displacement vector of the filler phase within the inter-frame time interval is directly obtained, and this centroid displacement vector is determined as the second displacement vector. Specifically, from the at least two simulated snapshots obtained in step S1, the first simulated snapshot with the earlier timestamp and the second simulated snapshot with the later timestamp are determined. The centroid coordinates of the filler phase in the first simulated snapshot are calculated respectively. and the centroid coordinates in the second frame simulation snapshot The formula for calculating the centroid coordinates is the same as the formula for calculating the centroid coordinates in the steps above. The second displacement vector of the filler phase during the inter-frame time interval. The calculation formula is: ; when > At that time, the overall displacement field of the filler phase during the inter-frame time interval is calculated, and the centroid translational component is decomposed from the overall displacement field. This centroid translational component is then determined as the second displacement vector. The specific decomposition method is as follows: First, the set of atomic coordinates of the filler phase in the first frame simulation snapshot is obtained. and the set of atomic coordinates in the second frame simulation snapshot The least squares method is used to fit the rigid body motion transformation of the filler phase from the first frame to the second frame. This rigid body motion transformation is described by the translation vector T and the rotation matrix Q. The optimal translation vector T and the optimal rotation matrix Q are obtained by solving the following minimization problem: ; The minimization problem can be solved using singular value decomposition or quaternion methods, which are well-known rigid body registration algorithms in this field. The optimal translational vector T obtained by the solution is the centroid translational component decomposed from the global displacement field. This centroid translational component T is determined as the second displacement vector of the filler phase within the inter-frame time interval, i.e.: ; Through the above decomposition, the second displacement vector only reflects the true translational trend of the packing phase centroid, eliminating the contribution of the in-situ rotation caused by thermal disturbance of slender packing or sheet packing to the displacement signal.

[0046] S4. Based on the difference between the first displacement vector and the second displacement vector, calculate the local stability index, which characterizes the degree of decoupling between the filler phase and the bulk polymer chain segments within the three-dimensional grid region.

[0047] Furthermore, in S4, the step of calculating the local stability index, which characterizes the degree of decoupling between the filler phase and the bulk polymer chain segments within the three-dimensional grid region, specifically includes: For the current three-dimensional mesh region, calculate the spatial vector difference between the second displacement vector and the first displacement vector to obtain the difference displacement vector; Calculate the vector magnitude of the difference displacement vector, and calculate the vector magnitude of the first displacement vector; Calculate the ratio of the vector magnitude of the difference displacement vector to the vector magnitude of the first displacement vector, and determine this ratio as the local stability index of the three-dimensional mesh region.

[0048] Specifically, for the current 3D mesh region, obtain the first displacement vector output from the above steps. Second displacement vector Calculate the spatial vector difference between the second displacement vector and the first displacement vector to obtain the difference displacement vector. The calculation formula is as follows: ; The difference displacement vector The physical meaning of this value lies in the fact that it characterizes the deviation between the actual displacement of the filler phase and the displacement of the background thermal motion of the polymer matrix represented by the background reference frame within the same inter-frame time interval. If the movement of the filler phase is completely synchronized with the thermal motion of the polymer matrix, then this difference displacement vector should approach a zero vector. If the filler phase exhibits directional migration independent of the polymer matrix, then this difference displacement vector will have a significant non-zero modulus.

[0049] Calculate the vector magnitudes of the difference displacement vector and the first displacement vector separately. The vector magnitude is calculated as the square root of the sum of the squares of the vector components. (The vector magnitude of the difference displacement vector is...) The calculation formula is: ; in, , , They are the difference displacement vectors. The components in the x, y, and z coordinate directions.

[0050] Vector magnitude of the first displacement vector The calculation formula is: ; in, , , The first displacement vectors are respectively The components in the x, y, and z coordinate directions.

[0051] The ratio of the magnitude of the difference displacement vector to the magnitude of the first displacement vector is calculated, and this ratio is determined as the local stability index of the 3D mesh region. The local stability index is denoted as... The calculation formula is as follows: ; in, To prevent extremely small constants with a denominator of zero, a typical value is taken as... When the vector magnitude of the first displacement vector When it is extremely small, The introduction of this method can prevent numerical overflow caused by division operations.

[0052] Local stability index It has clear physical interpretability. When When the value approaches zero, it indicates that the modulus of the displacement vector of the filler phase is basically consistent with the modulus of the displacement vector of the background reference system, and their directions are also basically consistent. The movement of the filler phase is highly coupled with the thermal movement of the polymer matrix, and the dispersion state of this local region is a true steady state. When the value is significantly greater than zero, it indicates that the filler phase has an independent displacement relative to the polymer matrix. The direction and amplitude of this displacement are decoupled from the background thermal motion, and this local region shows a tendency for directional migration of the filler phase, which is a potential risk area for secondary agglomeration. When the value is greater than or equal to one, it indicates that the independent displacement amplitude of the filler phase has exceeded the baseline thermal motion amplitude of the polymer matrix, and the dispersion state of this local area is in a highly unstable metastable state.

[0053] S5. Calculate the local stability index of all three-dimensional grid regions, generate the global stability assessment result, and determine whether the packing dispersion is in a metastable state based on the global stability assessment result.

[0054] Furthermore, in S5, the steps for generating global stability assessment results specifically include: Collect all calculated local stability indices for the three-dimensional mesh regions to form a set of local stability indices; Calculate the statistical characteristic values ​​of the set of local stability indices; The statistical eigenvalues ​​are output as the results of the global stability assessment.

[0055] Furthermore, in S5, the step of determining whether the filler dispersion is in a metastable state based on the global stability assessment results specifically includes: Obtain the global stability assessment results and compare them with the preset stability threshold; If the global stability assessment result is less than the preset stability threshold, the packing dispersion is determined to be in a true steady state, and the first assessment conclusion characterizing the dispersion stability is output. If the global stability assessment result is greater than or equal to the preset stability threshold, the packing dispersion is determined to be in a metastable state, and a second assessment conclusion characterizing the dispersion instability is output.

[0056] Specifically, all the 3D mesh regions generated in step S1 are traversed. Using the unique spatial index assigned in the previous steps, the data storage location corresponding to each 3D mesh region is accessed one by one, and the calculated local stability indices are collected. The local stability indices of all 3D mesh regions are summarized to form a set of local stability indices, denoted as […]. Where i is the index of the three-dimensional mesh region, The total number of 3D mesh regions, its value is equal to .

[0057] For the set of local stability indices Perform statistical analysis to calculate the statistical characteristic values ​​of the set. These statistical characteristic values ​​are used to characterize the overall distribution level and extreme cases of the local stability index within the simulated box space from different perspectives. The statistical characteristic values ​​can be at least one of the following statistical measures: the first statistical characteristic value is the arithmetic mean, denoted as... The calculation formula is as follows: ; The arithmetic mean reflects the overall average level of the local stability index for all three-dimensional mesh regions.

[0058] The second statistical characteristic value is the median, denoted as... It is defined as the set of local stability indices. The median is the value in the middle position after being sorted by numerical value. The median is not sensitive to extreme outliers in the set of local stability indices and can more robustly reflect the central tendency of the data.

[0059] The third statistical characteristic value is the maximum value, denoted as... The calculation formula is as follows: ; The maximum value reflects the stability of the local region with the most severe motion decoupling in the simulated box space, and is indicative of the identification of local high-risk regions.

[0060] The fourth statistical characteristic value is the high percentile value. A preset percentile P is selected, typically the 95th or 99th percentile, to set the local stability indices. After sorting the values ​​in ascending order, the value corresponding to the Pth percentile is taken as the statistical characteristic value, denoted as . High percentile values ​​can eliminate the interference of a very small number of anomalous meshes, while capturing the behavioral characteristics of the top few percent of local regions with the worst stability in the system. In practical applications, one or more of the above-mentioned arithmetic mean, median, maximum value, or high percentile value can be selected as statistical characteristic values ​​according to the type of filler and the required evaluation accuracy.

[0061] The statistical characteristic values ​​calculated in the above steps are output as the global stability assessment result. The global stability assessment result is denoted as... If multiple statistical characteristic values ​​were calculated in the above steps, then It can be an evaluation result vector composed of multiple statistics, or it can select the most sensitive statistic as the final scalar output according to preset rules. Typically, high percentile values ​​can be used. As a result of global stability assessment This focuses on monitoring the behavior of the most unstable local areas within the monitoring system.

[0062] Obtain the global stability assessment results output by the above steps. The global stability assessment result is compared with the preset stability threshold. Compare the results. Preset stability threshold. The value of is determined as follows: Several standard simulation systems with known dispersion stability are selected as calibration samples, and the global stability assessment results of each calibration sample are calculated. The optimal boundary value that can correctly distinguish between stable and unstable systems is used as the preset stability threshold. Typically, the preset stability threshold is... The value range is from 0.2 to 0.5.

[0063] If the comparison result is If the packing dispersion is deemed to be in a true steady state, it indicates that the movement of the packing phase within the entire simulation system is highly synchronized with the background thermal motion of the polymer matrix represented by the background reference frame, and there is no significant independent directional migration trend. Therefore, the packing dispersion structure is thermodynamically stable. Under this determination, a first assessment conclusion characterizing the dispersion stability is output. The specific form of the first assessment conclusion may include: outputting a text message to the graphical user interface stating "High confidence level of dispersion structure, current dispersion state is stable"; displaying a green indicator indicating stability in the assessment report interface; or sending a command signal to the simulation control module to terminate the current simulation calculation.

[0064] If the comparison result is If the particle size distribution is not found to be stable, the filler dispersion is determined to be in a metastable state. This indicates that although the total potential energy of the system may have converged, the filler phase exhibits independent directional migration relative to the polymer matrix, and the filler dispersion structure is in a thermodynamically metastable state. There is a risk of secondary aggregation during subsequent resting or long-term relaxation. Under this determination, a second assessment conclusion characterizing the dispersion instability is output. The specific form of the second assessment conclusion may include: outputting a text warning message to the graphical user interface stating "Local slip pseudo-steady state detected; current energy convergence is an illusion; it is recommended to extend the relaxation time"; displaying a red warning indicator in the assessment report interface; or sending a prompt to the simulation control module indicating that the simulation time needs to be extended or the simulation parameters adjusted.

[0065] Finally, it should be noted that the above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art can still modify the technical solutions described in the foregoing embodiments or make equivalent substitutions for some of the technical features. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A simulation method for assessing filler dispersion in a high polymer material, characterized in that, Includes the following steps: S1. Obtain at least two simulation snapshots generated during the system potential energy convergence phase of the molecular dynamics simulation, and divide the simulation box space into a three-dimensional grid region with a preset side length. S2. For each of the three-dimensional grid regions, extract the polymer chain segments in the grid region that are at a distance greater than a preset cutoff radius from the surface of the filler phase as the background reference system, and calculate the first displacement vector of the background reference system in the inter-frame time interval based on the at least two simulated snapshots. S3. For the filler phase within the three-dimensional grid region, calculate the second displacement vector of the filler phase within the inter-frame time interval based on the at least two simulated snapshots; S4. Based on the difference between the first displacement vector and the second displacement vector, calculate the local stability index, which characterizes the degree of decoupling between the filler phase and the bulk polymer chain segments within the three-dimensional grid region. S5. Calculate the local stability index of all three-dimensional grid regions, generate a global stability assessment result, and determine whether the packing dispersion is in a metastable state based on the global stability assessment result.

2. The simulation method for evaluating filler dispersion in polymer materials according to claim 1, characterized in that, In S1, the step of obtaining at least two simulation snapshots generated during the system potential energy convergence phase of the molecular dynamics simulation specifically includes: Read the original trajectory file output by the molecular dynamics simulation and extract the energy curve data of the total potential energy of the system as a function of time; In the energy curve data, identify the continuous time period in which the fluctuation amplitude of the total potential energy of the system is continuously less than a preset energy threshold, and mark this continuous time period as the potential energy convergence stage of the system. From the original trajectory file, extract all simulation snapshots whose timestamps are within the potential energy convergence phase of the system to form a candidate snapshot set; From the candidate snapshot set, at least two simulated snapshots whose time interval satisfies the preset inter-frame interval condition are selected as the at least two simulated snapshots output.

3. The simulation method for evaluating filler dispersion in polymer materials according to claim 1, characterized in that, In step S1, the step of dividing the simulated box space into three-dimensional mesh regions with preset side lengths specifically includes: Read any one of the at least two simulated snapshots and obtain the box length in the three dimensions of the simulated box space in that simulated snapshot. Obtain the preset grid side length parameters, calculate the ratio of the box length to the grid side length parameters in the three dimensions, and round each ratio up to obtain the number of grid divisions in the three dimensions. Based on the number of mesh divisions in the three dimensions, the simulation box is divided into multiple three-dimensional mesh regions by equal intervals in the three dimensions. Assign a unique spatial index identifier to each of the three-dimensional mesh regions.

4. The simulation method for evaluating filler dispersion in polymer materials according to claim 1, characterized in that, In step S2, the step of extracting polymer chain segments within the grid region that are at a distance greater than a preset cutoff radius from the filler phase surface as the background reference specifically includes: For the current three-dimensional mesh region, obtain the monomer coordinates of all polymer chain segments within the three-dimensional mesh region, and obtain the surface atomic coordinates of all filler phase atoms within the three-dimensional mesh region; For each monomer of the polymer chain segment, calculate the shortest spatial distance from the monomer to the surface atoms of the filler phase; The shortest spatial distance is compared with the preset cutoff radius, and individual units whose shortest spatial distance is greater than the preset cutoff radius are selected to form a candidate individual set; Determine whether the proportion of the number of monomers in the candidate monomer set to the total number of polymer chain segment monomers in the three-dimensional grid region is greater than a preset proportion threshold. If the ratio is greater than the preset threshold, the polymer chain segment to which the candidate monomer set belongs is extracted as the background reference system of the three-dimensional grid region; If the ratio is less than or equal to the preset ratio threshold, then the parent coarse-grained grid region to which the three-dimensional grid region belongs is retrieved, and the polymer chain segments in the parent coarse-grained grid region that satisfy the condition that the shortest spatial distance is greater than the preset cutoff radius are extracted as the background reference system of the three-dimensional grid region.

5. A simulation method for evaluating filler dispersion in polymer materials according to claim 1, characterized in that, In S2, the step of calculating the first displacement vector of the background reference frame within the inter-frame time interval based on the at least two simulated snapshots specifically includes: For the current 3D mesh region, determine the first simulation snapshot with the earlier timestamp and the second simulation snapshot with the later timestamp from the at least two simulation snapshots; Obtain the first centroid coordinates of the background reference system in the first frame simulation snapshot, and obtain the second centroid coordinates of the background reference system in the second frame simulation snapshot; Calculate the spatial vector difference between the second centroid coordinates and the first centroid coordinates, and determine the spatial vector difference as the first displacement vector of the background reference system within the inter-frame time interval.

6. The simulation method for evaluating filler dispersion in polymer materials according to claim 1, characterized in that, In S3, the step of calculating the second displacement vector of the filler phase within the inter-frame time interval based on the at least two simulated snapshots specifically includes: For the filler phase within the current three-dimensional mesh region, calculate the radius of gyration tensor of the filler phase, and determine the anisotropy of the filler phase based on the eigenvalues ​​of the radius of gyration tensor. Determine whether the anisotropy degree is greater than a preset anisotropy degree threshold; If the anisotropy is less than or equal to the preset anisotropy threshold, then the centroid displacement vector of the filler phase within the inter-frame time interval is obtained, and the centroid displacement vector is determined as the second displacement vector. If the anisotropy is greater than the preset anisotropy threshold, the overall displacement field of the filler phase during the inter-frame time interval is calculated, the centroid translation component is decomposed from the overall displacement field, and the centroid translation component is determined as the second displacement vector.

7. The simulation method for evaluating filler dispersion in polymer materials according to claim 1, characterized in that, In step S4, the step of calculating the local stability index, which characterizes the degree of decoupling between the filler phase and the bulk polymer chain segments within the three-dimensional grid region, specifically includes: For the current three-dimensional mesh region, calculate the spatial vector difference between the second displacement vector and the first displacement vector to obtain the difference displacement vector; Calculate the vector magnitude of the difference displacement vector, and calculate the vector magnitude of the first displacement vector; Calculate the ratio of the vector magnitude of the difference displacement vector to the vector magnitude of the first displacement vector, and determine this ratio as the local stability index of the three-dimensional mesh region.

8. A simulation method for evaluating filler dispersion in polymer materials according to claim 1, characterized in that, In step S5, the step of generating the global stability assessment result specifically includes: Collect all the calculated local stability indices of the three-dimensional mesh regions to form a set of local stability indices; Calculate the statistical characteristic values ​​of the set of local stability indices; The statistical characteristic values ​​are output as the global stability assessment results.

9. A simulation method for evaluating filler dispersion in polymer materials according to claim 1, characterized in that, In step S5, the step of determining whether the filler dispersion is in a metastable state based on the global stability assessment result specifically includes: Obtain the global stability assessment result and compare the global stability assessment result with a preset stability threshold; If the global stability assessment result is less than the preset stability threshold, the filler dispersion is determined to be in a true steady state, and a first assessment conclusion characterizing the dispersion stability is output. If the global stability assessment result is greater than or equal to the preset stability threshold, the filler dispersion is determined to be in a metastable state, and a second assessment conclusion characterizing the dispersion instability is output.

Images