CN116680962A - Method for predicting surface quality and residual stress of formed part under selective laser melting technology - Google Patents

Abstract

The invention relates to the technical field of metal additive manufacturing, in particular to a method for predicting the surface quality and residual stress of a formed part under a selective laser melting technology, which is used for solving the problem that the surface quality and the residual stress of the formed part under the selective laser melting technology cannot be predicted simultaneously in the prior art, and comprises the following steps: acquiring an initial powder particle distribution model in a selective laser melting process; dividing a computing grid, dispersing the whole material domain into object points and initializing the computing grid; searching a background grid in the initialized calculation grid; making a classification marking strategy, and classifying and marking the background grid unit and the substance points by combining a distribution model; calculating and updating particle temperature of the object particles based on the result of the classification mark, determining flow behavior of the fluid in the fluid region under the drive of heat, and solving residual stress and solid deformation of the solid region; the surface quality of the shaped part in the selective laser melting technology is predicted according to particle temperature, flow behavior and solid deformation.

Description

Method for predicting surface quality and residual stress of formed part under selective laser melting technology

Technical Field

The invention relates to the field of metal additive manufacturing, in particular to a method for predicting surface quality and residual stress of a formed part under a selective laser melting technology.

Background

In recent years, additive manufacturing technology has attracted widespread attention in the domestic and foreign industries. The additive manufacturing technology is a revolutionary manufacturing technology which adopts high-energy heat beams to melt materials layer by layer and accumulate the materials layer by layer and can realize near net shape forming of any complex shape. Unlike traditional turning and other material reducing and equal material producing technology, the material adding process obtains complicated part configuration in layer-by-layer smelting and accumulating mode and has the advantages of short period, low cost, saving in material, etc. Among them, selective laser melting technology is one of the common metal additive manufacturing technologies. However, the parts prepared by the Selective Laser Melting (SLM) technology at present have the problem of difficult control of thermal stress and dimensional deviation, and are particularly characterized by rough and uneven surfaces of the parts and warp deformation caused by thermal stress. The existing prediction method for the surface quality and the residual stress of the molded part by the selective laser melting technology can only predict the surface quality or the residual stress singly, can not predict the surface quality or the residual stress simultaneously, and lacks consideration of the interaction between the surface quality and the residual stress.

Disclosure of Invention

Aiming at the defects in the prior art, the invention provides a method for predicting the surface quality and residual stress of a formed part under a selective laser melting technology.

In order to achieve the above purpose, the method for predicting the surface quality and the residual stress of a formed part under the selective laser melting technology provided by the invention comprises the following steps: modeling a powder bed laying process by using a discrete element method to obtain an initial powder particle distribution model in a selective laser melting process; dividing a calculation grid by using preset calculation parameters, dispersing the whole material domain into object points, and initializing the calculation grid by using the object points; searching a background grid in the initialized calculation grid according to the quality field of the background grid node in the object point method; formulating a classification marking strategy, and classifying and marking the background grid cells and the substance points in the background grid by combining the distribution model; calculating and updating particle temperature of the object particles based on the result of the classification mark, determining flow behavior of fluid in the fluid region under the drive of heat, and solving residual stress and solid deformation of the solid region; the surface quality of the shaped part under the selective laser melting technology is predicted according to the particle temperature, the flow behavior and the solid deformation. The invention can simultaneously predict the surface quality and residual stress of the formed part under the selective laser melting technology.

Optionally, the classification marking strategy includes the steps of:

classifying and marking the background grid cells by utilizing physical information carried by neighbor grid cells of the background grid cells;

and classifying and marking the material points according to the physical information carried by the material points and the distribution model.

Furthermore, classifying and marking the background grid units and the substance points is convenient for respectively obtaining the temperature change of the metal material, the flow behavior of the fluid region, the solid deformation and the residual stress of the solid region, so as to further realize the prediction of the surface quality and the residual stress of the formed part.

Optionally, the classifying and marking the background grid unit by using the physical information carried by the background grid unit and the neighboring grid units of the background grid unit includes the following steps:

marking a background grid cell which does not contain object points as an empty cell, otherwise marking the background grid cell as a metal material cell;

marking the metal material units of which the neighbor grid units are the empty units as surface grid units;

marking a metal material unit with a temperature less than the solidus temperature as a solid unit;

The metallic material units having a temperature between the solidus temperature and the liquidus temperature are labeled as paste-like zone units.

Further, the metal material unit participates in calculation of particle temperature to obtain temperature change of the metal material; the surface grid unit and the pasty area unit participate in the calculation of particle velocity and particle position so as to obtain the flowing behavior of the fluid under the drive of heat; the solid units participate in the calculation of residual stress and solid deformation.

Optionally, the classifying and marking the material points according to the physical information carried by the material points and the distribution model includes the following steps:

marking material points in the powder particles as unmelted material points and marking the rest material points as solid material points according to the distribution model;

marking the unmelted points of material having a temperature greater than the liquidus temperature as fluid particles, and marking the unmelted points of material having a temperature between the solidus temperature and the liquidus temperature as paste-like region particles;

the material points with the temperature exceeding the solidus temperature are named as intermediate material points, the solid stress carried by the intermediate material points with the temperature not lower than the liquidus temperature is set to zero and marked as fluid material points, the solid stress carried by the intermediate material points with the temperature between the solidus temperature and the liquidus temperature is set to zero and marked as pasty area material points, and the intermediate material points with the temperature lower than the solidus temperature are marked as solid material points.

Optionally, the metal material units other than the solid unit are non-solid units, the non-solid units include the pasty area unit, the area formed by the non-solid units is the fluid area, and the area formed by the solid units is the solid area;

the method for calculating and updating the particle temperature of the object particles based on the result of the classification mark, determining the flow behavior of the fluid in the fluid region under the drive of heat, and solving the residual stress and the solid deformation of the solid region comprises the following steps:

determining a heat transfer solution domain within the background grid based on the locations of the object points, and calculating and updating a particle temperature for each object point in the background grid cell within the heat transfer solution domain;

solving the particle velocity and the particle position of the fluid particles and the paste region particles, and determining the flow behavior of the fluid in the fluid region under the heat driving according to the particle velocity and the particle position;

and solving residual stress at the solid particles and total point force of background grid nodes in the solid region by adopting a standard mass point method, so as to solve solid deformation of the solid region.

Optionally, the determining a heat transfer solution domain in the background grid based on the positions of the object particles, and calculating and updating the particle temperature of each object particle in the background grid unit in the heat transfer solution domain includes the steps of:

determining whether a new metal material unit appears in a calculation domain, and initializing the temperature of the new metal material unit;

determining the heat transfer solving domain by using the positions of object points in all the metal material units, circulating all the metal material units in the heat transfer solving domain, solving a transient heat conduction equation of the metal material units by adopting a finite difference method to obtain a discrete format equation of the transient heat conduction equation, and obtaining a predicted temperature increment of the metal material units by utilizing the discrete format equation, wherein the transient heat conduction equation and the discrete format equation respectively satisfy the following relations:

，

，

wherein ,for the material density->For the specific heat capacity of the material, D is differentiated, T is the temperature, T is the time variable, ++>Is a metal material volume fraction->Is metal material density->Is latent heat of phase change->For fluid volume fraction, ++>For gradient operator- >Is heat conduction coefficient>Vaporization heat loss caused by vaporization of metal material, +.>For heat loss caused by convection heat exchange->For radiating heat loss>Is a volumetric heat source->Temperature probe solution for the next temperature of the metallic material element, < >>For the predicted value of the temperature at the next moment of the metallic material unit,/->For the current temperature of the metallic material element, < >>For the pre-treatment ofTemperature increment, I/O>For the predicted fluid volume fraction of the next moment of the metallic material unit, +.>A current fluid volume fraction for the metallic material unit;

adding the predicted temperature increment and the current temperature of the metal material unit to obtain the predicted value;

if the predicted value does not span the temperature interval between the solidus temperature and the liquidus temperature, directly updating the lattice core temperature of the metal material unit by using the predicted value;

if the predicted value spans the temperature interval, iteratively updating the lattice core temperature by using a semi-implicit iteration method;

and directly giving the cell core temperature to the object points in the metal material unit, and further updating the particle temperature of each object point in the background grid unit.

Furthermore, the particle temperature can reflect the temperature of the metal material in the processing area in real time, so that the temperature change engineering of the metal material in the processing area can be obtained in real time by calculating the particle temperature in real time, and the surface quality of a formed part can be predicted.

Optionally, the solving the particle velocity and the particle position of the fluid particles and the paste-like region particles and determining the flow behavior of the fluid in the fluid region under thermal drive according to the particle velocity and the particle position comprises the steps of:

dividing each surface grid cell into a plurality of sub-cells, creating a virtual influence domain of object points in the surface grid cells, determining the occupancy rate of the object points in the surface grid cells to the sub-cells through the virtual influence domain, and taking the occupancy rate as the grid core volume fraction of the surface grid cells;

solving the surface force on the surface grid unit by using the grid core volume fraction and adopting a staggered derivative method, and simultaneously calculating the recoil pressure existing on the surface grid unit where vaporization and evaporation occur;

calculating and updating particle velocities of the fluid particles and the paste-like region particles based on the surface force and the recoil pressure, and further calculating and updating particle positions of the fluid particles and the paste-like region particles using the particle velocities;

determining a flow behavior of the fluid in the fluid region under thermal drive based on the particle velocity and the particle location.

Optionally, the calculating the recoil pressure existing on the surface grid unit where vaporization and evaporation occur by using the grid core volume fraction and adopting an interlaced derivative method to solve the surface force on the surface grid unit includes the following steps:

calculating a grid normal vector of the surface grid unit at the grid center according to the grid center volume fraction;

interpolating the grid core volume fraction to each node through a shape function N to obtain a node volume fraction;

calculating the gradient of the node volume fraction, and further obtaining a node normal vector;

calculating the lattice center curvature of the surface grid unit at the lattice center by using the node normal vector and the gradient of the shape function N;

calculating the surface force according to the lattice center normal vector and the lattice center curvature, and simultaneously calculating the recoil pressure by using the lattice center normal vector, wherein the surface force and the recoil pressure respectively meet the following relations:

，

，

wherein ,for said surface force of the surface grid cell J +.>For the tension coefficient>For the curvature of the lattice center of the surface lattice unit J, < >>D represents the differential operation for the lattice normal vector of the surface lattice unit J, T is the temperature, < >>For gradient operator->For fluid density- >Is metal material density->For the density of the ambient gas>Is a metal material volume fraction->For the recoil pressure of the surface grid cell J, +.>For recoil pressure coefficient, +.>For reference pressure +.>For evaporating latent heat->Is the molar mass of the metallic material, +.>For evaporating temperature, ++>Is an ideal gas constant.

Optionally, the calculating and updating the particle velocities of the fluid particles and the paste-like region particles based on the surface force and the recoil pressure, and further calculating and updating the particle positions of the fluid particles and the paste-like region particles using the particle velocities includes the steps of:

splitting an N-S equation containing a Darcy damping term by introducing an intermediate velocity field to obtain an intermediate velocity field of a node in the non-solid unit;

solving the lattice core pressure of the non-solid unit at the lattice core by utilizing the intermediate speed field, and interpolating the lattice core pressure through a linear function to obtain the pressure gradient at each node in the non-solid unit;

updating the junction speeds of the junctions in the non-solid unit using the intermediate velocity field and the pressure gradient;

and calculating the acceleration of each node in the non-solid unit by using the surface force, the recoil pressure and the node speed, and calculating and updating the particle speeds of the fluid particles and the paste-like region particles according to the acceleration, so as to update the particle positions of the fluid particles and the paste-like region particles by using the particle speeds.

Optionally, calculating the strain change rate at the solid finger points by means of the velocity gradient of the solid particles by using a standard mass point method, and further calculating the residual stress at the solid particles by means of an elastoplastic constitutive equation;

solving the total point force of the background grid nodes in the solid area by adopting a discrete solving mode of a momentum equation in a standard mass point method;

and updating the speed and the position of the solid particles by adopting a standard mass point method according to the total point force, so as to obtain the solid deformation of the solid region.

Further, the residual stress at the solid particles is the predicted residual stress of the formed part.

In summary, the temperature change of the metal material unit is obtained by calculating the particle temperature of each object particle in the background grid unit, the flow behavior of the fluid in the fluid area under the heat driving is determined by calculating the particle speeds and particle positions of the fluid particle and the paste area particle, the solid deformation is obtained by calculating the total dot force of the background grid node in the solid area, meanwhile, the residual stress of the formed piece is directly calculated, and the surface quality of the formed piece is predicted by utilizing the temperature change of the metal material unit, the flow behavior of the fluid in the fluid area under the heat driving and the solid deformation of the solid area, so that the surface quality and the residual stress of the formed piece under the selective laser melting technology are predicted simultaneously.

Drawings

In order to more clearly illustrate the technical solutions of the embodiments of the present application, the drawings that are needed in the embodiments will be briefly described below, it being understood that the following drawings only illustrate some embodiments of the present application and therefore should not be considered as limiting the scope, and other related drawings may be obtained according to these drawings without inventive effort for a person skilled in the art.

FIG. 1 is a schematic flow chart of a method for predicting surface quality and residual stress of a formed part by a selective laser melting technology according to an embodiment of the application;

FIG. 2 is a schematic diagram of virtual domains of object points according to an embodiment of the present application;

FIG. 3 is a schematic diagram of spatial information for solving gradients in accordance with an embodiment of the present application.

Detailed Description

Specific embodiments of the application will be described in detail below, it being noted that the embodiments described herein are for illustration only and are not intended to limit the application. In the following description, numerous specific details are set forth in order to provide a thorough understanding of the present application. However, it will be apparent to one of ordinary skill in the art that: no such specific details are necessary to practice the application. In other instances, well-known circuits, software, or methods have not been described in detail in order not to obscure the application.

Throughout the specification, references to "one embodiment," "an embodiment," "one example," or "an example" mean: a particular feature, structure, or characteristic described in connection with the embodiment or example is included within at least one embodiment of the invention. Thus, the appearances of the phrases "in one embodiment," "in an embodiment," "one example," or "an example" in various places throughout this specification are not necessarily all referring to the same embodiment or example. Furthermore, the particular features, structures, or characteristics may be combined in any suitable combination and/or sub-combination in one or more embodiments or examples. Moreover, those of ordinary skill in the art will appreciate that the illustrations provided herein are for illustrative purposes and that the illustrations are not necessarily drawn to scale.

It should be noted in advance that in an alternative embodiment, the meaning and number of the same symbol or alphabet appearing in all formulas are the same, except where separate descriptions are made; in addition, in this embodiment, "current" is "nth time" and "next" is "n+1th" time.

In an alternative embodiment, referring to fig. 1, the present invention provides a method for predicting surface quality and residual stress of a formed part by selective laser melting, the method comprising the steps of:

S1, modeling a powder bed laying process by using a discrete element method to obtain an initial powder particle distribution model in a selective laser melting process.

Specifically, in this embodiment, this step is to obtain the initial distribution and morphology of powder particles in the selective laser melting process, that is, the initial distribution model of powder particles, so as to perform classification marking on the material points later.

Further, this step is prior art and is not described in detail herein.

S2, dividing a calculation grid for the calculation domain by using preset calculation parameters, dispersing the whole material domain into object points, and initializing the calculation grid by using the object points.

Specifically, in this embodiment, the calculation grid is divided according to preset calculation parameters, and the whole material domain is discretized into object particles, and then the physical quantity carried by the grid unit and the physical quantity carried by the particles are initialized. Reference is made specifically to the standard particle method, which is not described in detail herein.

Further, the preset calculation parameters comprise the length of the grid unit and the number of initial substance points in the unit.

S3, searching the background grid in the initialized calculation grid according to the quality field of the background grid node in the object point method.

And S4, formulating a classification marking strategy, and classifying and marking the background grid units and the substance points in the background grid by combining the distribution model.

Specifically, in this embodiment, the classification marking policy specifically includes the following steps:

and F1, classifying and marking the background grid cells by utilizing physical information carried by neighbor grid cells of the background grid cells.

Wherein, F1 comprises the following steps:

f11 marks the background grid cells that do not contain the object particles as empty cells, otherwise as metallic material cells.

F12 marks the metallic material cells of which the neighbor grid cells are the empty cells as surface grid cells.

F13 marks the metallic material unit having a temperature less than the solidus temperature as a solid unit.

F14 marks the metallic material unit having a temperature between the solidus temperature and the liquidus temperature as a paste-like region unit.

And F2, classifying and marking the material points according to the physical information carried by the material points and the distribution model.

Wherein, F2 comprises the following steps:

f21, marking the unmelted material points with the temperature being higher than the liquidus temperature as fluid particles, and marking the unmelted material points with the temperature being between the solidus temperature and the liquidus temperature as pasty area particles.

F22, designating a material point with a temperature exceeding the solidus temperature as an intermediate material point, setting solid stress carried by the intermediate material point with a temperature not less than the liquidus temperature as a fluid material point, setting solid stress carried by the intermediate material point with a temperature between the solidus temperature and the liquidus temperature as a pasty region material point, and marking the intermediate material point with a temperature less than the solidus temperature as a solid material point.

Further, for convenience of the following description, this embodiment further makes the following description: the other metal material units except the solid unit are non-solid units, the non-solid units comprise pasty area units, the area formed by the non-solid units is a fluid area, and the area formed by the solid units is a solid area.

And S5, calculating and updating the particle temperature of the object particles based on the result of the classification mark, determining the flow behavior of the fluid in the fluid region under the drive of heat, and solving the residual stress and the solid deformation of the solid region.

Wherein, S5 further comprises the following steps:

s51, a heat transfer solution domain is determined based on the positions of the object particles in the background grid, and the particle temperature of each object particle in the background grid unit is calculated and updated in the heat transfer solution domain.

Wherein, S51 further comprises the following steps:

s511, determining whether a new metal material unit appears in the calculation domain, and initializing the temperature of the new metal material unit.

Specifically, in this embodiment, when the metal material is processed by using the selective laser melting technology, since the fluid in the molten pool has fluidity, the mass point in the metal material unit will also move, when the mass point moves to the empty unit to change the empty unit into the metal material unit, a new metal material unit appears, and at this time, the temperature of the newly appearing metal material unit needs to be initialized, so that the newly appearing metal material unit can participate in the process of calculating and updating the mass point temperature, and then the mass point temperature of each mass point in the background grid unit can be accurately calculated and updated.

S512, determining the heat transfer solving domain by using the positions of the material points in all the metal material units, circulating all the metal material units in the heat transfer solving domain, solving a transient heat conduction equation of the metal material unit by adopting a finite difference method to obtain a discrete format equation of the transient heat conduction equation, and obtaining the predicted temperature increment of the metal material unit by utilizing the discrete format equation.

Specifically, in this embodiment, all metal material units including the object points are heat transfer solving domains. The transient heat conduction equation satisfies the following relationship:

wherein ,for the material density->For the specific heat capacity of the material, D is differentiated, T is the temperature, T is the time variable, ++>Is a metal material volume fraction->Is metal material density->Is latent heat of phase change->For fluid volume fraction, ++>For gradient operator->Is heat conduction coefficient>Vaporization heat loss caused by vaporization of metal material, +.>For heat loss caused by convection heat exchange->For radiating heat loss>Is a volumetric heat source.

The heat loss due to evaporation and convection and the radiation heat loss respectively satisfy the following relations:

，

，

，

wherein ,for vaporization heat loss coefficient, < >>M is the molar mass of the metal material, R is the ideal gas constant, for the latent heat of evaporation, ++>For reference pressure +.>For evaporating temperature, ++>Is the specific heat capacity of the metal material +.>For the density of the ambient gas>Is the specific heat capacity of the ambient gas->Is a convection heat transfer coefficient>For reference temperature->Is Stefan-Boltzmann constant, & gt>Is the emissivity.

Volumetric heat sourceIs the sum of the volumetric heat sources of all the metallic material units. For the metallic material unit J, its volumetric heat source +. >The method comprises the following steps:

，

，

wherein ,is a surface heat source of the metal material unit J +.>For unit height, +.>Is the height of the surface of the metal material, +.>Is the cell height of the metallic material unit J, < ->Is the absorptivity of metal material to laser>For the laser power to be high,and r is the radial distance from the space point to the center of the laser spot.

Further, a discrete format equation obtained by solving a transient heat conduction equation by adopting a finite difference method satisfies the following relationship:

，

wherein ,temperature probe solution for the next temperature of the metallic material element,>for the predicted value of the next temperature of the metallic material element,/->For the current temperature of the metallic material element, +.>For predicting temperature increase, < > in->Predicted fluid volume fraction for next moment of metallic material element,/->The predicted temperature increment +.>。

And S513, adding the predicted temperature increment and the current temperature of the metal material unit to obtain the predicted value.

Specifically, in the present embodiment, the predicted value of the temperature at the next moment of the metal material unit satisfies the following relationship:

and S514, if the predicted value does not cross the temperature interval between the solidus temperature and the liquidus temperature, directly updating the lattice core temperature of the metal material unit by using the predicted value.

Specifically, in the present embodiment, if the predicted value of the next temperature of the metal material unit does not span the temperature interval between the solidus temperature and the liquidus temperature, it is considered that the solid-liquid phase does not occur in the metal material unit, and at this time, the predicted value of the next temperature of the metal material unit is directly used as the core temperature of the next moment of the metal material unit to update the current core temperature of the metal material unit.

And S515, if the predicted value spans the temperature interval, iteratively updating the grid core temperature by using a semi-implicit iteration method.

Specifically, in this embodiment, if the predicted value of the next temperature of the metal material unit spans the temperature interval between the solidus temperature and the liquidus temperature, it is considered that the metal material in the metal material unit may generate solid-liquid phase, and a new semi-implicit iteration method is adopted to calculate the next temperature of the lattice center of the metal material unit. The calculation principle of the new semi-implicit iteration method is as follows: first, whether the temperature in the metal material unit is increased or decreased is determined according to the positive or negative of the predicted temperature increase. The phase change latent heat only slows down the temperature change in one iteration calculation and does not change the trend of the temperature change, so that the initial upper limit and the initial lower limit of the temperature change can be determined according to the solidus temperature and the liquidus temperature of the metal material. First, predict times at a given temperature Number of digitsNext, the dichotomy is used to proceed +.>The secondary temperature trial results in a predicted temperature, which is specifically as follows:

，

，

，

wherein ,is->Predicted temperature from secondary temperature probe, +.>Is->Predicted temperature from secondary temperature probe, +.>Is->Predicted temperature from secondary temperature probe, +.>Is->Lower limit of temperature variation during sub-temperature probing, < ->Is->Upper limit of temperature variation during sub-temperature probing, < >>Is->Lower limit of temperature variation during sub-temperature probing, < ->Is->Upper limit of temperature variation during sub-temperature probing, < >>For metal material units at->Temperature heuristic solution of the current temperature at the time of the sub-temperature heuristic, < >>For metal material units at->Temperature heuristic solution of the current temperature at the time of temperature heuristic.

Further, it willPerforming secondary iterative calculation by taking the predicted temperature obtained by secondary temperature prediction as the initial value of Newton iterative method, and when +.>When the temperature obtained by Newton iteration calculation is used as the lattice temperature of the next moment of the metal material unit, otherwise, the temperature obtained by Newton iteration calculation is fed back to be used as the initial value of the Newton iteration method to continue the iteration calculation until the iteration is reachedUpper limit of number of calculations->。

Further, the method comprises the steps of,to converge the residual error->The temperature calculated for newton's iteration method.

S516, the cell core temperature is directly given to the object points in the metal material unit, so that the particle temperature of each object point in the background grid unit is updated.

Specifically, in this embodiment, the metal material unit is actually a background grid unit, and the material points are all present, so that the material point in the metal material unit is directly given the next temperature of the lattice center of the metal material unit, so that the particle temperature of each material point in the background grid unit can be updated.

S52, solving the particle speeds and the particle positions of the fluid particles and the paste region particles, and determining the flow behavior of the fluid in the fluid region under the heat driving according to the particle speeds and the particle positions.

Wherein S52 is performed in the fluid region, S52 further comprises the steps of:

s521, dividing each surface grid cell into a plurality of sub-cells, creating a virtual influence domain of material points in the surface grid cells, determining the occupancy rate of the material points in the surface grid cells to the sub-cells through the virtual influence domain, and taking the occupancy rate as the volume fraction of the surface grid cells.

Specifically, in this embodiment, referring to fig. 2, the area surrounded by the black solid line is a surface grid unit, each surface grid unit is divided into 16 sub-units, where the black dots are fluid particles or pasty area particles, and the dashed boxes with the black dots as geometric centers are virtual influence areas of the fluid particles or pasty area particles.

And judging whether the sub-unit is occupied by the fluid particles or the pasty area particles by calculating whether the grid center coordinates of the sub-unit are positioned in the virtual influence domain of the fluid particles or the pasty area particles, and finally obtaining the grid center volume fraction of the surface grid unit by calculating the occupancy rate of the fluid particles or the pasty area particles to the sub-unit.

The occupancy of fluid or paste-like area particles to a sub-cell is the ratio of all fluid and paste-like area particles occupying all sub-cells in a surface grid cell.

S522, solving the surface force on the surface grid unit by using the grid core volume fraction and adopting an interlaced derivative method, and simultaneously calculating the recoil pressure existing on the surface grid unit where vaporization and evaporation occur.

Wherein S522 further includes the following steps:

And S5221, calculating the lattice normal vector of the surface lattice unit at the lattice center according to the lattice center volume fraction.

Specifically, in the present embodiment, the cell normal vector of the surface mesh cell at the cell center satisfies the following relationship:

，

，

wherein ,for the lattice normal vector of the surface lattice cell J at the lattice center, +.>Cell volume fraction for construction using the central difference method>Gradient solver of>For any scalar physical quantity, use +.>Replacement->，、/> and />The length of the surface grid cells in x, y and z directions, respectively, +.>、/> and />Base vectors on x, y and z axes, respectively, i, j and k represent cell indices in the background grid cell in the x, y and z directions, respectively,/->、、/>、/>、/>、/>Is a different scalar physical quantity; />For gradient operators, i.e.)>。

Further, please refer to (a) in fig. 3, in the calculationWhen (I)>、/>、/>、、/>、/>The cell volume fraction for the different surface grid cells is shown as black circles in fig. 3 (a).

And S5222, interpolating the lattice core volume fraction to each node through a shape function N to obtain the node volume fraction.

Specifically, in this embodiment, the node volume fraction satisfies the following relationship:

wherein ,for the volume fraction of node I, +. >To influence the surface of the junction IA set of grid cells>The value of the shape function for the surface mesh unit J core at the junction I of the surface mesh unit.

S5223, calculating the gradient of the node volume fraction, and further obtaining a node normal vector.

Specifically, in this embodiment, the normal vector of the node satisfies the following relationship:

wherein ,is the normal vector of the node at node I, +.>Node volume fraction for construction using the center difference method>Gradient solver of>、/>、/>、/>、/>、/>The position of (b) is shown as a black dot in (b) of 3.

Further, referring to fig. 3 (b), the and calculation may be usedThe same method calculates->。

S5224, calculating the lattice center curvature of the surface lattice unit at the lattice center by using the node normal vector and the gradient of the shape function N.

Specifically, in the present embodiment, the cell center curvature of the surface mesh cell at the cell center satisfies the following relationship:

wherein ,is the curvature of the lattice center of the surface lattice unit J, < >>Is a set of nodes in the surface mesh unit J.

S5225, calculating the surface force according to the lattice center normal vector and the lattice center curvature, and simultaneously calculating the recoil pressure by using the lattice center normal vector.

Specifically, in the present embodiment, the surface force and the recoil pressure satisfy the following relationships, respectively:

，/>

，

wherein ,for the surface force of the surface grid cell J +.>For the tension coefficient>For fluid density->Is a metal material volume fraction->Recoil pressure for surface grid cell J +.>Is the recoil pressure coefficient.

Further, in calculating the surface force of the surface mesh unit, the temperature gradient of the surface mesh unit at the lattice centerCan use and calculate +.>Calculated by the same method.

And S523, calculating and updating the particle speeds of the fluid particles and the pasty region particles based on the surface force and the recoil pressure, and further calculating and updating the particle positions of the fluid particles and the pasty region particles by using the particle speeds.

Wherein S523 further comprises the following steps:

s5231, splitting an N-S equation containing Darcy damping terms by introducing an intermediate velocity field to obtain the intermediate velocity field of the nodes in the non-solid unit.

Specifically, in this embodiment, the N-S equation containing the Darcy damping term is decoupled based on the chord projection method.

The N-S equation containing the Darcy damping term is first split into the following two parts by introducing an intermediate velocity field:

，

，

wherein ,for the time step +.>For the intermediate velocity field introduced at node I in the non-solid element,/ >For the current node speed at node I, +.>Acceleration of gravity, ++>For shear viscosity coefficient, C is Darcy damping term,>for the current recoil pressure on the surface grid cell J, < >>For the current surface force on the surface grid cell J +.>For the next moment at node I in the non-solid element, < >>Is the cell core pressure at the next moment of the non-solid cell core.

And then, ignoring the Gecko pressure and the Darcy damping term to obtain an intermediate velocity field of the node I in the non-solid unit, namely, the intermediate velocity field introduced into the node I in the non-solid unit is:

wherein ,for the current recoil pressure at node I, +.>Is the current surface force at node I. /> and />By putting-> and />Mapping to node I results, and specific methods are referred to in the art and will not be described in detail herein.

Further, the method comprises the steps of,

，

，

wherein ,is the junction speed of junction I in the non-solid unit, < >>The velocity of the solid carried by node I at the fluid-solid interface, and (2)>For particle sets that have an influence on the junction I, < +.>Mass of object point->Is the speed of solid particles>Is the shape function value at the junction I in the metallic material unit where the substance point p is located. />

S5232, solving the lattice core pressure of the non-solid unit at the lattice core by using the intermediate speed field, and interpolating the lattice core pressure through a linear function to obtain the pressure gradient at each node in the non-solid unit.

Specifically, in the present embodiment, the pressure poisson equation can be obtained using the intermediate velocity field obtained in S5231:

wherein ,is the cell pressure of the non-solid unit J at the next moment.

In the pressure poisson equation, the second derivative term of the lattice heart pressure can be discretized into the following by adopting a three-dimensional center differential format:

wherein ,discrete form of the second derivative term of the lattice heart pressure in three-dimensional central differential format, L is characteristic length of non-solid unit, i, j and k respectively represent unit indexes along x, y and z directions in background grid unit, and->The next moment of the centroid pressure for the non-solid cell that is indexed in the x, y and z directions.

The poisson equation can then be written as a system of linear equations containing seven pressure unknowns:

wherein A represents a coefficient matrix, p is a vector of unknown lattice point pressures, and B represents a vector of intermediate velocity field divergences. By solving the linear equation set, the lattice center pressure of the lattice center in the next non-solid unit J can be obtained。

Further, it willPerforming linear function interpolation to obtain node pressure of each node in the non-solid unit at the next moment, wherein the method comprises the following steps of:

wherein ,The node pressure at the next moment of the node I.

Then according toUsage and computation->The pressure gradient +.about.at node I is calculated in the same way>。

S5233 updating the junction speed of each junction in the non-solid unit using the intermediate velocity field and the pressure gradient.

Specifically, in the present embodiment, the intermediate velocity field obtained in S5231 is usedAnd S5232 calculated pressure gradient +.>The node speed of the node I in the non-solid unit at the next moment can be calculated, the node speed of the node I in the non-solid unit is updated, and the node speed of the node I in the non-solid unit at the next moment meets the following relation:

wherein ,is the junction speed of the junction I at the next moment in the non-solid unit.

S5234, calculating the acceleration of each node in the non-solid unit by using the surface force, the recoil pressure and the node speed, and calculating and updating the particle speeds of the fluid particles and the paste-like region particles according to the acceleration, so as to update the particle positions of the fluid particles and the paste-like region particles by using the particle speeds.

Specifically, in the present embodiment, the acceleration of each node in the non-solid unit satisfies the following relationship:

wherein ,is the current acceleration of each node I in the non-solid unit.

Further, the particle velocities of the fluid particles and the paste-region particles satisfy the following relationship:

wherein ,for the next dot position of the fluid dot p or paste-like region dot p +.>The current particle position for a fluid particle p or paste region particle p; />The value range of the custom parameter is generally 0-1, and the custom parameter has better robustness when the value is 0.03. The particle velocities of the fluid particles and paste-like particles can be updated by this relationship.

Further, the particle positions of the fluid particles and paste-area particles satisfy the following relationship:

wherein ,for the dot position of the fluid dot p or paste-like region dot p at the next moment +.>The current particle position is for fluid particle p or paste region particle p.

S525, determining the flow behavior of the fluid in the fluid region under thermal drive according to the particle velocity and the particle position.

Since the fluid and paste-like regions are located in the fluid region, the behavior of fluid flow in the fluid region under thermal drive can be determined by the particle velocities and particle positions of the fluid and paste-like regions.

And S53, solving residual stress at the solid particles and total point force of background grid nodes in the solid region by adopting a standard substance particle method, and further solving solid deformation of the solid region.

Wherein, S53 further comprises the following steps:

s531, calculating the strain change rate at the solid particles by using a standard mass point method and the velocity gradient of the solid particles, and further calculating the residual stress at the solid particles by using an elastoplastic constitutive equation.

Specifically, in this embodiment, the rate of change of strain at the solid particlesThe following relationship is satisfied:

wherein ,is the velocity gradient of solid particles>Representation pair->Transpose, don't care>Is the linear expansion coefficient>Is the temperature change rate of solid particles +.>Is a second order unit tensor.

Then the residual stress at the solid particle can be calculated by means of elastoplastic constitutive equation, and the plastic flow followsFlow theory, its yield conditions are:

，

，

wherein ,is the cauchy stress in the solid region and is also the residual stress tensor at the solid particle; />Is equivalent to plastic strain>Is the yield stress. />Is a deflection stress tensor->The strain change rate, the stress rate, and the plastic flow direction of the second invariant of (a) satisfy:

，

，

，

wherein the residual strain tensor of the solid regionCan be made of->Obtained by time integration,/->Is the elastic strain rate tensor, +.>Is the plastic strain rate tensor; />For stress rate- >Is the first derivative of cauchy stress with respect to time; />Is plastic flow direction>Is the equivalent residual stress at the solid particles.

Further, in the plastic solving process, the yield surface of Mises criterion is a cylindrical surface in the main stress space, the reinforced yield surface is a geometric approximation thereof, and the method is thatThe plane is circular, and the normal vector is radial, so that the plastic correction is performed by adopting a radial return method.

Furthermore, the equivalent residual stress of the solid particles in the solid region can be obtained by solving the stepsTensor of residual stress->Residual Strain tensor->Wherein the residual stress tensor is convenient for observing the residual stress rule of the solid particles in all directions in the solid region. By equivalent residual stress->Tensor of residual stress->Residual Strain tensor->The prediction of the residual stress of the formed part can be realized.

S532, solving the total point force of the background grid nodes in the solid area by adopting a discrete solving mode of a momentum equation in a standard mass point method.

Specifically, in the present embodiment, it is first necessary to calculate the fluid-to-solid pressure by the following relationship:

，

，

wherein ,for the pressure of the fluid on the fluid-solid interface against the solid junction I, +. >Is the cell core pressure at the non-solid cell J cell core,>is the metal volume fraction value at the non-solid unit J cell center, +.>Is the volume of the non-solid unit J +.>Is the mass at the non-solid unit junction I.

The specific relational expression for solving the total point force of the background grid nodes in the solid area by adopting the discrete solving mode of the momentum equation in the standard mass point method is as follows:

wherein ,is an external force at the background grid node I in the solid area, < ->Is the internal force at the background mesh node I in the solid region,/>Is the total dot force of the background grid node I in the solid region.

Further, the method comprises the steps of, and />The following relationships are satisfied:

，

，/>

wherein ,，/>is the residual stress tensor of solid particle p +.>Is the volume represented by solid particle p.

And S533, updating the speed and the position of the solid particles by adopting a standard mass particle method according to the total node force, so as to obtain the solid deformation of the solid region.

Specifically, in this embodiment, the speed and position of the solid particles are updated according to the node force by using the standard mass point method as the prior art, and will not be described in detail herein.

S6, predicting the surface quality of the formed part under the selective laser melting technology according to the particle temperature, the flow behavior and the solid deformation.

Specifically, in this embodiment, the temperature change process of the metal material in the processing region may be obtained in real time according to the particle temperature, and the flow process of the metal material in the molten pool and the deformation of the unmelted metal material and the solidified region may be obtained in real time according to the flow behavior of the fluid in the fluid region under the heat drive. Thus, the surface quality of the shaped article can be predicted by the particle temperature, the flow behavior of the fluid in the fluid region, and the solid deformation of the solid region.

Further, since the residual stress at the solid particles is directly calculated in S53, the present embodiment can simultaneously realize the prediction of the surface quality and the residual stress of the formed article.

Further, in other alternative embodiments, after performing the operations S1 to S6 once, the operation frequency threshold may be set to repeatedly perform the steps S3 to S6, so that the surface quality and the residual stress of the formed article may be continuously predicted, a material domain described by a series of material points may be obtained, the forming morphology of the cladding channel may be easily observed, and further the surface quality of the formed article may be quantitatively evaluated by using quantitative indexes such as surface roughness. Meanwhile, the residual stress of the cladding channel and the substrate area can be obtained from the residual stress carried on the object points.

It should be noted that, in some cases, the actions described in the specification may be performed in a different order and still achieve desirable results, and in this embodiment, the order of steps is merely provided to make the embodiment more clear, and it is convenient to describe the embodiment without limiting it.

In summary, the temperature change of the metal material unit is obtained by calculating the particle temperature of each object particle in the background grid unit, the flow behavior of the fluid in the fluid area under the heat driving is determined by calculating the particle speeds and particle positions of the fluid particle and the paste area particle, the solid deformation is obtained by calculating the total dot force of the background grid node in the solid area, meanwhile, the residual stress of the formed piece is directly calculated, and the surface quality of the formed piece is predicted by utilizing the temperature change of the metal material unit, the flow behavior of the fluid in the fluid area under the heat driving and the solid deformation of the solid area, so that the surface quality and the residual stress of the formed piece under the selective laser melting technology are predicted simultaneously.

Finally, it should be noted that: the above embodiments are only for illustrating the technical solution of the present invention, and not for limiting the same; although the invention has been described in detail with reference to the foregoing embodiments, it will be understood by those of ordinary skill in the art that: the technical scheme described in the foregoing embodiments can be modified or some or all of the technical features thereof can be replaced by equivalents; such modifications and substitutions do not depart from the spirit of the invention, and are intended to be included within the scope of the appended claims and description.

Claims (10)

1. The method for predicting the surface quality and residual stress of the formed part under the selective laser melting technology is characterized by comprising the following steps of:

modeling a powder bed laying process by using a discrete element method to obtain an initial powder particle distribution model in a selective laser melting process;

dividing a calculation grid by using preset calculation parameters, dispersing the whole material domain into object points, and initializing the calculation grid by using the object points;

searching a background grid in the initialized calculation grid according to the quality field of the background grid node in the object point method;

formulating a classification marking strategy, and classifying and marking the background grid cells and the substance points in the background grid by combining the distribution model;

calculating and updating particle temperature of the object particles based on the result of the classification mark, determining flow behavior of fluid in the fluid region under the drive of heat, and solving residual stress and solid deformation of the solid region;

the surface quality of the shaped part under the selective laser melting technology is predicted according to the particle temperature, the flow behavior and the solid deformation.

2. The method for predicting surface quality and residual stress of a formed part by selective laser melting according to claim 1, wherein the classification marking strategy comprises the steps of:

Classifying and marking the background grid cells by utilizing physical information carried by neighbor grid cells of the background grid cells;

and classifying and marking the material points according to the physical information carried by the material points and the distribution model.

3. The method for predicting surface quality and residual stress of a formed part under a selective laser melting technology according to claim 2, wherein the classifying and marking the background grid cells by using physical information carried by the background grid cells and neighbor grid cells of the background grid cells comprises the following steps:

marking a background grid cell which does not contain object points as an empty cell, otherwise marking the background grid cell as a metal material cell;

marking the metal material units of which the neighbor grid units are the empty units as surface grid units;

marking a metal material unit with a temperature less than the solidus temperature as a solid unit;

the metallic material units having a temperature between the solidus temperature and the liquidus temperature are labeled as paste-like zone units.

4. A method for predicting surface quality and residual stress of a molded part by selective laser melting according to claim 3, wherein said classifying and marking said material points according to physical information carried by said material points and said distribution model comprises the steps of:

Marking material points in the powder particles as unmelted material points and marking the rest material points as solid material points according to the distribution model;

marking the unmelted points of material having a temperature greater than the liquidus temperature as fluid particles, and marking the unmelted points of material having a temperature between the solidus temperature and the liquidus temperature as paste-like region particles;

the material points with the temperature exceeding the solidus temperature are named as intermediate material points, the solid stress carried by the intermediate material points with the temperature not lower than the liquidus temperature is set to zero and marked as fluid material points, the solid stress carried by the intermediate material points with the temperature between the solidus temperature and the liquidus temperature is set to zero and marked as pasty area material points, and the intermediate material points with the temperature lower than the solidus temperature are marked as solid material points.

5. The method for predicting the surface quality and residual stress of a formed part under the selective laser melting technology according to claim 4, wherein the method comprises the following steps of:

the metal material units except the solid unit are non-solid units, the non-solid units comprise pasty area units, the area formed by the non-solid units is the fluid area, and the area formed by the solid units is the solid area;

The method for calculating and updating the particle temperature of the object particles based on the result of the classification mark, determining the flow behavior of the fluid in the fluid region under the drive of heat, and solving the residual stress and the solid deformation of the solid region comprises the following steps:

determining a heat transfer solution domain within the background grid based on the locations of the object points, and calculating and updating a particle temperature for each object point in the background grid cell within the heat transfer solution domain;

solving the particle velocity and the particle position of the fluid particles and the paste region particles, and determining the flow behavior of the fluid in the fluid region under the heat driving according to the particle velocity and the particle position;

and solving residual stress at the solid particles and total point force of background grid nodes in the solid region by adopting a standard mass point method, so as to solve solid deformation of the solid region.

6. The method of claim 5, wherein determining a heat transfer solution domain based on the locations of the object particles in the background grid, and calculating and updating the particle temperature of each object particle in the background grid cell in the heat transfer solution domain comprises the steps of:

Determining whether a new metal material unit appears in a calculation domain, and initializing the temperature of the new metal material unit;

determining the heat transfer solving domain by using the positions of object points in all the metal material units, circulating all the metal material units in the heat transfer solving domain, solving a transient heat conduction equation of the metal material units by adopting a finite difference method to obtain a discrete format equation of the transient heat conduction equation, and obtaining a predicted temperature increment of the metal material units by utilizing the discrete format equation, wherein the transient heat conduction equation and the discrete format equation respectively satisfy the following relations:

，

，

wherein ,for the material density->For the specific heat capacity of the material, D is differentiated, T is the temperature, T is the time variable, ++>Is a metal material volume fraction->Is metal material density->Is latent heat of phase change->For fluid volume fraction, ++>For the purpose of the gradient operator,is heat conduction coefficient>Vaporization heat loss caused by vaporization of metal material, +.>For heat loss due to convective heat transfer,for radiating heat loss>Is a volumetric heat source->Temperature probe solution for the next temperature of the metallic material element, < >>For the predicted value of the temperature at the next moment of the metallic material unit,/- >For the current temperature of the metallic material element, < >>For the predicted temperature increase,/->Is the metal materialPredicted fluid volume fraction at next cell time, +.>A current fluid volume fraction for the metallic material unit;

adding the predicted temperature increment and the current temperature of the metal material unit to obtain the predicted value;

if the predicted value does not span the temperature interval between the solidus temperature and the liquidus temperature, directly updating the lattice core temperature of the metal material unit by using the predicted value;

if the predicted value spans the temperature interval, iteratively updating the lattice core temperature by using a semi-implicit iteration method;

and directly giving the cell core temperature to the object points in the metal material unit, and further updating the particle temperature of each object point in the background grid unit.

7. The method of claim 6, wherein said solving for the particle velocity and particle position of the fluid particles and paste-like region particles and determining the flow behavior of the fluid in the fluid region under thermal drive based on the particle velocity and particle position comprises the steps of:

Dividing each surface grid cell into a plurality of sub-cells, creating a virtual influence domain of object points in the surface grid cells, determining the occupancy rate of the object points in the surface grid cells to the sub-cells through the virtual influence domain, and taking the occupancy rate as the grid core volume fraction of the surface grid cells;

solving the surface force on the surface grid unit by using the grid core volume fraction and adopting a staggered derivative method, and simultaneously calculating the recoil pressure existing on the surface grid unit where vaporization and evaporation occur;

calculating and updating particle velocities of the fluid particles and the paste-like region particles based on the surface force and the recoil pressure, and further calculating and updating particle positions of the fluid particles and the paste-like region particles using the particle velocities;

determining a flow behavior of the fluid in the fluid region under thermal drive based on the particle velocity and the particle location.

8. The method for predicting surface quality and residual stress of a formed part by selective laser melting according to claim 7, wherein said calculating the recoil pressure existing on the surface grid unit in which vaporization and evaporation occur while solving the surface force on the surface grid unit by using the cell volume fraction and using an interleaved derivative method comprises the steps of:

Calculating a grid normal vector of the surface grid unit at the grid center according to the grid center volume fraction;

interpolating the grid core volume fraction to each node through a shape function N to obtain a node volume fraction;

calculating the gradient of the node volume fraction, and further obtaining a node normal vector;

calculating the lattice center curvature of the surface grid unit at the lattice center by using the node normal vector and the gradient of the shape function N;

calculating the surface force according to the lattice center normal vector and the lattice center curvature, and simultaneously calculating the recoil pressure by using the lattice center normal vector, wherein the surface force and the recoil pressure respectively meet the following relations:

，

，

wherein ,for said surface force of the surface grid cell J +.>For the tension coefficient>For the curvature of the lattice center of the surface lattice unit J, < >>D represents the differential operation for the lattice normal vector of the surface lattice unit J, T is the temperature, < >>For gradient operator->For fluid density->Is metal material density->For the density of the ambient gas>Is a metal material volume fraction->For the recoil pressure of the surface grid cell J, +.>For recoil pressure coefficient, +.>For reference pressure +.>For evaporating latent heat->Is the molar mass of the metallic material, +. >For evaporating temperature, ++>Is an ideal gas constant.

9. The method of claim 8, wherein the calculating and updating the particle velocities of the fluid particles and the paste particles based on the surface force and the recoil pressure, and further using the particle velocities to calculate and update the particle positions of the fluid particles and the paste particles comprises the steps of:

splitting an N-S equation containing a Darcy damping term by introducing an intermediate velocity field to obtain an intermediate velocity field of a node in the non-solid unit;

solving the lattice core pressure of the non-solid unit at the lattice core by utilizing the intermediate speed field, and interpolating the lattice core pressure through a linear function to obtain the pressure gradient at each node in the non-solid unit;

updating the junction speeds of the junctions in the non-solid unit using the intermediate velocity field and the pressure gradient;

and calculating the acceleration of each node in the non-solid unit by using the surface force, the recoil pressure and the node speed, and calculating and updating the particle speeds of the fluid particles and the paste-like region particles according to the acceleration, so as to update the particle positions of the fluid particles and the paste-like region particles by using the particle speeds.

10. The method for predicting surface quality and residual stress of a shaped article by selective laser melting according to claim 9, wherein said solving for residual stress at said solid particles and total point forces of background grid nodes in said solid region by standard mass point method, and further solving for solid deformation of said solid region comprises the steps of:

calculating the strain change rate at the solid finger points by means of the velocity gradient of the solid point by using a standard mass point method, and further calculating the residual stress at the solid point by means of an elastoplastic constitutive equation;

solving the total point force of the background grid nodes in the solid area by adopting a discrete solving mode of a momentum equation in a standard mass point method;

and updating the speed and the position of the solid particles by adopting a standard mass point method according to the total point force, so as to obtain the solid deformation of the solid region.