Numerical Simulation of Thermal Field in Mass Concrete With Pipe Water Cooling

Water pipe cooling is mainly used to control temperature in the construction of mass concrete structures. It is important to reveal how to accurately stimulate the temperature field of mass concrete under action of this water pipe cooling. This paper presents a new method for this purpose. In this method, the contact surface of the water pipe and the concrete is used as the heat dissipation surface into the control equation and the composite Multiquadrics radial basis function (MQ-RBF) and low-order linear polynomial combination are used to discrete the spatial domain. The heat dissipation surface of the water pipe is included in the boundary conditions so that there is no need to build the refined water pipe modeling. This new method not only reduces the calculation cost but also ensures calculation accuracy. Through four calculation examples, this paper show that the algorithm has advantages in the numerical simulation of the concrete temperature field with water pipe cooling.


INTRODUCTION
Mass concrete can release a lot of hydration heat during the construction process, especially in its early stage. In this early stage of concrete construction, a higher temperature is generated through the violent chemical reaction and the poor thermal conductivity. Generally, the highest temperature on the surface of the mass concrete bearing platform is about 30°C, and the highest temperature inside is approximately 50°C. The surface is directly contact with the air to dissipate heat. Thus, the temperature difference between the inside and outside of the concrete becomes larger. As a result, it generates more significant tensile stress. For the safety of the structure, it is necessary to control the concrete temperature in the early stage of construction.
It is very difficult to accurately simulate the temperature field of the concrete around the cooling water pipe. The reason is that it is a typical problem of the heat conduction between circulating fluids and solids in small calibers . The diameter of the water pipe is much smaller than the size of the concrete structure. In addition, the temperature gradient around the water pipe is also large.
There are three main approaches to solving this problem. The first type is the equivalent adiabatic temperature rise method proposed by Zhu [1] to simulate the impact of water pipe cooling. The algorithm principle is to use the cooling water pipe as a heat loss term and to only consider the temperature field in the average sense of concrete. Although this method only needs to draw a simple grid to calculate the temperature field in the cooling influence of the water pipe, the drastic changes in the temperature near the water pipe were ignored. It cannot obtain an accurate temperature field around the water pipe. The second type is named a refine grid algorithm. This algorithm considers the actual existence of water pipes, and the temperature in the boundary surface-contact surface between concrete and water pipers as a known one, and arranges dense grids near the water pipes. It can accurately simulate the temperature gradient near the water pipe by methods of refining the grid, but the calculation and modelling cost is too high to achieve in the actual modelling. In this paper, the solution obtained by the fine grid is called a refined solution. The third type is to construct a special water pipe unit, whose purpose is to reduce the calculation scale and ensure a certain accuracy. Many people put forward different ways to solve this problem. Kim [2] regards water pipes as line elements and uses common nodes to simulate the heat exchange between water pipes and concrete. However, the water pipes in this method can only be arranged on grid nodes, and the calculation uses Lagrange units. In order to obtain accurate solutions, it is necessary to refine the grid around the water pipe. Zuo [3] et al. used an extended finite element method (FEM) to solve the water pipe cooling problem, but at present, he can only simulate the water pipe cooling problem in the two-dimensional state. It is noticeable that he did not give a method for the temperature field simulation in the three-dimensional state. Chen [4] proposed to place the water pipe inside the concrete unit, to set up a "water pipe embedded unit," and put the contact surface between the water pipe and concrete-the heat dissipation surface-into the governing equation. However, the unit around the water pipe is still linear. To accurately simulate the temperature around the water pipe, a refined mesh is needed.
In this paper, a meshless interpolation method based on the global weak-form RPIM is used to simulate the heat exchange process between water pipes and concrete. The RPIM was proposed to solve the singularity problem of the moment matrix in the polynomial PIM. Traditionally, the shape function of PIM was constructed by using a polynomial basis function. This method has a significant advantage in dealing with essential boundary conditions because its shape function has the δ function property. But in the process of calculating the shape function, its moment matrix is very prone to singularity. In order to improve this method, Liu [24] proposed the RPIM. The RPIM solves the problem of matrix singularity by coupling the polynomial basis function with the radial basis function, which retains the characteristics of the polynomial with the δ function property. Therefore, this paper uses the RPIM to simulate the temperature field of the mass concrete containing cooling water pipes.
Reference [29] uses the localized radial basis function collocation method (LRBFCM) to simulate the temperature field of the cooling water pipe. Although it also uses RBF to discretize the calculation domain, the method of [1] is used to deal with the cooling effect of water pipes, and the temperature field obtained is still an average temperature field so that the accurate temperature field around the water pipe cannot be calculated. In this paper, based on the literature [4], RPIM is used to simulate the concrete temperature field with cooling water pipes. This method not only ensures the accuracy near the water pipe but also reduces the calculation scale. Apart from that, in order to further improve the accuracy of the algorithm, adaptive shape parameters and the Two-Side Equilibration Method [30] are used.

Governing Equation
The heat transfer differential equation of concrete is shown in Eq. 1, and the calculation model is shown in Figure 1.
In this equation T is the unknown temperature; τ is the time; α is the concrete thermal diffusivity coefficient; θ is the adiabatic temperature rise.
The water pipes discussed in this article are of non-metallic material. Basis on [4], the inner boundary of the non-metallic water pipe is in contact with water, so the temperature is equal to the water temperature T w . The outside of the water pipe is in contact with the inner surface of the concrete, and the temperature is equal to the temperature T c of the inner surface of the concrete. Thus, the heat flux q from concrete to water pipe can be expressed as Among them: β′ is the equivalent heat transfer coefficient of water pipe and concrete, and h R 0 − r 0 is the thickness of water pipe. When h → 0, at this time β′ → ∞, the contact surface between the water pipe and the concrete is equivalent to the first type of boundary condition. λ 1 is thermal conductivity coefficient of non-metallic water pipes; R 0 is the outer diameter of the water pipe; r 0 is the inner diameter of the water pipe, and T w is the water temperature of a given function. The method for calculating the water temperature is referred to Zhu [31].
Then the initial conditions and boundary conditions of the concrete with cooling water pipe can be expressed as the following formula Frontiers in Physics | www.frontiersin.org August 2021 | Volume 9 | Article 716859 In this formula the initial temperature of the concrete is expressed in T 0 ; and T a is the atmospheric temperature; λ is the thermal conductivity coefficient of concrete; and β is heat transfer coefficient of concrete to air. Furthermore, Γ 2 represents the concrete adiabatic boundary. Γ 3 represents the boundary between concrete and air, and Γ 1 represents the boundary between water pipe and concrete.

Integration Scheme
Since the RPIM cannot convert the integrals in the solution domain into the sum of the integrals of the discrete units, this paper uses the background grid integral method to integrate the problem domain. The background grid integration method is to divide the problem domain into n regular background grids, and then use Gaussian integration to calculate numerically in each background grid, and finally superimpose the background grid integrals. In this paper, we call each background grid a background unit.
The background grid is used to integrate the function, and the formula is as follow In the above formula, N c is the number of background grids, and Ω k is the unit field of the k th background grid. N r is the number of Gaussian integration points in Ω k , w r is the rth Gauss weight coefficient, and |J rk | is the Jacobian matrix of the k th background grid at the integration points x r .

Water Pipe Cooling Scheme
This paper uses the method [4] proposed to deal with the cooling effect of water pipes. The idea is to use the contact surface between concrete and water pipe as the boundary condition for heat dissipation and include the heat dissipation surface into the governing Eq. 1. The heat dissipation coefficient of the contact surface is estimated according to the thickness of the water pipe and the thermal conductivity. The solution of the three-dimensional transient temperature field of concrete with cooling water pipes is equivalent to the extreme value of the following functional In the equation V e , S e , and S′ are the domain of the background unit, the heat release surface of concrete and air, and the heat dissipation surface of the water pipe and concrete respectively.
The temperature of any point in the element can be obtained by linear interpolation of the element nodes.
Among them: T is the overall node temperature; n is the number of field nodes in the support domain; T e is the field node temperature in the support domain and _ T is the derivative of temperature time. Substituting Eqs 8, 9 into Eq. 7, Simplify the Eq. 10, we can obtain Where H + G is global thermal transfer matrix, F is global temperature load matrix, and R is thermal capacity matrix. According to the variational principle, the variation of Eq. 11 is δΠ 0, and we can obtain: where [c] is the element selection matrix and the element coefficient matrix is the following expression The last term in the above Eqs 16, 17 is the correction term for the cooling water pipe. In the formula β β/cρ, c is the specific heat capacity of concrete, and ρ is the density of concrete.
During the construction of large-volume concrete, since the spacing of the cooling water pipes is much smaller than its length, and the heat of the concrete is mainly transmitted perpendicular to the water pipe, the integral of the water pipe contact surface can be converted into a one-dimensional integral along the direction of the water pipe. The expression of the correction term can be rewritten as Where l is the water pipe length in the background unit, and T w is the temperature of the water flow.

Time Discretization Scheme
The difference method is used to discretise the time term in Eq. 12, Firstly the time domain [0, t] is divided into nt time steps. Then the time domain becomes [τ 0 , τ 1 , /, τ nt ] and the step size of each step is Δτ τ i+1 − τ i . Assuming Δτ temperature changes linearly so Taking s 1 as the backward difference method and asking the above Eqs 21, 22 into Eq. 12, we can obtain By combining the initial and boundary conditions, the Eq. 23 can be used to solve the system of equations to obtain the field node temperature array at each time.

RADIAL POINT INTERPOLATION METHOD
The RPIM is formed by adding a set of PIM to the radial basis function (RBF), which overcomes the problem of matrix singularity caused by polynomial basis functions. The calculation principle of the RPIM is to take a set of field nodes in a small local domain for a calculation point and use a combination of an RBF and a polynomial to fit the value of the calculation point on the field node. RPIM interpolation can be expressed as: Where T(x) is the numerical solution of the calculation point. R i (x) is the RBF. m is the number of RBFs. p j (x)is a monomial of space coordinates, k is the number of polynomial basis functions.
x is the calculation point. a j and b j are unknown coefficients.
In order to obtain the unknown coefficients in the above Eq. 24, the RBF and the polynomial basis function are firstly determined. This paper uses the MQ-BRF method and first-order polynomial basis functions. The function expression is as follows: The expression of the MQ-RBFs function is: Frontiers in Physics | www.frontiersin.org August 2021 | Volume 9 | Article 716859 4 Polynomial basis function expression of first degree: x − y i is the distance between the calculation point x and the field node y i in the local domain. Then r i x − y i is used indicate the distance between them. w 0 and c are shape parameters of MQ-RBF. This paper obtains these parameters c by trying a series of data. l is the average distance of field nodes in the local domain. In order to obtain the unknown coefficients a i and b j , a local domain needs to be formed at the calculation point x, and RPIM uses the field nodes in the local domain to form n linear equations. Then Eq. 24 can be expressed in a matrix as the function value T s {T 1 , T 2 , ..., T n } T .  The moment matrix of RBFs is equal to R: The polynomial moment matrix is: Eq. 27 is written in matrix form: At this time, the unknown coefficients a and b are obtained, and the unknown coefficients a and b are brought back to Eq. 27 to obtain: Where Φ T (x) {ϕ 1 , ϕ 2 , /, ϕ n } is the shape function of RPIM. The derivative of the numerical solutions are as follows: Here z is a differential operator. This paper adopts two methods to improve the accuracy so that these make the numerical solution have better accuracy. Firstly, it uses adaptive shape parameters w in each supported domain Ω s . The formula is as follows: Where w 0 stands for the initial shape parameter. x is the calculation point in the support domain and x i is the field node in the support domain. Additionally, for the matrix in Eq. 12, this paper uses Liu [30] method to reduce the condition number of the matrix, thereby improving the stability of the algorithm.

NUMERICAL EXPERIMENTS
In this section, four examples are used to experimentally examine RPIM in terms of accuracy, stability, and convergence. To analyse the error of the numerical solution, relative error (RE), and the relative root mean square error (RMSE) are selected as the error-index.
T(x i ) is the analytic solution or the exact solution of the FEM at the calculation point x i , and T(x i ) is the numerical solution of the RPIM. N c is the total number of test points.

Example 1
We verify the effectiveness of the algorithm through a twodimensional calculation example. It supposes that the thin plate is a square and the side length is 1.6 m. The concrete domain uses 1,089 field nodes with an interval of 0.05 m. The adiabatic temperature rise of concrete is θ(τ) 15.32 × (1 − e −0.4τ 0.66 ) and thermal conductivity coefficient is α 0.1m 2 /d. The outer diameter of the water pipe is R 0 0.016m. The thickness is h 0.002m. The thermal conductivity of the water pipe is λ 1.66kJ/(m · h ·°C). The initial temperature is T 0 20°C and the initial water temperature is T w 20°C. The time step is Δτ 0.02day. Figure 2 shows the distribution of concrete temperature fields at different times. It can be seen that the results of the proposed method are basically consistent with those of the FEM. When simulating example 1, the FEM used 4,176 elements. The relative error is shown in Figure 3 also verifies the event, and the maximum error is near the water pipe and remains at about 2.3%. Other changes with time, the algorithm also maintains good stability. Through this example, the algorithm can accurately and effectively simulate the transient heat conduction problem of concrete with a cooling pipe in the rectangular region.

Example 2
The example two considers the temperature field of cubic concrete under a single water pipe. Selecting the concrete block whose calculation domain is Ω {(x, y, x)|0 < x, y, z < 1.5m}. The concrete domain uses 1,375 field nodes, of which 1,331 (11 × 11 × 11) are evenly arranged with a spacing of 0.15 m. Adding four nodes    around the water pipe is as shown in the picture of Figure 4A. The initial temperature of the concrete is 20°C, and the concrete surface is adiabatic. The material parameters of concrete and water pipes are shown in Table 1. The adiabatic temperature rise of concrete is θ(τ) 21.66 × (1 − e −0.49τ 0.56 ). The outer diameter of the water pipe is R 0 0.016m, the thickness is h 0.002m. The water temperature along the water pipe is obtained by linear interpolation of the inlet and outlet water temperature.   The flow velocity is q w 1.0m 3 /h and the Initial water temperature is T w 5°C. Concrete continues to flow for 20 days. The time step is Δτ 0.1day. The refined solution uses 11760 Element and is refined around the water pipe, as shown in the picture of Figure 4B.
In this case, selecting point 1(0.6,1.5,0.6), point 2(0.45,1.2,0.45), point 3(0.75,0.75,1.5), three points are as the test points. Figure 5A shows the numerical solution and the FEM fine solution of the test point over time. It can be concluded from the figure that the fine solution of the test point is basically consistent with the numerical solution. Figure 5B shows that the relative error between the fine solution and the numerical solution of the three test points is below 1E-2. Especially in the near water pipe point1 position, it still maintains better accuracy. Therefore, it is proved that the proposed water pipe embedding method is accurate and effective in combination with the RPIM method under the three-dimensional model.
The Example two also verifies the convergence of the algorithm. Table 2 shows the RMSE values of the three test points at different time steps. It can be seen from the table that as the step size Δτ decreases, the RMSE also decreases. It shows that the algorithm is gradually converging. In addition, when the number of nodes in the field is the same, when the time step is less than Δτ 0.2day, the RMSE values are basically equal, indicating that the numerical solution obtained by RPIM has converged when Δτ 0.2day. This calculation example takes into account the water temperature change along the route and the heat conduction of concrete containing the cooling water pipe is discussed under the third boundary condition. The result shows that the algorithm has good properties in precision, stability, and convergence under the above conditions.

Example 3
This example studies the influence of flow velocity, water temperature, and pipe diameter on the concrete temperature field based on example 2. Except for water flow velocity, water temperature, and pipe inner diameter, other parameters are the same as those in calculation example 2. The water flow time is 20 days.
Tables 3-5 respectively list the flow velocity, water temperature, and the change of the parameters of the pipe inner diameter.
In this example, when discussing the influence of parameter changes on the temperature field, the highest temperature of concrete is selected for discussion. Table 3 discusses the influence of cooling water pipe flow velocities changes on the temperature field and selects four different flow velocities. Figure 6A shows that the temperature of the concrete is greatly reduced after the cooling water pipe is added, and the maximum temperature is only about 30°C while the temperature field without adding water pipes still keeps rising after the 20 th day and is above 40°C. Under the latter three conditions, the temperature curve did not change much, and the maximum temperature difference was only 0.3°C, which indicates that when the flow velocities reached a certain level, increasing the flow velocities did not significantly improve the cooling effect. Table 4 shows the changes in water temperature parameters. Figure 6B shows the maximum temperature value after the water temperature changes. The maximum concrete temperature is 30.1°C under 5°C water temperature. And the maximum concrete temperature is 32.7°C under 20°C water temperature. It is not particularly obvious the difference in temperature between the highest temperatures from 5°C to 20°C. But after the 4th day, the water temperature becomes lower and then the faster temperature of the concrete will drop. The temperature difference of the highest temperature on the 20 th day is 11.1°C. It shows that the water temperature has a great influence on the change of concrete temperature. But, when the water temperature is too low, a large temperature gradient will be generated near the water pipe, which will occur a large local tensile stress and cause damage the structure. Therefore, the cooling water temperature should not be too low during actual construction. Table 5 shows the Change in r 0 of the water pipe. Figure 6C shows that changing the r 0 of the water pipe still has a greater impact on the temperature of the concrete. When r 0 gradually becomes smaller, the heat dissipation capacity of the water pipe becomes worse. The reason is that when the water pipe r 0 becomes smaller, its heat transfer coefficient β′ becomes smaller, and the amount of heat transferred per unit condition decreases, resulting in poor heat dissipation effect. The effective contact area of other water pipes with water and the decrease in the flow rate per unit time of the water pipe are also factors that cause the deterioration of the heat dissipation capacity of the water pipe. Therefore, if the r 0 of the water pipe increases, then the temperature of the concrete will drop accordingly.
The RMSE values in Tables 3-5 have always been maintained under 1E-2, which verified that this method has good accuracy, stability, and robustness for the heat transfer of the threedimensional concrete model below the third kind of boundary conditions under the condition of water flow.
Example 4 Based on the two calculation examples 2 and 3, we have added multiple water pipes to the concrete model in example four and have increased the boundary conditions of the heat exchange between the concrete surface and the air. The size of the concrete model is Ω {(x, y, z)|0 < x < 3m, 0 < y < 2m, 0 < z < 1m}. Arranged three water pipes in the direction parallel to the Y axis, it is as shown in Figure 7A. In this example, the model has a total of 2,046 (31 × 11 × 6) field nodes, of which 18 (6 × 3) water pipes are arranged. Figure 7B shows the layout of field nodes and boundary conditions. It is worth noting that the material parameters of concrete and water pipes are the same as those in Example 2, such as the water flow velocity and the initial water temperature. The boundary temperature of Γ 3 is T a 30°C. The initial temperature of concrete and water pipe is 20°C, and the other boundaries are insulated. At the same time, the water supply time is also time 20 days.
For the three-dimensional model, we choose point 1 (0.0, 2.0, 1.0), point 2 (0.6, 2.0, 0.6), point 3 (1.8, 0.8, 0.5) three points and the FEM fine solution for discussion. (The FEM calculation includes 17520 elements). Figure 8A shows the numerical solution and fine solution of RPIM under the time-varying three test points. Consequently, the result is that the numerical solution of the algorithm can still maintain high accuracy under the condition of multi-water pipe. It can also be concluded from the relative error graph in Figure 8B that the error between the test point and the fine solution is below 1E-2, especially at point 2 near the water pipe. In actual engineering calculations, the accuracy of the algorithm is sufficient. Table 6 shows the comparison between the numerical solution and the fine solution of the test points under different total numbers of field nodes. And Table 6 shows that as the total number of field nodes increases, the RMSE gradually decreases, which indicates that the convergence of the algorithm is also better under multiple water pipes. Furthermore, in terms of calculation time, the more the number of field nodes, the more time it takes. But compared with the FEM algorithm, it only takes 220 s to reach an accuracy of 1E-3.
This example also compares the maximum temperature difference ΔT between the inside and the surface of the concrete. The Maximum temperature difference including pipes is ΔT 19.4°C. According to the engineering experience, the safety value of the temperature difference between the inside and surface of the concrete should be less than 25°C, so the water pipe layout in this paper is in line with the actual project. In conclusion, for the temperature field problem of concrete with cooling water pipes, the method proposed in this paper has high accuracy for temperature simulation, and it is valuable to apply it to actual engineering.

CONCLUSION
This paper uses the global weak-form radial basis function method and the method of 'embedded water pipe' to simulate the temperature field of concrete containing cooling water pipes. By incorporating the contact surface between the water pipe and the concrete as the heat dissipation surface into the governing equation, this reduces the modelling of water pipes. Using a Two-Side Equilibration Method and adaptive shape parameter method to further improve the accuracy and stability of the algorithm, and the discussion of four examples, it can be concluded that the algorithm can provide accurate prediction results. In conclusion, the method proposed in this paper has advantages in simulating the temperature field of hydration heat of concrete containing cooling water pipes.
The following conclusions can be drawn: 1) This method not only reduces the calculation cost, but also maintains the calculation accuracy.
2) The method is simple and convenient applying to actual projects, which still has good accuracy.
3) The method has good performance in stability, accuracy, and convergence.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.