A method for calculating implicit heat transfer of moving particles using a multi-time-layer difference scheme

Abstract

A multi-time-layer differential scheme for implicit heat transfer calculation using moving particles is presented, with the following specific steps: 1. Before calculation, virtual particles are arranged at the boundary to implement boundary layout, using virtual particles to arrange the first, second, and third types of boundaries for heat transfer calculation, as well as symmetric and periodic geometric boundary conditions; 2. During heat transfer calculation, the physical and mathematical relationships between different particles are given according to the input boundary conditions; 3. The distances between particle i and its neighboring particle j from the boundary surface formed by the boundary particles are calculated, and boundary conditions are applied; 4. The particle properties are given according to the initial input temperature, and the calculation time step length is given according to the stability condition; 5. The coefficient matrix and source terms of the multi-time-layer scheme of the heat diffusion equations are combined to solve the energy equation for temperature; 6. A consistent solution is performed using relevant solving algorithms to solve for the next temperature and determine whether convergence has occurred; 7. After all MPS calculations for a single time step are completed, the next time step is initiated.

Description

Technical Field

[0001] This invention relates to the field of computational technology of the moving semi-implicit particle method, specifically to a moving particle implicit heat transfer calculation method using a multi-time-layer difference scheme. Background Technology

[0002] The moving semi-implicit particle method (MSP) is a Lagrangian meshless method developed from smooth particle hydrodynamics. Compared to traditional mesh methods, it avoids the inherent defects of mesh distortion and eliminates dissipation terms in the momentum equation based on the Lagrangian discretization scheme, thus improving the accuracy of numerical calculations. Due to its Lagrangian discretization method, it has unique advantages in capturing free surface and interface motion, phase changes, and solid deformation. It has significant advantages in dealing with multi-component, multi-phase problems of phase change melting, bubble dynamics, fluid-structure interaction, droplet behavior, and boiling behavior, which are of great interest to traditional mesh methods. Currently, the MSP is widely used internationally in severe reactor core meltdown accidents to study phenomena of great concern in the field, such as the generation, migration, solidification, and repositioning of molten material, the interaction between molten material and concrete, zirconium-water reaction, and the interaction between molten material and coolant. For any numerical calculation method, heat transfer calculation is an indispensable part. However, the current semi-implicit moving particle method has problems with accurately applying the heat transfer boundary, and it also has problems with small particles generating small time steps and large-scale cases taking a long time to compute. These issues need to be addressed through other measures. Summary of the Invention

[0003] To further improve the computational speed and accuracy of moving semi-implicit particle method simulation analysis, this invention aims to provide a multi-time-layer difference scheme for calculating implicit heat transfer in moving particles. This invention uses an implicit solution method, eliminating the stringent time step limitations of the explicit portion in the original semi-implicit method. It also enhances the original semi-implicit method's ability to handle situations where the time step is too small due to stability strategy limitations caused by small particles arranged in a small computational domain, special physical properties, and other factors. Simultaneously, it meets the requirement for more accurate results by using a multi-time-layer difference scheme to solve the implicit equations, obtaining more accurate results than the single-time-layer scheme.

[0004] To achieve the above objectives, the present invention adopts the following technical solution:

[0005] A method for calculating implicit heat transfer of moving particles using a multi-time-layer difference scheme, comprising the following steps:

[0006] Step 1: Arrangement of the computational domain and its boundaries. Before computation, arrange the particles of the computational domain according to its extent, and then arrange the boundary-related virtual particles according to each boundary of the computational domain. This includes three types of heat transfer boundary conditions: the first type of temperature boundary condition, the second type of heat flux boundary condition, and the third type of convective heat transfer boundary condition; and two types of geometric boundary conditions used to simplify the computational domain: symmetric boundaries and periodic boundaries. The specific functions of each boundary are as follows:

[0007] The first type of temperature boundary condition is applicable to heat transfer calculations under isothermal walls; the second type of heat flux boundary condition is applicable to surfaces with constant power input heat, among which the adiabatic boundary is a special type of heat flux boundary condition; the third type of convective heat transfer boundary condition is applicable to convective heat transfer between the heat transfer surface and the external medium; geometric boundary conditions are used for symmetrically arranged computational domains and for periodic boundaries used in infinitely large or periodic computational domains. By using geometric boundary conditions, only a portion of the computational domain needs to be arranged, thereby simplifying calculations, reducing computational costs, and accelerating calculation speed.

[0008] In particle arrangement, all three heat transfer boundary arrangements require three layers of virtual particles to be arranged outside the computational domain to satisfy the support domain. The additional virtual particles not only satisfy the support domain of the particles in the computational domain, but also serve as supplementary nodes for the realization of heat transfer boundary conditions.

[0009] The two geometric boundary conditions supplement the support domain by mapping based on geometric relationships. Therefore, after defining the geometric boundary, the distance from the boundary is mapped into the outside of the boundary as supplementary particles. These mappings are one-to-one, which means that the position parameters are determined, so there is no need to arrange virtual particles. The geometric boundary needs to be arranged according to the overall simulation content and the requirements of simplifying the computational domain. The arrangement of periodic boundaries requires declaring that the two boundaries are periodic boundaries. Symmetrical boundaries cannot be arranged on curved surfaces.

[0010] Step 2: Construct the relationship between virtual particles and internal particles based on the boundary layout;

[0011] During the calculation, the temperature values ​​of the virtual particles are not directly added to the equations formed by the thermal equations as source terms for solving. Instead, relationships between the particles are constructed through different selected boundaries, and all particles are listed in the thermal equations. The resulting large matrix for solving is then combined with the relationships between the particles at the boundaries. Specifically, particles are added near the boundaries: for heat transfer boundaries, virtual particles are used; for geometric boundaries, mappings generated by geometric relationships are used. These added particles are then replaced with equations using selected boundary conditions. There are five relationships in total for three types of heat transfer boundaries and two types of geometric boundaries, as follows:

[0012] For the first type of temperature boundary condition, the temperature relationship between particles inside the computational domain and virtual particles in neighboring domains, as well as the distances of both types of particles from the boundary, are as follows:

[0013]

[0014] In the formula

[0015] i — an internal particle of the computational domain

[0016] j — Particles within the neighborhood of the internal particle. dummy — Refers to the set of virtual particles that supplement the boundary under the following boundary conditions: temperature boundary condition, heat flux boundary condition, convective heat transfer boundary condition, symmetric boundary condition, or periodic boundary condition.

[0017] j,dummy — Virtual particle j within the neighborhood of internal particle i.

[0018] T Dummy —Temperature at which virtual particles are used to fill the gaps where particles are missing at the boundary

[0019] T Boundary —Given boundary temperature

[0020] T Internal —Internal particle temperature of the computational domain

[0021] d j,dummy — The distance of virtual particle j within the neighborhood of internal particle i from the boundary.

[0022] d i — Distance of internal particle i from the boundary

[0023] For the second type of heat flux boundary condition and the third type of convective heat transfer boundary condition, the additional node method is used to perform Taylor expansion at the boundary, ensuring that virtual particles outside the boundary and particles inside the computational domain satisfy the applied boundary conditions. The particles inside the computational domain and virtual particles in neighboring domains have the following relationship:

[0024]

[0025]

[0026] In the formula

[0027] O(h 2 — Taylor expansion remainder

[0028] —Temperature gradient relationship at the boundary, For the second type of boundary For the setting of convective heat transfer boundaries, i.e., there is a relationship

[0029]

[0030] In the formula

[0031] T ∞ Ambient temperature

[0032] h—convective heat transfer coefficient

[0033] k—thermal conductivity of the material

[0034] When calculating the convection boundary, the physical property parameters required for the boundary arrangement need to be input in advance. For constant physical properties, the boundary arrangement can be directly solved by solving the virtual particle temperature outside the boundary. For polynomials with given physical properties and temperatures, algorithms for solving nonlinear problems are used for the solution.

[0035] Under symmetric and periodic boundaries, there exists a one-to-one mapping between particles at the boundary. The relationship between particles with one-to-one mapping at the geometric boundary is as follows:

[0036] T Symmetric,j′ =T Internal,j Formula (5)

[0037] T Periodic,j′ =T Internal,j Formula (6)

[0038] In the formula

[0039] T Internal,j —The temperature of particle j within the neighborhood of particle i, and particle j is within the computational domain.

[0040] T Symmetric,j′ —The temperature of particle j′ within the neighborhood of particle i, but particle j′ is not within the computational domain, and there is a mapping relationship between it and particle j due to the existence of the symmetry boundary.

[0041] T Periodic,j′ —The temperature of particle j′ within the neighborhood of particle i, but particle j′ is not within the computational domain, and there is a mapping relationship between the periodic boundary and particle j.

[0042] Step 3: Calculate the distance between the particle and the boundary;

[0043] Formulas (1), (2), (3) and (4) show that the application of the convective heat transfer boundary uses the Taylor expansion method based on the interparticle distance, so it is necessary to calculate the distance from the particles in the computational domain to the heat transfer boundary.

[0044] At a given convective heat transfer boundary, boundary nodes are generated based on particles within the computational domain. Boundary node curves and surface equations are then fitted using numerical methods. The generation principle for boundary nodes is to use the first layer of particles at the boundary, extrapolating the normal vector of each particle by half a particle radius to obtain a boundary node. After acquiring the position information of all boundary nodes, a polynomial curve or surface is fitted. In the two-dimensional simulation, the fitted boundary node curve is represented as:

[0045] z(x) = A + Bx + Cx 2 Formula (7)

[0046] In the 3D simulation, the fitted surface for the boundary nodes is written as:

[0047] z(x,y)=A+Bx+Cy+Dx 2 +Ey 2 +Fxy formula(8)

[0048] In the formula:

[0049] z(x) — the fitted curve function

[0050] z(x, y) — the fitted surface function

[0051] A, B, C, D, E, and F are all coefficients of the fitted curve and surface.

[0052] In particle computation near the boundary, given a heat transfer boundary, it is necessary to obtain the particle distance to the boundary, that is, to solve the particle distance to the boundary, which is transformed into solving the x and y coordinate distance fitting boundary curve, i.e., a two-dimensional or surface, i.e., a three-dimensional distance problem.

[0053] Step 4: Update the physical properties based on the input initial temperature or the temperature of the previous time step; the time step length for a single calculation process is:

[0054]

[0055] In the formula:

[0056] Δt — the length of the time step in a single computation process

[0057] S—Dissipation number, a constant that controls the calculation to prevent divergence, determined according to the selected time difference scheme; l0—Initial particle arrangement distance in the particle method, also known as particle radius.

[0058] α max — The maximum thermal diffusivity of all particles in the computational domain at the previous or initial time step. Thermal diffusivity is the ratio of thermal conductivity to the product of density and specific heat.

[0059] Step 5: After obtaining the distance between the particle and the boundary, accurate heat transfer boundary conditions can be applied. Simultaneously, after obtaining the time step length of a single computational process, the heat transfer problem can be solved. In the particle method, solving the heat transfer problem is equivalent to solving the diffusion equation.

[0060]

[0061] In the formula

[0062] T—Particle temperature within the computational domain

[0063] —Internal heat source term, pre-defined based on specific calculation conditions.

[0064] ρ — density, given as a function of temperature

[0065] c — heat capacity, given as a function of temperature.

[0066] For each internal particle i that needs to be solved within the computational domain, a discrete equation can be written, that is, the governing equation is discretized as:

[0067]

[0068] In the formula:

[0069] d — the number of dimensions in the simulation

[0070] λ — A coefficient that makes the particle interaction model consistent with the analytical solution of the Gaussian function.

[0071] n0 — Initial particle number density

[0072] w(r ij ,r e — The weight function used in the particle method, r ij ,r e The independent variable of the weight function

[0073] r ij —Distance between internal particle i and neighboring particle j

[0074] r e —Given particle size

[0075] — This represents the solution for the temperature of internal particle i, where the superscript n+1 represents the result to be solved at the next time step.

[0076] — This represents the solution for the temperature of internal particle i, where the superscript n represents the known result of the previous time step.

[0077] — This represents solving for the temperature of particle j within the neighborhood of particle i, where the superscript n represents the known result of the previous time step.

[0078] k ij — Thermal conductivity, used in the calculation, is the harmonic mean of the internal particle i and its neighboring particles j.

[0079]

[0080] In the formula:

[0081] k i — Solving for the thermal conductivity of internal particle i

[0082] k j — Solving for the thermal conductivity of particle j in the neighborhood of particle i.

[0083] The difference method using a two-level time scheme implicitly defines the control equations. The multi-time-level difference scheme for the control equations is as follows:

[0084]

[0085] In the formula:

[0086] — Solve for the temperature of particle j within the neighborhood of particle i, where the superscript n+1 represents the result to be solved at the next time step.

[0087] f is a constant that controls the time difference scheme. f = 0 is an explicit difference scheme, f = 0.5 is a Crank-Nicholson difference scheme, and f = 1 is an implicit difference scheme. The difference between 0 and 0.5 is a partially explicit difference scheme, and the difference between 0.5 and 1 is a partially implicit difference scheme.

[0088] In explicit difference schemes, the coefficient S used to control the time step in step 4 needs to be less than 1. In the Crank-Nicholson difference scheme, the coefficient S used to control the time step needs to be less than 2. For implicit difference schemes, there is no limit to the size of the time step.

[0089] In various schemes, except for f=0 (i.e., explicit calculation does not require solving implicit equations), all other difference schemes require solving a system of equations for all particles within the computational domain, and then solving the matrix. When the neighboring particle j of the internal particle i is a supplementary particle j, dummy, based on the temperature relationship between the two particles (i.e., the boundary condition), the boundary particle representing the boundary condition is converted into a particle term within the computational domain. These terms are then placed into the diagonal and source terms of the matrix and solved consistently. The specific operation is as follows: For the first type of boundary, the matrix assembly process is as follows: If the neighboring particle i of the internal particle i has a supplementary particle j, dummy representing the first type of temperature boundary condition, the supplementary particle j, dummy is separated from the particles within the computational domain, i.e.:

[0090]

[0091] Based on formula (14), the relationship between the virtual particle j,dummy in the neighborhood domain and the i particle to be solved under the first type of heat transfer boundary is supplemented and written as:

[0092]

[0093] The operation is the same for the second and third types of heat transfer boundaries, and the boundary condition relations are used to replace the virtual particles j,dummy supplemented in the neighboring domain:

[0094]

[0095] For symmetric and periodic boundary conditions, since there is a mapping relationship between the supplementary particle j′ in the neighborhood and the j particle in the neighborhood of the i particle being solved, the boundary relationship is also substituted into the equation, and the symmetric boundary conditions are:

[0096]

[0097] The periodic boundaries are:

[0098]

[0099] Meanwhile, the geometric boundary is used to supplement the particle j′ term in the neighboring domain, and can be solved by merging it with the particle j term in the neighboring domain that is in the computational domain according to the mapping relationship;

[0100] Step 6: Solving the diffusion equation using a multi-time-layer difference scheme

[0101] After obtaining the difference scheme for heat transfer calculation in step 5, equations 15, 16, 17, and 18 are used. Since an equation can be written for each particle i in the computational domain, the equations formed by all particles in the solution domain are combined to obtain the heat diffusion equation set, which can be simplified to matrix form Ax = b, as follows:

[0102]

[0103] In the formula

[0104] b n —Source term value, where the subscript n represents the source term of the equation for the nth particle, b1 refers to the source term of the equation written for the 1st particle, b a The source term of the equation written for particle a —The temperature value to be solved. The subscript n indicates that the temperature value represents the temperature of particle n, and the superscript n+1 indicates that the temperature value A to be solved in the next time step. mn — The elements in the matrix are formed by the subscript mn, which represents the coefficient between particle m and particle n, that is, the coefficient before the temperature term of particle n in the diffusion equation for particle m; there are three types of coefficients in the coefficient matrix: diagonal elements, elements within neighboring particles, and elements outside neighboring particles.

[0105] The coefficient of the nth diagonal element is A nn For particle number 1, there is A 11 :

[0106]

[0107] For the element 'a' inside the neighboring particle and the element 'n' outside the neighboring particle, we have:

[0108]

[0109] A 1n =0 formula (22)

[0110] The source term is:

[0111]

[0112] In the moving particle method, during the calculation of the smallest unit, i.e., the particle, only the interaction between neighboring particles within the neighborhood is considered. Therefore, in the thermal diffusion equations, this is reflected in the fact that only the neighboring particle terms within the neighborhood are non-zero in each row, while all other terms are zero. Thus, the coefficient matrix of the thermal diffusion equations is sparse. At the same time, the neighborhoods of all particles are of equal size, and they are all neighbors of each other. Therefore, the coefficient matrix of the equations is symmetric. Considering that after substituting the boundary condition relations and simplifying by rearranging terms, only the diagonal terms of the coefficient matrix and the source terms of the thermal diffusion equations change, the solution of the boundary conditions does not cause a change in the symmetry of the thermal diffusion equations.

[0113] Since the coefficient matrix of the heat diffusion equations is a large symmetric sparse matrix, a large sparse symmetric matrix solver and a compatible solution algorithm, such as the conjugate gradient method, the conjugate gradient square algorithm, the biconjugate gradient method, the stable biconjugate gradient method, the conjugate residual method, the symmetric Langios algorithm, or the stable biconjugate residual algorithm, can be used to iteratively solve the equations.

[0114] The convergence criterion for the analytical solution of the thermal diffusion equations is:

[0115] ||r2<||b2ε formula (24)

[0116] In the formula:

[0117] ||r||2——The L2 norm of the residual vector, where the residual vector is calculated as r=b-Ax

[0118] ||b||2——The 2-norm of the source terms in the thermal diffusion equations

[0119] ε—a very small number, used to represent the constant of the convergence criterion for iteration. If the residual under the solution vector x after one iteration satisfies the above convergence condition, then it is considered to have converged at this time step, and the calculation of the next time step is performed.

[0120] Step 7: Advance the time step;

[0121] After a time step calculation is completed and the iterative solution converges, it is necessary to determine whether the calculation termination condition is met. If the calculation termination condition is not met, return to steps 4, 5, 6, and 7 to continue to the next time step until the calculation terminates.

[0122] In summary, steps 1 to 3 input the boundary conditions of the specific problem to be calculated, determine the given parameters, generate particles within the computational domain, and form the boundaries of the computational domain by arranging particles and boundary nodes. During the neighbor particle search process, the distance of each particle at the boundary is determined, and the relationship between the boundary virtual particles and the particles within the computational domain is given. Step 4 updates the initial physical properties or the physical properties calculated in the previous time step, and gives a single time step for subsequent calculations that meet the stability conditions. Step 5 performs differential calculus on the governing equations of heat transfer calculation to form a large sparse symmetric coefficient matrix that can be solved. Step 6 solves the matrix formed in the previous step to obtain the temperature of each smallest computational unit particle in the next time step. The final step 7 determines whether the calculation meets the termination condition. If the termination condition is not met, steps 4, 5, 6, and 7 are returned to continue the next time step until the calculation terminates.

[0123] Compared with existing technologies, the method of the present invention has the following advantages:

[0124] 1) This method optimizes the heat transfer calculation module in current mobile semi-implicit particle method calculations. It uses more accurate heat transfer boundaries and geometry editing. Building upon the original explicit calculation method, it employs an implicit scheme with multiple time-level differences to solve the equations for all particles simultaneously. The explicit and implicit proportions in the calculation process are specified by inputting a constant f, and boundary conditions are consistently added to the coefficient matrix. While addressing the inherent errors in the boundary arrangement of the original mobile particle method, it uses a more accurate boundary arrangement method and significantly increases the time step size of a single calculation implicitly, while still maintaining stability. Therefore, this method enables stable and accurate large-time-step calculations.

[0125] 2) Currently, no researchers have proposed an algorithm that optimizes heat transfer calculations using the moving semi-implicit particle method by incorporating implicit computation with boundary conditions through multi-time-layer difference schemes. This invention proposes an advanced and accurate boundary condition arrangement method (three heat transfer boundaries and two geometric boundaries). This more accurate boundary arrangement avoids errors caused by current boundary condition placement methods. Furthermore, this invention includes most boundary condition forms required for simulation calculations, exhibiting high accuracy and versatility. Employing implicit computation, it allows for larger time steps to advance the calculation process, provided stability conditions permit, resulting in faster computation speed while maintaining computational stability. Therefore, this invention offers better stability and computational speed. Attached Figure Description

[0126] Figure 1 This is a schematic diagram of the calculation process for heat transfer calculation using the multi-time-layer difference scheme of the meshless moving particle method according to the present invention.

[0127] Figure 2 This is a schematic diagram of particles with various boundary arrangements during the initial calculation of this invention. Figure a is a schematic diagram of the heat transfer boundary arrangement, including internal particles, boundary particles, and boundary (virtual) particles; Figure b is a schematic diagram of the geometric boundary arrangement, showing the mapping relationship generated under the geometric relationship between particles.

[0128] Figure 3 This is a schematic diagram of the boundary node arrangement for the boundary in this invention.

[0129] Figure 4 This is an iterative convergence curve of one-dimensional plate heat conduction under different algorithms, demonstrating the calculation results of this invention.

[0130] Figure 5 This is an iterative convergence curve of the thermal conductivity of a square plate under different algorithms, demonstrating the calculation results of this invention.

[0131] Figure 6This is a one-dimensional flat plate heat conduction final temperature distribution cloud map demonstrating the calculation results of this invention.

[0132] Figure 7 This is a cloud map showing the final temperature distribution of the square plate's thermal conductivity, demonstrating the calculation results of this invention. Detailed Implementation

[0133] The present invention will now be described in further detail with reference to the accompanying drawings and specific embodiments.

[0134] like Figure 1 As shown, the present invention provides a method for calculating implicit heat transfer of moving particles using a multi-time-layer difference scheme, comprising the following steps:

[0135] Step 1: Arrangement of the computational domain and its boundaries. Before computation, arrange the particles of the computational domain according to its extent, and then arrange the boundary-related virtual particles according to each boundary of the computational domain. This includes three types of heat transfer boundary conditions: the first type of temperature boundary condition, the second type of heat flux boundary condition, and the third type of convective heat transfer boundary condition; and two types of geometric boundary conditions used to simplify the computational domain: symmetric boundaries and periodic boundaries. The specific functions of each boundary are as follows:

[0136] The first type of temperature boundary condition is applicable to heat transfer calculations under isothermal walls; the second type of heat flux boundary condition is applicable to surfaces with constant power input heat, among which the adiabatic boundary is a special type of heat flux boundary condition; the third type of convective heat transfer boundary condition is applicable to convective heat transfer between the heat transfer surface and the external medium; geometric boundary conditions are used for symmetrically arranged computational domains and for periodic boundaries used in infinitely large or periodic computational domains. By using geometric boundary conditions, only a portion of the computational domain needs to be arranged, thereby simplifying calculations, reducing computational costs, and accelerating calculation speed.

[0137] In particle arrangement, all three heat transfer boundary arrangements require three layers of virtual particles to be arranged outside the computational domain to satisfy the supporting domain, such as... Figure 2 (a), where the straight line is the boundary, the point on the boundary is the boundary node, and the area above the boundary is the internal particle. In the moving particle method, the interaction between particles is based on the neighbor domain. There are internal particles with missing support domains near the boundary. At this time, the particles outside the boundary (light-colored particles at the bottom of the boundary) are added. The additional virtual particles not only satisfy the support domain of the particles in the computational domain, but also serve as supplementary nodes for the realization of heat transfer boundary conditions.

[0138] The two geometric boundary conditions, based on the mapping generated by the geometric relationship, supplement the support domain. Therefore, after defining the geometric boundary, the distance from the boundary is mapped into the outside of the boundary as supplementary particles. These mappings are one-to-one, meaning the position parameters are determined, thus eliminating the need to arrange virtual particles. Figure 2(b) shows the arrangement of periodic and symmetric boundaries. Since the particles missing in the neighborhood of the internal particles at the boundary are all subsets of the internal particles due to the geometric relationships, there is no need to arrange virtual particles. The light-colored virtual particles in the figure are only drawn for explanation. The arrangement of geometric boundaries needs to be based on the overall simulation content and the requirements of simplifying the computational domain. The arrangement of periodic boundaries requires that the two boundaries be declared as periodic boundaries. Symmetric boundaries cannot be arranged on curved surfaces.

[0139] Step 2: Construct the relationship between virtual particles and internal particles based on the boundary layout;

[0140] During the calculation, the temperature values ​​of the virtual particles are not directly added to the equations formed by the thermal equations as source terms for solving. Instead, relationships between the particles are constructed through different selected boundaries, and all particles are listed in the thermal equations. The resulting large matrix for solving is then combined with the relationships between the particles at the boundaries. Specifically, particles are added near the boundaries: for heat transfer boundaries, virtual particles are used; for geometric boundaries, mappings generated by geometric relationships are used. These added particles are then replaced with equations using selected boundary conditions. There are five relationships in total for three types of heat transfer boundaries and two types of geometric boundaries, as follows:

[0141] For the first type of temperature boundary condition, the temperature relationship between particles inside the computational domain and virtual particles in neighboring domains, as well as the distances of both types of particles from the boundary, are as follows:

[0142]

[0143] In the formula

[0144] i — an internal particle of the computational domain

[0145] j — Particles within the neighborhood of the internal particle. dummy — Refers to the set of virtual particles that supplement the boundary under the following boundary conditions: temperature boundary condition, heat flux boundary condition, convective heat transfer boundary condition, symmetric boundary condition, or periodic boundary condition.

[0146] j,dummy — Virtual particle j within the neighborhood of internal particle i.

[0147] T Dummy —Temperature at which virtual particles are used to fill the gaps where particles are missing at the boundary

[0148] T Boundary —Given boundary temperature

[0149] T Internal —Internal particle temperature of the computational domain

[0150] d j,dummy— The distance of virtual particle j within the neighborhood of internal particle i from the boundary.

[0151] d i — Distance of internal particle i from the boundary

[0152] For the second type of heat flux boundary condition and the third type of convective heat transfer boundary condition, the additional node method is used to perform Taylor expansion at the boundary, ensuring that virtual particles outside the boundary and particles inside the computational domain satisfy the applied boundary conditions. The particles inside the computational domain and virtual particles in neighboring domains have the following relationship:

[0153]

[0154]

[0155] In the formula

[0156] O(h 2 — Taylor expansion remainder

[0157] —Temperature gradient relationship at the boundary, For the second type of boundary For the setting of convective heat transfer boundaries, i.e., there is a relationship

[0158]

[0159] In the formula

[0160] T ∞ Ambient temperature

[0161] h—convective heat transfer coefficient

[0162] k—thermal conductivity of the material

[0163] When calculating the convection boundary, the physical property parameters required for the boundary arrangement need to be input in advance. For constant physical properties, the boundary arrangement can be directly solved by solving the virtual particle temperature outside the boundary. For polynomials with given physical properties and temperatures, algorithms for solving nonlinear problems are used for the solution.

[0164] Under symmetric and periodic boundaries, there exists a one-to-one mapping between particles at the boundary. The relationship between particles with one-to-one mapping at the geometric boundary is as follows:

[0165] T Symmetric,j′ =T Internal,j Formula (5)

[0166] T Periodic,j′ =T Internal,j Formula (6)

[0167] In the formula

[0168] TInternal,j —The temperature of particle j within the neighborhood of particle i, and particle j is within the computational domain.

[0169] T Symmetric,j′ —The temperature of particle j′ within the neighborhood of particle i, but particle j′ is not within the computational domain, and there is a mapping relationship between it and particle j due to the existence of the symmetry boundary.

[0170] T Periodic,j′ —The temperature of particle j′ within the neighborhood of particle i, but particle j′ is not within the computational domain, and there is a mapping relationship between the periodic boundary and particle j.

[0171] Step 3: Calculate the distance between the particle and the boundary;

[0172] Formulas (1), (2), (3) and (4) show that the application of the convective heat transfer boundary uses the Taylor expansion method based on the interparticle distance, so it is necessary to calculate the distance from the particles in the computational domain to the heat transfer boundary.

[0173] At a given convective heat transfer boundary, boundary nodes are generated based on particles within the computational domain. Boundary node curves and surface equations are then fitted using numerical methods. The generation principle for boundary nodes is to use the first layer of particles at the boundary, extrapolating the normal vector of each particle by half a particle radius to obtain a boundary node. After acquiring the position information of all boundary nodes, a polynomial curve or surface is fitted. In the two-dimensional simulation, the fitted boundary node curve is represented as:

[0174] z(x) = A + Bx + Cx 2 Formula (7)

[0175] In the 3D simulation, the fitted surface for the boundary nodes is written as:

[0176] z(x,y)=A+Bx+Cy+Dx 2 +Ey 2 +Fxy formula(8)

[0177] In the formula:

[0178] z(x) — the fitted curve function

[0179] z(x, y) — the fitted surface function

[0180] A, B, C, D, E, and F are all coefficients of the fitted curve and surface.

[0181] like Figure 3 As shown, the curve function fitted to the boundary nodes is given in the form of a polynomial function, which makes it easier to calculate the distance from the inner particles to the boundary.

[0182] In particle computation near the boundary, given a heat transfer boundary, it is necessary to obtain the particle distance to the boundary, that is, to solve the particle distance to the boundary, which is transformed into solving the x and y coordinate distance fitting boundary curve, i.e., a two-dimensional or surface, i.e., a three-dimensional distance problem.

[0183] Step 4: Update the physical properties based on the input initial temperature or the temperature of the previous time step; the time step length for a single calculation process is:

[0184]

[0185] In the formula:

[0186] Δt — the length of the time step in a single computation process

[0187] S—Dissipation number, a constant that controls the calculation to prevent divergence, determined according to the selected time difference scheme; l0—Initial particle arrangement distance in the particle method, also known as particle radius.

[0188] α max — The maximum thermal diffusivity of all particles in the computational domain at the previous or initial time step. Thermal diffusivity is the ratio of thermal conductivity to the product of density and specific heat.

[0189] Step 5: After obtaining the distance between the particle and the boundary, accurate heat transfer boundary conditions can be applied. Simultaneously, after obtaining the time step length of a single computational process, the heat transfer problem can be solved. In the particle method, solving the heat transfer problem is equivalent to solving the diffusion equation.

[0190]

[0191] In the formula

[0192] T—Particle temperature within the computational domain

[0193] —Internal heat source term, pre-defined based on specific calculation conditions.

[0194] ρ — density, given as a function of temperature

[0195] c — heat capacity, given as a function of temperature.

[0196] For each internal particle i that needs to be solved within the computational domain, a discrete equation can be written, that is, the governing equation is discretized as:

[0197]

[0198] In the formula:

[0199] d — the number of dimensions in the simulation

[0200] λ — A coefficient that makes the particle interaction model consistent with the analytical solution of the Gaussian function.

[0201] n0 — Initial particle number density

[0202] w(r ij ,r e — The weight function used in the particle method, r ij ,r e The independent variable of the weight function

[0203] r ij —Distance between internal particle i and neighboring particle j

[0204] r e —Given particle size

[0205] — This represents the solution for the temperature of internal particle i, where the superscript n+1 represents the result to be solved at the next time step.

[0206] — This represents the solution for the temperature of internal particle i, where the superscript n represents the known result of the previous time step.

[0207] — This represents solving for the temperature of particle j within the neighborhood of particle i, where the superscript n represents the known result of the previous time step.

[0208] k ij — Thermal conductivity, used in the calculation, is the harmonic mean of the internal particle i and its neighboring particles j.

[0209]

[0210] In the formula:

[0211] k i — Solving for the thermal conductivity of internal particle i

[0212] k j — Solving for the thermal conductivity of particle j in the neighborhood of particle i.

[0213] The difference method using a two-level time scheme implicitly defines the control equations. The multi-time-level difference scheme for the control equations is as follows:

[0214]

[0215] In the formula:

[0216] — Solve for the temperature of particle j within the neighborhood of particle i, where the superscript n+1 represents the result to be solved at the next time step.

[0217] f is a constant that controls the time difference scheme. f = 0 is an explicit difference scheme, f = 0.5 is a Crank-Nicholson difference scheme, and f = 1 is an implicit difference scheme. The difference between 0 and 0.5 is a partially explicit difference scheme, and the difference between 0.5 and 1 is a partially implicit difference scheme.

[0218] In explicit difference schemes, the coefficient S used to control the time step in step 4 needs to be less than 1. In the Crank-Nicholson difference scheme, the coefficient S used to control the time step needs to be less than 2. For implicit difference schemes, there is no limit to the size of the time step.

[0219] In various schemes, except for f=0 (i.e., explicit calculation does not require solving implicit equations), all other difference schemes require solving a system of equations for all particles within the computational domain, and then solving the matrix. When the neighboring particle j of the internal particle i is a supplementary particle j, dummy, based on the temperature relationship between the two particles (i.e., the boundary condition), the boundary particle representing the boundary condition is converted into a particle term within the computational domain. These terms are then placed into the diagonal and source terms of the matrix and solved consistently. The specific operation is as follows: For the first type of boundary, the matrix assembly process is as follows: If the neighboring particle i of the internal particle i has a supplementary particle j, dummy representing the first type of temperature boundary condition, the supplementary particle j, dummy is separated from the particles within the computational domain, i.e.:

[0220]

[0221] Based on formula (14), the relationship between the virtual particle j,dummy in the neighborhood domain and the i particle to be solved under the first type of heat transfer boundary is supplemented and written as:

[0222]

[0223] The operation is the same for the second and third types of heat transfer boundaries, and the boundary condition relations are used to replace the virtual particles j,dummy supplemented in the neighboring domain:

[0224]

[0225] For symmetric and periodic boundary conditions, since there is a mapping relationship between the supplementary particle j′ in the neighborhood and the j particle in the neighborhood of the i particle being solved, the boundary relationship is also substituted into the equation, and the symmetric boundary conditions are:

[0226]

[0227] The periodic boundaries are:

[0228]

[0229] Meanwhile, the geometric boundary is used to supplement the particle j′ term in the neighboring domain, and can be solved by merging it with the particle j term in the neighboring domain that is in the computational domain according to the mapping relationship;

[0230] Step 6: Solving the diffusion equation using a multi-time-layer difference scheme

[0231] After obtaining the difference scheme for heat transfer calculation in step 5, equations 15, 16, 17, and 18 are used. Since an equation can be written for each particle i in the computational domain, the equations formed by all particles in the solution domain are combined to obtain the heat diffusion equation set, which can be simplified to matrix form Ax = b, as follows:

[0232]

[0233] In the formula

[0234] b n —Source term value, where the subscript n represents the source term of the equation for the nth particle, b1 refers to the source term of the equation written for the 1st particle, b a The source term of the equation written for particle a —The temperature value to be solved. The subscript n indicates that the temperature value represents the temperature of particle n, and the superscript n+1 indicates that the temperature value A to be solved in the next time step. mn — The elements in the matrix are formed by the subscript mn, which represents the coefficient between particle m and particle n, that is, the coefficient before the temperature term of particle n in the diffusion equation for particle m; there are three types of coefficients in the coefficient matrix: diagonal elements, elements within neighboring particles, and elements outside neighboring particles.

[0235] The coefficient of the nth diagonal element is A nn For particle number 1, there is A 11 :

[0236]

[0237] For the element 'a' inside the neighboring particle and the element 'n' outside the neighboring particle, we have:

[0238]

[0239] A 1n =0 formula (22)

[0240] The source term is:

[0241]

[0242] In the moving particle method, during the calculation of the smallest unit, i.e., the particle, only the interaction between neighboring particles within the neighborhood is considered. Therefore, in the thermal diffusion equations, this is reflected in the fact that only the neighboring particle terms within the neighborhood are non-zero in each row, while all other terms are zero. Thus, the coefficient matrix of the thermal diffusion equations is sparse. At the same time, the neighborhoods of all particles are of equal size, and they are all neighbors of each other. Therefore, the coefficient matrix of the equations is symmetric. Considering that after substituting the boundary condition relations and simplifying by rearranging terms, only the diagonal terms of the coefficient matrix and the source terms of the thermal diffusion equations change, the solution of the boundary conditions does not cause a change in the symmetry of the thermal diffusion equations.

[0243] Since the coefficient matrix of the heat diffusion equations is a large symmetric sparse matrix, a large sparse symmetric matrix solver and a compatible solution algorithm, such as the conjugate gradient method, the conjugate gradient square algorithm, the biconjugate gradient method, the stable biconjugate gradient method, the conjugate residual method, the symmetric Langios algorithm, or the stable biconjugate residual algorithm, can be used to iteratively solve the equations.

[0244] The convergence criterion for the analytical solution of the thermal diffusion equations is:

[0245] ||r2<||b2ε formula (24)

[0246] In the formula:

[0247] ||r||2——The L2 norm of the residual vector, where the residual vector is calculated as r=b-Ax

[0248] ||b||2——The 2-norm of the source terms in the thermal diffusion equations

[0249] ε—a very small number, a constant used to represent the convergence criterion of the iteration. If the residual under the solution vector x after one iteration satisfies the above convergence condition, then it is considered convergent at this time step, and the calculation is performed for the next time step, such as... Figure 4 and Figure 5 The diagram shows the iterative convergence curves for a single time step for two basic examples, a semi-infinite flat plate and a square plate. It can be seen that the conjugate gradient method, the conjugate residual method, and the conjugate squaring method have good convergence, among which the conjugate squaring method has the fastest convergence speed.

[0250] Step 7: Advance the time step;

[0251] After a time step calculation is completed and the iterative solution converges, it is necessary to determine whether the calculation termination condition is met. If the termination condition is not met, the process returns to steps 4, 5, 6, and 7 to continue to the next time step until the calculation terminates. Figure 6 and Figure 7The results of the thermal conductivity calculation for a semi-infinite flat plate and a square plate are shown, along with temperature distribution contour maps. The calculation terminates when the steady-state result is reached. The results agree well with the analytical solution.

[0252] In summary, steps 1 to 3 input the boundary conditions of the specific problem to be calculated, determine the given parameters, generate particles within the computational domain, and form the boundaries of the computational domain by arranging particles and boundary nodes. During the neighbor particle search process, the distance of each particle at the boundary is determined, and the relationship between the boundary virtual particles and the particles within the computational domain is given. Step 4 updates the initial physical properties or the physical properties calculated in the previous time step, and gives a single time step for subsequent calculations that meet the stability conditions. Step 5 performs differential calculus on the governing equations of heat transfer calculation to form a large sparse symmetric coefficient matrix that can be solved. Step 6 solves the matrix formed in the previous step to obtain the temperature of each smallest computational unit particle in the next time step. The final step, step 7, determines whether the calculation meets the termination condition. If the termination condition is not met, steps 4, 5, 6, and 7 are returned to continue the next time step until the calculation terminates. This invention provides a more accurate boundary condition arrangement. The above steps can meet the requirements of calculations in most heat transfer fields.

Claims

1. A method for calculating implicit heat transfer of moving particles using a multi-time-layer difference scheme, characterized in that: The steps are as follows: Step 1: Arrangement of the computational domain and its boundaries. Before computation, arrange the particles of the computational domain according to its extent, and then arrange the boundary-related virtual particles according to each boundary of the computational domain. This includes three types of heat transfer boundary conditions: the first type of temperature boundary condition, the second type of heat flux boundary condition, and the third type of convective heat transfer boundary condition; and two types of geometric boundary conditions used to simplify the computational domain: symmetric boundaries and periodic boundaries. The specific functions of each boundary are as follows: The first type of temperature boundary condition is applicable to heat transfer calculations under isothermal walls; the second type of heat flux boundary condition is applicable to surfaces with constant power input heat, among which the adiabatic boundary is a special type of heat flux boundary condition; the third type of convective heat transfer boundary condition is applicable to convective heat transfer between the heat transfer surface and the external medium; geometric boundary conditions are used for symmetrically arranged computational domains and for periodic boundaries used in infinitely large or periodic computational domains. By using geometric boundary conditions, only a portion of the computational domain needs to be arranged, thereby simplifying calculations, reducing computational costs, and accelerating calculation speed. In particle arrangement, all three heat transfer boundary arrangements require three layers of virtual particles to be arranged outside the computational domain to satisfy the support domain. The additional virtual particles not only satisfy the support domain of the particles in the computational domain, but also serve as supplementary nodes for the realization of heat transfer boundary conditions. The two geometric boundary conditions supplement the support domain by mapping based on geometric relationships. Therefore, after defining the geometric boundary, the distance from the boundary is mapped into the outside of the boundary as supplementary particles. These mappings are one-to-one, which means that the position parameters are determined, so there is no need to arrange virtual particles. The geometric boundary needs to be arranged according to the overall simulation content and the requirements of simplifying the computational domain. The arrangement of periodic boundaries requires declaring that the two boundaries are periodic boundaries. Symmetrical boundaries cannot be arranged on curved surfaces. Step 2: Construct the relationship between virtual particles and internal particles based on the boundary layout; During the calculation, the temperature values ​​of the virtual particles are not directly added to the equations formed by the thermal equations as source terms for solving. Instead, relationships between the particles are constructed through different selected boundaries, and all particles are listed in the thermal equations. The resulting large matrix for solving is then combined with the relationships between the particles at the boundaries. Specifically, particles are added near the boundaries: for heat transfer boundaries, virtual particles are used; for geometric boundaries, mappings generated by geometric relationships are used. These added particles are then replaced with equations using selected boundary conditions. There are five relationships in total for three types of heat transfer boundaries and two types of geometric boundaries, as follows: For the first type of temperature boundary condition, the temperature relationship between particles inside the computational domain and virtual particles in neighboring domains, as well as the distances of both types of particles from the boundary, are as follows: In the formula i — an internal particle of the computational domain j — Particle j in the neighborhood of the internal particle dummy refers to the set of virtual particles that supplement the boundary under the following boundary conditions: temperature boundary condition (Type I), heat flux boundary condition (Type II), convective heat transfer boundary condition (Type III), symmetric boundary condition, or periodic boundary condition. j,dummy — Virtual particle j within the neighborhood of internal particle i. T Dummy —Temperature at which virtual particles are used to fill the gaps where particles are missing at the boundary T Boundary —Given boundary temperature T Internal —Internal particle temperature of the computational domain d j,dummy — The distance of virtual particle j within the neighborhood of internal particle i from the boundary. d i — Distance of internal particle i from the boundary For the second type of heat flux boundary condition and the third type of convective heat transfer boundary condition, the additional node method is used to perform Taylor expansion at the boundary, ensuring that virtual particles outside the boundary and particles inside the computational domain satisfy the applied boundary conditions. The particles inside the computational domain and virtual particles in neighboring domains have the following relationship: In the formula O(h 2 — Taylor expansion remainder —Temperature gradient relationship at the boundary, For the second type of boundary For the setting of convective heat transfer boundaries, i.e., there is a relationship In the formula T ∞ Ambient temperature h—convective heat transfer coefficient k—thermal conductivity of the material When calculating the convection boundary arrangement, the physical property parameters required for the boundary arrangement need to be input in advance. For constant physical properties, the boundary arrangement can be directly solved by solving the virtual particle temperature outside the boundary. For polynomials with given physical properties and temperatures, algorithms for solving nonlinear problems are used for the solution. Under symmetric and periodic boundaries, there exists a one-to-one mapping between particles at the boundary. The relationship between particles with one-to-one mapping at the geometric boundary is as follows: T Symmetric,j′ =T Internal,j Formula (5) T Periodic,j′ =T Internal,j Formula (6) In the formula T Internal,j —The temperature of particle j within the neighborhood of particle i, and particle j is within the computational domain. T Symmetric,j′ —The temperature of particle j′ within the neighborhood of particle i, but particle j′ is not within the computational domain, and there is a mapping relationship between it and particle j due to the existence of the symmetry boundary. T Periodic,j′ —The temperature of particle j′ within the neighborhood of particle i, but particle j′ is not within the computational domain, and there is a mapping relationship between the periodic boundary and particle j. Step 3: Calculate the distance between the particle and the boundary; Formulas (1), (2), (3) and (4) show that the application of the convective heat transfer boundary uses the Taylor expansion method based on the interparticle distance, so it is necessary to calculate the distance from the particles in the computational domain to the heat transfer boundary. At a given convective heat transfer boundary, boundary nodes are generated based on particles within the computational domain. Boundary node curves and surface equations are then fitted using numerical methods. The generation principle for boundary nodes is to use the first layer of particles at the boundary, extrapolating the normal vector of each particle by half a particle radius to obtain a boundary node. After acquiring the position information of all boundary nodes, a polynomial curve or surface is fitted. In the two-dimensional simulation, the fitted boundary node curve is represented as: z(x) = A + Bx + Cx 2 Formula (7) In the 3D simulation, the fitted surface for the boundary nodes is written as: z(x,y) = A + Bx + Cy + Dx 2 + Ey 2 + Fxy Formula (8) In the formula: z(x) — the fitted curve function z(x, y) — the fitted surface function A, B, C, D, E, and F are all coefficients of the fitted curve and surface. In particle computation near the boundary, given a heat transfer boundary, it is necessary to obtain the particle distance to the boundary, that is, to solve the particle distance to the boundary, which is transformed into solving the x and y coordinate distance fitting boundary curve, i.e., a two-dimensional or surface, i.e., a three-dimensional distance problem. Step 4: Update the physical properties based on the input initial temperature or the temperature of the previous time step; the time step length for a single calculation process is: In the formula: Δt — the length of the time step in a single computation process S—Dissipation number, a constant that controls the calculation to prevent divergence, determined according to the selected time difference scheme; l0—Initial particle arrangement distance in the particle method, also known as particle radius. α max — The maximum thermal diffusivity of all particles in the computational domain at the previous or initial time step. Thermal diffusivity is the ratio of thermal conductivity to the product of density and specific heat. Step 5: After obtaining the distance between the particle and the boundary, accurate heat transfer boundary conditions can be applied. Simultaneously, after obtaining the time step length of a single computational process, the heat transfer problem can be solved. In the particle method, solving the heat transfer problem is equivalent to solving the diffusion equation. In the formula T—Particle temperature within the computational domain —Internal heat source term, pre-defined based on specific calculation conditions. ρ — density, given as a function of temperature c — heat capacity, given as a function of temperature. For each internal particle i that needs to be solved within the computational domain, a discrete equation can be written, that is, the governing equation is discretized as: In the formula: d — the number of dimensions in the simulation λ — A coefficient that makes the particle interaction model consistent with the analytical solution of the Gaussian function. n0 — Initial particle number density w(r ij ,r e — The weight function used in the particle method, r ij ,r e The independent variable of the weight function r ij —Distance between internal particle i and neighboring particle j r e —Given particle size T i n+1 — This represents the solution for the temperature of internal particle i, where the superscript n+1 represents the result to be solved at the next time step. T i n — This represents the solution for the temperature of internal particle i, where the superscript n represents the known result of the previous time step. — This represents solving for the temperature of particle j within the neighborhood of particle i, where the superscript n represents the known result of the previous time step. k ij — Thermal conductivity, used in the calculation, is the harmonic mean of the internal particle i and its neighboring particles j. In the formula: k i — Solving for the thermal conductivity of internal particle i k j — Solving for the thermal conductivity of particle j in the neighborhood of particle i. The difference method using a two-level time scheme implicitly defines the control equations. The multi-time-level difference scheme for the control equations is as follows: In the formula: — Solve for the temperature of particle j within the neighborhood of particle i, where the superscript n+1 represents the result to be solved at the next time step. f is a constant that controls the time difference scheme. f = 0 is an explicit difference scheme, f = 0.5 is a Crank-Nicholson difference scheme, and f = 1 is an implicit difference scheme. The difference between 0 and 0.5 is a partially explicit difference scheme, and the difference between 0.5 and 1 is a partially implicit difference scheme. In explicit difference schemes, the coefficient S used to control the time step in step 4 needs to be less than 1. In the Crank-Nicholson difference scheme, the coefficient S used to control the time step needs to be less than 2. For implicit difference schemes, there is no limit to the size of the time step. In various schemes, except for f=0 (i.e., explicit calculation does not require solving implicit equations), all other difference schemes require solving a system of equations for all particles within the computational domain, and then solving the matrix. When the neighboring particle j of the internal particle i is a supplementary particle j, dummy, based on the temperature relationship between the two particles (i.e., the boundary condition), the boundary particle representing the boundary condition is converted into a particle term within the computational domain. These terms are then placed into the diagonal and source terms of the matrix and solved consistently. The specific operation is as follows: For the first type of boundary, the matrix assembly process is as follows: If the neighboring particle i of the internal particle i has a supplementary particle j, dummy representing the first type of temperature boundary condition, the supplementary particle j, dummy is separated from the particles within the computational domain, i.e.: Based on formula (14), the relationship between the virtual particle j,dummy in the neighborhood domain and the i particle to be solved under the first type of heat transfer boundary is supplemented and written as: The operation is the same for the second and third types of heat transfer boundaries, and the boundary condition relations are used to replace the virtual particles j,dummy supplemented in the neighboring domain: For symmetric and periodic boundary conditions, since there is a mapping relationship between the supplementary particle j′ in the neighborhood and the j particle in the neighborhood of the i particle being solved, the boundary relationship is also substituted into the equation, and the symmetric boundary conditions are: The periodic boundaries are: Meanwhile, the geometric boundary is used to supplement the particle j′ term in the neighboring domain, and can be solved by merging it with the particle j term in the neighboring domain that is in the computational domain according to the mapping relationship; Step 6: Solving the diffusion equation using a multi-time-layer difference scheme After obtaining the difference scheme for heat transfer calculation in step 5, equations 15, 16, 17, and 18 are used. Since an equation can be written for each particle i in the computational domain, the equations formed by all particles in the solution domain are combined to obtain the heat diffusion equation set, which can be simplified to matrix form Ax = b, as follows: In the formula b n —Source term value, where the subscript n represents the source term of the equation for the nth particle, b1 refers to the source term of the equation written for the 1st particle, b a The source term of the equation written for particle a —The temperature value to be solved. The subscript n indicates that the temperature value represents the temperature of particle n, and the superscript n+1 indicates that the temperature value A to be solved in the next time step. mn — The elements in the matrix are formed by the subscript mn, which represents the coefficient between particle m and particle n, that is, the coefficient before the temperature term of particle n in the diffusion equation for particle m; there are three types of coefficients in the coefficient matrix: diagonal elements, elements within neighboring particles, and elements outside neighboring particles. The coefficient of the nth diagonal element is A nn For particle number 1, there is A 11 : For the element 'a' inside the neighboring particle and the element 'n' outside the neighboring particle, we have: A 1n = 0 Formula (22) The source term is: In the moving particle method, during the calculation of the smallest unit, i.e., the particle, only the interaction between neighboring particles within the neighborhood is considered. Therefore, in the thermal diffusion equations, this is reflected in the fact that only the neighboring particle terms within the neighborhood are non-zero in each row, while all other terms are zero. Thus, the coefficient matrix of the thermal diffusion equations is sparse. At the same time, the neighborhoods of all particles are of equal size, and they are all neighbors of each other. Therefore, the coefficient matrix of the equations is symmetric. Considering that after substituting the boundary condition relations and simplifying by rearranging terms, only the diagonal terms of the coefficient matrix and the source terms of the thermal diffusion equations change, the solution of the boundary conditions does not cause a change in the symmetry of the thermal diffusion equations. Since the coefficient matrix of the heat diffusion equations is a large symmetric sparse matrix, a large sparse symmetric matrix solver and a compatible solution algorithm, such as the conjugate gradient method, the conjugate gradient square algorithm, the biconjugate gradient method, the stable biconjugate gradient method, the conjugate residual method, the symmetric Langios algorithm, or the stable biconjugate residual algorithm, can be used to iteratively solve the equations. The convergence criterion for the analytical solution of the thermal diffusion equations is: ||r||2<||b||2ε Formula (24) In the formula: ||r||2——The L2 norm of the residual vector, where the residual vector is calculated as r=b-Ax ||b||2——The 2-norm of the source terms in the thermal diffusion equations ε—a very small number, used to represent the constant of the convergence criterion for iteration. If the residual under the solution vector x after one iteration satisfies the above convergence condition, then it is considered to have converged at this time step, and the calculation of the next time step is performed. Step 7: Advance the time step; After a time step calculation is completed and the iterative solution converges, it is necessary to determine whether the calculation termination condition is met. If the calculation termination condition is not met, return to steps 4, 5, 6, and 7 to continue to the next time step until the calculation terminates. In summary, steps 1 to 3 input the boundary conditions of the specific problem to be calculated, determine the given parameters, generate particles within the computational domain, and form the boundaries of the computational domain by arranging particles and boundary nodes. During the neighbor particle search process, the distance of each particle at the boundary is determined, and the relationship between the boundary virtual particles and the particles within the computational domain is given. Step 4 updates the initial physical properties or the physical properties calculated in the previous time step, and gives a single time step for subsequent calculations that meet the stability conditions. Step 5 performs differential calculus on the control equations for heat transfer calculation to form a large sparse symmetric coefficient matrix that can be solved. Step 6 solves the matrix formed in the previous step to obtain the temperature of each smallest computational unit particle in the next time step. The final step 7 determines whether the calculation meets the termination condition. If the termination condition is not met, steps 4, 5, 6, and 7 are returned to continue the next time step until the calculation terminates.

Images