A ceramic injection molding multiphase flow calculation method, device and medium based on lattice boltzmann and discrete element method

Abstract

The present application belongs to the technical field of ceramic injection molding, and particularly relates to a ceramic injection molding multiphase flow calculation method, equipment and medium based on lattice Boltzmann and discrete element method. The present application establishes a two-way coupling model of fluid and particles, calculates local porosity in real time in simulation, and corrects fluid evolution equation, so as to accurately characterize the heterogeneous flow characteristics of high solid content feedstock. The method can simulate the micro behaviors such as particle migration and agglomeration, predict the density gradient distribution of green body, and handle the flow mutation and blockage risk under complex mold geometry. The present application simulates the interaction between the binder and the powder in the feedstock, realizes accurate prediction of heterogeneous flow, shear-induced migration and phase separation, and also realizes accurate capture of the interaction between solid and liquid phases, flow front evolution and micro defect formation mechanism in the injection molding process, providing reliable simulation basis for process optimization and quality control.

Description

Technical Field

[0001] This invention belongs to the field of ceramic injection molding technology, specifically relating to a method, equipment, and medium for calculating multiphase flow in ceramic injection molding based on the lattice Boltzmann and discrete element method. Background Technology

[0002] Ceramic injection molding, as a near-net-shape forming technology, is widely used in manufacturing ceramic parts with complex shapes, high dimensional accuracy, and specific microstructures. The core of this process lies in mixing ceramic powder with an organic binder to form a flowable feedstock, which is then injected into the mold cavity under high pressure.

[0003] The simulation accuracy of fluid flow under high solid content conditions is insufficient. Ceramic powder accounts for a high volume proportion in ceramic injection molding feedstock, which belongs to a multiphase flow system with high solid content. Existing simulation technologies mostly use fluid evolution equations for low solid content scenarios without making specific corrections for the characteristics of dense particle packing. This makes it impossible to accurately quantify the obstructive effect of particle packing on fluid flow, resulting in the simulated fluid velocity and shear rate not matching the actual working conditions, and making it difficult to capture fluid flow anomalies under high solid content.

[0004] Existing technologies do not take into account the dynamic changes in the adhesion force and fluid-structure interaction force between particles when calculating particle motion parameters. This leads to deviations in the updating of particle translational velocity, rotational velocity, and position and orientation, making it impossible to accurately reproduce the migration, collision, and aggregation behavior of particles during injection molding. Consequently, it affects the accuracy of subsequent porosity calculations and equation corrections.

[0005] Existing simulation technologies often use structured meshes to discretize mold boundaries, which makes it difficult to accurately fit the contours of complex molds and easily leads to boundary approximation errors. Existing technologies output flow field or particle distribution data, but cannot predict potential defects such as uneven density inside the green preform in advance; they also lack the ability to capture abnormal phenomena such as local rheological changes and particle blockage during injection molding. Summary of the Invention

[0006] This invention provides a calculation method for multiphase flow in ceramic injection molding based on lattice Boltzmann and discrete element method. It simulates the interaction between binder and powder in the feed and realizes accurate prediction of phenomena such as heterogeneous flow, shear-induced migration and phase separation. It provides quantitative basis for optimizing mold design, injection parameters and predicting sintering defects.

[0007] The methods include: S1: Establish a fluid-particle bidirectional coupling model to simulate the flow of ceramic feed in the mold; S2: Based on the fluid-particle bidirectional coupling model, input the initial ceramic powder particle parameters, binder fluid parameters and mold cavity boundary conditions, and perform transient simulation calculations on the injection molding filling process of ceramic feeding to obtain the evolution data of particle position, velocity and fluid flow state. S3: During the simulation calculation, for each time step, based on the current fluid flow state obtained in step S2, the fluid-structure interaction force acting on each solid particle is calculated using the momentum exchange method. S4: Based on the fluid-structure interaction force, interparticle contact force and gravity calculated in step S3, and using Newton's second law, calculate the translational velocity and rotational velocity of each solid particle in the next time step, and update its position and orientation. S5: Based on the updated positions of all solid particles in step S4, calculate the local porosity at each spatial location within the simulation area; and based on the local porosity, correct the evolution equation of the lattice Boltzmann method in step S2 to calculate fluid flow under high solid content conditions. S6: Based on the simulation results of steps S2 to S5, the final spatial distribution of solid particles in the mold cavity when the injection molding process is completed is statistically analyzed, and the density gradient distribution inside the green blank is predicted based on the uniformity of the particle distribution. S7: During the simulation process in step S2, the geometric boundary conditions of the micro or complex mold shape are processed by the lattice Boltzmann method. Combined with the force and motion calculations in steps S3 to S5, local rheological abrupt changes and potential particle blockage locations are captured, and finally the flow field, particle distribution and defect prediction results are output for evaluating the injection molding quality.

[0008] According to another embodiment of this application, an electronic device is provided, including a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the program to implement the steps of the ceramic injection molding multiphase flow calculation method based on the lattice Boltzmann and discrete element method.

[0009] According to another embodiment of this application, a storage medium is also provided, on which a computer program is stored, which, when executed by a processor, implements the steps of the ceramic injection molding multiphase flow calculation method based on the lattice Boltzmann and discrete element method.

[0010] As can be seen from the above technical solutions, the present invention has the following advantages: This invention provides a multiphase flow calculation method for ceramic injection molding based on the lattice Boltzmann method and the discrete element method. It establishes a fluid-particle bidirectional coupling model, using LBM for the fluid phase and DEM for the solid particle phase. The LBGK evolution equation and JKR contact force model are explicitly selected to adapt to the characteristics of the ceramic feedstock. Based on the momentum exchange mechanism, the method sets rules for selecting fluid grid points on the particle surface, momentum transfer scale transformation coefficients, and coupling timing, and prioritizes the fluid-solid interaction calculations. This enables real-time feedback of the bidirectional fluid-solid interaction, capturing the discrete motion, collision, and aggregation behavior of particles, making the modeling results more consistent with the multiphase characteristics of the ceramic feedstock, and resulting in a more reliable simulation foundation.

[0011] Based on particle location, the local porosity of each fluid grid point is calculated, and outliers are removed through boundary smoothing. Then, combining porosity data and the Ergun drag coefficient model, a drag coefficient correction value is calculated to determine the porosity-related force term, which is integrated into the original LBM evolution equation and replaces the original external force term. This allows the fluid evolution equation to adapt in real-time to the particle packing state under high solids content, accurately quantifying the hindering effect of particles on fluid flow and improving the accuracy of fluid flow simulation under high solids content conditions.

[0012] At each time step, the momentum exchange method is used to select fluid grid points within 5 layers of the particle surface, calculate the momentum difference, and accumulate them to obtain the fluid-structure interaction force, thus verifying abnormal data. Integrating the fluid-structure interaction force, interparticle contact force, and gravity, and based on Newton's second law, the Euler forward integration method is used to calculate particle motion parameters. After updating the position and attitude, the boundary and abnormal data are verified. This method can characterize the six-degree-of-freedom motion characteristics of particles, match the dynamic changes of the fluid-structure state in real time, reduce the deviation between particle motion parameters and position updates, and realistically reproduce the migration and aggregation states of particles during injection molding.

[0013] After the simulation terminates, the mold cavity is divided into three-dimensional sub-regions, the final particle distribution is statistically analyzed, the solid volume fraction is verified, and the green body density gradient is predicted based on the actual packing density of the ceramic powder. Fluid shear rate, porosity, and particle packing density are monitored in real time, and thresholds are set to capture rheological abrupt changes and particle blockage. After the simulation terminates, all simulation data are integrated and organized into a visualization package according to evaluation indicators. This allows for the early prediction of potential defects such as uneven density within the green body and the real-time identification of process anomalies during injection molding. Attached Figure Description

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

[0015] Figure 1Flowchart of the calculation method for multiphase flow in ceramic injection molding; Figure 2 Flowchart of an embodiment of the calculation method for multiphase flow in ceramic injection molding; Figure 3 This is a schematic diagram of an electronic device. Detailed Implementation

[0016] The multiphase flow calculation method for ceramic injection molding based on the lattice Boltzmann and discrete element method provided by this invention establishes a direct interaction model between fluid and particles at the mesoscale through bidirectional coupling of LBM-DEM. This enables the capture of shear-induced particle migration occurring at high shear rates during injection molding, as well as solid-liquid phase separation at corners and gates. This method can obtain a dynamic distribution map of particles in the flow field, thereby identifying potential powder-deficient or powder-rich regions in advance.

[0017] The ceramic injection molding of this invention is commonly used to manufacture miniature parts or parts with complex internal cavities. This invention utilizes the inherent advantages of the lattice-Boltzmann method in parallel computing and handling complex boundaries, combined with the discrete element method for accurate description of particle contact mechanics, to more realistically reproduce the clogging, bridging effect, and wall slip behavior of ceramic powder in narrow flow channels, thereby improving the prediction accuracy of microscopic defects.

[0018] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0019] like Figure 1 and Figure 2 As shown, the method includes the following steps: S1: Establish a fluid-particle bidirectional coupling model for simulating the flow of ceramic feed in a mold, wherein the fluid phase is modeled using the lattice Boltzmann method and the solid particle phase is modeled using the discrete element method, and define the momentum exchange mechanism between the fluid phase and the solid particle phase.

[0020] Step S1 specifically includes the following methods: S11. Establish and initialize a solution containing a lattice Boltzmann solver and a discrete element method solver. Define a computational mesh covering the mold cavity for the lattice Boltzmann solver, configure the fluid relaxation time τ according to formula (1), and back-calculate and confirm the value of τ according to formula (2) by using the set kinematic viscosity ν of the binder.

[0021] The spatiotemporal evolution of the binder fluid distribution function is described as follows: ; In the formula, f i For position x, time t, along discrete velocity direction c i The fluid particle distribution function characterizes the microscopic state of the binder fluid. x represents the spatial coordinates, corresponding to the grid point positions in the fluid computational grid. c i Δt is the discrete velocity vector, representing the allowed direction of motion of fluid particles within the lattice. Δt is the time step in the LBM calculation. τ is the relaxation time, controlling the rate at which the fluid distribution function tends towards local equilibrium, and is related to the fluid viscosity. It is: the local equilibrium distribution function, which is determined by the local density and the macroscopic velocity.

[0022] This is the component of the external force term in the discrete velocity direction i, used to introduce macroscopic forces such as gravity, pressure gradient, or particle interaction.

[0023] Establish the relationship between macroscopic physical properties and microscopic parameters: ; ν is the kinematic viscosity of the binder fluid, a macroscopic physical property. is the lattice sound velocity, a constant related to the selected discrete velocity model.

[0024] The initial position, particle size distribution, density, and material parameters of the ceramic powder particles are imported into the discrete element method solver, and the equations of motion are set according to formulas (6) and (7). A shared data exchange interface is established for the two solvers.

[0025] The method for describing the positional changes of ceramic powder particles is as follows: ; In the formula, m k Let dv be the mass of the kth ceramic powder particle. k / dt is the translational acceleration of the k-th particle. The vector sum of all contact forces acting on the k-th particle, F f,k Let g be the force exerted by the fluid on the k-th particle, and g be the gravitational acceleration vector.

[0026] The method for describing the angle change of ceramic powder particles is as follows: ; In the formula, I k Let be the moment of inertia of the k-th ceramic powder particle. Let be the angular acceleration of the k-th particle. The resultant torque produced by the contact force on the k-th particle. Let be the torque generated by the fluid force on the k-th particle.

[0027] In an exemplary embodiment, the lattice Boltzmann solver generates an Eulerian background mesh based on the three-dimensional geometry of the mold cavity, the size of which determines the spatial resolution of the simulation. The discrete element method solver creates a set containing all discrete ceramic powder particles, each particle being treated as an independent entity with mass, moment of inertia, and initial state of motion. The determination of the kinematic viscosity ν and relaxation time τ established by formula (2) is a key preliminary operation to ensure the physical accuracy of the fluid dynamics behavior, directly transforming the macroscopic properties of the material into microscopic simulation parameters.

[0028] S12. In the computational grid, identify all grid cells occupied by solid particles and mark them as fluid-structure interaction interfaces.

[0029] At each simulation time step, after the lattice Boltzmann solver completes the evolution of the fluid distribution function according to formula (1), at the coupling interface, the resultant force and resultant torque of the fluid acting on the surface of each particle are calculated according to the momentum exchange method shown in formula (9).

[0030] ; In the formula, F f,k Let be the total force exerted by the fluid on the k-th particle. For all boundary grid points x associated with the particle k surface b Summation, To sum over all discrete velocity directions i pointing into the particle, This is the distribution function from the fluid towards the lattice points at the particle boundary before the collision. The distribution function of the fluid rebounding from the particle boundary after a collision. Let Δx be the discrete velocity vector, and Δx be the grid spacing.

[0031] For regions with high solids content, the fluid evolution is locally porosity corrected according to formula (11), and the particle drag force after local porosity correction is calculated according to formulas (12) and (13).

[0032] ; In the formula, This is an additional force term, related to the local porosity ϵ, used to characterize the additional distributed resistance to fluid flow generated by solid particles under high solid content conditions.

[0033] ; In the formula, U represents the drag force acting on the particles, β is the drag coefficient, and is a function of local porosity, relative velocity, etc. f vp represents the local fluid velocity, and vp represents the particle velocity.

[0034] ; In the formula, ϵ is the local porosity, μ is the hydrodynamic viscosity, and d p Where ρ is the characteristic diameter of the ceramic powder particles, and ρ is the fluid density, as in Formula 4. This represents the relative velocity between the fluid and the particles.

[0035] In one exemplary embodiment, the identification of the coupling interface is a real-time geometric search process that determines which fluid grid points are cut or occupied by the moving particle surfaces. The momentum exchange method of Equation (9) is performed on these interface grid points, directly statistically analyzing the change in fluid molecule momentum caused by the presence of particles, and then integrating to obtain the total fluid force on the particles. Meanwhile, Equation (11) introduces an additional term related to the local particle concentration in the lattice Boltzmann evolution. This force term explicitly expresses the additional volumetric impediment effect of the particle swarm on the fluid flow, which is a necessary modification to the standard model under the sparsity assumption. Equations (12) and (13) provide a drag force calculation method based on the idea of ​​local averaging, whose coefficients strongly depend on the local pore structure composed of particle swarms.

[0036] S13. The forces calculated in step S12 and the fluid forces based on the modified drag force are transmitted to the discrete element method solver.

[0037] The discrete element method solver combines the fluid forces with the calculated interparticle contact forces and gravity, substitutes them into equations (6) and (7), and integrates to update the velocity, angular velocity, position, and orientation of all ceramic powder particles. The updated particle position and velocity information is fed back to the lattice Boltzmann solver for re-identification of the coupling interface in the next time step, thereby closing the bidirectional coupling loop.

[0038] In one exemplary embodiment, the order, content, and timing of data exchange between the parallel fluid solver and particle solver are specified. Fluid force data transmitted from the lattice Boltzmann solver is added to the resultant force of each corresponding particle in the discrete element solver. The discrete element solver then solves a set of coupled equations (6) and (7), updating the particle motion state using an explicit time integration method. The updated particle position information is immediately fed back to the fluid solver. This feedback allows the solid boundary representation in the fluid mesh to be updated in real time, thereby affecting the evolution of equation (1) in the next cycle, the interface grid point determination in equation (9), and the recalculation of porosity in equation (11).

[0039] S2: Based on the fluid-particle bidirectional coupling model, input the initial ceramic powder particle parameters, binder fluid parameters and mold cavity boundary conditions, and perform transient simulation calculations on the injection molding filling process of ceramic feeding to obtain the evolution data of particle position, velocity and fluid flow state.

[0040] In an exemplary embodiment, based on the bidirectional coupling model established by S1, three types of core initial parameters are input in batches to ensure that the parameters are consistent with the actual ceramic injection molding process.

[0041] The parameters of ceramic powder particles include: particle size distribution, particle density, interparticle adhesion coefficient, and elastic recovery coefficient.

[0042] The binder fluid parameters include kinematic viscosity, density, and surface tension. The relaxation time parameter is adjusted by formula (2) to match the binder viscosity in the model with the actual binder formulation.

[0043] ; The boundary conditions of the mold cavity include: the three-dimensional geometric model of the mold, the surface roughness of the cavity, the inlet injection pressure, the inlet flow rate, and the simulation termination conditions.

[0044] The transient simulation employs parallel computing with a fixed time step. The selection of the time step balances computational accuracy and efficiency, ensuring that subtle changes in particle motion and fluid flow can be captured. During the simulation, the evolution data for each time step is recorded in real time, including: the three-dimensional spatial coordinates of the ceramic particles, translational velocity, rotational velocity, and inter-particle contact state.

[0045] Based on the macroscopic density field of the binder fluid involved in formula (4), the macroscopic velocity field of formula (5), the fluid distribution function, and the shear rate, all evolution data are classified and archived according to time steps to form a simulation database, and the validity of the data is verified synchronously.

[0046] Calculating macroscopic density and velocity from microscopic distribution functions: ; Formula 4 is the statistical formula for macroscopic fluid density. At any spatial point, the macroscopic density of the fluid can be obtained by simply summing the microscopic distribution functions along all possible directions of motion at that point. This is one of the fundamental operations for extracting macroscopic variables from discrete, statistical particle distribution functions in continuum mechanics.

[0047] ; Equation 5 is used to calculate macroscopic fluid velocity. The principle of Equation 5 is to multiply the distribution function in all directions by their respective discrete velocity vectors and then sum the results; the result is proportional to the fluid's momentum density. An additional external force correction term is added to the formula to ensure the accuracy of the macroscopic velocity definition in the presence of an external force field, which is a way for the lattice Boltzmann method to achieve a correct physical response.

[0048] S3: During the simulation, for each time step, based on the current fluid flow state obtained in step S2, the fluid-structure interaction force acting on each solid particle is calculated using the momentum exchange method.

[0049] Step S3 specifically includes the following steps: S31: Identify the fluid computation grid points covering the surface of each solid particle and obtain the evolution data of the fluid flow state at the previous and current times of the fluid computation grid points.

[0050] In an exemplary embodiment, for each solid particle, based on its current position and radius, the surrounding fluid is traversed to calculate grid points. It is determined whether the grid point is within the influence range of the particle surface, i.e., the distance from the grid point center to the particle center is less than the sum of the particle radius and half the grid step size. For grid points that meet the condition, the distribution function of the previous moment and the distribution function of the current moment obtained in step S2 are extracted, and the distribution function value of each discrete velocity direction is recorded.

[0051] S32: For each identified fluid calculation grid point, calculate the momentum difference between the fluid and the particle boundary before and after the collision at the grid point according to formula (9), accumulate the momentum differences of all relevant grid points and perform scale transformation to obtain the preliminary fluid-structure interaction force.

[0052] In one exemplary embodiment, for each identified fluid grid point, the distribution functions of the previous time step and the current time step are substituted into: ; The momentum difference in each discrete velocity direction is calculated, which is the distribution function multiplied by the discrete velocity vector. The momentum differences in all discrete velocity directions are summed to obtain the momentum change at that lattice point. The momentum changes of all relevant lattice points are accumulated and multiplied by the scale factor of the lattice step size and time step to obtain the initial fluid-structure interaction force acting on the particle. Here, the force is calculated from microscopic data of fluid evolution, without the need for empirical assumptions, ensuring consistency between the force calculation and fluid dynamics.

[0053] S33: Based on the calculated local porosity, the initial fluid-structure interaction force is corrected by calculating the drag coefficient after high solids content correction, and the final fluid-structure interaction force is obtained.

[0054] In an exemplary embodiment, the calculation method for the local porosity calculated in step S2 is as follows: ; In the formula, ϵ(x) is the local porosity at position x. δ(x−xi) is the Dirac delta function, used to determine whether particle i affects position x. xi is the center position of the i-th particle. V k This refers to the volume of a single particle.

[0055] Substitute it into formula (13) to calculate the drag coefficient under high solid content.

[0056] Substituting the drag coefficient into formula (12), the drag correction value corresponding to the relative velocity between the particles and the fluid is calculated. The correction value is then superimposed with the preliminary fluid-structure interaction force obtained in step S32 to obtain the final fluid-structure interaction force. Here, the corrected fluid-structure force better matches the actual flow resistance in a high-solids-content powder bed, avoiding errors caused by the low-concentration assumption and improving the accuracy of the coupling calculation.

[0057] S4: Based on the fluid-structure interaction force, interparticle contact force and gravity calculated in step S3, and using Newton's second law, calculate the translational velocity and rotational velocity of each solid particle in the next time step, and update its position and orientation.

[0058] In an exemplary embodiment, the fluid-structure interaction force matrix output in step S3 is retrieved, the drag force and lift force corresponding to each particle are extracted, and other forces acting on each particle are calculated: the contact force between particles is calculated by the JKR contact force model formula (8), and the adjacent particles around each particle are traversed to determine whether the distance between particles is less than the sum of the radii of the two particles.

[0059] Describe the collisions between particles and the adhesion caused by the binder: ; In the formula, F normal This refers to the normal contact force between particles. The equivalent elastic modulus is calculated from the elastic modulus and Poisson's ratio of the two particles. The equivalent radius is calculated from the radii of the two particles, δ is the normal overlap between the particles, and Δγ is the surface energy parameter, which characterizes the strength of the adhesion effect caused by the binder or the material itself.

[0060] If the contact conditions are met, substitute into formula (8) to calculate the normal contact force and decompose to obtain the tangential contact force; the particle's own weight is calculated based on the particle's mass and gravitational acceleration, and its direction is vertically downward.

[0061] The fluid-structure interaction force, interparticle contact force, and gravity are vector superimposed to obtain the resultant force on each particle. This resultant force is then substituted into the particle translation equation (6), and the particle's translational acceleration is calculated by dividing the resultant force by the particle's mass. Using the Euler forward integral method, the particle's translational velocity (current velocity + acceleration × time step) is calculated for the next time step. The particle's displacement is then obtained by multiplying the translational velocity by the time step, and vector superimposed with the particle's current three-dimensional spatial coordinates to update the particle's position.

[0062] Furthermore, the contact torque generated by the tangential contact force between particles and the fluid torque generated by the fluid-structure interaction force are vector-superimposed to obtain the resultant torque on each particle. Substituting this into the particle rotation equation formula 7, the angular acceleration is calculated by dividing the resultant torque by the particle's moment of inertia. Using the same integration method, combined with the current angular velocity, the angular velocity of the next time step is calculated. The rotation angle of the particle is obtained by multiplying the angular velocity by the time step, and the spatial attitude of the particle is updated.

[0063] Furthermore, after the update is completed, it is checked whether the particle position exceeds the boundary of the mold cavity. If it does, it is adjusted to the inside of the boundary. The particle velocity and angular velocity are checked for abnormality. Unreasonable data is removed and recalculated. Based on formulas (6), (7), and (8), which reflect the effects of all forces on the particle, the accuracy of the motion parameter calculation is ensured, and the rapid iterative update of particle motion parameters can be realized.

[0064] S5: Based on the updated positions of all solid particles in step S4, calculate the local porosity at each spatial location within the simulation area; and based on the local porosity, correct the evolution equation of the lattice Boltzmann method in step S2 to calculate the fluid flow under high solid content conditions.

[0065] S6: Based on the simulation results of steps S2 to S5, the final spatial distribution of solid particles in the mold cavity when the injection molding process is completed is statistically analyzed, and the density gradient distribution inside the green blank is predicted based on the uniformity of the particle distribution.

[0066] S7: During the simulation process in step S2, the geometric boundary conditions of the micro or complex mold shape are processed by the lattice Boltzmann method. Combined with the force and motion calculations in steps S3 to S5, local rheological abrupt changes and potential particle blockage locations are captured, and finally the flow field, particle distribution and defect prediction results are output for evaluating the injection molding quality.

[0067] In an exemplary embodiment, during the transient simulation in step S2, mold boundary processing and anomaly capture operations are performed simultaneously. For irregular cavities and fine channels, unstructured meshes are used to discretize the mold cavity boundaries to ensure accurate fit to the mold contour. The boundary grid points are divided into three categories: solid wall boundaries, inlet boundaries, and outlet boundaries. The solid wall boundaries adopt rebound boundary conditions. Combined with LBGK evolution equation formula 1 and local equilibrium distribution function formula 3, the migration and collision of the fluid distribution function are handled to avoid fluid momentum leakage at the boundaries.

[0068] Furthermore, the inlet and outlet boundaries are set with the macroscopic velocity and density of the fluid based on the injection molding process parameters to ensure that the boundary conditions are consistent with the actual working conditions.

[0069] Furthermore, the fluid-structure interaction force from step S3, the particle motion data from step S4, the local porosity from step S5, and the calculation results of the modified LBM equation are used to monitor the fluid shear rate, local porosity, and particle packing density at each time step in real time. Monitoring thresholds are set, and any threshold reached is marked as a region of sudden local rheological change. If the porosity remains above 0.25 for multiple time periods and the particle velocity approaches zero, it is marked as a potential particle blockage location.

[0070] After the simulation terminates, the flow field data from step S2, the final particle distribution from step S4, the density gradient prediction results from step S6, and the abnormal area records from step S7 are integrated and categorized according to injection molding quality evaluation indicators. This generates flow field cloud maps, particle density distribution maps, density gradient maps, and abnormal area annotation maps, which are then output. This differentiated processing of different types of boundaries ensures that boundary conditions closely match actual injection molding conditions. Real-time anomaly capture can be completed synchronously during the simulation, and the setting of monitoring thresholds aligns with the characteristics of ceramic injection molding processes, avoiding misjudgments or omissions.

[0071] In one embodiment of the present invention, in step S5, calculating the local porosity of each spatial location within the simulation area based on the updated positions of all solid particles in step S4 specifically includes the following method: S51: Traverse all fluid calculation grid points within the simulation area, extract the spatial location and particle radius of all solid particles after the update in step S4, set a local influence range for each fluid calculation grid point, determine whether there are solid particles within the range and mark them.

[0072] In an exemplary embodiment, all preset fluid calculation grid points within the simulation area are scanned one by one, and the three-dimensional coordinates of each grid point are recorded. At the same time, the latest spatial position coordinates of all ceramic powder solid particles and the radius parameter of each particle are extracted from the output of step S4.

[0073] Furthermore, for each fluid calculation grid point, combined with the simulated grid step size, its local influence range is set as a spherical region with the grid center as the center and a radius equal to 1.5 times the grid step size. By calculating the straight-line distance from the center of each solid particle to the center of the grid point, and comparing it with the sum of the local influence range radius and the particle radius, it is determined whether the solid particle is within the local influence range of the grid point. All solid particles that meet the conditions are marked, and their position and radius information are stored.

[0074] Traversing all fluid calculation grid points here ensures no omissions within the simulation area. The reasonable setting of the local influence range matches the local characteristics of fluid flow. Geometric distance judgment can accurately screen relevant particles, eliminate particles that do not interact with the fluid at that grid point, and reduce unnecessary calculations.

[0075] S52: For each fluid calculation grid point, calculate the total volume of all marked solid particles within its local influence range, and then calculate the total volume of that local influence range to obtain the solid volume fraction of that grid point.

[0076] In an exemplary embodiment, for each fluid computation grid point where particles have been marked, the volume of each particle is calculated using the sphere volume formula based on the radius parameters of the marked solid particles. The volumes of all solid particles within the local influence range of that grid point are summed to obtain the total volume of solid particles.

[0077] Furthermore, based on the preset local influence range parameters, the total volume of the spherical influence region is calculated. The sum of the solid particle volumes obtained by accumulation is divided by the total volume of the local influence range to obtain the volume ratio of solid particles at the fluid calculation grid point, i.e., the solid volume fraction. This value is recorded and associated with the corresponding grid point coordinates.

[0078] This embodiment uses the sphere volume formula to calculate the particle volume, which closely matches the actual shape of ceramic powder particles. The calculation logic of the volume summation ratio is direct and has no additional empirical assumptions. The obtained solid volume fraction accurately reflects the local particle distribution.

[0079] S53: Substitute the solid volume fraction obtained in step S52 into formula (10), calculate the local porosity of each fluid calculation grid point by subtracting the solid volume fraction from 1, and perform boundary smoothing on the calculation results to remove outliers.

[0080] In an exemplary embodiment, for each fluid calculation grid point, the solid volume fraction of that grid point calculated in step S52 is retrieved and directly substituted into formula (10). According to the calculation rules defined in the formula, the solid volume fraction is subtracted from 1 to obtain the local porosity value corresponding to that grid point.

[0081] Furthermore, after completing the porosity calculation for all fluid calculation grid points, the calculation results are smoothed using the neighborhood averaging method. The porosity values ​​of the eight neighboring grid points around each grid point are selected, the average value is calculated, and the abnormal porosity of that grid point is replaced. Finally, continuous local porosity distribution data at each spatial location within the simulation area are obtained.

[0082] In this embodiment, the calculation is performed by directly substituting into formula (10), strictly following the core formula of the initial technical solution. The calculation results are consistent and traceable. Neighborhood smoothing can eliminate outliers that occur during the calculation process. The obtained continuous porosity data can be directly used to correct the LBM evolution equation in subsequent steps without additional data processing.

[0083] In one embodiment of the present invention, in step S5, based on the local porosity, the evolution equation of the lattice Boltzmann method described in step S2 is modified to calculate the fluid flow under high solids content conditions, specifically including the following methods: S511: Retrieve the local porosity data of each fluid calculation grid point in the simulation area obtained in step S5, and simultaneously extract the fluid distribution function, relaxation time, discrete velocity vector, and current fluid density and velocity of the LBM model in step S2.

[0084] In an exemplary embodiment, the coordinates of each fluid calculation grid point within the simulation area are matched one by one, and the smoothed local porosity values ​​from step S5 are retrieved to establish the correspondence between porosity and grid point coordinates and LBM model parameters. From the simulation cache data in step S2, the fluid distribution function of the current time step, the relaxation time parameter obtained by formula (2), the preset discrete velocity vector, and the current fluid macroscopic density and macroscopic velocity field data obtained by formulas (4) and (5) are extracted. All parameters are classified and archived according to grid point coordinates to ensure that the porosity parameters of each grid point correspond one-to-one with the core parameters of LBM. Invalid data with mismatched coordinates are removed, and a complete set of parameters is retained for subsequent correction calculations.

[0085] Archiving parameters by grid point coordinates ensures targeted corrections, and synchronously retrieving existing LBM parameters eliminates the need for re-initialization, reducing computational redundancy. Eliminating invalid data avoids correction errors caused by parameter mismatches, while retaining a complete set of parameters ensures the continuity of subsequent correction calculations, aligning with the localized differences in fluid flow under high solids content.

[0086] S512: Based on the retrieved local porosity, combined with the porosity data in formula (10) and the Ergun drag coefficient model in formula (13), the drag coefficient correction value under high solid content conditions is calculated, and then substituted into formula (11) to determine the porosity additional force term in the corrected LBGK evolution equation.

[0087] Step S512 specifically includes the following steps: S5121: Batch retrieve the local porosity of each fluid calculation grid point calculated in step S5, substitute it into formula (10) to verify the validity of the data, divide it into three intervals according to the porosity value: high solid content dense phase, transition phase, and dilute phase, and set the corresponding drag coefficient correction weight for each interval.

[0088] In an exemplary embodiment, the local porosity values ​​of all fluid calculation grid points after boundary smoothing in step S5 are retrieved in batches. The porosity of each grid point is substituted into formula (10) one by one to check whether the value is within a reasonable range of 0-1. If the porosity is less than 0, it is corrected to 0; if it is greater than 1, it is corrected to 1 to ensure that the data conforms to the physical definition of formula (10).

[0089] Furthermore, based on the high solids content process characteristics of ceramic injection molding, the porosity is divided into three intervals: porosity ≤ 0.3 is the high solids content dense phase region, 0.3 < porosity < 0.5 is the transition phase region, and porosity ≥ 0.5 is the dilute phase region. Drag coefficient correction weights are assigned to each of the three intervals: the weight for the dense phase region is set to 1.0, the weight for the transition phase region is linearly interpolated based on porosity, and the weight for the dilute phase region is set to 0.5. A correction weight mapping table is established by mapping the porosity intervals, correction weights, and fluid grid coordinates one-to-one, and the calibrated porosity data and weight parameters are stored.

[0090] In this embodiment, local porosity is defined as the fraction of space that the fluid can occupy, with a value in the range of 0-1. The verification operation is to ensure that the data conforms to physical laws and avoid invalid data leading to subsequent calculation deviations. The core feature of ceramic injection molding is high solid content. The particle packing density and fluid resistance differences in different porosity ranges are graded according to porosity and a differentiated correction weight is set. This is based on the variation law of fluid resistance under high solid content, so that the drag coefficient correction can accurately match the particle distribution state in different regions, which is consistent with the logic of formula (13) Ergun model adapting to high-concentration particle beds.

[0091] S5122: Substitute the effective porosity of each grid point into the formula (13) Ergun drag coefficient model to calculate the basic drag coefficient, and then multiply it by the correction weight of the corresponding interval to obtain the drag coefficient correction value under high solid content conditions. Remove abnormal data of correction value and supplement interpolation.

[0092] In an exemplary embodiment, for each fluid calculation grid point, the effective porosity value after calibration in step S5121 is retrieved and substituted into the Ergun drag coefficient model of formula (13) to calculate the basic drag coefficient of the grid point under the current porosity, and to determine the contribution ratio of the two items in formula (13) to the viscous force-dominated drag and the inertial force-dominated drag, respectively.

[0093] Furthermore, based on the porosity range of the grid point, the corresponding correction weight is retrieved from the correction weight mapping table. The basic drag coefficient is multiplied by the correction weight to obtain the corrected drag coefficient value for the grid point under high solids content conditions. The correction values ​​of all grid points are iterated, and a correction value threshold is set. Correction values ​​exceeding the threshold are judged as abnormal data. The average of the correction values ​​of four adjacent grid points is used to obtain interpolated data, which replaces the abnormal values, ensuring that the drag coefficient correction values ​​of all fluid grid points are continuous and reasonable, and finally forming a complete drag coefficient correction value matrix.

[0094] S5123: Substitute the drag coefficient correction value and local porosity data of each grid point into formula (11), and combine the discrete velocity vector retrieved in step S511 to calculate the porosity additional force term for each fluid grid point and each discrete velocity direction, and complete the amplitude and direction calibration of the additional force term.

[0095] In an exemplary embodiment, the discrete velocity vector extracted in step S511 is retrieved, and combined with the drag coefficient correction value of each grid point obtained in step S5122 and the local porosity data calibrated in step S5121, it is substituted into formula (11) to calculate the porosity additional force term amplitude value corresponding to each discrete velocity direction for each fluid calculation grid point.

[0096] Furthermore, based on the direction of the discrete velocity vector, the direction of the additional force term is calibrated to ensure that the direction of the additional force term is opposite to the direction of fluid flow, thus conforming to the physical characteristics of fluid resistance. The additional force terms for all discrete velocity directions at each grid point are summarized, and the resultant force amplitude of the additional force terms is calculated. If the resultant force amplitude exceeds a preset threshold, the amplitude of the additional force terms for each discrete velocity direction is scaled proportionally to ensure that the additional force terms are adapted to the current macroscopic physical quantities of the fluid. Finally, the calibrated porosity additional force term for each fluid grid point and each discrete velocity direction is output for integration of the evolution equations.

[0097] In this embodiment, the calculation of the additional force term is combined with the discrete velocity vector to ensure compatibility with the core parameters of the LBM model. The direction calibration conforms to the physical characteristics of fluid resistance, and the amplitude calibration can match the current macroscopic state of the fluid. The accurate calculation of each grid point and each discrete velocity direction allows the additional force term to truly reflect the local fluid resistance, providing an accurate and suitable core correction term for the establishment of the modified LBM evolution equation. At the same time, no additional adjustment of model parameters is required, reducing computational complexity.

[0098] S513: Integrate the porosity additional force term obtained in step S512 into the original LBGK evolution equation used in step S2, replacing the original uncorrected evolution equation, to obtain the corrected LBM evolution equation adapted to the high solids content environment, which is used for fluid flow calculation in the next time step.

[0099] Step S513 specifically includes the following steps: S5131: Analyze the original LBGK evolution equation used in step S2, clarify the parameter positions and calculation sequence of the migration term, collision term and original external force term, and calibrate the compatibility of the porosity additional force term obtained in step S512 with the time sequence and scale parameters.

[0100] In an exemplary embodiment, the original LBGK evolution equation corresponding to formula (1) used for fluid evolution calculation in step S2 is retrieved, and the migration term, collision term and original external force term in the equation are decomposed one by one. The mathematical expression, discrete velocity vector, relaxation time and migration term are marked for each term, and the external force term is executed first, collision term is executed later and external force term is superimposed synchronously.

[0101] Furthermore, retrieve the porosity additional force term matrix after calibration in step S512, and combine it with the LBM model lattice step size and time step size parameters extracted in step S511 to calibrate the amplitude scale of the additional force term, ensuring that the dimensions of the additional force term are consistent with the original external force term in formula (1).

[0102] The calculation sequence of the additional force term is matched, the calculation order of the additional force term, the migration term, and the collision term is clarified, the additional force term is superimposed synchronously with the original external force term, and it is substituted after the collision term is calculated and before the migration term is executed. The correspondence between the additional force term in each discrete velocity direction and the corresponding discrete velocity distribution function in formula (1) is distinguished, the additional force term with abnormal adaptation is marked and recalibrated, and the parameters and timing matching table of the equation term after disassembly are retained.

[0103] S5132: Integrate the calibrated porosity-added force term into the original equation, replace the original external force term that did not consider porosity, and adopt a differentiated integration method for the grid points at the mold boundary and the grid points inside the simulation area to establish an integrated parameter mapping table.

[0104] In an exemplary embodiment, based on the timing matching table in step S5131, an integration operation of the porosity additional force term is performed for each fluid calculation grid point: for grid points inside the simulation region, the calibrated porosity additional force term is directly used to replace the original external force term in formula (1) to ensure that the calculation logic of the additional force term, migration term, and collision term is consistent, and the expression of the replaced external force term is recorded.

[0105] Furthermore, after the grid points are determined by the mold cavity boundary, the porosity additional force term is attenuated by 0.8-0.9 times the amplitude and then replaced with the original external force term. The attenuation coefficient is set based on the grid density of the mold boundary. The denser the grid, the smaller the attenuation coefficient, in order to adapt to the constraint characteristics of fluid flow at the boundary.

[0106] After integration, an integrated parameter mapping table is established for each fluid calculation grid point, recording the grid point coordinates, the magnitude of the additional force term, the attenuation coefficient, and the expression of the replaced external force term. The integrated equation terms are verified one by one, and problems such as expression errors and missing parameters caused by integration errors are eliminated to ensure that the integration operation of each grid point is accurate and error-free, forming the corrected LBGK evolution equation.

[0107] Furthermore, the fluid flow states of the internal grid points and the boundary grid points differ in the simulation region. The fluid flow in the internal grid points is not constrained by the mold, which can reflect the hindering effect of the porosity-added force term. The boundary grid points are constrained by the solid wall of the mold, and the fluid flow resistance comes not only from particle accumulation but also from the friction of the mold wall. If the added force term is not attenuated, the fluid resistance at the boundary will be over-calculated, deviating from the actual working condition. The relationship between the attenuation coefficient and the grid density is based on the interaction strength between the fluid at the boundary grid points and the mold wall. The denser the grid, the stronger the boundary constraint, and the smaller the attenuation of the added force term, which fits the core logic of formula (1) for handling complex boundary conditions.

[0108] S5133: Verify the mathematical integrity and conservation of the integrated equations, update the LBM equation library of the simulation system, synchronize the calculation parameters of the next time step of the associated step S2, and ensure that the corrected equations can be directly called.

[0109] In an exemplary embodiment, the mathematical integrity and physical conservation of the preliminary modified LBGK evolution equation obtained in step S5132 are verified: In terms of mathematical integrity, the expressions of the migration term, collision term, and modified external force term in the equation are checked to see if they are complete and if any parameters are missing, ensuring that the equation has no syntax errors. 10-15 fluid calculation grid points are randomly selected, and the grid point parameters are substituted into the modified equation to calculate the fluid distribution function. Then, the function is substituted into formulas (4) and (5) to statistically calculate the macroscopic density and velocity field, verifying the conservation of mass and momentum. If the conservation requirements are not met, the process returns to step S5131 to recalibrate the amplitude and timing of the additional force term. After verification, the modified LBM evolution equation is updated to the equation library of the simulation system, overwriting the original unmodified formula (1). The calculation parameters of the next time step in step S2 are synchronously associated, the integrated version of the equation and the corresponding time step are marked, and an equation update log is generated to ensure that the fluid flow calculation in the next time step can directly call the modified equation. At the same time, the verification data and update log are retained for easy anomaly investigation.

[0110] In one embodiment of the present invention, S6 specifically includes the following: S61. Extract the three-dimensional coordinates of all solid particles from the final particle position data output in step S4. Based on the geometry of the mold cavity, establish a three-dimensional analysis grid that is similar to or slightly larger than the original LBM fluid grid scale to cover the entire cavity. Assign each particle to its corresponding analysis grid cell according to its coordinates, and accumulate the number of particles in each cell.

[0111] In one exemplary embodiment, after the simulation is completed, the discrete element solver outputs the coordinate position (x, y, z) of each particle at the final time step.

[0112] To perform spatial statistics, the continuous space is discretized. An analysis grid independent of the fluid computation grid is established. The choice of the analysis grid size is a trade-off: too small a grid size will result in many cells without particles and high statistical noise; too large a grid will smooth out meaningful local variations.

[0113] Furthermore, the grid size is set to 3 to 5 times the average particle diameter. This ensures that each cell contains a number of particles to obtain statistical significance while maintaining a certain spatial resolution. All particles are traversed, and their corresponding analysis grid cell index is quickly determined based on their coordinates. The particle count or volume accumulation is then added to that cell. This process is similar to the volume projection in step S51, but the purpose is different. Here, the focus is more on counting and simple volume accumulation to prepare for calculating the solid fraction.

[0114] This embodiment divides the continuous cavity space into discrete cubic units, classifying the continuously distributed particles into these units. This is similar to creating a three-dimensional histogram, where the height of each histogram bar represents the number of particles or the total volume falling into that spatial region.

[0115] This embodiment transforms thousands of discrete particle coordinates into a regular three-dimensional array through meshing, which is suitable for calculating the mean, standard deviation, gradient, and visualization.

[0116] S62. For each analysis grid cell, calculate the solid volume fraction φlocal based on the number of particles contained within it and the average volume of a single particle.

[0117] Specifically, the solid volume fraction of the analytical grid cell is the sum of the volumes of all particles within the cell divided by the volume of that analytical grid cell. Simultaneously, using formula (10), the porosity of the analytical grid cell is defined as εlocal = 1 - φlocal. By traversing all analytical grid cells, the spatial distribution field of the solid phase fraction φ(x) and its corresponding porosity field ε(x) within the entire cavity are generated.

[0118] In one exemplary embodiment, for each analysis grid cell, the program reads its cumulative total particle volume. If the particle size distribution in the simulation is narrow, the total volume can be approximated by multiplying the average particle volume by the number of particles.

[0119] If the particle size distribution is wide, the actual volume of each particle is used during accumulation in step S61. Dividing this total volume by the volume of the grid cell yields the local solid volume fraction φlocal for the analytical grid cell. This value reflects the degree of filling of the analytical grid cell space. According to the definition of porosity, the local porosity εlocal = 1 - φlocal.

[0120] S63. Based on the solid fractional spatial distribution field φ(x) generated in step S62, calculate the spatial gradient ∇φ.

[0121] In three-dimensional space, the magnitude of the gradient |∇φ| characterizes the severity of local density changes. The entire cavity is divided into several macroscopic regions, and the mean and standard deviation of φ(x) within each macroscopic region are calculated. By comparing the differences in mean values ​​between different macroscopic regions and the standard deviations within each region, the non-uniformity of green body density at both macroscopic and microscopic scales is quantitatively assessed, and a density gradient distribution map is generated.

[0122] In an exemplary embodiment, the solid fractional field φ(x) is numerically differentiated, such as by using the central difference formula to calculate the gradient components of each grid cell in the x, y, and z directions, thereby obtaining the gradient vector ∇φ and its magnitude |∇φ|.

[0123] The presence of a large gradient modulus indicates a drastic change in density between adjacent regions, corresponding to potential defect interfaces such as density abrupt change zones. To assess overall and local homogeneity, a larger analysis region is defined. For example, along the path of the melt flowing from the gate to the end, the cavity is divided into 10 equal segments, and the average value of φ for all mesh elements in each segment is calculated.

[0124] If these average values ​​show a monotonically increasing or decreasing trend from the gate to the end, it indicates the presence of a macroscopic density gradient. Calculating the standard deviation of φ within each segment allows for the assessment of the microscopic uniformity within that segment. Areas with large standard deviations indicate significant density fluctuations even within small ranges, potentially indicating poor green body quality. Finally, the |∇φ| field, segmented average values, and standard deviations are output as contour plots, line graphs, or bar charts, and can be compared with the maximum allowable density variation rate to provide a conclusive judgment on whether the area is acceptable or at risk.

[0125] It should be understood that the sequence number of each step in the above embodiments does not imply the order of execution. The execution order of each process should be determined by its function and internal logic, and should not constitute any limitation on the implementation process of the embodiments of the present invention.

[0126] like Figure 3 As shown, this application also provides an electronic device, including a memory 402, a processor 401, an input device 403, an output device 404, and a computer program stored in the memory and executable on the processor 401. When the processor 401 executes the program, it implements the steps of a multiphase flow calculation method for ceramic injection molding based on the lattice Boltzmann and discrete element method.

[0127] In this embodiment of the invention, the input device 403 may be a mouse and keyboard, or a voice input device. The output device 404 may be a display screen, a display panel, or a touch screen. The processor 101 may be implemented using at least one of an application-specific integrated circuit, a microcontroller, a microprocessor, or an electronic unit designed to perform the functions described herein; in some cases, such an implementation may be implemented in a controller. For software implementations, implementations such as processes or functions may be implemented with separate software modules that allow the performance of at least one function or operation. The software code may be implemented by a software application (or program) written in any suitable programming language, and the software code may be stored in memory and executed by the controller.

[0128] The present invention also provides a storage medium storing a computer program thereon, wherein the computer program, when executed by a processor, implements the steps of the ceramic injection molding multiphase flow calculation method based on the lattice Boltzmann and discrete element method.

[0129] The above description of the disclosed embodiments enables those skilled in the art to make or use the invention. Various modifications to these embodiments will be readily apparent to those skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of the invention. Therefore, the invention is not to be limited to the embodiments shown herein, but is to be accorded the widest scope consistent with the principles and novel features disclosed herein.

Claims

1. A method for calculating multiphase flow in ceramic injection molding based on lattice Boltzmann and discrete element method, characterized in that the method... include: S1: Establish a fluid-particle bidirectional coupling model to simulate the flow of ceramic feed in the mold; S2: Based on the fluid-particle bidirectional coupling model, input the initial ceramic powder particle parameters, binder fluid parameters and mold cavity boundary conditions, and perform transient simulation calculations on the injection molding filling process of ceramic feeding to obtain the evolution data of particle position, velocity and fluid flow state. S3: During the simulation calculation, for each time step, based on the current fluid flow state obtained in step S2, the fluid-structure interaction force acting on each solid particle is calculated using the momentum exchange method. S4: Based on the fluid-structure interaction force, interparticle contact force and gravity calculated in step S3, and using Newton's second law, calculate the translational velocity and rotational velocity of each solid particle in the next time step, and update its position and orientation. S5: Based on the updated positions of all solid particles in step S4, calculate the local porosity at each spatial location within the simulation area; and based on the local porosity, correct the evolution equation of the lattice Boltzmann method in step S2 to calculate fluid flow under high solid content conditions. S6: Based on the simulation results of steps S2 to S5, the final spatial distribution of solid particles in the mold cavity when the injection molding process is completed is statistically analyzed, and the density gradient distribution inside the green blank is predicted based on the uniformity of the particle distribution. S7: During the simulation process in step S2, the geometric boundary conditions of the micro or complex mold shape are processed by the lattice Boltzmann method. Combined with the force and motion calculations in steps S3 to S5, local rheological abrupt changes and potential particle blockage locations are captured, and finally the flow field, particle distribution and defect prediction results are output for evaluating the injection molding quality.

2. The method for calculating multiphase flow in ceramic injection molding based on lattice Boltzmann and discrete element method according to claim 1, characterized in that, S1 specifically includes the following methods: S11. Establish and initialize the lattice Boltzmann solver and the discrete element method solver, configure the fluid and particle parameters, and establish a data exchange interface for the two solvers. S12. Identify fluid-structure interaction interfaces in the computational grid, calculate the force exerted by the fluid on the particles based on the momentum exchange method, and correct the force in the high solid content region. S13. The fluid force is transmitted to the discrete element method solver to update the particle motion state, and the updated particle information is fed back to the lattice Boltzmann solver, closing the bidirectional coupling loop.

3. The method for calculating multiphase flow in ceramic injection molding based on lattice Boltzmann and discrete element method according to claim 1, characterized in that, S3 specifically includes the following methods: S31: Identify the fluid computation grid points associated with the surface of each solid particle and obtain the fluid distribution function data at the grid points; S32: Based on the acquired fluid distribution function data, the momentum change of each fluid grid point is calculated by the momentum exchange method, and the initial fluid-structure interaction force acting on the particles is obtained by summing them up. S33: Based on the drag coefficient calculated from the local porosity, the initial fluid-structure interaction force is corrected to obtain the final fluid-structure interaction force.

4. The method for calculating multiphase flow in ceramic injection molding based on lattice Boltzmann and discrete element method according to claim 1, characterized in that, In S5, the local porosity at each spatial location within the simulation region is calculated based on the updated positions of all solid particles from step S4, specifically in the following ways: S51: Traverse all fluid calculation grid points within the simulation area, extract the spatial location and particle radius of all solid particles after the update in step S4, set a local influence range for each fluid calculation grid point, determine whether there are solid particles within the range and mark them. S52: For each fluid calculation grid point, calculate the total volume of all marked solid particles within its local influence range, and then calculate the total volume of the local influence range to obtain the solid volume fraction of that grid point; S53: Calculate the solid volume fraction obtained in step S52, and then calculate the local porosity of each fluid calculation grid point by subtracting the solid volume fraction from 1. Perform boundary smoothing on the calculation results to remove outliers.

5. The method for calculating multiphase flow in ceramic injection molding based on lattice Boltzmann and discrete element method according to claim 4, characterized in that, In S5, based on the local porosity, the evolution equation of the lattice Boltzmann method described in step S2 is modified to calculate the fluid flow under high solids content conditions. This specifically includes the following methods: S511: Retrieve the local porosity data of each fluid calculation grid point in the simulation area obtained in step S5, and extract the fluid distribution function, relaxation time, discrete velocity vector, and current fluid density and velocity from the LBM model in step S2; S512: Based on the retrieved local porosity, combined with porosity data and Ergun drag coefficient model, calculate the corrected value of drag coefficient under high solid content conditions, and then determine the porosity additional force term in the corrected LBGK evolution equation. S513: Integrate the porosity additional force term obtained in step S512 into the original LBGK evolution equation used in step S2, replacing the original uncorrected evolution equation, to obtain the corrected LBM evolution equation adapted to the high solids content environment, which is used for fluid flow calculation in the next time step.

6. The method for calculating multiphase flow in ceramic injection molding based on lattice Boltzmann and discrete element method according to claim 5, characterized in that, SS512 specifically includes the following methods: S5121: Batch retrieve the local porosity of each fluid calculation grid point calculated in step S5, verify the validity of the data, divide the three intervals of high solid content dense phase, transition phase and dilute phase according to the porosity value, and set the corresponding drag coefficient correction weight for each interval. S5122: Substitute the effective porosity of each grid point into the Ergun drag coefficient model to calculate the basic drag coefficient, and then multiply it by the correction weight of the corresponding interval to obtain the drag coefficient correction value under high solid content conditions. Remove abnormal data of correction value and supplement interpolation. S5123: Combine the drag coefficient correction value of each grid point and the local porosity data with the discrete velocity vector retrieved in step S511 to calculate the porosity additional force term for each fluid grid point and each discrete velocity direction, and complete the amplitude and direction calibration of the additional force term.

7. The method for calculating multiphase flow in ceramic injection molding based on lattice Boltzmann and discrete element method according to claim 5, characterized in that, S513 specifically includes the following methods: S5131: Analyze the original LBGK evolution equation used in step S2, clarify the parameter positions and calculation sequence of the migration term, collision term and original external force term, and calibrate the compatibility of the porosity additional force term obtained in step S512 with the time sequence and scale parameters. S5132: Integrate the calibrated porosity-added force term into the original equation, replace the original external force term that did not consider porosity, and adopt a differentiated integration method for the grid points at the mold boundary and the grid points inside the simulation area to establish an integrated parameter mapping table. S5133: Verify the mathematical integrity and conservation of the integrated equations, update the LBM equation library of the simulation system, synchronize the calculation parameters of the next time step of the associated step S2, and ensure that the corrected equations can be directly called.

8. The method for calculating multiphase flow in ceramic injection molding based on lattice Boltzmann and discrete element method according to claim 1, characterized in that, S6 specifically includes the following methods: S61. Assign the three-dimensional coordinates of the solid particles at the final moment to the analysis mesh cells covering the mold cavity, and count the number of particles in each cell. S62. Based on the number and volume of particles in each analytical grid cell, calculate the solid volume fraction of the analytical grid cell and convert it into porosity to form a spatial distribution field of solid phase fraction and porosity. S63. Calculate the spatial gradient based on the solid fractional distribution field, divide the macroscopic region and statistically analyze the mean and standard deviation of each region to quantify the non-uniformity of the density distribution.

9. An electronic device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that, When the processor executes the program, it implements the steps of the ceramic injection molding multiphase flow calculation method based on the lattice Boltzmann and discrete element method as described in any one of claims 1 to 8.

10. A storage medium having a computer program stored thereon, characterized in that, When the computer program is executed by the processor, it implements the steps of the ceramic injection molding multiphase flow calculation method based on the lattice Boltzmann and discrete element method as described in any one of claims 1 to 8.

Images