Thermal analysis of mass concrete structures with cooling pipes by an adaptive moving Kriging method

Abstract

Numerical simulation of the temperature field in massive concrete with embedded cooling water pipes is of great engineering significance for temperature control, crack prevention, and optimization of cooling strategies. To obtain an accurate temperature field, conventional finite element approaches often require local mesh refinement in regions with steep thermal gradients near the pipe wall, which makes model construction cumbersome. To overcome the dependence on mesh generation, this study couples a local weak-form Moving Kriging (MK) interpolation method with a moving mesh partial differential equation (MMPDE) approach for the thermal analysis of pipe cooled concrete. By exploiting the smoothness and differentiability of the Gaussian basis function, the proposed MK approximation is able to represent the continuous temperature field with sharp gradients in the vicinity of cooling pipes. Meanwhile, the MMPDE method constructs a metric tensor from the Hessian of the temperature field, thereby automatically driving nodes to concentrate in regions where the numerical solution varies rapidly. As a result, strong-gradient-induced numerical oscillations and excessive errors can be effectively suppressed without a significant increase in the degrees of freedom. Numerical results demonstrate that, compared with conventional fixed node distributions, the proposed algorithm achieves higher accuracy and improved stability.

1 Introduction

Mass concrete structures are among the most common forms of civil infrastructure and are widely used in hydraulic dams, large bridges, and underground constructions. During hardening, the substantial heat released by cement hydration induces pronounced internal temperature differences and associated thermal stresses [1, 2], which constitute a primary cause of temperature-induced cracking. Once such cracks develop, they not only jeopardize the long-term durability of the structure but can also markedly reduce its load-carrying capacity and may even trigger sudden brittle failure, leading to losses in life and property [3]. In current engineering practice, embedded cooling water pipes are frequently adopted to mitigate thermal stresses in mass concrete. However, an improperly designed cooling scheme can readily generate extremely steep temperature gradients in the vicinity of the cooling pipes, resulting in localized internal cracking or even through-cracks, thereby compromising structural safety [4–6]. Therefore, accurate prediction and effective regulation of the temperature evolution in mass concrete are fundamental to preventing detrimental thermal cracking and ensuring structural integrity.

At present, temperature-field simulations for concrete containing cooling pipes in practical engineering design are predominantly based on equivalent approaches. In such methods, the cooling pipe is represented as a “negative heat source,” and only the average cooling effect is considered [7]. Although this strategy is computationally simple and efficient, it neglects the intense temperature variations near the pipe wall and thus cannot accurately resolve the local temperature field around the pipe [8–11]. To achieve a more faithful description of the near-pipe thermal behavior, engineering analyses commonly employ element-based numerical methods, among which the finite element method (FEM) is the most widely used [12–16, 46]. Owing to its rigorous theoretical foundation, FEM solutions are generally reliable; moreover, the availability of mature FEM software packages has substantially facilitated routine implementation in engineering computations.

Two major challenges arise when FEM is used to analyze pipe-cooled concrete. First, obtaining accurate numerical solutions typically requires local mesh refinement in the vicinity of the cooling pipes; for large-scale hydraulic projects, this leads to enormous computational and storage demands [17]. Second, FEM meshing is time-consuming and often requires substantial manual intervention, sometimes accounting for more than 90% of the overall analysis workflow [18, 19]. Therefore, meshfree methods, which only require nodes to be placed in the computational domain and on its boundary without connecting the nodes through fixed elements, have attracted increasing attention from researchers. Compared with FEM, meshfree methods may simplify preprocessing and are often more amenable to automation in problems where repeated local remeshing is required; they have therefore been widely applied across diverse fields [20–23].

Smoothed particle hydrodynamics (SPH) discretizes a continuum into a set of particles carrying physical quantities such as mass and energy [24]. Within a Lagrangian framework, the macroscopic behavior of the continuum is described through particle–particle interactions; SPH is particularly suited to problems involving large deformations, free surfaces, moving boundaries, and complex multiphase flows. For quasi-static problems, crack propagation, and solid mechanics analyses with complex geometries, the element-free Galerkin (EFG) method based on moving least squares (MLS) approximation has been extensively employed [25]. The Moving Kriging (MK) method constructs shape functions from Kriging interpolation theory, where the key principle is to satisfy the optimality conditions of unbiased estimation and minimum variance, thereby enabling high-accuracy approximation of physical fields [26]. This property is advantageous for complex fields exhibiting sharp variations, such as the steep temperature gradients near embedded cooling pipes. Moreover, owing to the Kronecker-delta property, essential boundary conditions can be imposed in MK in a manner similar to FEM, avoiding special boundary treatments and simplifying the imposition of essential boundary conditions.

Wang et al. [27] introduced Moving Kriging interpolation into a meshfree Galerkin scaled boundary framework and achieved improved accuracy and convergence for steady-state heat conduction. Chen et al. [28] combined MK with a local Petrov–Galerkin formulation to solve transient heat conduction problems. Despite its favorable accuracy, MK still requires a large number of nodes around each cooling pipe to resolve the high temperature gradients. This often results in ill-conditioned system matrices, reduced computational efficiency, and numerical instability. Consequently, adaptive node redistribution provides an effective avenue for improvement. Existing adaptive strategies can be broadly classified into three categories: (i) h-adaptivity, which locally increases or decreases the number of nodes based on an error indicator [29, 30]; (ii) p-adaptivity, which raises or lowers the order of the approximation functions in regions requiring refinement [31]; and (iii) r-adaptivity, which relocates nodes according to solution features [32]. In this study, an r-type moving mesh partial differential equation (MMPDE) approach [33, 34] is introduced to alleviate the loss of accuracy and deterioration of stability that may occur when nodes are locally concentrated.

This study develops a hybrid algorithm that couples MK interpolation with the MMPDE method, denoted as MM–MK, to accurately simulate the temperature field in mass concrete containing embedded cooling water pipes. The proposed strategy decouples adaptive node redistribution from the physical-field solution: MMPDE acts as an independent dynamic node-placement module that generates an optimally uniform node distribution in the computational domain under a prescribed tensor metric by solving an optimization problem during time stepping. Based on the resulting node set, the temperature field is then analyzed using MK interpolation, which possesses the Kronecker-delta property. Within this framework, the mesh produced by MMPDE serves only as a coordinate-mapping tool that links the physical and computational domains, rather than directly participating in the discretization of the governing equations. The effectiveness of the proposed method is verified through numerical experiments on three representative benchmark problems.

The remainder of this paper is organized as follows. Section 2 presents the governing equations for concrete with embedded cooling pipes. Section 3 describes the core numerical methodology that integrates MK interpolation with the MMPDE technique. Section 4 provides comprehensive numerical investigations and comparative discussions to validate the proposed framework. Finally, Section 5 summarizes the main findings of this study.

2 Mathematical model

As shown in Figure 1, to investigate the cooling effect of embedded water pipes in mass concrete, we consider a model in which a cooling pipe passes through the concrete body. Cooling water circulates inside the pipe and removes the heat released during cement hydration, thereby suppressing the peak temperature and mitigating the risk of temperature-induced cracking. Based on this physical process, and neglecting heat transfer along the pipe-axis direction, this section derives the governing heat-conduction equation for mass concrete with embedded cooling pipes. The concrete is assumed to be isotropic.

FIGURE 1

The transient heat-conduction equation of the model is given bywhere denotes the temperature at a point in the concrete domain , as a function of spatial coordinates and time . Here, is the adiabatic temperature-rise rate of the concrete, is the adiabatic temperature rise, and is the thermal diffusivity of the concrete.

The initial condition of the concrete is prescribed as a known temperature . The boundaries corresponding to the top surface, bottom surface, left surface, right surface, and the pipe wall of the concrete are denoted by , , , , and , respectively. To solve the governing equation given in Equations 1, the initial condition and boundary conditions of the concrete domain must be specified, as given in Equations 2–5.

The initial condition isthe Dirichlet boundary condition isthe Neumann boundary condition isand under the Robin boundary condition,where is the convective heat-transfer coefficient between the concrete surface and the ambient air, is the ambient temperature, and is the thermal conductivity of concrete.

3 Numerical scheme

3.1 Moving Kriging interpolation

Let denote a physical field defined on the problem domain . The domain is discretized by field nodes. For an arbitrary point , suppose that its local support domain contains field nodes, and the corresponding nodal values are . Then, for any evaluation point in the computational domain. The MK interpolation approximation, the associated shape-function construction, and the derivative expressions are given in Equations 6–13.where are the MK shape functions, , and .

The shape-function vector constructed by the MK interpolation is given byThe matrices and are defined aswhere is a set of polynomial basis functions, is the number of polynomial terms, and is the identity matrix. In this study, linear polynomials are employed as the basis functions; for the two-dimensional case,Accordingly, the matrices and can be written as

The correlation vector is given bywhere is the correlation function between the point and node ; its choice directly affects the accuracy and smoothness of the interpolation. Common choices include the Gaussian, exponential, and multi-quadrics functions. In this work, the Gaussian correlation function is adopted because it is infinitely differentiable and can provide highly smooth interpolation, which is well suited for continuously varying and smooth physical fields such as temperature. The Gaussian correlation function is defined aswhere is the distance between nodes and , and the shape parameter is determined using the leave-one-out cross-validation (LOOCV) method [35, 36].

The partial derivatives of the shape functions can be expressed aswhere denotes .

3.2 Moving mesh partial differential equation

The core idea of the MMPDE method is to improve numerical accuracy by controlling the tensor metric , while keeping the total number of degrees of freedom unchanged. In particular, is constructed such that an error estimator is redistributed as uniformly as possible over the computational domain, thereby achieving an equidistribution of numerical errors. Let and denote the physical domain and a simplicial mesh on it, respectively. For an element , denotes its physical coordinate. Similarly, let and denote the computational domain and a simplicial mesh on it, respectively. For an element , denotes its computational coordinate. These two simplicial meshes share the same number of nodes and the same number of elements [37]. For each element, there exists an invertible affine mapping , whose Jacobian matrix is denoted by [38]. Let the edge matrices of and be and , respectively. In this work, the metric tensor is defined aswhere is the Hessian matrix of the numerical solution , and is the identity matrix. This construction is motivated by the minimization of the -norm of the linear interpolation error, and it can effectively capture second-order variations of the solution. Since Equation 14 is computationally expensive in practical applications, an element-wise averaged metric defined in Equation 15, is often used:The MMPDE method interprets a general nonuniform mesh as a uniform mesh under the metric . Such meshes satisfy the equidistribution condition [39] and the alignment condition [40] in The normalization terms used in the equidistribution and alignment conditions are defined in Equation 18.where and denote the areas of the elements and , respectively. Equation 16 is the equidistribution condition, which controls the element sizes, whereas Equation 17 is the alignment condition, which controls the element shapes and orientations. To generate an optimal adaptive mesh that satisfies these two conditions as closely as possible, the MMPDE method minimizes a mesh functional via a variational formulation. The mesh functional is typically written aswhich can be regarded as a Riemann-sum approximation of the corresponding continuous functional. When and , the functional admits a minimizer [41–43]. However, the system resulting from minimizing Equation 19 is highly nonlinear and is difficult to solve directly. Therefore, the mesh evolution is commonly driven by the MMPDEwhere is a parameter that controls the responsiveness of the mesh to changes in the metric , and ensures that Equation 20 is invariant under transformations induced by . The derivative can be obtained using the differential calculus of scalar matrix functions. Equation 20 can be rewritten in the vertex-wise form given by Equation 21.where is the set of elements associated with the vertex , is the velocity of the vertex in element , and denotes the local index of the vertex in element .

The local nodal moving velocity can be defined asThe edge matrices are given in Equation 23, and the derivatives required for the local velocity calculation are expressed in Equations 24, 25.

The expressions of and areand the nodal velocity computed from Equation 22, the nodal positions at time can be updated accordingly, which effectively drives nodes toward regions where the numerical solution exhibits strong variations.

3.3 Local weak-form discretization for heat transfer problems

The problem domain and its boundary are discretized by field nodes. For an arbitrary node , , the local weak form of Equation 1 is constructed over its local integration subdomain . In the present study, , where , with denoting the nodal spacing and being a scaling parameter, typically taken in the range . The local weak form can be written aswhere denotes the weight function. Applying the divergence theorem to Equation 26 to reduce the order of derivatives yields

The Heaviside piecewise function is adopted as the weight function. With further simplification the Equation 27, one obtainswhere , is the internal boundary of the integration subdomain, is the part of the local boundary intersecting the Neumann boundary, and is the part intersecting the essential boundary. Using the Heaviside function eliminates the need to traverse the entire computational domain, thereby significantly reducing the computational cost and simplifying the numerical implementation [44]. For the temperature-field problem of mass concrete with embedded cooling pipes, the temperature approximation depends on time and the spatial coordinate , and can be expressed asandwhere is the shape-function matrix, while and are the column vectors of nodal temperatures and their time derivatives, respectively. Substituting Equations 29, 30 into Equation 28 yields the assembled system equation in Equation 31, with the component matrices and vector defined in Equation 32.whereConsequently, a linear system of ordinary differential algebraic equations with time as the independent variable is obtained. In this work, the time domain is discretized using the backward difference scheme.

3.4 Adaptive refinement of MM–MK algorithm

When the MM–MK algorithm is applied to transient problems, the adaptive procedure is carried out at selected physical time steps to track the evolution of the temperature field. Within each physical time step, a staggered coupling strategy is employed: the MMPDE is solved to update the node distribution, the temperature solution is transferred to the updated nodes, and the temperature field is then recomputed using the MK meshless method. In the present work, node adaptation is performed once per physical time step. For steady-state problems, the same procedure is embedded in an outer adaptive iteration framework until the relative changes in both the numerical solution and the node positions satisfy a prescribed convergence criterion. Within each adaptive iteration, the MMPDE is solved using a pseudo-time stepping procedure for mesh relaxation. The overall procedure is summarized in Algorithm 1.

4 Numerical examples

We employ the relative error norm and the root-mean-square error (RMS) as error indicators, as defined in Equations 33, 34:where and denote the exact and numerical values at the -th node, respectively; and is the total number of sampling nodes.

4.1 Example 1: convection–diffusion problem

We first consider a two-dimensional steady convection–diffusion problem. The accuracy of uniform and adaptive node distributions is compared to assess the effectiveness of the proposed MM–MK method. The computational domain is , which is discretized by 441 nodes .

The governing equation of the steady convection–diffusion problem iswhere coefficients in Equation 35 are given by , where is a prescribed constant diffusion coefficient. The boundary condition is given in Equation 36.with being the entire boundary of the domain.

The analytical solution is given in Equation 37.

To validate the effectiveness of the proposed algorithm, Table 1 reports the RN of the MM–MK method under uniform and adaptive node distributions. For a more comprehensive comparison, Table 1 also includes the results reported by Kaennakham [45] using other algorithms for the same problem. It should be noted that when , the problem becomes convection-dominated. In this regime, the pronounced convection gradients significantly increase the difficulty of the numerical solution.

As can be seen from Table 1, when , all methods listed in the table, regardless of whether uniform or adaptive nodes are employed, produce numerical solutions that agree well with the analytical solution, and the error remains below . When the case turns convection-dominated , the computational difficulty increases gradually as decreases. Nevertheless, the MM–MK method remains more accurate than ACNAA. In particular, for with uniform nodes, the accuracy of MM–MK is approximately 300 times that of the No-ACNAA method. Under adaptive nodes, the accuracy of MM–MK is about 8 times that of ACNAA. However, for the extremely convection-dominated case , the adaptive MM-MK result still yields RN = 1.05, indicating that although the method remains substantially more accurate than the comparison methods, the current nodal resolution is still insufficient for quantitatively accurate prediction in this regime. Moreover, by comparing the MK-based results under uniform and adaptive node distributions, the accuracy of MM–MK is approximately 37 times that of the fixed-node MK scheme. The results indicate that, compared with the global RBF collocation approach, the MK method based on a local weak form is more amenable to high-resolution capturing. Moreover, the adaptive node redistribution driven by the MMPDE method can further enhance the accuracy of the MK scheme in strong-gradient physical fields without increasing the degrees of freedom.

Figure 2 further illustrates the improvement achieved by adaptive node relocation over the fixed-node scheme at , where the exact solution is shown in Figure 2A. As can be clearly observed from the adaptive node distribution in Figure 2B, nodes migrate toward regions with large convection gradients, which is consistent with the underlying assumption of the adaptive strategy. Compared with the fixed-node numerical solution in Figure 2C, the adaptive solution in Figure 2D provides a noticeably better approximation. Together with Table 1 and Figure 2, these results demonstrate that the MM–MK algorithm can effectively enhance the quality of the numerical solution.

FIGURE 2

In addition, to examine the convergence behavior of the proposed method, Table 2 reports the RN obtained with different numbers of nodes. As the number of nodes increases from to , the relative error decreases monotonically from to , indicating the consistency of the numerical scheme. This trend also suggests that the method enters an asymptotic convergence regime once the discretization becomes sufficiently refined.

4.2 Example 2: one pipe with exact solution

To assess the numerical accuracy of the MM–MK algorithm for heat-conduction problems, a single-pipe cooling case is considered. The benchmark problem in Equation 38 describes the temperature field around a cooling pipe. In this scenario, the heat exchange between the pipe wall and the surrounding medium is intense, producing a steep temperature gradient near the wall. Consequently, the temperature varies sharply in space, which makes this example well suited for evaluating the capability of the proposed method in capturing strong-gradient regions. The computational domain is a square , where denotes time, is the pipe-center coordinate, the pipe radius is , and in Equation 38 is the exact solution. In this example, the initial temperature in the computational domain is set to 50, and the boundary conditions are obtained by directly differentiating Equation 38.

To determine reasonable choices of the total number of nodes and the time step for the MM–MK algorithm in pipe-cooling simulations, Table 3 reports the RMS error at for different node numbers and time steps . An analysis of the results in Table 3 shows that, when and , the error reaches , which meets the desired accuracy for engineering computations. Further reducing the time step to yields no substantial improvement, indicating that is sufficient to eliminate the dominant influence of temporal discretization error. With , the RMS error decreases monotonically as increases; however, the marginal improvement becomes much smaller when increases from 696 to 852. Balancing computational cost and accuracy, and are adopted in the following computations.

To highlight the advantage of adaptive node redistribution in MM–MK, Figure 3 compares the RMS errors at obtained by uniform and adaptive node distributions for different total numbers of nodes . It is evident that, under the same degrees of freedom, the adaptive strategy consistently outperforms the uniform distribution. When the uniform node number becomes sufficiently large, further refinement provides little improvement; by contrast, the adaptive strategy continues to reduce the error, indicating that newly added degrees of freedom are preferentially allocated to the error-dominant regions, thereby improving the accuracy under the same number of degrees of freedom. This result demonstrates the necessity and advantage of adaptive node relocation for strong-gradient heat-conduction problems.

FIGURE 3

Figure 3 indicates that MM–MK achieves good overall accuracy for heat-conduction simulations. However, for engineering pipe-cooling applications, more attention is typically paid to the temperature distribution in the strong-gradient region around the pipe. Figure 4 presents the temperature fields of the numerical and analytical solutions at , together with the distribution of the error RN. As shown in Figures 4A,B, the numerical and analytical temperature fields exhibit highly consistent patterns, and the isotherms nearly coincide, confirming the high approximation accuracy over the computational domain. Moreover, Figure 4C shows that the error is mainly concentrated in the strong-gradient region near the cooling pipe. The maximum relative error at the pipe is , while the error rapidly decays to a negligible level away from the pipe. This behavior is consistent with the physical characteristics of heat conduction and further verifies the accuracy of MM–MK in resolving strong-gradient regions.

FIGURE 4

The above results validate the accuracy of the proposed algorithm under a symmetric pipe configuration. To further evaluate the applicability and robustness of MM–MK under asymmetric strong-gradient conditions, Figure 5 compares the node distribution, numerical solution, and analytical solution at for an asymmetric pipe layout. Although the asymmetric arrangement leads to an evidently off-center temperature field, the numerical and analytical solutions in Figures 5C,D remain in close agreement in terms of the overall pattern. The node relocation in Figures 5A,B shows that nodes cluster toward the pipe region, indicating strong temperature-gradient variations there, which is consistent with the adaptive strategy. Furthermore, the RN error distribution in Figure 6 suggests that the error remains localized in the pipe neighborhood and decays rapidly away from the pipe without pronounced numerical oscillations or error diffusion. These observations demonstrate that MM–MK remains effective for asymmetric pipe-cooling cases and retains good robustness. Overall, the results confirm the stability and reliability of MM–MK in simulating strong-gradient heat-conduction problems.

FIGURE 5

FIGURE 6

4.3 Example 3: multiple pipes with pipe cooling

To examine the stability and accuracy of the proposed algorithm in cases where no analytical solution is available, the third example considers the temperature field in a mass-concrete block containing multiple cooling pipes. The computational domain is The Robin boundary condition imposed on is given in Equation 39.where is the convective heat-transfer coefficient between concrete and air, is the thermal conductivity of concrete, is the ambient air temperature, and is time. The pipe radius is , the initial temperature of concrete is , and the adiabatic temperature rise of concrete is given by . In this example, the top surface is exposed to air, the and boundaries are maintained at , and the is insulated. The pipe wall is in contact with cooling water, the pipe-wall temperature is prescribed by Equation 40.where the initial water temperature is .

For concrete components with embedded cooling pipes, accurately resolving the temperature field in the vicinity of the pipe is of particular importance. To quantify the numerical accuracy, the RMS error is evaluated using all nodes located within of the pipes at . A refined finite-element solution is adopted as the reference solution, where the mesh is locally refined around the pipes; the finite-element model has 2272 degrees of freedom and uses a time step of . As shown in Figure 7, for a fixed number of pipes, the RMS error decreases as the number of discretization nodes increases, indicating that the accuracy of MM-MK improves monotonically with increasing nodal resolution. Moreover, when , the RMS error for different pipe configurations remains below , demonstrating good adaptability of the proposed method with respect to varying pipe layouts.

FIGURE 7

To further examine the numerical stability of the proposed MM-MK method, additional tests were performed using larger time step sizes beyond those reported in Table 3. The numerical errors obtained with different time step sizes are presented in Table 4. As shown in Table 4, the MM-MK method consistently produces smaller numerical errors than the fixed-node MK method for all tested time step sizes. In particular, when larger time steps are used, the error of the fixed-node MK method increases more significantly, while the MM-MK method maintains relatively smaller errors. These results indicate that the MM-MK method is less sensitive to the time step size and exhibits improved numerical stability.

Furthermore, the temporal behavior of the error is illustrated in Figure 8, which presents the RN evolution over time for the four-pipe configuration. The error remains consistently low and exhibits no oscillatory behavior, indicating favorable temporal stability. In addition, Figure 9 shows the time histories of the maximum concrete temperature for the case without cooling pipes and for cases with different numbers of pipes. The results indicate that pipe cooling can significantly reduce the maximum temperature of the concrete. As the number of pipes increases, the peak temperature decreases progressively; in particular, using four pipes yields a noticeably lower than using two pipes. This indicates that increasing the number of cooling pipes can effectively improve the cooling performance and suppress the peak temperature rise in mass concrete. However, a larger number of pipes may also intensify local temperature gradients, which could increase the risk of thermal cracking. Therefore, in practical cooling design, the number of cooling pipes should be selected by balancing peak temperature reduction and temperature gradient control according to the specific thermal control requirements of the structure.

FIGURE 8

FIGURE 9

5 Conclusion

This study is primarily motivated by the strong temperature gradients in the vicinity of cooling pipes, for which conventional mesh-based methods may suffer from excessive local refinement, deteriorated numerical stability, and an overly large computational scale. To address these issues, an adaptive node strategy based on coupling the MK meshfree temperature solver with the MMPDE method is developed. The proposed coupling enables stable and accurate simulations of the strong-gradient temperature field around cooling pipes while keeping the total number of degrees of freedom unchanged.

A hybrid meshfree framework, referred to as MM–MK, is proposed. In this framework, the MMPDE method is treated as an independent node-generation module that performs adaptive node relocation in each time iteration based on a tensor metric. On the resulting node set, the heat-transfer problem of pipe cooling is solved using the MK interpolation that satisfies the Kronecker-delta property. In this manner, the proposed method alleviates the loss of accuracy and the poor stability induced by local over-refinement in strong-gradient regions.

In the first benchmark example, a convection-dominated strong-gradient problem is employed to verify the computational accuracy and stability of the proposed MM–MK method. In the second example, both symmetric and asymmetric pipe layouts are examined, and no pronounced numerical oscillations or error diffusion are observed, demonstrating the stability and reliability of the method for strong-gradient heat-conduction problems associated with pipe cooling. For multi-pipe configurations, the error near the pipes decreases continuously as the number of nodes increases.

Once the nodal resolution reaches a moderate level, the method maintains good accuracy in most tested cases and shows robust error reduction trends across different pipe-number configurations. Moreover, the error remains at a low level throughout the time evolution without oscillatory behavior, indicating good adaptability to various pipe layouts as well as favorable temporal stability.

Overall, the proposed MM–MK meshfree approach provides a numerical strategy that balances stability and accuracy for simulating the temperature field of mass concrete with embedded cooling pipes, and it offers an effective reference for cooling-system design and thermal crack-control analysis.

References

• 1.

  NguyenDHNguyenVTLuraPDaoVT. Temperature-stress testing Machine–a state-of-the-art design and its unique applications in concrete research. Cement and Concrete Composites (2019) 102:28–38. 10.1016/j.cemconcomp.2019.04.019

  • CrossRef
  • Google Scholar

• 2.

  LiLDaoVLuraP. Autogenous deformation and coefficient of thermal expansion of early-age concrete: initial outcomes of a study using a newly-developed temperature stress testing machine. Cement and Concrete Composites (2021) 119:103997. 10.1016/j.cemconcomp.2021.103997

  • CrossRef
  • Google Scholar

• 3.

  ZhuHHuYLiQMaR. Restrained cracking failure behavior of concrete due to temperature and shrinkage. Construction Building Mater (2020) 244:118318. 10.1016/j.conbuildmat.2020.118318

  • CrossRef
  • Google Scholar

• 4.

  ChenGDingXCaiMQiaoW. Analytical solution for temperature field of nonmetal cooling pipe embedded in mass concrete. Appl Therm Eng (2019) 158:113774. 10.1016/j.applthermaleng.2019.113774

  • CrossRef
  • Google Scholar

• 5.

  LiQChenGZhuF. Simulation of thermal field in mass concrete with cooling pipes based on the isogeometric analysis method. Front Phys (2024) 12:1338718. 10.3389/fphy.2024.1338718

  • CrossRef
  • Google Scholar

• 6.

  HongYLinJChengAH-DWangYChenW. 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

  • CrossRef
  • Google Scholar

• 7.

  ZhuB. The equivalent heat conduction equation of pipe cooling in mass concrete considering influence of external temperature. J Hydraul Eng (2003) 34 (3):49–54.

  • Google Scholar

• 8.

  TasriASusilawatiA. Effect of material of post-cooling pipes on temperature and thermal stress in mass concrete. Structures (Elsevier) (2019) 20:204–12. 10.1016/j.istruc.2019.03.015

  • CrossRef
  • Google Scholar

• 9.

  HongYLinJVafaiK. Thermal effect and optimal design of cooling pipes on mass concrete with constant quantity of water flow. Numer Heat Transfer, A: Appl (2020) 78:619–35. 10.1080/10407782.2020.1805222

  • CrossRef
  • Google Scholar

• 10.

  Van TranMChauVNNguyenPH. Prediction and control of temperature rise of massive reinforced concrete transfer slab with embedded cooling pipe. Case Stud Construction Mater (2023) 18:e01817. 10.1016/j.cscm.2022.e01817

  • CrossRef
  • Google Scholar

• 11.

  YanFLuoZGengYLiX. Study on air-pipe cooling effect and cooling strategy of mass reinforced concrete wall. Scientific Rep (2025) 15:10988. 10.1038/s41598-025-87033-4

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 12.

  ChengJLiTLiuXZhaoL. 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:487–508. 10.1002/nag.2409

  • CrossRef
  • Google Scholar

• 13.

  LiuX-h.DuanYZhouWChangX. Modeling the piped water cooling of a concrete dam using the heat-fluid coupling method. J Eng Mech (2013) 139:1278–89. 10.1061/(asce)em.1943-7889.0000532

  • CrossRef
  • Google Scholar

• 14.

  GuoWWangFWuYHanFYuCLiuJet alA non-equilibrium micromechanics-based thermo-hydro-mechanical model for freezing/thawing in saturated cementitious materials: from elasticity to elastic-plasticity. Cement Concrete Res (2023) 173:107267. 10.1016/j.cemconres.2023.107267

  • CrossRef
  • Google Scholar

• 15.

  GuoWWangFWuYJiangJXuW. New insights into freezing behavior of saturated and air-entrained porous media via a micromechanics-based thermo-hydro-mechanical model. Water Resour Res (2023) 59:e2022WR034211. 10.1029/2022wr034211

  • CrossRef
  • Google Scholar

• 16.

  GuoWWuYJiaMMaZJiangJXuW. Micromechanics-based thermo-hydro-mechanical model for air-entrained porous materials subjected to freezing-thawing cycles. J Rock Mech Geotechnical Eng (2025) 17 (7):4413–4428. 10.1016/j.jrmge.2024.11.055

  • CrossRef
  • Google Scholar

• 17.

  StrouboulisTZhangLBabuškaI. Assessment of the cost and accuracy of the generalized fem. Int Journal Numerical Methods Engineering (2007) 69:250–83. 10.1002/nme.1750

  • CrossRef
  • Google Scholar

• 18.

  YuY-YLinYJiZ-S. New method for ship finite element method preprocessing based on a 3d parametric technique. J Marine Science Technology (2009) 14:398–407. 10.1007/s00773-009-0058-1

  • CrossRef
  • Google Scholar

• 19.

  WangYYuYLinY. Isogeometric analysis with embedded stiffened shells for the hull structural mechanical analysis. J Mar Sci Technology (2022) 27:786–805. 10.1007/s00773-021-00868-0

  • CrossRef
  • Google Scholar

• 20.

  NguyenVPRabczukTBordasSDuflotM. Meshless methods: a review and computer implementation aspects. Mathematics Computers Simulation (2008) 79:763–813. 10.1016/j.matcom.2008.01.003

  • CrossRef
  • Google Scholar

• 21.

  PatelVGRachchhNV. Meshless method–review on recent developments. Mater Today: Proceedings (2020) 26:1598–603. 10.1016/j.matpr.2020.02.328

  • CrossRef
  • Google Scholar

• 22.

  DaxiniSDPrajapatiJM. Parametric shape optimization techniques based on meshless methods: a review. Struct Multidisciplinary Optimization (2017) 56:1197–214. 10.1007/s00158-017-1702-8

  • CrossRef
  • Google Scholar

• 23.

  ShojaeiABoroomandBMossaibyF. A simple meshless method for challenging engineering problems. Eng Computations (2015) 32:1567–600. 10.1108/ec-06-2014-0131

  • CrossRef
  • Google Scholar

• 24.

  VakilhaMSaghatchiRAlexiadisAYildizMShadlooM. A fully explicit incompressible smoothed particle hydrodynamics approach for modeling transient heat transfer and thermo-capillary flows. Comput and Fluids (2024) 269:106112. 10.1016/j.compfluid.2023.106112

  • CrossRef
  • Google Scholar

• 25.

  YuZHeYTanFTongD. Three-dimensional numerical manifold method for thermal stress calculation of mass concrete. Construction Building Mater (2025) 491:142769. 10.1016/j.conbuildmat.2025.142769

  • CrossRef
  • Google Scholar

• 26.

  HouDWangLYanJLiewKM. Vibration analysis of a strain gradient plate model via a mesh-free moving kriging interpolation method. Eng Anal Boundary Elem (2022) 135:156–66. 10.1016/j.enganabound.2021.11.014

  • CrossRef
  • Google Scholar

• 27.

  WangFLinGZhouY-h.ChenD-h.Element-free galerkin scaled boundary method based on moving kriging interpolation for steady heat conduction analysis. Eng Anal Boundary Elem (2019) 106:440–51. 10.1016/j.enganabound.2019.05.027

  • CrossRef
  • Google Scholar

• 28.

  ChenLLiewKM. A local petrov-galerkin approach with moving kriging interpolation for solving transient heat conduction problems. Comput Mech (2011) 47:455–67. 10.1007/s00466-010-0553-6

  • CrossRef
  • Google Scholar

• 29.

  WuWLiuY-LZhangA-MLiuM. An h-adaptive local discontinuous galerkin method for second order wave equation: applications for the underwater explosion shock hydrodynamics. Ocean Eng (2022) 264:112526. 10.1016/j.oceaneng.2022.112526

  • CrossRef
  • Google Scholar

• 30.

  MoreiraCACaicedoMACerveraMChiumentiMBaigesJ. A multi-criteria h-adaptive finite-element framework for industrial part-scale thermal analysis in additive manufacturing processes. Eng Comput (2022) 38:4791–813. 10.1007/s00366-022-01655-0

  • CrossRef
  • Google Scholar

• 31.

  HeXHuangJYangDZhouYHuangX. An efficient dynamic p-adaptive discontinuous galerkin method for seismic wave simulations. IEEE Trans Geosci Remote Sensing (2025) 63. 10.1109/TGRS.2025.3600174

  • CrossRef
  • Google Scholar

• 32.

  OmellaÁJPardoD. r-adaptive deep learning method for solving partial differential equations. Comput and Mathematics Appl (2024) 153:33–42. 10.1016/j.camwa.2023.11.005

  • CrossRef
  • Google Scholar

• 33.

  YangYYangQDengYHeQ. Moving sampling physics-informed neural networks induced by moving mesh pde. Neural Networks (2024) 180:106706. 10.1016/j.neunet.2024.106706

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 34.

  SubichCJ. A robust moving mesh method for spectral collocation solutions of time-dependent partial differential equations. J Comput Phys (2015) 294:297–311. 10.1016/j.jcp.2015.04.003

  • CrossRef
  • Google Scholar

• 35.

  RippaS. An algorithm for selecting a good value for the parameter c in radial basis function interpolation. Adv Comput Mathematics (1999) 11:193–210. 10.1023/a:1018975909870

  • CrossRef
  • Google Scholar

• 36.

  CavorettoRDe RossiA. An adaptive loocv-based refinement scheme for rbf collocation methods over irregular domains. Appl Mathematics Lett (2020) 103:106178. 10.1016/j.aml.2019.106178

  • CrossRef
  • Google Scholar

• 37.

  LopezMShontzSMHuangW. A parallel variational mesh quality improvement method for tetrahedral meshes based on the mmpde method. Computer-Aided Des (2022) 148:103242. 10.1016/j.cad.2022.103242

  • CrossRef
  • Google Scholar

• 38.

  TannahillCWanJ. Mm-admm: implicit integration of mmpdes in parallel. Comput and Mathematics Appl (2023) 141:67–79. 10.1016/j.camwa.2023.03.019

  • CrossRef
  • Google Scholar

• 39.

  HuangW. Mathematical principles of anisotropic mesh adaptation. Commun Comput Phys (2006) 1:276–310. 10.4208/cicp.2006.v1.p276

  • CrossRef
  • Google Scholar

• 40.

  ZhangFHuangWLiXZhangS. Moving mesh finite element simulation for phase-field modeling of brittle fracture and convergence of newton’s iteration. J Comput Phys (2018) 356:127–49. 10.1016/j.jcp.2017.11.033

  • CrossRef
  • Google Scholar

• 41.

  HuangWKamenskibLSibH. Mesh smoothing: an mmpde moving mesh approach. In: Proceedings of the 24th international meshing roundtable (2015). p. 1–5.

  • Google Scholar

• 42.

  HuangWKamenskiL. A geometric discretization and a simple implementation for variational mesh generation and adaptation. J Comput Phys (2015) 301:322–37. 10.1016/j.jcp.2015.08.032

  • CrossRef
  • Google Scholar

• 43.

  DassiFKamenskiLFarrellPSiH. Tetrahedral mesh improvement using moving mesh smoothing, lazy searching flips, and rbf surface reconstruction. Computer-Aided Des (2018) 103:2–13. 10.1016/j.cad.2017.11.010

  • CrossRef
  • Google Scholar

• 44.

  AtluriSNShenS. The basis of meshless domain discretization: the meshless local petrov–galerkin (mlpg) method. Adv Comput Mathematics (2005) 23:73–93. 10.1007/s10444-004-1813-9

  • CrossRef
  • Google Scholar

• 45.

  KaennakhamSChuathongN. An automatic node-adaptive scheme applied with a rbf-collocation meshless method. Appl Mathematics Comput (2019) 348:102–25. 10.1016/j.amc.2018.11.066

  • CrossRef
  • Google Scholar

• 46.

  DongXGuoWDengJGePLiuZ. A dual-scale damage model for freezing–thaw degradation in cement-based materials under thermo-hydro-mechanical coupling. J Build Eng (2026) 120:115519. 10.1016/j.jobe.2026.115519

  • CrossRef
  • Google Scholar