Radiation heat transfer simulation method suitable for multi-angle parallel computing

Abstract

The invention discloses a radiation heat transfer simulation method suitable for multi-angle parallel computing, which comprises the following steps of: carrying out numerical discretization in each radiation transmission direction, mutually independent solving processes in each direction, and independently carrying out computing according to the number of CPU (Central Processing Unit) cores of a computer; meanwhile, the solving algorithm of the radiation intensity in the spatial domain is optimized, the original overall calculation of the spatial domain is optimized into successive calculation in the minimum calculation unit, the calculation efficiency is greatly improved, and the method is suitable for radiation transmission in various complex media and can provide technical support for thermal system analysis in the industrial field.

Description

Technical Field

[0001] The present invention belongs to the technical field of radiation heat transfer simulation, and in particular relates to a radiation heat transfer simulation method suitable for multi-angle parallel computing. Background Art

[0002] Radiative heat transfer, one of the three fundamental heat transfer mechanisms, is widely present in high-temperature and vacuum environments, playing a particularly important role in aerospace, nuclear energy, and energy engineering. Compared to conduction and convection, radiative heat transfer is independent of the medium and can be performed efficiently in a vacuum, thus playing a leading role in many extreme operating conditions. Radiation simulation technology, a key computational tool for simulating radiative transfer processes, has been widely applied in a variety of technical fields, including spacecraft thermal control, solar photovoltaic systems, laser transmission, and thermal management of nuclear reactors. Any object with a temperature above absolute zero (0K) emits thermal radiation, the intensity and characteristics of which depend on the object's temperature and emissivity. Based on the emissivity characteristics of different objects, radiating surfaces are classified as black bodies, gray bodies, and non-gray bodies. During radiative transfer, the absorption and scattering properties of the medium significantly affect the transfer of radiative energy, resulting in an uneven distribution of radiation in space. This poses a challenge for high-precision radiative heat transfer analysis.

[0003] Currently, the numerical methods for studying radiation intensity in three-dimensional space can be mainly divided into two categories: one is the numerical method based on ray tracing, which simulates the transmission of radiation energy by calculating spatial geometric relationships. Taking the patent "Surface-to-Surface Radiation Heat Transfer Algorithm Suitable for Large-Scale Parallel Computing" (patent number: CN109948259A) proposed by Niu Hongpan et al. in 2019 as an example, this method simulates the propagation of a large number of radiation rays in space and calculates the transfer process of radiation heat flux between surfaces, thereby solving the surface temperature. Although this method can provide higher calculation accuracy, due to the need to simulate a large number of rays and the large amount of calculation, it usually leads to a long calculation time, which is difficult to meet the real-time and large-scale computing requirements in engineering applications.

[0004] Another type of method is numerical methods based on solving partial differential equations. These methods use the physical principles of radiation intensity transmission in a medium to calculate the radiation transfer process and implement radiation simulation. Common computational methods in this category include the finite element method (FEM), the finite volume method (FVM), and the spectral element method (SEM). The element differential method (EDM) is a hybrid numerical method that combines the finite element method with the boundary element method (BEM). This method discretizes the computational domain into isoparametric elements and applies the finite element method to each isoparametric element, effectively simulating the radiation intensity distribution in three-dimensional space and calculating various radiation characteristic parameters. In "Direct Solution of Discrete Coordinate Radiation Transfer Equations by Chebyshev Collocation Point Spectral Method," Li (2010) used the spectral method to discretize the radiation transfer process into a matrix equation and successfully simulated the radiation transfer process in a three-dimensional rectangular furnace model. Although many numerical methods require complex geometric models to be discretized into detailed grids for calculation, as the industry's requirements for computing speed and efficiency continue to increase, there is an urgent need to design efficient numerical methods that can provide high-precision calculations and adapt to large-scale parallel computing needs. Summary of the Invention

[0005] In order to overcome the shortcomings of the existing technology, the present invention provides a radiation heat transfer simulation method suitable for multi-angle parallel computing. Numerical discretization is performed in each radiation transmission direction. The solution processes in each direction are independent of each other and are calculated independently according to the number of computer CPU cores. At the same time, the solution algorithm for radiation intensity in the spatial domain is optimized, and the original overall calculation of the spatial domain is optimized into successive calculations within the smallest computing unit, which greatly improves the computing efficiency. It is suitable for radiation transmission in a variety of complex media and can provide technical support for thermal system analysis in the industrial field.

[0006] The technical solutions adopted by the present invention to solve the technical problems are as follows:

[0007] Step 1: Determine the number of discrete radiation intensity angles, the number of grids in three-dimensional space, and the number of nodes within the grid based on the expected level of refinement of the simulation target, and determine the physical parameters included in the radiation transfer equation;

[0008] Step 2: Discretize the spatial computational domain according to the expected number of discrete grids and the number of nodes in the grid, and the spatial domain is decomposed into multiple units; the units are hexahedral units, and adjacent units are connected by shared common faces. The units can be passed through by direction vectors. According to the direction of the direction vector passing through the common face, the unit located at the starting position of the direction vector is the upstream unit, and the unit located at the end position of the direction vector is the downstream unit; sort the unit calculation order; continue to discretize nodes in the unit, and the nodes are divided into boundary points, internal points and connection points according to their positions in the unit; the boundary points are nodes distributed on the computational model spatial domain, the internal points are nodes distributed inside the unit, and the connection points are nodes distributed at the connection of two or more units.

[0009] Step 3: Discretize the spatial solid angle according to the expected number of discrete angles, and the radiation transfer equation is discretized into each The equation of direction, where Ω is the total variable of direction, is the circular angle variable, θ is the zenith angle variable, m is the serial number of the circular angle variable, n is the serial number of the zenith angle variable, and the direction indicated by the mth circular angle and the nth zenith angle is the direction variable Ω m,n , the radiation intensity calculation processes in a single direction are independent of each other, and the unit calculation order is determined according to the angle value; the nature of the connection point between the upstream unit and the downstream unit in step 2 is changed, and the connection point is changed to an internal point in the upstream unit and a boundary point in the downstream unit;

[0010] Step 4: Set the parallel strategy according to the number of CPU cores in the computer and solve each The equation of direction, for each Ω m,n The spatial domain units are sorted in the direction, and each unit in each direction is used as the minimum calculation unit. The spatial domain of the minimum calculation unit is interpolated using the Lagrange interpolation basis function, and the radiation transfer equation coefficient matrix A of the minimum calculation unit is solved one by one according to the calculation order; according to Ω m,n The boundary conditions in the direction are assigned to the boundary points, and the source matrix S of the radiation transfer equation on the node is solved according to the node type in the spatial unit;

[0011] Step 5: Each The initial minimum calculation unit of the direction equation is assigned the radiation intensity value of the boundary point by the boundary conditions, and the remaining internal points are assigned by the source term matrix S in step 4, and combined to form B in the iteration matrix of the current minimum calculation unit, and the matrix AI=B is constructed to solve and obtain the radiation intensity I in the current minimum calculation unit; after the calculation of this minimum calculation unit is completed, the downstream unit is found according to the calculation unit sorting result of step 2, and the radiation intensity value is assigned to the boundary point of the downstream unit. The calculation method of the minimum calculation unit in step 5 is repeated according to the calculation order to complete the radiation intensity calculation of each direction and each spatial unit;

[0012] Step 6: When the radiation intensity calculation in the spatial domain in all directions is completed, solve the input radiation G in the spatial domain and determine whether the error of G is less than the convergence condition. If it converges, end the iteration; otherwise, return to step 5 and iterate until the input radiation G meets the convergence condition, and output the radiation intensity I in the calculation domain.

[0013] Furthermore, the physical parameters include absorption coefficient, scattering coefficient, scattering phase function, refractive index and wall emissivity.

[0014] Furthermore, the radiation transmission process expression of the radiation intensity in three-dimensional space is:

[0015]

[0016] Where, κ a is the absorption coefficient; κ s is the scattering coefficient; r is the spatial coordinate of the radiation intensity, Indicates that the spatial variable is r, the angular variable Ω is the partial derivative with respect to the spatial coordinates x, y, and z, I(r,Ω) represents the radiation intensity value when the spatial variable is r and the angular variable is Ω; Ω is the angular coordinate of the radiation intensity, including the zenith angle θ and the circumferential angle The calculation interval is [0,π] and [0,2π], expressed as:

[0017]

[0018] μ, η, ξ represent the projection values ​​of the unit vector of the angle variable Ω in the x, y, and z directions, respectively; i, j, k represent the unit vectors in the x, y, and z directions, respectively; Г is the angular redistribution term caused by the gradient refractive index, which includes the angular partial differential and is expressed as:

[0019]

[0020] α, β and γ are the spatial derivatives of the refractive index n:

[0021]

[0022] S is represented by:

[0023]

[0024] I b is the blackbody radiation intensity at a point in space, and its relationship with temperature T is I b =σT 4 , σ is the Stefan-Boltzmann constant; Ω′ is the other incident directions; Φ(Ω,Ω′) is the spectral scattering phase function;

[0025] Define the blackbody radiation force of the medium as E b, and the relationship with the blackbody radiation intensity is:

[0026] E b =πI b .

[0027] Furthermore, the input radiation G is the radiation integral of the space point, which is expressed as:

[0028] G(r)=∫ 4π I(r,Ω)dΩ.

[0029] Furthermore, the boundary conditions are determined by the simulation target. For a semi-transparent diffuse reflection wall, the boundary conditions are:

[0030]

[0031] Where, I b,w (.) represents the blackbody radiation intensity at the spatial variable r on the wall boundary condition, I w (.) represents the radiation intensity value at the spatial variable r and the directional variable Ω in the wall boundary condition. The subscript w represents the wall, ε w is the wall emissivity, n w represents the unit vector of the exterior normal to the wall.

[0032] Furthermore, in step 3, the angle segment constant quadrature method is used to divide Ω into For discrete angles Angular redistribution term Г m,n is represented as:

[0033]

[0034] Where, are the weights of the discrete angles, as well as They are obtained by the following recursive formulas:

[0035]

[0036] Furthermore, in step 3, the radiation transfer equation is discretized into each The direction equation is expressed as:

[0037]

[0038] Where Г m,n and S m,n Respectively expressed as:

[0039]

[0040] Furthermore, in step 4, the Lagrange interpolation basis function is used to discretize the calculation unit, and the shape function K in the calculation unit is v The Lagrange interpolation basis function is expressed by interpolation of nodes on the x, y, and z axes. The total function after the basis functions in the three directions are orthogonal is expressed as:

[0041]

[0042] Where, They represent the Lagrange interpolation basis functions in the x, y, and z directions respectively, and the subscript v is the isoparametric element number;

[0043] For radiation intensity:

[0044]

[0045] Among them, i, j, k represent the node numbers of the nodes in the isoparametric element along the x, y, and z axes. Respectively represent the number of nodes of the calculation unit in the x, y, and z axes; use the isoparametric element to calculate the spatial partial derivative of the radiation intensity:

[0046]

[0047] Where δ is the local coordinate within the isoparametric element; J x is the Jacobian matrix of the isoparametric elements on the x-axis.

[0048] Furthermore, in step 4, after the minimum calculation unit is discretized by angle, the spatial discretization equation is expressed as:

[0049]

[0050] The above equation is expressed in each node within the computational unit, and the left side of the equation at each node is extracted. The coefficients of the matrix A v m,n , extract the right side of the equation corresponding to each node as matrix B v m,n , the matrix A of the assembled nodes v m,n and B v m,n , get the direction Ω m,n The matrix calculation form in the spatial calculation domain is:

[0051] A m,n I m,n =B m,n

[0052] The beneficial effects of the present invention are as follows:

[0053] This invention provides a three-dimensional spatial radiation heat transfer numerical simulation method that can perform parallel calculations across angles. It performs numerical discretization in each radiation transmission direction, with independent solutions in each direction, calculated independently based on the number of CPU cores in the computer. Furthermore, this invention optimizes the algorithm for solving radiation intensity in the spatial domain, reducing the original overall spatial domain calculation to sequential calculations within the smallest computational unit. This significantly improves computational efficiency and is applicable to radiation transmission within a variety of complex media. It can provide technical support for thermal system analysis in the industrial field. BRIEF DESCRIPTION OF THE DRAWINGS

[0054] Figure 1 is a flow chart of the method of the present invention;

[0055] Figure 2 The names of the three types of nodes and the schematic diagram of the calculation transfer process after the method of the present invention is discretized in the spatial domain;

[0056] Figure 3 It is an upstream and downstream relationship diagram of several types of isoparametric elements after discretization according to the method of the present invention;

[0057] Figure 4 This is a schematic diagram of parallel calculation in different directions of the method of the present invention;

[0058] Figure 5 This is a radiation simulation diagram of the method of the present invention in a model with a discontinuous refractive index medium. DETAILED DESCRIPTION

[0059] The present invention will be further described below with reference to the accompanying drawings and examples.

[0060] In order to realize the full process simulation of the three-dimensional space radiation heat transfer process in the spatial domain, the present invention provides a radiation simulation method suitable for parallel computing, and optimizes the unit calculation process to increase the calculation speed to meet the needs of the industrial field for thermal radiation calculation.

[0061] To achieve the above object, the present invention provides the following technical solution: a radiation heat transfer simulation method suitable for multi-angle parallel computing, comprising the following steps:

[0062] Step 1: Determine the number of discrete radiation intensity angles, the number of grids in three-dimensional space, and the number of nodes within the grid based on the expected level of refinement of the simulation target, and determine the physical parameters included in the radiation transfer equation, including the absorption coefficient, scattering coefficient, scattering phase function, refractive index, and wall emissivity.

[0063] Step 2: Discretize the spatial computational domain according to the expected number of discrete grids and the number of nodes within the grid. The spatial domain is first decomposed into multiple units, and nodes are further discretized within the units. Nodes are divided into boundary points, internal points, and connection points based on their positions within the units. Boundary points are nodes distributed in the computational model spatial domain, internal points are nodes distributed within the units, and connection points are nodes distributed at the connection between two or more units.

[0064] Step 3: Discretize the spatial solid angle according to the expected discrete angle number, and the radiation transfer equation is discretized into each Directional equations: The radiation intensity calculation processes in a single direction are independent of each other. The order of unit calculation is determined based on the angle value. The nature of the connection point between the upstream unit and the downstream unit in step 2 is changed. The connection point is changed to an internal point in the upstream unit and to a boundary point in the downstream unit.

[0065] Step 4: Set the parallel strategy according to the number of CPU cores in the computer and solve each The equation of direction, for each Ω m,n The spatial domain units are sorted in the direction, and each unit in each direction is used as the minimum calculation unit. The spatial domain of the minimum calculation unit is interpolated using the Lagrange interpolation basis function, and the radiation transfer equation coefficient matrix A of the minimum calculation unit is solved one by one according to the calculation order; according to Ω m,n The boundary conditions in the direction are assigned to the boundary points, and the source matrix S of the radiation transfer equation on the node is solved according to the node type in the spatial unit;

[0066] Step 5: Each The initial minimum calculation unit of the direction equation is assigned radiation intensity values ​​at the boundary points by the boundary conditions. The remaining internal points are assigned values ​​by the original matrix S in step 4, and combined to form B in the iteration matrix of the current minimum calculation unit. The matrix AI = B is constructed and solved to obtain the radiation intensity I in the current minimum calculation unit. After the calculation of this minimum calculation unit is completed, the downstream unit is found according to the calculation unit sorting result in step 2, and the radiation intensity value is assigned to the boundary point of the downstream unit. The calculation method of the minimum calculation unit described in step 5 is repeated according to the calculation order to complete the radiation intensity calculation of each direction and each spatial unit.

[0067] Step 6: When the radiation intensity calculation in the spatial domain in all directions is completed, solve the input radiation G in the spatial domain and determine whether the error of G is less than the convergence condition. If it converges, end the iteration; otherwise, return to step 5 and iterate until the input radiation G meets the convergence condition, and output the radiation intensity I in the calculation domain.

[0068] The process described in the above steps is as follows Figure 1 As shown in Figure 2, the radiation transmission process of radiation intensity in three-dimensional space is expressed as follows:

[0069]

[0070] Where, κ a is the absorption coefficient; κ s is the scattering coefficient, defined in step 1; r is the spatial coordinate of the radiation intensity; Ω is the angular coordinate of the radiation intensity, including the zenith angle θ and the circumferential angle The calculation interval is [0,π] and [0,2π], expressed as:

[0071]

[0072] Where Г is the angular redistribution term caused by the gradient refractive index, which includes the angular partial differential and is expressed as:

[0073]

[0074] The spatial derivatives α, β, and γ of the refractive index n described in step 1 are:

[0075]

[0076] S described in step 4 is expressed as:

[0077]

[0078] In the above formula, I b is the blackbody radiation intensity at a point in space, and its relationship with temperature T is I b =σT 4 , σ is the Stefan-Boltzmann constant; Ω′ is the other incident directions; Φ is the spectral scattering phase function.

[0079] The blackbody radiation force of the medium is E b , and its relationship with the blackbody radiation intensity is:

[0080] E b =πI b

[0081] The input radiation G described in step 6 is the radiation integral of the spatial point, expressed as:

[0082] G(r)=∫ 4π I(r,Ω)dΩ

[0083] The boundary conditions described in step 5 are determined by the simulation target. Taking the semi-transparent diffuse reflection wall as an example, its boundary conditions are:

[0084]

[0085] Where, subscript w represents the wall, ε w is the wall emissivity, determined by step 1, nw represents the unit vector of the exterior normal to the wall.

[0086] After clarifying the composition of the radiation transfer equation required for the simulation of the radiation transfer process, such as Figure 2 As shown, step 2 performs spatial discretization, then step 3 performs angle discretization, and step 4 performs step-by-step discrete solution of the radiation transfer equation.

[0087] In step 3, the angle segment constant quadrature method is used to divide Ω into For discrete angles Angular redistribution term Г m,n Can be expressed as:

[0088]

[0089] Where, are the weights of the discrete angles, as well as They are obtained by the following recursive formulas:

[0090]

[0091]

[0092] In step 3, the radiative transfer equation is discretized into The direction equation can be expressed as:

[0093]

[0094] Where Г m,n and S m,n Respectively expressed as:

[0095]

[0096] In step 4, Lagrange interpolation basis function is used to calculate the angle Discretize the spatial domain of the direction and calculate the shape function K within the unit v Expressed as:

[0097]

[0098] Where L represents the Lagrangian shape function and the subscript v is the isoparametric element number.

[0099] The physical quantity at any point in the minimum calculation unit can be expressed using the value on the node. Take radiation intensity as an example:

[0100]

[0101] Where i, j, k represent the node numbers of the nodes in the isoparametric element along the x, y, and z axes. Therefore, the isoparametric element can be used to calculate the spatial partial derivative of the radiation intensity:

[0102]

[0103] Where δ is the local coordinate within the isoparametric element; J x is the Jacobian matrix of the isoparametric element on the x-axis, and the transformation method of the y-axis and the z-axis is the same as the above formula.

[0104] Figure 3 It expresses the upstream and downstream unit relationship of the unit connected by the minimum computing unit connection point. The connection point is associated with all adjacent isoparametric element information. Since the radiation intensity propagates in the form of rays, at the connection point, the radiation intensity information is transmitted from the upstream isoparametric element to the downstream isoparametric element. The present invention defines the number of upstream units adjacent to the connection point as N. up .

[0105] The equation construction of the connection point requires the use of all upstream unit information. The discrete equation in the smallest calculation unit described in step 4 is expressed as:

[0106]

[0107] Assemble the matrix A of nodes v m,n and B v m,n , we can get the direction Ω m,n The matrix calculation form in the spatial calculation domain is:

[0108] A m,n I m,n =B m,n

[0109] As described in step 5, the radiation intensity I calculated by the upstream minimum calculation unit is passed to the downstream unit, which uses it as a boundary condition to assemble the matrix A in its calculation process. v m,n and B v m,n , completing the calculation transfer.

[0110] Step 4 and Step 5 Figure 4 As shown, as described in step 6, the cycle from step 4 to step 6 is repeated independently in each parallel direction until the input radiation G in the calculation domain of the calculation space meets the convergence condition, and the output calculation domain can represent the radiation heat transfer result.

[0111] Example:

[0112] Figure 5Describes the calculation process of radiation transfer in a complex refractive index medium, and calculates the optical thickness of the medium in the box as τ = κ a H=1, the temperature of the bottom of the box z=0 is 1000K, ε w = 1. The rest of the walls are black cold walls. H at the center of the box 3 The / 16 area is a uniform refractive index medium with n=2, and the rest is a linear refractive index medium with n=5-4 ​​(z / H). The grid used in the calculation is N e =3×3×3, N n =3×3×3,

[0113] The simulation results show that the method provided by the present invention can quickly and effectively simulate the radiation heat transfer process of complex refractive index media in three-dimensional space. The simulation results are consistent with the expected goals, and the feasibility and effectiveness of the theory and method proposed in the present invention have been verified.

Claims

1. A radiation heat transfer simulation method suitable for multi-angle parallel computing, characterized in that: The steps include: Step 1: Determine the number of discrete radiation intensity angles, the number of grids in three-dimensional space, and the number of nodes within the grid based on the expected level of refinement of the simulation target, and determine the physical parameters included in the radiation transfer equation; Step 2: Discretize the spatial computational domain according to the expected number of discrete grids and the number of nodes in the grid, and the spatial domain is decomposed into multiple units; the units are hexahedral units, and adjacent units are connected by shared common faces. The units can be passed through by direction vectors. According to the direction of the direction vector passing through the common face, the unit located at the starting position of the direction vector is the upstream unit, and the unit located at the end position of the direction vector is the downstream unit; sort the unit calculation order; continue to discretize nodes in the unit, and the nodes are divided into boundary points, internal points and connection points according to their positions in the unit; the boundary points are nodes distributed on the computational model spatial domain, the internal points are nodes distributed inside the unit, and the connection points are nodes distributed at the connection of two or more units. Step 3: Discretize the spatial solid angle according to the expected number of discrete angles, and the radiation transfer equation is discretized into each The equation of direction, where Ω is the total variable of direction, is the circular angle variable, θ is the zenith angle variable, m is the serial number of the circular angle variable, n is the serial number of the zenith angle variable, and the direction indicated by the mth circular angle and the nth zenith angle is the direction variable Ω m,n , the radiation intensity calculation processes in a single direction are independent of each other, and the unit calculation order is determined according to the angle value; the nature of the connection point between the upstream unit and the downstream unit in step 2 is changed, and the connection point is changed to an internal point in the upstream unit and a boundary point in the downstream unit; Step 4: Set the parallel strategy according to the number of CPU cores in the computer and solve each The equation of direction, for each Ω m,n The spatial domain units are sorted in the direction, and each unit in each direction is used as the minimum calculation unit. The spatial domain of the minimum calculation unit is interpolated using the Lagrange interpolation basis function, and the radiation transfer equation coefficient matrix A of the minimum calculation unit is solved one by one according to the calculation order; according to Ω m,n The boundary conditions in the direction are assigned to the boundary points, and the source matrix S of the radiation transfer equation on the node is solved according to the node type in the spatial unit; Step 5: Each The initial minimum calculation unit of the direction equation is assigned the radiation intensity value of the boundary point by the boundary conditions, and the remaining internal points are assigned by the source term matrix S in step 4, and combined to form B in the iteration matrix of the current minimum calculation unit, and the matrix AI=B is constructed to solve and obtain the radiation intensity I in the current minimum calculation unit; after the calculation of this minimum calculation unit is completed, the downstream unit is found according to the calculation unit sorting result of step 2, and the radiation intensity value is assigned to the boundary point of the downstream unit. The calculation method of the minimum calculation unit in step 5 is repeated according to the calculation order to complete the radiation intensity calculation of each direction and each spatial unit; Step 6: When the radiation intensity calculation in the spatial domain in all directions is completed, solve the input radiation G in the spatial domain and determine whether the error of G is less than the convergence condition. If it converges, end the iteration; otherwise, return to step 5 and iterate until the input radiation G meets the convergence condition, and output the radiation intensity I in the calculation domain.

2. The radiation heat transfer simulation method suitable for multi-angle parallel computing according to claim 1, characterized in that: The physical parameters include absorption coefficient, scattering coefficient, scattering phase function, refractive index and wall emissivity.

3. The radiation heat transfer simulation method suitable for multi-angle parallel computing according to claim 2, characterized in that: The radiation transmission process expression of the radiation intensity in three-dimensional space is: Where, κ a is the absorption coefficient; κ s is the scattering coefficient; r is the spatial coordinate of the radiation intensity, Indicates that the spatial variable is r, the angular variable Ω is the partial derivative with respect to the spatial coordinates x, y, and z, I(r,Ω) represents the radiation intensity value when the spatial variable is r and the angular variable is Ω; Ω is the angular coordinate of the radiation intensity, including the zenith angle θ and the circumferential angle The calculation interval is [0,π] and [0,2π], expressed as: μ, η, ξ represent the projection values ​​of the unit vector of the angle variable Ω in the x, y, and z directions, respectively; i, j, k represent the unit vectors in the x, y, and z directions, respectively; Г is the angular redistribution term caused by the gradient refractive index, which includes the angular partial differential and is expressed as: α, β and γ are the spatial derivatives of the refractive index n: S is represented by: I b is the blackbody radiation intensity at a point in space, and its relationship with temperature T is I b =σT 4 , σ is the Stefan-Boltzmann constant; Ω′ is the other incident directions; Φ(Ω,Ω′) is the spectral scattering phase function; Define the blackbody radiation force of the medium as E b , and the relationship with the blackbody radiation intensity is: AND b =πI b 。 4. The radiation heat transfer simulation method suitable for multi-angle parallel computing according to claim 3 is characterized in that: The input radiation G is the radiation integral of the space point, expressed as: G(r)=∫ 4π I(r,Ω)dΩ。 5. The radiation heat transfer simulation method suitable for multi-angle parallel computing according to claim 4, characterized in that: The boundary conditions are determined by the simulation target. For a semi-transparent diffuse reflection wall, the boundary conditions are: Where, I b,w (.) represents the blackbody radiation intensity at the spatial variable r on the wall boundary condition, I w (.) represents the radiation intensity value at the spatial variable r and the directional variable Ω in the wall boundary condition. The subscript w represents the wall, ε w is the wall emissivity, n w represents the unit vector of the exterior normal to the wall.

6. The radiation heat transfer simulation method suitable for multi-angle parallel computing according to claim 5, characterized in that: In step 3, the angle segment constant quadrature method is used to divide Ω into For discrete angles Angular redistribution term Г m,n is represented as: Where, are the weights of the discrete angles, as well as They are obtained by the following recursive formulas:

7. The radiation heat transfer simulation method suitable for multi-angle parallel computing according to claim 6, characterized in that: In step 3, the radiation transfer equation is discretized into The direction equation is expressed as: Where Г m,n and S m,n Respectively expressed as:

8. The radiation heat transfer simulation method suitable for multi-angle parallel computing according to claim 7, characterized in that: In step 4, the Lagrange interpolation basis function is used to discretize the calculation unit, and the shape function K in the calculation unit is v The Lagrange interpolation basis function is expressed by interpolation of nodes on the x, y, and z axes. The total function after the basis functions in the three directions are orthogonal is expressed as: Where, They represent the Lagrange interpolation basis functions in the x, y, and z directions respectively, and the subscript v is the isoparametric element number; For radiation intensity: Among them, i, j, k represent the node numbers of the nodes in the isoparametric element along the x, y, and z axes. Respectively represent the number of nodes of the calculation unit in the x, y, and z axes; use the isoparametric element to calculate the spatial partial derivative of the radiation intensity: Where δ is the local coordinate within the isoparametric element; J x is the Jacobian matrix of the isoparametric elements on the x-axis.

9. The radiation heat transfer simulation method suitable for multi-angle parallel computing according to claim 8, characterized in that: In step 4, after the minimum calculation unit is discretized by angle, the spatial discretization equation is expressed as: The above equation is expressed in each node within the computational unit, and the left side of the equation at each node is extracted. The coefficients of the matrix A v m,n , extract the right side of the equation corresponding to each node as matrix B v m,n , the matrix A of the assembled nodes v m,n and B v m,n , get the direction Ω m,n The matrix calculation form in the spatial calculation domain is: A m,n I m,n =B m,n 。

Images