Simulation design method of graphite heater

Abstract

The application relates to the technical field of simulation design, and discloses a simulation design method of a graphite heater. The method comprises the following steps: establishing a space mapping database by extracting a notch corner node set and giving an anisotropic resistivity tensor, tracking a focal point self-enhancement index by node in electric-thermal bidirectional iterative coupling solution, and outputting a failure thermal cycle prediction number by pre-calibration mapping relationship. The application solves the problem that the existing method cannot quantitatively map the notch geometric parameters and the failure thermal cycle number of the heater in the design stage, and the problem that the life prediction and geometric optimization are mutually separated.

Description

Technical Field

[0001] This application relates to the field of simulation design technology, and in particular to a simulation design method for a graphite heater. Background Technology

[0002] Graphite heaters are widely used in high-temperature equipment such as vacuum sintering furnaces, single crystal growth furnaces, and ground thermal intensity tests for hypersonic vehicles. A typical structure consists of a cylindrical body with alternating through-slots at both ends along the axial direction, forming a sawtooth current path. High current and low voltage are applied to achieve resistance heating. To predict the temperature field distribution and structural reliability of the heater before manufacturing, existing technologies generally employ finite element simulation design methods. This involves establishing an electrothermal coupling model, assigning the graphite resistivity as a scalar function of global temperature to the global mesh nodes, solving the electric field equations and temperature field equations, outputting a global temperature field cloud map and heat flux density distribution, and combining this with response surface methodology to perform parametric scanning optimization of geometric parameters such as slot width and wall thickness.

[0003] However, the sawtooth slotted structure of the graphite heater causes the current to deflect at the slot end. This abrupt change in cross-sectional geometry at the deflection point results in a high density of current lines, creating a current concentration effect. The Joule heat density at this point is significantly higher than that in the straight section. Simultaneously, the resistivity of isostatically pressed graphite is anisotropic, with differences between the in-plane and thickness directions, and both exhibit nonlinear changes with temperature. Existing methods treat resistivity as a globally uniform temperature scalar function, failing to capture the localized temperature rise at the slot corner caused by current concentration and its impact on the resistivity distribution. This neglects the reverse effect of local resistivity deflection on current redistribution, resulting in the complete absence of the actual current concentration-local resistivity change-current redistribution positive feedback effect at the slot corner in the simulation model. Consequently, the peak Joule heat density at this location is systematically underestimated.

[0004] Because the aforementioned positive feedback effect cannot be quantitatively represented in existing simulations, engineers cannot determine whether a dangerous hot spot self-reinforcing behavior exists at a certain slot corner geometry during the design phase, nor can they establish a quantitative correlation between hot spot intensity and the actual working life of the heater—since hot spot intensity has not been extracted as a separately calculable index, a feedback loop between life prediction and geometric design cannot be formed. The direct consequence of this is that the slot fillet radius design of graphite heaters can only rely on engineering experience; the actual number of failure cycles of the heater is far lower than the simulation prediction value; and the failure location is inconsistent with the high-temperature region in the simulation temperature field, reflecting the systematic prediction inaccuracies of existing simulation methods under extreme operating conditions. Summary of the Invention

[0005] This application provides a simulation design method for graphite heaters, which solves the problem that existing graphite heater simulation design methods cannot quantify the positive feedback effect between anisotropic resistivity and current concentration at the slot corner, resulting in a systematic underestimation of hot spot intensity and inaccurate prediction of failure location. It also solves the problem that existing methods cannot establish a quantitative mapping between slot geometric parameters and the number of thermal cycles for heater failure during the design stage, resulting in a disconnect between lifetime prediction and geometric optimization.

[0006] This application provides a simulation design method for a graphite heater, the simulation design method for the graphite heater comprising: Step S1: Perform mesh refinement processing on the slot corner area of ​​the graphite heater, extract the nodes belonging to the slot corner area in the refined mesh, and obtain the slot corner node set; Step S2: Assign the anisotropic resistivity tensor to each node in the set of slot corner nodes, and establish a spatial mapping database of the resistivity of each node as a function of temperature. Step S3: Substitute the spatial mapping database into the electrothermal bidirectional iterative coupling solution. In each iteration, extract the Joule heat density of each node in the slot corner node set. Using the ratio of the Joule heat density increment of each node in two adjacent iterations to the global average Joule heat density increment, perform iterative convergence segment averaging on each node in the slot corner node set to obtain the hotspot self-enhancing index of each node. Step S4: Substitute the hot spot self-reinforcing index into the pre-calibrated mapping relationship between the hot spot self-reinforcing index and the number of failure thermal cycles to obtain the predicted number of failure thermal cycles for the graphite heater.

[0007] The technical solution provided in this application introduces a spatial positioning mechanism—a set of slot corner nodes—into the simulation design method for graphite heaters. This reduces the scope of refined simulation modeling from uniformly distributed across the entire domain to the critical current-deflection region at the slot corner. This allows subsequent assignment of anisotropic resistivity tensors and Joule heat density tracking calculations to be concentrated in the most physically active local regions, rather than uniformly distributing computational resources across flat segments that contribute little to hotspots. Furthermore, the spatial mapping database uses a key-value pair structure with node spatial coordinates as indexes and temperature-dependent anisotropic resistivity tensors as values. This ensures that the resistivity update of each node in each iteration strictly corresponds to the node's current temperature and the principal anisotropic direction. This achieves a gradual capture of the impact of local resistivity deflection on current redistribution in the bidirectional electrothermal iterative coupling solution—a level of modeling accuracy that is fundamentally impossible under existing assumptions of uniform resistivity across the entire domain.

[0008] The hotspot self-enhancing index, defined as the core computational quantity in this application, is characterized by the average ratio of the Joule heat density increment at the slot corner node to the global average Joule heat density increment during the iteration convergence phase. Its algorithmic contribution lies in transforming the strength of the positive feedback effect, which was originally hidden within the iteration process and difficult to directly observe from the temperature field contour map, into a dimensionless scalar index that can be quantified during the simulation phase. This allows the coupling strength between current concentration at the slot corner and anisotropic resistivity changes to be objectively measured during the design phase, rather than relying on engineering experience. Substituting this index into a pre-calibrated power function mapping relationship directly outputs the predicted number of failure thermal cycles, enabling the simulation design method for graphite heaters to quantitatively predict the heater's service life for the first time without physical prototyping. This integrates the geometric optimization of the slot fillet radius and the satisfaction of life constraints into a single closed simulation design process. Attached Figure Description

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

[0010] Figure 1 This is a schematic diagram of an embodiment of the simulation design method for a graphite heater in this application. Figure 2 This is a schematic diagram illustrating the change of the single-step self-enhancement ratio of each node in the slot corner node set under different fillet radii as a function of the number of iteration steps in the embodiments of this application; Figure 3 This is a schematic diagram of the global Joule heat density distribution of the sawtooth cross section of the graphite heater after the convergence of the electrothermal bidirectional iterative coupling solution in the embodiments of this application. Detailed Implementation

[0011] This application provides a simulation design method for a graphite heater. The terms "first," "second," "third," "fourth," etc. (if present) in the specification, claims, and accompanying drawings of this application are used to distinguish similar objects and are not necessarily used to describe a specific order or sequence. It should be understood that such data can be interchanged where appropriate so that the embodiments described herein can be implemented in a sequence other than that illustrated or described herein. Furthermore, the terms "comprising" or "having" and any variations thereof are intended to cover a non-exclusive inclusion; for example, a process, method, system, product, or apparatus that comprises a series of steps or units is not necessarily limited to those steps or units explicitly listed, but may include other steps or units not explicitly listed or inherent to such processes, methods, products, or apparatus.

[0012] For ease of understanding, the specific process of the embodiments of this application is described below. Please refer to [link / reference]. Figure 1 One embodiment of the simulation design method for graphite heaters in this application includes: Step S1: Perform mesh refinement processing on the slot corner area of ​​the graphite heater, extract the nodes belonging to the slot corner area in the refined mesh, and obtain the slot corner node set. Specifically, the set of nodes at the slot corner is a subset of nodes selected from the global finite element mesh of the graphite heater based on the criterion that the spatial coordinates of the nodes fall within the region near the arc surface of the slot corner. The cylindrical body of the graphite heater has through slots alternately opened from both ends along the axial direction. The current deflects at the end of the slot, and the abrupt change in cross-sectional geometry at this deflection results in denser current lines, with a significantly higher current density than in the straight section. At least eight layers of fine-grained mesh are set in the radial direction, with a minimum element side length of 0.05 mm, to ensure that the numerical solution of the potential gradient in this region achieves convergence accuracy and to avoid underestimation of the current concentration effect due to sparse meshing.

[0013] Step S2: Assign the anisotropic resistivity tensor to each node in the slot corner node set, and establish a spatial mapping database of the resistivity of each node as a function of temperature. Specifically, the spatial mapping database is a set of key-value pairs indexed by the spatial coordinates of nodes and their values ​​being the temperature-dependent anisotropic resistivity tensors at those nodes. For nodes within the slot corner node set, the resistivity tensor distinguishes between in-plane and thickness components, both described by piecewise polynomials with temperature as the independent variable. For straight-segment nodes outside the set, both components take the same value, degenerating into isotropic resistivity. The in-plane resistivity temperature fitting coefficient set is obtained by least-squares fitting after high-temperature testing of graphite samples from the same batch using the four-probe method. The thickness anisotropy ratio is taken as the average of the measured values ​​from the same batch, typically between 1.05 and 1.30. Binding these coefficients to the spatial coordinates of each node forms a complete database, which is updated in each iteration based on the current temperature field, ensuring that the conductivity tensor always corresponds to the real-time temperature.

[0014] Step S3: Substitute the spatial mapping database into the electrothermal bidirectional iterative coupling solution. In each iteration, extract the Joule heat density of each node in the slot corner node set. Using the ratio of the Joule heat density increment of each node in two adjacent iterations to the global average Joule heat density increment, perform iterative convergence segment averaging on each node in the slot corner node set to obtain the hotspot self-enhancing index of each node. Specifically, the hotspot self-enhancement index is the arithmetic mean of the ratio of the Joule heat density increment of each node in the slot corner node set to the global average Joule heat density increment in each iteration step during the latter half of the electrothermal iteration convergence. The Joule heat density is calculated by bilinearly shrinking the potential gradient vector of the current iteration step and the conductivity tensor of the node. The global average Joule heat density is the volume-weighted average of the Joule heat densities of all nodes. The iteration convergence segment is determined by adaptively backtracking from the last iteration step. The criterion is that the maximum standard deviation of the single-step self-enhancement ratio of all nodes in the slot corner node set is less than a preset stability threshold. Iteration steps in the early stage of iteration where the temperature of each node is not yet stable and the ratio is disturbed by the initial conditions are excluded. The interval where the increment of each node has tended to be stable is retained. The obtained mean can reflect the local amplification factor of the node near the steady state. When the self-reinforcing index of a hot spot at a node exceeds the preset safety threshold of 1.50, it indicates that the Joule heat density of the node continues to increase at a rate of more than 1.50 times the average growth rate of the entire region. This means that there is a positive feedback hot spot driven by the local deflection of anisotropic resistivity, and the node is included in the dangerous corner node set.

[0015] Step S4: Substitute the hot spot self-reinforcing index into the pre-calibrated mapping relationship between the hot spot self-reinforcing index and the number of failure thermal cycles to obtain the predicted number of failure thermal cycles for the graphite heater.

[0016] Specifically, the number of predicted failure thermal cycles is obtained by substituting the global maximum hot spot self-reinforcing index into a power function mapping relationship. This mapping relationship is calibrated using thermal cycling test data from four sets of standard specimens with different fillet radii. The hot spot self-reinforcing index for each set of specimens is used as the independent variable, and the number of thermal cycles at the first visible crack is used as the dependent variable. Two coefficients in the power function are then fitted. When the number of predicted failure thermal cycles is less than the preset target number of thermal cycles, the target number of thermal cycles is substituted back into the power function to obtain the inverse, yielding a target hot spot self-reinforcing index threshold that meets the life requirement. This threshold is then used to replace the preset safety threshold, and the process is repeated until the number of predicted failure thermal cycles meets the constraint. The current fillet radius is then output as the final design result.

[0017] In one specific embodiment, step S1 includes: A parametric three-dimensional geometric model is established for the cylindrical body of the graphite heater. The parametric three-dimensional geometric model includes the outer diameter, wall thickness, total axial height, through slots alternately opened from both ends of the cylinder along the axial direction, the spacing between adjacent slots, the width of a single slot, and the fillet radius at the corner of the slot opening. Based on the parametric three-dimensional geometric model, the entire domain is divided into hexahedral structured meshes. For the corner region of the slot, no less than eight layers of mesh are set in the radial direction with the center of the corner arc as the reference. Sparse meshes are set for the straight sections of the cylinder to obtain the global mesh model. Based on the global mesh model, all nodes are classified according to whether their location belongs to the corner area of ​​the slot, thus obtaining the set of slot corner nodes. Record the spatial coordinates, slot number, corner azimuth, and anisotropic principal direction vector of each node in the slot corner node set to obtain the slot corner node set carrying the node position attributes.

[0018] Specifically, in the parametric 3D geometric model, the outer diameter and wall thickness determine the radial cross-sectional area of ​​the cylinder, which in turn determines the cross-sectional width of the current path; the total axial height and the spacing between adjacent slots together determine the effective length of each serrated conductive path; the width of a single slot determines the interval between adjacent conductive segments; the fillet radius at the slot opening is a key variable parameter in this method, initially set at 0.5 mm. This value corresponds to the lower limit of the minimum machinable fillet radius of a milling cutter in actual machining, using this as a starting point to ensure that subsequent optimizations start from the extreme cases within the feasible range of the process. Through slots are alternately opened from both ends of the cylinder along the axial direction, forming a serrated current path, causing the current to bend at the end of the slot opening. The current lines are denser at the bend, and the current density is significantly higher than that of the straight section. This area is the slot opening corner region that this method focuses on modeling.

[0019] In the global mesh generation, at least 8 mesh layers are set radially around the corner region of the slot, with the center of the corner arc as the reference. The minimum side length of the element is 0.05 mm. This value is derived from the analysis of the ratio of typical wall thickness to corner radius of graphite heaters. When the corner radius is set to the minimum value of 0.5 mm, the single-layer thickness of the 8 radial mesh layers is approximately 0.06 mm, which allows the numerical solution of the potential gradient in this region to converge and avoids the underestimation of the current concentration effect. The mesh side length of the straight section of the cylinder is 0.5 mm, forming a contrast in density with the corner region, controlling the total number of nodes in the global region while ensuring the calculation accuracy at the corner. The anisotropic principal direction vector records the unit vectors of the in-plane direction and thickness direction of the graphite lattice at the node in the global coordinate system. The two are orthogonal and are used to project the in-plane and thickness resistivity components to the global coordinate system during subsequent resistivity tensor assembly. The corner azimuth angle records the circumferential position of the corner of the node relative to the axis of the cylinder, which is used to distinguish the geometric differences of nodes corresponding to different slots in the case of multiple slots.

[0020] In one specific embodiment, step S2 includes: Based on the anisotropic principal direction vectors recorded by each node in the slot corner node set carrying node position attributes, tensor assembly processing is performed on the in-plane resistivity component and the thickness resistivity component of each node to obtain the initial anisotropic resistivity tensor of each node. The measured high-temperature resistivity data of the same batch of graphite material samples under the four-probe method were processed by piecewise polynomial least squares fitting of the resistivity in the in-plane direction and the resistivity in the thickness direction, respectively, to obtain the temperature fitting coefficient group of resistivity in the in-plane direction and the anisotropy ratio in the thickness direction. Based on the in-plane resistivity-temperature fitting coefficient set and the thickness-direction anisotropy ratio, temperature-dependent assignment is performed on each component of the initial anisotropic resistivity tensor to obtain the temperature-dependent anisotropic resistivity tensor of each node. By indexing and binding the spatial coordinates of each node in the slot corner node set carrying node position attributes with the temperature-dependent anisotropic resistivity tensor, a spatial mapping database of the resistivity of each node as a function of temperature is obtained.

[0021] Specifically, the assembly of the anisotropic resistivity tensor is based on the anisotropic principal direction vector recorded at each node. The product of the in-plane direction unit vector and its own dyadic vector is multiplied by the in-plane direction resistivity component, and then the product of the thickness direction unit vector and its own dyadic vector is multiplied by the thickness direction resistivity component. The sum of these two terms yields a 3×3 symmetric resistivity matrix at that node. For straight segment nodes outside the slot corner node set, the thickness direction anisotropy ratio is taken as 1.00, and the resistivity tensor degenerates into an isotropic scalar multiplied by an identity matrix.

[0022] In the four-probe method for high-temperature resistivity testing, standard samples were cut from the same batch of isostatically pressed graphite. The test temperature range covered 300K to 2800K, with a test point taken every 100K. Two sets of samples were prepared independently for testing, one along the in-plane direction and the other along the thickness direction, resulting in two sets of discrete temperature and resistivity data sequences. A second-order polynomial was fitted to each of the in-plane resistivity discrete data sequences in the two temperature ranges of 300K to 1000K and 1000K to 2800K, respectively. The fit expressions converged using the least squares criterion. The fitted expressions are as follows: when hour:

[0023] when hour:

[0024] in This is absolute temperature, measured in Kelvin. The resistivity is in-plane, and the two temperature ranges together contain... There are a total of 6 fitting coefficients.

[0025] Taking a certain grade of isostatically pressed graphite material as an example, the discrete data of in-plane resistivity measured by the four-probe method are shown in the table below:

[0026] The coefficients obtained by fitting the above data using the least squares criterion are: Low temperature range ;High temperature section The sum of squared residuals are respectively goodness of fit All are greater than 0.999. Anisotropy ratio in the thickness direction. The arithmetic mean of the ratio of the resistivity in the thickness direction to the resistivity in the in-plane direction of the two sets of samples at each test temperature point was taken. This was obtained from the actual measurements of the above material. That is, resistivity in the thickness direction .

[0027] The spatial mapping database uses the spatial coordinate triplet of each node as the index key and the temperature-dependent anisotropic resistivity tensor and corresponding fitting coefficient set of that node as the value, storing it as a key-value pair structure. In the bidirectional iterative coupling solution of electrothermal energy, each iteration looks up a table using the node coordinates as the key according to the current temperature field. The current temperature is substituted into the piecewise polynomial of the corresponding node to obtain the real-time in-plane resistivity component. Then, it is multiplied by the thickness anisotropy ratio to obtain the thickness resistivity component. The resistivity tensor of the current iteration step is reconstructed by dyadic assembly. After inversion, the conductivity tensor is obtained. It is substituted into the steady-state current continuity equation, i.e., the divergence of the product of conductivity tensor and potential gradient is set to zero. Fixed potential boundary conditions at the electrode contact surface and electrically insulating boundary conditions on the outer surface are applied. After finite element discretization, a system of linear equations with node potential as unknowns is solved to obtain the global potential distribution. Then, the potential gradient is calculated from the node potential difference and node spacing. It is bilinearly condensed with the conductivity tensor to obtain the Joule heat density of each node, in watts per cubic meter. This is substituted into the transient heat conduction equation as the volume heat source term to solve the temperature field.

[0028] In one specific embodiment, step S3 involves substituting the spatial mapping database into the electrothermal bidirectional iterative coupled solution, and extracting the Joule heat density of each node in the slot corner node set in each iteration, including: Substitute the anisotropic resistivity tensor corresponding to the current temperature of each node in the spatial mapping database into the electric field solution equation, perform finite element analysis on the global potential distribution, and obtain the global potential distribution of the current iteration step. Based on the anisotropic resistivity tensor of each node in the global potential distribution and spatial mapping database, the Joule heat density of each node in the global domain is calculated and processed to obtain the global Joule heat density distribution of the current iteration step. The global Joule heat density distribution is substituted into the temperature field solution equation as a heat source term. The global temperature field is then solved by finite element method to obtain the updated global temperature field. The updated global temperature field triggers the update of the resistivity tensor of the spatial mapping database. The solution loop is returned until the maximum node-wise difference between two adjacent global temperature fields is less than 0.5 Kelvin, and the converged global temperature field is obtained. The Joule heat density of each node in the slot corner node set is indexed and stored in all iteration steps according to the node number and iteration sequence number to obtain the Joule heat density tracking matrix. Then, the global average Joule heat density of each iteration step is stored sequentially according to the iteration sequence number to obtain the global average Joule heat density sequence.

[0029] Specifically, the electric field equation is a steady-state current continuity equation, which physically signifies global current conservation under conditions of no free charge accumulation. Specifically, it is formed by taking the divergence of the product of the conductivity tensor and the potential gradient to zero. The conductivity tensor is obtained by inverting the resistivity tensor corresponding to the current temperature of each node in the spatial mapping database. The potential gradient is the vector of the first-order partial derivatives of the potential field with respect to spatial coordinates. In the boundary conditions, fixed potential values ​​are applied to the contact surfaces of the two electrodes, and an electrically insulating boundary condition with zero normal current density is applied to the outer surface of the heater and all slot surfaces. After finite element discretization, a system of linear equations is formed with nodal potentials as unknowns. The coefficient matrix is ​​assembled from the conductivity tensors of each element through shape function integration. Solving this system of linear equations yields the global nodal potential distribution. The Joule heat density of each node is calculated by bilinearly contracting the potential gradient vector and the conductivity tensor of that node, with units of watts per cubic meter. The global Joule heat density distribution is the set of all nodal Joule heat density values.

[0030] The temperature field solution equation is a transient heat conduction equation. The product of mass density and specific heat capacity multiplied by the partial derivative of temperature with respect to time equals the divergence of the product of thermal conductivity tensor and temperature gradient plus Joule heat density. Joule heat density is substituted node by node as a volume heat source term. Under vacuum heating conditions, only radiation boundary conditions are retained for the outer surface boundary conditions. The radiation heat flux density is calculated using the Stefan-Boltzmann law and surface emissivity. After each iteration of the temperature field solution, the new temperature values ​​of each node in the updated global temperature field trigger a lookup in the spatial mapping database, reassembles the resistivity tensor of each node, and proceeds to the electric field solution of the next iteration. The convergence criterion is that the maximum node-by-node difference between two adjacent global temperature fields is less than 0.5 Kelvin. This value corresponds to the minimum resolvable temperature difference of isostatically pressed graphite within the engineering allowable temperature error range. When the value is lower than this, the impact of continuing iteration on the Joule heat density distribution is below the numerical noise level. The Joule heat density tracking matrix stores the Joule heat density value of each node in each iteration step using the node number in the slot corner node set as the row index and the iteration number as the column index. The global average Joule heat density sequence uses the iteration number as the index, and each element is the volume-weighted arithmetic mean of the Joule heat density of each node in the global domain in that iteration step. The weight is the ratio of the volume of the control volume corresponding to each node to the total volume of the global domain.

[0031] Figure 2This diagram illustrates the variation of the single-step self-enhancement ratio of each node in the slot corner node set under different fillet radii with the number of iteration steps in this embodiment of the application. The horizontal axis represents the number of iteration steps, and the vertical axis represents the average single-step self-enhancement ratio. The four curves correspond to the simulation results for fillet radii of 0.5mm, 1.0mm, 1.5mm, and 2.0mm, respectively. The shaded area represents the step number interval of the iteration convergence segment, and the dashed line represents the preset safety threshold of 1.50. The average value of each curve that tends to stabilize within the convergence segment is the hotspot self-enhancement index of each node under the corresponding fillet radius. The curve corresponding to a fillet radius of 0.5mm stabilizes at 2.41 within the convergence segment, exceeding the safety threshold, and the corresponding node is included in the dangerous corner node set. The curve corresponding to a fillet radius of 1.5mm stabilizes at 1.47 within the convergence segment, below the safety threshold, and satisfies the hotspot safety constraint.

[0032] Figure 3 This is a schematic diagram of the global Joule heat density distribution of the sawtooth cross-section of the graphite heater after the convergence of the electrothermal bidirectional iterative coupling solution in the embodiments of this application. The grayscale in the figure represents the normalized value of the Joule heat density, and the darker the color, the higher the Joule heat density. The white area is the location of the through slot, and the dark circular concentrated area is the current concentration hot spot at the corner of the slot. The hot spot location corresponds one-to-one with the spatial coordinates of the dangerous corner node in the set of corner nodes at the slot. The peak value of the Joule heat density in the hot spot area is about 2.4 to 2.5 times the global average value, which verifies the existence of the positive feedback effect between the anisotropic resistivity and the current concentration at the corner of the slot.

[0033] In one specific embodiment, step S3, which uses the ratio of the Joule heat density increment of each node in two adjacent iterations to the global average Joule heat density increment, includes: Based on the Joule heat density tracking matrix, the Joule heat density of each node in the slot corner node set is subtracted between two adjacent iteration steps to obtain the stepwise Joule heat density increment sequence of each node. Based on the global average Joule heat density sequence, the difference between the global average Joule heat density of two adjacent iteration steps is processed to obtain the stepwise global average Joule heat density increment sequence. Divide the Joule heat density increment of each node in each iteration step in the progressive Joule heat density increment sequence by the global average Joule heat density increment of the corresponding iteration step in the progressive global average Joule heat density increment sequence. Perform progressive ratio calculation on each node in the slot corner node set in all iteration steps to obtain the single-step self-enhancing ratio sequence of each node.

[0034] Specifically, in the Joule heat density tracking matrix, each row corresponds to a node in the slot corner node set, and each column corresponds to an iteration step. The matrix elements are the Joule heat density values ​​of that node in that iteration step. The difference between adjacent columns is performed row-by-row, i.e., the element in the (k+1)th column is subtracted from the element in the corresponding row of the kth column, yielding the Joule heat density increment of each node between the kth and (k+1)th iteration steps. The increments of all iteration steps are arranged by iteration number to form a progressive Joule heat density increment sequence for each node, with a length equal to the total number of iteration steps minus one. In the global average Joule heat density sequence, the difference between adjacent elements is performed, i.e., the (k+1)th element is subtracted from the kth element, yielding the global average Joule heat density increment for the kth iteration step. The increments of all iteration steps are arranged by iteration number to form a progressive global average Joule heat density increment sequence, with the same length as the progressive Joule heat density increment sequence for each node.

[0035] Each element of the single-step self-enhancing ratio sequence for each node is obtained by dividing the increment value of that node in the corresponding iteration step in the stepwise Joule heat density increment sequence by the global average increment value of the same iteration step in the stepwise global average Joule heat density increment sequence. When the absolute value of a certain iteration step in the stepwise global average Joule heat density increment sequence is less than the regularization term, the regularization term is used to replace that value in the division operation. The regularization term is set to one hundred millionth of a watt per cubic meter. This value is lower than the minimum resolvable value of Joule heat density under the operating conditions of the graphite heater and does not affect the ratio calculation result of the effective iteration steps. The physical meaning of the single-step self-enhancing ratio sequence for each node is the amplification factor of the rate of change of Joule heat density of that node relative to the global average rate of change in each iteration step. When this value is greater than 1 in a certain iteration step, it indicates that the growth rate of Joule heat density of that node in that step is faster than the global average level, reflecting that the positive feedback between current concentration and local resistivity change has produced a net enhancement effect on that node in that iteration step.

[0036] In the computational implementation, the starting step number of the iterative convergence segment was previously determined by taking the quartile of the convergence iteration steps, i.e., fixing the latter half of the total iteration steps as the iterative convergence segment. This method could obtain a stable mean of the ratio sequence within the convergence segment with low computational cost when the total number of iteration steps was relatively stable. However, as the actual number of convergence iteration steps corresponding to the global temperature field under different fillet radii becomes inconsistent, the method of fixing the latter half of the total steps as the starting step number of the convergence segment changes proportionally with the total number of iteration steps. This introduces a systematic bias related to the total number of iteration steps into the hotspot self-enhancing index calculated under different fillet radii, masking the true moment when the single-step self-enhancing ratio sequence of each node tends to stabilize. To eliminate this bias, this application changes the method of determining the starting step number of the iterative convergence segment from truncating it according to a fixed proportion of the total number of iteration steps to adaptive backtracking based on the standard deviation of the ratio sequence node by node. This ensures that the delineation of the convergence segment is determined only by the stability of the ratio sequence itself, and no longer systematically shifts with the difference in the total number of iteration steps.

[0037] In one specific embodiment, step S3 involves iteratively averaging the convergence segments of each node in the slot corner node set to obtain the hotspot self-reinforcing index of each node, including: Based on the actual number of convergence iterations corresponding to the convergence global temperature field, backtracking step by step from the last iteration step, for each candidate starting step number, calculate the maximum value of the node-by-node standard deviation of the single-step self-enhancement ratio of all nodes in the slot corner node set within the interval from the candidate starting step number to the last iteration step, and take the earliest candidate starting step number that satisfies the maximum value being less than the preset stability threshold and the interval length being no less than the preset minimum window number of steps as the starting number of the iteration convergence segment, thus obtaining the step number interval of the iteration convergence segment; Based on the step number interval of the iterative convergence segment, extract the single-step self-enhancement ratio subsequence of each node within the step number interval of the iterative convergence segment from the single-step self-enhancement ratio sequence of each node, and perform arithmetic mean processing on each value in the single-step self-enhancement ratio subsequence to obtain the hotspot self-enhancement index of each node. The hotspot self-reinforcement index is compared with a preset safety threshold node by node, and nodes whose hotspot self-reinforcement index exceeds the preset safety threshold are selected to obtain a set of dangerous corner nodes.

[0038] Specifically, the actual number of convergence iterations is determined by the iteration termination number corresponding to the convergence global temperature field, and is denoted as the total number of steps. The step number range for the iterative convergence segment is determined using the following adaptive method: for each candidate starting step number... ( from (decreasing progressively forward by -1), calculate the node-by-node standard deviation of the single-step self-enhancing ratio sequence of all nodes in the slot corner node set within the interval [k, K-1]:

[0039] in For nodes In the The single-step self-enhancement ratio of the iteration step, Let $i$ be the arithmetic mean of the single-step self-enhancement ratio of node $i$ in the interval [k, K-2]. Take the maximum standard deviation of all nodes. ,when And interval length When, continue to expand the candidate starting step number forward; when When the iteration stops backtracking, $k+1$ is taken as the starting number of the iteration convergence segment step number interval, and $K-1$ is taken as the ending number. A preset stability threshold is set. Set the value to 0.05, and preset the minimum window step size. Take 5. The value is based on the fact that when the standard deviation is less than 0.05, the gradual fluctuation amplitude of the single-step self-enhancement ratio of each node is less than 5% of its mean, which is sufficient to reflect the true local amplification factor near the steady state. A value of 5 ensures that the convergence segment contains a sufficient number of effective sampling steps to eliminate numerical noise. This adaptive criterion allows the determination of the convergence segment to be entirely determined by the actual stability of the ratio sequence during the iteration process, avoiding the systematic bias caused by the difference in convergence speed between different models due to fixed proportion truncation.

[0040] Subsequences are extracted from the single-step self-enhancing ratio sequence of each node according to the step number interval of the iteration convergence segment. The arithmetic mean of all elements in the subsequence is calculated to obtain the hotspot self-enhancing index of the node. This index is a dimensionless scalar, representing the average amplification factor of the Joule heat density growth rate relative to the global average growth rate during the iteration process near the steady state of the node.

[0041] The preset safety threshold is 1.50. This value is based on the fact that when the hot spot self-reinforcing index does not exceed 1.50, the peak Joule heat density at the corner does not exceed 3.0 times the average value of the entire region. This corresponds to the upper limit of the maximum acceptable local heat flux density deviation for isostatic graphite under long-term stable working conditions below 2000 degrees Celsius in engineering practice. Exceeding this multiple will significantly accelerate the initiation of intergranular cracks due to the thermal stress concentration of the graphite material at the corner. The hot spot self-reinforcing index of all nodes in the corner node set is compared with 1.50 node by node. Nodes with a hot spot self-reinforcing index strictly greater than 1.50 are identified as dangerous corner nodes. Their node numbers, spatial coordinates, and corresponding hot spot self-reinforcing index values ​​are extracted to form a dangerous corner node set. When the number of nodes in this set is zero, it indicates that all corner nodes under the current fillet radius meet the safety constraints. Otherwise, the degree of danger needs to be judged based on the difference between the maximum hot spot self-reinforcing index in the dangerous corner node set and the preset safety threshold, driving the reverse optimization of the fillet radius.

[0042] In one specific embodiment, step S4 includes: Graphite standard samples of the same batch with different fillet radii were subjected to thermal cycling tests under target working conditions. The number of thermal cycles when cracks first appeared on each sample was recorded. The data points of the hot spot self-reinforcing index and the number of thermal cycles corresponding to each sample were subjected to power function least squares fitting to obtain the mapping relationship between the hot spot self-reinforcing index and the number of failure thermal cycles. Based on the hotspot self-reinforcement index of each node in the dangerous corner node set, the maximum value of the hotspot self-reinforcement index of all nodes in the dangerous corner node set is taken to obtain the maximum hotspot self-reinforcement index of the entire domain. The maximum hotspot self-reinforcement index of the entire domain is input into the mapping relationship between the hotspot self-reinforcement index and the number of failure thermal cycles. Power function evaluation is then performed to obtain the predicted number of failure thermal cycles. The number of predicted failure thermal cycles is compared with the preset target number of thermal cycles. When the number of predicted failure thermal cycles is less than the preset target number of thermal cycles, the preset target number of thermal cycles is solved in reverse to obtain the target hotspot self-reinforcing index threshold. The preset safety threshold is replaced with the target hotspot self-reinforcing index threshold and the process is restarted. When the number of predicted failure thermal cycles is not less than the preset target number of thermal cycles, the current fillet radius and the number of predicted failure thermal cycles are output as the final design result.

[0043] Specifically, the target operating parameters for the thermal cycling test include the heating rate, peak temperature, holding time, and cooling rate. The heating rate is set at 50 degrees Celsius per minute, the peak temperature is consistent with the actual peak operating temperature of the designed graphite heater, the holding time is the actual holding time under operating conditions, and the cooling method is natural cooling with a cooling rate of approximately 20 degrees Celsius per minute. Four groups of standard graphite samples from the same batch were selected with fillet radii of 0.5 mm, 1.0 mm, 1.5 mm, and 2.0 mm, with at least three samples in each group. Each group of samples underwent thermal cycling tests under the aforementioned conditions. The cumulative number of thermal cycles completed when the first visible crack exceeding 0.1 mm in length appeared on the surface of each sample was recorded. The arithmetic mean of multiple samples in the same group was taken as the failure thermal cycle count for that group. For each group of samples, the corresponding hot spot self-reinforcing index was calculated according to the aforementioned steps. This index, along with the number of failure thermal cycles, formed four data point pairs. With the hot spot self-reinforcing index as the independent variable and the number of failure thermal cycles as the dependent variable, the fitting form was a power function. That is, the number of failure thermal cycles was equal to the coefficient A multiplied by the negative B power of the hot spot self-reinforcing index, where A and B are both positive real numbers. By taking the logarithm of the four data points, the problem was transformed into a linear regression problem. The values ​​of A and B were obtained by solving the least squares criterion, thus completing the calibration of the mapping relationship.

[0044] The maximum hotspot self-reinforcing index is obtained by taking the maximum value from the hotspot self-reinforcing indices of each node in the set of dangerous corner nodes. The basis for taking the maximum value is that the failure is dominated by the most dangerous node. Substituting the maximum value into the power function yields the predicted number of failure thermal cycles. This number represents the number of thermal cycles that the heater is expected to experience its first failure under the target operating condition under the current fillet radius design. The predicted number of failure thermal cycles is compared with the preset target number of thermal cycles, which is given by the design requirements. When the predicted number of failure thermal cycles is less than the preset target number of thermal cycles, the preset target number of thermal cycles is substituted into the power function, and the corresponding hotspot self-reinforcing index threshold is calculated in reverse. That is, the preset target number of thermal cycles is divided by the coefficient A and the negative B root is taken to obtain the target hotspot self-reinforcing index threshold. After replacing the preset safety threshold of 1.50 with this threshold, the selection step from the set of dangerous corner nodes is re-executed, driving the fillet radius to further increase until the predicted number of failure thermal cycles is not less than the preset target number of thermal cycles. At this point, the current fillet radius and the predicted number of failure thermal cycles are output as the final design result.

[0045] In one specific embodiment, the mapping relationship between the hot spot self-reinforcing index and the number of failure thermal cycles is determined as follows: thermal cycling tests are conducted on four groups of graphite standard samples from the same batch, with fillet radii of 0.5 mm, 1.0 mm, 1.5 mm, and 2.0 mm, respectively. The data points of the hot spot self-reinforcing index and the number of thermal cycles corresponding to each group of samples are subjected to power function least squares fitting. The minimum sum of squared fitting residuals is used as the convergence condition to obtain the power function coefficient group in the mapping relationship between the hot spot self-reinforcing index and the number of failure thermal cycles.

[0046] The process involves comparing the predicted number of failure thermal cycles with the preset target number of thermal cycles. This includes: based on the mapping relationship between the hotspot self-reinforcement index and the number of failure thermal cycles, the preset target number of thermal cycles is substituted in reverse to obtain the target hotspot self-reinforcement index threshold corresponding to the preset target number of thermal cycles; the difference between the target hotspot self-reinforcement index threshold and the maximum hotspot self-reinforcement index in the entire domain is calculated. When the difference is less than 0, the preset safety threshold is replaced with the target hotspot self-reinforcement index threshold and the process is re-executed; when the difference is not less than 0, the current fillet radius and the predicted number of failure thermal cycles are output as the final design result.

[0047] The selection of four fillet radii (0.5mm, 1.0mm, 1.5mm, and 2.0mm) is based on the feasible range of fillet radii in the actual processing of isostatic graphite heaters. 0.5mm is the lower limit of milling cutter processing, and 2.0mm is the upper limit of structure under the constraints of wall thickness and groove width. Four equally spaced sampling points ensure the interpolation accuracy of the power function fitting within this interval. The hot spot self-reinforcing index corresponding to each group of samples is calculated using the aforementioned method. It forms four data point pairs with the failure thermal cycle number of each group. The logarithm of both ends of the data point pairs is taken to transform the power function relationship into a linear relationship. The slope and intercept are solved using the linear least squares method. After inverse transformation, the power function coefficients A and B are obtained. The convergence condition is that the sum of squares of the fitting residuals reaches its minimum value in the logarithmic space, that is, the sum of squares of the differences between the logarithmic transformation values ​​of all four groups of data points and the predicted values ​​of the fitted straight line reaches its minimum. At this time, the values ​​of A and B are uniquely determined, and the calibration of the power function coefficient group is completed.

[0048] The specific operation of the reverse substitution process is as follows: divide the preset target number of thermal cycles by the coefficient A, take the negative B root of the quotient (i.e., take the natural logarithm of the quotient, divide by negative B, and then take the exponent), to obtain the target hotspot self-reinforcing index threshold. The physical meaning of this threshold is the maximum allowable hotspot self-reinforcing index that can just meet the preset target number of thermal cycles under the current mapping relationship. The difference between the target hotspot self-reinforcing index threshold and the global maximum hotspot self-reinforcing index is calculated. The difference is the target hotspot self-reinforcing index threshold minus the global maximum hotspot self-reinforcing index. When the difference is less than 0, it indicates that the current global maximum hotspot self-reinforcing index exceeds the threshold for meeting the lifetime requirement. After replacing the preset safety threshold of 1.50 with the target hotspot self-reinforcing index threshold, the process is re-executed from the dangerous corner node set screening stage, driving the fillet radius to increase until the global maximum hotspot self-reinforcing index falls within the target hotspot self-reinforcing index threshold. When the difference is not less than 0, the global maximum hotspot self-reinforcing index under the current fillet radius has met the lifetime constraint. The current fillet radius and the predicted number of failure thermal cycles are output as the final design result.

[0049] 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 the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features. Such modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. A simulation design method for a graphite heater, characterized in that, The method includes: Step S1: Perform mesh refinement processing on the slot corner area of ​​the graphite heater, extract the nodes belonging to the slot corner area in the refined mesh, and obtain the slot corner node set; Step S2: Assign the anisotropic resistivity tensor to each node in the set of slot corner nodes, and establish a spatial mapping database of the resistivity of each node as a function of temperature. Step S3: Substitute the spatial mapping database into the electrothermal bidirectional iterative coupling solution. In each iteration, extract the Joule heat density of each node in the slot corner node set. Using the ratio of the Joule heat density increment of each node in two adjacent iterations to the global average Joule heat density increment, perform iterative convergence segment averaging on each node in the slot corner node set to obtain the hotspot self-enhancing index of each node. Step S4: Substitute the hot spot self-reinforcing index into the pre-calibrated mapping relationship between the hot spot self-reinforcing index and the number of failure thermal cycles to obtain the predicted number of failure thermal cycles for the graphite heater.

2. The simulation design method for a graphite heater according to claim 1, characterized in that, Step S1 includes: A parametric three-dimensional geometric model is established for the cylindrical body of the graphite heater. The parametric three-dimensional geometric model includes the outer diameter, wall thickness, total axial height, through slots alternately opened from both ends of the cylinder along the axial direction, the spacing between adjacent slots, the width of a single slot, and the radius of the fillet at the corner of the slot opening. Based on the parametric three-dimensional geometric model, the entire domain is divided into a hexahedral structured mesh. For the corner region of the slot, no less than eight layers of mesh are set in the radial direction with the center of the corner arc as the reference. For the straight section of the cylinder, a sparse mesh is set to obtain the global mesh model. Based on the global mesh model, all nodes are classified according to whether their location belongs to the corner area of ​​the slot, resulting in a set of corner nodes. For each node in the slot corner node set, record its spatial coordinates, slot number, corner azimuth angle, and anisotropic principal direction vector to obtain the slot corner node set carrying node position attributes.

3. The simulation design method for a graphite heater according to claim 2, characterized in that, Step S2 includes: Based on the anisotropic principal direction vector recorded by each node in the slot corner node set carrying node position attributes, tensor assembly processing is performed on the in-plane resistivity component and the thickness resistivity component of each node to obtain the initial anisotropic resistivity tensor of each node. The measured high-temperature resistivity data of the same batch of graphite material samples under the four-probe method were processed by piecewise polynomial least squares fitting of the resistivity in the in-plane direction and the resistivity in the thickness direction, respectively, to obtain the temperature fitting coefficient group of resistivity in the in-plane direction and the anisotropy ratio in the thickness direction. Based on the in-plane resistivity temperature fitting coefficient set and the thickness direction anisotropy ratio, temperature-dependent assignment is performed on each component of the initial anisotropic resistivity tensor to obtain the temperature-dependent anisotropic resistivity tensor of each node. The spatial coordinates of each node in the slot corner node set carrying node position attributes are indexed and bound to the temperature-dependent anisotropic resistivity tensor to obtain a spatial mapping database of the resistivity of each node as a function of temperature.

4. The simulation design method for a graphite heater according to claim 1, characterized in that, In step S3, the spatial mapping database is substituted into the electrothermal bidirectional iterative coupling solution. In each iteration, the Joule heat density of each node in the slot corner node set is extracted, including: Substitute the anisotropic resistivity tensor corresponding to the current temperature of each node in the spatial mapping database into the electric field solution equation, perform finite element analysis on the global potential distribution, and obtain the global potential distribution of the current iteration step. Based on the global potential distribution and the anisotropic resistivity tensor of each node in the spatial mapping database, the Joule heat density of each node in the global domain is calculated and processed to obtain the global Joule heat density distribution of the current iteration step. The global Joule heat density distribution is substituted as a heat source term into the temperature field solution equation. The global temperature field is then solved by finite element method to obtain the updated global temperature field. The updated global temperature field is then used to trigger the resistivity tensor update of the spatial mapping database. The solution loop is returned until the maximum node-wise difference between two adjacent global temperature fields is less than 0.5 Kelvin, thus obtaining the converged global temperature field. The Joule heat density of each node in the slot corner node set is indexed and stored in all iteration steps according to the node number and iteration sequence number to obtain the Joule heat density tracking matrix. Then, the global average Joule heat density of each iteration step is stored sequentially according to the iteration sequence number to obtain the global average Joule heat density sequence.

5. The simulation design method for a graphite heater according to claim 4, characterized in that, Step S3, which uses the ratio of the Joule heat density increment of each node in two adjacent iterations to the global average Joule heat density increment, includes: Based on the Joule heat density tracking matrix, the Joule heat density of each node in the slot corner node set is subtracted between two adjacent iteration steps to obtain the stepwise Joule heat density increment sequence of each node. Based on the global average Joule heat density sequence, the difference between the global average Joule heat density of two adjacent iteration steps is processed to obtain the stepwise global average Joule heat density increment sequence. Divide the Joule heat density increment of each node in the progressive Joule heat density increment sequence of each node in each iteration step by the global average Joule heat density increment of the corresponding iteration step in the progressive global average Joule heat density increment sequence, and perform progressive ratio calculation on each node in the slot corner node set in all iteration steps to obtain the single-step self-enhancing ratio sequence of each node.

6. The simulation design method for a graphite heater according to claim 5, characterized in that, In step S3, the nodes in the slot corner node set are subjected to iterative convergence segment averaging to obtain the hotspot self-reinforcing index of each node, including: Based on the actual number of convergence iterations corresponding to the convergence global temperature field, backtracking step by step from the last iteration step, for each candidate starting step number, calculate the maximum value of the node-by-node standard deviation of the single-step self-enhancement ratio of all nodes in the slot corner node set within the interval from the candidate starting step number to the last iteration step, and take the earliest candidate starting step number that satisfies the maximum value being less than the preset stability threshold and the interval length being no less than the preset minimum window number of steps as the starting number of the iteration convergence segment, thus obtaining the step number interval of the iteration convergence segment; Based on the step number interval of the iterative convergence segment, extract the single-step self-enhancement ratio subsequence of each node within the step number interval of the iterative convergence segment from the single-step self-enhancement ratio sequence of each node, and perform arithmetic mean processing on each value in the single-step self-enhancement ratio subsequence to obtain the hotspot self-enhancement index of each node. The hotspot self-reinforcement index is compared with a preset safety threshold node by node, and nodes whose hotspot self-reinforcement index exceeds the preset safety threshold are selected to obtain a set of dangerous corner nodes.

7. The simulation design method for a graphite heater according to claim 6, characterized in that, Step S4 includes: Graphite standard samples of the same batch with different fillet radii were subjected to thermal cycling tests under target working conditions. The number of thermal cycles when cracks first appeared on each sample was recorded. The data points of the hot spot self-reinforcing index and the number of thermal cycles were subjected to power function least squares fitting to obtain the mapping relationship between the hot spot self-reinforcing index and the number of failure thermal cycles. Based on the hotspot self-reinforcement index of each node in the set of dangerous corner nodes, the maximum value of the hotspot self-reinforcement index of all nodes in the set of dangerous corner nodes is taken to obtain the maximum hotspot self-reinforcement index of the entire domain. The global maximum hotspot self-reinforcement index is input into the mapping relationship between the hotspot self-reinforcement index and the number of failure thermal cycles, and power function evaluation is performed to obtain the predicted number of failure thermal cycles. The predicted number of failure thermal cycles is compared with the preset target number of thermal cycles. When the predicted number of failure thermal cycles is less than the preset target number of thermal cycles, the preset target number of thermal cycles is solved in reverse to obtain the target hotspot self-reinforcement index threshold. The preset safety threshold is replaced with the target hotspot self-reinforcement index threshold and the process is restarted. When the predicted number of failure thermal cycles is not less than the preset target number of thermal cycles, the current fillet radius and the predicted number of failure thermal cycles are output as the final design result.

8. The simulation design method for a graphite heater according to claim 7, characterized in that, The mapping relationship between the hotspot self-enhancing index and the number of failure thermal cycles is defined in the following way: Thermal cycling tests were conducted on four groups of graphite standard samples from the same batch, with fillet radii of 0.5 mm, 1.0 mm, 1.5 mm, and 2.0 mm, respectively. The hot spot self-reinforcing index and the number of thermal cycles corresponding to each group of samples were subjected to power function least squares fitting. The minimum sum of squared fitting residuals was used as the convergence condition to obtain the power function coefficient group in the mapping relationship between the hot spot self-reinforcing index and the number of failure thermal cycles.

9. The simulation design method for a graphite heater according to claim 7, characterized in that, The process of comparing the predicted number of failure thermal cycles with the preset target number of thermal cycles includes: Based on the mapping relationship between the hotspot self-reinforcement index and the number of failure thermal cycles, the preset target number of thermal cycles is substituted in reverse to obtain the target hotspot self-reinforcement index threshold corresponding to the preset target number of thermal cycles. The difference between the target hotspot self-reinforcement index threshold and the global maximum hotspot self-reinforcement index is calculated. When the difference is less than 0, the preset safety threshold is replaced with the target hotspot self-reinforcement index threshold and the process is restarted. When the difference is not less than 0, the current fillet radius and the number of failure thermal cycle predictions are output as the final design result.

Images