A ceramic 3D printing sintering reverse compensation method based on measured radial shrinkage field and SOVS constitutive constraint

Abstract

The present application relates to a kind of ceramic 3D printing sintering reverse compensation methods based on measured radial shrinkage field and SOVS constitutive constraint, which includes the following steps: obtaining target three-dimensional model, preparing un-compensated ceramic sample and carrying out defatting-sintering, model space registration before and after sintering, constructing measured radial shrinkage field, generating initial reverse compensation model, establishing finite element model under the SOVS sintering constitutive constraint, inverting SOVS material parameters on measured data, compensation model sintering simulation verification, iterative update reverse compensation model, output final compensation printing model, after compensation model is determined, further actual printing compensation piece and sintering according to same process and re-measurement, to measured geometric error as the final verification basis of compensation effect.The present application solves the problems existing in the present ceramic 3D printing sintering compensation, such as too simple shrinkage rate assumption, insufficient local topography compensation accuracy, finite element simulation and measured topography disconnection, lack of physical constraint of compensation parameters and the like.

Description

Technical Field

[0001] This invention relates to the field of ceramic additive manufacturing technology, and in particular to a method for reverse compensation of ceramic 3D printing sintering based on measured radial shrinkage field and SOVS constitutive constraints. It is applicable to debinding-sintering shrinkage compensation of ceramic materials such as zirconia ceramics, alumina ceramics, silicon carbide ceramics, and silicon nitride ceramics after photopolymerization molding, slurry direct writing molding, or other ceramic additive manufacturing processes. Background Technology

[0002] Ceramic materials possess high hardness, high wear resistance, high corrosion resistance, high temperature stability, and good chemical inertness, making them widely used in aerospace, biomedicine, electronic devices, energy equipment, and civil engineering particle material simulation. With the development of additive manufacturing technology, photopolymer ceramic 3D printing, slurry direct writing, selective laser sintering, and other ceramic forming methods can prepare complex structures that are difficult to obtain using traditional molding, slip casting, or machining methods.

[0003] However, ceramic 3D printed parts typically require debinding and high-temperature sintering processes. The printed green body often contains resin, dispersant, photosensitive monomers, or other organic components. During the debinding stage, these organic components are removed. During the sintering stage, neck growth, pore closure, grain growth, and densification occur between particles, ultimately resulting in a high-density ceramic component. These processes inevitably cause significant volume shrinkage and morphological changes.

[0004] Existing ceramic 3D printing sintering compensation methods mainly adopt the following categories: The first type of method is empirical scaling, which involves scaling up the original 3D model during the modeling stage based on the overall linear shrinkage rate obtained from historical experiments. For example, if the overall linear shrinkage rate after sintering is 20%, the model is scaled up uniformly by a fixed ratio. This method is simple, but it implicitly assumes that the shrinkage is uniform, isotropic, and independent of local geometry. However, in actual ceramic components, the wall thickness, curvature, porosity, printing direction, photopolymerization depth, local constraint state, temperature gradient, and grain growth degree may all vary in different regions, resulting in significant spatial non-uniformity of sintering shrinkage. For irregularly shaped particles, complex porous structures, thin-walled structures, structures with sharp corners, and simulated particle models used in soil mechanics, overall scaling up is difficult to guarantee the accuracy of the contour after sintering.

[0005] The second type of method is unidirectional or tridirectional shrinkage rate compensation, which involves measuring the average shrinkage rate in the X, Y, and Z directions respectively and scaling it according to the direction. This method takes into account a certain degree of anisotropy, but it still cannot describe the shrinkage differences between different local areas in the same direction, and it cannot handle local morphological deviations caused by local curvature, sharp corners, local wall thickness, temperature field differences, and support contact.

[0006] The third type of method is finite element forward simulation, which establishes a ceramic sintering shrinkage model to predict the temperature, stress, and deformation fields during the sintering process. This type of method can reveal the sintering stress and shrinkage mechanism, but most applications remain at the level of forward analysis of sintering stress or deformation trends. It has not yet established a closed loop with the CT-measured three-dimensional morphology before and after sintering, nor has it used the measured local shrinkage field back to compensate the STL model.

[0007] The fourth type of method is data-driven or machine learning compensation, which involves training a model using historical samples to predict compensation coefficients. However, simple machine learning models are easily limited by the number of samples, the model's generalization ability, and the coverage of training data. When the ceramic slurry, sintering regime, model scale, or printing parameters are changed, the prediction results are prone to inaccuracy, and there is a lack of constraints from the sintering physical mechanism.

[0008] Based on the above problems, it is necessary to provide a ceramic 3D printing sintering compensation method that can utilize measured three-dimensional morphology data before and after sintering, combined with the SOVS sintering constitutive model, to achieve local radial reverse compensation and iterative convergence.

[0009] The existing technology has the following drawbacks: (1) Existing proportional compensation or triaxial compensation cannot accurately describe the local non-uniform shrinkage behavior of ceramic sintering.

[0010] (2) Existing finite element sintering simulations are mostly used for forward prediction of sintering stress or overall deformation, lacking parameter inversion and compensation optimization mechanisms coupled with CT measured local shrinkage fields.

[0011] (3) Existing compensation methods are difficult to unify microscopic material parameters, sintering temperature history, relative density evolution, grain size evolution and macroscopic STL model geometric compensation.

[0012] (4) Existing machine learning compensation methods are prone to becoming black box predictions, lack SOVS constitutive constraints, and are difficult to guarantee reliability on new material systems and new geometric models.

[0013] (5) For ceramic 3D printing models that need to maintain particle morphology characteristics, such as soil mechanics particle test, triaxial compression test, and single particle crushing test, the existing overall compensation method is difficult to guarantee the accuracy of particle edge angle, sphericity, local concave and convex features and real morphology after sintering. Summary of the Invention

[0014] This invention provides a reverse compensation method for ceramic 3D printing sintering based on measured radial shrinkage field and SOV constitutive constraints, aiming to solve the problems existing in ceramic 3D printing sintering compensation, such as overly simplistic shrinkage rate assumptions, insufficient accuracy of local morphology compensation, disconnect between finite element simulation and measured morphology, and lack of physical constraints on compensation parameters.

[0015] This invention provides a method for reverse compensation in ceramic 3D printing sintering based on measured radial shrinkage field and SOVS constitutive constraints. The method includes the following steps: Step S10: Obtain the target three-dimensional model of the ceramic component to be prepared, and perform format unification, coordinate normalization and basic network check on the target three-dimensional model. Step S20: Print the uncompensated ceramic green body directly according to the target three-dimensional model, or initially enlarge the target three-dimensional model according to the preliminary empirical shrinkage rate, and print the ceramic green body according to the enlarged model; Step S30: The ceramic green body is processed according to a preset degreasing-sintering regime to obtain a sintered ceramic component; wherein, the preset temperature-time curve is used as the process input condition for subsequent finite element sintering simulation. Step S40: Perform CT scans or three-dimensional morphology scans on the ceramic green body before sintering and the ceramic component after sintering, respectively, and perform segmentation, noise reduction, hole repair, surface reconstruction and mesh processing on the scan data to obtain the model before sintering and the model after sintering. Step S50: Spatial registration is performed on the pre-sintering model, post-sintering model, and target model to establish the point correspondence between the three. The registration methods include centroid alignment, principal inertial axis alignment, feature point registration, iterative nearest point registration, or a combination thereof. After registration, the three models are resampled to the same common reference surface, and the correspondence is established using nearest triangle projection and centroid coordinate interpolation. Abnormal points whose corresponding distance or normal angle exceeds a preset threshold are removed or rematched. Step S60: Based on the point correspondence after registration, calculate the radial change of the corresponding position before and after sintering relative to the model reference point, obtain the measured radial shrinkage field characterizing the local non-uniform shrinkage behavior of the ceramic component, and perform outlier removal and neighborhood smoothing on the measured radial shrinkage field. Step S70: Construct a vertex-level or region-level reverse compensation coefficient field based on the measured radial contraction field, perform non-uniform radial compensation on the target model, generate an initial compensation model, and perform mesh repair, patch quality check, self-intersection detection, normal unification and print adaptation processing on the initial compensation model. Step S80: Import the initial compensation model into Abaqus to establish a ceramic sintering finite element model; the ceramic sintering finite element model includes a transient heat conduction model and a sintering shrinkage mechanics model. The sintering shrinkage mechanics model is constrained by the SOVS sintering constitutive model and is used to describe the coupling relationship between temperature history, relative density evolution, grain size evolution, viscosity change, sintering stress and volume shrinkage. Step S90: Using the measured radial shrinkage field, relative density evolution data, sintered geometric error and grain size test data, the material parameters in the SOVS sintering constitutive model are inverted so that the finite element simulation results obtained by using the inverted material parameters match the measured results. Step S100: Import the initial compensation model or the updated compensation model into the ceramic sintering finite element model, use the SOVS material parameters obtained by inversion to perform sintering simulation, obtain the simulated sintered model, and calculate the geometric error between the simulated sintered model and the target model. Step S110: Based on the maximum error and root mean square error between the simulated sintered model and the target model, determine whether the current compensation model meets the preset convergence condition; if the preset convergence condition is met, proceed to step S120; if the preset convergence condition is not met, proceed to step S130. Step S120: Output the final compensation model that satisfies the preset convergence condition, and perform patch repair, normal unification, slice adaptation, and printing path check on the final compensation model to obtain the final STL model for ceramic 3D printing; print the compensation ceramic green body according to the final compensation model, and prepare the compensation ceramic component using the same debinding-sintering process as in the parameter calibration stage; perform CT scanning or three-dimensional morphology retesting on the compensation ceramic component, and use the measured maximum error and root mean square error between the retested model and the target model as the verification results of the final compensation effect; Step S130: Based on the geometric error field between the simulated sintered model and the target model, update the reverse compensation coefficient field and the compensation model, and re-check the mesh quality and print adaptation of the updated compensation model. Then return to step S100 to re-perform the finite element sintering simulation and geometric error evaluation until the preset convergence condition is met or the preset maximum number of iterations is reached.

[0016] A further technical solution of the present invention is that step S10 includes: A zirconia ceramic particle model used for soil mechanics particle mechanics tests was selected as the target model. The model has an irregular surface and a preset edge angle; wherein the edge angle is 0.00, 0.05, 0.10, 0.15 or 0.20; the target model is input in STL format and imported into MATLAB or Python programs for face inspection, deletion of duplicate vertices, repair of non-manifold edges, unification of normals, and recording of vertex set and triangle face set.

[0017] A further technical solution of the present invention is that the degreasing-sintering stage in step S30 includes a heating and holding stage from 25°C to 200°C, from 200°C to 300°C, from 300°C to 550°C, and from 550°C to 1050°C; the sintering stage includes a heating stage from 1050°C to 1480°C, a holding stage at 1480°C for 120 minutes, and a cooling stage to 150°C.

[0018] A further technical solution of the present invention is that step S40 includes: CT scans or 3D morphology scans were performed on the printed and dried ceramic green bodies, but before degreasing, to obtain the 3D morphology data of the ceramic green bodies before sintering, which was recorded as the pre-sintering model. After sintering, the sintered ceramic component is subjected to another CT scan or 3D morphology scan to obtain the post-sintering model. When using CT scanning, the scanned images undergo thresholding, noise filtering, pore repair, surface reconstruction, and meshing to obtain STL models of the ceramic components before and after sintering. If the CT reconstruction results contain internal pore information, the porosity distribution is further extracted for subsequent SOVS material parameter inversion and mapped to the initial relative density field ρ at the finite element integration points. r0 (x); When using three-dimensional topography scanning, point cloud or surface mesh data is obtained, and through denoising, point cloud registration, surface reconstruction and mesh repair, a pre-sintering model and a post-sintering model consistent with the target model coordinate system are obtained.

[0019] A further technical solution of the present invention is that step S50 includes: Model before sintering Sintered model and target model When importing into a MATLAB or Python program, first perform centroid translation on each model to make their centroids coincide. For a discrete 3D mesh model, its centroid is represented as:

[0020] in, For the model centroid, For the number of vertices in the model, For the first The coordinates of each vertex are determined. After centroid alignment, the model is initially rotated and registered using the principal axis of inertia method to ensure consistency in the main geometric directions. For irregular ceramic particle models, the iterative nearest point algorithm is further used for fine registration to establish the initial correspondence between the pre-sintering model, the post-sintering model, and the target model. Subsequently, the three models are resampled to the same common reference surface, and a one-to-one correspondence is established using nearest triangle projection and centroid coordinate interpolation, while outliers with distances or normal angles exceeding a preset threshold are removed. After registration, the corresponding points in the pre-sintering model are denoted as:

[0021] The corresponding points in the sintered model are denoted as:

[0022] The corresponding point in the target model is denoted as:

[0023] in, .

[0024] A further technical solution of the present invention is that step S60 includes: Using the model's centroid, geometric center, printing coordinate origin, or a preset reference point as the radial reference point, denoted as:

[0025] For the The radial distance between the corresponding points before sintering is:

[0026] The radial distance after sintering is:

[0027] No. The radial shrinkage rate at each point is defined as:

[0028] in, For the first Measured radial shrinkage rate at each point. The unit radial direction vector of each point is:

[0029] The corresponding radial contraction displacement is expressed as:

[0030] Radial shrinkage rate of all surface points or radial contraction displacement Together, they constitute the measured radial shrinkage field. To reduce the influence of scanning noise and local registration error, the measured radial shrinkage field is smoothed by neighborhood weighted average. The radial shrinkage field is applicable to components that are approximately star-shaped relative to the reference point and whose sintering deformation is mainly radial. When there is significant tangential displacement, warping, or non-star-shaped region, this radial compensation method is not directly used. Instead, three-dimensional displacement field compensation or local reference region is divided for the component.

[0031]

[0032] in, The smoothed radial shrinkage rate. For the first The set of neighborhood points of each vertex. , where is the neighborhood weight coefficient.

[0033] A further technical solution of the present invention is that step S70 includes: Based on the smoothed radial contraction field, calculate the first... The reverse compensation coefficients for each vertex of the target model:

[0034] Where, β i The compensation coefficient at the i-th vertex is calculated precisely based on the measured shrinkage rate. To avoid excessive iterative compensation leading to local geometric distortion, subsequent error updates use a relaxation coefficient μ to control the step size, with a value range of: 0 < μ ≤ 1 Based on the compensation coefficient, the vertices of the target model are locally radially magnified to obtain the compensated vertex coordinates:

[0035] in, For the compensation model, the first The coordinates of each vertex; or, expressed in the form of radial displacement compensation:

[0036] in:

[0037] The initial inverse compensation model is formed by the topological relationships between all compensated vertices and the original triangular facets. ; After generating the initial compensation model, the compensation model is subjected to quality checks on triangular facets, hole repair, self-intersection detection, non-manifold edge repair, normal consistency, mesh smoothing, and facet simplification to ensure that the compensation model meets the requirements for ceramic 3D printing slicing and molding.

[0038] A further technical solution of the present invention is that step S80 includes: Initial compensation model Import the data into Abaqus and establish a ceramic sintering finite element model. This model includes a transient heat conduction model and a sintering shrinkage mechanics model. For transient heat conduction analysis, a three-dimensional heat conduction element is used, with material density, thermal conductivity, and specific heat capacity input. The debinding-sintering process curve is defined as the furnace temperature amplitude curve or temperature boundary condition. The transient heat conduction control equation is: when a significant porosity gradient is detected by the CT (Cellular Temperature Coefficient), ... Assign corresponding integration points; when the printing direction causes significant anisotropy, define direction-dependent viscosity through local material orientation. Without setting the above spatial field, this model is only applicable to homogeneous, approximately isotropic billets.

[0039]

[0040] in, For material density, For specific heat capacity, For temperature, For time, Thermal conductivity, The heat source term is used; the temperature field of the ceramic component during sintering is obtained through transient heat conduction analysis, and this temperature field is used as the predefined temperature field for subsequent sintering shrinkage mechanics analysis; in the sintering shrinkage mechanics model, C3D8R elements or equivalent three-dimensional solid elements are used to define the elastic modulus, Poisson's ratio, coefficient of thermal expansion, and user CREEP subroutine; the total strain rate during sintering is expressed as:

[0041] in, The total strain rate, e For elastic strain rate, th For thermal strain rate, vp The viscoplastic or creep strain rate under deviatoric stress. s The strain rate is the volumetric shrinkage strain rate caused by sintering densification; the deviation viscoplastic strain rate is taken as... The volume shrinkage strain rate tensor is taken ,in And the compression is positive. The thermal strain rate is:

[0042] in, The coefficient of thermal expansion is For the rate of temperature change, It is a second-order unit tensor; The relationship between porosity and relative density satisfies:

[0043] in, Porosity The relative density is given. The relative density and the sintering volume shrinkage strain satisfy the mass conservation relationship:

[0044] in, The initial relative density, For sintering volume shrinkage strain; this application uniformly adopts the strain convention of tensile strain being positive and compressive strain being negative, and adopts... Furthermore, the hydrostatic pressure is defined as positive during compression; the viscosity during sintering is expressed in Arrhenius form.

[0045] in, For viscosity, The coefficient before exponential growth of the viscosity term. The activation energy for viscous flow. The gas constant is... The absolute temperature is given; the shear viscosity and bulk viscosity are expressed as follows:

[0046]

[0047] in, Shear viscosity, This is the volume viscosity; in this embodiment, we take... , Scope of application is And in numerical calculations, let To avoid singular bulk viscosity; sintering stress is expressed as:

[0048] in, For sintering stress, The specific surface energy of the material, The average grain size is used in this embodiment. The grain size evolution is expressed by the grain growth kinetic equation as follows:

[0049] in, This is the current grain size. This is the initial grain size. The grain growth index. This is the pre-exponential coefficient for grain growth. The activation energy for grain growth is given by: The volumetric shrinkage strain during sintering densification is expressed as:

[0050] Where p is the hydrostatic pressure according to the positive convention of compression. For volume viscosity; increment at each time step Within, the increase in strain due to sintering volume shrinkage is:

[0051] In the Abaqus subroutine, a complete sintering constitutive equation is established using Fortran. The established viscoelastic constitutive equation is defined in the CREEP subroutine of Abaqus to form a complete SOVS sintering constitutive model. The CREEP subroutine reads temperature, time, stress, state variables, and material parameters at each integration point and at each time increment, calculates relative density, porosity, grain size, viscosity, shear viscosity, bulk viscosity, sintering stress, and sinter shrinkage strain increment, and returns the updated state variables to the Abaqus solver. Specifically, the deviation viscoplastic strain increment is written into DECRA(1) instead of directly into the deviation strain tensor. According to the deviation viscoplastic strain rate relationship, let

[0052] in, This is the Mises equivalent stress, corresponding to QTILD in the CREEP subroutine. Shear viscosity, The time increment is used; Abaqus converts the equivalent deviation viscoplastic strain increment into a deviation viscoplastic strain tensor based on the Mises flow direction. DESWA(1) writes the volumetric sintering shrinkage strain increment, i.e.

[0053] in, To compress the equivalent pressure stress to positive, corresponding to P in the CREEP subroutine, For sintering stress, The value is the bulk viscosity; the negative sign indicates that sintering produces volume shrinkage under the sign convention of positive tensile strain. The material definition simultaneously enables CREEP, LAW=USER and SWELLING, LAW=USER, and updates STATEV according to the aforementioned state evolution equation in each time increment. When using explicit creep integral, i.e., LEXIMP=0, only DECRA(1) and DESWA(1) need to be defined, and a sufficiently small time increment is used to control the calculation accuracy and stability; when using implicit creep integral, i.e., LEXIMP=1, consistent derivatives corresponding to DECRA(2)~DECRA(5) and DESWA(2)~DESWA(5) should also be provided according to the Abaqus interface requirements to ensure iterative convergence. Reducing the time increment cannot replace the derivatives required for implicit integral.

[0054] The CREEP subroutine shall at least set the following state variables:

[0055] in, Linear shrinkage rate, Let be the volumetric shrinkage rate. Under the assumption of isotropic shrinkage, , .

[0056] A further technical solution of the present invention is that step S90 includes: Based on experimental radial shrinkage field, relative density evolution curves, post-sintering geometric errors, and grain size test results, SOVS material parameters are inverted; the set of parameters to be inverted is denoted as: θ = [η0, Q] v D, Q g ] Where η0 is the pre-exponential coefficient of the viscosity term, Q v Q is the activation energy for viscous flow, D is the pre-exponential coefficient for grain growth, and Q is the activation energy for viscous flow. g The activation energy for grain growth; the initial relative density ρ r0 The initial grain size g0 is directly determined experimentally and used as a known initial condition. The specific surface energy γ is preferentially determined by independent experiments or reliable literature, and γ is included in the inversion parameter set only when independent constraint data are available. The objective function for material parameter inversion is expressed as:

[0057] in, Let be the overall objective function. This is the relative density evolution error term. For the radial contraction field error term, This is the geometric morphology error term after sintering. This is the grain size error term. , , , These are the weighting coefficients for each error term; the relative density evolution error term is defined as follows: Each error term is dimensionless according to the measurement variance, tolerance error, or characteristic scale, and , The radial error term and the geometric shape error term are corrected for correlation or weighted by one method to avoid repeated weighting of the same surface data.

[0058]

[0059] in, To use parameters The simulation results of the first Relative density at a given moment For the first The experimental calibration of relative density at each time point, The relative density sampling point number is given. The radial contraction field error term is defined as:

[0060] in, The simulation results are as follows: Radial contraction displacement of a surface point The first obtained from CT scan or three-dimensional topographic scan Measured radial contraction displacement at one surface point Let be the feature length of the target model. The geometric shape error term after sintering is defined as:

[0061] in, The simulated sintering result is the first... Coordinates of a surface point For the actual measurement of the sintering stage Coordinates of surface points; the grain size error term is defined as:

[0062] in, The simulated grain size after sintering The sintered grain size is obtained through scanning electron microscopy or other microscopic characterization methods; the material parameter inversion problem is expressed as:

[0063] in, The optimal SOVS parameters obtained through inversion. The parameter range is defined. During the inversion process, genetic algorithms, particle swarm optimization, Bayesian optimization, Levenberg-Marquardt algorithm, gradient descent algorithm, or a combination of multiple algorithms are used to minimize the objective function. When the initial experimental data is insufficient, the relative density evolution error term and the radial shrinkage field error term are used for preliminary inversion. After obtaining the CT model and grain size test results after sintering, the geometric morphology error term and the grain size error term after sintering are introduced for secondary correction.

[0064] A further technical solution of the present invention is that steps S100-S130 include: SOVS material parameters obtained by inversion For the initial compensation model Perform sintering simulation to obtain the simulated sintered model. The simulated sintered model is spatially registered with the target model, and errors are calculated. The geometric error of each vertex is defined as:

[0065] in, To compensate for the model after sintering simulation, the first vertex coordinates The coordinates of the corresponding vertex in the target model; the maximum error is:

[0066] The root mean square error is:

[0067] The current compensation model is output as the final compensation model only when both the maximum error and the root mean square error meet the preset accuracy threshold; the successful convergence condition is:

[0068]

[0069]

[0070] in, The maximum error threshold, Let δ be the root mean square error threshold, δ be the error change rate threshold between two adjacent iterations, and k be the iteration number. If the error change rate between two adjacent iterations is not greater than δ but the accuracy threshold is not met, it is only considered as iteration stall, and the calculation is continued after reducing the relaxation coefficient μ, re-inverting the parameters, or adjusting the mesh; it is not considered as successful convergence. If the error does not meet the accuracy condition, the compensation model is updated according to the radial error between the simulated sintered model and the target model. In the k-th iteration, the compensation radius is updated as follows:

[0071] in, For the first During the nth iteration The compensation radius of each vertex, For the target model The radius of each vertex For the first The model after the second simulation sintering The radius of each vertex The relaxation coefficient is used; the updated compensation vertex coordinates are:

[0072] Alternatively, the error displacement field can be used for updating:

[0073] After updating, mesh quality checks, finite element sintering simulations, and error assessments are performed again until the aforementioned accuracy conditions are met. If the accuracy conditions are still not met after reaching the maximum number of iterations, a non-convergence flag is output, and the parameter inversion or compensation step size adjustment steps are returned. This model is not considered a final qualified compensation model. The compensation model that meets the accuracy conditions is actually printed, degreased and sintered according to the same procedure, and retested to verify the compensation effect with measured geometric errors.

[0074] 1. This invention uses a measured radial shrinkage field instead of a single overall shrinkage rate, which can describe the local non-uniform shrinkage during the ceramic sintering process and improve the dimensional accuracy of ceramic components with complex morphology after sintering.

[0075] 2. This invention incorporates real geometric data obtained from CT scans or three-dimensional morphology scans before and after sintering into the compensation model, so that the compensation coefficients are derived from measured morphological changes, rather than relying solely on empirical shrinkage rates.

[0076] 3. This invention introduces the SOVS sintering constitutive model into the compensation iteration, and improves the physical interpretability of the compensation method by predicting shrinkage constrained by relative density, porosity, grain size, viscosity and sintering stress.

[0077] 4. This invention uses material parameter inversion to adapt the finite element model to current ceramic slurry, printing equipment, debinding process, and sintering process, avoiding the increase in error caused by directly applying empirical compensation coefficients between different material systems.

[0078] 5. This invention can realize a closed-loop iteration of "actual measurement - inversion - simulation - compensation - verification", which can continuously reduce the geometric deviation between the simulation sintering model and the target model.

[0079] 6. This invention can perform local reverse compensation on STL vertices, and is applicable to particle models and freeform ceramic components that are approximately star-shaped and dominated by radial contraction.

[0080] 7. This invention can preserve the angularity, surface unevenness and irregular morphology of ceramic particle models, and is suitable for triaxial compression tests in soil mechanics, single particle crushing tests and mechanical behavior studies of particle aggregates.

[0081] 8. The machine learning module of this invention is an optional initial prediction branch, which can reduce the number of iterations by utilizing historical data and ensure reliability through SOVS physical constraints. Attached Figure Description

[0082] Figure 1 This is the overall flowchart of the ceramic 3D printing sintering reverse compensation method based on the measured radial shrinkage field and SOVS constitutive constraints of the present invention.

[0083] Figure 2 A schematic diagram of spatial registration of the pre-sintering model, the post-sintering model, and the ideal model.

[0084] Figure 3 This is a radial error cloud diagram between the pre-sintering model, the post-sintering model, and the ideal model.

[0085] Figure 4 This is a comparison chart of the STL model before and after reverse compensation.

[0086] Figure 5 This is a contour map of the displacement field for reverse compensation.

[0087] Figure 6 The flowchart shows the finite element sintering simulation and iterative convergence process for the compensation model.

[0088] Figure 7 This is a radial error cloud diagram between the simulated sintered model and the ideal model. Detailed Implementation

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

[0090] This invention proposes a reverse compensation method for ceramic 3D printing sintering based on measured radial shrinkage field and SOVS constitutive constraints. The aim is to establish an experimentally driven, physically constrained, geometrically realizable and iteratively optimized ceramic particle compensation method. Under the conditions of fixed material system, fixed printing process and fixed sintering curve, the radial error between the sintered particle model and the ideal model is reduced round by round, thereby improving the accuracy of ceramic particle morphology control.

[0091] The technical solution adopted in this invention mainly involves obtaining real three-dimensional models of ceramic components before and after sintering through CT scanning or three-dimensional morphology scanning. After registration, the local radial shrinkage field is calculated and used as the basis for SOVS sintering constitutive parameter inversion and STL model reverse compensation. Through finite element simulation and iterative optimization, the compensated model can more closely approximate the target model after undergoing the same debinding-sintering process, thereby improving the dimensional accuracy, morphology retention, and experimental repeatability of the ceramic 3D printed parts.

[0092] Specifically, such as Figures 1 to 7 As shown, a preferred embodiment of the ceramic 3D printing sintering reverse compensation method based on measured radial shrinkage field and SOVS constitutive constraint of the present invention includes the following steps: S10. Obtain the target 3D model: Obtain the target 3D model of the ceramic component to be prepared, and perform format standardization, coordinate normalization, and basic mesh check on the model. The target model can be an STL, OBJ, PLY, or STEP format model, or a 3D mesh model obtained by CT reconstruction.

[0093] Step S20: Print the uncompensated ceramic green body directly according to the target three-dimensional model, or initially enlarge the target three-dimensional model according to the preliminary empirical shrinkage rate, and print the ceramic green body according to the enlarged model; Step S30: The ceramic green body is processed according to a preset degreasing-sintering regime to obtain a sintered ceramic component; wherein, the preset temperature-time curve is used as the process input condition for subsequent finite element sintering simulation. Step S40: Perform CT scans or three-dimensional morphology scans on the ceramic green body before sintering and the ceramic component after sintering, respectively, and perform segmentation, noise reduction, hole repair, surface reconstruction and mesh processing on the scan data to obtain the model before sintering and the model after sintering. Step S50: Spatial registration is performed on the pre-sintering model, post-sintering model, and target model to establish the point correspondence between the three. The registration methods include centroid alignment, principal inertial axis alignment, feature point registration, iterative nearest point registration, or a combination thereof. After registration, the three models are resampled to the same common reference surface, and the correspondence is established using nearest triangle projection and centroid coordinate interpolation. Abnormal points whose corresponding distance or normal angle exceeds a preset threshold are removed or rematched. Step S60: Based on the point correspondence after registration, calculate the radial change of the corresponding position before and after sintering relative to the model reference point, obtain the measured radial shrinkage field characterizing the local non-uniform shrinkage behavior of the ceramic component, and perform outlier removal and neighborhood smoothing on the measured radial shrinkage field. Step S70: Construct a vertex-level or region-level reverse compensation coefficient field based on the measured radial contraction field, perform non-uniform radial compensation on the target model, generate an initial compensation model, and perform mesh repair, patch quality check, self-intersection detection, normal unification and print adaptation processing on the initial compensation model. Step S80: Import the initial compensation model into Abaqus to establish a ceramic sintering finite element model; the ceramic sintering finite element model includes a transient heat conduction model and a sintering shrinkage mechanics model. The sintering shrinkage mechanics model is constrained by the SOVS sintering constitutive model and is used to describe the coupling relationship between temperature history, relative density evolution, grain size evolution, viscosity change, sintering stress and volume shrinkage. Step S90: Using the measured radial shrinkage field, relative density evolution data, sintered geometric error and grain size test data, the material parameters in the SOVS sintering constitutive model are inverted so that the finite element simulation results obtained by using the inverted material parameters match the measured results. Step S100: Import the initial compensation model or the updated compensation model into the ceramic sintering finite element model, use the SOVS material parameters obtained by inversion to perform sintering simulation, obtain the simulated sintered model, and calculate the geometric error between the simulated sintered model and the target model. Step S110: Based on the maximum error and root mean square error between the simulated sintered model and the target model, determine whether the current compensation model meets the preset convergence condition; if the preset convergence condition is met, proceed to step S120; if the preset convergence condition is not met, proceed to step S130. Step S120: Output the final compensation model that satisfies the preset convergence condition, and perform patch repair, normal unification, slice adaptation, and printing path check on the final compensation model to obtain the final STL model for ceramic 3D printing; print the compensation ceramic green body according to the final compensation model, and prepare the compensation ceramic component using the same debinding-sintering process as in the parameter calibration stage; perform CT scanning or three-dimensional morphology retesting on the compensation ceramic component, and use the measured maximum error and root mean square error between the retested model and the target model as the verification results of the final compensation effect; Step S130: Based on the geometric error field between the simulated sintered model and the target model, update the reverse compensation coefficient field and the compensation model, and re-check the mesh quality and print adaptation of the updated compensation model. Then return to step S100 to re-perform the finite element sintering simulation and geometric error evaluation until the preset convergence condition is met or the preset maximum number of iterations is reached.

[0094] In this embodiment, after compensation, the STL model undergoes triangular facet quality checks, hole repair, self-intersection detection, facet simplification, and mesh smoothing to ensure that the compensated model meets the manufacturability requirements for ceramic 3D printing.

[0095] This embodiment forms a closed-loop iteration of "actual measurement - inversion - simulation - compensation - re-simulation - verification". When the geometric error between the simulated sintered model and the target model meets the preset convergence condition, the final compensation model is output. "Verification" includes actual printing of the final compensation part, sintering with the same process, and CT or 3D morphology remeasurement. The verification sample or verification geometry does not participate in the aforementioned parameter inversion.

[0096] After compensation, the STL model undergoes triangular facet quality checks, hole repair, self-intersection detection, facet simplification, and mesh smoothing to ensure that the compensated model meets the manufacturability requirements for ceramic 3D printing.

[0097] The beneficial effects of the reverse compensation method for ceramic 3D printing sintering based on measured radial contraction field and SOVS constitutive constraint in this embodiment are: 1. This invention uses a measured radial shrinkage field instead of a single overall shrinkage rate, which can describe the local non-uniform shrinkage during the ceramic sintering process and improve the dimensional accuracy of ceramic components with complex morphology after sintering.

[0098] 2. This invention incorporates real geometric data obtained from CT scans or three-dimensional morphology scans before and after sintering into the compensation model, so that the compensation coefficients are derived from measured morphological changes, rather than relying solely on empirical shrinkage rates.

[0099] 3. This invention introduces the SOVS sintering constitutive model into the compensation iteration, and improves the physical interpretability of the compensation method by predicting shrinkage constrained by relative density, porosity, grain size, viscosity and sintering stress.

[0100] 4. This invention uses material parameter inversion to adapt the finite element model to current ceramic slurry, printing equipment, debinding process, and sintering process, avoiding the increase in error caused by directly applying empirical compensation coefficients between different material systems.

[0101] 5. This invention can realize a closed-loop iteration of "actual measurement - inversion - simulation - compensation - verification", which can continuously reduce the geometric deviation between the simulation sintering model and the target model.

[0102] 6. This invention can perform local reverse compensation on STL vertices, and is applicable to particle models and freeform ceramic components that are approximately star-shaped and dominated by radial contraction.

[0103] 7. This invention can preserve the angularity, surface unevenness and irregular morphology of ceramic particle models, and is suitable for triaxial compression tests in soil mechanics, single particle crushing tests and mechanical behavior studies of particle aggregates.

[0104] 8. The machine learning module of this invention is an optional initial prediction branch and is not a core necessary step. It can reduce the number of iterations by utilizing historical data and ensure reliability through SOVS physical constraints.

[0105] Furthermore, in this embodiment, step S10 specifically includes: A zirconia ceramic particle model used for soil mechanics particle mechanics tests is selected as the target model. The model has an irregular surface and a preset edge angle; wherein, the edge angle can be 0.00, 0.05, 0.10, 0.15 or 0.20; the target model is input in STL format, imported into MATLAB or Python programs for face inspection, deletion of duplicate vertices, repair of non-manifold edges, unification of normals, and recording of vertex set and triangle face set.

[0106] This embodiment uses zirconia ceramic material. Typical material parameters are shown in Table 1.

[0107] Table 1

[0108] Among the parameters mentioned above, the initial average grain size g0 and the initial relative density ρ r0 Directly calibrated by experiments; initial billet density according to Calculation, where This is the theoretical density of the actual material, determined through true density testing or supplier data. If the material systems are different, the experimentally measured value of the corresponding material or the inverse value constrained by independent data should be used.

[0109] Furthermore, in this embodiment, the degreasing-sintering stage in step S30 includes heating and holding stages from 25°C to 200°C, from 200°C to 300°C, from 300°C to 550°C, and from 550°C to 1050°C; the sintering stage includes heating from 1050°C to 1480°C, holding at 1480°C for 120 minutes, and cooling to 150°C.

[0110] In this embodiment, a continuous degreasing-sintering process is adopted, and the specific temperature-time-relative density relationship is shown in Table 2.

[0111] Table 2

[0112] Furthermore, in this embodiment, step S40 specifically includes: Model acquisition and 3D reconstruction before and after sintering: CT scans or 3D morphology scans are performed on the printed and dried ceramic green bodies, but before degreasing, to obtain the 3D morphology data of the ceramic green bodies before sintering, which is denoted as the pre-sintering model. After sintering, the sintered ceramic component is subjected to another CT scan or 3D morphology scan to obtain the post-sintering model. When using CT scanning, the scanned images undergo thresholding, noise filtering, pore repair, surface reconstruction, and meshing to obtain STL models of the ceramic components before and after sintering. If the CT reconstruction results contain internal pore information, the porosity distribution is further extracted for subsequent SOVS material parameter inversion and mapped to the initial relative density field ρ at the finite element integration points. r0 (x); When using three-dimensional topography scanning, point cloud or surface mesh data is obtained, and through denoising, point cloud registration, surface reconstruction and mesh repair, a pre-sintering model and a post-sintering model consistent with the target model coordinate system are obtained.

[0113] Furthermore, in this embodiment, step S50 specifically includes: Spatial registration of the pre-sintering model, post-sintering model, and target model: The pre-sintering model... Sintered model and target model When importing into a MATLAB or Python program, first perform centroid translation on each model to make their centroids coincide. For a discrete 3D mesh model, its centroid is represented as:

[0114] in, For the model centroid, For the number of vertices in the model, For the first The coordinates of each vertex are determined. After centroid alignment, the model is initially rotated and registered using the principal axis of inertia method to ensure consistency in the main geometric directions. For irregular ceramic particle models, the iterative nearest point algorithm is further used for fine registration to establish the initial correspondence between the pre-sintering model, the post-sintering model, and the target model. Subsequently, the three models are resampled to the same common reference surface, and a one-to-one correspondence is established using nearest triangle projection and centroid coordinate interpolation, while outliers with distances or normal angles exceeding a preset threshold are removed. After registration, the corresponding points in the pre-sintering model are denoted as:

[0115] The corresponding points in the sintered model are denoted as:

[0116] The corresponding point in the target model is denoted as:

[0117] in, .

[0118] Furthermore, in this embodiment, step S60 specifically includes: Calculation of measured radial contraction field: The model's centroid, geometric center, printed coordinate origin, or preset reference point are used as radial reference points, denoted as:

[0119] For the The radial distance between the corresponding points before sintering is:

[0120] The radial distance after sintering is:

[0121] No. The radial shrinkage rate at each point is defined as:

[0122] in, For the first Measured radial shrinkage rate at each point. The unit radial direction vector of each point is:

[0123] The corresponding radial contraction displacement is expressed as:

[0124] Radial shrinkage rate of all surface points or radial contraction displacement Together, they constitute the measured radial shrinkage field. To reduce the influence of scanning noise and local registration error, the measured radial shrinkage field is smoothed by neighborhood weighted average. The radial shrinkage field is applicable to components that are approximately star-shaped relative to the reference point and whose sintering deformation is mainly radial. When there is significant tangential displacement, warping, or non-star-shaped region, this radial compensation method is not directly used. Instead, three-dimensional displacement field compensation or local reference region is divided for the component.

[0125]

[0126] in, The smoothed radial shrinkage rate. For the first The set of neighborhood points of each vertex. , where is the neighborhood weight coefficient.

[0127] Furthermore, in this embodiment, step S70 specifically includes: Initial reverse compensation model generation: Based on the smoothed radial contraction rate field, calculate the first... The reverse compensation coefficients for each vertex of the target model:

[0128] Where, β i The compensation coefficient at the i-th vertex is calculated precisely based on the measured shrinkage rate. To avoid excessive iterative compensation leading to local geometric distortion, subsequent error updates use a relaxation coefficient μ to control the step size, with a value range of: 0 < μ ≤ 1 The vertex coordinates of the target model are obtained by locally radially magnifying the vertex based on the compensation coefficient:

[0129] in, For the compensation model, the first The coordinates of each vertex; or, expressed in the form of radial displacement compensation:

[0130] in:

[0131] The initial inverse compensation model is formed by the topological relationships between all compensated vertices and the original triangular facets. ; After generating the initial compensation model, the compensation model is subjected to quality checks on triangular facets, hole repair, self-intersection detection, non-manifold edge repair, normal consistency, mesh smoothing, and facet simplification to ensure that the compensation model meets the requirements for ceramic 3D printing slicing and molding.

[0132] Furthermore, in this embodiment, step S80 specifically includes: Establish an Abaqus finite element model under SOVS sintering constitutive constraints: The initial compensation model... Import the data into Abaqus and establish a ceramic sintering finite element model. This model includes a transient heat conduction model and a sintering shrinkage mechanics model. For transient heat conduction analysis, a three-dimensional heat conduction element is used, with material density, thermal conductivity, and specific heat capacity input. The debinding-sintering process curve is defined as the furnace temperature amplitude curve or temperature boundary condition. The transient heat conduction control equation is: when a significant porosity gradient is detected by CT, ρ... r0 (x) Assign the corresponding integration point; when the printing direction causes significant anisotropy, define the direction-dependent viscosity through the local material direction. Without setting the above spatial field, this model is only used for homogeneous, approximately isotropic billets.

[0133]

[0134] in, For material density, For specific heat capacity, For temperature, For time, Thermal conductivity, The heat source term is used; the temperature field of the ceramic component during sintering is obtained through transient heat conduction analysis, and this temperature field is used as the predefined temperature field for subsequent sintering shrinkage mechanics analysis; in the sintering shrinkage mechanics model, C3D8R elements or equivalent three-dimensional solid elements are used to define the elastic modulus, Poisson's ratio, coefficient of thermal expansion, and user CREEP subroutine; the total strain rate during sintering is expressed as:

[0135] in, The total strain rate, e For elastic strain rate, th For thermal strain rate, vp The viscoplastic or creep strain rate under deviatoric stress. s The strain rate is the volumetric shrinkage strain rate caused by sintering densification; the deviation viscoplastic strain rate is taken as... The volume shrinkage strain rate tensor is taken ,in And the compression is positive. The thermal strain rate is:

[0136] in, The coefficient of thermal expansion is For the rate of temperature change, It is a second-order unit tensor; The relationship between porosity and relative density satisfies:

[0137] in, Porosity The relative density is given. The relative density and the sintering volume shrinkage strain satisfy the mass conservation relationship:

[0138] in, The initial relative density, For sintering volume shrinkage strain; this application uniformly adopts the strain convention of tensile strain being positive and compressive strain being negative, and adopts... Furthermore, the hydrostatic pressure is defined as positive during compression; the viscosity during sintering is expressed in Arrhenius form.

[0139] in, For viscosity, The coefficient before exponential growth of the viscosity term. The activation energy for viscous flow. The gas constant is... The absolute temperature is given; the shear viscosity and bulk viscosity are expressed as follows:

[0140]

[0141] in, Shear viscosity, This is the volume viscosity; in this embodiment, we take... , Scope of application is And in numerical calculations, let To avoid singular bulk viscosity; sintering stress is expressed as:

[0142] in, For sintering stress, The specific surface energy of the material, The average grain size is used in this embodiment. The grain size evolution is expressed by the grain growth kinetic equation as follows:

[0143] in, This is the current grain size. This is the initial grain size. The grain growth index. This is the pre-exponential coefficient for grain growth. The activation energy for grain growth is given by: The volumetric shrinkage strain during sintering densification is expressed as:

[0144] Where p is the hydrostatic pressure according to the positive convention of compression. For volume viscosity; increment at each time step Within, the increase in strain due to sintering volume shrinkage is:

[0145] In the Abaqus subroutine, a complete sintering constitutive equation is established using Fortran. The established viscoelastic constitutive equation is defined in the CREEP subroutine of Abaqus to form a complete SOVS sintering constitutive model. The CREEP subroutine reads temperature, time, stress, state variables, and material parameters at each integration point and at each time increment, calculates relative density, porosity, grain size, viscosity, shear viscosity, bulk viscosity, sintering stress, and sinter shrinkage strain increment, and returns the updated state variables to the Abaqus solver. Specifically, the deviation viscoplastic strain increment is written into DECRA(1) instead of directly into the deviation strain tensor. According to the deviation viscoplastic strain rate relationship, let

[0146] in, This is the Mises equivalent stress, corresponding to QTILD in the CREEP subroutine. Shear viscosity, The time increment is used; Abaqus converts the equivalent deviation viscoplastic strain increment into a deviation viscoplastic strain tensor based on the Mises flow direction. DESWA(1) writes the volumetric sintering shrinkage strain increment, i.e.

[0147] in, To compress the equivalent pressure stress to positive, corresponding to P in the CREEP subroutine, For sintering stress, The value is the bulk viscosity; the negative sign indicates that sintering produces volume shrinkage under the sign convention of positive tensile strain. The material definition simultaneously enables CREEP, LAW=USER and SWELLING, LAW=USER, and updates STATEV according to the aforementioned state evolution equation in each time increment. When using explicit creep integral, i.e., LEXIMP=0, only DECRA(1) and DESWA(1) need to be defined, and a sufficiently small time increment is used to control the calculation accuracy and stability; when using implicit creep integral, i.e., LEXIMP=1, consistent derivatives corresponding to DECRA(2)~DECRA(5) and DESWA(2)~DESWA(5) should also be provided according to the Abaqus interface requirements to ensure iterative convergence. Reducing the time increment cannot replace the derivatives required for implicit integral.

[0148] The CREEP subroutine shall at least set the following state variables:

[0149] in, Linear shrinkage rate, Let be the volumetric shrinkage rate. Under the assumption of isotropic shrinkage, , .

[0150] Furthermore, in this embodiment, step S90 specifically includes: SOVS material parameter inversion: Based on experimental radial shrinkage field, relative density evolution curve, post-sintering geometric error, and grain size test results, SOVS material parameters are inverted; the set of parameters to be inverted is denoted as: θ = [η0, Q] v D, Q g ] Where η0 is the pre-exponential coefficient of the viscosity term, Q v Q is the activation energy for viscous flow, D is the pre-exponential coefficient for grain growth, and Q is the activation energy for viscous flow. g The activation energy for grain growth; the initial relative density ρ r0 The initial grain size g0 is directly determined experimentally and used as a known initial condition. The specific surface energy γ is preferentially determined by independent experiments or reliable literature, and γ is included in the inversion parameter set only when independent constraint data are available. The objective function for material parameter inversion is expressed as:

[0151] in, Let be the overall objective function. This is the relative density evolution error term. For the radial contraction field error term, This is the geometric morphology error term after sintering. This is the grain size error term. , , , These are the weighting coefficients for each error term; the relative density evolution error term is defined as follows: each error term is dimensionless according to the measurement variance, tolerance error, or characteristic scale, and , The radial error term and the geometric shape error term are corrected for correlation or weighted by one method to avoid repeated weighting of the same surface data.

[0152]

[0153] in, To use parameters The simulation results of the first Relative density at a given moment For the first The experimental calibration of relative density at each time point, The relative density sampling point number is given. The radial contraction field error term is defined as:

[0154] in, The simulation results are as follows: Radial contraction displacement of a surface point The first obtained from CT scan or three-dimensional topographic scan Measured radial contraction displacement at one surface point Let be the feature length of the target model. The geometric shape error term after sintering is defined as:

[0155] in, The simulated sintering result is the first... Coordinates of a surface point For the actual measurement of the sintering stage Coordinates of surface points; the grain size error term is defined as:

[0156] in, The simulated grain size after sintering The sintered grain size is obtained through scanning electron microscopy or other microscopic characterization methods; the material parameter inversion problem is expressed as:

[0157] in, The optimal SOVS parameters obtained through inversion. The parameter range is defined. During the inversion process, genetic algorithms, particle swarm optimization, Bayesian optimization, Levenberg-Marquardt algorithm, gradient descent algorithm, or a combination of multiple algorithms are used to minimize the objective function. When the initial experimental data is insufficient, the relative density evolution error term and the radial shrinkage field error term are used for preliminary inversion. After obtaining the CT model and grain size test results after sintering, the geometric morphology error term and the grain size error term after sintering are introduced for secondary correction.

[0158] Furthermore, in this embodiment, steps S100-S130 specifically include: Compensation model sintering simulation, error judgment and iterative update: using SOVS material parameters obtained by inversion For the initial compensation model Perform sintering simulation to obtain the simulated sintered model. The simulated sintered model is spatially registered with the target model, and errors are calculated. The geometric error of each vertex is defined as:

[0159] in, To compensate for the model after sintering simulation, the first vertex coordinates The coordinates of the corresponding vertex in the target model; the maximum error is:

[0160] The root mean square error is:

[0161] The current compensation model is output as the final compensation model only when both the maximum error and the root mean square error meet the preset accuracy threshold; the successful convergence condition is:

[0162]

[0163]

[0164] in, The maximum error threshold, Let δ be the root mean square error threshold, δ be the error change rate threshold between two adjacent iterations, and k be the iteration number. If the error change rate between two adjacent iterations is not greater than δ but the accuracy threshold is not met, it is only considered as iteration stall, and the calculation is continued after reducing the relaxation coefficient μ, re-inverting the parameters, or adjusting the mesh; it is not considered as successful convergence. If the error does not meet the accuracy condition, the compensation model is updated according to the radial error between the simulated sintered model and the target model. In the k-th iteration, the compensation radius is updated as follows:

[0165] in, For the first During the nth iteration The compensation radius of each vertex, For the target model The radius of each vertex For the first The model after the second simulation sintering The radius of each vertex The relaxation coefficient is used; the updated compensation vertex coordinates are:

[0166] Alternatively, the error displacement field can be used for updating:

[0167] After updating, mesh quality checks, finite element sintering simulations, and error assessments are performed again until the aforementioned accuracy conditions are met. If the accuracy conditions are still not met after reaching the maximum number of iterations, a non-convergence flag is output, and the parameter inversion or compensation step size adjustment steps are returned; the model is not considered a final qualified compensation model. For compensation models that meet the accuracy conditions, actual printing, debinding and sintering according to the same procedure, and retesting are performed to verify the compensation effect using measured geometric errors. The beneficial effects of the reverse compensation method for ceramic 3D printing sintering based on measured radial shrinkage field and SOVS constitutive constraint of this invention are: 1. This invention uses a measured radial shrinkage field instead of a single overall shrinkage rate, which can describe the local non-uniform shrinkage during the ceramic sintering process and improve the dimensional accuracy of ceramic components with complex morphology after sintering.

[0168] 2. This invention incorporates real geometric data obtained from CT scans or three-dimensional morphology scans before and after sintering into the compensation model, so that the compensation coefficients are derived from measured morphological changes, rather than relying solely on empirical shrinkage rates.

[0169] 3. This invention introduces the SOVS sintering constitutive model into the compensation iteration, and improves the physical interpretability of the compensation method by predicting shrinkage constrained by relative density, porosity, grain size, viscosity and sintering stress.

[0170] 4. This invention uses material parameter inversion to adapt the finite element model to current ceramic slurry, printing equipment, debinding process, and sintering process, avoiding the increase in error caused by directly applying empirical compensation coefficients between different material systems.

[0171] 5. This invention can realize a closed-loop iteration of "actual measurement - inversion - simulation - compensation - verification", which can continuously reduce the geometric deviation between the simulation sintering model and the target model.

[0172] 6. This invention can perform local reverse compensation on STL vertices, and is applicable to particle models and freeform ceramic components that are approximately star-shaped and dominated by radial contraction.

[0173] 7. This invention can preserve the angularity, surface unevenness and irregular morphology of ceramic particle models, and is suitable for triaxial compression tests in soil mechanics, single particle crushing tests and mechanical behavior studies of particle aggregates.

[0174] 8. The machine learning module of this invention is an optional initial prediction branch, which can reduce the number of iterations by utilizing historical data and ensure reliability through SOVS physical constraints.

[0175] The above description, in conjunction with specific preferred embodiments, provides a further detailed explanation of the present invention. It should not be construed that the specific implementation of the present invention is limited to these descriptions. For those skilled in the art, various simple deductions or substitutions can be made without departing from the concept of the present invention, and all such modifications and substitutions should be considered within the scope of protection of the present invention.

Claims

1. A method for reverse compensation in ceramic 3D printing sintering based on measured radial shrinkage field and SOVS constitutive constraints, characterized in that, The method includes the following steps: Step S10: Obtain the target three-dimensional model of the ceramic component to be prepared, and perform format unification, coordinate normalization and basic network check on the target three-dimensional model. Step S20: Print the uncompensated ceramic green body directly according to the target three-dimensional model, or initially enlarge the target three-dimensional model according to the preliminary empirical shrinkage rate, and print the ceramic green body according to the enlarged model; Step S30: The ceramic green body is processed according to a preset degreasing-sintering regime to obtain a sintered ceramic component; wherein, the preset temperature-time curve is used as the process input condition for subsequent finite element sintering simulation. Step S40: Perform CT scans or three-dimensional morphology scans on the ceramic green body before sintering and the ceramic component after sintering, respectively, and perform segmentation, noise reduction, hole repair, surface reconstruction and mesh processing on the scan data to obtain the model before sintering and the model after sintering. Step S50: Spatial registration is performed on the pre-sintering model, post-sintering model, and target model to establish the point correspondence between the three. The registration methods include centroid alignment, principal inertial axis alignment, feature point registration, iterative nearest point registration, or a combination thereof. After registration, the three models are resampled to the same common reference surface, and the correspondence is established using nearest triangle projection and centroid coordinate interpolation. Abnormal points whose corresponding distance or normal angle exceeds a preset threshold are removed or rematched. Step S60: Based on the point correspondence after registration, calculate the radial change of the corresponding position before and after sintering relative to the model reference point, obtain the measured radial shrinkage field characterizing the local non-uniform shrinkage behavior of the ceramic component, and perform outlier removal and neighborhood smoothing on the measured radial shrinkage field. Step S70: Construct a vertex-level or region-level reverse compensation coefficient field based on the measured radial contraction field, perform non-uniform radial compensation on the target model, generate an initial compensation model, and perform mesh repair, patch quality check, self-intersection detection, normal unification and print adaptation processing on the initial compensation model. Step S80: Import the initial compensation model into Abaqus to establish a ceramic sintering finite element model; the ceramic sintering finite element model includes a transient heat conduction model and a sintering shrinkage mechanics model. The sintering shrinkage mechanics model is constrained by the SOVS sintering constitutive model and is used to describe the coupling relationship between temperature history, relative density evolution, grain size evolution, viscosity change, sintering stress and volume shrinkage. Step S90: Using the measured radial shrinkage field, relative density evolution data, sintered geometric error and grain size test data, the material parameters in the SOVS sintering constitutive model are inverted so that the finite element simulation results obtained by using the inverted material parameters match the measured results. Step S100: Import the initial compensation model or the updated compensation model into the ceramic sintering finite element model, use the SOVS material parameters obtained by inversion to perform sintering simulation, obtain the simulated sintered model, and calculate the geometric error between the simulated sintered model and the target model. Step S110: Based on the maximum error and root mean square error between the simulated sintered model and the target model, determine whether the current compensation model meets the preset convergence condition; if the preset convergence condition is met, proceed to step S120; if the preset convergence condition is not met, proceed to step S130. Step S120: Output the final compensation model that satisfies the preset convergence condition, and perform patch repair, normal unification, slice adaptation, and printing path check on the final compensation model to obtain the final STL model for ceramic 3D printing; print the compensation ceramic green body according to the final compensation model, and prepare the compensation ceramic component using the same debinding-sintering process as in the parameter calibration stage; perform CT scanning or three-dimensional morphology retesting on the compensation ceramic component, and use the measured maximum error and root mean square error between the retested model and the target model as the verification results of the final compensation effect; Step S130: Based on the geometric error field between the simulated sintered model and the target model, update the reverse compensation coefficient field and the compensation model, and re-check the mesh quality and print adaptation of the updated compensation model. Then return to step S100 to re-perform the finite element sintering simulation and geometric error evaluation until the preset convergence condition is met or the preset maximum number of iterations is reached. Step S40 specifically includes: CT scans or 3D morphological scans were performed on the printed and dried ceramic green bodies, but before degreasing, to obtain the 3D morphological data of the ceramic green bodies before sintering, which was recorded as the pre-sintering model. After sintering, the sintered ceramic component is subjected to another CT scan or 3D morphology scan to obtain the post-sintering model. When using CT scanning, the scanned images undergo thresholding, noise filtering, pore repair, surface reconstruction, and meshing to obtain STL models of the ceramic components before and after sintering. If the CT reconstruction results contain internal pore information, the porosity distribution is further extracted for subsequent SOVS material parameter inversion and mapped to the initial relative density field of the finite element integration points. When using 3D topography scanning, point cloud or surface mesh data is obtained, and through denoising, point cloud registration, surface reconstruction and mesh repair, a pre-sintering model and a post-sintering model consistent with the coordinate system of the target model are obtained. Step S50 includes: Model before sintering Sintered model and target model When importing into a MATLAB or Python program, first perform centroid translation on each model to make their centroids coincide; for a discrete 3D mesh model, its centroid is represented as: in, For the model centroid, For the number of vertices in the model, For the first The coordinates of each vertex are determined. After centroid alignment, the model is initially rotated and registered using the principal axis of inertia method to ensure consistency in the main geometric directions. For irregular ceramic particle models, the iterative nearest point algorithm is further used for fine registration to establish the initial correspondence between the pre-sintering model, the post-sintering model, and the target model. Subsequently, the three models are resampled to the same common reference surface, and a one-to-one correspondence is established using nearest triangle projection and centroid coordinate interpolation, while outliers with distances or normal angles exceeding a preset threshold are removed. After registration, the corresponding points in the pre-sintering model are denoted as: The corresponding points in the sintered model are denoted as: The corresponding point in the target model is denoted as: in, ; Step S60 includes: Using the model's centroid, geometric center, printing coordinate origin, or a preset reference point as the radial reference point, denoted as: For the The radial distance between the corresponding points before sintering is: The radial distance after sintering is: No. The radial shrinkage rate at each point is defined as: in, For the first The measured radial shrinkage rate at point n, the th The unit radial direction vector of each point is: The corresponding radial contraction displacement is expressed as: Radial shrinkage rate of all surface points or radial contraction displacement Together, they constitute the measured radial shrinkage field. To reduce the influence of scanning noise and local registration error, the measured radial shrinkage field is smoothed by neighborhood weighted average. The radial shrinkage field is applicable to components that are approximately star-shaped relative to the reference point and whose sintering deformation is mainly radial. When there is significant tangential displacement, warping, or non-star-shaped region, this radial compensation method is not directly used. Instead, three-dimensional displacement field compensation or local reference region is divided for the component. in, The smoothed radial shrinkage rate. For the first The set of neighborhood points of each vertex. These are the neighborhood weight coefficients; Step S70 includes: Based on the smoothed radial contraction field, calculate the first... The reverse compensation coefficients for each vertex of the target model: Where, β i The compensation coefficient at the i-th vertex is calculated precisely based on the measured shrinkage rate. To avoid excessive iterative compensation leading to local geometric distortion, subsequent error updates use a relaxation coefficient μ to control the step size, with a value range of: 0 < μ ≤ 1 Based on the compensation coefficient, the vertices of the target model are locally radially magnified to obtain the compensated vertex coordinates: in, For the compensation model, the first The coordinates of each vertex; or, expressed in the form of radial displacement compensation: in: The initial inverse compensation model is composed of the topological relationships between all compensated vertices and the original triangular facets. ; After generating the initial compensation model, the compensation model is subjected to triangular facet quality checks, hole repair, self-intersection detection, non-manifold edge repair, normal consistency, mesh smoothing, and facet simplification to ensure that the compensation model meets the requirements for ceramic 3D printing slicing and forming. Steps S100-S130 include: SOVS material parameters obtained by inversion For the initial compensation model Sintering simulation was performed to obtain the simulated sintered model. The simulated sintered model is spatially registered with the target model, and errors are calculated. The geometric error of each vertex is defined as: in, To compensate for the model after sintering simulation, the first vertex coordinates The coordinates of the corresponding vertex in the target model; the maximum error is: The root mean square error is: The current compensation model is output as the final compensation model only when both the maximum error and the root mean square error meet the preset accuracy threshold; the successful convergence condition is: in, The maximum error threshold, Let δ be the root mean square error threshold, δ be the error change rate threshold between two adjacent iterations, and k be the iteration number. If the error change rate between two adjacent iterations is not greater than δ but the accuracy threshold is not met, it is only considered as iteration stall, and the calculation is continued after reducing the relaxation coefficient μ, re-inverting the parameters, or adjusting the mesh; it is not considered as successful convergence. If the error does not meet the accuracy condition, the compensation model is updated according to the radial error between the simulated sintered model and the target model. In the k-th iteration, the compensation radius is updated as follows: in, For the first During the nth iteration The compensation radius of each vertex, For the target model The radius of each vertex For the first The model after the second simulation sintering The radius of each vertex The relaxation coefficient is used; the updated compensation vertex coordinates are: Alternatively, the error displacement field can be used for updating: After updating, perform mesh quality checks, finite element sintering simulations, and error assessments again until the aforementioned accuracy conditions are met. If the accuracy conditions are still not met after reaching the maximum number of iterations, output a non-convergence flag and return the parameter inversion or compensation step size adjustment steps. Do not consider the model as the final qualified compensation model. For compensation models that meet the accuracy conditions, actually print them, degrease and sinter them according to the same procedure, and retest them to verify the compensation effect with actual measured geometric errors.

2. The method for reverse compensation of ceramic 3D printing sintering based on measured radial shrinkage field and SOVS constitutive constraint according to claim 1, characterized in that, Step S10 includes: A zirconia ceramic particle model used for soil mechanics particle mechanics tests was selected as the target model. The model has an irregular surface and a preset edge angle; wherein the edge angle is 0.00, 0.05, 0.10, 0.15 or 0.20; the target model is input in STL format and imported into MATLAB or Python programs for face inspection, deletion of duplicate vertices, repair of non-manifold edges, unification of normals, and recording of vertex set and triangle face set.

3. The method for reverse compensation of ceramic 3D printing sintering based on measured radial shrinkage field and SOVS constitutive constraint as described in claim 2, characterized in that, The degreasing-sintering stage in step S30 includes heating and holding stages from 25°C to 200°C, from 200°C to 300°C, from 300°C to 550°C, and from 550°C to 1050°C; the sintering stage includes heating from 1050°C to 1480°C, holding at 1480°C for 120 minutes, and cooling to 150°C.

4. The method for reverse compensation of ceramic 3D printing sintering based on measured radial shrinkage field and SOVS constitutive constraint according to claim 1, characterized in that, Step S80 includes: Initial compensation model Import the data into Abaqus and establish a ceramic sintering finite element model. This model includes a transient heat conduction model and a sintering shrinkage mechanics model. For transient heat conduction analysis, a three-dimensional heat conduction element is used, with material density, thermal conductivity, and specific heat capacity input. The debinding-sintering process curve is defined as the furnace temperature amplitude curve or temperature boundary condition. The transient heat conduction control equation is: when a significant porosity gradient is detected by the CT (Cellular Temperature Coefficient), ... Assign corresponding integration points; when the printing direction causes significant anisotropy, define the direction-dependent viscosity through the local material direction; when no spatial field is set, this model is only used for homogeneous, approximately isotropic billets; in, For material density, For specific heat capacity, For temperature, For time, Thermal conductivity, The heat source term is used; the temperature field of the ceramic component during sintering is obtained through transient heat conduction analysis, and this temperature field is used as the predefined temperature field for subsequent sintering shrinkage mechanics analysis; in the sintering shrinkage mechanics model, C3D8R elements or equivalent three-dimensional solid elements are used to define the elastic modulus, Poisson's ratio, coefficient of thermal expansion, and user CREEP subroutine; the total strain rate during sintering is expressed as: in, The total strain rate, e For elastic strain rate, th For thermal strain rate, vp The viscoplastic or creep strain rate under deviatoric stress. s The strain rate is the volumetric shrinkage strain rate caused by sintering densification; the deviation viscoplastic strain rate is taken as... The volume shrinkage strain rate tensor is taken ,in And the compression is positive; the thermal strain rate is: in, The coefficient of thermal expansion is... For the rate of temperature change, It is a second-order unit tensor; The relationship between porosity and relative density satisfies: in, Porosity The relative density is given; the relative density and the sintering volume shrinkage strain satisfy the mass conservation relationship: in, The initial relative density, The strain represents the volumetric shrinkage strain during sintering; where tension is positive and compression is negative, and the strain convention is adopted. Furthermore, the hydrostatic pressure is defined as positive during compression; the viscosity during sintering is expressed in Arrhenius form. in, For viscosity, The coefficient before exponential growth of the viscosity term. The activation energy for viscous flow. The gas constant is Temperature; shear viscosity and bulk viscosity are expressed as follows: in, Shear viscosity, Let be the bulk viscosity; take , Scope of application is And in numerical calculations, let To avoid singular bulk viscosity; sintering stress is expressed as: in, For sintering stress, The specific surface energy of the material, The average grain size is taken as... Grain size evolution is expressed by the grain growth kinetic equation as follows: in, This is the current grain size. This is the initial grain size. The grain growth index. This is the pre-exponential coefficient for grain growth. The activation energy for grain growth is given by: The volumetric shrinkage strain during sintering densification is expressed as: Where p is the hydrostatic pressure according to the positive convention of compression. For volume viscosity; increment at each time step Within, the increase in strain due to sintering volume shrinkage is: In the Abaqus subroutine, a complete sintering constitutive equation is established using Fortran. The established viscoelastic constitutive equation is defined in the CREEP subroutine of Abaqus to form a complete SOVS sintering constitutive model. The CREEP subroutine reads temperature, time, stress, state variables, and material parameters at each integration point and at each time increment, calculates relative density, porosity, grain size, viscosity, shear viscosity, bulk viscosity, sintering stress, and sinter shrinkage strain increment, and returns the updated state variables to the Abaqus solver. Specifically, the deviation viscoplastic strain increment is written into DECRA(1) instead of directly into the deviation strain tensor. According to the deviation viscoplastic strain rate relationship, let in, This is the Mises equivalent stress, corresponding to QTILD in the CREEP subroutine. Shear viscosity, The time increment; Abaqus converts the equivalent deviation viscoplastic strain increment into a deviation viscoplastic strain tensor based on the Mises flow direction; DESWA(1) writes the volumetric sintering shrinkage strain increment, i.e. in, To compress the equivalent pressure stress to positive, corresponding to P in the CREEP subroutine, For sintering stress, For volumetric viscosity; the negative sign indicates that sintering produces volume shrinkage under the sign convention of positive tensile strain; the material definition simultaneously enables CREEP, LAW=USER and SWELLING, LAW=USER, and updates STATEV according to the state evolution equation in each time increment; when using explicit creep integral, i.e. LEXIMP=0, only DECRA(1) and DESWA(1) need to be defined, and a sufficiently small time increment is used to control the calculation accuracy and stability; when using implicit creep integral, i.e. LEXIMP=1, consistent derivatives corresponding to DECRA(2)~DECRA(5) and DESWA(2)~DESWA(5) should also be provided according to the Abaqus interface requirements to ensure iterative convergence, and reducing the time increment cannot replace the derivative required by implicit integral; The CREEP subroutine shall at least set the following state variables: in, Linear shrinkage rate, The volumetric shrinkage rate; under the assumption of isotropic shrinkage, , .

5. The method for reverse compensation of ceramic 3D printing sintering based on measured radial shrinkage field and SOVS constitutive constraint according to claim 4, characterized in that, Step S90 includes: Based on experimental radial shrinkage field, relative density evolution curves, post-sintering geometric errors, and grain size test results, SOVS material parameters are inverted; the set of parameters to be inverted is denoted as: θ = [η0,Q v ,D,Q g ] Where η0 is the pre-exponential coefficient of the viscosity term, Q v Q is the activation energy for viscous flow, D is the pre-exponential coefficient for grain growth, and Q is the activation energy for viscous flow. g The activation energy for grain growth; the initial relative density ρ r0 The initial grain size g0 is directly determined experimentally and used as a known initial condition. The specific surface energy γ is preferentially determined by independent experiments or reliable literature, and γ is included in the inversion parameter set only when independent constraint data are available. The objective function for material parameter inversion is expressed as: in, Let be the overall objective function. This is the relative density evolution error term. For the radial contraction field error term, This is the geometric morphology error term after sintering. This is the grain size error term. , , , These are the weighting coefficients for each error term; the relative density evolution error term is defined as follows: each error term is dimensionless according to the measurement variance, tolerance error, or characteristic scale, and , The radial error term and the geometric shape error term are corrected for correlation or weighted by one method to avoid repeated weighting of the same surface data. in, To use parameters The simulation results of the first Relative density at a given moment For the first The experimental calibration of relative density at each time point, The relative density sampling point number; the radial contraction field error term is defined as: in, The simulation results are as follows: Radial contraction displacement of a surface point The first obtained from CT scan or three-dimensional topographic scan Measured radial contraction displacement at one surface point The feature length of the target model; the geometric shape error term after sintering is defined as: in, The simulated sintering result is the first... Coordinates of a surface point For the actual measurement of the sintering stage Coordinates of surface points; the grain size error term is defined as: in, The simulated grain size after sintering The sintered grain size is obtained through scanning electron microscopy or other microscopic characterization methods; the material parameter inversion problem is expressed as: in, The optimal SOVS parameters obtained through inversion. The parameter range is defined. During the inversion process, genetic algorithms, particle swarm optimization, Bayesian optimization, Levenberg-Marquardt algorithm, gradient descent algorithm, or a combination of multiple algorithms are used to minimize the objective function. When the initial experimental data is insufficient, the relative density evolution error term and the radial shrinkage field error term are used for preliminary inversion. After obtaining the CT model and grain size test results after sintering, the geometric morphology error term and the grain size error term after sintering are introduced for secondary correction.

Images