Intelligent control method for blow molding process parameters based on digital twinning

Abstract

The present application relates to a digital twin-based intelligent regulation method for blow molding process parameters, aiming at the problem that the phase change state of the melt in the traditional process is difficult to identify in time and space with high precision, a data-driven digital twin modeling method based on array sensor space-time anchor point identification, phase change event intelligent screening and multi-dimensional dynamic template reverse locking is proposed. The method accurately deploys high-response array sensors, real-time collects and maps heat flow signals, relies on multi-stage filtering and feature extraction algorithm to identify effective phase change anchor points, realizes local adaptive correction of material dynamic parameters; and with the help of anchor point confidence dynamic evaluation and redundant point cooperative mechanism, the time and space alignment accuracy and robustness of virtual and real models are improved. The invention can significantly improve the melt phase change synchronous monitoring, parameter tracing and virtual-real integrated control capability of the blow molding process, and is conducive to improving the intelligent manufacturing level.

Description

Technical Field

[0001] This invention relates to the field of intelligent monitoring of advanced manufacturing processes and physical data-driven digital twin modeling technology, and in particular to a method for intelligent control of blow molding process parameters based on digital twins. Background Technology

[0002] Currently, with the rapid development of intelligent manufacturing and the Industrial Internet, digital twin technology has received widespread attention in the fields of intelligent monitoring and parameter optimization of manufacturing processes. Especially in high-temperature, high-speed molding processes such as blow molding, digital twin models are becoming a crucial technological foundation for achieving real-time equipment status perception, adaptive process control, and product quality feedforward control. Most existing digital twin models adopt an architecture combining physical-based modeling and data-driven correction. Their technological development trends include: leveraging real-time data feedback from multiple sensors, periodic batch training to assist model recalibration, utilizing surrogate models for multi-scale state prediction, and introducing adaptive algorithms such as reinforcement learning for dynamic process parameter tuning. The industry generally considers improving the dynamic state synchronization capability of digital twin models, reducing modeling and update latency, and enhancing the model's adaptability to changing process conditions as key research directions.

[0003] Currently, most representative digital twin dynamic calibration technologies rely on large-capacity sensor data within a preset acquisition period for global error identification and surrogate model retraining. For example, by periodically collecting full data on plastic melt temperature, pressure, and flow rate distribution, a deep learning model is used to offline train the state mapping between the virtual twin and the physical production line; when the error exceeds a threshold, a full model parameter reassessment is triggered periodically or in batches; or an online reinforcement learning algorithm based on deviation feedback is used to adaptively adjust key process parameters. These technologies are mostly used in typical manufacturing processes such as thermo-fluid-structure interaction and complex deformation field dynamics, and are effective for application scenarios with stable processes and well-defined parameter variation patterns. They are widely used in digital twin modeling and process optimization of continuous molding equipment such as injection molding, extrusion, and blow molding.

[0004] However, such technologies generally suffer from the following prominent problems and technical shortcomings: First, the periodic full data backhaul and global proxy model retraining bring huge communication load and computational pressure. Especially in the scenario where resources are limited at the edge of the actual production line, it is difficult to achieve millisecond-level model updates and dynamic adaptation of high-frequency parameters. The model state synchronization often has significant delays, which affects the real-time performance of process control.

[0005] Secondly, current dynamic calibration mechanisms based on deviation thresholds or periodic triggers often fail to accurately match the high-frequency, high-disturbance blow molding process. When physical process parameters are frequently fine-tuned or cross-batch process drift occurs, the synchronization between the virtual model and the physical equipment on the time axis and physical state axis rapidly decreases, leading to a decline in model prediction accuracy, affecting process consistency and product stability.

[0006] Furthermore, while general data-driven paradigms (such as reinforcement learning and global retraining) possess a certain degree of adaptability, their lack of deep coupling with the intrinsic phase transition physical processes of blow molding results in a lack of strong physical interpretability and regional specificity in the model update mechanism. This makes it impossible to drive efficient local calibration of the virtual model solely for key physical events, and it is difficult to balance state synchronization and computational resource consumption while maintaining extremely low update latency.

[0007] In the actual complex blow molding process environment, parameters such as material properties, mold structure, and cooling rate change significantly. If the update delay of the virtual twin model is too large or the global synchronization effect is lagging, it will inevitably cause the adjustment of process parameters to be hindered, the product size and performance to be unstable, and affect production consistency and production line yield.

[0008] Therefore, there is an urgent need for an innovative method for dynamically updating digital twin models. This method should effectively break through the existing general technical paradigms that rely on periodic full data and suffer from large global model synchronization lags. It should fully utilize the unique melt phase transition critical events in the blow molding process as dynamic anchor points, achieving high-precision, localized, and millisecond-level updates of the virtual model to the physical entity's state using only a low-overhead spatiotemporal anchoring mechanism. This ensures efficient response under limited resources at the edge and can adaptively handle complex operating conditions such as batch switching and environmental fluctuations. This will provide key technical support for intelligent blow molding and phase transition-oriented physical data-driven digital twin modeling, significantly improving process consistency and final product quality. It will also provide a practical and feasible technical solution to address industry pain points such as the insufficient process adaptability of existing intelligent manufacturing digital twin models. Summary of the Invention

[0009] This application provides a method for intelligent control of blow molding process parameters based on digital twins, which aims to solve one of the problems or issues of the prior art mentioned in the background.

[0010] The intelligent control method for blow molding process parameters based on digital twins provided in this application specifically includes: S1: Acquire the heat flow signals collected by the array sensors and record the spatial coordinate labels corresponding to each sensor to form a temperature field dataset.

[0011] S2: Perform time-domain differentiation on the temperature field dataset to extract the heat flow inflection point features caused by the jump in specific heat capacity and the abrupt change in thermal conductivity, and generate a candidate event set.

[0012] S3: Based on the confidence level of each event in the candidate event set, select valid anchor events, define the occurrence time of the valid anchor events as timestamp anchors and the occurrence location as spatial coordinate anchors, and construct spatiotemporal anchors.

[0013] S4: Match the corresponding target template from the preset template library based on the material properties and working parameters in the spatiotemporal anchor point.

[0014] S5: Using the spatiotemporal anchor point as a constraint, reverse locking is performed on the simulation step size, mesh node displacement, and material constitutive parameter combination in the target template to generate parameter remapping results.

[0015] S6: Based on the time offset standard deviation of phase transition events at the same location within the subsequent three consecutive blow molding cycles, calculate the weight coefficient of the spatiotemporal anchor point and construct the decay function to generate a weighted anchor point confidence index.

[0016] S7: If the weighted anchor confidence index is lower than the preset threshold, the collaborative verification process of adjacent redundant measurement points is initiated; otherwise, a simplified rheological model is used in conjunction with the parameter remapping results to perform interpolation and deduction, generating a virtual process state sequence to regulate the blow molding process parameters.

[0017] S8: Compare the virtual process state sequence with the real-time running state of the physical entity to identify the deviation, and update the internal state variables of the digital twin module based on the comparison results to complete the dual adaptive alignment of the virtual blow molding process and the physical process on the time axis and state axis.

[0018] The intelligent control method for blow molding process parameters based on digital twins provided in this application has the following beneficial effects: (1) By deploying array sensors in the key temperature measurement area of ​​the blow molding equipment mold, the transient heat flow inflection point signal caused by the sudden change in specific heat capacity and thermal conductivity when the melt front passes through is accurately captured. This physical phenomenon is identified as a "phase change anchoring event", and a unique spatiotemporal anchor point is constructed with its timestamp and spatial coordinates. This effectively overcomes the response lag and synchronization inaccuracy problems caused by the reliance on periodic batch data sampling or offline training of proxy models in traditional digital twin systems. Compared with the high computational overhead mechanism of using deviation threshold to trigger global relearning or reinforcement learning online parameter tuning in the existing technology, this solution uses the intrinsic phase change behavior of materials as an endogenous synchronization beacon, which significantly improves the real-time performance and physical consistency of virtual-real mapping. Under the premise of not needing frequent full model reconstruction, millisecond-level model state convergence is achieved, which greatly reduces the computational load and communication latency at the edge end. The overall model update latency is stably controlled within 1.5% of a single blow molding cycle, which effectively ensures the dynamic tracking capability and operational stability of the digital twin system under complex working conditions.

[0019] (2) A pre-set library of phase change dynamics templates matching different material grades, wall thickness ranges and cooling rates is introduced. After a real phase change anchoring event is detected, the corresponding template is automatically activated. The simulation step size, mesh node displacement and material constitutive parameter combination of the melt front in the virtual model are locked in reverse with the time and space information of the actual event as constraints. Only the local subdomain directly affected by the anchor point is incrementally remapped, avoiding the high resource consumption and computational redundancy caused by global model reconstruction in traditional methods. At the same time, combined with the "attenuation function", the weight of the phase change event at the same position in the subsequent multiple blow molding cycles is dynamically adjusted according to the time offset standard deviation. When the time offset is found to be increasing, the weight is automatically reduced and the collaborative verification mechanism of the adjacent redundant measurement points is started. This enhances the robustness of the system to uncertain factors such as process drift and sensor degradation, significantly improves the adaptive accuracy and reliability of model parameters in cross-batch production, and solves the model mismatch problem caused by environmental disturbance or material fluctuation in long-term operation.

[0020] (3) A simplified rheological model is used for interpolation between two phase change anchoring events. Its computational complexity decreases linearly with the anchor point interval, further optimizing the system resource utilization rate and ensuring that a high level of simulation continuity and process integrity can still be maintained under limited computing power. This design deeply couples the physical essence of the strong coupling between thermo-mechanical-phase change in the blow molding process, transforming the material phase change behavior into a dynamic synchronization mechanism with low communication overhead and high interpretability. It gets rid of the dependence on general artificial intelligence paradigms such as reinforcement learning or big data-driven agent modeling, and builds a lightweight, highly cohesive, and self-evolving digital twin architecture for specific manufacturing scenarios. It not only improves the flexibility and scalability of system deployment, but also provides efficient and reliable process monitoring and prediction capabilities for resource-constrained production line edge environments, and has good engineering implementation value and industrial application prospects.

[0021] The aforementioned technologies work together to construct a high-precision, low-latency, and robust digital twin synchronization framework for the blow molding process, driven by physical mechanisms. This framework enables a fundamental shift from passive response to active anchoring, from global recalculation to local correction, and from static modeling to dynamic adaptation, significantly improving the real-time performance, accuracy, and sustainable operation capabilities of the digital twin system in complex manufacturing scenarios. Attached Figure Description

[0022] Figure 1 This is the main flowchart of a method for intelligent control of blow molding process parameters based on digital twins.

[0023] Figure 2 This is a sub-flowchart of a method for intelligent control of blow molding process parameters based on digital twins.

[0024] Figure 3This is another sub-flowchart of the intelligent control method for blow molding process parameters based on digital twins. Detailed Implementation

[0025] Embodiments of the present invention are described in detail below, examples of which are shown in the accompanying drawings, wherein the same or similar reference numerals denote the same or similar elements or elements having the same or similar functions throughout. The embodiments described below with reference to the accompanying drawings are exemplary and are only used to explain the present invention, and should not be construed as limiting the present invention.

[0026] The following disclosure provides many different embodiments or examples for implementing different structures of the invention. To simplify the disclosure, specific examples of components and arrangements are described below. Of course, these are merely examples and are not intended to limit the invention. Furthermore, reference numerals and / or letters may be repeated in different examples; such repetition is for simplification and clarity and does not in itself indicate a relationship between the various embodiments and / or arrangements discussed.

[0027] like Figure 1 As shown, this application provides a method for intelligent control of blow molding process parameters based on digital twins, specifically including: S1: Acquire the heat flow signals collected by the array sensors and record the spatial coordinate labels corresponding to each sensor to form a temperature field dataset.

[0028] S2: Perform time-domain differentiation on the temperature field dataset to extract the heat flow inflection point features caused by the jump in specific heat capacity and the abrupt change in thermal conductivity, and generate a candidate event set.

[0029] S3: Based on the confidence level of each event in the candidate event set, select valid anchor events, define the occurrence time of the valid anchor events as timestamp anchors and the occurrence location as spatial coordinate anchors, and construct spatiotemporal anchors.

[0030] S4: Match the corresponding target template from the preset template library based on the material properties and working parameters in the spatiotemporal anchor point.

[0031] S5: Using the spatiotemporal anchor point as a constraint, reverse locking is performed on the simulation step size, mesh node displacement, and material constitutive parameter combination in the target template to generate parameter remapping results.

[0032] S6: Based on the time offset standard deviation of phase transition events at the same location within the subsequent three consecutive blow molding cycles, calculate the weight coefficient of the spatiotemporal anchor point and construct the decay function to generate a weighted anchor point confidence index.

[0033] S7: If the weighted anchor confidence index is lower than the preset threshold, the collaborative verification process of adjacent redundant measurement points is initiated; otherwise, a simplified rheological model is used in conjunction with the parameter remapping results to perform interpolation and deduction, generating a virtual process state sequence to regulate the blow molding process parameters.

[0034] S8: Compare the virtual process state sequence with the real-time running state of the physical entity to identify the deviation, and update the internal state variables of the digital twin module based on the comparison results to complete the dual adaptive alignment of the virtual blow molding process and the physical process on the time axis and state axis.

[0035] Step S1: Acquire the heat flow signals collected by the array sensors and record the spatial coordinate labels corresponding to each sensor to form a temperature field dataset. Specifically, this includes: S1.1: Based on the thermodynamic distribution characteristics of the blow molding mold cavity, the key temperature measurement area is determined, and an array of sensors is deployed in the key temperature measurement area according to a preset spatial topology to obtain a set of physical sensing nodes capable of capturing millisecond-level temperature fluctuations.

[0036] Based on the thermodynamic distribution characteristics inside the blow molding die cavity, a three-dimensional steady-state and transient thermal field distribution model is constructed using known process conditions and mold material thermal conductivity data. The region with the most significant temperature gradient change and sensitivity to phase transition critical points along the melt front's path during blow molding is identified as the key temperature measurement area. For this key temperature measurement area, node positions are selected according to a preset spatial topology structure. This structure is constructed based on the die cavity geometry and the main direction of heat flow distribution, ensuring sensor coverage at different flow branches and wall thickness variations. The installation feasibility of the selected node positions is analyzed, comprehensively considering the matching degree between the machining hole positions and the sensor response curves. Locations that meet the installation conditions and ensure stable heat signal transmission to the sensing surface are selected as the final deployment points. A sensor selection algorithm is used to select high-response micro-thermop models with response times less than the specified upper limit and sensitivity coefficients meeting the minimum resolution requirements from the sensor performance database, based on the requirement to capture millisecond-level temperature fluctuations. The sensor placement spacing is then calculated in conjunction with the spatial distribution of the node positions. Multi-objective optimization was performed on the deployment spacing and thermal field coverage. An iterative algorithm was used to adjust the local position of each sensor to eliminate signal redundancy or void areas, ensuring that the deployed physical sensor node set forms an equivalent continuous sampling field in space. By binding the sensor and node positions, the thermal field simulation results from the previous step were transformed into a list of sensor deployment coordinates and performance parameters, realizing the physical deployment output of the high-response micro thermopile array in the key temperature measurement area and the ability to capture millisecond-level temperature fluctuations.

[0037] For example, in the blow molding cavity of a daily chemical product packaging bottle, the wall thickness distribution exhibits a sharp change from the neck to the bottle body transition zone. Thermal field simulation results show that the maximum circumferential temperature gradient in this transition zone reaches 45℃ / mm, and the peak transient rate of change corresponds to a flow front arrival time of approximately 0.25s. Six nodes distributed circumferentially within this region were selected as key temperature measurement locations, with a node spacing of 12mm. Installation feasibility analysis indicates that four nodes can be directly embedded into the mold's inner wall through slots, while the remaining two nodes require installation via drilling holes in the outer wall. From the candidate sensor database, models with a response time not exceeding 50ms and a sensitivity coefficient of 0.018V / ℃ were selected. After optimizing the layout spacing and thermal field coverage, the positions of two edge nodes were adjusted to eliminate sampling gaps, ultimately forming a physical sensor node set with 97% coverage. The millisecond-level temperature fluctuation signals detected by this set during actual blow molding were used to construct the raw melt temperature field data for subsequent steps, achieving a significant improvement in phase transition critical point capture in the verification experiment.

[0038] S1.2: The array sensor is used to monitor the flow process of the plastic melt in the mold in real time, and the continuous heat flow signal caused by the advancement of the melt front is collected to generate raw voltage sequence data containing time dimension change characteristics.

[0039] Based on the array of sensors already deployed in the critical temperature measurement area of ​​the blow molding die, the sampling trigger signals of each sensor are synchronously activated to start the data acquisition channel through a hardware interrupt mechanism, ensuring that the millisecond-level time accuracy is not deviated due to software polling delay.

[0040] The analog voltage signal output by the sensor is amplitude modulated by a dedicated low-noise amplification module and quantized by a high-precision analog-to-digital converter. The sampling accuracy is set to 16 bits to ensure complete recording of the details of transient heat flow changes.

[0041] The quantized voltage sampling sequence is buffered and sorted according to a preset time synchronization flag to form a continuous time series matrix. The sorting process adopts a clock cycle alignment method to eliminate time drift between different sampling channels.

[0042] Dynamic range compression is performed on the data in each channel of the time series matrix to improve the stability of subsequent feature extraction.

[0043] The normalized time series data is stored as a raw voltage series data file and structured and encapsulated according to sensor ID and sampling time index to ensure that the data can be directly accessed by downstream processing modules.

[0044] Through the above-mentioned acquisition, quantization, synchronization, normalization and encapsulation processing methods, the set of physical sensor nodes obtained in the previous step is transformed into raw voltage sequence data containing time dimension change characteristics, so as to achieve high-precision and low-latency recording of heat flow signals caused by the advancement of the plastic melt front.

[0045] S1.3: Based on the installation geometric parameters of each sensor in the physical sensing node set, perform spatial coordinate calibration processing to assign a unique three-dimensional spatial coordinate label to each physical sensing node, so as to form a static spatial index library characterizing the physical location of the sensor.

[0046] For the set of physical sensor nodes that have been deployed and have acquired heat flow signals, the input conditions are the installation geometric parameters of each array sensor, including its radial distance, axial distance, and circumferential angle relative to the reference coordinate system of the blow molding mold cavity, as well as the origin and orientation definitions of the mold reference coordinate system. Based on the installation geometric parameters, a 3D spatial calibration algorithm is used to establish a coordinate mapping matrix from the mold reference coordinate system to the global digital twin model coordinate system, ensuring accurate correspondence between physical measurement points and simulation mesh nodes. During the calibration process, the actual installation deviation of each sensor is verified using 3D contour scanning results, obtaining the deviation vector and correcting the original geometric parameters to eliminate positional inconsistencies caused by manufacturing tolerances and installation errors. The corrected geometric parameters are then substituted into the coordinate transformation formula: in, This represents the sensor's three-dimensional coordinate vector in the global coordinate system. This represents the three-dimensional transformation matrix from the local coordinate system to the global coordinate system of the mold. This is the corrected 3D position vector of the sensor in the local coordinate system. This transformation operation is performed on each sensor, and the result is recorded, forming a unique 3D spatial coordinate label. The 3D spatial coordinate labels of all sensors are written into a static spatial index in numerical order, and the index is bound one-to-one with the set of physical sensor nodes. The coordinate labels enable spatial retrieval and association of subsequent data. Through spatial coordinate calibration, the physical acquisition nodes from the previous step are transformed into data indexes with unique 3D coordinate identifiers, achieving a precise correspondence between the raw data and the physical location, laying the foundation for the spatiotemporal fusion of subsequent temperature data.

[0047] For example, 12 high-response micro-thermopile sensors are deployed on a mold in a blow molding production line. The origin of the mold cavity reference coordinate system is located at the center of the bottom surface, with the axial direction along the height of the bottle. The radial and circumferential parameters are obtained through three-dimensional scanning. The position of sensor #5 in the local coordinate system is measured as (radial distance 16.25mm, axial distance 85.40mm, circumferential angle 135°). The manufacturing tolerance causes an actual position deviation vector of (-0.15mm, +0.20mm, -0.5°), and the corrected position vector is (16.10mm, 85.60mm, 134.5°). The transformation matrix T from the mold to the global coordinate system is calibrated and measured as: [0.866, 0.500,0.000;0.500,0.866,0.000;0.000,0.000,1.000] Convert the corrected local position vector to rectangular coordinates (mm) and substitute them into the formula to obtain P=( 10.31mm, (18.20mm, 85.60mm) were used as 3D spatial coordinate labels and stored in the index library. Verification showed that the calibration error was controlled within ±0.05mm, ensuring the high-precision matching capability of the static spatial index library in subsequent spatiotemporal fusion.

[0048] S1.4: Perform thermoelectric conversion calculations on the original voltage sequence data, and map the voltage amplitude to temperature values ​​based on the sensitivity coefficient of the micro thermopile sensor to generate a time-series temperature data stream with physical meaning.

[0049] S1.5: Based on the static spatial index library, the time series temperature data stream is spatiotemporally fused and bound with the corresponding three-dimensional spatial coordinate labels to construct a temperature field dataset containing complete time series information and location information, which serves as the sole data input object for subsequent phase transition feature extraction.

[0050] Step S2: Perform time-domain differentiation on the temperature field dataset to extract the heat flow inflection point features caused by the jump in specific heat capacity and abrupt change in thermal conductivity, and generate a candidate event set. Specifically, this includes: S2.1: Perform sliding window denoising processing on the time series data of each channel in the temperature field dataset, and use the adaptive Kalman filter algorithm to remove high-frequency electromagnetic interference and sensor thermal noise to generate a clean melt temperature time series signal with smooth characteristics.

[0051] It should be noted that a phase transition event refers to the instantaneous process during which a melt releases latent heat due to crystallization or glass transition during cooling, resulting in an inflection point in the heat flux signal; that is, the physical event corresponding to the heat flux inflection point characteristic. Each candidate event is assigned a confidence level, which comprehensively reflects three dimensions: signal-to-noise ratio (the ratio of the inflection point amplitude to the background noise), temporal consistency (whether the rate of temperature change before and after the inflection point conforms to the same physical process), and the degree of matching with the standard phase transition kinetic template. The higher the confidence level of a candidate event, the more likely it is to be a true phase transition critical point.

[0052] For the time series signals of each channel in the temperature field dataset, a sliding data truncation mechanism with a fixed number of sampling points as the window length is established, and the sequence in each window is used as an independent data segment to input the denoising processing algorithm.

[0053] The extracted data segments are input into the state estimation model of the adaptive Kalman filter. During the filter initialization phase, the initial values ​​of the state transition matrix and the observation matrix are set, and the process noise covariance matrix and the measurement noise covariance matrix are configured according to the noise variance parameters calibrated by the sensor.

[0054] By using a sliding window to update the estimated values ​​of process noise and measurement noise in real time, and by recursively calculating the Kalman gain coefficient, the filter can be adaptively adjusted to the signal at different time periods, thereby enhancing the ability to suppress non-stationary noise.

[0055] Within each window, iterative calculations of state prediction and observation update are performed, and the filter weights are calculated using the following Kalman gain formula: in, To predict the state covariance matrix, For the observation matrix, To measure the noise covariance matrix, For matrix transpose, For Kalman gain.

[0056] The observed and predicted values ​​of each sampling point within the window are weighted and fused according to the Kalman gain coefficient to generate a smoothed temperature signal data segment, which is then spliced ​​together during the sliding process to form a complete clean melt temperature time sequence signal.

[0057] By using adaptive Kalman filtering and sliding window segmentation, the temperature field data sequence containing high-frequency electromagnetic interference and sensor thermal noise is transformed into a clean temperature time series signal with continuity and smoothness, thereby improving the stability and accuracy of subsequent phase transition characteristic calculations.

[0058] For example, five micro-thermopil sensor nodes are arranged in the wall thickness variation area of ​​the blow molding mold cavity. The sampling frequency is set to 200Hz, the window length is selected as 40 sampling points, the initial value of the process noise covariance matrix is ​​set to 0.005, and the initial value of the measurement noise covariance matrix is ​​set to 0.01. In the initial filtering stage, the Kalman gain is calculated according to the above formula. Dynamically adjusts to the window size as it slides. Within each window, the peak fluctuation of the temperature signal is reduced to 40% of the original peak value after filtering, the signal curve shows a smooth trend, and the signal-to-noise ratio is significantly improved compared to the initial state. The clean temperature time series signal generated after complete sliding window processing can stably output temperature gradient and heat flux curvature characteristics in subsequent differential operations, ensuring a significant improvement in the accurate identification capability of phase change events.

[0059] S2.2: Perform central differential operation based on the clean melt temperature time sequence signal to calculate the temperature change rate at adjacent sampling times and generate an instantaneous temperature gradient sequence characterizing the melt heating or cooling rate.

[0060] Based on the clean melt temperature time-series signal obtained through adaptive Kalman filtering, a central difference operation is performed on the continuously sampled temperature data of each channel to form a numerical difference sequence that reflects the temperature changes at adjacent times.

[0061] The temperature value at each sampling point in the clean temperature time-series signal is calculated using a functional difference between the temperature values ​​at its immediate and adjacent sampling points, employing the central difference formula. Perform the calculation, where and These are the temperature values ​​at the previous and next sampling points, respectively. The sampling interval time. This represents the instantaneous temperature gradient.

[0062] The calculated instantaneous rate of change values ​​are arranged into a temperature rate of change vector in channel order, and a corresponding timestamp label is retained in each vector element to ensure time series traceability.

[0063] Outlier removal is performed on the temperature change rate vector. By setting upper and lower thresholds for the change rate, rate spikes that are physically impossible are filtered out to ensure the physical rationality of the gradient data.

[0064] The filtered temperature change rate vector is mapped to the corresponding spatial coordinate index to form an instantaneous temperature gradient sequence with spatiotemporal binding characteristics, which serves as the input for subsequent second-order differentiation and phase transition inflection point extraction.

[0065] By using central differential operation and spatiotemporal mapping, the temperature time sequence signal of the clean melt from the previous step is transformed into an instantaneous temperature gradient sequence with clear physical meaning and spatial identification, thereby achieving precise quantification of the melting heating or cooling rate.

[0066] For example, in the production line of blow-molded bottles for daily chemical products, sampling interval time is used. Set to 0.05 seconds, the sensor recorded consecutive temperature values ​​of 210.3℃, 211.0℃, and 211.5℃. Using the second sample as the calculation center point, the data before and after were substituted into the central difference formula, yielding an instantaneous temperature gradient of 12℃ / second at that sampling moment. This gradient value was then bound to the timestamp at 0.05 seconds and associated with the spatial coordinates (35mm, 120mm, 5mm) of the sensor channel on the mold wall, forming a data record in the instantaneous temperature gradient sequence. During continuous sampling, this method stably outputs the gradient sequence, with the gradient curve exhibiting a clear peak in the phase transition region. This provides a highly physically reliable input for subsequent phase transition inflection point detection, significantly improving the accuracy of the heating and cooling rates during dynamic model updates.

[0067] S2.3: Perform second-order time-domain differential processing on the instantaneous temperature gradient sequence, extract the acceleration component of the temperature change rate, and generate a transient heat flux curvature feature vector that reflects the intensity of the jump in specific heat capacity and the abrupt change in thermal conductivity.

[0068] For the instantaneous temperature gradient sequence data generated by the preceding sub-step S2.2, a second-order time-domain differential operation is performed for each time sampling point. The central difference method is used to construct the second-order derivative estimation formula. The temperature gradient difference results of adjacent sampling points are used as input to form the acceleration component sequence of the temperature change rate.

[0069] The formula for estimating the second derivative is defined as follows: in, Represents the instantaneous temperature gradient. The sampling time interval, For transient heat flux curvature, the formula uses [missing information] in the denominator. Power operations ensure the correct dimensionality of physical quantities.

[0070] The obtained temperature change rate acceleration component is normalized to eliminate the dimensional differences between different sensing nodes, and the normalization coefficient is set to be calculated based on the maximum rate of change amplitude of the historical temperature curve of each node.

[0071] Curvature calculation is performed on the normalized acceleration component sequence. The temperature change rate acceleration is regarded as a second-order geometric property of the curve. The curvature intensity is quantified by obtaining the absolute value, and a transient heat flow curvature feature vector is formed.

[0072] The curvature feature vectors are sorted by time index and mapped to three-dimensional spatial coordinate labels to ensure that subsequent zero-crossing detection and extreme value search can be located to the accurate spatial position, thereby realizing the physical characterization of the intensity of thermal capacity jump and thermal conductivity abrupt change.

[0073] By using the above-mentioned second-order time-domain differentiation and curvature solution processing method, the instantaneous temperature gradient sequence of the previous step is transformed into a transient heat flow curvature feature vector with physical meaning, thereby improving the ability to quantitatively analyze the critical transient phase change during blow molding.

[0074] For example, in a blow molding machine for daily chemical product packaging bottles, the sampling frequency of the micro thermopile sensor is set to 2000Hz, and the instantaneous temperature gradient sequence is in °C / ms. The sampling time interval is 0.5ms. For three adjacent points with temperature gradient values ​​of 0.8℃ / ms, 1.1℃ / ms, and 0.9℃ / ms, respectively, substituting into the above second-order time-domain differential formula, the acceleration component value is obtained as 1.6℃ / ms². Comparing this value with the historical maximum rate of change amplitude of the node, 2.0℃ / ms², the normalization coefficient is 0.8, and the absolute value of the normalized curvature intensity value is 1.28. This curvature intensity value is mapped to the corresponding sensor spatial coordinate label (12.5mm, 8.0mm, 3.2mm), combined with the time index of 125ms, to form a curvature feature vector unit and stored in the feature library. In this scenario, the curvature intensity value is significantly higher than the critical curvature reference value of 0.9 for phase transition, which verifies that the melt is in the critical region of specific heat capacity jump and thermal conductivity abrupt change, ensuring the accuracy and sensitivity of subsequent phase transition event identification.

[0075] S2.4: Based on the transient heat flux curvature feature vector, execute the zero-crossing detection and extreme value search algorithm to identify abnormal fluctuation points where the curvature sign flips or the amplitude exceeds the preset phase transition threshold, and generate a preliminary phase transition inflection point index list containing potential phase transition moments.

[0076] A point-by-point scanning process is performed on the transient heat flux curvature feature vector, which reflects the intensity of the jump in specific heat capacity and the abrupt change in thermal conductivity, to establish a curvature symbol sequence and record the positions of symbol changes. The positions of symbol changes are input into a zero-crossing detection algorithm, which determines the directional critical point of curvature change based on symbol flipping. An extreme value search algorithm is performed on the curvature feature vector to calculate the curvature amplitude at each sampling point and compare it with a preset phase transition threshold. An abnormal fluctuation point set is established using the threshold determination results, containing all sampling points whose amplitude exceeds the preset phase transition threshold. A temporal intersection operation is performed on the zero-crossing detection results and the abnormal fluctuation point set to select sampling points that simultaneously satisfy both symbol flipping and amplitude exceeding the limit as potential phase transition moment candidates. The time index of the candidate sampling points is written into a preliminary phase transition inflection point index list for subsequent spatial coordinate and timestamp binding.

[0077] By combining zero-crossing detection and extreme value search, the curvature feature vector generated in the previous step is transformed into a preliminary phase transition inflection point index list containing potential phase transition moments, thereby achieving precise positioning of the critical phase transition state.

[0078] For example, in one blow molding cycle, the curvature feature vector has a length of 2000 sampling points. The zero-crossing detection algorithm sets the sign change judgment window to 3 sampling points, and the extreme value search sets the preset phase transition threshold to 0.85, with the amplitude unit being W / (m²·K). The curvature feature vector is scanned point by point. The sign reversal judgment occurs at points 453, 927, and 1624, corresponding to the critical events of the curvature change direction reversing. The extreme value search detects a peak amplitude of 1.12 between points 450 and 460, a peak amplitude of 0.97 between points 922 and 930, and a peak amplitude of 1.05 between points 1620 and 1628, all exceeding the threshold. The temporal intersection results retain points 453, 927, and 1624 as potential phase transition moment candidates and write them into the preliminary phase transition inflection point index list. The list is then traced back to its spatial coordinates in subsequent steps to pinpoint the actual physical location, significantly improving the reliability and accuracy of phase transition event identification.

[0079] S2.5: Based on the preliminary phase transition inflection point index list, backtrack the corresponding spatial coordinate labels and timestamp information, aggregate the discrete inflection points into independent data units with spatiotemporal attributes, and generate a set of candidate events representing the critical state of phase transition.

[0080] like Figure 2 As shown, step S3 involves: filtering out valid anchoring events based on the confidence level of each event in the candidate event set, defining the occurrence time of the valid anchoring event as a timestamp anchor point and the occurrence location as a spatial coordinate anchor point, and constructing a spatiotemporal anchor point. Specifically, this includes: S3.1: Based on the transient heat flux inflection point amplitude of each event in the candidate event set, calculate the ratio of the current event signal strength to the root mean square value of the background noise, and generate a signal-to-noise ratio index.

[0081] S3.2: Using the signal-to-noise ratio index as an input variable, execute a time-series consistency verification algorithm based on a sliding window, compare the temperature change rate trend within a preset time window before and after the current event, and generate a time-series consistency coefficient.

[0082] Based on the signal-to-noise ratio (SNR) calculated in S3.1, this SNR is selected as the input variable for the time-series consistency verification algorithm. Time-domain pattern analysis is performed on the temperature change rate sequence of a single event in the candidate event set. The temperature change rate curve of the current event is extracted from the temperature change rate curves within the time windows containing a preset number of sampling points before and after it, ensuring that the start and end boundaries of the extracted segments are strictly aligned on the time axis. A sliding window structure is used to perform point-by-point interpolation on the extracted segments, calculating the temperature difference for each corresponding sampling point. The squared values ​​of the temperature differences at each sampling point are then summed to form the overall deviation within the window. The following normalized correlation coefficient formula is used to quantify the overall deviation within the window: in, The rate of temperature change at the i-th sampling point within the current event window. The correlation coefficient is the temperature change rate of the corresponding sampling point within the preceding window. The reciprocal of this correlation coefficient is used as the initial value of the time-series consistency coefficient. This initial value is updated through recursive calculations within the sliding window, generating a time-series consistency curve covering both windows before and after the event. Linear smoothing is performed on the time-series consistency curve to eliminate occasional spike noise. The mean of the smoothed curve is then used as the output time-series consistency coefficient for the current event. Through a sliding window-based time-series consistency verification process, the signal-to-noise ratio index from the previous step is transformed into a time-series consistency coefficient reflecting the physical rationality of the event, achieving accurate matching and reliability quantification of candidate phase transition events in the time domain.

[0083] For example, in a test on a blow molding production line, the temperature sampling frequency was 5kHz, the preset sliding window length was 250 sampling points, and the signal-to-noise ratio (SNR) was calculated to be 36.2 using S3.1. Data from sampling points 112 to 362 were extracted from the temperature change rate sequences before and after the current event. After performing point-by-point interpolation, the sum of squares of the temperature difference was calculated to be 18.6, and the sum of squares of the temperature change rate before the current event was 26.3. Substituting these values ​​into the above formula yielded the correlation coefficient, the reciprocal of which was 1.413. This reciprocal was used as an initial value for recursion within the window range, resulting in a smoothed mean temporal consistency coefficient of 1.408. This value is higher than the system's preset physical reasonableness threshold of 1.35, indicating a significant improvement in the consistency of the event's temporal change pattern, making it a valid input parameter for subsequent multi-dimensional feature matching operations.

[0084] S3.3: Perform multi-dimensional feature matching operation based on the time sequence matching coefficient and the preset phase transition dynamics feature template, calculate the Euclidean distance between the current event feature vector and the standard phase transition feature vector, and generate a phase transition matching score.

[0085] The execution objects for feature vector matching based on temporal consistency coefficients and pre-set phase transition dynamics feature templates are event data units that have passed consistency verification in the candidate event set and the standard feature vector set in the phase transition feature template library. The curvature feature vector of each candidate event is normalized with its corresponding temporal consistency coefficient to form a dimensionless multidimensional physical feature vector. For the normalized multidimensional feature vector, sub-components corresponding to melt propagation velocity, surface crystallinity change rate, and local contraction strain gradient are extracted, and Euclidean distance is calculated between these sub-components and the standard template feature vector in the same dimensional subspace.

[0086] The Euclidean distance is compared with the template matching tolerance threshold, and the distance is inversely mapped to a matching score. All candidate events are ranked according to the matching score, and a sequence of matching indexes characterizing the accuracy of event categories is generated. Through the above multi-dimensional feature matching operation, the physical features extracted in the previous steps are quantitatively compared with the standard template to form scoring data that can be used for confidence screening, achieving precise quantification of matching judgment and event classification.

[0087] For example, on a blow molding production line for plastic bottles used in daily chemical products, the mold cooling rate is set to 15℃ / s, the material grade is PET1101, and after normalization, the melt propagation speed component of the candidate event feature vector is 0.78, the surface crystallinity change rate component is 0.65, and the local shrinkage strain gradient component is 0.52; the corresponding components of the standard template are 0.80, 0.60, and 0.50, respectively. The Euclidean distance is calculated, yielding D1=0.0583. The matching score is calculated to be 0.945, significantly higher than the preset confidence threshold of 0.90, thus identifying the event as a high-accuracy phase transition anchoring event, which is subsequently included in the valid event subset. Performing the above matching process ensures that the spatiotemporal synchronous anchoring benchmark of the virtual blow molding process has an accurate physical feature correspondence, significantly improving the correctness and stability of the model's dynamic updates.

[0088] S3.4: Perform a threshold filtering operation on the candidate event set based on the phase transition matching degree score, remove invalid events with scores lower than a preset confidence threshold, retain high-scoring events and mark them as valid anchoring events, and generate a cleaned subset of valid anchoring events.

[0089] Numerical analysis is performed on the phase transition matching degree score sequence obtained by multidimensional feature matching operation, and the matching degree score corresponding to each candidate event is extracted as the input parameter for the filtering operation.

[0090] A score comparison matrix is ​​constructed based on the input parameters of the filtering operation, and a preset threshold is encapsulated as a constant threshold vector. The vector is then compared element by element with the score comparison matrix to determine the candidate event index set that is below the threshold.

[0091] The corresponding spatiotemporal attribute records are located in the original candidate event set using the candidate event index set, and such events are marked as invalid to remove their data entries.

[0092] Generate a list of high-scoring events from the remaining candidate events that have not been removed, and add an effective anchor tag field to the event data structure to explicitly mark the high-scoring events.

[0093] The high-scoring events after being marked are sorted in a secondary order according to their timestamps and spatial coordinate distribution characteristics, forming an output structure of an effective anchored event subset that has been optimized by the index.

[0094] By using threshold filtering and labeling, the phase transition matching score result from the previous step is transformed into an effective subset of anchoring events with spatiotemporal consistency and high signal purity, thereby providing reliable data for subsequent spatiotemporal anchor point construction.

[0095] For example, in a blow molding production line for plastic bottles used in daily chemical products, the heat flow signal sequences collected by 12 micro-thermopil sensors deployed within the mold were matched using Euclidean distance feature matching, resulting in a matching score range of 0.45 to 0.92. During the screening process, a preset threshold was set to 0.68. A score comparison matrix was constructed, and each element was compared with the threshold vector, resulting in a candidate event index set of 4 events with scores below the threshold. After the elimination operation, the remaining 8 high-scoring events were retained and marked as valid anchoring events. The screening intensity ratio was calculated using the following formal formula: in, The total number of candidate events. The number of events below the threshold is used. Substituting the total number of candidate events (12) and the number of events below the threshold (4), the screening intensity ratio is 0.67, indicating a significant increase in the proportion of high-scoring events. The effective anchoring event subset output after secondary sorting is applied in the subsequent spatiotemporal synchronization steps. During the verification process, the time difference between the virtual model and the physical measurement points remained within ±0.05 seconds, and the spatial deviation was less than 0.8 mm, achieving a significant reduction in model update latency and a significant improvement in process consistency.

[0096] S3.5: Extract the occurrence time data and sensor spatial coordinate labels of each event in the effective anchored event subset, encapsulate the occurrence time as a timestamp anchor parameter and the spatial coordinates as a spatial coordinate anchor parameter, and combine them to generate the spatiotemporal anchor point for spatiotemporal synchronization.

[0097] Step S4: Based on the material properties and working parameters in the spatiotemporal anchor point, match the corresponding target template from the preset template library. Specifically, this includes: S4.1: Based on the material grade code and mold cooling rate value obtained from the spatiotemporal anchor point parsing, perform a multi-dimensional index retrieval operation on the pre-set phase transformation dynamics template library to generate an initial template candidate set containing candidate template IDs and corresponding working condition matching degrees.

[0098] It should be noted that the pre-set template library refers to the set of phase change dynamics templates stored in the digital twin system in advance, corresponding to different blow molding process stages (such as melt extrusion, blow molding, and cooling and shaping). Each template contains the material constitutive relationship, thermophysical parameters and simulation control conditions for the corresponding stage. It can be quickly called up according to the remapping results of local subdomain parameters to realize the rapid deduction of the virtual process.

[0099] Based on the material grade code and mold cooling rate value obtained from spatiotemporal anchor point parsing, the two are combined as input parameters for multidimensional retrieval to determine the index dimension required to perform matching retrieval operation in the template library.

[0100] For the material grade code dimension, the template library material property index table is called to extract candidate template index branches containing the corresponding material thermophysical parameters, phase transformation characteristic curves and melt rheological coefficients.

[0101] For the mold cooling rate dimension, the template library working condition parameter index table is called to extract candidate template index branches that include cooling rate range, cooling medium properties and cooling distribution uniformity.

[0102] Cartesian product operations are performed using the material grade code index branch and the cooling rate index branch to generate all possible template combination indices, and the matching degree is calculated based on the input parameters in the spatiotemporal anchor point.

[0103] The matching degree is calculated using the normalized difference ratio method, and the matching degree value is obtained through the following formula: in, For matching degree, These are the measured values ​​for the dimensions corresponding to the template to be matched. As a reference value in the spatiotemporal anchor point, This indicates taking the absolute value.

[0104] For the matching degree of multi-dimensional parameter combinations, a weighted summation is performed to generate an overall matching degree index. The weight allocation is determined based on the sensitivity ranking of each parameter to template selection.

[0105] The overall matching degree index is compared with the preset matching degree threshold. Template IDs that meet the threshold conditions are retained to form an initial template candidate set, and corresponding matching degree metadata is attached as input for the S4.2 call.

[0106] Through the above-mentioned multidimensional index retrieval and weight matching degree calculation processing method, the spatiotemporal anchor point of the previous step is transformed into an initial template candidate set containing candidate template IDs and working condition matching degrees, so as to realize the rapid and accurate positioning of phase transformation dynamic templates that conform to specific material properties and cooling conditions in the template library.

[0107] For example, in the blow molding process of plastic bottles for daily chemical products, the material grade code resolved from the spatiotemporal anchor is "PET-1100", and the mold cooling rate is 8.5 K / s. In the template library, 20 template index branches containing PET series grades are retrieved for the material dimension, each branch having different specific heat capacity and Tg values; for the cooling rate dimension, 10 template index branches covering cooling rates of 6-10 K / s are retrieved, including different cooling media and flow distributions. 200 combined indices are generated using Cartesian product, and the matching degree is calculated according to the above formula. The reference values ​​in the spatiotemporal anchor are the material specific heat capacity of 1.15 kJ / (kg·K) and the cooling rate of 8.5 K / s, while the measured values ​​for the corresponding dimensions of the template to be matched are the corresponding parameter values ​​in the template. The multidimensional matching degree weights were set to 0.6 for the material dimension and 0.4 for the cooling rate dimension. After weighted summation, the optimal matching degree of 0.92 was obtained, corresponding to the template ID T-PET-1100C8, which was retained in the initial template candidate set. The output set contains 20 template IDs and matching degree metadata that meet the threshold of 0.85. These are used in subsequent steps to call the melt front advance velocity curves of these templates for time-series alignment verification.

[0108] S4.2: Using the candidate template IDs in the initial template candidate set, call the melt front advance velocity calibration curve data stored in the database, and perform time-series alignment calculation with the actual melt arrival time recorded by the spatiotemporal anchor point to generate a time-series deviation quantification index characterizing the degree of consistency in the time dimension.

[0109] Based on the candidate template IDs in the initial template candidate set, the melt front propulsion velocity calibration curve data stored in the database are retrieved to form the theoretical propulsion behavior dataset for the corresponding template. Preprocessing is performed on the actual melt arrival times recorded at spatiotemporal anchor points, unifying the sampling format and correcting the units of the raw time data to eliminate timescale differences arising across devices or batches. A time-series matching relationship is established between the theoretical propulsion behavior dataset and the actual melt arrival times, and alignment operations are used to synchronize the two time series according to the node positions corresponding to the spatial coordinate labels. Deviation calculation is performed on the synchronized time series using the following time-series deviation quantification formula: in, This is the actual melt arrival time series. This is the arrival time series corresponding to the theoretical propulsion speed curve. The number of spatial nodes participating in the comparison is specified. The calculated Δt index is standardized to ensure horizontal comparison under different batches and material properties. The standardized Δt index is output as a time-series deviation quantification index and stored in the template matching evaluation cache for subsequent template selection. By calling the database calibration curve and combining it with the actual time data of the unique spatiotemporal anchor point, time-series deviation quantification is performed to evaluate the temporal dimension consistency between theoretical propulsion characteristics and measured propulsion behavior.

[0110] For example, on a production line, there are a total of 20 spatial nodes within the mold cavity. The initial template candidate set contains 3 template IDs, and the corresponding database calibration curves have average propulsion speeds of 0.35 m / s, 0.38 m / s, and 0.40 m / s, respectively. The actual melt arrival time recorded by the spatiotemporal anchor points is standardized to milliseconds (ms), and the data range for the 20 spatial nodes is between 350 ms and 420 ms. After calculating the arrival time series of the theoretical propulsion speed curve for template ID=2, it is aligned with the actual arrival time series according to the node positions to obtain the difference sequence for each node. Substituting the difference sequence into the formula to calculate Δt, if the Σ difference is 1200 ms and the number of nodes n=20, then Δt is 60 ms. The Δt index is normalized according to the maximum and minimum difference intervals to obtain a normalized deviation index of 0.15. This deviation index is stored in the matching evaluation cache and is used to determine whether it meets the tolerance range in the subsequent screening stage. Through quantitative evaluation of this deviation, the template’s consistency in the time dimension is significantly improved, which meets the requirements of high-precision dynamic modeling.

[0111] S4.3: Select the preferred template subset that meets the preset tolerance range based on the time series deviation quantification index, and extract the surface crystallinity evolution curve parameters associated with each template in the preferred template subset. Combine the ambient temperature sensor readings of the current blow molding cycle to perform thermal history correction processing to generate an environmentally compensated surface crystallinity correction evolution curve.

[0112] Based on the time-series deviation quantification index output from step S4.2, each candidate template in the initial template candidate set is numerically filtered according to its deviation value and a preset tolerance range to form a preferred template subset that meets the time-matching accuracy requirements. For each template in the preferred template subset, its stored surface crystallinity evolution curve parameter set is read, including the initial crystallinity value, maximum crystallinity, crystallization rate constant, and curve shape factor. Based on the continuous temperature data collected by the ambient temperature sensor during the current blow molding cycle, the average ambient temperature and temperature fluctuation amplitude of the cycle are calculated as external input conditions for thermal history correction. A thermal history compensation function is introduced, and a parameterized nonlinear mapping method is used to calculate the correction amount of the ambient temperature change on the crystallization rate constant. The correction formula is as follows: in, The original crystallization rate constant of the template curve. For temperature sensitivity coefficient, The average ambient temperature for the current period. The reference ambient temperature is used for the material under template calibration conditions. The corrected rate constant replaces the original parameters of the template curve, and the entire curve is recalculated to generate an environmentally compensated surface crystallinity correction evolution curve. By incorporating temperature fluctuation amplitude, a curve smoothing filter is used to eliminate short-term high-frequency noise, thereby improving the stability of the correction curve in the digital twin model. Through this processing method, the optimal template selected by the time-series deviation quantification index is transformed into crystallinity evolution data that includes the influence of ambient temperature, achieving adaptive adjustment of the thermal history of the virtual blow molding process.

[0113] For example, for a blow-molded bottle made of PP-245 material, the average ambient temperature is kept stable at 22°C during one blow molding cycle, with temperature fluctuations not exceeding 0.5°C. Preferably, the crystallization rate constant k in the template subset is 0.018 s². -1 The reference ambient temperature Tref is 20℃, and the temperature sensitivity coefficient α is experimentally calibrated to be 0.0005s. -1 / ℃. Substituting the above parameters into the correction formula, the corrected rate constant k' is obtained as 0.019s. -1 The original rate constant was then replaced in the template, and the surface crystallinity correction curve was regenerated. During the curve generation process, a three-point weighted moving average filter was applied to eliminate instantaneous noise within 0.3℃. The resulting correction curve showed that the peak time of melt crystallinity appeared significantly earlier in the numerical simulation, the matching degree of the virtual blow molding model with the physical process was improved, the time deviation of melt state prediction was shortened to half of the original, and the process adaptability of the blow molding process was significantly improved.

[0114] S4.4: Based on the surface crystallinity correction evolution curve after environmental compensation and the spatial coordinate anchor point position information locked by spatiotemporal anchor points, the local contraction strain gradient sub-model in the preferred template subset is subjected to spatial grid node mapping verification to generate a local contraction strain gradient verification result with spatial topological consistency.

[0115] Based on the environmentally compensated surface crystallinity correction evolution curve and the spatial coordinate anchor point position information locked by the spatiotemporal anchor point, the grid node coordinate parameters of each local contraction strain gradient sub-model in the preferred template subset are extracted into two-dimensional or three-dimensional coordinate arrays. Spatial consistency preprocessing is then performed on these coordinate arrays to eliminate coordinate mapping deviations caused by differences in node order during template construction. For each time sampling point in the environmentally compensated surface crystallinity correction evolution curve, the corresponding spatial position parameters recorded in the spatiotemporal anchor point are called, and the expected contraction strain gradient value sequence at that location is calculated using a spatial interpolation function. This sequence is then used as a reference curve input for the verification process. For each local contraction strain gradient sub-model in the preferred template subset, a node-by-node comparison operation is performed based on the actual contraction strain gradient data sequence of its grid nodes and the reference curve. The deviation vector is calculated, and the mean square value of the deviation is quantified as a spatial mapping error index. The formula for calculating the mean square value is as follows: in, This represents the actual shrinkage strain gradient value. For reference curve values, The number of nodes is specified. A threshold judgment operation is performed on the spatial mapping error index, discarding template instances whose errors exceed the preset tolerance limit, and recording the spatial topology consistency verification results of the remaining templates. The node mapping relationship table, error index values, and reference curve fitting degree from the verification results are encapsulated into a structured output object for subsequent weighted scoring and fusion. Through the spatial grid node mapping verification process, the environmentally compensated surface crystallinity correction evolution curve is effectively mapped to the actual position locked by the unique spatiotemporal anchor point, achieving consistency verification of the local contraction strain gradient sub-model in the spatial topology.

[0116] S4.5: Combining the quantitative index of time-series deviation, the evolution curve of surface crystallinity correction, and the verification results of local shrinkage strain gradient, execute the weighted scoring fusion algorithm to determine the optimal matching object from the preferred template subset, so as to generate the final target template for reverse locking operation.

[0117] Using multi-source evaluation data from a subset of optimized templates as input, the temporal deviation quantification index, environmentally compensated surface crystallinity correction evolution curve parameters, and local shrinkage strain gradient verification results are extracted for each template. The temporal deviation quantification index is normalized to map its value range to a consistent range with other evaluation parameters for score fusion calculation. A time-adaptability scoring matrix is ​​constructed based on the normalized temporal deviation quantification index, and the fitting error value of the surface crystallinity correction evolution curve is converted into a crystallinity consistency scoring matrix according to preset evaluation weights. Using the spatial topology matching rate parameter from the local shrinkage strain gradient verification results, a strain matching scoring matrix is ​​constructed to ensure that the evaluation of the spatial structure consistency dimension has equal weight with the evaluation of the time adaptability and material consistency dimensions. A weighted scoring fusion algorithm is executed on the above three scoring matrices, using a multi-dimensional vector weighted summation method to calculate a single comprehensive score value, as shown in the following formula: in, Indicates the timing suitability score. Indicates the crystallinity consistency score. Indicates the strain matching score. , , The weighting coefficients are preset according to the importance of the process. The comprehensive scores are sorted in descending order, and each template in the preferred template subset is ranked. The template with the highest comprehensive score is selected as the final optimal matching object. Through the above weighted scoring fusion processing method, the multi-dimensional matching results of the previous step are transformed into a unique target template, achieving optimal template selection accuracy and suitability.

[0118] For example, during the blow molding cycle of PET material, the material grade code is "PET-GF20", the mold cooling rate is 15℃ / min, and the preferred template subset contains 3 candidate templates. The time series deviation quantification indices are 0.8, 0.6, and 0.7, respectively; the surface crystallinity correction evolution curve fitting error values ​​after environmental compensation are 0.05, 0.08, and 0.07, respectively; and the spatial topological matching rates corresponding to the strain matching score are 0.92, 0.85, and 0.88, respectively. With the weighting coefficients set to w_1=0.4, w_2=0.35, and w_3=0.25, the comprehensive score value is calculated by substituting into the formula. For the first template, the overall score is 0.4·0.8+0.35·0.95+0.25·0.92=0.86; for the second template, it is 0.4·0.6+0.35·0.92+0.25·0.85=0.78; and for the third template, it is 0.4·0.7+0.35·0.93+0.25·0.88=0.82. After ranking by overall score, the first template has the highest score, so the system automatically locks this template as the final target template and performs a reverse locking operation using this template in the subsequent S5 step. The verification results show that the deviation between the virtual blow molding process and the measured state is significantly reduced, the temperature field fitting error is reduced to within 0.03, and the model's state convergence speed and process adaptability are significantly improved.

[0119] like Figure 3 As shown, step S5 involves using the spatiotemporal anchor point as a constraint to perform reverse locking on the simulation step size, mesh node displacement, and material constitutive parameter combination in the target template, generating parameter remapping results. Specifically, this includes: S5.1: Based on the spatial coordinate anchor points contained in the spatiotemporal anchor points, locate the corresponding reference mesh node in the finite element mesh model of the target template, and calculate the set of neighboring mesh nodes directly covered by the thermal conduction and diffusion effect of the reference mesh node according to the preset thermal influence radius algorithm, and generate a local subdomain mesh topology structure containing the reference mesh node and neighboring mesh nodes.

[0120] It should be noted that the aforementioned reverse locking refers to using the measured time and spatial location of the phase transition provided by the spatiotemporal anchor point as boundary conditions to iteratively adjust the simulation parameters (including simulation step size, mesh node displacement, and material constitutive parameters) in the target template until the temperature field output by the template simulation is consistent with the measured value at that spatiotemporal anchor point. Specifically, if the measured time of arrival of the melt front at a certain spatial coordinate anchor point in the blow molding die cavity is t0, then the reverse locking process forces the finite element simulation to reach the phase transition temperature at that spatial coordinate anchor point at a time equal to t0, thereby inversely correcting the thermophysical parameters in the simulation model.

[0121] Based on the spatial coordinate anchor point data resolved from the spatiotemporal anchor points, the finite element mesh model index interface of the target template is invoked to locate the reference mesh node in the 3D mesh that perfectly matches the spatial coordinates, and a unique mesh number and physical property index for this node are generated. The thermal influence radius algorithm is numerically calculated for the located reference mesh node. The thermal diffusivity is determined based on the material's thermal conductivity, density, and current cooling rate, and the thermal conduction influence radius of the node in the current time step is calculated using Fourier's law of heat conduction. The thermal influence radius is compared with the spatial topology data of the finite element mesh, and all neighboring mesh nodes with a distance less than or equal to the thermal influence radius within the 3D Euclidean distance are selected, generating a spatial distance matrix for this set. The mesh topology analysis module is invoked to connect the reference mesh node and the set of neighboring mesh nodes into a local subdomain mesh topology structure, recording the connectivity between nodes, element boundary conditions, and heat flow transmission paths for each node within this structure. Through a geometric correction process, this local subdomain mesh topology structure is mapped one-to-one with the node indices in the finite element global mesh model to ensure that the subsequent local parameter remapping process can accurately locate the physically affected region. By using a processing method based on thermal influence radius calculation and mesh topology filtering, the spatial location of spatiotemporal anchor points is transformed into local subdomain mesh topology data containing the reference node and the set of neighboring nodes, thereby achieving the modeling accuracy control target limited to the direct influence range of the anchor points.

[0122] For example, in a blow molding production batch, the spatial coordinates of the spatiotemporal anchor points are (125.4mm, 87.2mm, 36.5mm), the material's thermal conductivity is 0.245W / (mm·K), its density is 1.38g / cm³, and the cooling rate is measured to be 0.75K / s. The reference mesh node numbered N_482 is located in the finite element mesh model. According to Fourier's law, the formula for calculating the thermal diffusivity α is: in For the thermal conductivity of the material, For material density, Given the specific heat capacity of the material, and substituting the value of 1.25 J / (g·K), we obtain the thermal diffusivity a≈0.142 mm² / s. At the current time step Δt=0.05 s, the thermally affected radius R is calculated as follows: The calculated result is R≈2.67mm. Based on the 3D Euclidean distance matrix of the local mesh, 16 neighboring nodes with a distance of less than or equal to 2.67mm from the reference node were selected, forming a local subdomain mesh topology structure of the reference node and the 16 neighboring nodes. The connectivity information and spatial coordinates between nodes were marked in this structure. After performing geometric correction, this local subdomain perfectly matches the global mesh index table. Subsequent calculations of local temperature residuals and reverse locking of material parameters can be accurately applied to this subdomain, achieving a significant improvement in the adaptability of the virtual model to local areas in this batch of production.

[0123] S5.2: For each grid node in the local subdomain grid topology, extract the predicted temperature field data at the current simulation step size, and perform a difference operation between the predicted temperature field data and the measured phase transition temperature data corresponding to the timestamp anchor point recorded in the spatiotemporal anchor point to generate a local temperature residual vector.

[0124] In the constructed local subdomain mesh topology, the reference mesh node and its neighboring mesh nodes are used as operation objects. The temperature field prediction data corresponding to the current simulation step size in the target template is called to form a node temperature prediction matrix. Based on the timestamp anchor point in the spatiotemporal anchor point, the measured phase transition temperature value at the corresponding position of the anchor point is retrieved, and this value is mapped to the local subdomain mesh node according to the node index to form the node measured temperature vector.

[0125] A difference operation is performed on the corresponding elements of the predicted values ​​and measured temperature vectors at each node in the temperature prediction matrix. The output of the difference operation is then reorganized into a local temperature residual vector according to the node spatial order and stored in a high-precision floating-point data object to ensure the numerical stability of the subsequent sensitivity matrix construction. Through the above difference operation, the state difference between the virtual model and the physical entity is quantified at the thermal field parameter level, realizing a structured characterization of the deviation.

[0126] For example, at an ambient temperature of 40℃, a simulation prediction calculation with a step size of 0.002s was performed on a blow-molded part made of PP-270 material. The local subdomain mesh contains one reference node and eight neighboring nodes. The target template predicted the node temperatures as follows: reference node 189.4℃, neighboring node ranges from 187.1℃ to 191.0℃. The spatiotemporal anchor point records the measured phase transition temperature of the reference node as 188.7℃, and the corresponding measured values ​​of the neighboring nodes range from 186.5℃ to 190.2℃. The difference was calculated sequentially, with ΔT of 0.7℃ for the reference node and ΔT values ​​of the neighboring nodes distributed as follows: Between 0.4℃ and 0.8℃, a local temperature residual vector of length 9 is formed. This vector is directly used as input in the subsequent construction of the sensitivity matrix of the adjoint variable method, ensuring that the correction direction during optimization iteration is consistent with the physical deviation, thereby significantly improving the correction accuracy of the local subdomain thermophysical parameters.

[0127] S5.3: Based on the local temperature residual vector, a sensitivity matrix for the material constitutive parameters is constructed using the adjoint variable method, and the material thermal property correction coefficients that can eliminate the local temperature residual vector are calculated in reverse iteration using the gradient descent optimization algorithm, generating an updated combination of material constitutive parameters including the corrected specific heat capacity parameter and the corrected thermal conductivity parameter.

[0128] Based on the numerical input of the local temperature residual vector, the initial values ​​of the material constitutive parameters in the target template are called as the calculation basis to establish the coupled system of the state equation and error equation required by the adjoint variable method.

[0129] The adjoint equation of the temperature field control equation after finite element discretization is constructed. The local temperature residual vector is used as the error weight coefficient. The influence of each material thermophysical parameter on the temperature field error at the current time step is compared to generate a complete material parameter sensitivity matrix.

[0130] For each element of the sensitivity matrix, gradient calculation is performed to map the partial derivative of the temperature field error to the influence coefficients of the two physical property parameters, specific heat capacity and thermal conductivity, which serve as the search direction benchmark for the optimization algorithm.

[0131] During gradient descent optimization, a learning rate and a convergence threshold are set, and the material constitutive parameters are iteratively updated. The update amount for a single iteration is calculated using the following formula: in, For material parameter update amount, For learning rate, This is the gradient vector constructed based on the sensitivity matrix.

[0132] The updated values ​​obtained from the iterative calculation are superimposed with the initial values ​​of the current material parameters to form a new combination of specific heat capacity and thermal conductivity parameters. This combination is then substituted into the adjoint equation to execute the next round of iterations until the residual norm satisfies the convergence condition.

[0133] In the convergence determination stage, the following formula for calculating the residual norm is used: in, To predict temperature, For actual measured temperature, This represents the number of nodes.

[0134] By combining the adjoint variable method and the gradient descent optimization algorithm, the local temperature residual vector from the previous step is transformed into an updated combination of material constitutive parameters, including the corrected specific heat capacity parameter and the corrected thermal conductivity parameter, thereby achieving accurate dynamic correction of the material's thermal properties within the local subdomain.

[0135] For example, during the blow molding cycle, the maximum absolute value of the local temperature residual vector is 3.2℃, the initial value of specific heat capacity is set to 2200 J / (kg·K), the initial value of thermal conductivity is set to 0.25 W / (m·K), the learning rate η is set to 0.015, and the convergence threshold is 0.001. Based on the temperature field of a finite element mesh with a total of 64 nodes, after constructing the sensitivity matrix, the maximum component of the gradient vector is calculated to be 45.8, corresponding to an update amount of the specific heat capacity parameter. The update amount for the thermal conductivity parameter is 0.687. After 8 iterations, the residual norm decreased to 0.0009, satisfying the convergence condition. The output updated material constitutive parameter combination is a specific heat capacity of 2199.313 J / (kg·K) and a thermal conductivity of 0.239 W / (m·K). In subsequent explicit integration of the local mesh, this combination significantly improved the fit between the virtual model temperature field and the physical entity temperature field, keeping the model update delay within 1.5% of a single cycle.

[0136] S5.4: Based on the updated material constitutive parameter combination, adaptively recalculate the explicit integration time step of each grid node in the local subdomain grid topology to determine the maximum allowable simulation step size that satisfies the numerical stability condition, and generate a local optimized simulation step size sequence that adapts to the current phase transition dynamics interval.

[0137] The input conditions include the updated material constitutive parameter combination output from step S5.3, and the current simulation step size and thermodynamic state data of each grid node in the local subdomain mesh topology. Based on the updated material constitutive parameter combination, the corrected specific heat capacity and thermal conductivity parameters are extracted as the physical input for time step recalculation. For each local grid node, the stability criterion in the explicit integration method is called to combine the corrected material parameters with the local temperature gradient values ​​on the node to calculate the thermal diffusivity, which characterizes the constraint of heat conduction velocity on the integration step size. The stability criterion formula is used: in, To the maximum allowable simulation step size, The stability coefficient, The thermal diffusivity is calculated based on the corrected specific heat capacity and thermal conductivity. The local grid node spacing is used. The thermal diffusivity in the above formula is obtained through... The calculation yielded, where This is the corrected thermal conductivity. For material density, The calculated specific heat capacity is used. This calculation is performed on all local nodes to generate a node-level maximum allowable simulation step size array. A global minimum selection operation is performed on the node-level step size array to ensure numerical stability of the entire local subdomain during the integral iteration process. The selection results are divided into stages according to the phase transition dynamics interval to generate a locally optimized simulation step size sequence, which is matched with the phase transition critical point time information to ensure that the time step automatically shrinks within the critical interval and moderately widens within the stable interval. Through the above processing method, the results of the previous step are transformed into adaptive time step size data that meets the premise of numerical stability, realizing the accurate evolution of the virtual blow molding model in the local subdomain.

[0138] For example, within a local subdomain of the mold shoulder of a blow-molded chemical bottle, the mesh node spacing... The corrected thermal conductivity is 0.002m. The material density is 0.22 W / (m·K). 970 kg / m³, specific heat capacity The thermal diffusivity is 1800 J / (kg·K). The calculation is 1.27 × 10 -7 When the speed is m² / s and the stability coefficient C is 0.5, the maximum allowable simulation step size is... 1.59×10 -9 The result is then globally compared with the calculated step sizes of other nodes. For this region, a step size of 1.59 × 10⁻⁶ is selected. -9 s is used as the critical interval step size, which is relaxed to 2.0 × 10⁻⁶ in the stationary interval. -9 s. It has been verified that the step size sequence significantly improves the simulation accuracy of local phase transition in the shoulder under numerical stability conditions and effectively controls the model update delay.

[0139] S5.5: Write the local optimization simulation step size sequence, the updated material constitutive parameter combination, and the grid node displacement after correction based on the measured position into the corresponding storage address of the target template to complete the local subdomain parameter remapping of the region directly affected by the spatiotemporal anchor point and generate the parameter remapping result.

[0140] Step S6: Based on the time offset standard deviation of phase transition events at the same location within three consecutive blow molding cycles, calculate the weighting coefficient of the spatiotemporal anchor point and construct a decay function to generate a weighted anchor point confidence index. Specifically, this includes: S6.1: Obtain the occurrence time data of valid anchoring events recorded at the same spatial coordinate anchor point within three consecutive blow molding cycles, and perform time-series alignment processing on the occurrence time data to generate an original time series set containing time offsets.

[0141] S6.2: Perform standard deviation calculation based on the time offset in the original time series set to quantify the dispersion of phase change events on the time axis and generate a phase change time offset standard deviation index that characterizes process stability.

[0142] S6.3: Using the phase transition time offset standard deviation index as an input variable, and substituting it into a preset nonlinear decay mapping function for numerical transformation processing, in order to establish an inverse correlation between time fluctuation and trust level, and generate initial anchor point weight coefficients.

[0143] The standard deviation of the phase transition time offset is obtained as a numerical input, and the parameter table of the preset nonlinear decay mapping function in the digital twin module is called.

[0144] The input indicators are numerically bound to the control variables of the function, and the validity interval conditions within the function's domain are checked to ensure the stability of the mapping calculation.

[0145] Under the premise of meeting the conditions, the index value is substituted into the core formula structure of the decay mapping function. This formula adopts an exponential decay model to establish an inverse correlation between time fluctuation and trust level.

[0146] in, The initial anchor point weight coefficients, For the theoretical maximum weight, It is a natural constant. The decay rate factor, This represents the standard deviation of the time offset for the phase transition event.

[0147] The exponential calculation unit is executed to multiply the attenuation factor by the input index, take the negative of the result, and then perform exponential evaluation to obtain the attenuation ratio.

[0148] Multiply the attenuation ratio by the theoretical maximum weight to calculate the corresponding initial anchor weight coefficient.

[0149] The calculation output is written to the weight register of the current spatiotemporal anchor point, which is then used by the weighted fusion algorithm.

[0150] By using an exponential nonlinear mapping function, time fluctuation indicators are transformed into quantified initial anchor weight coefficients, thereby enabling the construction of basic data for dynamic evaluation of the model anchor credibility.

[0151] For example, under a certain blow molding cycle condition, the standard deviation of the time offset of the phase transition event at the same position in three consecutive cycles is 0.85 seconds, the preset theoretical maximum weight is 1.0, and the decay rate factor is set to 0.9. Substituting this standard deviation into the formula, the exponential calculation part is... The calculated attenuation ratio is approximately 0.461. Multiplying this value by 1.0 yields an initial weight coefficient of 0.461. This coefficient is used as a core weight input in the subsequent weighted fusion operation in S6.4. After smoothing, it is applied to determine the credibility of the anchor point. In this scenario, the virtual blow molding process's dependence on this anchor point is significantly reduced, avoiding the risk of state mismatch caused by anchor points with large time fluctuations directly driving the model.

[0152] S6.4: Perform a weighted fusion operation based on the initial anchor point weight coefficient and the environmental interference factor of the current blow molding cycle to eliminate the disturbance of occasional noise to the evaluation results and generate a smoothed dynamic adjustment weight value.

[0153] The input conditions are the initial anchor point weight coefficient generated in the previous step S6.3 and the environmental interference factor data collected in real time during the current blow molding cycle, including multi-dimensional disturbance indicators such as ambient temperature fluctuation value, equipment vibration acceleration value, and instantaneous deviation of air supply pressure.

[0154] Normalization is performed based on the initial anchor point weight coefficient and the environmental disturbance factor vector to convert each disturbance component into a dimensionless standard value according to its own physical quantity dimension and maintain the consistency of the numerical range.

[0155] The normalized environmental disturbance factor vector is multiplied by the preset weight matrix and the initial anchor point weight coefficient to generate a weighted noise coefficient that characterizes the overall intensity of environmental disturbance.

[0156] An exponential smoothing filtering algorithm is used to smooth the weighted noise coefficients. In the filtering operation, the dynamic sampling interval of the current blow molding cycle is used as the time constant to adapt to the stability requirements of different cycles.

[0157] The difference correction operation is performed between the smoothed weighted noise coefficient and the initial anchor weight coefficient. Through this weighted fusion processing method, the initial anchor weight coefficient of the previous step is transformed into a dynamically adjusted weight value after environmental disturbance compensation, thereby suppressing the impact of environmental noise on the evaluation results and improving stability.

[0158] S6.5: Based on the dynamically adjusted weight value, encapsulate a structured data object containing timestamp identifiers and confidence values ​​to form a standardized evaluation result that can be called by subsequent verification processes, and generate a dynamically adjusted weighted anchor confidence index.

[0159] Step S7: If the weighted anchor point confidence index is lower than a preset threshold, a collaborative verification process for adjacent redundant measurement points is initiated; otherwise, a simplified rheological model combined with the parameter remapping results is used for interpolation and deduction to generate a virtual process state sequence for controlling blow molding process parameters. Specifically, this includes: S7.1: Obtain the weighted anchor confidence index and the preset threshold judgment result, and perform logical branch routing processing based on the judgment result. If the index is lower than the threshold, generate a nearby redundant measurement point collaborative anchoring verification instruction. If the index is higher than the threshold, generate a simplified rheological model interpolation and deduction enable signal to establish the virtual state update strategy of the current blow molding cycle.

[0160] It should be noted that the simplified rheological model is a lightweight thermo-rheological coupling computational model adapted to the real-time simulation requirements of digital twins. It takes the local subdomain parameters (including material constitutive parameters, thermophysical parameters, and mesh node displacements) corrected after inverse locking as input, and constructs a computational framework based on the phase change dynamics templates of the corresponding stages in the pre-set template library. The model uses the spatiotemporal anchor points corresponding to two effective phase change events as boundary constraints. Under the premise of keeping the basic control equations of melt rheology and heat transfer unchanged, it constructs a hybrid computational architecture of high-precision local simulation and low-latency global interpolation by performing explicit iteration only on the local subdomains affected by the anchor points and using field distribution interpolation for the remaining regions. This significantly reduces the real-time computational load of the digital twin system while ensuring the physical consistency of the evolution trend of the melt temperature field and rheological field.

[0161] It should be noted that the collaborative verification process refers to a closed-loop process that uses the aforementioned spatiotemporal anchor points as a benchmark to compare the consistency between the virtual state sequence of the digital twin model and the actual state data of the physical entity. It achieves dynamic correction and collaborative synchronization of virtual and real states through the linkage of confidence assessment, parameter remapping, and model interpolation. For example, the collaborative verification process first assesses the confidence of the collected spatiotemporal anchor points to select valid anchor points; then, it remaps the local subdomain parameters based on the valid anchor points; finally, it calls a simplified rheological model to generate a virtual state sequence and compares it with the next collected spatiotemporal anchor points to complete the closed-loop verification.

[0162] S7.2: In response to the collaborative anchoring verification command of neighboring redundant measurement points, the heat flow signal sequence of neighboring redundant measurement points is obtained. The spatiotemporal correlation weighted fusion algorithm is used to perform consistency comparison and error compensation processing on multi-source phase change events, and a high-precision spatiotemporal anchor point identifier after collaborative correction is generated to eliminate the model mismatch risk caused by single-point measurement drift.

[0163] The heat flow signal sequence collected from nearby redundant measuring points is acquired and loaded into the multi-source phase change event comparison buffer. The spatiotemporal matching index is established using the three-dimensional spatial coordinates of each measuring point and the sampling timestamp. Each measuring point sequence is resampled to a uniform sampling interval according to the time axis to ensure the comparability of subsequent algorithms.

[0164] The correlation coefficient matrix is ​​calculated for multi-source signals based on the spatiotemporal matching index. A spatial correlation threshold is set and abnormal measurement points with low correlation are removed to retain data sources with high physical consistency, thereby reducing the impact of noise interference on anchor point correction.

[0165] A weighted fusion operation is performed on the retained data source signals. The weights are calculated using the following formula: in, For the first The weight of each neighboring measurement point The denominator is the spatial correlation coefficient, and the sum of the correlation coefficients of all retained measurement points is used to ensure weight normalization.

[0166] The fused signal sequence is compared with the current low-confidence anchor point signal. The degree of single-point drift is quantified by zero-crossing detection and amplitude deviation statistics. An error compensation vector is constructed and used to correct the timestamp and spatial coordinate parameters of the anchor point.

[0167] The compensated anchor point parameters are encapsulated into structured data objects and marked as high-precision spatiotemporal anchor point identifiers after collaborative correction. Through the above processing method, the results of the previous step are transformed into high-confidence anchor point data that can be called for subsequent virtual blow molding process state updates, thereby achieving the expected technical effect of eliminating model mismatch caused by single-point measurement drift.

[0168] For example, on a plastic bottle blow molding production line, one main measuring point and three adjacent redundant measuring points are arranged inside the mold cavity. The sampling frequency of the high-response micro-thermop sensor is 200Hz. The spatial correlation coefficients of the redundant measuring points are 0.92, 0.88, and 0.85, respectively. A spatial correlation threshold of 0.8 is set, and measuring points below the threshold are removed, retaining all three sets of data sources. According to the formula, the weights of the three measuring points are calculated to be 0.34, 0.32, and 0.31, respectively. The difference between the fused signal and the main measuring point signal shows a time offset of 12ms and an amplitude deviation of 0.5℃. By constructing an error compensation vector, the anchor point timestamp is corrected to the original value minus 12ms, and the spatial coordinates are corrected to the original coordinates plus a displacement of 0.2mm along the melt front direction. The final high-precision spatiotemporal anchor point identifier significantly improves the consistency between the virtual model and the physical process in subsequent state updates.

[0169] S7.3: Responding to the simplified rheological model interpolation enable signal, obtain the local subdomain parameter remapping results and the state boundary conditions of the preceding effective anchoring events, call the simplified rheological model to perform nonlinear state interpolation operations, and generate intermediate melt rheological field data between the two phase transition anchoring events to maintain the temporal continuity of the virtual process.

[0170] When responding to the simplified rheological model interpolation enable signal, the local subdomain parameter remapping result is used as the main input variable, and the state boundary conditions of the preceding effective anchoring events are analyzed to establish the temporal and spatial boundary constraints of the interpolation computational domain. A simplified rheological model running instance is retrieved, and the interpolation time step is set according to the time start and end points locked by the boundary conditions. The coefficients of the nonlinear rheological equation are determined based on the material thermophysical parameters contained in the remapping result. State extrapolation is performed at each moment within the interpolation computational domain. Using the finite volume method combined with explicit time integration, the state variables of each grid node under the dual coupling of temperature and velocity fields are predicted, forming a nonlinear transition sequence from the state of the adjacent anchor point to the state at the intermediate moment. A fidelity control factor is introduced to correct the property constraints of the predicted sequence, ensuring that the intermediate rheological field data generated by interpolation remains consistent under thermodynamic and hydrodynamic conditions. Through the above processing method, the local subdomain parameter remapping result of the previous step is transformed into intermediate melt rheological field data covering the period between two phase change anchoring events, achieving the technical effect of continuous evolution of the virtual process with controlled delay.

[0171] For example, in the blow molding process of daily chemical product packaging bottles, the local subdomain parameter remapping results include material correction parameters with a specific heat capacity of 1.95 kJ / (kg·K) and a thermal conductivity of 0.245 W / (m·K). The state boundary conditions recorded by the preceding effective anchoring event are the melt front position and temperature field distribution corresponding to the start time of 0.35 s and the end time of 0.50 s. The interpolation time step is set to 0.005 s. When calling the simplified rheological model, the equation coefficients are calculated based on the combination of specific heat capacity and thermal conductivity. The rheological equation is: ,in Thermal conductivity, For density, For specific heat capacity, This represents the first derivative of temperature with respect to time. In this example, assuming a density of 955 kg / m³, the temperature change rate at each step is obtained by substituting the parameters and used for time integration. The finite volume method is applied to update the state of each grid node between 0.35 s and 0.50 s, generating intermediate melt rheological data containing temperature and velocity fields. This data is then corrected for thermodynamic consistency by a fidelity control factor of 0.85. This output data is used for subsequent incremental grid state updates in S7.4. Verification results show that the difference between the state changes of the virtual process and the physical measurements is significantly reduced within this interval, successfully achieving temporal continuity and low-delay response in the virtual process.

[0172] S7.4: Receive high-precision spatiotemporal anchor point identifiers or intermediate melt rheological field data after collaborative correction, map the above data to digital twin grid nodes based on discrete time step synchronization mechanism, perform incremental update operation of state variables, and generate a virtual process state sequence containing complete spatiotemporal evolution characteristics, so as to realize the virtual model's millisecond-level following of physical entities.

[0173] The system receives high-precision spatiotemporal anchor point markers or intermediate melt rheological field data after co-calibration, parses the included timestamp parameters and 3D spatial coordinate labels, and locates the corresponding target mesh node set in the digital twin finite element mesh model. Based on the discrete time step synchronization mechanism, it calculates the time difference between the current frame time of the virtual simulation engine and the timestamp of the high-precision spatiotemporal anchor point marker, and performs index positioning on the simulation step sequence according to the time difference to identify the time step interval that needs to be mapped to the state. For the target mesh node set, it loads the melt rheological field state variables generated by co-calibration or interpolation, including temperature field distribution, viscosity field distribution, crystallinity distribution, and local strain distribution. Using a spatial interpolation matching algorithm, it performs consistency correction between the spatial resolution of the input state variables and the resolution of the digital twin mesh model to ensure that the node spacing and coordinate system of the state mapping are completely aligned. For the spatially consistent state variables, it performs an incremental update operation, directly accumulating the changes to the existing state values ​​of the digital twin mesh nodes, and writes the updated state values ​​into the simulation engine's runtime cache to ensure a smooth transition of the model state in terms of temporal continuity. By incrementally updating the grid node states, the collaborative correction data or interpolation results from the previous step are transformed into a virtual process state sequence containing complete temporal evolution and spatial coupling characteristics, enabling the virtual model to dynamically follow physical entities at the millisecond level.

[0174] Step S8: Compare the virtual process state sequence with the real-time operating state of the physical entity, and update the internal state variables of the digital twin module based on the comparison result, thus completing the dual adaptive alignment of the virtual blow molding process and the physical process on the time axis and state axis. Specifically, this includes: S8.1: Acquire the grid node displacement data and material constitutive parameter combination data in the virtual process state sequence, and simultaneously collect the sensor measurement data of the corresponding spatiotemporal coordinates in the real-time operation state of the physical entity to form a paired state dataset containing virtual predicted values ​​and physical measured values.

[0175] S8.2: Perform multi-scale residual calculation processing based on the paired state dataset, and use a spatial interpolation algorithm to map the virtual predicted values ​​to the sampling grid of the physical measured values ​​to generate a state deviation vector field that represents the difference between the virtual blow molding process state sequence and the real-time operating state of the physical entity.

[0176] Based on the virtual predicted values ​​and physical measured values ​​contained in the paired state dataset, a multi-scale residual calculation framework is invoked to construct differential analysis chains at macroscopic, local, and microscopic scales, forming a coarse-to-fine deviation quantification structure. At the macroscopic scale, the global displacement distribution of the virtual predicted value grid nodes and the displacement distribution of the physical measured value sensor sampling locations are extracted, and a preliminary deviation matrix is ​​generated through global statistical differencing. At the local scale, region-weighted differencing is performed on significantly anomalous regions in the macroscopic deviation matrix, using weighting coefficients to improve the deviation calculation accuracy of edge nodes or thermally affected region nodes. At the microscopic scale, parameter-by-parameter residual calculations are performed on the predicted and measured material constitutive parameters of locally anomalous nodes, forming a material parameter deviation submatrix as a supplement.

[0177] Macroscopic, local, and microscopic deviation matrices are synthesized into a multi-scale deviation vector according to a preset ratio. A spatial interpolation algorithm is then used to map the virtual predicted values ​​to the sampling grid positions of the physically measured values. The interpolation process employs a three-dimensional spline interpolation method, adjusting the interpolation weights based on the Euclidean distance between nodes to ensure consistency between the mapped data and the physical sampling points. Within the interpolated sampling grid, the deviation values ​​of each node are recalculated, constructing a state deviation vector field covering the entire sampling space. The data structure of this vector field includes the three-dimensional coordinates, deviation value, and scale source label for each grid node. This processing method transforms the paired state dataset from the previous step into a complete state deviation vector field, achieving a spatially resolved representation of the differences between the virtual blow molding process state sequence and the physical operating state.

[0178] For example, during a certain blow molding cycle, the virtual predicted value mesh consists of 500 three-dimensional nodes, with node displacement ranging from 0.2 mm to 1.5 mm, and the predicted specific heat capacity of the material ranging from 1.2 to 1.8 kJ / (kg·K). The physical measured values ​​are sampled from 120 sensors distributed on the inner wall of the mold, with measured node displacement ranging from 0.25 mm to 1.6 mm, and measured specific heat capacity ranging from 1.25 to 1.85 kJ / (kg·K). The formula for calculating the macroscopic deviation matrix is: in The total number of nodes. For virtual predicted value vectors, This represents a vector of physically measured values. After macroscopic-scale deviation calculations, the mean nodal deviation is 0.05 mm. In local-scale deviation calculations, the weighting coefficient for specific heat capacity deviation within the heat-affected zone is set to 2, increasing the mean deviation to 0.07 kJ / (kg·K). In microscopic-scale deviation calculations, the deviations of specific heat capacity and thermal conductivity are calculated. Spatial interpolation uses three-dimensional spline interpolation, with the interpolation error controlled within 0.002. After mapping, a state deviation vector field containing 120 sampling points is generated. Execution results show that this state deviation vector field accurately characterizes the microscopic differences between the virtual and physical processes on the time and state axes, significantly improving the stability of subsequent modal decomposition and trend deviation index calculations.

[0179] S8.3: Perform mode decomposition on the state deviation vector field, extract the error feature spectrum containing high-frequency noise components and low-frequency drift components, and perform smooth estimation of the low-frequency drift components based on the Kalman filter algorithm to generate a trend deviation index for characterizing systematic model mismatch.

[0180] S8.4: Perform backpropagation optimization operation based on the trend deviation index, and use the adjoint variable method to calculate the sensitivity matrix of the internal state variables of the digital twin module to the trend deviation index, so as to generate a state update gradient vector containing the correction values ​​of grid node displacement and the correction values ​​of material constitutive parameters.

[0181] Based on the trend deviation index input, the sensitivity analysis process for the internal state variables of the digital twin module is initiated, with the calculation object set as the combination of grid node displacements and material constitutive parameters. The adjoint variable method is used to solve the partial derivative matrix of the internal state variables with respect to the trend deviation index. A numerical strategy combining finite difference and inverse integration is employed to ensure the analytical accuracy of the sensitivity matrix in both the spatiotemporal dimensions. Column vector normalization is performed on the sensitivity matrix to eliminate the influence of differences in the dimensions of different physical quantities on the iterative convergence. A state update gradient objective function is constructed, mapping the sensitivity matrix and trend deviation index to the gradient space through element-wise multiplication. The gradient vector is defined as the set of corrections for each state variable, and the amplitude of each component is adjusted according to its importance weight. For the material constitutive parameter components, a thermo-mechanical-phase transition coupling constraint check is performed to restrict the gradient components within the physically feasible region to avoid numerical divergence caused by parameter mutations. The above chained calculation generates a state update gradient vector containing corrections for grid node displacements and material constitutive parameters.

[0182] By using the adjoint variable method in chain derivation and gradient constraint processing, the trend deviation index of the previous step is transformed into serializable state update gradient data, thereby achieving adaptive alignment between the virtual blow molding process and the physical process at the level of internal state variables.

[0183] For example, in the application of a blow molding production line for plastic bottles used in daily chemical products, the trend deviation index is taken from the average displacement residuals of three key nodes in area C of the mold, in millimeters, with values ​​of 0.12, 0.15, and 0.09, respectively. The sensitivity matrix elements are obtained by solving using the adjoint variable method. The displacement component sensitivities are approximately 2.5, 3.1, and 2.0; the material specific heat capacity sensitivities are 0.008, 0.010, and 0.006; and the thermal conductivity sensitivities are 0.004, 0.005, and 0.003. After normalization, the sensitivity matrix and deviation vector are multiplied element-wise to generate a gradient vector: the displacement corrections are 0.30, 0.465, and 0.18 mm, the material specific heat capacity corrections are 0.00096, 0.0015, and 0.00054 J / (g·℃), and the thermal conductivity corrections are 0.00048, 0.00075, and 0.00027 W / (m·℃). When applying physical feasible region constraints, the thermophysical parameters are verified to be within the allowable range for the HDPE material grade, ensuring that the correction values ​​do not cause numerical instability in the model. This gradient vector significantly reduces the difference between the model-predicted displacement and the measured displacement in subsequent iterations, improves the consistency between the material's thermophysical properties and the physical process, and achieves rapid convergence of the virtual blow molding process in terms of material constitutive and displacement states.

[0184] S8.5: Based on the state update gradient vector, perform incremental iterative updates on the internal state variables of the digital twin module, and reload the corrected mesh node displacements and material constitutive parameters into the simulation engine to complete the dual adaptive alignment of the virtual blow molding process and the physical process on the time axis and state axis.

[0185] For those skilled in the art, various other corresponding changes and modifications can be made based on the technical solutions and concepts described above, and all such changes and modifications should fall within the protection scope of the claims of this invention.

[0186] Unless otherwise defined, the technical or scientific terms used herein shall have the ordinary meaning as understood by one of ordinary skill in the art to which this application pertains. The terms “first,” “second,” “third,” and similar terms used in this patent application specification and claims do not indicate any order, quantity, or importance, but are merely used to distinguish different components. Similarly, the terms “an” or “a” and similar terms do not indicate a quantity limitation, but rather indicate the presence of at least one. The terms “comprising” or “including” and similar terms mean that the element or object preceding “comprising” or “including” encompasses the element or object listed following “comprising” or “including” and its equivalents, and do not exclude other elements or objects. The “multiple” mentioned in the embodiments of this application refers to two or more. A and / or B indicate three possibilities: A; B; and A and B.

[0187] The above description is merely an exemplary embodiment of this application, but the scope of protection of this application is not limited thereto. Any person skilled in the art can easily conceive of various equivalent modifications or substitutions within the technical scope disclosed in this application, and such modifications or substitutions should all be covered within the scope of protection of this application. Therefore, the scope of protection of this application should be determined by the scope of the claims.

Claims

1. A method for intelligent control of blow molding process parameters based on digital twins, characterized in that, include: S1: Acquire the heat flow signals collected by the array sensors and record the spatial coordinate labels corresponding to each sensor to form a temperature field dataset; S2: Perform time-domain differentiation on the temperature field dataset to extract the heat flow inflection point features caused by the jump in specific heat capacity and the abrupt change in thermal conductivity, and generate a candidate event set. S3: Based on the confidence level of each event in the candidate event set, select valid anchor events, define the occurrence time of the valid anchor events as timestamp anchor points and the occurrence location as spatial coordinate anchor points, and construct spatiotemporal anchor points; S4: Match the corresponding target template from the preset template library based on the material properties and working parameters in the spatiotemporal anchor point; S5: Using the spatiotemporal anchor point as a constraint, reverse locking is performed on the simulation step size, mesh node displacement, and material constitutive parameter combination in the target template to generate parameter remapping results; S6: Based on the time offset standard deviation of phase transition events at the same location in the subsequent three consecutive blow molding cycles, calculate the weighting coefficient of the spatiotemporal anchor point and construct the attenuation function to generate a weighted anchor point confidence index. S7: If the weighted anchor confidence index is lower than the preset threshold, the collaborative verification process of adjacent redundant measurement points is initiated; otherwise, a simplified rheological model is used in conjunction with the parameter remapping results to perform interpolation and deduction, generating a virtual process state sequence to regulate the blow molding process parameters.

2. The intelligent control method for blow molding process parameters based on digital twins according to claim 1, characterized in that, Step S7 is followed by step S8: S8: Compare the deviation between the virtual process state sequence and the real-time running state of the physical entity, update the internal state variables of the digital twin module based on the comparison result, and complete the dual adaptive alignment of the virtual blow molding process and the physical process on the time axis and state axis.

3. The intelligent control method for blow molding process parameters based on digital twins according to claim 1, characterized in that, Step S3 specifically includes: Based on the transient heat flux inflection point amplitude of each event in the candidate event set, the ratio of the current event signal strength to the root mean square value of the background noise is calculated to generate a signal-to-noise ratio index. Using the signal-to-noise ratio index as an input variable, a time-series consistency verification algorithm based on a sliding window is executed to compare the temperature change rate trend within a preset time window before and after the current event and generate a time-series consistency coefficient. Based on the time series matching coefficient and the preset phase transition dynamics feature template, a multi-dimensional feature matching operation is performed to calculate the Euclidean distance between the current event feature vector and the standard phase transition feature vector, and a phase transition matching score is generated. Based on the phase transition matching score, a threshold filtering operation is performed on the candidate event set to remove invalid events with scores lower than a preset confidence threshold, retain high-scoring events and mark them as valid anchoring events, and generate a cleaned subset of valid anchoring events. Extract the occurrence time data and sensor spatial coordinate labels of each event in the effective anchored event subset, encapsulate the occurrence time as a timestamp anchor parameter and the spatial coordinates as a spatial coordinate anchor parameter, and combine them to generate the spatiotemporal anchor point used for spatiotemporal synchronization.

4. The intelligent control method for blow molding process parameters based on digital twins according to claim 1, characterized in that, Step S5 specifically includes: Based on the spatial coordinate anchor points contained in the spatiotemporal anchor points, the corresponding reference mesh node is located in the finite element mesh model of the target template, and the set of neighboring mesh nodes directly covered by the thermal conduction and diffusion effect of the reference mesh node is calculated according to the preset thermal influence radius algorithm, thereby generating a local subdomain mesh topology structure containing the reference mesh node and neighboring mesh nodes. For each grid node within the local subdomain grid topology, the predicted temperature field data at the current simulation step size is extracted, and the difference between the predicted temperature field data and the measured phase transition temperature data corresponding to the timestamp anchor point recorded in the spatiotemporal anchor point is calculated to generate a local temperature residual vector. Based on the local temperature residual vector, a sensitivity matrix for the material constitutive parameters is constructed using the adjoint variable method. Then, the material thermal property correction coefficients that can eliminate the local temperature residual vector are calculated in reverse iteration using the gradient descent optimization algorithm, generating an updated combination of material constitutive parameters that includes the corrected specific heat capacity parameter and the corrected thermal conductivity parameter. Based on the updated material constitutive parameter combination, the explicit integration time step of each grid node in the local subdomain grid topology is adaptively recalculated to determine the maximum allowable simulation step size that satisfies the numerical stability condition, and a locally optimized simulation step size sequence adapted to the current phase transition dynamics range is generated. The local optimization simulation step sequence, the updated material constitutive parameter combination, and the grid node displacements corrected based on measured positions are written together into the corresponding storage address of the target template to complete the local subdomain parameter remapping of the region directly affected by the spatiotemporal anchor point and generate the parameter remapping result.

5. The intelligent control method for blow molding process parameters based on digital twins according to claim 1, characterized in that, The temperature field dataset includes: deploying array sensors in key temperature measurement areas according to a preset spatial topology to collect continuous heat flow signals caused by the advance of the melt front in real time and generate raw voltage sequence data; performing spatial coordinate calibration on the sensors to build a static spatial index library, and mapping the raw voltage sequence to temperature values ​​through thermoelectric conversion to generate a time series temperature data stream; and then spatiotemporally fusing and binding the time series temperature data stream with three-dimensional spatial coordinate labels to obtain the temperature field dataset.

6. The intelligent control method for blow molding process parameters based on digital twins according to claim 1, characterized in that, The candidate event set includes: obtaining a clean temperature time series signal and calculating an instantaneous temperature gradient sequence by performing sliding window denoising and Kalman filtering on the temperature field dataset; extracting the transient heat flux curvature feature vector by performing second-order time-domain differentiation on the instantaneous temperature gradient sequence; and generating the candidate event set by identifying abnormal fluctuation points through zero-crossing detection and extreme value search.

7. The intelligent control method for blow molding process parameters based on digital twins according to claim 1, characterized in that, The parameter remapping result includes: locating the reference grid node based on the spatial coordinate anchor point and determining the local subdomain grid topology; calculating the local temperature residual between the virtual model and the physical entity; using the adjoint variable method and gradient descent optimization to obtain the material thermal property correction coefficient; adaptively determining the local optimization simulation step size sequence; and then writing the simulation step size, the updated material constitutive parameters, and the corrected grid node displacement into the target template to generate the parameter remapping result.

8. The intelligent control method for blow molding process parameters based on digital twins according to claim 1, characterized in that, The weighted anchor confidence index includes: extracting the occurrence time of the phase transition event at the same spatial coordinate anchor position within three consecutive blow molding cycles, calculating the standard deviation of the time offset between each cycle; constructing a decay function as a weighting coefficient based on the standard deviation, and multiplying the weighting coefficient by the original confidence level to obtain the weighted anchor confidence index.

9. The intelligent control method for blow molding process parameters based on digital twins according to claim 1, characterized in that, The decay function is an exponential decay function, whose decay rate constant is positively correlated with the time offset standard deviation, and the weighting coefficient decreases monotonically with the increase of the blow molding cycle number.

10. The intelligent control method for blow molding process parameters based on digital twins according to claim 1, characterized in that, The simplified rheological model uses the spatiotemporal boundary determined by two adjacent spatiotemporal anchor points as constraints, and employs the finite volume method to solve the temperature field control equations containing parameters of thermal conductivity, specific heat capacity, and density. It also performs continuous nonlinear interpolation derivation of the melt state between the two anchor points.

Images