Modeling method for thermal stress evolution of fire-resistant glass

Abstract

The present application relates to the technical field of fireproof glass thermodynamic simulation and structure safety analysis, in particular to a model simulation method for thermal stress evolution of fireproof glass, comprising: establishing a three-dimensional model of a glass base layer, a laminated layer and a coating layer of the fireproof glass and carrying out mesh division, and setting material models respectively; establishing a corresponding relationship between thermal physical parameters of each layer material and temperature and dynamically updating; setting a contact unit between layer interfaces and using a cohesive force model to represent interlayer mechanical properties; applying a temperature rising boundary condition on a fire surface, and carrying out thermal-structure bidirectional coupling calculation; monitoring the first principal stress and switching the solving mode to calculate the crack propagation path when the crack threshold is reached; combining displacement time history with fire resistance test data to correct model parameters and output simulation results; the present application realizes unified solution of the whole process of heat transfer, deformation, debonding, cracking and failure.

Description

Technical Field

[0001] This invention relates to the field of thermodynamic simulation and structural safety analysis technology for fireproof glass, specifically a model simulation method for the thermal stress evolution of fireproof glass. Background Technology

[0002] Fire-resistant glass is widely used in energy storage cabin observation windows, building fire-resistant partitions, and industrial high-temperature protection. Under fire conditions, it not only undergoes temperature transfer but also experiences various responses such as thermal expansion mismatch, interlaminar debonding, thermal bending deformation, and glass cracking. To ensure the integrity and fire resistance of fire-resistant glass under unilateral flame impact or continuous heating conditions, it is usually necessary to analyze its thermal failure process during the product design stage and use numerical simulation methods to pre-evaluate its fire resistance. However, existing methods for analyzing the heat resistance of fire-resistant glass often focus on calculating a single temperature field or employ simplified thermo-structural sequential analysis, treating the multi-layered structure, including the glass substrate, interlayer, and coating layer, as a homogeneous material. While these methods can provide preliminary temperature distribution or deformation results, under high-temperature conditions, the material's thermophysical parameters change significantly with temperature, and interlayer interfaces may experience damage, debonding, and air gap evolution. The glass substrate may also develop cracks after stress concentration. If fixed material parameters, constant interface constraints, and a single solution mode are still used, it is not only difficult to accurately characterize the entire process of fire-resistant glass from heating and deformation to cracking failure, but it is also prone to large deviations in the predicted results of failure time, crack location, and displacement response, thus affecting the accuracy of the fire resistance performance assessment of fire-resistant glass. Summary of the Invention

[0003] To address the aforementioned technical problems, this invention provides a model simulation method for the thermal stress evolution of fire-resistant glass. Specifically, the technical solution of this invention includes: S1. In the simulation space, establish three-dimensional models of the glass base layer, the laminated layer and the coating layer of the fireproof glass respectively, mesh the three-dimensional models and set the corresponding material models respectively. S2. Establish the correspondence between the thermal properties of each layer of material and temperature, and update the thermal properties based on the temperature calculation results of the mesh nodes during the simulation process; S3. Contact units are set at the interface between the glass substrate and the laminated layer, and at the interface between the laminated layer and the coating layer, and the interlayer mechanical properties of the contact units are set using the cohesive force model. S4. Apply a heating boundary condition to the fire-exposed surface on the outside of the glass substrate, and use implicit statics to perform a two-way thermo-structural coupling calculation to obtain the temperature distribution, thermal bending deformation and thermal stress distribution of each layer of material. S5. During the bidirectional coupling calculation, monitor the first principal stress of the glass base mesh. When the first principal stress reaches the preset crack initiation threshold, switch the solution method from implicit statics to explicit dynamics and use extended finite element method to calculate the crack propagation path. When the calculation shows that the crack penetrates the entire glass base and the surrounding contact units are completely debonded, the fireproof glass is determined to be in failure. S6. Obtain the displacement time history data obtained from the simulation, and obtain the test displacement data measured in advance through the fire resistance test of fireproof glass. Compare the displacement time history data with the test displacement data, correct the model parameters according to the comparison results, and output the simulation results.

[0004] Optionally, in step S1, the corresponding material models are set, including: A linear elastic brittle fracture model was established for the glass substrate; A generalized Maxwell viscoelastic model is set for the interlayer; an anisotropic model is set for the coating layer.

[0005] Optionally, in step S2, the thermophysical parameters include elastic modulus, thermal conductivity, specific heat capacity, and coefficient of thermal expansion. For glass substrates, the relationship between the elastic modulus and temperature is a non-linear decay function: ,in, For temperature The elastic modulus below, Reference temperature The initial elastic modulus, This is the temperature softening coefficient. It is a non-linear exponent; Step S2 also includes: using the glass softening point temperature as the phase transition temperature, and when the temperature of a mesh node is greater than or equal to the softening point temperature, switching the material model of the corresponding mesh to a viscous fluid model.

[0006] Optionally, in step S3, the interlayer mechanical properties of the contact elements are set using a cohesive force model, including: The bonding and debonding process at the interface was simulated using a bilinear traction-separation criterion. When the equivalent shear stress at the interface reaches the critical shear strength, the interface is considered to have begun to be damaged. When the relative slip displacement of the interface reaches the failure displacement, the stiffness of the contact element is reduced to zero, and the adjacent layers are determined to be completely debonded.

[0007] Optionally, in step S4, the thermal-structural bidirectional coupling calculation includes: Within the current increment step, the temperature field solver calculates heat conduction and heat radiation based on the thermal conductivity to obtain the temperature distribution of each layer of material. The structural field solver calculates thermal bending deformation and thermal stress distribution based on temperature distribution, combined with the coefficient of thermal expansion and interface state; The calculated thermal bending deformation and the air gap formed based on the interface debonding are obtained. The grid node coordinates and interlayer thermal resistance are updated, and the updated grid node coordinates and interlayer thermal resistance are input into the temperature field solver for the next incremental step.

[0008] Optionally, in step S5, when the first principal stress reaches the preset crack initiation threshold, the solution method is switched from implicit statics to explicit dynamics, and the extended finite element method is used to calculate the crack propagation path, including: The tensile strength that decreases with temperature is used as the crack initiation threshold; When the first principal stress of a certain grid cell is greater than or equal to the tensile strength that decays with temperature, the solution method for that grid cell and its adjacent grids within a preset topological order that reflects the extent of grid space expansion is switched from implicit statics to explicit dynamics. When the first principal stress is less than the tensile strength that decreases with temperature, the implicit static solution method is maintained. The discontinuity of the crack surface is simulated by the step function of the extended finite element method, the stress singularity at the crack tip is simulated by the asymptotic crack tip function, and the crack propagation path is calculated by explicit time integration.

[0009] Optionally, in step S6, the displacement time history data is compared with the experimental displacement data, and the model parameters are corrected based on the comparison results, including: Obtain the glass center displacement time history curve from the displacement time history data obtained from the simulation, and compare it with the measured curve of the fire resistance test of fireproof glass, which is used as the test displacement data. Calculate the peak deflection error between the displacement time history data and the test displacement data; when the peak deflection error exceeds the preset deflection error threshold, use the gradient descent method to correct the initial thermal expansion coefficient or the viscoelastic relaxation time of the interlayer until the error is within the preset range. When the error of the peak deflection is less than or equal to the preset deflection error threshold, the current model parameters are confirmed to be uncorrectable and the simulation results are output.

[0010] Compared with the prior art, the present invention has the following beneficial effects: This invention establishes three-dimensional models of the glass base layer, interlayer, and coating layer of fire-resistant glass and performs mesh generation. It combines differentiated material models for each layer, dynamically updates thermophysical parameters with temperature, uses interface cohesive contact elements, employs thermo-structural bidirectional coupling calculations, and switches from implicit statics to explicit dynamics after crack initiation, using extended finite element method to track crack propagation paths. Simultaneously, it compares simulated displacement time histories with fire resistance test displacement data to correct model parameters. This allows for the unified solution of heat transfer, thermal expansion mismatch, interlayer debonding, thermal bending deformation, crack propagation, and eventual failure of fire-resistant glass during fire exposure within a single timeframe. This invention is particularly suitable for evaluating the fire resistance performance of multi-layer fire-resistant glass in scenarios such as observation windows in energy storage cabins, fire-resistant partitions in buildings, and high-temperature protection in industrial applications. It has the advantages of improving the accuracy of predicting failure time, crack location, and displacement response, and enhancing the reliability of simulation verification of fire-resistant glass fire resistance capabilities. Attached Figure Description

[0011] To more clearly illustrate the technical solutions in the embodiments of this application and the prior art, the accompanying drawings used in the description of the embodiments and the prior art will be briefly introduced below: Figure 1 This is a flowchart illustrating the model simulation method for the thermal stress evolution of fireproof glass provided in the embodiments of this application. Detailed Implementation

[0012] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to specific embodiments.

[0013] Model simulation methods for the thermal stress evolution of fireproof glass include: S1. In the simulation space, establish three-dimensional models of the glass base layer, the laminated layer and the coating layer of the fireproof glass respectively, mesh the three-dimensional models and set the corresponding material models respectively. S2. Establish the correspondence between the thermal properties of each layer of material and temperature, and update the thermal properties based on the temperature calculation results of the mesh nodes during the simulation process; S3. Contact units are set at the interface between the glass substrate and the laminated layer, and at the interface between the laminated layer and the coating layer, and the interlayer mechanical properties of the contact units are set using the cohesive force model. S4. Apply a heating boundary condition to the fire-exposed surface on the outside of the glass substrate, and use implicit statics to perform a two-way thermo-structural coupling calculation to obtain the temperature distribution, thermal bending deformation and thermal stress distribution of each layer of material. S5. During the bidirectional coupling calculation, monitor the first principal stress of the glass base mesh. When the first principal stress reaches the preset crack initiation threshold, switch the solution method from implicit statics to explicit dynamics and use extended finite element method to calculate the crack propagation path. When the calculation shows that the crack penetrates the entire glass base and the surrounding contact units are completely debonded, the fireproof glass is determined to be in failure. S6. Obtain the displacement time history data obtained from the simulation, and obtain the test displacement data measured in advance through the fire resistance test of fireproof glass. Compare the displacement time history data with the test displacement data, correct the model parameters according to the comparison results, and output the simulation results.

[0014] This embodiment provides a model simulation mechanism for the thermal stress evolution of fireproof glass, such as... Figure 1 As shown; specifically, taking the fire resistance verification of composite fireproof glass for observation windows of lithium battery energy storage compartments as the unified main scenario, the entire process of fire-induced temperature rise, interlayer debonding, glass cracking and eventual failure is digitally reproduced. The observation window is located on the side wall of the energy storage compartment and is used by personnel to observe the situation inside the compartment in the early stages of a fire. Therefore, it is required not only to be fire-resistant, but also to maintain its integrity for a period of time when the flame is impacted from one side. To this end, a three-dimensional model of the glass base layer, the laminated layer and the coating layer is established in the simulation space, and the three-layer structure is discretized separately. Then, incremental coupling is performed between the thermal field and the structural field to solve the problem. Specifically, in the modeling stage, a fireproof glass model with a length of 1200mm and a width of 800mm can be built first based on the actual sample size; for example, assuming that the thickness of the glass base layer is 6mm, the thickness of the laminate layer is 2mm, and the thickness of the coating layer is 0.1mm, it can be divided into base layer units G1 and G2, laminate layer unit J1, and coating layer unit M1 along the thickness direction. The model is further divided into three regions along its length: left, middle, and right, thus forming local meshes such as G1-L, G1-C, and G1-R. The purpose of this is not to limit the number of elements, but to ensure that each element can be assigned an independent state variable during the solution process. It should be noted that, in order to conform to the conventional terminology of finite element analysis, the terms such as element, basic element, and adjacent element appearing in this embodiment strictly refer to the smallest independent computational micro-element formed after the three-dimensional model is meshed, i.e., the mesh or mesh element. Similarly, the monitoring and calculation of the physical state of the grid is essentially a numerical analysis of the discrete units and their corresponding grid nodes. The two have the same meaning in the technical solution. The glass base layer is used to bear the main mechanical constraints, the laminated layer is used for buffering and interlayer bonding under thermal shock, and the coating layer is used to change the thermal radiation absorption and reflection characteristics. Therefore, the three layers are described by different materials. Furthermore, during the thermal property update stage, a temperature-driven parameter table can be established for each mesh node. For example, for node N1, its initial temperature is 20℃, and the corresponding elastic modulus, thermal conductivity, specific heat capacity, and coefficient of thermal expansion are taken as initial values. When the temperature of N1 is calculated to rise to 180℃ in the next increment step, the corresponding update value is read from the preset function or data table and the update value is filled back into the adjacent cells containing N1. If element G1-C consists of four nodes with node temperatures of 160℃, 170℃, 180℃, and 190℃, then the average node temperature of 175℃ can be used as the basis for updating the thermal properties of this element in this step. Alternatively, the integration point temperature update method can be used, which updates directly at the integration point level when the solver supports it. This process ensures that the material properties are not fixed constants, but evolve dynamically with the temperature rise of the fire. During the interlayer interface treatment stage, contact units are inserted between the glass substrate and the interlayer, and between the interlayer and the coating layer. For ease of understanding, the interface between the substrate and the interlayer can be discretized into three contact segments: C1, C2, and C3. Initially, the normal stiffness and tangential stiffness of the three segments are both effective values, representing continuous interlayer bonding. When the difference in thermal expansion leads to an increase in the relative slip of the middle segment C2, the cohesive model will reduce the local load-bearing capacity according to its traction-separation response; if C2 further reaches complete failure, the interface force transmission capacity corresponding to this segment will drop to the preset lower limit of load-bearing capacity, thereby forming a local thermal resistance change in the subsequent thermal field, and manifesting as additional bending induced by debonding in the structural field. During the coupled calculation phase, a temperature rise boundary condition is applied to the fire-exposed surface outside the glass substrate. A standard fire resistance curve or a measured temperature rise curve of the thermal runaway flame inside the energy storage chamber can be used. During the solution process, the temperature of each node in the current step is first calculated by the temperature field, and then the temperature distribution is read by the structural field. The thermal bending and thermal stress are calculated in combination with thermal expansion and interface state. The structural field outputs the geometric position after deformation and the amount of interface opening, and the temperature field then corrects the radiation distance and interlayer thermal resistance accordingly. Taking a simplified incremental step as an example: In step 10, the temperature of the center node on the exposed side rises to 420℃, while the temperature of the center node on the unexposed side is 180℃. Based on this, the structural field calculates the center deflection to be 8mm. In step 11, due to local cracking at the central interface, a 0.3mm air gap is formed. The temperature field uses this air gap as a new thermal resistance layer, causing a change in the heating rate of the unexposed side. The structural field then continues to update the deflection and stress. By repeating this cycle, the time-history-based temperature, deformation, and stress results can be obtained. During the crack initiation and failure determination stage, the first principal stress of each mesh in the glass substrate is continuously monitored. For ease of deduction, it is assumed that the first principal stress of the three key elements G1-L, G1-C, and G1-R at a certain moment is 32MPa, 51MPa, and 40MPa, respectively, and the preset crack initiation threshold is 48MPa. Then only G1-C meets the crack initiation condition. At this time, the solution method is switched from implicit statics based on overall steady-state tracking to explicit dynamics, which is more stable for rapid crack propagation, and the characterization of the propagating finite element crack is initiated near G1-C. If subsequent calculations show that the crack extends from the center to the left and right edges and eventually penetrates the thickness of the glass substrate, and the surrounding contact units C1, C2, and C3 have completely detached, then the sample can be determined to be faulty. If the crack only stays on the local surface and does not penetrate, or if the substrate is cracked but the surrounding interface still maintains some constraint, then the subsequent incremental step calculations will continue instead of terminating prematurely. In addition, as a supplementary explanation, exception handling logic also needs to be set in the actual solution process; if the temperature growth rate of a certain incremental step exceeds the preset temperature rise threshold, causing the implicit solution to fail to converge, the time step can be automatically reduced, for example, from 1 second to 0.2 seconds; if it still fails to converge after reduction, the state of the previous stable step is saved and the calculation is switched to local explicit transition. If some interface elements exhibit negative volume or abnormal penetration due to mesh distortion, local remapping can be performed to project the state variables of the abnormal elements onto newly generated neighboring meshes, thus avoiding global computation interruption. If the experimental displacement curve is missing data for a certain time period, the existing data segments are first compared by aligning them with the time, and the missing intervals are marked as not participating in the correction, so as to ensure that subsequent parameter updates are not interfered with by invalid data. In the fire-resistant design project of the energy storage compartment, engineers simulated an observation window rated for 60 minutes of fire resistance. The sample was installed in a steel frame with the fire-exposed surface facing the inside of the compartment. After the simulation started, the first 15 minutes showed rapid heating of the fire-exposed surface and slow heating of the unexposed surface, with the deflection in the middle increasing from 0 to 6 mm. After 20 minutes, the laminated layer and the coating layer partially debonded, forming a zone of thermal resistance change. At 33 minutes, the first principal stress of the middle glass base unit reached the crack initiation condition, and the system automatically entered the crack propagation calculation. After 35 minutes, the crack extended to the edge but did not penetrate, and the sample maintained its integrity. After 41 minutes, the crack penetrated and the surrounding interface failed. The system determined that the configuration did not meet the target fire resistance time. The simulation center displacement time history curve was compared with the furnace test curve second by second. If the simulated peak deflection was 28 mm while the actual measurement was 31 mm, the parameter correction process was triggered and recalculated until the curve error fell within the allowable range. The purpose of this step is to incorporate the heat transfer, deformation, debonding and fracture of multi-layer fireproof glass under fire conditions into the same time frame for unified solution, so as to achieve traceable reproduction of failure mechanism and verifiable output of fire resistance of engineering sample.

[0015] Furthermore, in step S1, corresponding material models are set, including: a linear elastic brittle fracture model for the glass substrate; a generalized Maxwell viscoelastic model for the laminated layer; and an anisotropic model for the coating layer.

[0016] This embodiment provides a modeling mechanism for the differentiated material responses of multilayer fireproof glass. Specifically, in the aforementioned energy storage cabin observation window scenario, if all three layers are approximated with a uniform isotropic elastic material, although a rough temperature field can be obtained, the hysteretic deformation of the interlayer is often underestimated at high temperatures, and the thermal conductivity difference of the coating layer along the in-plane and thickness directions cannot be correctly reflected, which leads to deviations in the cracking time. Therefore, this embodiment sets up material models that are more in line with actual working conditions for the base layer, interlayer, and coating layer respectively. Specifically, the glass substrate can adopt a linear elastic brittle fracture model. The reason is that the deformation of conventional fireproof glass substrate in the elastic stage before cracking is lower than the preset deformation threshold. However, after the tensile stress exceeds the strength threshold, the crack propagation rate increases nonlinearly and the plastic dissipation energy is lower than the preset energy threshold. For ease of explanation, it is assumed that the elastic modulus of a substrate unit is 70 GPa and the Poisson's ratio is 0.22 at 20℃. When the temperature rise has not yet entered the softening zone, the unit only produces a small elastic strain under a tensile stress of 40MPa; when its first principal stress approaches the tensile limit, it does not release energy through continuous yielding, but directly enters the brittle fracture process controlled by the crack criterion; this makes the behavior of high stiffness before cracking and sudden drop in load after cracking in the structural field closer to actual observation. For the interlayer, a generalized Maxwell viscoelastic model can be used. Under fire conditions in the energy storage compartment, the interlayer does not fail immediately, but exhibits the characteristic of being constrained for a short time and gradually relaxing over a long period. To conduct finite element numerical simulation, it can be simplified into two parallel branches: one fast relaxation branch and one slow relaxation branch. It is assumed that at 100°C, the fast relaxation branch decays from the initial modulus to 30% within 10 seconds, and the slow relaxation branch decays from the initial modulus to 60% within 300 seconds. If the central region experiences sustained shear displacement at 200 seconds, the interlayer can initially suppress the relative slippage between the two panes of glass. However, as time accumulates, the stress will be gradually released, thereby inducing stress redistribution between the layers. This type of model can more realistically represent the viscoelastic decay under the continuous action of a fire. For the coating layer, an anisotropic model can be used; in practice, the functional coatings on the surface of fireproof glass often have different thermal or mechanical responses along the coating surface and perpendicular to the coating surface; for ease of understanding, it is assumed that the in-plane thermal conductivity of the coating layer is 1.8 W / (m·K) and the thermal conductivity in the thickness direction is 0.9 W / (m·K). Under the same fire conditions, heat is more likely to diffuse along the coating surface rather than quickly penetrate the thickness direction. If an isotropic single-value thermal conductivity is still used, the heat flow diffusion path will be oversimplified, causing the prediction of the temperature difference between the central hot spot and the edge to be distorted. The anisotropic model can characterize this difference through the material parameter matrix in different directions, without requiring complex formula derivation in the explanation stage. In one alternative implementation, if complete relaxation spectrum data of the interlayer cannot be obtained in the actual project, a low-order generalized Maxwell approximation model with two or three branches can be used first, and parameter regression can be performed through subsequent displacement time history correction. If the coating layer is extremely thin and separate meshing would result in an excessively large element length-to-thickness ratio, it can be treated as an equivalent shell layer or interface functional layer, as long as the treatment can still reflect the directional thermal characteristics. If the base layer has initial conditions such as tempered prestress, a prestress field can be superimposed on the initial conditions without changing the basic modeling idea of ​​linear elasticity before cracking and brittle fracture after cracking. In the same energy storage cabin observation window sample, the engineers first used a three-layer unified elastic model to perform trial calculations and found that the deflection at the simulation center was 12 mm at 30 min and 18 mm at the furnace test, with the error value exceeding the preset allowable range. Further investigation revealed that the main reason was that the equivalent stiffness of the interlayer was set too high, failing to reflect high-temperature relaxation; at the same time, the thermal diffusion of the coating layer was treated too evenly, resulting in a small thermal gradient in the middle; after changing to a model of brittle base layer, viscoelastic interlayer, and anisotropic coating layer, the simulated deflection at 30 min increased to 17.2 mm, and the temperature difference distribution between the exposed and unexposed surfaces was closer to the experimental curve. The purpose of this step is to enable the intrinsic responses of different functional layers to be expressed in a targeted manner, thereby achieving a more accurate thermo-mechanical coupling simulation basis that closely approximates the actual failure process.

[0017] Furthermore, in step S2, the thermophysical parameters include elastic modulus, thermal conductivity, specific heat capacity, and coefficient of thermal expansion; For glass substrates, the relationship between the elastic modulus and temperature is a non-linear decay function. This function uses a reference temperature Initial elastic modulus Temperature softening coefficient and nonlinear exponent Calculated temperature elastic modulus below ; Step S2 also includes: using the glass softening point temperature as the phase transition temperature, and when the temperature of a mesh node is greater than or equal to the softening point temperature, switching the material model of the corresponding mesh to a viscous fluid model.

[0018] This embodiment provides a temperature-dependent material update mechanism for the high-temperature softening stage. Specifically, in the aforementioned basic scheme, if only the constant parameters at room temperature are used for solving, the early thermal deflection can be described. However, as the fire in the energy storage compartment continues to develop and the temperature approaches the glass softening zone, the stiffness of the base layer will be significantly overestimated, resulting in cracks appearing too late or even incorrect prediction of the long-term integrity of the glass. Therefore, this embodiment further introduces dynamic updates of thermophysical parameters as temperature changes and switches the material state after reaching the softening point. Among them, temperature softening coefficient and nonlinear exponent The modulus scatter data was obtained by conducting high-temperature mechanical property tests on glass substrate samples at different temperature levels in advance, and then the nonlinear decay function was calibrated by curve fitting using the least squares method. Specifically, the updated objects include at least the elastic modulus, thermal conductivity, specific heat capacity, and coefficient of thermal expansion. Taking glass substrate as an example, the decrease of elastic modulus with temperature can be described by nonlinear decay. For ease of understanding, assuming the reference temperature is 20℃, the initial elastic modulus is 70GPa, the temperature softening coefficient is a value within the first preset positive range, and the nonlinear exponent is a value greater than 1, then at 100℃, the modulus decrease rate is less than the first preset rate threshold, while near 400℃ and 500℃, the decrease rate is greater than the second preset rate threshold. This reflects that glass does not soften linearly from the beginning, but its stiffness decreases rapidly after entering a specific high-temperature region; after writing this function into the material update module, the corresponding modulus can be automatically obtained by reading the current temperature of the unit for each incremental step. It should be noted that, to prevent excessively high temperatures from causing the decay term to exceed 1 and resulting in a negative modulus that violates physical laws, the update logic must be configured as follows: Set a lower limit for the residual modulus; when the modulus result calculated by the above nonlinear decay function is lower than the lower limit, the system automatically truncates and maintains it at the lower limit to ensure the physical rationality of the material stiffness data in the extreme high temperature range close to the softening point. To conduct a microscopic simulation, assume that the average temperatures of the three basic units B1, B2, and B3 at a certain moment are 80℃, 260℃, and 520℃, respectively. For B1, its elastic modulus may still remain above 90% of its initial value; for B2, the modulus has dropped to about 70% of its initial value; and for B3, the modulus may be only 40% or even lower. At the same time, thermal conductivity and specific heat capacity can also be updated by looking up tables according to temperature. For example, the specific heat capacity of B2 is higher than that of B1, indicating that its heat absorption capacity is enhanced; the increase in the coefficient of thermal expansion means that under the same temperature difference, B3 is more prone to thermal expansion mismatch than B1; this distributed update method with different units and different parameters at the same time is an important source of subsequent thermal stress differentiation. Furthermore, the glass softening point temperature is used as the state switching threshold. It should be noted that the instruction manual and phase transition temperature here do not refer to microscopic phase transitions in the crystallographic sense, but are artificially set algorithm switching judgment thresholds based on the macroscopic stress performance in engineering, in order to achieve a smooth dynamic transition from elastic solid to significantly viscous flow. When the temperature of a certain grid node reaches or exceeds the softening point, its grid is no longer treated as a traditional solid elastic matrix, but is switched to a viscous fluid model. For example, if the softening point is set to 620℃, and the temperatures of the four nodes of unit B3 are 610℃, 625℃, 630℃, and 635℃ respectively, then B3 can be identified as entering the softening state according to the node ratio or integral point criterion. After entering the softening state, the way the unit is subjected to shear and bending changes. It is no longer mainly elastic recovery, but instead exhibits viscous flow and macroscopic structural instability deformation. This can explain the phenomenon that some furnace tests have not yet formed typical brittle cracks, but the glass has obviously sagged. In terms of algorithm integration, state flags can be used for switching. For example, each unit can be assigned a state code of 0, 1, or 2 to represent the normal solid state, the transitional softened state, and the viscous flow dynamics, respectively. If the average temperature of a unit is lower than the softening point minus the hysteresis band, it remains in the solid state. If it enters the hysteresis region near the softening point, it can be marked as a transitional state, and its elastic modulus can be reduced, and its damping or viscosity parameters increased. If two consecutive incremental steps are both above the softening point, it switches to the flow dynamics to prevent the solid-liquid state from jumping back and forth due to numerical fluctuations. This approach reflects the physical process and is also beneficial to numerical stability. In one alternative implementation, if the node temperature distribution of a certain unit is extremely uneven, for example, half of the nodes in the same unit are below the softening point and half are above the softening point, then the integral point can be used as the primary method for determination, or the unit can be locally refined before determination, to avoid the average temperature masking local hot spots; if the glass material used in the project does not provide a clear softening point, an initial estimate can be given based on the standard value of similar materials, and then reverse-checked by the sag and displacement mutation points in the fire resistance test; if the time step becomes extremely sensitive after switching, a lower viscosity limit can be introduced in the softening range to avoid numerical divergence; At the 38th minute of the simulation of the energy storage compartment observation window, the temperature of the central fire zone was significantly higher than that of the edge; at this time, the temperature of the central unit B3 reached 628℃, while that of the edge unit B1 was 410℃; accordingly, the system switched B3 to a viscous flow dynamic, while B1 remained solid, and B2 was in a transitional softened state; due to the step-like decrease in stiffness in the central part and the aggravation of thermal expansion mismatch, the central deflection rapidly increased from 22mm to 29mm within 90 seconds; if the entire domain solid treatment is still used, the sudden increase in displacement in this section cannot be correctly captured; The purpose of this step is to uniformly express the nonlinear softening and post-softening flow dynamics of the base material at high temperatures, thereby achieving continuous simulation of the stages before and after high-temperature instability, sagging, and cracking.

[0019] Further, in step S3, the interlayer mechanical properties of the contact elements are set using a cohesive force model, including: The bonding and debonding process at the interface was simulated using a bilinear traction-separation criterion. When the equivalent shear stress at the interface reaches the critical shear strength, the interface is considered to have begun to be damaged. When the relative slip displacement of the interface reaches the failure displacement, the stiffness of the contact element is reduced to zero, and the adjacent layers are determined to be completely debonded.

[0020] This embodiment provides an interface damage mechanism for describing interlayer adhesion degradation. Specifically, in the aforementioned scheme, even if the material models of the base layer, the interlayer, and the coating layer are distinguished, if only a never-separating bonding relationship is adopted between the layers, the impact of interlayer peeling, local bulging debonding, and local air gap evolution on heat conduction and stress path will be underestimated in the middle and late stages of an energy storage tank fire. Therefore, this embodiment adopts a cohesive force model at the interface and uses a bilinear traction-separation criterion to express the entire process from complete adhesion to complete debonding. Specifically, the bilinear traction-separation relationship can be simplified into two stages: the first stage is the linear rising stage, which means that when the interface separates or slips at a small displacement, the traction force increases with the increase of the relative displacement, indicating that the interface is still transmitting the load; the second stage is the damage and degradation stage, which means that when the traction force reaches its peak, the interface stiffness begins to decrease until the relative displacement reaches the failure displacement and the interface loses its bearing capacity. It should be noted that, in order to maintain consistency between the terminology and physical representation throughout the text, the traction force in this embodiment when explaining the bilinear traction-separation criterion is specifically mapped to equivalent shear stress in the current thermal expansion stress mode dominated by interlayer slippage; the peak value of the corresponding traction force is strictly equivalent to the equivalent shear stress reaching the critical shear strength; the aforementioned relative displacement also refers to the relative slip displacement within the interface surface. For ease of explanation, it is assumed that the base layer-insulated layer interface is discretized into three segments: I1, I2, and I3, all with a critical shear strength of 1.2 MPa. However, due to the higher temperature in the middle, the relative slip of I2 increases the fastest. When the equivalent shear stress of I2 first reaches 1.2 MPa, the system marks it as the beginning of damage. If the relative slip displacement of I2 continues to increase from 0.05 mm to the preset failure displacement of 0.18 mm, its stiffness drops to zero, corresponding to complete debonding. In numerical implementation, a damage variable D can be set for each contact unit, with a value ranging from 0 to 1. Initially, D=0, indicating no damage. After damage begins, D gradually increases. When D=1, it indicates complete failure. The specific parameter evolution process is as follows: If the damage variables of I1, I2, and I3 at a certain moment are 0.1, 0.7, and 0, respectively, it indicates that the middle segment has significantly degraded, while the two sides still remain bonded. At this point, the force transmission capacity of region I2 in the structural field is weakened, which will cause local deflection concentration; in the temperature field, a small opening may appear at the corresponding position of I2, changing the interlayer heat transfer path; as the calculation progresses, if I2 reaches D=1, a through thermal resistance gap may be formed in the middle, while I1 and I3 can still provide residual constraints to the two sides. Furthermore, the equivalent shear stress can be synthesized from the tangential components of the interface, while the relative slip displacement is one of the main controlling variables of failure evolution. The advantage of this setting is that, under fire conditions, interlaminar failure is not always dominated by normal opening, but in most cases, it is triggered by in-plane shear slip caused by thermal expansion mismatch. Especially when the observation window is constrained by the steel frame boundary, the thermal expansion difference between the middle and the edge is more likely to form a tangential slip peak. Therefore, controlling interface degradation in conjunction with shear strength and slip failure displacement is more in line with engineering practice. In one alternative implementation, if a local contact segment experiences short-term compression, i.e., the normal displacement is negative while the tangential slippage is large, shear damage can still be allowed to occur, and it is not forced to remain intact simply because it has not opened up. If some interface data is difficult to test directly, initial values ​​can be obtained through small-sample peeling tests or high-temperature shear tests, and then back-calibrated using full-window displacement curves. If the calculation shows that the interface is repeatedly damaged and recovered in adjacent incremental steps, an irreversible damage evolution rule can be adopted, which means that once D increases, it will no longer regress, in order to better reflect the actual irreversible characteristics of thermally induced debonding. At the 26th minute after the observation window of the energy storage compartment was exposed to fire, the central segment I2 near the coating layer first showed large tangential slip due to concentrated temperature difference, and the equivalent shear stress reached the critical value, and the system recorded that it entered the damage state; by the 29th minute, the slip displacement of I2 reached the failure threshold, and the central coating layer and the interlayer separated locally; at the 31st minute, the central segment of the base layer-interlayer interface also began to be damaged, forming a double interface degradation; this change caused the thermal resistance of the central area to increase, the local temperature difference to steepen, and promoted the initiation of cracks in the central part of the base layer in the following minutes; The purpose of this step is to explicitly incorporate the evolution process of interlayer bonding from complete adhesion to gradual debonding into the simulation, thereby achieving an accurate characterization of the changes in heat transfer paths and stress redistribution caused by interface degradation.

[0021] Furthermore, in step S4, the thermal-structural bidirectional coupling calculation includes: Within the current increment step, the temperature field solver calculates heat conduction and heat radiation based on the thermal conductivity to obtain the temperature distribution of each layer of material. The structural field solver calculates thermal bending deformation and thermal stress distribution based on temperature distribution, combined with the coefficient of thermal expansion and interface state; The calculated thermal bending deformation and the air gap formed based on the interface debonding are obtained. The grid node coordinates and interlayer thermal resistance are updated, and the updated grid node coordinates and interlayer thermal resistance are input into the temperature field solver for the next incremental step.

[0022] This embodiment provides a bidirectional coupling mechanism for geometric deformation feedback. Specifically, in the aforementioned interface damage scheme, if only unidirectional coupling is performed, i.e., temperature is calculated first and then stress is calculated without feeding the deformation results back to the thermal field, although the stress distribution at a certain moment can be obtained, it cannot reflect the closed-loop process of glass bulging, interface opening, thermal resistance change, and temperature redifferentiation. For single-sided fire-affected components such as the observation window of the energy storage compartment, this closed loop is particularly critical. Therefore, this embodiment forms bidirectional coupling by gradually interacting with the thermal field and the structural field. Specifically, in the current increment step, the temperature field solver calculates the temperature distribution based on the thermal conductivity, radiation boundary, and current geometry; the structural field reads this temperature distribution and, in conjunction with the thermal expansion coefficient and interface damage state, calculates the thermal bending deformation and thermal stress; then, the newly obtained nodal displacements and interface openings are returned to the temperature field to update the mesh nodal coordinates and interlayer thermal resistance. For ease of deduction, we can assume that there are two nodes, P1 and P2, at the center position; initially, the distance between them is 2.0 mm, corresponding to the original thickness of the interlayer; if the structural field calculation results in a local opening to 2.4 mm due to debonding, the temperature field in the next step will no longer follow the original heat conduction path, but will be re-solved according to the series thermal resistance of the interlayer + local air gap. Furthermore, the thermal resistance update can be handled using simplified rules; assuming that when a certain interface segment does not debond, its equivalent thermal resistance is R1; after a local air gap of 0.2mm occurs, the equivalent thermal resistance increases to R2; if the air gap expands to 0.8mm, it can be further increased to R3. In this process, considering that the heat transfer due to thermal radiation will increase exponentially within the thin air gap formed by the opening of the interface under the high temperature of a fire, if the series thermal resistance is calculated solely based on the extremely low thermal conductivity of the air medium, it will lead to abnormal blockage of local heat flow and cause numerical distortion. Therefore, when the temperature field solver receives the air gap thickness and updates the thermal resistance, it needs to simultaneously convert the radiative heat flux density between the two surfaces of the microscopic air gap into an equivalent thermal conductivity and compensate it into the comprehensive thermal conductivity model of the air gap. In specific conversion, the equivalent thermal conductivity is equal to the sum of the thermal conductivity of the air medium and the equivalent radiative thermal conductivity. The equivalent radiative thermal conductivity is obtained by calculating the amount of radiative heat exchange on both sides of the micro gap using the Stefan-Boltzmann law and dividing it by the ratio of the temperature difference between the two sides of the surface to the current thickness of the gap. Therefore, under the same heat flow, the heating rate of adjacent areas on the unexposed side will be slower, while the hot spots inside the exposed side may be more concentrated. At the same time, the update of geometric coordinates will also affect the radiative heat transfer path. For example, after the middle bulges out, its viewing angle coefficient with the surrounding frame changes, and the radiation receiving distribution changes accordingly. Although it is not required to expand the complete radiation formula in the embodiment, it should be clear that the coordinate update is not only for displaying the deformation diagram, but also directly participates in the next round of thermal calculation. Let's illustrate this with microscopic numerical values: Assume the first... At the end of the step, the central node displacement is 6 mm, and the I2 segment forms a 0.3 mm air gap; the temperature field is thus determined in the first step. In the first step, the interlayer thermal resistance of the I2 region was increased by 20%; the result was that... At the end of the step, the center temperature of the unexposed surface was revised from the originally predicted 210℃ to 198℃, while the hot spot near the center of the exposed surface rose to 462℃; based on this, the structural field was calculated that the stress peak shifted slightly from the original location to the adjacent unbonded area; it can be seen that debonding not only changes the local temperature, but also changes the position of the stress peak, thus affecting the subsequent judgment of the crack initiation area. In one optional implementation, if the displacement generated by a certain incremental step exceeds the preset displacement tolerance, causing the mesh distortion to exceed the threshold, geometric smoothing or local reconstruction can be performed first, and then the updated mesh can be transferred back to the thermal field; if the air gap formed by debonding is smaller than the preset micro gap threshold, for example, smaller than the preset lower limit of 0.02mm, it can be treated as an air gap layer without being treated separately, but can be equivalently treated by increasing the interface thermal resistance coefficient to avoid a surge in computation. If a segment is completely debonded and then comes into contact again under compression, the structural field can allow the normal contact to recover, but the thermal resistance characteristics after the damage will still be maintained in the thermal field, or the thermal conduction will be recalculated as closed contact, depending on the contact heat conduction model used. During the 24th to 32nd minute of the fire exposure at the observation window of the energy storage compartment, slight bulging first appeared in the central region, and local debonding of the interface of the laminated layer formed an air gap. At 1450 seconds, the system calculated that the central deflection was 10 mm and the local air gap was 0.15 mm. By 1560 seconds, due to the addition of new interlaminar thermal resistance in the thermal field, there was a temperature difference of 12°C between the center temperature of the unexposed surface and the unfeeded model, and the temperature gradient near the center increased, resulting in a faster increase in tensile stress in the edge transition zone of the structural field. Finally, the crack initiation location shifted from the absolute center to a region 10 mm off-center, which was more consistent with the crack initiation location in the experiment. The purpose of this step is to establish a closed-loop feedback between the thermal field, interface state, and geometric deformation, so as to achieve synchronous tracking of temperature redistribution and stress redistribution during the fire process.

[0023] Furthermore, in step S5, when the first principal stress reaches the preset crack initiation threshold, the solution method is switched from implicit statics to explicit dynamics, and the extended finite element method is used to calculate the crack propagation path, including: The tensile strength that decreases with temperature is used as the crack initiation threshold; When the first principal stress of a certain grid cell is greater than or equal to the tensile strength that decays with temperature, the solution method for that grid cell and its adjacent grids within a preset topological order that reflects the extent of grid space expansion is switched from implicit statics to explicit dynamics. When the first principal stress is less than the tensile strength that decreases with temperature, the implicit static solution method is maintained. The discontinuity of the crack surface is simulated by the step function of the extended finite element method, the stress singularity at the crack tip is simulated by the asymptotic crack tip function, and the crack propagation path is calculated by explicit time integration.

[0024] This embodiment provides a hybrid solution mechanism for the rapid crack propagation stage. Specifically, the aforementioned two-way coupling scheme can accurately characterize the thermal stress accumulation before crack initiation. However, if implicit static continuous solution is still used after crack initiation, problems such as difficulty in convergence, unclear crack path jumps, or frequent mesh re-drilling often occur. Especially in brittle substrates such as the observation window of the energy storage compartment, the rate of crack propagation after crack initiation is greater than the rate of thermal bending deformation in the early stage. Therefore, in this embodiment, the solution strategy is switched when the crack initiation condition is met, and the extended finite element method is used to track the crack path. Specifically, the crack initiation threshold is not taken as a fixed room temperature strength, but as a tensile strength that decays with temperature. The tensile strength that decays with temperature is obtained by interpolation from a pre-established tensile strength-temperature relationship table. This relationship table is constructed by conducting uniaxial tensile tests on the glass substrate material at a series of set characteristic high temperatures, measuring and recording the measured values ​​of the ultimate tensile strength under different temperature conditions. This is because even if the glass has not yet softened during the later stages of a fire, its tensile strength has already decreased significantly. For example, assuming that the temperatures of units K1, K2, and K3 are 300℃, 420℃, and 500℃ at a certain moment, their tensile strengths will decrease to 80%, 60%, and 45% of their room temperature values, respectively. If the first principal stresses of the three are 35MPa, 38MPa, and 34MPa respectively, then although the stress of K3 is not the highest, its strength decreases the most, and it may be the first to meet the crack initiation condition that the stress is greater than or equal to the strength. This can avoid the actual situation of missing the preferential cracking of the high-temperature weakened zone by judging solely by the stress peak value. Regarding the switching range, it is not necessary to change the entire model to explicit dynamics at the same time. Instead, the crack initiating element and its adjacent meshes within the preset topological order can be switched. For ease of understanding, if the preset topological order is 1, then adjacent elements sharing edges or faces with K2 are included in the explicit region. If the topological order is 2, then it is expanded outward by another ring. This approach maintains the computational stability of the crack zone while avoiding excessively small time steps and a surge in computational load caused by explicit solutions across the entire window. For example, if there are 25 base elements in the middle section, and only 9 of them are within the second-order neighborhood of K2, then the explicit region is first limited to these 9 elements, while the remaining regions are solved in the original manner or coupled to the boundary of the explicit region. Furthermore, the extended finite element method describes cracks using additional functions instead of relying on pre-made meshes along crack paths; among them, the step function is used to express the discontinuity of displacement on both sides of the crack surface, and the crack tip function is used to express the concentrated characteristics of the stress field near the crack tip. For the implementation of the disclosure, it can be understood as follows: when the crack passes through the interior of K2, there is no need to cut K2 into two new elements. Instead, by adding degrees of freedom, the two sides of the same element exhibit separate displacement fields. When the crack tip advances to the vicinity of the boundary of the adjacent element, the stress gradient near the tip is maintained by the crack tip enhancement function. The explicit time integral is responsible for gradually advancing the crack propagation within a very small time step, which is suitable for handling transient fracture. To conduct a microscopic deduction, assume that the current number is... When K2 satisfies the crack initiation condition, the system switches K2 and its first-order neighbors K1, K3, K4, and K5 to explicit regions; at the time when... At any moment, among them, For time steps, the crack extends from the center of K2 to the upper right and penetrates into K3; the... At that moment, the stress concentration at the interface between K3 and its adjacent element K6 above it intensifies, and the crack continues to extend towards the edge; If the crack length subsequently increases to penetrate the effective width of the base layer, and the interface elements near the crack zone have failed, then the failure determination is established; if the crack stops due to the obstruction of the local compressive stress zone, the explicit zone can be monitored for a period of time; when the crack no longer advances in several consecutive steps, the stable region far from the crack tip can be restored to the conventional solution as appropriate to reduce the consumption of computational resources. In addition, if multiple elements simultaneously meet the crack initiation conditions at a certain moment, the main crack source can be determined according to the principle of prioritizing the one with the smallest strength margin after temperature weakening, or multiple crack sources can be allowed to propagate in parallel; if energy reflection or displacement discontinuity occurs at the boundary between the explicit region and the surrounding implicit region, a transition zone element can be set, or compatible constraints can be used at the boundary. If the local time step is excessively reduced due to the small minimum unit size, local merging can be performed after the crack moves away from the area, or only the crack area of ​​the base layer can be refined while keeping the mesh of the adhesive layer and the coating layer coarser. If the experiment shows that the crack often starts from the corner defect, the corner micro-notch can also be introduced into the initial model without affecting the switching logic of this embodiment. 33 minutes after the observation window of the energy storage compartment was exposed to fire, the temperature of unit K2 in the upper center had risen to 465℃, and its tensile strength had decreased significantly. Although the first principal stress value of a certain unit at the edge was slightly higher, the stress-to-strength ratio of K2 reached the crack initiation threshold first. Therefore, the system preferentially established an explicit crack zone near K2. In the explicit integration time of 2 seconds, the crack extended to the upper right region along an oblique upward path and continued to advance along the banded region with a large temperature gradient. This path is consistent with the oblique cracking phenomenon from the upper center to the edge observed in the furnace test, rather than a simple horizontal or vertical crack. The purpose of this step is to adopt a solution strategy more suitable for brittle transient fracture after crack initiation, so as to achieve stable tracking of crack origin, propagation direction and penetration time.

[0025] Further, in step S6, the displacement time history data is compared with the experimental displacement data, and the model parameters are corrected based on the comparison results, including: Obtain the glass center displacement time history curve from the displacement time history data obtained from the simulation, and compare it with the measured curve of the fire resistance test of fireproof glass, which is used as the test displacement data. Calculate the peak deflection error between the displacement time history data and the test displacement data; when the peak deflection error exceeds the preset deflection error threshold, use the gradient descent method to correct the initial thermal expansion coefficient or the viscoelastic relaxation time of the interlayer until the error is within the preset range. When the error of the peak deflection is less than or equal to the preset deflection error threshold, the current model parameters are confirmed to be uncorrectable and the simulation results are output.

[0026] This embodiment provides a parameter correction mechanism for improving simulation credibility. Specifically, the aforementioned implementation methods have established a complete solution chain from heating and debonding to crack propagation. However, if the initial material values ​​come from a public manual or supplier's nominal values, rather than the measured values ​​of the corresponding batch of samples, the simulation results may still deviate from the fire resistance test. Especially in the scenario of the observation window of the energy storage compartment, differences in the adhesive layer formulation, coating process and installation pre-tightening will be reflected in the center displacement curve; therefore, this embodiment uses displacement time history comparison to adjust the model parameters. Specifically, the time history curve of the glass center displacement can be selected as the comparison object. The reason is that the center deflection can usually reflect the three aspects of temperature gradient, interlayer constraint and material softening, and it is easy to continuously obtain the data through the displacement gauge during the test. To facilitate the explanation of the specific numerical comparison process, it is assumed that the center displacement of the simulation curve at 10 minutes, 20 minutes, 30 minutes, and 35 minutes is 2mm, 9mm, 18mm, and 27mm, respectively, while the actual experimental measurements are 3mm, 11mm, 21mm, and 31mm, respectively. It can be seen that the simulation is generally too stiff. Further calculation of the peak deflection error shows that if the simulated peak value is 27mm and the experimental peak value is 31mm, the calculated relative error is 12.9%, which is higher than the preset deflection error threshold of 10%. Therefore, parameter correction is required. Parameter correction can be achieved using the gradient descent method, prioritizing the adjustment of the initial thermal expansion coefficient or the viscoelastic relaxation time of the interlayer, which have a greater impact on displacement. For example, if the current value of the thermal expansion coefficient is slightly too small, the thermal bending driving force will be underestimated. For example, if the current value of the relaxation time is too long, the continuous constraint of the interlayer at high temperature will be overestimated. To illustrate the update process, the objective function can be defined as the square of the difference between the simulated peak deflection and the experimental peak deflection, plus a weighted sum of displacement differences at several key moments; specifically, the expression for the error objective function is: in, and The peak deflections are from the simulation and the experiment, respectively. and The first Simulation and experimental center displacement data at key moments, These are the weighting coefficients for the corresponding time points. The total number of key moments selected. These are the model parameters to be corrected; During the initial iteration, the system slightly increases the thermal expansion coefficient by one step and performs a fast solution again; if the error decreases, it continues to update in that direction; if the error increases, it decreases the step size or adjusts the relaxation time. Furthermore, to clarify the specific calculation process of the algorithm steps, the above-mentioned optimization process based on gradient descent is explicitly quantified here using the numerical difference method: since the finite element two-way coupled model is an implicit nonlinear system, the error objective function cannot be directly evaluated. Find the analytic derivative, therefore the parameter to be corrected is needed. For example, the coefficient of thermal expansion is used in the first... In the iteration, the system uses the ... The parameters to be corrected in the next iteration Using small step sizes as a baseline, first... Simulation calculation to obtain the first Finite difference gradient of the next iteration: This gradient is the value of the error objective function after increasing the parameters by a small step size. Compared with the current error objective function value Calculations show that, based on this finite difference gradient... According to the updated formula And combined with learning rate Calculate the first The parameters to be corrected in the next iteration Among them, the learning rate The dimension of the parameter to be corrected is Dimensions and gradient The ratio of dimensions is used to ensure the consistency of dimensions on both sides of the updated formula. This underlying numerical logic, which uses step size perturbation to find an approximate gradient and guide the descent optimization, is the very essence of the aforementioned incremental step size and update based on the error direction. For example, in the first round, the coefficient of thermal expansion is increased by 3%, and the peak error is reduced from 12.9% to 8.4%, which is within the allowable range, so the update stops. If it still does not meet the standard, a second round of correction is made to the relaxation time of the interlayer. Furthermore, gradient descent does not require adjusting only one parameter at a time; an alternating update strategy can also be used. For example, the first round corrects the thermal expansion coefficient, the second round corrects the fast relaxation branch time constant, and the third round fine-tunes the slow relaxation branch time constant. If an update causes the error to increase, the parameters are rolled back to the previous round and the learning step size is reduced. This avoids drastic parameter oscillations and allows convergence to a more reasonable set of parameters specific to the sample within a limited number of iterations. In addition, if there are noise spikes in the test displacement curve, smoothing can be performed first, and then the peak value and key moment points can be extracted; if the peak value error meets the requirements, but the deviation of the crack initiation time is still large, the crack initiation time error can be added as a second correction index. If the thermal expansion coefficient and relaxation time are adjusted but convergence is still not achieved, it indicates that the deviation may come from boundary constraints, coating anisotropy, or interface parameter settings. In this case, the multi-parameter joint correction mode can be switched. If the experiment only retains the displacement data for the first 30 minutes due to sensor failure, the material's early response can be fitted with the curve for that period first, and then indirectly verified by combining the later damage photos or fracture time. In the same energy storage cabin observation window project, the first version of the model used data from the material handbook, and the peak displacement of the simulation center was 4mm smaller than that of the furnace test. Based on this, the system determined that the equivalent stiffness of the sample in the model was greater than the actual physical stiffness. Therefore, the initial thermal expansion coefficient of the base layer was increased by 2%, and after recalculation, the displacement curve from 20 minutes to 35 minutes shifted upward as a whole; however, the peak value was still 2mm different. By shortening the relaxation time of the fast relaxation branch of the interlayer by 15%, the interlayer releases the constraint earlier. After solving again, the peak error is reduced to less than 5%, and the displacement growth inflection point near 33 minutes is closer to the actual measurement. The final corrected model is used for batch simulation screening of observation windows of different sizes in the same series. The purpose of this step is to use the experimental observation results in reverse to correct the model, thereby achieving closed-loop calibration between simulation parameters and actual sample conditions, and improving the reliability of subsequent engineering prediction results.

[0027] It should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are not intended to limit it. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can be made to the technical solutions of the present invention without departing from the spirit and scope of the technical solutions of the present invention.

Claims

1. A model simulation method for the thermal stress evolution of fireproof glass, characterized in that, include: S1. In the simulation space, establish three-dimensional models of the glass base layer, the laminated layer and the coating layer of the fireproof glass respectively, divide the three-dimensional models into meshes and set the corresponding material models respectively; S2. Establish the correspondence between the thermal properties of each layer of material and temperature, and update the thermal properties based on the temperature calculation results of the mesh nodes during the simulation process; S3. Contact units are provided at the interface between the glass substrate and the laminated layer, and at the interface between the laminated layer and the coating layer, and the interlayer mechanical properties of the contact units are set using a cohesive force model. S4. Apply a heating boundary condition to the fire-exposed surface outside the glass substrate, and use implicit statics to perform a two-way thermo-structural coupling calculation to obtain the temperature distribution, thermal bending deformation and thermal stress distribution of each layer of material. S5. During the bidirectional coupling calculation process, monitor the first principal stress of the glass base layer mesh. When the first principal stress reaches the preset crack initiation threshold, switch the solution method from implicit statics to explicit dynamics, and use extended finite element method to calculate the crack propagation path. When the calculation shows that the crack penetrates the entire glass base layer and the surrounding contact units are completely debonded, the fireproof glass is determined to be in failure. S6. Obtain the displacement time history data obtained from the simulation, and obtain the test displacement data measured in advance through the fire resistance test of fireproof glass. Compare the displacement time history data with the test displacement data, correct the model parameters according to the comparison results, and output the simulation results.

2. The model simulation method for the thermal stress evolution of fireproof glass according to claim 1, characterized in that, In step S1, the corresponding material models are set up, including: A linear elastic brittle fracture model is set for the glass substrate; A generalized Maxwell viscoelastic model is set for the adhesive layer; an anisotropic model is set for the coating layer.

3. The model simulation method for the thermal stress evolution of fireproof glass according to claim 1, characterized in that, In step S2, the thermophysical parameters include elastic modulus, thermal conductivity, specific heat capacity, and coefficient of thermal expansion. For the glass substrate, the elastic modulus has a non-linear decay function with respect to temperature: ,in, For temperature The elastic modulus below, Reference temperature The initial elastic modulus, This is the temperature softening coefficient. It is a non-linear exponent; Step S2 further includes: using the glass softening point temperature as the phase transition temperature, and when the temperature of a mesh node is greater than or equal to the softening point temperature, switching the material model of the corresponding mesh to a viscous fluid model.

4. The model simulation method for the thermal stress evolution of fireproof glass according to claim 1, characterized in that, In step S3, the interlayer mechanical properties of the contact unit are set using a cohesive force model, including: The bonding and debonding process at the interface was simulated using a bilinear traction-separation criterion. When the equivalent shear stress at the interface reaches the critical shear strength, the interface is considered to have begun to be damaged. When the relative slip displacement of the interface reaches the failure displacement, the stiffness of the contact element is reduced to zero, and the adjacent layers are determined to be completely debonded.

5. The model simulation method for the thermal stress evolution of fireproof glass according to claim 1, characterized in that, In step S4, the thermal-structural bidirectional coupling calculation includes: Within the current increment step, the temperature field solver calculates heat conduction and heat radiation based on the thermal conductivity to obtain the temperature distribution of each layer of material. The structural field solver calculates the thermal bending deformation and thermal stress distribution based on the temperature distribution, combined with the coefficient of thermal expansion and the interface state; The calculated thermal bending deformation and the air gap formed according to the interface debonding determination are obtained, the grid node coordinates and interlayer thermal resistance are updated, and the updated grid node coordinates and interlayer thermal resistance are input into the temperature field solver of the next incremental step.

6. The model simulation method for the thermal stress evolution of fireproof glass according to claim 1, characterized in that, In step S5, when the first principal stress reaches the preset crack initiation threshold, the solution method is switched from implicit statics to explicit dynamics, and the extended finite element method is used to calculate the crack propagation path, including: The tensile strength that decreases with temperature is used as the crack initiation threshold. When the first principal stress of a certain grid cell is greater than or equal to the tensile strength that decays with temperature, the solution method for that grid cell and its adjacent grids within a preset topological order that reflects the extent of grid space expansion is switched from implicit statics to explicit dynamics. When the first principal stress is less than the tensile strength that decreases with temperature, the implicit static solution method is maintained. The discontinuity of the crack surface is simulated by the step function of the extended finite element method, the stress singularity at the crack tip is simulated by the asymptotic crack tip function, and the crack propagation path is calculated by explicit time integration.

7. The model simulation method for the thermal stress evolution of fireproof glass according to claim 1, characterized in that, In step S6, the displacement time history data is compared with the experimental displacement data, and the model parameters are corrected based on the comparison results, including: Obtain the glass center displacement time history curve from the displacement time history data obtained from the simulation, and compare it with the measured curve of the fire resistance test of fireproof glass, which is the test displacement data. Calculate the peak deflection error between the displacement time history data and the test displacement data; when the peak deflection error exceeds the preset deflection error threshold, use the gradient descent method to correct the initial thermal expansion coefficient or the viscoelastic relaxation time of the interlayer until the error is within the preset range; When the error of the peak deflection is less than or equal to the preset deflection error threshold, it is confirmed that the current model parameters do not need to be corrected and the simulation results are output.

Images