A Method for Analyzing Fluid Temperature Distribution in a U-Shaped Medium-Deep Buried Pipe Heat Exchanger

Abstract

This invention discloses a method for analyzing fluid temperature distribution in a U-shaped medium-deep buried pipe heat exchanger, belonging to the field of ground source heat pump technology. This invention assumes that the fluids in the injection well, production well, and connecting well of the U-shaped medium-deep buried pipe heat exchanger satisfy a one-dimensional heat transfer equation along the flow direction. It also assumes that the backfill material and surrounding soil in the three wells satisfy a one-dimensional radial heat conduction perpendicular to the fluid flow direction, and ignores the heat capacity of the insulation material. Then, a one-dimensional transient heat transfer equation for the fluid is established, and the heat transfer process in the backfill material and soil is analyzed based on an analytical model, thereby establishing a semi-analytical heat transfer model. Finally, the fluid is divided into several nodes along the flow direction, the simulation time is divided into several time nodes, the relevant equations are discretized into algebraic equations, and the iterative method is used to calculate the distribution of fluid temperature with time and space. This invention features simple modeling, low computational load, and high computational accuracy, and has high promotional value.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of ground source heat pumps, and more particularly relates to a U-shaped middle-deep buried pipe heat exchanger fluid temperature distribution analysis method. BACKGROUND

[0002] The U-shaped middle-deep buried pipe heat exchanger is a buried pipe heat exchanger that can effectively extract middle-deep geothermal energy, and has the advantages of high heat exchange efficiency and small land occupation. The U-shaped middle-deep buried pipe heat transfer model is the theoretical basis for rock-soil thermal property testing, heat exchanger design, heat exchanger performance prediction and optimization, and has important research value. The existing technology mainly uses a numerical model to simulate the U-shaped middle-deep buried pipe and calculate the fluid temperature distribution.

[0003] In September 2020, Mingzhi Yu et al. disclosed an invention patent (application number: 202010485475.X) named "Simplified calculation method for heat exchange performance of U-shaped middle-deep buried pipe heat exchanger". The invention divides the numerical calculation area of the U-shaped middle-deep buried pipe heat exchanger into a descending pipe area, a horizontal pipe area, and an ascending pipe area, and each area can be solved in two-dimensional cylindrical coordinates. The originally complex three-dimensional unsteady heat transfer problem is simplified into two-dimensional unsteady problems in three different areas, reducing a large number of grid nodes and improving the calculation efficiency. Therefore, the U-shaped middle-deep buried pipe heat exchanger can be quickly calculated, and the calculation efficiency is greatly improved.

[0004] In February 2021, Li Chao et al. published an article named "Research on the influence of U-shaped deep buried pipe cementing layer on buried pipe heat exchange performance" in the Journal of Solar Energy, Vol. 42, No. 2, pp. 267-273. The article establishes a three-dimensional full-size numerical calculation model and verifies the model in combination with a U-shaped deep buried pipe project. On this basis, the influence of buried pipe cementing layers with different thermal conductivities on the heat exchange capacity of the buried pipe is simulated and analyzed. The results show that when the thermal conductivity of the cementing layer is less than and close to the thermal conductivity of the rock-soil layer, the buried pipe heat exchange intensity increases significantly; the influence of the thermal conductivity of the cementing layer on the buried pipe heat exchange is related to the size of the surrounding vertical rock-soil thermal conductivity, and has no obvious relationship with the vertical temperature distribution of the stratum.

[0005] In August 2022, the article "Research on Heat Extraction Performance of Medium-Deep U-shaped Buried Pipe Heat Exchanger Based on Longitudinal Layering of Rock-soil" by Lingling Bao et al. was published in the 37th volume of the Advances in Geophysics, pages 1371-1378. To explore the influence of different factors on the heat transfer performance of medium-deep U-shaped buried pipe heat exchanger under different geological structure layering, a numerical heat transfer model of medium-deep U-shaped buried pipe heat exchanger based on longitudinal layering of rock-soil was established. Using MATLAB software, the influence of circulating flow rate, inlet temperature, horizontal pipe length, thermal conductivity of insulation material, and thermal conductivity of backfill material on the heat transfer performance of U-shaped buried pipe heat exchanger was simulated using the alternating direction implicit method. The results show that: (1) When the circulating fluid flow rate increases from 10m 3 / h to 30m 3 / h, the heat extraction gradually increases, but increasing the inlet temperature leads to a decrease in heat extraction performance. To ensure efficient operation of the heat pump unit and geothermal heat exchanger, the selection of circulating flow rate and the adjustment of reasonable inlet temperature according to actual working conditions need to be considered; (2) Increasing the length of the horizontal pipe significantly improves the performance of the buried pipe heat exchanger, and the length of the horizontal pipe can be appropriately increased within the economic cost and technical limits; (3) After setting an insulation layer at the outlet of the uplink pipe, the heat extraction increases, but as the thermal conductivity decreases, the increase in heat extraction gradually decreases, so using low thermal conductivity insulation material can reduce the reverse heat transfer phenomenon; (4) When other parameters remain unchanged, increasing the thermal conductivity of the backfill material can significantly increase the outlet water temperature of the heat exchanger and the heat extraction, so high thermal conductivity backfill material can be used to improve the heat extraction capacity of the heat exchanger.

[0006] However, the above three documents all use numerical models to calculate the fluid temperature distribution of U-shaped medium-deep buried pipe heat exchangers, which have the disadvantages of complex modeling and large amount of calculation. SUMMARY

[0007] 1. Problems to be solved

[0008] In order to overcome the shortcomings of the existing U-shaped medium-deep buried pipe heat exchanger fluid temperature distribution calculation method, such as complex modeling and large amount of calculation, the present application provides a U-shaped medium-deep buried pipe heat exchanger fluid temperature distribution analysis method, which calculates the fluid temperature distribution with space and time variation by establishing a semi-analytical heat transfer model. The present application considers the transient heat transfer process of fluid, backfill and rock-soil, and also considers the actual situation of the change of rock-soil thermal properties with depth. It is suitable for fluid temperature distribution calculation under the conditions of given inlet fluid temperature and inlet heat flow, and has the advantages of simple modeling, small amount of calculation and high precision.

[0009] 2. Technical scheme

[0010] In order to solve the above problems, the technical scheme adopted by the present application is as follows:

[0011] A U-shaped middle-deep ground heat exchanger fluid temperature distribution analysis method of the present application assumes that the fluid in the injection well 8, the connecting well 10 and the production well 9 of the U-shaped middle-deep ground heat exchanger satisfies the one-dimensional heat transfer equation along the flow direction, assumes that the backfill 3 in the three wells and the surrounding rock-soil 7 satisfy the one-dimensional radial heat conduction perpendicular to the fluid flow direction, and ignores the heat capacity of the thermal insulation material 6 in the three wells. Then, the U-shaped middle-deep ground heat exchanger can be simplified as Figure 1 Figure 2 . Then, the one-dimensional transient heat transfer equation of the fluid is established, and the heat transfer process in the backfill and the rock-soil is analyzed based on the analytical model, and then a semi-analytical heat transfer model is established. Finally, the fluid is divided into several nodes along the fluid flow direction, the simulation time is divided into several time nodes, the related equations are discretized into algebraic equations, and the iterative method is used to calculate the distribution of fluid temperature with time and space.

[0012] The establishment process of the semi-analytical heat transfer model of the U-shaped middle-deep ground heat exchanger is as follows:

[0013] After simplifying the U-shaped middle-deep ground heat exchanger, the internal fluid energy equation can be expressed as

[0014]

[0015] r i (z) - the inner radius of the injection pipe, the connecting pipe or the production pipe, which can vary with z;

[0016] C f - the volumetric specific heat capacity of the fluid;

[0017] r o (z) - the outer radius of the injection pipe, the connecting pipe or the production pipe, which can vary with z;

[0018] C p - the volumetric specific heat capacity of the injection pipe, the connecting pipe and the production pipe wall;

[0019] T f (z, t) - the fluid temperature, which varies with z and t;

[0020] t - time;

[0021] V - the volumetric flow rate of the fluid;

[0022] z - the coordinate along the fluid flow direction;

[0023] L t - the total length of the U-shaped middle-deep ground heat exchanger, and L t = L1 + L2 + L3, where L1, L2 and L3 are the lengths of the injection pipe, the connecting pipe and the production pipe, respectively;​

[0024] q p (z, t) - heat flux transferred to the fluid through the pipe, which varies with z and t:

[0025]

[0026] T g (z, t) - temperature of the inner surface of the backfill;

[0027] R fg (z) - thermal resistance between the inner surface of the backfill and the fluid:

[0028]

[0029] λ p - thermal conductivity of the injection pipe, connection pipe, and production pipe;

[0030] λ ins - thermal conductivity of the insulation material;

[0031] r ins - outer radius of the insulation material;

[0032] L ins - length of the insulation material;

[0033] h(z) - convective heat transfer coefficient between the fluid and the inner wall of the pipe;

[0034] The backfill and the rock soil both satisfy one-dimensional radial heat conduction, so the temperature T g (z, t) of the inner surface of the backfill can be expressed as follows:

[0035]

[0036] t m - the mth time node;

[0037] T0(z) - initial temperature;

[0038] t i - the ith time node;

[0039] G(z, t) - G function of the composite medium model, related to the thermal properties of the rock soil, used to analyze the one-dimensional radial heat conduction process in the backfill and the rock soil;

[0040] At the inlet of the U-shaped deep geothermal pipe, the boundary conditions are as follows:

[0041]

[0042] T in (t) - inlet fluid temperature;

[0043] Q out (t) - heat output power.

[0044] The initial conditions are as follows:

[0045] T f (z,t)| t＝0 = T0(z) (6)

[0046] The above equations (1) to (6) constitute a complete description of the heat transfer process of the U-shaped middle-deep ground heat exchanger.

[0047] The fluid temperature distribution solving process of the U-shaped middle-deep ground heat exchanger is as follows:

[0048] The fluid in the injection pipe, the connecting pipe and the production pipe is divided into N1, N2 and N3 segments along the z direction respectively, and the total number of fluid nodes in the three pipes is N = N1+N2+N3+1, and the z coordinate of the arbitrary nth fluid node along the z direction is set as z n .

[0049] The simulation time t t is equally divided into M time nodes, and the arbitrary mth time node is t m = mΔt, wherein Δt is the difference value of any two adjacent time nodes, and Δt = t t / M.

[0050] For the arbitrary mth time node and the nth (2≤n≤N) fluid node, the equation (1) is discretized, and the equations (2) and (4) are connected to obtain:

[0051]

[0052] Wherein, A(z n ), B(z n ), C(z n-1 ) and D(z n-1 ,t m ) are intermediate variables:

[0053]

[0054]

[0055]

[0056]

[0057] Wherein, the calculation formula of q p (z n-1 ,t m ) is as follows:

[0058] q p (z n-1 ,t m )=C(z n-1 )[T0(0.5z n +0.5z n-1 )-D(z n-1 ,t m )-0.5T f (z n ,t m )-0.5T f (z n-1 ,t m )] (12)

[0059] Equations (5) and (6) can be discretized as

[0060]

[0061] T f (z n ,t0)=T0(z n ),(1≤n≤N) (14)

[0062] Among the above intermediate variables, A(z n ), B(z n ) and C(z n-1 ) can be calculated according to equations (8), (9) and (10) respectively, and D(z n-1 ,t m ) can be calculated according to equation (11) and q m-1 (z,t) of t1-t p time nodes, wherein q p (z,t) can be calculated according to equation (12), and it is worth noting that q p (z,t0) of t0 time node is 0. Based on the fluid temperature distribution of t0 time node in equation (14), equations (7) and (13) are solved by iteration method, and then all fluid node temperatures of t1 time node are calculated, and q p (z,t1) of all fluid nodes are calculated according to equation (12); based on the fluid temperature distribution of t1 time node and q p (z,t1) of all fluid nodes, equations (7) and (13) are solved by iteration method, and then all fluid node temperatures of t2 time node are calculated, and q p (z,t2) of all fluid nodes are calculated according to equation (12); based on the fluid temperature distribution of t2 time node and q p (z,t1), q p(z, t2), the equations (7) and (13) are solved by iteration method, then all fluid node temperatures of t3 time node are calculated, and q of all fluid nodes are calculated according to equation (12) p (z, t3); in this way, fluid temperature distribution of all time nodes can be calculated.

[0063] 3. Advantages

[0064] Compared with the prior art, the present application has the following advantages:

[0065] (1) The fluid temperature distribution analysis method of the U-shaped middle-deep ground buried pipe heat exchanger has the advantages of simple modeling, small calculation amount, etc., because only one-dimensional grid division is performed on the fluid and the number of space nodes is small, so the method can be conveniently applied to engineering practice.

[0066] (2) The fluid temperature distribution analysis method of the U-shaped middle-deep ground buried pipe heat exchanger adopts the analysis method to analyze one-dimensional transient heat transfer in the backfill and the rock-soil, and adopts the numerical method to simulate one-dimensional transient heat transfer in the fluid, so that not only the fluid heat capacity, the pipe wall heat capacity and the backfill heat capacity are considered, but also the change of the rock-soil thermal properties with depth and the change of the heat flow transferred to the fluid through the pipe with time and depth are considered, so that the calculation precision is high enough to meet the precision requirement of engineering practice. BRIEF DESCRIPTION OF DRAWINGS

[0067] Figure 1 Fig. 1 is a schematic diagram of the U-shaped middle-deep ground buried pipe heat exchanger;

[0068] Figure 2 Fig. 2 is a schematic diagram of the U-shaped middle-deep ground buried pipe heat exchanger simplified by the present application;

[0069] Figure 3 Fig. 3 is a node division of the fluid of the U-shaped middle-deep ground buried pipe heat exchanger according to the present application;

[0070] Figure 4 Fig. 4 is a comparison diagram of the outlet fluid temperature distribution calculation value, the experimental result and the results of other models in embodiment 1 of the present application;

[0071] Figure 5 Fig. 5 is a comparison diagram of the outlet fluid temperature distribution calculation value, the numerical model result and the results of other models in embodiment 2 of the present application;

[0072] Figure 6 Fig. 6 is a comparison diagram of the fluid temperature distribution calculation value, the numerical model result and the results of other models in embodiment 2 of the present application;

[0073] Figure 7 Fig. 7 is a comparison diagram of the fluid temperature distribution along the flow direction in embodiment 2 of the present application.

[0074] Explanation of the labels in the diagram:

[0075] 1. Fluid; 2. Injection pipe; 3. Backfill material; 4. Connecting pipe; 5. Production pipe; 6. Insulation material; 7. Rock and soil; 8. Injection well; 9. Production well; 10. Connecting well. Detailed Implementation

[0076] To further understand the content of this invention, a detailed description of the invention will be provided in conjunction with the accompanying drawings and embodiments.

[0077] Example 1

[0078] Figure 1 This is a schematic diagram of a U-shaped medium-deep buried pipe heat exchanger. Figure 2 This is a simplified schematic diagram of the U-shaped medium-deep buried pipe heat exchanger according to the present invention. Figure 3 This invention relates to the node division of fluid in a U-shaped medium-deep buried pipe heat exchanger. This embodiment calculates the outlet fluid temperature distribution over 0-72 hours for a specific U-shaped medium-deep buried pipe heat exchanger and compares the calculated values ​​with experimental results and other models (analysis models based on the infinite linear heat source theory). The parameters of the U-shaped medium-deep buried pipe heat exchanger tested in the field are shown in Table 1.

[0079] Table 1. Parameters of U-shaped medium-deep buried pipe heat exchangers tested in the field.

[0080]

[0081] Combination Figures 1-3 For injection well 8 and connecting well 10, the inner radius r of injection pipe 2, connecting pipe 4 or production pipe 5 i (z) = 0.0622m, the outer radius r of injection pipe 2, connecting pipe 4, or output pipe 5 o (z) = 0.0699m; for producing well 9, r i (z) = 0.0797m, r o (z) = 0.0889m. The inner radius r of the backfill material 3. g (z) is

[0082]

[0083] The initial temperature T0(z) is

[0084]

[0085] T sur —Surface temperature;

[0086] a—Geothermal gradient;

[0087] θ — the angle between the connecting tube and the output tube.

[0088] For the aforementioned U-shaped medium-deep buried pipe heat exchanger, this embodiment divides the fluid in the injection pipe and connecting pipe into 34 and 25 equal segments along the z-direction, respectively, i.e., N1 = 34 and N2 = 25; the fluid in the output pipe is divided into 56 segments along the z-direction, wherein the uninsulated segment and the insulated segment are each divided into N... 31 =44 and N 32 =12 segments, i.e., N3 = 56. Therefore, the total number of fluid nodes is N = N1 + N2 + N3 + 1 = 116, and the z-coordinate of any nth fluid node is...

[0089]

[0090] Due to simulation time t t If the time interval is 72 hours, it can be divided into 1080 equal time points, i.e., M = 1080, then Δt = t t / M=240s.

[0091] Then, calculate A(z) according to equations (8), (9) and (10) respectively. n B(z) n ) and C(z) n-1 ):

[0092]

[0093]

[0094]

[0095] Based on the fluid temperature distribution at time node t0 in equation (14), equations (7) and (13) are solved using the Gaussian iteration method, and then the temperatures of all fluid nodes at time node t1 are calculated. Finally, the q of all fluid nodes is calculated according to equation (12). p (z,t1), i.e., q p (z1,t1),q p (z2,t1)…q p (z N-1 ,t1); Fluid temperature distribution based on time node t1 and q of all fluid nodes p (z,t1), calculate D(z) according to equation (11). n-1 The Gaussian iteration method is used to solve equations (7) and (13) at time t2, and then the temperature of all fluid nodes at time t2 is calculated. The q of all fluid nodes is calculated according to equation (12). p (z,t2), i.e., q p (z1,t2),q p (z2,t2)…q p (z N-1,t2); Fluid temperature distribution based on time node t2 and q of all fluid nodes p (z,t1), q p (z,t2), calculate D(z) according to equation (11). n-1 The Gaussian iteration method is used to solve equations (7) and (13) at time t3, and then the temperature of all fluid nodes at time t3 is calculated. The q of all fluid nodes is calculated according to equation (12). p (z,t3), i.e., q p (z1,t3),q p (z2,t3)…q p (z N-1 ,t3); and so on, all temperatures of 116 fluid nodes at 1080 time nodes can be calculated.

[0096] The outlet fluid temperature (i.e., T) at all time points calculated in this invention f (z N Comparison of ,t) with experimental results and other models, for example Figure 4 As shown, it is evident that the calculated results of this invention generally follow the experimental trend, while the calculated results of the analytical model based on the infinite linear heat source theory show a significant difference from the experimental trend, especially over a shorter time period. Furthermore, the outlet fluid temperature calculated by this invention shows a high degree of agreement with the experimental results. Additionally, comparisons at some time points are shown in Table 2. It is clear that the calculation accuracy of this invention is far superior to that of the analytical model based on the infinite linear heat source theory.

[0097] Table 2 compares the outlet fluid temperature calculated by this invention at some time points with experimental results and other models.

[0098]

[0099] Example 2

[0100] Example 1 addresses the condition of a known inlet fluid temperature, while this example addresses the condition of a known thermal output power Q. out (t) Conditions. The parameters of the U-shaped medium-deep buried pipe heat exchanger in this embodiment are basically the same as those in Table 1, except that the inlet fluid temperature T in this embodiment is different. in (t) is unknown, while Q out (t) is known, and Q out (t)=100L t =574200W. This invention is used to calculate the outlet fluid temperature distribution over 120 days and the fluid temperature distribution at different time points, and the calculation results are compared with a numerical heat transfer model and other models (analysis models based on the infinite linear heat source theory).

[0101] Due to simulation time t t If the time period is 120 days, it can be divided into 2880 time points, i.e., M = 2880, then Δt = t t / M=3600s. The method for mesh generation of the fluid in this embodiment is the same as in Embodiment 1, and the related calculation steps are also basically the same as in Embodiment 1. The difference lies in A(z) n ) and C(z) n-1 The calculation results are different:

[0102]

[0103]

[0104] The comparison results of the outlet fluid temperature distribution are as follows: Figure 5 and Figure 6 As shown. The outlet fluid temperature calculated by this invention shows a high degree of agreement with the numerical model results, and the deviation between the two is generally less than 0.04℃ over the 120-day simulation period. The calculation results based on the infinite linear heat source theory show a lower degree of agreement with the numerical model results, and the deviation between the two is greater than 1℃ within one day, while the deviation over 120 days is 0.13℃. A comparison of the fluid temperature distribution along the flow direction is shown below. Figure 7 As shown, the fluid temperature distributions calculated by this invention at three different time points (10 hours, 80 hours, and 120 days) are in good agreement with the numerical model results. However, the calculation results of the analytical model based on the infinite linear heat source theory show a lower agreement with the numerical model results, especially over shorter time periods. These results demonstrate that the calculation accuracy of this invention is high and meets the requirements of engineering practice.

Claims

1. A method for analyzing fluid temperature distribution of a U-shaped middle-deep ground heat exchanger, characterized in that: Assuming that the fluid in the U-type middle-deep ground heat exchanger satisfies the one-dimensional heat transfer equation along the flow direction, the backfill and the surrounding rock-soil satisfy the one-dimensional radial heat conduction perpendicular to the fluid flow direction, and the heat capacity of the thermal insulation material is ignored; then, the one-dimensional transient heat transfer equation of the fluid is established, and the heat transfer process in the backfill and the rock-soil is analyzed based on the analytical model, and then the semi-analytical heat transfer model is established; finally, the fluid is divided into several nodes along the flow direction, the simulation time is divided into several time nodes, the related equations are discretized into algebraic equations, and the iterative method is used to calculate the distribution of the fluid temperature with time and space; The semi-analytical heat transfer model of the U-type middle-deep ground heat exchanger is as follows: The one-dimensional transient heat transfer equation of the fluid is r i (z) - inner radius of the injection pipe, connection pipe or production pipe; C f - the volumetric specific heat capacity of the fluid; r o (z) - outer radius of the injection pipe, connection pipe or production pipe; C p - volumetric specific heat capacity of the injection pipe, connection pipe and production pipe walls; T f (z, t) - fluid temperature; t——time; V——volume flow rate of the fluid; z——coordinate along the flow direction of the fluid; L t - the total length of the U-shaped middle-deep buried pipe, and L t = L1+L2+L3, wherein L1, L2 and L3 are the lengths of the injection pipe, the connecting pipe and the production pipe, respectively; q p (z, t) - heat flow transferred to the fluid by the tube: T g (z, t) - backfill interior surface temperature; R fg (z) Thermal resistance between the inner surface of the backfill and the fluid: λ p - thermal conductivity of the injection pipe, connection pipe and production pipe; λ ins - thermal conductivity of the insulating material; r ins - the outer radius of the thermal insulation material; L ins — length of the thermal insulation material; h(z)——convective heat transfer coefficient between the fluid and the inner wall of the pipe; The backfill and the rock soil satisfy one-dimensional radial heat conduction, and the inner surface temperature T g (z, t) can be expressed as follows: t m - mth time node; T0(z)——initial temperature; t i - the i-th time node; G(z,t)——G function of the composite medium model, related to the thermal properties of the rock-soil, used to analyze the one-dimensional radial heat conduction process in the backfill and the rock-soil; At the inlet of the U-type middle-deep ground heat exchanger, the boundary conditions are as follows: T in (t) — inlet fluid temperature; Q out (t) - thermal output power; The initial conditions are as follows: T f (z,t)| t＝0 = T0(z) (6).

2. The method according to claim 1, wherein the U-shaped middle-deep ground heat exchanger is characterized in that: The fluids in the injection pipe, the connecting pipe and the production pipe are divided into N1, N2 and N3 segments along the z direction respectively, the total number of fluid nodes in the three pipes is N=N1+N2+N3+1, and the z coordinate of the nth fluid node along the z direction is set as z n .

3. The method according to claim 2, wherein: The simulation time t t is equally divided into M time nodes, then any mth time node is t m = mΔt, where Δt is the difference value of any two adjacent time nodes, and Δt = t t / M.

4. The method according to claim 3, wherein: For any mth time node and nth fluid node, 2≤n≤N, equation (1) is discretized, and equations (2) and (4) are simultaneously solved to obtain where A(z n ), B(z n ), C(z n-1 ) and D(z n-1 ,t m ) are intermediate variables: where q p (z n-1 ,t m ) is calculated as follows: q p (z n-1 ,t m ) = C(z n-1 )[T0(0.5z n + 0.5z n-1 ) - D(z n-1 ,t m ) - 0.5T f (z n ,t m ) - 0.5T f (z n-1 ,t m )] (12) Equations (5) and (6) can be discretized as T f (z n ,t0)=T0(z n ),(1≤n≤N) (14) The intermediate variables A(z n ), B(z n ) and C(z n-1 ) can be calculated according to equations (8), (9) and (10) respectively, and D(z n-1 ,t m ) can be calculated according to equation (11) and q p (z,t) at time nodes t1-t m-1 , wherein q p (z,t) can be calculated according to equation (12), q p (z,t0)=0 at time node t0; based on the fluid temperature distribution at time node t0 in equation (14), equations (7) and (13) are solved by iteration to calculate all fluid node temperatures at time node t1, and q p (z,t1) of all fluid nodes is calculated according to equation (12); based on the fluid temperature distribution at time node t1 and q p (z,t1) of all fluid nodes, equations (7) and (13) are solved by iteration to calculate all fluid node temperatures at time node t2, and q p (z,t2) of all fluid nodes is calculated according to equation (12); based on the fluid temperature distribution at time node t2 and q p (z,t1), q p (z,t2) of all fluid nodes, equations (7) and (13) are solved by iteration to calculate all fluid node temperatures at time node t3, and q p (z,t3) of all fluid nodes is calculated according to equation (12); and so on, the fluid temperature distribution at all time nodes can be calculated.

Images