G-code based comprehensive thermal risk assessment and local thermal compensation method for 3D printing

Abstract

This invention presents a comprehensive thermal risk assessment and local thermal compensation method for 3D printing based on G-code. It uses the G-code of material extrusion 3D printing as the processing object, constructing the manufacturing path state and printing time difference of each printing layer by parsing the G-code. Based on historical material sampling points within the current printing layer, a planar thermal risk value is calculated. This planar thermal risk value is then continuously corrected according to the local coverage of multiple layers below the current printing layer to obtain a comprehensive thermal risk value. This comprehensive thermal risk value is mapped to a thermal compensation coefficient, and the local extrusion amount and fan control amount are corrected based on the thermal compensation coefficient, outputting corrected G-code. This invention simultaneously considers the path execution order specified by the G-code, printing time difference, directional heat dissipation differences, and the material coverage of multiple layers below. Furthermore, the assessment results can be directly converted into local process parameters and written into the G-code.

Description

Technical Field

[0001] This invention relates to the field of additive manufacturing technology, specifically to a method for comprehensive thermal risk assessment and local thermal compensation in 3D printing based on G-code. Background Technology

[0002] 3D printing typically begins by slicing a 3D model into layers using slicing software, generating G-code that executes layer by layer in a predetermined order. The G-code contains the printing path, execution order, and corresponding process parameters for each layer. In material extrusion 3D printing, molten material is extruded and deposited sequentially along the path specified by the G-code, and the thermal state of local areas continuously changes as the path is executed. The historical thermal impact on the current printing location depends not only on the presence of nearby printed material but also on the execution order of adjacent paths, printing timestamps, local path direction, the main cooling direction of the equipment, and the material coverage of the layers already printed below the current layer.

[0003] Existing methods for assessing and controlling the thermal state of 3D printing mainly include online monitoring methods based on temperature sensors or infrared imaging, numerical calculation methods for temperature fields based on finite element analysis, machine learning prediction methods based on historical samples, and approximate judgment methods based on local geometric distance, material distribution density, or fixed process thresholds. Online monitoring methods can provide feedback adjustments based on the actual printing status; numerical calculation methods for temperature fields can obtain relatively complete information on the thermal process; methods based on geometry or empirical rules have relatively low computational complexity and are easy to integrate with slicing software or G-code post-processing programs.

[0004] However, existing methods still have certain limitations when directly using G-code for pre-printing local thermal risk assessment and compensation. Online monitoring methods typically require additional temperature detection devices to obtain detection information during the printing process; finite element temperature field calculations usually have high requirements for material parameters, boundary conditions, and computational resources; machine learning-based prediction results are affected by training samples and applicable operating conditions. Methods with lower computational overhead, such as geometric distance, local density, or fixed threshold methods, often lack a comprehensive characterization of the actual path execution timing, historical heat input attenuation, and directional heat dissipation differences. For locations with similar local geometric distances but different printing timings or heat dissipation conditions, the above methods may yield similar risk assessment results.

[0005] Furthermore, even if a local thermal risk assessment result is obtained, if the correspondence between the assessment result and the specific printing position in the original G-code cannot be established, and the local process parameter correction amount and the writing position and method of the control instructions cannot be further determined, the assessment result is still difficult to directly use for printing equipment execution. Especially when the compensation amount at adjacent positions changes continuously, it is also necessary to deal with problems such as sudden changes in local process parameters, excessive control instructions, and the impact of parameter modifications on subsequent execution timing. Therefore, a complete processing method is needed that can construct printing timing information and path space information from G-code, calculate the comprehensive thermal risk value, convert the comprehensive thermal risk value into local process parameter correction amounts, and output corrected G-code. Summary of the Invention

[0006] To address the problem that existing methods for assessing local thermal risks before printing cannot simultaneously consider the path execution order specified by G-code, printing time difference, directional heat dissipation differences, and the coverage of continuous multi-layer materials below, and that the assessment results cannot be directly converted into local process parameters and written into G-code, this invention provides a comprehensive thermal risk assessment and local thermal compensation method for 3D printing based on G-code.

[0007] This invention uses the G-code of material extrusion 3D printing as the processing object. By parsing the G-code, it constructs the manufacturing path state and printing time difference of each printing layer. Based on the historical material sampling points in the current printing layer, it calculates the planar thermal risk value. Then, it continuously corrects the planar thermal risk value according to the local coverage of the continuous multiple layers below the current printing layer to obtain a comprehensive thermal risk value. The comprehensive thermal risk value is mapped to a thermal compensation coefficient, and the local extrusion amount and fan control amount are corrected according to the thermal compensation coefficient, and the corrected G-code is output.

[0008] To solve the above-mentioned technical problems, the present invention adopts the following technical solution: In a first aspect, the present invention provides a 3D printing local thermal compensation method based on G-code, comprising the following steps: S1: Analyze the G-code of material extrusion 3D printing to identify the manufacturing segment where material deposition occurs and the auxiliary motion segment where material deposition does not occur; discretize the manufacturing segment to construct a set of manufacturing paths for each printing layer, and set printing time markers for material sampling points according to the length, feed speed and execution order of each motion segment; wherein, the auxiliary motion segment is not written into the material distribution map, but its execution time is included in the printing time markers.

[0009] S2: For the current sampling point to be evaluated, select historical material sampling points from the current printing layer whose printing timestamps are earlier than the current sampling point to be evaluated and are located within a preset time window; determine the time decay weight based on the printing time difference between the historical material sampling points and the current sampling point to be evaluated; determine the local direction reference for the current sampling point to be evaluated, wherein, when the main direction of local cooling of the equipment can be obtained, the main direction of local cooling of the equipment is used as the local direction reference, and when the main direction of local cooling of the equipment cannot be obtained, the local path tangential direction of the manufacturing section where the current sampling point to be evaluated is located is used as the local direction reference; determine the directional thermal influence weight based on the planar distance and relative orientation between the historical material sampling points and the current sampling point to be evaluated, and perform weighted accumulation of the thermal influence of the historical material sampling points to obtain the planar thermal risk value of the current sampling point to be evaluated.

[0010] S3: Statistically analyze the local material coverage of the current sampling point under the current printing layer in multiple consecutive layers, and set the depth attenuation weight according to the distance between each lower printing layer and the current printing layer to obtain the interlayer heat transfer score; continuously correct the planar thermal risk value according to the interlayer heat transfer score to obtain the comprehensive thermal risk value.

[0011] S4: Normalize the comprehensive thermal risk value and obtain the thermal compensation coefficient through continuous mapping; smooth the thermal compensation coefficient along the execution order of G code; based on the smoothed thermal compensation coefficient, correct at least one of the extrusion amount and fan control amount of the corresponding manufacturing section; reconstruct the G code based on the correction result and output the corrected G code.

[0012] Secondly, the present invention provides a comprehensive thermal risk assessment method for 3D printing based on G-code, the assessment method comprising the following: S1: Obtain the G-code for material extrusion 3D printing; parse the motion segments in each printing layer according to the execution order of the G-code; identify the manufacturing segments where material deposition occurs and the auxiliary motion segments where material deposition does not occur; discretize the manufacturing segments to obtain material sampling points arranged according to the execution order of the G-code; and construct the manufacturing path set and material distribution map for each printing layer; determine the cumulative execution time based on the length, feed speed, and execution order of each motion segment; and set printing time markers for the material sampling points; wherein, the auxiliary motion segments do not form material occupation, but their execution time is included in the cumulative execution time; S2: Select the material sampling points sequentially as the current sampling points to be evaluated. In the current printing layer, select historical material sampling points located in the local planar neighborhood of the current sampling point to be evaluated, whose printing timestamps are earlier than the current sampling point to be evaluated, and whose printing time difference is within a preset time window. Determine the time decay weight based on the printing time difference. Determine the local orientation reference of the current sampling point to be evaluated. Determine the directional thermal influence weight based on the planar distance, relative orientation, and local orientation reference between the historical material sampling points and the current sampling point to be evaluated. Then, perform a weighted summation based on the material occupancy state, time decay weight, and directional thermal influence weight of the historical material sampling points to obtain the planar thermal risk value of the current sampling point to be evaluated. S3: Simultaneously, based on the occupied status of the printed material at the current sampling point to be evaluated, determine the interlayer heat transfer status of the area and determine the interlayer heat transfer score; use the interlayer heat transfer score to correct the planar thermal risk value and obtain the comprehensive thermal risk value.

[0013] Furthermore, the method for determining the local direction reference is as follows: when the main direction of local cooling of the equipment can be obtained, the main direction of local cooling of the equipment is used as the local direction reference; when the main direction of local cooling of the equipment cannot be obtained, the local path tangential direction of the manufacturing segment where the current sampling point to be evaluated is located is used as the local direction reference.

[0014] Compared with the prior art, the present invention has the following beneficial effects: 1. This invention constructs a material distribution map and printing time stamp by parsing G-codes, and includes the execution time of the auxiliary motion segment in the printing time stamp, so that the thermal risk assessment corresponds to the actual path execution sequence specified by the G-codes.

[0015] 2. In planar thermal risk assessment, this invention comprehensively considers the printing time difference, relative position, and directional heat dissipation differences of historical material sampling points, and prioritizes the main local cooling direction of the equipment as the local directional reference for the directional thermal influence core. When the main local cooling direction of the equipment cannot be obtained, the local path tangential direction of the manufacturing section where the current sampling point to be assessed is located is used as an alternative local directional reference. The thermal influence of historical material sampling points is weighted and accumulated, and then corrected using the material coverage of multiple consecutive layers below the current printing layer. This invention can distinguish local locations with similar geometric distances but different path thermal histories and heat dissipation conditions.

[0016] 3. This invention employs continuous interlayer correction and thermal compensation mapping, and smooths the thermal compensation coefficient along the G-code execution order, which can reduce the impact of sudden changes in compensation amount on the printing process.

[0017] 4. This invention can convert the comprehensive thermal risk assessment results into local extrusion amount and fan control amount, and output corrected G-code, which can be used as a G-code post-processing method to be combined with existing material extrusion 3D printing processes. Attached Figure Description

[0018] Figure 1 A schematic diagram of the comprehensive thermal risk assessment and local thermal compensation method.

[0019] Figure 2 Material distribution diagram of the current layer to be evaluated (layer k).

[0020] Figure 3 Material distribution diagram of the first layer (k-1 layer) below the current layer to be evaluated.

[0021] Figure 4 Material distribution diagram of the second layer (k-2) below the current layer to be evaluated.

[0022] Figure 5 Material distribution diagram of the third layer (k-3) below the current layer to be evaluated.

[0023] Figure 6 Material distribution diagram of the fourth layer (k-4th layer) below the current layer to be evaluated.

[0024] Figure 7 Material distribution diagram of the fifth layer (k-5th layer) below the current layer to be evaluated. Detailed Implementation

[0025] The present invention will be further explained below with reference to the embodiments and accompanying drawings, but this is not intended to limit the scope of protection of this application.

[0026] Example 1 This embodiment takes material extrusion 3D printing as an example. The input includes a G-code file generated by slicing software, as well as necessary process parameters such as target linewidth and discrete sampling step size. The method flow is as follows: Figure 1 As shown, firstly, the manufacturing path and execution sequence of each printing layer are extracted from the G code to construct a set of manufacturing paths and printing time stamps; then, the planar thermal risk value, interlayer heat transfer score, comprehensive thermal risk value and thermal compensation coefficient are calculated in sequence, and the local process parameters are corrected and the G code is reconstructed based on the thermal compensation coefficient.

[0027] I. G-code Parsing and Basic Data Construction: The G-code output by the slicing software is parsed line by line. The currently valid coordinate mode, extrusion mode, feed rate, position coordinates, and extrusion amount are updated according to the actual execution order of the G-code. Each motion segment is determined based on the motion relationship between adjacent positions, and the starting coordinates, ending coordinates, corresponding layer number, feed rate, extrusion amount changes, and execution order of each motion segment are obtained.

[0028] Based on the planar displacement and extrusion rate changes of each motion segment, the motion segment where actual material deposition occurs is identified and denoted as the manufacturing segment; the idle stroke segment, retraction segment, lifting segment, transposition segment, and other motion segments that do not result in actual material deposition are denoted as auxiliary motion segments. For the k-th printing layer, the set of manufacturing paths consisting of all manufacturing segments in that printing layer is denoted as... Auxiliary motion segments are not included in the manufacturing path set. It does not create material occupation, but its execution time is still included in the cumulative execution time of G-code.

[0029] Set along the manufacturing path According to the preset step size Discrete sampling is performed to obtain material sampling points arranged according to the execution order of the G-code. Each material sampling point retains its planar coordinates, corresponding layer number, manufacturing segment, execution order, and local path direction. The preset step size... Based on the target line width In accordance with the predetermined calculation resolution setting, in this embodiment, a resolution not greater than the target line width is selected. The positive value. The target linewidth. This refers to the target material forming width corresponding to a single manufacturing path within the current printing layer plane, which is determined by the slicing parameters.

[0030] A two-dimensional computational grid is established for each printing layer. For any grid center in the k-th layer of the two-dimensional computational grid... According to its set of manufacturing paths The shortest distance is used to construct the material distribution map of the k-th layer according to the following formula. : (1) in, Indicates the center of the pixel grid ( ) to the Layer manufacturing path set The shortest distance; This indicates the presence of deposited material at that location. This indicates that there is no deposited material at this location.

[0031] To ensure that material sampling points near the edge of the manufacturing area can establish a complete local planar neighborhood, when establishing the two-dimensional computational grid, a blank boundary region with a width not less than the effective radius of the directional thermally affected core is set outside the manufacturing path; since no material deposition occurs in the grid points within the blank boundary region, their material occupancy state is set to 0.

[0032] Simultaneously, based on the motion quantity, feed rate, and execution sequence of each manufacturing segment and auxiliary motion segment, the start time, end time, and duration of each motion segment are calculated sequentially. For discrete material sampling points on the manufacturing segment, a linear interpolation is performed between the start time and end time of the manufacturing segment based on the relative position of the material sampling point on its respective manufacturing segment to obtain the printing time stamp of the material sampling point. Although the auxiliary motion section does not generate material sampling points, its execution time is included in the start time of the subsequent manufacturing section to reflect the actual cooling time interval formed by the idle stroke, retraction, lifting and transposition processes.

[0033] Therefore, each material sampling point corresponds to its layer number, planar coordinates, printing time, path direction, and material occupancy status.

[0034] II. Selection of Historical Material Sampling Points and Time Decay Based on the aforementioned material distribution map and printing timestamps, any material sampling point P in the k-th layer is selected as the current sampling point to be evaluated. Within the local planar neighborhood of the current sampling point P, points with printing timestamps earlier than [previous timetamp] are selected. The material sampling points Q constitute the candidate historical material sampling point set.

[0035] Let the printing timestamps of the current sampling point P to be evaluated and the historical material sampling point Q be respectively... and The time difference between their printing is: (2) in, The printing time of the current sampling point P to be evaluated. The printing time for historical material sampling point Q. When When Q did not complete deposition earlier than P, it does not participate in the thermal impact calculation at the current point; when When Q is included in the historical material sampling point set; when At that time, it is considered that the residual thermal influence of the sampling point on the current point is already weak, and it is not included in the calculation. This is a preset time window; you can set it yourself according to the actual situation.

[0036] For historical material sampling points that meet the time condition, a time decay weight is used to represent the decrease in thermal effect as the time interval increases: (3) in, This is the time decay coefficient. The shorter the time interval, the lower the decay rate. The larger the value; the longer the time interval, The smaller.

[0037] III. Directional Thermally Affected Core Construction After obtaining the set of historical material sampling points, the directional thermal influence of each historical material sampling point on the current sampling point to be evaluated is further calculated. Let the historical material sampling points be... Relative to the current sampling point to be evaluated Planar offset vector for: (4) This embodiment uses a directional thermally affected kernel for representation. right Thermal influence weight: (5) in, This is the kernel amplitude coefficient, used to adjust the overall weight of the directional thermally affected kernel. This is the covariance matrix, used to control the extent of thermal influence in different directions. (Superscript) This represents the transpose operation. The covariance matrix can be represented as: (6) in, It is a two-dimensional rotation matrix. This refers to the local orientation angle near the current sampling point to be evaluated. The influence of the major axis direction on the scale, The minor axis direction influences the scale. The major axis direction represents the direction where heat is relatively difficult to dissipate, while the minor axis direction represents the direction where the thermal effect decays faster. The local direction angle is determined based on the main local cooling direction of the equipment, and is taken as the angle between the axis containing the main local cooling direction of the equipment and the positive direction of the global X-axis; when the equipment does not have a local cooling direction or the main local cooling direction of the equipment cannot be obtained, the angle is taken as the angle between the local path tangent direction of the manufacturing section where the current sampling point to be evaluated is located and the positive direction of the global X-axis.

[0038] When it is necessary to express the difference in heat dissipation on both sides of the current point, Different values ​​can be taken on both sides of the local principal direction. For example, when the historical material sampling point is located on the positive side of the local principal direction, different values ​​can be taken. When located on the negative side of the local principal direction, take The positive and negative sides are determined based on the sign of the projection component of the planar offset vector of the historical material sampling point relative to the current sampling point to be evaluated in the local principal direction. and It can be preset according to the equipment cooling direction and path direction. This process is used to enable the directional thermal influence core to reflect the differences in local thermal influence in different directions.

[0039] IV. Planar Thermal Risk Value calculate After obtaining the time decay weight and directional thermal influence kernel, the historical material sampling points are weighted and accumulated to obtain the planar thermal risk value of the current sampling point P to be evaluated: (7) in, This indicates that the sample point P to be evaluated is located within the local planar neighborhood and satisfies the following conditions: A collection of historical material sampling points; Sampling points for historical materials The material occupies a state; As a time decay weight; This represents the directional thermal influence nucleus. From this formula, it can be seen that only historical material sampling points that have already been deposited, printed earlier than the current point, and are located within a local planar neighborhood will contribute to the planar thermal risk value of the current point. When no historical material sampling point meets the printing time and spatial neighborhood conditions, the planar thermal risk value of the current sampling point to be evaluated is set to zero.

[0040] V. Interlayer heat transfer score and overall thermal risk value Based solely on the planar thermal risk value of the current layer, it is difficult to reflect the impact of the material below the current point on heat transfer and support. Therefore, this embodiment further statistically analyzes the material coverage within the local statistical window corresponding to the N consecutive layers below the current sampling point. When the actual number of printed layers below the current printed layer is less than the preset number N, only the actual printed layers below are statistically analyzed. The depth attenuation weight can be normalized according to the actual number of layers involved in the statistics. No interlayer correction is performed for the first layer.

[0041] Suppose the current sampling point to be evaluated is located in the k-th layer, and the local coverage of the n-th layer below it is... The depth decay weight is Then the interlayer heat transfer score for: (8) in, Indicates the first The local material coverage rate within the local statistical window corresponding to the current sampling point to be evaluated in the layer is obtained by the ratio of the number of grids with a material occupancy state of 1 within the local statistical window to the total number of grids within the local statistical window. This represents the depth decay weight of the nth layer below. Layers closer to the current layer have a greater impact on heat transfer at the current point; layers farther away have a smaller impact. The depth decay weight can be calculated using: (9) in, This is the depth attenuation coefficient, and .

[0042] Let the risk gain coefficient be Based on the interlayer heat transfer score of the current sampling point to be evaluated, the planar thermal risk value is continuously corrected. A lower interlayer heat transfer score indicates less local coverage of the formed material below the current sampling point, a relatively weaker heat transfer path to the downward solid material, and a greater degree of correction applied to the planar thermal risk value. Interlayer heat transfer correction amount. for: (10) The higher the interlayer heat transfer score, the smaller the interlayer heat transfer correction; when the interlayer heat transfer score reaches its upper limit, the interlayer heat transfer correction is zero. (Comprehensive thermal risk value) for: (11) in, This represents the overall thermal risk value of the sampling points currently to be evaluated.

[0043] VI. Generation of Thermal Compensation Coefficient The range of the overall thermal risk value may differ depending on the model and the printing layer. To facilitate unified control in the future, the overall thermal risk value is normalized: (12) in, The upper bound for normalization; This is the normalized thermal risk index.

[0044] Subsequently, the normalized thermal risk index is input into the smoothing mapping function to obtain the thermal compensation coefficient: (13) in, This represents the thermal compensation coefficient corresponding to the current sampling point to be evaluated. The maximum compensation strength; The mapping slope parameter; This is to compensate for the risk threshold that begins to increase significantly. The higher the thermal risk index, the larger the thermal compensation coefficient; the lower the thermal risk index, the smaller the thermal compensation coefficient.

[0045] To reduce abrupt changes in the thermal compensation coefficients between adjacent material sampling points, a one-dimensional moving average can be applied to the thermal compensation coefficient sequence along the G-code execution order. Let the material sampling points arranged in the G-code execution order be... The moving average window length is an odd number. Then the material sampling point The thermal compensation coefficient after smoothing is denoted as And calculate according to the following formula: (14) Where L is the number of material sampling points included in the moving average window. For material sampling points located near the start or end of the path and unable to form a complete moving average window, L can be reduced accordingly. Subsequent local process parameter corrections use a smoothed thermal compensation coefficient. When no smoothing is performed, let .

[0046] VII. Local Process Parameter Correction and G-code Reconstruction After obtaining the thermal compensation coefficient, it is mapped to local process parameters. This embodiment uses extrusion rate and fan control parameters as examples for illustration.

[0047] If we map the thermal compensation coefficient to the extrusion rate, let the original extrusion rate increment be... The maximum extrusion reduction ratio is The corrected increment of extrusion amount is: (15) When the overall thermal risk at the current point is high, The increased local extrusion amount is appropriately reduced after correction, thereby reducing local material accumulation and heat input. If the heat compensation coefficient is mapped to the fan control value, let the reference fan control value be... The maximum fan control value is The corrected fan control variable is: (16) When the overall thermal risk at the current point is high, As the fan speed increases, the fan control parameters are correspondingly increased, thereby enhancing localized cooling.

[0048] For situations where different locations within the same manufacturing section correspond to different thermal compensation coefficients, the manufacturing section is divided into several sub-segments along its direction. The representative value of the thermal compensation coefficient for each sub-segment is used for parameter correction. During G-code reconstruction, the original path geometry remains unchanged; only the extrusion rate or fan control commands within the manufacturing section are locally rewritten.

[0049] After completing the one-dimensional moving average processing described above, to avoid excessive control commands, a minimum difference threshold for parameter updates is further set. New control commands are only generated when the difference between the old and new parameters reaches or exceeds this threshold; otherwise, the previous control parameters are used. For the idle stroke segment, retraction segment, lifting segment, and transposition segment, no material deposition-type thermal compensation correction is applied, but their execution time is still retained in the printing time series for the aforementioned printing time difference calculation.

[0050] Thus far, this embodiment has completed the process from G-code parsing, comprehensive thermal risk value calculation, thermal compensation coefficient generation to correcting G-code output.

[0051] Example 2 This embodiment uses the material extrusion additive manufacturing process as an example to construct a six-layer test model. For example... Figures 2 to 7 As shown in the figure, white sampling points represent sampling points where the material occupancy state is 0, gray sampling points represent sampling points where the material occupancy state is 1, and black point A is the sampling point to be evaluated in this example.

[0052] The test model uses a 15×15 two-dimensional grid with a grid spacing of [missing information]. The grid coordinates are denoted as x=1, 2, ..., 15 from left to right and y=1, 2, ..., 15 from top to bottom; layer 6 is designated as the current layer k to be evaluated; layers 5 through 1 correspond to layers (k-1) through (k-5) below, respectively. The current sampling point to be evaluated is point A, and its grid coordinates are: The grid coordinates (x, y) and the actual printed coordinates (X, Y) satisfy the following relationship: .

[0053] This layer has 27 material sampling points, and the printing order is as follows: (4,5)→(5,5)→…→(11,5)→(11,6)→(11,7)→(11,8)→(10,8)→(9,8)→A(8,8)→(7,8)→(6,8)→(5,8)→(5,9)→(5,10)→(5,11)→(6,11)→…→(12,11).

[0054] The continuous manufacturing path of the (k-1)th layer below is (8,8)→(8,9), with 2 material sampling points; the continuous manufacturing path of the (k-2)th layer below is (8,8)→(9,8)→(9,9)→(8,9), with 4 material sampling points; the continuous manufacturing path of the (k-3)th layer below is (9,7)→(9,8)→(8,8)→(7,8)→(7,9)→(8,9)→(9,9), with 7 material sampling points; the continuous manufacturing path of the (k-1)th layer below is (8 ...9)→(9,9), with 7 material sampling points; the continuous manufacturing path of the (k-1)th layer below is (8,8)→(8,9)→(9,9)→(9,9), with 7 material sampling points; the continuous manufacturing path of the (k-1)th layer below is (8,8)→(8,9)→(9,9)→(9,9), with 7 material sampling points; the continuous manufacturing path of the (k-1)th layer below is (8,8)→(8,9)→(9,9)→(9,9), with 7 material sampling points; the continuous manufacturing path of the (k-1)th layer below is (8,8)→(9,9)→(8,9)→(9,9), with 7 material sampling points; the continuous manufacturing path of the (k-1)th layer below is (8,8)→(9,9)→(8,9)→(9,9), with 7 material sampling points; the continuous manufacturing path of the (k-1)th layer below is (8,8)→(9,9)→(8,9)→(9,9)→(9,9), with 7 material sampling The continuous manufacturing path for layer 4 is (9,7)→(8,7)→(7,7)→(7,8)→(8,8)→(9,8)→(9,9)→(8,9)→(7,9)→(6,9)→(6,8), with 11 material sampling points. The continuous manufacturing path for layer k-5 below is (6,7)→⋯→(10,7)→(10,8)→⋯→(6,8)→(6,9)→⋯→(10,9), with 15 material sampling points. Except for the locations along these paths, the material occupancy state of all other grid locations in each layer is 0.

[0055] The feed speed of the manufacturing section is 5.00 mm / s; the feed speed of the auxiliary motion section is 100 mm / s. A filament with a diameter of 1.75 mm, a line width of 0.25 mm, and a layer height of 0.20 mm is used. This embodiment employs a relative extrusion mode. This represents the relative extrusion increment for each manufacturing segment. For a manufacturing segment with a length of 0.25 mm, the relative extrusion increment is:

[0056] Therefore, in G code, take Each layer first moves to the starting point of the continuous path using the G0 command, then converts the adjacent coordinates into the G1 forward extrusion command. The original manufacturing segment near point A is M83; relative extrusion mode is enabled. G1 X2.750 Y2.000 E0.00520 F300 G1 X2.500 Y2.000 E0.00520 G1 X2.250 Y2.000 E0.00520 G1 X2.000 Y2.000 E0.00520 ; Point A The distance between adjacent sampling points is 0.25mm, therefore the printing time interval between adjacent sampling points is: Taking the printing time of the first material sampling point in the current layer as zero, the printing time of point A is: .

[0057] The 13 sampling points located before point A that have already completed material deposition constitute the historical material sampling point set of point A, denoted as H(A), i.e.:

[0058] The sampling points (7,8) to (12,11) after point A have not yet been printed and do not belong to H(A). The sampling point in H(A) farthest from point A is (4,5), and its planar distance is... Since the radius of the local plane neighborhood is less than 1.75 mm, all 13 historical material sampling points are within the local plane neighborhood, and the material occupancy state I(Q) of each point is 1. There is no additional empty travel segment before point A, so the G0 instruction does not add time to the time difference between the 13 sampling points and point A.

[0059] In this embodiment, a preset time window is used. Take 15.0s, time decay coefficient Take 0.50. The printing time difference and time decay formulas can be calculated using formulas (2) and (3). For example, point Its printing time is 0.60s, therefore:

[0060]

[0061] In this embodiment, the axis of the local cooling main direction of the device coincides with the global X-axis. Point A is defined as the positive side along the +X direction and the negative side along the -X direction; therefore, the local direction angle θ is set to 0°. The positive side of the local main direction, i.e., the +X side, is the relatively unfavorable side for heat dissipation, and the influence scale along the major axis is taken as... The negative side of the local principal direction, i.e., the -X side, is the side with relatively favorable heat dissipation. The influence of the major axis direction on the scale is taken as follows. The minor axis direction affects the uniformity of scale. Therefore, the two-dimensional rotation matrix is: For historical material sampling points located on the positive side of the local principal direction of point A or having the same X-coordinate as point A, i.e. The kernel amplitude coefficient A is uniformly set to 0.15 for all historical material sampling points; when historical material sampling points are located on the positive and negative sides of the local principal direction, only the influence scale of the major axis direction is changed; the influence scale of the major axis on the positive side of the local principal direction is... The negative side major axis affects the scale. The minor axis affects the scale. Substituting into equation (6), we get:

[0062]

[0063] For historical material sampling points located on the negative side of the local principal direction of point A, i.e. The major axis direction affects the scale. The influence of the minor axis direction on the scale is still taken as .at this time:

[0064]

[0065] For point Its planar offset vector relative to point A is: Substituting the above plane offset vector and the corresponding covariance matrix into equation (5), we get:

[0066] Therefore, point For point A .

[0067] Calculating other sampling points in the historical material sampling point set in the same way yields:

[0068] Then, from equation (7), the planar thermal risk value of point A can be obtained:

[0069] Next, the interlayer heat transfer score at point A is calculated. A local statistical window centered at point A with a radius of 2 grids is defined as follows: Up to 10 Up to 10, there are a total of 25 grid locations. The number of material sampling points in this window for layers k-1 to k-5 below are 2, 4, 7, 11, and 15 respectively. Therefore, the local coverage rates of each lower layer are as follows:

[0070] The number of consecutive layers N participating in the statistics below is 5, and the depth attenuation coefficient is... Take 0.75. According to equation (9), the depth attenuation weights corresponding to the five layers can be calculated as follows:

[0071] Substituting into equation (8), we get:

[0072] In this embodiment, the risk gain coefficient Take 0.65. and Substituting into equation (10), we get:

[0073] Will and Substituting into equation (11), we get:

[0074] Comprehensive thermal risk value normalized upper bound Taking 1.50, therefore:

[0075] Maximum compensation strength Take 0.85, mapping slope Taking 10 as the risk threshold where compensation begins to significantly increase. Take 0.55. Substituting into the aforementioned smoothing mapping function, we obtain the thermal compensation coefficient for point A:

[0076] During the G-code refactoring phase, a one-dimensional moving average of 15 sampling points is applied to the thermal compensation coefficient sequence obtained in the current layer according to the path order. Point A is the 14th material sampling point in the current path, and the thermal compensation coefficients of the 7th to 21st material sampling points are averaged. After obtaining the thermal compensation coefficients of these 15 sampling points through the same calculation process, the smoothed thermal compensation coefficient of point A is:

[0077] Maximum reduction ratio of extrusion volume Take 0.20. Substitute the smoothed thermal compensation coefficient into equation (15):

[0078] Reference fan control quantity Take 127, the maximum fan control value. Take 255. According to the aforementioned fan control correction formula:

[0079] The control value is set to a step length of 5, rounded to the nearest multiple of 5, therefore the fan control value is quantized to 170. This value differs from the baseline fan control value of 127 by 43, reaching the minimum difference threshold of 5 for parameter updates. Therefore, a new fan control instruction is written before the manufacturing segment corresponding to point A. The corresponding corrected G-code is: G1 X2.250 Y2.000 E0.00520 M106 S170 G1 X2.000 Y2.000 E0.00487; Point A corresponds to the manufacturing section. When the local high-heat risk area corresponding to point A ends, and the target fan control quantity calculated by the subsequent manufacturing segment is restored to the reference fan control quantity 127, write instruction M106 S127 before the corresponding subsequent manufacturing segment to restore the reference fan control state.

[0080] This completes the calculation of the printing time stamp for point A, the screening of historical material sampling points, the calculation of time decay weight, the calculation of directional thermal influence kernel, the calculation of planar thermal risk value, the calculation of interlayer heat transfer score, the calculation of interlayer heat transfer correction, the calculation of comprehensive thermal risk value, the calculation and smoothing of thermal compensation coefficient, the correction of local process parameters, and the reconstruction of G-code. Similarly, for the remaining material sampling points in the current layer, based on their printing time, the set of historical material sampling points, and the material coverage within the corresponding local statistical window of the N consecutive layers below, the calculations are performed sequentially according to the above method to obtain the corrected local process parameters for the manufacturing segment corresponding to each material sampling point, and the reconstruction of G-code is completed.

[0081] Any aspects not covered in this invention are applicable to existing technologies.

Claims

1. A comprehensive thermal risk assessment method for 3D printing based on G-code, characterized in that, The evaluation method includes the following: S1: Obtain the G-code for material extrusion 3D printing; parse the motion segments in each printing layer according to the execution order of the G-code; identify the manufacturing segments where material deposition occurs and the auxiliary motion segments where material deposition does not occur; discretize the manufacturing segments to obtain material sampling points arranged according to the execution order of the G-code; and construct the manufacturing path set and material distribution map for each printing layer; determine the cumulative execution time based on the length, feed speed, and execution order of each motion segment; and set printing time markers for the material sampling points; wherein, the auxiliary motion segments do not form material occupation, but their execution time is included in the cumulative execution time; S2: Select the material sampling points sequentially as the current sampling points to be evaluated. In the current printing layer, select historical material sampling points located in the local planar neighborhood of the current sampling point to be evaluated, whose printing timestamps are earlier than the current sampling point to be evaluated, and whose printing time difference is within a preset time window. Determine the time decay weight based on the printing time difference. Determine the local orientation reference of the current sampling point to be evaluated. Determine the directional thermal influence weight based on the planar distance, relative orientation, and local orientation reference between the historical material sampling points and the current sampling point to be evaluated. Then, perform a weighted summation based on the material occupancy state, time decay weight, and directional thermal influence weight of the historical material sampling points to obtain the planar thermal risk value of the current sampling point to be evaluated. S3: Simultaneously, based on the occupied status of the printed material at the current sampling point to be evaluated, determine the interlayer heat transfer status of the area and determine the interlayer heat transfer score; use the interlayer heat transfer score to correct the planar thermal risk value and obtain the comprehensive thermal risk value.

2. The evaluation method according to claim 1, characterized in that, The process of establishing the material distribution information for each layer is as follows: The G-code is parsed line by line, and the currently valid coordinate mode, extrusion mode, feed speed, position coordinates and extrusion amount are updated according to the actual execution order of the G-code. Each motion segment is determined based on the position coordinate changes between adjacent motion commands, and the starting coordinates, ending coordinates, corresponding printing layer, feed speed, extrusion amount changes and execution order of each motion segment are obtained. Based on the planar displacement of each motion segment and the change in extrusion amount determined according to the current extrusion mode, the motion segment with planar displacement and positive change in extrusion amount is identified as the manufacturing segment, and the motion segment that does not form actual material deposition is identified as the auxiliary motion segment; the auxiliary motion segment includes at least one of the following: idle stroke segment, retraction segment, lifting segment, and transposition segment; Discrete sampling is performed along the manufacturing section at a preset step size to obtain material sampling points arranged in the order of G-code execution. The preset step size is determined based on the target linewidth and a predetermined calculation resolution, and is not greater than the target linewidth. The manufacturing segment or the material sampling point is projected onto the two-dimensional computing grid corresponding to the current printing layer. Based on the distance between the grid position in the two-dimensional computing grid and the manufacturing segment and the target line width, the material occupancy status of each grid position is determined, and a material distribution map of the current printing layer is constructed. Based on the motion length, feed speed, execution sequence, and control delay of each manufacturing segment and auxiliary motion segment, the execution start time, execution end time, and duration of each motion segment are determined. Then, based on the relative position of the material sampling point in its respective manufacturing segment, interpolation is performed between the execution start time and execution end time of the manufacturing segment to obtain the printing time mark of each material sampling point.

3. The evaluation method according to claim 2, characterized in that, The auxiliary motion segment does not form material sampling points or material occupancy states; the duration of the auxiliary motion segment is included in the cumulative execution time of the G-code and the execution start time of the subsequent manufacturing segment, and is used as the cooling time between adjacent manufacturing segments to participate in the calculation of the printing time difference and time decay weight between the historical material sampling point and the current sampling point to be evaluated, so that the auxiliary motion segment only affects the local cooling process and thermal history interval, and does not contribute to the thermal input of material deposition.

4. The evaluation method according to claim 1, characterized in that: Let the printing timestamps of the current sampling point P to be evaluated and the historical material sampling point Q be respectively... and The time difference between their printing is: (2) The time decay weight of historical material sampling point Q on the current sampling point P to be evaluated for: (3) in, For the preset time window, This is the time decay coefficient; The planar offset vector of the historical material sampling point Q relative to the current sampling point P to be evaluated is: ， The directional thermal influence weight uses a directional thermal influence kernel. And expressed as: ， ， Where A is the kernel amplitude coefficient, Let covariance matrix be the variance matrix. It is a two-dimensional rotation matrix. The local orientation reference is the local orientation angle relative to the preset coordinate axis. The scale is influenced by the major axis direction along the local directional reference. The minor axis direction perpendicular to the local direction reference influences the scale; The planar thermal risk value of the current sampling point to be evaluated for: ， in, It is located within the local planar neighborhood of the current sampling point P to be evaluated and satisfies A collection of historical material sampling points, The material occupancy state of historical material sampling point Q; The method for determining the local direction reference is as follows: when the main direction of local cooling of the equipment can be obtained, the main direction of local cooling of the equipment is used as the local direction reference; when the main direction of local cooling of the equipment cannot be obtained, the local path tangential direction of the manufacturing section where the current sampling point to be evaluated is located is used as the local direction reference.

5. The evaluation method according to claim 1, characterized in that, The specific process of step S3 is as follows: Calculate the local material coverage of the current sampling point to be evaluated in the continuous multiple layers below the current printing layer, determine the depth attenuation weight based on the distance between each lower printing layer and the current printing layer, and determine the interlayer heat transfer score based on the local material coverage and the depth attenuation weight. The planar thermal risk value is corrected based on the interlayer heat transfer score to obtain a comprehensive thermal risk value.

6. The evaluation method according to claim 5, characterized in that, The process for determining the interlayer heat transfer score is as follows: Suppose the current sampling point P to be evaluated is located in the k-th printing layer. Taking the current sampling point to be evaluated as the center, a local statistical region is established in the kn-th printing layer below the current printing layer. The ratio of the number of grids with a material occupancy state of 1 in the local statistical region to the total number of grids in the local statistical region is determined as the local material coverage rate corresponding to the kn-th printing layer. ; The interlayer heat transfer score for: ， Where N is the number of consecutive printed layers below the current printed layer that are included in the statistics. The depth decay weight is the weight corresponding to the nth lower printed layer, and the depth decay weight decreases as the distance between the lower printed layer and the current printed layer increases; The depth attenuation weight is adopted as follows: ， in, This is the depth attenuation coefficient, and ; When the number of existing printed layers below the current printed layer is less than N, only the existing printed layers below are counted, and normalized according to the depth attenuation weights of the layers actually participating in the count; no interlayer heat transfer correction is performed for the first printed layer.

7. The evaluation method according to claim 6, characterized in that: Let the risk gain coefficient be The interlayer heat transfer correction amount is determined based on the interlayer heat transfer score. for: ， The comprehensive thermal risk value of the current sampling point to be evaluated for: 。 8. A local thermal compensation method for 3D printing based on G-code, characterized in that, The comprehensive thermal risk value of the current sampling point to be evaluated is obtained by using the comprehensive thermal risk assessment method according to any one of claims 1 to 7. The comprehensive thermal risk value is normalized, and a thermal compensation coefficient is obtained based on this for local thermal compensation.

9. The compensation method according to claim 8, characterized in that, The comprehensive thermal risk value is normalized to obtain the normalized thermal risk index. : ， in, This serves as the normalized upper bound for the comprehensive thermal risk value. The normalized thermal risk index is input into a monotonically continuous smooth mapping function to obtain the thermal compensation coefficient corresponding to the current sampling point to be evaluated. : ， in, The maximum compensation strength; The mapping slope parameter; To compensate for the risk threshold that has begun to increase significantly; The material sampling points are arranged according to the G-code execution order, and the thermal compensation coefficient sequence is processed by a one-dimensional moving average along the G-code execution order to obtain the smoothed thermal compensation coefficient. Based on the smoothed thermal compensation coefficient, local thermal compensation is performed on at least one of the extrusion amount and fan control amount of the manufacturing segment corresponding to the current sampling point to be evaluated.

10. The 3D printing local thermal compensation method based on G-code according to claim 9, characterized in that, The local thermal compensation also includes G-code reconstruction, which includes: The smoothed thermal compensation coefficients are mapped to the manufacturing segments in the original G-code corresponding to each material sampling point; When there are differences in the thermal compensation coefficients at different locations within the same manufacturing segment, the manufacturing segment is divided into multiple sub-segments based on the location, magnitude, or range of change of the thermal compensation coefficients, and the local process parameters of the corresponding sub-segments are determined using the representative values ​​of the thermal compensation coefficients within each sub-segment. When mapping the smoothed thermal compensation coefficient to the extrusion rate, let the original extrusion rate increment be... The maximum extrusion reduction ratio is The corrected increment of extrusion amount is: ， When mapping the smoothed thermal compensation coefficient to the fan control quantity, let the reference fan control quantity be... The maximum fan control value is The corrected fan control variable is: ， in, The thermal compensation coefficient after smoothing; While maintaining the geometry and execution order of the original printing path, insert or rewrite at least one of the extrusion quantity control instructions and fan control instructions before or within the corresponding manufacturing segment; Set a minimum difference threshold for parameter updates. Only when the difference between the new local process parameter and the current local process parameter reaches or exceeds the minimum difference threshold for parameter updates will a new control instruction be written. Otherwise, the current local process parameter will be used. No material deposition-type thermal compensation correction is applied to the auxiliary motion segment, but the execution time of the auxiliary motion segment is retained; when the local high-heat risk area ends and the local process parameters corresponding to the subsequent manufacturing segment are restored to the reference value, a control instruction to restore the reference local process parameters is written before the subsequent manufacturing segment; The G-code is reconstructed based on the corrected local process parameters and corresponding control instructions, and the corrected G-code is output.

Images