ORIGINAL RESEARCH article

Front. Phys., 24 May 2024

Sec. Mathematical Physics

Volume 12 - 2024 | https://doi.org/10.3389/fphy.2024.1338718

Simulation of thermal field in mass concrete with cooling pipes based on the isogeometric analysis method

  • College of Mechanics and Materials, Hohai University, Nanjing, China

Abstract

As the water pipe cooling system is widely applied to controlling temperature in mass concrete structures, the precise simulation of the temperature field in mass concrete with cooling pipes embedded is meaningful. This paper presents an isogeometric analysis (IGA) with NURBS for heat transfer in mass concrete with consideration of the cooling pipe. The proposed method not only achieves the same level of accuracy with fewer nodes but also eliminates the time-consuming process of mesh in the traditional FEM. The coarsest parameter space which depicts small pipe and large concrete precisely is constructed to create an efficient model for numerical computation. In addition, the unique k-refinement in IGA is supposed to be the most appropriate encryption mechanism, and the knot insertion vector for effective refinement is calculated by considering the characteristics of temperature gradient distribution around the cooling pipes. In addition, a different calculation parameter has been discussed to show the stability and flexibility of the IGA. The obtained numerical results demonstrate the accuracy and efficiency of the proposed scheme in the simulation of transient temperature fields in concrete structures with cooling systems.

1 Introduction

The embedded cooling pipe has been extensively used in engineering structures with mass concrete due to the significant effect on heat removal and crack control. It can effectively reduce the temperature difference between the surface and inner of the concrete structures, especially in the early stage of cement hydration [13]. However, the temperature gradient around the small pipes cannot be neglected. It generates remarkable thermal stress and increases the possibility of cracking [4]. Therefore, it is necessary to predict the temperature gradient and thermal field in concrete with cooling pipes.

How to improve efficiency in thermal field numerical simulation while maintaining accuracy? Various numerical methods have been developed to solve this problem, including the finite element method (FEM) [58], virtual element method (VEM) [9], composite element method (CEM) [1012], boundary element method (BEM) [13], and singular boundary method (SBM) [14]. Among the methods mentioned earlier, the finite element method is the most widely applied method to simulate the cooling system. It is generally employed to discretize the calculation domain into non-overlapping elements and take the polynomial basis functions as unknown fields. As a result, it is hard to exactly represent the small pipes in the large size of concrete. In addition, the dense grids in the FEM around the tiny pipes lead to a computationally intensive process and are time-consuming. Research workers create numerous improved approaches based on the FEM. Zuo [15] proposed an extended finite element method (XFEM) which is only used in the two-dimensional water pipe cooling problem. To avoid intensively meshing the small pipes, Kim [16] and Liu [17] chose to simply neglect the practical size and shape of the pipes by using a line element; as a consequence, the location of the pipes needed to leech on the node inconveniently. Liu considered the cooling pipe as a line element, and Chen [18] meshed the concrete by neglecting the location of the pipes and treated the pipes as an inner negative heat source after that. However, the accuracy of the temperature field still depends on the size of the element around the cooling pipe. To improve the effectiveness of calculations while ensuring accuracy, it is necessary to consider the typical temperature gradient distribution around the pipes and refine the grids selectively. In other words, dense elements are needed in regions with sharp temperature gradients, whereas sparse elements are suitable for other regions. However, the process of element partitioning is extremely time-consuming in methods like the FEM.

Instead of the FEM, this study introduces the isogeometric analysis (IGA) method to simulate the heat transfer process in mass concrete with cooling pipes. The IGA method, first introduced by Hughes [19], employs smooth and higher-order splines (such as non-uniform rational B-splines, T-splines [19], and LR-splines [20]) as basis functions to represent both the geometric space and approximation of solution fields. In the IGA context, the non-uniform rational B-splines (NURBSs) are most widely adopted to solve partial differential equations. There are several attractive characteristics [21] in NURBS-based IGA, such as great geometric accuracy, higher-order continuity, sparse matrix, and automatic mesh generation. Consequently, the IGA has been successfully applied in many research fields, for example, elasticity, shape optimization, crack propagation, and fluid flows [2228]. In addition, many research workers utilize the IGA in heat conduction areas. Attracted by the unique quality, R. Duvigneau [29] summarized the advantages of the IGA by comparing it with the traditional FEM through the thermal conduction test case. Chen [30] solved the Dirichlet boundary processing error caused by the lack of interpolation of the NURBS basis function. Yu [31] improved the accuracy of the steady-state heat conduction problem with locally refined adaptive IGA, and An [32] combined IGA and BEM to solve the same problem. To reduce the problem dimension and obtain highly accurate results for the flux values on the boundaries, Özgür applied the IGA-BEM to the steady-state heat conduction with volumetric heat source and non-linear boundary conditions [33].

To the best of the author’s knowledge, the IGA has not been used in the temperature simulation of mass concrete with cooling pipes. This work attempts to simulate the mass concrete temperature fields with the IGA method for the following reasons: first, the tiny pipes and large concrete can be accurately represented, even with a coarse grid. Second, the higher-order basis functions are more suitable for the non-linear temperature gradient than the traditional polynomial interpolation [34]. Third, it is convenient that there are no more efforts needed for mesh or pretreatment. Finally, the flexible refinement of IGA is appropriate for the different grid density requirements of complex temperature gradients and complicated surrounding conditions.

In this paper, IGA based on the NURBS is exploited to simulate the temperature field of two-dimensional mass concrete with one cooling pipe, as the preparation for cubic with multiple cooling pipes. The parameter space which describes small pipe and large concrete precisely is presented to generate an efficient model for calculation. The advantages of the suggested parameter space are shown by comparing it with others. Taking the correspondence between temperature gradient and distance from the pipe into account, the area in which temperature gradient changes dramatically around the cooling pipe is suggested, and the knot insertion vector for grid encryption is proposed. Furthermore, the situations with different calculating parameters such as water temperature, pipe thickness, and time step are investigated. The environmental conditions around concrete are also considered to demonstrate the flexibility of IGA.

The contents of this paper are as follows: Section 2 introduces the fundamental equations for transient heat conduction. Section 3 describes the IGA method and presents the parameter space utilized in mass concrete. Section 4 demonstrates the accuracy and performance of the simulation with several numerical examples. Finally, the conclusion is presented in Section 5.

2 Governing equation

The simplified mathematical model for the typical two-dimensional mass concrete with cooling pipes embedded is shown in Figure 1.

FIGURE 1

As can be seen in Figure 1A, there is a cooling pipe assumed as a negative heat resource [18] at the center of the concrete domain , and the inner boundary is in contact with the pipe. In addition, the boundary denotes the concrete surface in contact with air, and the boundary is treated as absolute heat insulation at the bottom. It is supposed that the concrete is isotropic and homogeneous, and the governing equation in this problem is given below:

In Eq. 1, represents the unknown temperature in any location at time, is the concrete thermal diffusivity, and are the specific heat capacity, density, and the heat conductivity, respectively. is the given function of adiabatic temperature rise. Moreover, the initial conditions and boundary conditions are shown in the following formulations.where the initial temperature of concrete is supposed to be a known function as Eq. 2 and represent the temperature of cooling water and air, respectively. Based on Chen [18], is the coefficient of equivalent heat convection between the water pipe and concrete, and are the thermal conductivity coefficient and thickness of the non-metallic water pipe, respectively. is the heat convection coefficient from the concrete to the air. It should be noted that only the third boundary condition (Robin boundary condition) is considered in this study.

According to the governing equation and boundary conditions mentioned in Eqs 35, the functional appeared by applying Galerkin’s method and integration by parts.

In Eq. 6, and . Boundary is the inner boundary in contact with pipe and boundary denotes the concrete surface in contact with air. The solution of the two-dimensional transient temperature field of the mass concrete with cooling pipe embedded is equivalent to the extreme value of the functional above.

3 Isogeometric analysis

In this section, we first give a brief introduction to NURBS as the foundation of the IGA. After that, IGA formulation for thermal field in mass concrete with cooling pipe is to be discussed. Finally, the most suitable IGA calculation model and refinement for numerical simulations are presented.

3.1 B-spline basis function and refinement

Different from the FEM, the IGA generates mesh by constructing parametric surface (volumes, in terms of three dimension) with CAD basis function like B-spline or NURBS. It provides some unique advantages to the isogeometric analysis, such as geometric exactness and simple refinements.

As NURBS is built from B-spline, it is necessary to introduce the definition and character of B-spline. The B-spline is defined in parametric space named knot vector , which is constructed with a non-decreasing set of knots , where and are the number of basis functions (as well as control points) and the polynomial order, respectively. The B-spline basis function could be defined recursively asand

Based on the B-spline basis function described in Eqs. 7, 8, the NURBS basis function is defined as follows:

Moreover, the basis function of the NURBS surface is given bywhere is the weight, as well as .

There are several crucial properties that emerged from the basis functions above. The first one is that each basis function is non-negativity and constitutes partition of unity. The second main property to note is compact support, that is, mentioned in Eq. 9 is non-zero only in the knot span . This means that the stiffness matrix in numerical calculation is a large sparse matrix, which is similar to their FEM counterparts. Finally, every order basis function has continuous derivatives at knot which is repeated times in the knot span. As a result, the continuity at any knot value could be selected by controlling the number of the knot value. In addition, the basis is interpolatory at a knot when the multiplicity is . It is worth noting that this is the most distinctive property of the IGA.

Given NURBS basis functions, , and corresponding weighted control points , a piecewise-polynomial NURBS curve, is obtained by

Similarly, a tensor product NURBS surface is defined by

As shown in Eqs 11, 12, the NURBS curve (surface) is approximated by a control polygon consisting of control points and the control polygon is not unique. As a result, the choice of control points should be based on the geometric features of the objective curve like endpoints and inflection points. The weights at control points reflect the distance from the curve to the control polygon.

Similar to the FEM, appropriate refinement should be taken to ensure the computational accuracy; however, a coarse mesh could guarantee the geometric exactness in the isogeometric approach. Three main means are used frequently to enrich the basis function, namely, knot insertion, order elevation, and k-refinement. Here, a brief introduction to the approach is given.

The first means is knot insertion by which a new knot is inserted in a given one-dimensional knot vector Ξ and a new extended knot vector is created. It will increase the number of basis functions to , as well as the control points, and the mesh will be tessellated. The new weighted control points can obtained by Eqs 13, 14where

Unfortunately, the continuity of the curve based on is increased to at the new knot. To keep the continuity of knot unchanged, the multiplicity of the new knot should stay consistent with other knots, that is, the second means, k-refinement. The last one is order elevation, which is much similar to the p-refinement used in the FEM. It has no influence on the existing mesh, but the number of basis functions and control points in each mesh will be increased.

3.2 IGA formulation for thermal field

In the IGA method, the tensor products of the NURBS basis function are defined in the parametric space as Eq. 10. To apply the IGA, we first introduce the transformation between parametric space and physical space , which can be written as follows

In the non-zero mesh , based on Eq. 15, the approximation of the transient temperature field for this problem can be expressed aswhere is the temperature of control point at time , is the element basis function matrix in the mesh, and is the element control point temperature matrix. Substituting Eq. 16 into the function (6) and using the variational method to make , we can obtainwhere is the element thermal capacity matrix, is the element thermal transfer matrix, and is the element temperature load matrix. The element boundaries are in contact with the pipe and the air. is the computational element domain in . are the Jacobian determinants for transformation . Benefiting from the local support of the NURBS, the matrix mentioned above is a sparse matrix which is similar to the FEM. Based on Eqs 1720, we assemble the element matrices to obtain the global equation.where

To discretize the time, the difference method is employed in Eqs 21, 22 . It is supposed that the temperature changes linearly at each time step after the time domain has been divided into time steps. As the initial temperature of concrete is given, we assume that temperature is known before the calculation of . Then, we obtainand

As the solution based on the backward difference method is unconditionally convergent, we take into the equation above and then substituting Eqs 23, 24 into Eq. 21, we can obtain

The transient temperature field of mass concrete with the cooling pipe embedded is obtained through Eq. 25.

3.3 Mathematical model in IGA

The different size between the large scale of concrete and the small diameter of water pipe makes it hard to ensure the geometric accuracy in this problem, which greatly affects the results of numerical simulation. In the IGA, the parameter space is able to depict geometry exactly at all levels of discretization. As a result, this method is expected to acquire improved accuracy as compared with less geometrically precise methods. Moreover, the most time-consuming part is bypassed as the process of discretization is accomplished when the parameter space is determined. However, the parameter space is non-unique to model the same physical space, and the different parameter space has a significant impact on the accuracy and efficiency. This makes the suitable choice of parameter space important in the IGA.

The strategy of this study is that minimal control points are taken to represent geometry in detail and the refinement based on the temperature gradient in the mass concrete is implemented on those control points. The typical mathematical model for the two-dimensional mass concrete can be seen in Figure 1A, and it also can be assumed as an infinite plate with a circular hole. For simplicity of the expression, the parameter space is only constructed with a quarter of the plate. To depict the circular hole as water pipe, the polynomial order in direction is at least quadratic, and linear is the minimum order in direction. At least four control points are required in direction for this mathematical model, and there is a repeated control point in the upper left corner which will generate inhomogeneous mesh and make the element abnormal after refinement although Hughes [35] has given the advice solution. The coarsest mesh with four control points in direction is illustrated in Figures 2A, B, in which the control point is denoted as the green one. The corresponding inhomogeneous encryption mesh can also be seen under the arrows.

FIGURE 2

For the purpose of a well-proportioned mesh, the recommended parameter space is given in Table 1, in which five control points are adopted in direction. As shown in Figure 2, the performance of the three parameter spaces is identical in the coarsest mesh but is entirely different after refinement, and the harmonious refinement is available only in the suggested parameter space as exemplified in Figure 2C. Conveniently, the proposed parameter space is appropriate for the different diameter of water pipe , as well as concrete size . More details are available in Table 1.

TABLE 1

Control point

Recommended parameter space and control points.

The mechanisms of NURBS refinement also play an important part in the construction of parameter space. J.A. Cottrell [34] found that the k-refinement method achieves a significant improvement in accuracy in the problems of structural vibrations over the classical -continuous p-refinement method. R. Duvigneau [29] proved that quadratic functions result in a significant error reduction compared with the linear functions in the heat conduction problem with the same degree of freedom, whereas less accuracy improved compared with the cubic functions. However, the CPU time required for quadratic functions is considerably reduced compared to the cubic functions. Chen [30] figured out that the quadratic functions in the IGA have almost an equal convergence rate with the cubic Lagrange basis function in transient heat conduction of solid medium. To further improve the computational efficiency and accuracy of the IGA method, the order elevation is discarded for refinement in this paper, and the quadratic function is adopted in both and directions. It not only increases the accuracy of the recommended parameter space but also guarantees the continuity of basis function in the concrete domain. In addition, the knot insertion generates elements with an inhomogeneous distribution of control points although the grid is the same as k-refinement. In summary, the k-refinement is considered to be the most appropriate for this problem among the three main kinds of mechanisms mentioned in Section 3.1. To ensure the continuity in the knot and uniform distribution of the control point, new knots should be inserted with the multiplicity of polynomial order. Argumentation is given in example 1.

At last, the new knot vector should be established flexibly according to the temperature gradient in the practical condition. It is assumed that there is a cylindrical region named A around the water pipe in which the temperature gradient is perpendicular to the cooling pipe surface and has no relationship with the boundary conditions. This region is time-varying and the higher temperature gradient appears near the water pipe. In consequence, the dense elements are needed in region A and the parameter space is only refined in direction, whereas sparse elements are arranged in region B. The choice of the cylindrical region A is crucial. The CEM [12] roughly took the cylindrical region with a radius of 0.12 m from the water pipe as the large gradient regions according to engineering experience regardless of the different material parameters. The radial basis function finite difference method (RBF-FD) [36] only discussed the ARMS error with different radius values of the region. Zhu [37] has figured out the relationship between temperature gradient and distance from the pipe in time with step . The formulation is given below:In the Eq. 26, represents the heat transfer from the water pipe to the concrete through the inner face in contact with the pipe, and it can be found in Eq. 28. is the equivalent temperature gradient at the inner face and is given as Eq. 27. is the temperature of the concrete contacted with the pipe and is the external radius of the pipe. is the average temperature increment of the concrete between and . is the length of a section of the water pipe. are the given hydration heat parameters of concrete. The schematic plot and relationship are shown in Figure 3. It can be seen that the temperature gradient varies dramatically from to 0.2 , but slowly at remaining areas. Consequently, radius of the region to be refined can be calculated from selected , and the different element densities are adopted to the corresponding regions after that. Based on the method mentioned above, the new knot vector and proposed parameter space are constructed and the mathematical model in the IGA is established.

FIGURE 3

4 Numerical experiments

In this section, three examples are given to examine the accuracy, stability, and sensibility of the proposed parameter space in the IGA. To evaluate the numerical results, the relative error (RE) and the relative root mean square error (RMSE) are applied, and the definition is expressed as Eqs 29, 30.where denotes the exact solution based on the model as shown in Figure 1B in the FEM at point . is the numerical solution of the IGA. is the total number of test points.

4.1 Example 1

There is a two-dimensional square concrete with length and a cooling pipe in the center of the concrete with external radius . The thickness of the pipe is . The initial temperature of concrete and water are and , respectively. The main focus of this example is to examine the validity of the recommending parameter space and show guidance of the selection refinement, so only adiabatic boundary conditions are considered. The whole process lasts 20 days with a time step of 0.02 days. The material parameters are shown in Table 2, and the adiabatic temperature increase of the concrete is .

TABLE 2

Material parameterConcreteWaterCooling pipe
Density ()24201000
Specific heat capacity 0.894.187
Thermal conductivity coefficient 9.01.66

Material parameters of concrete and water pipe, for example 1.

Based on the formulation mentioned in part 3.3, we take as region A at the beginning, and the grid size of the region A is . The domain of the concrete is divided into four elements in direction, and each size of the elements is , as the initial model shown in Figure 4A. To illustrate the rationality of the parameter space and the need for refinement, the model with repeated control points and uniform refinement of the parameter space is displayed as the comparison in Figure 4B. The model of parameter space with recommended control points and uniform refinement is displayed in Figure 4C as well. There are three test points selected with the coordinates of P1 (−0.4, 0), P2 (−0.8, 0), and P3 (−0.8, 0.8). It is worth noting that the origin of coordinates is the center of the water pipe in this paper.

FIGURE 4

Table 3 shows the numerical solutions of different models with almost the same number of degrees of freedom in the IGA and the exact solution in the FEM with 1256 numbers of nodes. The temperature of the three test points at different times is given. By comparing the RMSE of the first two models, it is clear that the repeated control points result in malformed meshes that lead to low accuracy. From the results obtained by the two models with the same number of degrees of freedom, it can be found that encryption based on temperature gradients can increase the accuracy with any additional effort. Moreover, the good accuracy of the proposed parameter space can be found and in that the RMSE can reach up to 8.22E-3, even though only 297 control points are used in the IGA. These results suggest that the construction of the proposed parameter space is reasonable and necessary in this paper.

TABLE 3

Temperature (°C)Degree of freedomP1P2P3RMSE
t = 2.26 dayt = 10.88 dayt = 4.20 dayt = 15.02 dayt = 16.36 day
Model as Figure 4B26120.012118.137022.711318.423617.95241.98E-1
Model as Figure 4C29724.968622.757126.660321.637021.40353.34E-2
Initial model as Figure 4A29725.232723.475926.944122.415122.17998.22E-3
FEM125625.287023.700227.034122.673222.4473

RMSE of numerical solution varies with the model under different test points and time, for example 1.

To verify the convergence of the initial model and the rationality of encryption approach selection, different means of refinement established on the initial model are displayed in Figures 4D–F. For concise expression, we named k-refinement in direction for one time as , and knot insertion (similar to the h-refinement used in the FEM) in direction for two times as , and so on. The grid of knot insertion is not given for the reason that it has totally the same grid as k-refinement. As the refinement of the initial model is already based on temperature gradient, the subsequent encryption on the initial model adopts uniform encryption. Three test points are chosen with the coordinates of Q1 (−0.08, 0), Q2(−0.2, 0), and Q3 (−0.4, 0).

Figure 5 shows the relative error of the time variable for different refinement means at the given test points. In addition, the total time consumed is provided in Table 4. Incidentally, these tests are performed using an Intel Core i5 2.4 GHz. It can be seen that the coarse grid applied in the initial model results in a relative error increment and finally reaches a 1.6E-2 error. After k-refinement in direction, the accuracy is significantly improved and achieves a relative error of 4E-3, which is supposedly sufficient for practical engineering. In contrast to , the error is only unnoticeably decreased in , but the time consumed is doubled. Compared with , has a similar performance in not only the accuracy but also the time consumption. Furthermore, the same phenomenon occurs between and . From the consequences of , it has been proven that the temperature gradient changes in pipe radial and refinement in direction is totally ineffective in this case. It is worth noting that the time consumption is dramatically reduced compared with the FEM in almost the same degree of freedom.

FIGURE 5

TABLE 4

InitialFEM
Degree of freedom297561108911054296931256
Total time (s)7.2215.1931.7832.7913.5629.29129.63

Total time and numbers of freedom used in different models, for example 1.

From the discussion aforementioned, the parameter space for is supposed to be the suitable refinement in this condition, and it might be useful to determine the mechanism in other conditions. The temperature of these test points in the time domain is shown in Figure 6A, and it is obvious that the transient temperature and temperature variation regularity of the proposed scheme are identical to those of the FEM. In addition, the relative error of the test point is shown in Figure 6B. It can be seen that the relative error is nearly 3E-3 and converges with time finally. All the results demonstrate the accuracy and efficiency of the proposed scheme.

FIGURE 6

4.2 Example 2

The main focus of this example is to figure out the sensitivity of the IGA to thickness, water temperature, and time step. Concrete length and pipe external radius are and , respectively. The initial temperatures of concrete and cooling water are and , respectively. The adiabatic temperature increase of the concrete is , and the thickness of the water pipe is . The remaining parameters are not mentioned as they are the same as in example 1. The mechanism of refinement is referred to applied in example 1. The refinement interval used in region A is , whereas incremental intervals , , and are used in region B. The FEM model with 1792 elements and 1888 nodes is treated as the exact solution, and the points R1 (−0.08, 0), R2 (−0.3, 0), R3 (−0.5, 0), and R4 (−1, 0) are selected as the test points.

To display the influence of water pipe thickness, the numerical results with different thicknesses of the proposed scheme and the FEM in the test points are enumerated in Table 5. It can be found that the thickness has a considerable influence on the temperature, and the effect of the cooling pipe decreases as the thickness increases. The RMSE converges to a level around 1E-3 although the heat convection coefficient for the surface contacted with water pipe is crudely represented by the coefficient of equivalent heat convection.

TABLE 5

Thickness ()MethodR1R2R3R4RMSE (E)
t = 1.96 dayt = 3.54 dayt = 13.68 dayt = 15.82 dayt = 7.42 day
0.002IGA17.939723.594122.697224.090127.30572.11–3
FEM18.001923.626822.756124.144127.3349
0.004IGA18.913724.136123.782025.084127.76493.90–3
FEM19.023124.195823.901425.193227.8172
0.006IGA19.665624.558224.675825.900928.12574.68–3
FEM19.794624.629524.827126.038928.1882

Comparisons between temperature obtained from the proposed scheme and the FEM in test points with different thicknesses and time, for example 2.

Table 6 lists the comparison between the temperature obtained from the IGA with different water temperatures and those obtained from the FEM at the test points. It can be seen that colder water reduces the temperature observably at each test point and time. However, the increased temperature difference between water and concrete produces large temperature gradient which results in the increase in RMSE and the risk of cracking.

TABLE 6

Water temperature ()MethodR1R2R3R4RMSE (E)
t = 2.96 dayt = 4.54 dayt = 14.68 dayt = 16.82 dayt = 8.42 day
5IGA15.789222.458819.788221.308025.92243.09–3
FEM15.877322.505519.862121.376425.9641
10IGA18.427723.938822.396423.725227.17252.25–3
FEM18.496023.975222.456823.780927.2053
15IGA21.066325.418725.004626.142428.42271.54–3
FEM21.114825.444925.051526.185428.4465

Temperature and RMSE in test points with different water temperatures and times for IGA and FEM, for example 2.

Table 7 displays the temperature with different time steps by the proposed scheme compared with the FEM. The stability of time step for the proposed scheme has been proven from the results above. With increasing time steps, the temperature increases slightly in both the IGA and FEM, and the RMSE remains unchanged. From the numerical solutions listed in the three tables above, the accuracy and stability of the given model in IGA have been demonstrated.

TABLE 7

Time step (day)MethodR1R2R3R4RMSE (E)
t = 1.6 dayt = 6.4 dayt = 8.4 dayt = 14.4 dayt = 16.8 day
0.02IGA17.689024.137923.970124.598424.47562.23–3
FEM17.748324.180624.018424.649624.5306
0.16IGA17.693324.130423.962824.593024.47132.23–3
FEM17.752624.173124.011124.644224.5263
0.4IGA17.725524.096123.929324.568824.45302.23–3
FEM17.784924.138823.977624.619824.5078

Temperature and RMSE in those test points mentioned above with different time steps and times for IGA and FEM, for example 2.

4.3 Example 3

In the last example, the temperature field for the proposed scheme is illustrated comprehensively in the concrete domain instead of the test points, and the boundary in contact with air is considered as a simulation of the practical engineering situation. Heat convection coefficient for boundary is , and the temperature is . The concrete bottom is on the fundamental which is considered as the insulation panel. The material parameters are identical to those in example 1. In the previous examples, the mechanisms of refinement are only used in direction, as the other boundary conditions are adiabatic. To take boundary into account, uniform k-refinement is applied in direction for one time on the basis of model employed in example 1, and grids near the insulation panel remain unchanged to reduce the amount of computation.

To further verify the accuracy of the IGA, Figure 7 presents the temperature field at several time histories by the IGA compared with the FEM. We can see that the distribution of the temperature field is almost the same in the IGA and FEM. At the beginning of the experiment, the temperature of the concrete around the surface increases although the center is affected by the cooling pipe. With the heat of hydration proceeding, the temperature in most areas of the concrete is increased remarkably, and the temperature gradient around the pipe is changed sharply. When it comes to the time of 20 days, the temperature around the cooling pipe decreases after the peak of the hydration reaction. Despite the process being complicated, a good consistency can be observed between the temperature field of the offered scheme and the FEM. The slight difference in the bottom of the concrete is caused by the coarse grids used in boundary .

FIGURE 7

Figure 8 displays the absolute relative error correspondingly. It can be found that the relative error steadily converges to a level around E-3 in most regions over several time histories although the largest error near 4E-2 is concentrated around the cooling pipe. This phenomenon is mostly caused by the approximate representation of the cooling pipe effect and the large temperature gradient around the pipe. All the results in this example demonstrate the accuracy of the IGA in the entire concrete, as well as the flexibility to refine according to the different boundaries.

FIGURE 8

5 Conclusion

In this paper, the IGA method based on NURBS is introduced to simulate the transient temperature field in massive concrete with cooling pipe embedded. As a preparation for the multi-pipes in three-dimensional simulation, only a classical model of simplified two-dimensional concrete with one pipe is considered, where the cooling pipe is considered as a negative heat resource.

Based on the aforementioned classical model, the most concise parameter space that can precisely describe the simplified model is provided. Different mechanisms for encrypting parameter space are discussed, and the unique k-refinement has been proven to be the most efficient method. To improve the efficiency under the premise of ensuring accuracy, the characteristics of temperature gradient in this problem is considered and different levels of refinement are applied in constructing the parameter space. Furthermore, the parameter space also can be refined flexibly according to the actual boundary conditions, which is expected to reduce the time consumption. The construction of parameter space is not only to avoid the cumbersome and useless refinement but also to leap over the time-consuming process of the mesh.

Numerical examples are presented to verify the validity and applicability of the proposed scheme. The rationality of the constructed parameter space and the accuracy of the IGA are demonstrated through numerical results displayed in the first example. The stability for different calculating parameters is also indicated in the second example. In addition, the flexibility and applicability of the method are shown in the last example. As a result, the proposed parameter space in the IGA method is deemed to possess great potential in simulating the thermal field in mass concrete with a cooling pipe system.

Statements

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material; further inquiries can be directed to the corresponding author.

Author contributions

QL: writing–original draft. GC: writing–review and editing. FZ: writing–review and editing and resources.

Funding

The author(s) declare that no financial support was received for the research, authorship, and/or publication of this article.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The reviewer DL declared a shared affiliation with the author(s) to the handling editor at the time of review.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

  • 1.

    KwakYHWalewskiJSleeperDSadatsafaviH. What can we learn from the Hoover Dam project that influenced modern project management. Int J Project Manage (2014) 32(2):25664. 10.1016/j.ijproman.2013.04.002

  • 2.

    SeoT-SKimS-SLimC-K. Experimental study on hydration heat control of mass concrete by vertical pipe cooling method. J Asian Architecture Building Eng (2015) 14(3):65762. 10.3130/jaabe.14.657

  • 3.

    TasriASusilawatiA. Corrigendum to “Effect of cooling water temperature and space between cooling pipes of post-cooling system on temperature and thermal stress in mass concrete” [J. Build. Eng. 24 (2019) 100731]. J Building Eng (2019) 26:100922. 10.1016/j.jobe.2019.100922

  • 4.

    YuXChenJXuQZhouZ. Research on the influence factors of thermal cracking in mass concrete by model experiments. KSCE J Civil Eng (2017) 22(8):290615. 10.1007/s12205-017-2711-2

  • 5.

    ChengJLiTCLiuXZhaoLH. A 3D discrete FEM iterative algorithm for solving the water pipe cooling problems of massive concrete structures. Int J Numer Anal Methods Geomechanics (2016) 40(4):487508. 10.1002/nag.2409

  • 6.

    ZhuB. The equivalent heat equation of pipe cooling in mass concrete considering influence of external temperature. J Hydraulic Eng (2003) 3(3):4954. 10.3321/j.issn:0559-9350.2003.03.010

  • 7.

    ZhuYXuZHeJet alA calculation method for solving temperature field of mass concrete with cooling pipes. J Yangtze River Scientific Res Inst (2003) 20(2):1922. 10.1063/1.4982444

  • 8.

    ConceiçãOJFariaRAzenhaMet alA new method based on equivalent surfaces for simulation of the post-cooling in concrete arch dams during construction. Eng Structures (2020) 209. 10.1016/j.engstruct.2019.109976

  • 9.

    VaccaGAn H1-conforming virtual element for Darcy and Brinkman equations. Math Models Methods Appl Sci (2017) 28(01):15994. 10.1142/s0218202518500057

  • 10.

    ChenSHSuPShahrourI. Composite element algorithm for the thermal analysis of mass concrete. Int J Numer Methods Heat Fluid Flow (2011) 21(4):43447. 10.1108/09615531111123100

  • 11.

    ZhongRHouG-PQiangS. An improved composite element method for the simulation of temperature field in massive concrete with embedded cooling pipe. Appl Therm Eng (2017) 124:140917. 10.1016/j.applthermaleng.2017.06.124

  • 12.

    ShengQChaoWZhenyangZ. A new composite element algorithm for the temperature filed of concrete containing cooling pipe. SHUILI XUEBAO (2015) 46(6):73948. 10.13243/j.cnki.slxb.20140910

  • 13.

    SimõESISimõESNTadeuAReisMVasconcellosCMansurW. Experimental validation of a frequency domain BEM model to study 2D and 3D heat transfer by conduction. Eng Anal Boundary Elem (2012) 36(11):168698. 10.1016/j.enganabound.2012.05.007

  • 14.

    WeiXChenBChenSYinS. An ACA-SBM for some 2D steady-state heat conduction problems. Eng Anal Boundary Elem (2016) 71:10111. 10.1016/j.enganabound.2016.07.012

  • 15.

    ZuoZHuYLiQLiuG. An extended finite element method for pipe-embedded plane thermal analysis. Finite Elem Anal Des (2015) 102-103:5264. 10.1016/j.finel.2015.05.002

  • 16.

    KimJKKimKHYangJK. Thermal analysis of hydration heat in concrete structures with pipe-cooling system. Comput Structures (2001) 79:16371. 10.1016/s0045-7949(00)00128-0

  • 17.

    LiuXZhangCChangXZhouWChengYDuanY. Precise simulation analysis of the thermal field in mass concrete with a pipe water cooling system. Appl Therm Eng (2015) 78:44959. 10.1016/j.applthermaleng.2014.12.050

  • 18.

    GuorongCWentaoXYunYet alComputation method for temperature field of mass concrete containing cooling water pipes. Chin J Comput Phys (2012) 29:4116. 10.19596/j.cnki.1001-246x.2012.03.015

  • 19.

    BazilevsYCaloVMCottrellJAEvansJHughesTLiptonSet alIsogeometric analysis using T-splines. Comput Methods Appl Mech Eng (2010) 199(5-8):22963. 10.1016/j.cma.2009.02.036

  • 20.

    DokkenTLycheTPettersenKF. Polynomial splines over locally refined box-partitions. Comput Aided Geometric Des (2013) 30(3):33156. 10.1016/j.cagd.2012.12.005

  • 21.

    BazilevsYDa VeigaLCottrellJet alIsogeometric analysis: approximation, stability and error estimates for h-refined meshes. Math Models Methods Appl Sci (2006) 16(7):103190. 10.1142/S0218202506001455

  • 22.

    YoonMHaS-HChoS. Isogeometric shape design optimization of heat conduction problems. Int J Heat Mass Transfer (2013) 62:27285. 10.1016/j.ijheatmasstransfer.2013.02.077

  • 23.

    GillebaartEDe BreukerR. Low-fidelity 2D isogeometric aeroelastic analysis and optimization method with application to a morphing airfoil. Comput Methods Appl Mech Eng (2016) 305:51236. 10.1016/j.cma.2016.03.014

  • 24.

    FischerPKlassenMMergheimJSteinmannPMüllerR. Isogeometric analysis of 2D gradient elasticity. Comput Mech (2010) 47(3):32534. 10.1007/s00466-010-0543-8

  • 25.

    AkkermanIBazilevsYKeesCEFarthingM. Isogeometric analysis of free-surface flow. J Comput Phys (2011) 230(11):413752. 10.1016/j.jcp.2010.11.044

  • 26.

    BhardwajGSinghIVMishraBKBuiT. Numerical simulation of functionally graded cracked plates using NURBS based XIGA under different loads and boundary conditions. Compos Structures (2015) 126:34759. 10.1016/j.compstruct.2015.02.066

  • 27.

    BuffaASangalliGVazquezR. Isogeometric analysis in electromagnetics_ B-splines approximation. Comput Methods Appl Mech Eng (2010) 199(17–20):114352. 10.1016/J.CMA.2009.12.002

  • 28.

    WangJZhangHZhengYYeH. High-order NURBS elements based isogeometric formulation for swellable soft materials. Comput Methods Appl Mech Eng (2020) 363:112901. 10.1016/j.cma.2020.112901

  • 29.

    DuvigneauRAn introduction to isogeometric analysis with application to thermal conduction. INSTITUT NATIONAL DE RECHERCHE EN INFORMATIQUE ET EN AUTOMATIQUE (2012).

  • 30.

    ChenT. Imposing dirichlet boundary conditions with point collocation method in isogeometric analysis. J Mech Eng (2012) 48(5):15764. 10.3901/jme.2012.05.157

  • 31.

    YuTChenBNatarajanSBuiTQ. A locally refined adaptive isogeometric analysis for steady-state heat conduction problems. Eng Anal Boundary Elem (2020) 117:11931. 10.1016/j.enganabound.2020.05.005

  • 32.

    AnZYuTBuiTQWangCTrinhNA. Implementation of isogeometric boundary element method for 2-D steady heat transfer analysis. Adv Eng Softw (2018) 116:3649. 10.1016/j.advengsoft.2017.11.008

  • 33.

    GüMüŞÖCBaranoğluBÇetinB. Isogeometric and NURBS-enhanced boundary element formulations for steady-state heat conduction with volumetric heat source and nonlinear boundary conditions. Eng Anal Boundary Elem (2022) 145:299309. 10.1016/j.enganabound.2022.09.021

  • 34.

    CottrellJAHughesTJRRealiA. Studies of refinement and continuity in isogeometric structural analysis. Comput Methods Appl Mech Eng (2007) 196(41-44):416083. 10.1016/j.cma.2007.04.007

  • 35.

    HughesTJRCottrellJABazilevsY. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. Comput Methods Appl Mech Eng (2005) 194(39-41):413595. 10.1016/j.cma.2004.10.008

  • 36.

    HongYLinJChengAHDWangYChenW. Thermal analysis of heat transfer in pipe cooling concrete structure by a meshless RBF-FD method combined with an indirect model. Int J Therm Sci (2020) 152:106296. 10.1016/j.ijthermalsci.2020.106296

  • 37.

    ZhuZQiangSChenW. A new method solving the temperature field of concrete around cooling pipes. Comput Concrete (2013) 11(5):44162. 10.12989/cac.2013.11.5.441

Summary

Keywords

isogeometric analysis, mass concrete, cooling pipe system, temperature field, heat transfer, temperature gradient

Citation

Li Q, Chen G and Zhu F (2024) Simulation of thermal field in mass concrete with cooling pipes based on the isogeometric analysis method. Front. Phys. 12:1338718. doi: 10.3389/fphy.2024.1338718

Received

15 November 2023

Accepted

22 April 2024

Published

24 May 2024

Volume

12 - 2024

Edited by

Teoman Özer, Istanbul Technical University, Türkiye

Reviewed by

Dong Lei, Hohai University, China

Fuzheng Gao, Shandong University, China

Ziqi Yu, Toyota Research Institute of North America, United States

Updates

Copyright

*Correspondence: Qingwen Li,

Disclaimer

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.

Outline

Figures

Cite article

Copy to clipboard


Export citation file


Share article

Article metrics