Simulation of chemo-thermo-mechanical problems in cement-based materials with Peridynamics

A chemo-thermo-mechanical problem is solved using a peridynamic approach to investigate crack propagation in non-reinforced concrete at early-age. In the present study, the temperature evolution and the variation of the hydration degree in conjunction with the mechanical behaviour of cement-based materials are examined. Firstly, a new peridynamic model is introduced to solve fully coupled chemo-thermal problems by satisfying thermal equilibrium condition and hydration law simultaneously and then the effects of the chemo-thermal analysis are imposed in the mechanical framework to investigate all the interactions. The proposed approach is used to solve 2D chemo-thermo-elastic problems and then it is applied to investigate the fracture of concrete structures. Additionally, we examine the accuracy of the method by comparing the crack paths, temperature and hydration degree with those achieved by applying other numerical methods and the experimental data available in the literature. A good agreement is obtained between all sets of results.


List of symbols A
Cross-sectional area A T Ratio of the maximum value of the heat production rate to the latent hydration heat a Hydration constant determined by interpolating experimental data b Hydration constant determined by interpolating experimental data b Body Cement hydration and shrinkage is a phenomenon characterized by inhomogeneous stresses occurring in plain and non-reinforced concrete. Concrete is widely used in infrastructures such as airfields, slabs and dams [1]. A significant disadvantage of concrete is its sensitivity with respect to the early-age cracking in aforementioned structures. Hence, a basic understanding of material performance due to early-age cracking is required in the design of concrete structures [1,2]. Due to the fast change of mechanical properties such as strength of the material [3,4], prediction of the behaviour of concrete is still not a trivial task. In fact, during cement hydration process, the mechanical properties of the material strongly tie in with the physical and chemical processes. In this multi-physics problem all the thermal, mechanical and chemical effects should be considered [2]. Therefore, it is still a challenging issue to precisely predict the early-age behaviour of the cement based material, especially when discontinuities such as cracks occur. Numerous numerical and experimental studies have been so far performed to analyse the mechanical behaviour of early-age non-reinforced concrete. Among the experimental studies, Kolver et al. [5] and Qstergaard et al. [6] have investigated the tensile creep of concrete at early-age. In [7], a restrained shrinkage test has been developed to characterize the early-age mechanical behaviour. Moreover, the type of cement and effects of curing temperature on the shrinkage properties have been investigated in [8]. The emergence of cracks due to autogenous shrinkage in high performance and slag cement concrete have been investigated in [9,10], respectively. Later, early age hygro-thermo-chemo-mechanical reactions have been studied in [11][12][13][14][15][16][17][18][19].
For early age behaviour of concrete the numerical studies by [20][21][22][23] should be mentioned here. Moreover, in [24][25][26], cracking risk analysis is performed by using smeared-crack and crack band models, respectively. Other numerical approaches such as microplane approach in [27] or lattice approach in [28] are used to investigate the early-age concrete behaviour. Soon afterwards, Neguyen et al. [1,29] have proposed a phase field model to investigate complex fracture induced by early-age shrinkage and hydration heat in cement based materials. Other numerical methods in this field can be traced in [30][31][32]. However, due to the complexity of the hydration process, available computational approaches cannot precisely describe how complex crack patterns nucleate and evolve in concrete structures. The majority of the proposed computational techniques are yet able to merely predict the crack propagation induced by cement hydration qualitatively [11,12] and their application to 3D problems often turns out to be very complicated. Another drawback of these techniques in the literature is to use ad hoc modifications and simple assumptions. In some of these techniques, a large number of tests to characterize the model inputs are required [12,33]. Moreover, gradients/divergence terms in the governing equation could be a source of inconsistency with the onset of material discontinuity [34].
The main purpose of the present study is the solution of chemo-thermo-mechanical problems in cement based materials, including crack initiation and propagation, using a peridynamic based approach. In fact, the peridynamic approach offers some advantages [67,68], in terms of accuracy and efficiency, which can be highly important for the coupled chemothermo-mechanical problems conducted in the present study. The proposed model includes, the effects of hydration reactions on the mechanical properties of material. The model exhibits highly accurate results, with respect to experiments in the literature, not only in terms of displacement but also in terms of crack path, temperature and hydration degree. One of the major advantages of the proposed approach is that there is no need to re-mesh the model when there are discontinuities such as cracks. These features culminate in a robust model that performs properly even in situations where the nucleation and propagation of cracks with complex patterns occurs.
The paper is structured as follows. Section 2 provides the fundamentals of all the interactions taken into account during the hydration process. A peridynamic based mathematical model is developed in section 3, while in sect. 4, we explain the step by step of the procedure of the algorithm. Model validation by using experimental data is shown in Sect. 5. In Sect. 6, the accuracy of the method is investigated by several examples. Section 7 concludes the paper.

Fundamentals
Let us consider a body which in an initial configuration X occupies a 2D space (see Fig. 1). By applying chemo-thermo-mechanical conditions on X or on its boundary, oX, four state variables including displacement, uðx; tÞ, temperature Tðx; tÞ, hydration degree aðx; tÞ and the damage level /ðx; tÞ vary over the body. By calculating these four variables, one can determine the state of the whole system. The various fields are all interacting in this multi-physics model (see Fig. 2). In Fig. 2, each arrow represents the interaction direction e.g. chemical and thermal analysis influence the mechanical analysis because variations in temperature and hydration affect the displacements and therefore may generate cracks. On the other hand crack propagation determines a variation of thermal conductivity and of hydration which is the only influence mechanics has on the thermal and hydration fields. This implies that, fracture analysis is the direct result of the mechanical analysis and we use the fracture effect by using a history-dependent function which is discussed later in the mechanical kernel function.
In order to analyse this multi-physics model through an appropriate scheme, we decouple it into two main problems as follows: 1. Fully coupled chemo-thermal problem 2. Chemo-thermo-mechanical coupled problem In the next section, we explain how Peridynamics is used to solve the aforementioned problems.

Mathematical model
This section shows the equations adopted to study the multi-physics model using the peridynamic approach.

Fully coupled chemo-thermal problem
Based on Fourier's law, the thermodynamic equilibrium of thermal flux in a concrete specimen due to hydration reaction is governed by: in which, q and c v are the density of the material and the specific heat capacity. Q 1 stands for the potential heat of the hydration reaction and a is the hydration degree which describes the level of the hydration reaction between cement and water. a takes a value between 0 and 1, 0 6 a 6 1. a ¼ 1 means full hydration between cement and water takes place. Furthermore, T and k represent the temperature and the thermal conductivity respectively. Considering a constant value for thermal conductivity, Eq. (1) takes the following form: Bobaru and Duangpanya [44] firstly proposed the application of the Bond Based Peridynamics (BB-PD) to heat conduction problem. The heat flux is exchanged between a material point x (the source node) and the family points located within its integration domain; i.e, x 0 2 HðxÞ as shown in Fig. 1. Hence, the transient form of the BB-PD thermal diffusion equation is expressed by: in which h s ðx; tÞ is the heat source per unit volume and the kernel function, f h , is given by Chen and Bobaru [69]: where In BB-PD, two material points have the relative position n ¼ x À x 0 known as bond. We assume that the heat flux depends on the temperature difference between the points and it occurs over a bond connecting two material points. In Eq. (4), c TH stands for the thermal micro-conductivity in peridynamics correlated with the standard conductivity. This correlation can be obtained by equating the classical thermal potential to the peridynamic thermal potential at a material point with a specific horizon size, d, (the details can be found in [45,46]). The so-called horizon specifies the size of the region where non-local interactions occur. In Eq. (4), n is an integer, normally selected to be 0, 1, or 2 [44][45][46]. Based on [69], different values of n in the representation of the kernel function have been selected in the past. For example, a constructive approach to derive the specific form of the kernel function leads to n ¼ 2 (see [44,45]). The value n ¼ 1 is used in [46], and this value has been used most often for the mechanical peridynamic formulation [70]. A comprehensive study on different values of n can be found in [69]. In this paper, n ¼ 1 is used as suggested in [46]. Furthermore, the thermal micro-modulus c TH can be represented in terms of the thermal conductivity k for 1D, 2D and 3D cases as: where t and A stand for thickness of the body (2D case) and the cross-sectional area allocated to each material point, respectively. In the coupled chemo-thermal problem, the heat source per unit volume h s ðx; tÞ is due to hydration reaction and Eq. (2) can be rewritten as: In Eq. (7), we assume that the evolution of heat release due to the hydration process is obtained by the Arrhenius law [23]. The heat released during hydration reaction is defined by: where A T represents the ratio of the maximum value of the heat production rate to the latent hydration heat for a normalized definition of the hydration function [1,29,32], the damage variable / affects the heat during hydration reaction as assumed in [29]. Moreover, E a is the activation energy characterizing the rate of heat generation and R ¼ 8:314 Â 10 À3 kJ K À1 mol À1 is the ideal gas constant. In Eq. (8), / stands for the damage level of point x which is described in the next section where the fracture analysis effects are introduced in the chemo-thermomechanical problem. In fact, / in Eq. (8) takes a value between 0 and 1. Moreover, f ðaÞ is the evolution function providing the normalized heat production rate in terms of the hydration degree. Generally, there are many ways to choose this specific function e.g. a piecewise linear or exponential approximation of the experimental data [1,20,29,33,71]. In the present study, the power form of the evolution function, f ðaÞ, for the normalized heat production rate is considered as presented in [1,29,32]: The constants a, b and c in Eq. (9) are determined by interpolating experimental data [29]. The relation between the aforementioned constants is expressed as [20,29,32,71]: We can discretize Eq. (7) in different ways [72]. In this contribution, the discretization is performed by a meshfree approach originally introduced in [36]; this simple discretization have been commonly applied in the literature (see Fig. 3). In 2D cases, the domain is discretized by a grid of points known as nodes. We assign a volume equal to Dx 2 t to each node, where t is the specimen thickness. In this way, each node is located in a center of a square cell as depicted in Fig. 3. Each node x i interacts with all the nodes within its integration domain as shown in Fig. 3. Hereinafter, x i and x j are called the source node and the family node, respectively. The horizon is considered as d ¼ m Dx, where m is the ratio between d and the grid size Dx. In the usual way, d, m and Dx are considered uniform in the whole domain.
Moreover, one can discretize the time into instants as t 1 ; t 2 ; . . .; t n .
Therefore, by choosing the one-point Gauss quadrature rule for the integration in space, one can express the discretized form of Eq. (7) as: where V j stands for the volume of node x j and N represents the total number of family nodes. Finally, bðx j À x i Þ is the volume correction factor which determines the portion of V j that falls within the interaction domain of source node x i . In this paper, bðx j À x i Þ is determined as suggested in [73].
To solve Eq. (11), different time integration schemes can be taken into account. In this paper, we employ a forward difference time marching; in this way, after obtaining the field variables (a n i and T n i ) at time instant t n , one is able to advance to the next time step via set of coupled equations as follows: The unknown values in Eq. (12), are T nþ1 i and a nþ1 i . Various techniques can be used to solve Eq. (12), in this study, the Newton Raphson (N-R) method is adopted. In Eq. (12), Dt CT is the incremental time step related to the chemo-thermal analysis, to obtain a stable solution we must choose it within a proper range. The upper bound of the range can be expressed by [36,46]:

Boundary conditions
Chemo-thermal boundary conditions of Dirichlettype, must be applied on corresponding boundary surfaces represented by oX D . For this type of conditions, a temperature field should be assigned to each node located within a distance Dx=2 close to oX D (Fig. 1); however, Neumann-type boundary conditions is imposed in the form of convection, heat flux and radiation. To impose heat flux boundary conditions, firstly we evaluate the rate of heat entering into the body from the boundary surface, then we must convert the rate of flowing heat _ Q to a volumetric valueQ (heat generation per unit volume) by [45]: in which, q is the heat flux, S is the area to which the heat flux is imposed and V f is the volume of the boundary region. By using Eq. (14), one can determine the heat flux portion of each node on the corresponding boundary. The convection boundary condition (Neumann), in which a heat transfer between the surface of a body and its surrounding occurs, is expressed by [45,46]: in which h is the convective heat transfer coefficient. Moreover, T s ðx; tÞ and T a ðx; tÞ are the body surface and air temperature, respectively. In this way, volumetric heat generation due to the convection boundary condition is defined by [45]: For the chemical part, a Neumann-type boundary condition is assigned to the hydration degree (a); furthermore, the initial value of the hydration degree needs to be prescribed on the whole domain.
In the next subsection, we illustrate the strategy to cope with the mechanical part of solution, using Peridynamics. The same discretization scheme in space is used in the mechanical solution. Therefore, we describe how to incorporate the chemo-thermal effects into the mechanical solution.

Chemo-thermo-mechanical coupled problem
In this subsection, the basic concepts in addition to the main ideas of the peridynamic model for fracture are presented. Let us consider again a body, which an initial configuration X, occupies a 2D space (see Fig. 4). Based on BB-PD, the equation of motion is given by: where f stands for the pairwise force function of each bond and similar to the chemo-thermal problem, two material points interact only within a finite distance, d. Moreover, q, u and b are the mass density, displacement field and body force density, respectively. In Eq. (17), the position of point x at time t is denoted by yðx; tÞ ¼ x þ uðx; tÞ. In Fig. 4, the relative displacement vector between two points is defined as g ¼ u 0 À u. By using a constitutive law, with the Prototype Microelastic Brittle (PMB) model [36], the pairwise force function is expressed by: is the thermal expansion coefficient and s denotes the relative elongation (stretch) of a bond, defined by [35]: In Eq. (18), the average temperature of a bond is defined by: where T 0 represents the reference temperature, T and T 0 indicate the temperature at material points x and x 0 respectively. Similarly, the average level of hydration degree is expressed as: where a 0 represents the reference hydration degree, a and a 0 indicate the hydration degree at material points x and x 0 respectively. The material constant j, in Eq. (18) denotes the evaluation of autogenous shrinkage when hydration degree is greater than the mechanical percolation threshold a au and the : h i þ represents the positive operator [29]. To take into account the shrinkage effect of hydration phenomenon, j Eq. (18) has been introduced. Moreover, lðn; tÞ is a history-dependent function which takes the value of 0 for broken bonds and 1 for active bonds. c ME ðaÞ stands for the micro-modulus defined in Eq. (22). xðnÞ determines the degree of interactions between points. Taking into account the elastic deformation energy and xðnÞ ¼ 1, one is able write the mechanical micro-modulus in terms of the horizon radius d and of the Young's modulus E and as: c ME ðaÞ ¼ 9EðaÞ p td 3 ; Plane stress 48EðaÞ 5p td 3 ; Plane strain 12EðaÞ Moreover, the critical stretch, s 0 ðaÞ, is defined in terms of the critical energy release rate of the material, G 0 ðaÞ [36,44]: Due to the processes that occur in the hardening cement paste, its mechanical properties are strongly affected by the hydration value a [1,29], since in the previous formula the mechanical parameters depend on the hydration value. The Young's modulus, EðaÞ and the critical energy release rate, GðaÞ 0 , are given by [1,29]: where E 1 and G 01 are the final Young's modulus and final (when the hydration is completed) energy release rate, respectively. a a E and a a G 0 are two parameters in the range of (0, 1) related with the hydration value of the material a. Based on [24,29,74], a a E , a a G 0 are defined by: where a E and a G 0 are two threshold hydration degree constants, which express the hydration degree when material starts having both the mechanical properties, EðaÞ and G 0 ðaÞ greater than zero. We break a bond and remove it from further calculations when: Moreover, the damage level of point x at time t is given by: where / represents the ratio of the number of broken bonds to the total number of bonds initially connected to point x and its value is between 0 and 1. The case / ¼ 0 means that there is no broken bonds between the point and all surrounding points within its horizon, while / ¼ 1 represents a complete damaged state. It is well known that PBM constitutive law may lead to overstimations of the peak-load [75][76][77][78] and to difficulties in the capability to describe the post-peak behaviour. A possible strategy for reducing this effect is to adopt bi or tri-linear constitutive laws or nonuniform material resistance as indicated in [79]. Similar to the chemo-thermal part, to discretize Eq. (17), one can adopt the one-point Gauss quadrature rule as: In the present contribution, a quasi-static solution is obtained using the dynamic relaxation method firstly introduced in [80]. As suggested in [81,82], the equation of motion of each node in Eq. (29) is assembled, into a global system of equations and therefore one can introduce an artificial damping to the system by: K € UðX; tÞ þ CK _ UðX; tÞ ¼ FðU; U 0 ; X; X 0 Þ ð30Þ in which C stands for the damping coefficient and K is the fictitious diagonal density matrix. Furthermore, displacement vector U and the positions vector, X in Eq. (30) are given by: Moreover, F stands for the summation of the external and internal forces, whose its i-th component can be calculated: where M is the total number of nodes, and the density matrix, K, is defined according to [81]. For the damping ratio C, values in the range of 10 6 À 10 7 kg/m 3 s are used, as suggested in [82]. The time integration is performed by adopting a central difference explicit integration approach. Having found the displacement and acceleration of each node x i at t n , ðu n i ; € u n i Þ, displacements and velocities at t nþ1 ¼ t n þ Dt ME can be calculated by: In this way, the solver is able to advance to the next time step by: in which Dt ME represents the constant mechanical time step. For the first time iteration, the velocity at t 1=2 is obtained by:  For the dynamic relaxation method, as recommended in [81] the time step is taken as Dt ME ¼ 1.
4 The step by step procedure of the method 1. Inputting the constant parameters.
2. Pre-processing: applying the chemical, thermal and mechanical initial conditions. 3. Chemo-thermal analysis for step ðn þ 1Þ and calculate T nþ1 and a nþ1 by Eq. (12  To validate the present numerical method, we simulate a concrete cubic sample which is taken from an experimental study [1,83]. The dimensions of the cubic sample is 150 Â 150 Â 150 mm. The concrete properties are taken from [1,83] and are given in Table 1. The testing sample is equipped with thermal sensors located at the top-center of the sample and the whole system is stored in a polystyrene box with a thickness of 50 mm. The chamber temperature is 20 C and its humidity is 50%. The initial hydration degree all over the slab is a 0 ¼ 0:05. Furthermore, the initial temperature of the concrete slab is T 0 ¼ 20 C and the final temperature is T 1 ¼ 30 C. The model is discretized uniformly by a grid spacing of Dx ¼ Dy ¼ Dz ¼ 0:075 m and a total number of 9261 nodes. The incremental time step is taken as Dt ¼ 1 s and the simulation is carried out up to 600 h. The time history of the Young's modulus is represented in Fig. 5. There is an acceptable agreement between the Young's modulus of the three experimental tests [1,83] and that of the peridynamic model. We also investigate the thermal behaviour of the system. The temperature history is shown in Fig. 6. It can be observed that the maximum temperature reaches up to 38 C at the concrete age of 18 h in the peridynamic model and it is in a remarkable agreement with the experimental measurements [1,83]. Furthermore, the accumulated heat generated due to the hydration process in the peridynamic model is compared to the experimental test in Fig. 7. A very good agreement is achieved.

Shrinkage phenomenon
In order to get a deep insight into the chemo-thermomechanical behavior, the shrinkage phenomenon of concrete is examined. Shrinkage is the change of volume of the concrete due to the change of the hydration degree and it is one of the significant aspects of the early-age concrete. The experimental set up (rectangular slabs) is described in [1,83]. The dimensions of the slabs are 1000 Â 100 mm and the related numerical model is solved in a plane strain condition. C20/25 concrete is used in the experimental tests and its properties are taken from [1,83] (see Table 1). The aim of this test is to measure the autogenous shrinkage strain. To this end, the moisture exchange with the ambient medium is prevented. The ambient temperature was measured during the test in [1,83] as depicted in Figure 8 and used into the numerical model to impose a convection boundary   13 Example II: the evolution of the a, d, g damage level /, b, e, h temperature T ( C) and c, f, i hydration degree a during the hydration process using the present approach condition on all the boundaries [1,83]. On the left boundary of the concrete slab, both x and y-displacements are fixed, while on all other boundaries both x and y-displacements are free. The right side of the slab is in touch with a displacement transducer which records all displacement within a range of 2 mm.The slab is discretized uniformly by a grid spacing of Dx ¼ Dy ¼ 0:005 m and a total number of 4221 nodes. The incremental time step is taken as Dt ¼ 10 s and the simulation is carried out up to 28 days. In Fig. 9, the development of shrinkage strain is presented. Using the peridynamic model, the shrinkage strain reaches À184 lm=m after 28 days which is in a very good agreement with the two experimental tests in [1,83]. Up to here, we presented simple models to investigate the variation of the elastic modules, the evolution of the temperature and the shrinkage strain in early-age concrete. Moreover, we compared them to Fig. 14 Example II: the evolution of the phase field d, temperature T ( C) and hydration degree a during the hydration process using the phase field method [29] the experimental data available in the literature. In the rest of this study, we investigate more complex problems by three other examples.
6 Numerical examples 6.1 Example I: chemo-thermo-mechanical analysis of a concrete slab In this example, the temperature variation of a concrete slab, due to the early-age behaviour of concrete, is investigated. Using this example, we aim at validating the fully chemo-thermal reactions via a peridynamic approach. This example has been investigated in [31,32] by using FEM. The authors used two different types of elements known as hybrid and conform, then they compared their results with the experimental data in [84]. Heat transfer occurs across the slab thickness only. This problem is solved in a plane strain condition. The slab with the dimensions of 1:0 Â 0:35 m is considered here (see Fig. 10.) Two thermal sensors, A and B, are placed at 0.05 m from the bottom and the top surfaces respectively. The same material parameters used in [29,31] are assumed to solve this example. These parameters are given in Table 2. The convective boundary conditions have coefficients h c ¼ 7:5 W/(m 2 K) (convection with air) and h c ¼ 4:5 W/(m 2 K) (convection with soil) at the bottom-end (P1-P2) and the top-end (P3-P4) respectively. The environment air temperature represented in Fig. 11 is changing during the simulation. The initial hydration degree all over the slab is a 0 ¼ 0:05. Furthermore, the initial temperature of the concrete slab and the surface of the supporting soil is T 0 ¼ 17:5 C and T 0 Soil ¼ 17 C respectively. The isolated boundary condition is applied to the right (P2-P3) and the left (P1-P4) sides of the slab (see Fig. 10). Furthermore, the x-displacements are fixed in the left and the right sides of the slab while y-displacements are free [29,31]. On the bottom-end (P1-P2), y-displacements are fixed, while the x-displacements are free [29,31]. On the top side of the slab (P3-P4), both x and y-displacements are free. The slab is discretized uniformly by a grid spacing of Dx ¼ Dy ¼ 0:05 m and a total number of 14,271 nodes. The incremental time step is taken as Dt ¼ 10 s. Two nodes x A ¼ ð0:5; 0:05Þ m and x B ¼ ð0:5; 0:3Þ m are chosen to investigate the time history of the temperature. The time history of the temperature for nodes x A and x B are shown in Fig. 11. The results for both of the nodes are in an acceptable agreement with the FEM models (conform and hybrid) in [31,32]. Moreover, the peridynamic temperature time history is in a good agreement with the measured data in [84]. In this example, we investigate the chemo-thermomechanical crack propagation in a homogeneous circular concrete specimen. A square hole is considered at the centre of this specimen as shown in Fig. 12. This example has been solved in [29] using the phase field method. The material parameters in [29] are considered in this example and are listed in Table 3. The radius of the circular specimen and the size of the square hole are taken as r ¼ 0:5 m and A convection boundary condition with the convection coefficient of h ¼ 8 W/(m 2 K) is applied over the system as shown in Fig. 12. The ambient temperature is assumed to be constant during the whole interaction process and is taken as T a ¼ 15 C. To impose the mechanical boundary conditions, the left and the right sides of the square hole are fixed along the x-direction while the y-displacements are free. Moreover, the top and the bottom edges of the square hole are fixed along the y-direction while the x-displacements are free.
The slab is discretized uniformly by a grid spacing of Dx ¼ Dy ¼ 0:005 m and a total number of 25156 nodes. The incremental time step is taken as Dt ¼ 10 s for a total of 72000 steps. The contour plots of the damage, temperature field and hydration degree for three different time steps are depicted in Figure 13. The comparison of Fig. 13, with the phase field results taken from [29] in Fig. 14, shows that there is an acceptable agreement between the peridynamic and phase field approaches. As shown in Figure 13g, the damage evolution in the peridynamic model is absolutely symmetric which shows there is significantly low round-off error in our simulations. The comparison of Figs. 13 and 14 is apparently not very good, because the phase field model crack evolves earlier and completely breaks the sample, while the peridynamic model shows smaller damage values. However, if we impose a very small perturbation in the damage value of only one node we obtain damage evolution contour plots after 200 hours in Fig. 15a similar to the phase field results in Fig. 14. Since the model has a symmetrical geometry and we apply symmetric loads and constrains the numerical results must be symmetric; however, in some numerical codes due to the cumbersome computations a round-off error occurs in the computations which could lead to non-symmetric results. For further validation of the present method, two nodes x A ¼ ð0:2; 0Þ m and x B ¼ ð0:5; 0Þ m see (Fig. 12) are chosen to investigate the time history of the temperature and hydration degree. The time history of the temperature and hydration degree for nodes x A and x B are shown in Figs. 16a, b respectively. The results for both of the nodes are in an acceptable agreement with the phase field models in [29].
6.3 Example III: comparison of the crack shape in a concrete beam In this example, crack initiation and propagation in a non-reinforced concrete C20/25 beam is investigated by using the proposed peridynamic model. The results will be compared to the experimental tests and the phase field model simulation results [1]. As mentioned in [1], C20/25 is risky for the durability of the structures and it is prone to loss of performance due to cracking induced by early-age hydration. Damage and fracture may take place in a critical working condition as mentioned in [1]. Figure 17 shows the geometry  and the boundary conditions of the specimen. The central part of the H-shape specimen is bracing, while stiff steel tubes are arranged parallel and symmetrically to the long sides of the middle part. Four critical zones with a high risk of cracking emerges at four corners of the specimen as depicted in Fig. 17a. The material parameters in [1] are considered in this example and listed in Table 1. The initial temperature is assumed to be 20 C [1]. For the sake of simplicity, only half of the specimen is modelled in our simulations due to its symmetric shape. The steel tubes and plates are made of AISI 4130 Alloy steel; its properties are taken from Nguyen et al. [1] and are listed in   19 Example III: the distribution of the temperature T ( C) a using the present approach b using the phase field approach [1] and hydration degree a c using the present approach d using the phase field approach [1] after 24 h Fig. 20 Crack morphology obtained by a present approach b and c experiment in [1] steel tube are zero, while their y-displacements are free. It is assumed that there is a perfect cohesion between the concrete and the steel interfaces so that the displacements along x-directions at both surfaces BC and FG are fixed while their ydisplacements are free. In all other boundaries (AB, CD, FE, GH), both x and y-displacements are free. 4. The initial conditions T 0 ¼ 20 C and a 0 ¼ 0:01 are taken for the whole structure.
The slab is discretized uniformly by a grid spacing of Dx ¼ Dy ¼ 0:00425 m and a total number of 10599 nodes. The incremental time step is taken as Dt ¼ 10 s for a total of 86400 steps. The contour plots of the damage are depicted in Fig. 18 for three different time steps.
The comparison of the crack paths in Fig. 18a, c, e, and the results obtained by the phase field model in Fig. 18b, d, f, shows that there is a very good agreement between the peridynamic and phase field approaches. Moreover, the temperature field and hydration degree at time t ¼ 240 hours are represented in Fig. 19.
The comparison of Fig. 19a, c, and the temperature field and hydration degree contour-plots obtained by the phase field model in Fig. 19b, d, indicates that there is an acceptable agreement between the peridynamic and phase field approaches. In order to compare the crack paths quantitatively, one may analyse the crack angles as shown in Fig. 20. This angle for the experimental test [1] is 145:8 and 156:63 at corner 1 and 3 respectively (see Fig. 20). This angle in the present study is absolutely symmetric and is equal to 147 (see Fig. 20). The crack angle obtained using the proposed approach is in a good agreement with the experimental test.

Conclusions
In this paper, a new technique is developed in a peridynamic framework to investigate the complex fracture mechanism induced by the chemo-thermomechanical reactions at early ages in a cement based structure. Using this method, one is able to predict complex crack patterns without any a priori hypothesis on crack onset and geometry. The accuracy of the model is examined through problems such as chemo-thermo-mechanical analysis of a concrete slab, hydration reaction in a cement ring structure and crack propagation over a concrete beam. Hydration degree, temperature and crack paths are compared to the experimental data and to numerical solutions produced by the FEM, hybrid FEM and phase field methods. A proper agreement is achieved between all sets of results. The new computational method is capable of producing symmetric results. More importantly, the crack angle observed by the present method is in a good agreement with the experimental data.
Acknowledgements Authors acknowledge the support they received from MIUR under the research project PRIN2017-DEVISU and from University of Padua under the research projects BIRD2018 NR.183703/18 and BIRD2017 NR.175705/ 17.
Funding Open access funding provided by Politecnico di Milano within the CRUI-CARE Agreement.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.