Ordinary state-based peridynamic modelling for fully coupled thermoelastic problems

An ordinary state-based peridynamic model is developed for transient fully coupled thermoelastic problems. By adopting an integral form instead of spatial derivatives in the equation of motion, the developed model is still valid at discontinuities. In addition, the ordinary state-based peridynamic model eliminates the limitation on Poisson’s ratio which exists in bond-based peridynamics. Interactions between thermal and structural responses are also considered by including the coupling terms in the formulations. These formulations are also cast into their non-dimensional forms. Validation of the new model is conducted by solving some benchmark problems and comparing them with other numerical solutions. Thin plate and block under shock-loading conditions are investigated. Good agreements are obtained by comparing the thermal and mechanical responses with those obtained from boundary element method and finite element solutions. Subsequently, a three-point bending test simulation is conducted by allowing crack propagation. Then a crack propagation for a plate with a pre-existing crack is investigated under pressure shock-loading condition. Finally, a numerical simulation based on the Kalthoff experiment is conducted in a fully coupled manner. The crack propagation processes and the temperature evolutions are presented. In conclusion, the present model is suitable for modelling thermoelastic problems in which discontinuities exist and coupling effects cannot be neglected.


Undeformed
Deformed bond-based PD model is limited to have a fixed Poisson's ratio, i.e. 1/3 for two-dimensional problems and 1/4 for three-dimensional problems [48]. Considering this constraint, a generalized PD model which eliminates the aforementioned limitations is necessary for thermomechanical problems. The purpose of this paper is to develop a PD model overcoming the bond-based PD limitation on material properties for fully coupled thermoelastic problems. Therefore, ordinary state-based PD theory is employed and its formulations are derived in the realm of fully coupled thermomechanics. Then both dimensional and non-dimensional forms of the coupled equations in both deformation and thermal fields are provided and subsequently validated. Temperature changes and displacements of an isotropic thin plate and a block under various types of loading conditions are investigated. The validation is conducted via comparing the results with those from ANSYS and BEM solutions. The remarkable agreement indicates the capability of the proposed PD model to accurately predict the mechanical and thermal responses in a fully coupled manner. In the final part of this paper, crack propagation patterns are predicted for three-point bending test, Kalthoff problem. In addition, crack propagation for a plate with a pre-existing crack subjected to a pressure loading is conducted. The crack growth paths and temperature distributions are presented, and crack patterns are compared with the predictions from pure mechanical analysis. Hence, the effect of coupling term on crack growth is estimated and discussed.

Fully coupled peridynamic thermomechanics
The peridynamics falls into the category of non-local theories, i.e. particles interact with each other within a horizon, δ, as shown in Fig 1. Each material point is identified by its location represented by coordinate x in an undeformed state. The body region is R, and the interaction domain of material point x is denoted by H x . The other material points in H x , i.e. x , are called the family members of x. Furthermore, y and y represent the positions of x and x in a deformed configuration. In the deformed configuration, t − t is the force function exerted from point x to point x. Similarly, t − t is the force function acting at material point x from point x. In bond-based PD theory, the pairwise force functions are equal in magnitude and also parallel to their relative position in the deformed state, i.e. t//t and |t| = t [48,49]. Therefore, Poisson's ratio is 1/3 for two-dimensional (2D) problems and 1/4 for three-dimensional (3D) problems [15]. On the other hand, the state-based force functions are still parallel to their relative position in the deformed state, i.e. t//t . But they do not have to be equal in magnitude in ordinary state-based PD model, i.e. |t| = t or |t| = t . Thus, ordinary state-based formulation overcomes the limitation on Poisson's ratio. Hence, in this study ordinary state-based PD theory is adopted for fully coupled thermoelastic problems. Regarding the fully coupled thermomechanical analysis, the ordinary state-based PD model has not been explicitly provided, and only the expression in bond-based theory is available in the published literature [43,44,50]. Therefore, in this study, ordinary state-based fully coupled PD thermoelastic equations with explicit formulations of PD parameters are provided. Then these equations are cast into their corresponding non-dimensional forms.

Peridynamic thermoelastic formulations
The general form of fully coupled thermoelasticity is provided by Oterkus et al. [43]. The equation of motion can be shown as The integration on the right-hand side of Eq. (1) represents the total PD force density acting on point x. Thereby, b (x, t) represents the volumetric body force at material point x. In Eq. (1) t u − u, x − x, t and t u − u , x − x , t are PD force functions that belong to material point x and x . The ordinary state-based expressions for PD force functions are presented as provided in Chapter 4 of [49]. and where B 1 and B 2 are PD auxiliary parameters; those can be defined as for 2D for 3D (6) where T represents the temperature change of point x with respect to reference temperature, T = (x, t)− 0 . Similarly, T is the temperature change of point x . The dilatations, θ for point x and θ for point x , are defined as The PD bond stretch, s, between material points x and x can be defined as The PD auxiliary parameter, , is defined as The relationship between the PD material parameters, i.e. a, b and d, and classical material parameters are listed as provided in Chapter 4 of [49]; The equation of motion for ordinary state-based formulation including the effect of temperature change can be written as On the other hand, the heat conduction equation in the fully coupled thermomechanical PD model is [43,44] where h s (x, t) is the rate of heat generation per unit volume. In the above equation, κ is defined as PD microconductivity with its definitions listed as [37,39,43] The second term in Eq. (12) represents the deformation coupling effect on temperature. The time rate of change of stretch extension,ė x − x , can be defined aṡ The physical meaning and theory foundation of the PD thermal modulus are fully discussed by Oterkus et al. [43]. In this paper, the same derivation approach is adopted. The initial form of the ordinary state PD force function for point x is shown in Eqs. (2) and (4). In another form, being similar to the derivation conducted by Oterkus et al. [43], the PD force function can be divided into two parts as [43] The first part on the right-hand side includes only the structural deformation, and the second part is related to temperature effect. In Eq. (15), k is called the modulus state [51]; the term BT represents the effects of thermal state on deformation. Substituting the expression in Eq. (4) to Eq. (2) PD force function can be obtained as By plugging the dilatation term in Eq. (7a) the PD force function becomes After rewriting the PD force function by splitting into pure mechanical and thermal part Eq. (17) becomes By comparing Eqs. (15) and (18), local thermal modulus of bond between x and x can be obtained as with B = β y −y |y −y| [43]. Substituting the PD parameters provided in Eq. (10) as Then the explicit forms of the local thermal modulus, β, for different dimensions can be obtained as Furthermore, if the bond-based restriction, a = 0, is applied as explained in Chapter 4 of [49], Eq. (21) will reduce to its bond-based form, where β b = 2δbα or β b = 1/2 (cα) with c is the bond-based PD material constant [43]. The bond-based PD form of thermomechanical formulations is provided in Appendix A.

Non-dimensional form of peridynamics thermoelastic formulations
The governing equation can be put in a non-dimensional form by using non-dimensional variables [52]. Therefore, the fully coupled PD equations are cast into their non-dimensional forms by adopting the approach proposed by Sackman [53] and Oterkus et al. [43].
Regarding the heat conduction equation, the diffusivity is defined as the characteristic length/time quantity For the equation of motion, the characteristic length/time is the elastic wave speed where λ and μ are Lame's constants. Combining the characteristic length/timescale leads to characteristic length and time as follows As explained in [43] following non-dimensional forms can be used Length-related variables: Displacement: Stretch: Velocity-related variables: And the temperature: By using the above non-dimensional parameters and substituting the peridynamic parameters listed in Eqs. (10), (13) and (21) into Eqs. (11) and (12), the non-dimensional form of the fully coupled equations can be achieved by utilizing the non-dimensional parameters given in Eqs. (22)(23)(24)(25)(26)(27)(28)(29)(30) as 1D analysiṡ In the above equations, the non-dimensional variables are denoted with an over score. The parametersT anḋ T are the rate of temperature changes at material pointx andx , respectively. The non-dimensional coupling coefficient, , measures the strength of coupling effect on temperature distribution due to deformation. It can be defined as [43] = β 2 cl 0

Peridynamic failure criteria
Since PD equations are formulated without any spatial derivatives, PD theory is suitable to be applied for failure prediction. Once the stretch between material points exceeds the critical stretch value, s c , the bond will be broken and will be removed permanently. At the same time, the force between these two points becomes zero. The critical stretch value for bond failure is related to the critical energy release rate G c as provided in Chapter 6 of [49]; A history-dependent damage function χ x, x , t is implemented for each interaction between the material points. The value of the function χ x, x , t will be set to be zero when the bond is broken.
The local damage at a point represents the weighted ratio of the number of the broken interactions to the total number of interaction. Therefore, crack propagation path can be presented by the local damage value as [48];

Numerical implementation
The developed PD thermomechanical model is discretized for numerical implementation. Thus, the discretized form of equation of motion provided in Eq. (11) can be written as where i represents the point of interest, and j represents the family members of point i.
The discretized form of dilatation provided in Eq. (7a) becomes Here, N i represents the number of family members for point i.
On the other hand, the discretized form heat equation provided in Eq. (12) can be written as The local damage at point i provided in Eq. (37) can be written as [49] ϕ ( Explicit time integration is used to find the temperature, velocity and displacement profile at each time step [43,49]. The numerical procedure is provided in Fig. 2. In Fig. 2, N t represents the total number of time step, n represents the current number of the time step, and N node represents the total number of PD points.

Numerical results for verification problems
In this section, peridynamic simulations are conducted by using the proposed model. The validity of the fully coupled thermomechanical equations is established by comparing the PD solutions with previously considered BEM and ANSYS solutions. Firstly, a dimensionless isotropic plate is imposed with three types of loadings, i.e. pressure shock loading, thermal shock loading and their combination. The results from PD solution are in agreement with the ones from an existing BEM solution. As a next verification problem, a block is investigated with a temperature boundary condition. Good agreements are obtained by comparing with PD and ANSYS solutions. Hence, the non-dimensional form of the equations is validated for both twoand three-dimensional problems. Secondly, a dimensional isotropic square plate is separately subjected to a tension pressure shock loading and a combination of compression and tension pressure shock loading. The temperature and displacement responses from PD solution coincide very well with ANSYS solution. Thus, present PD model is thus validated via these numerical simulations. Subsequently, proposed model is further employed for failure prediction including fully coupled analysis. A plate with a pre-existing crack is studied under a pressure shock-loading condition. The temperature and structural responses without crack propagation are verified against ANSYS solution. Then, the crack propagation is simulated and its path is compared with the one from pure mechanical simulation. In this way, the coupling term effect on crack propagation is estimated and analysed.

Plate subjected to shock loading
The validity of the non-dimensional thermoelastic PD equations for 2D problems is established by constructing PD solutions for an existing BEM solution provided by Hosseini-Tehrani and Eslami [54]. The same geometry model, boundary conditions and loading conditions are adopted as in [54]. As shown in Fig. 3, a thin plate with a non-dimensional geometryL = 10,W = 10 and thicknessh = 1 is subject to a shock loading on the edge ofx = −L/2 and fixed on the edge ofx =L/2. The edges ofȳ =W /2 andȳ = −W /2 are traction free. Furthermore, atx = −L/2 the plate is subjected to temperature boundary condition and all other three edges are insulated. The Poisson's ratio is set to be 0.17. Regarding the PD discretization, the grid size is x = 0.05 and the horizon size is chosen as δ = 3.015 x. The uniform time step size is 5.0 × 10 −4 with total simulation timet total = 6. The boundary condition is implemented by using fictitious layers [43,44]. The applied loads are implemented on real boundary layer [49]. Since there is no heat source in this simulation case,h s = 0.
The initial conditions are:ūx withūx andūȳ representing the non-dimensional displacement components in the x and y directions, respectively. The shock-loading conditions are: Loading 2: Thermal shockP

Fig. 3 A thin plate subjected to shock-loading conditions
Loading 3: Combined pressure and thermal shock wheret represents the non-dimensional time. The applied pressure load is in the positive x direction, as illustrated in Fig. 3. Figures 4, 5 and 6 provide the thermal and mechanical responses att = 3 andt = 6 along the horizontal centreline of the plate for 3 different loading conditions. The coupling coefficient = 0 represents the uncoupled case, where the effect of deformation on temperature field is ignored. Figure 4a represents the temperature distribution when the plate is subjected to pure pressure shock loading (loading condition 1). As it can be seen from the figure, when = 0 no temperature change is observed. On the other hand, when = 0.1 both temperature drop and temperature rise are observed, which are induced by the applied pressure shock due to coupling effect. The magnitude of the temperature change is relatively small, within the range between -0.02 and 0.05. As the time progresses, the peak position of the temperature distribution moves towards the positive x direction. Figure 4b represents the dimensionless axial displacement along the horizontal centreline att = 3 andt = 6. The wave fronts at these two time points are clearly observed. There are slight differences between the displacement predictions from the coupled and uncoupled cases. As the time progresses, the difference becomes larger. Therefore, it could be inferred that due to coupling effect, the temperature change induced by deformation does affect the deformation. The same conclusion is obtained from the simulation cases with loading condition 2 as can be seen from Fig. 5. Due to the heating effect by the applied thermal loading, the plate experiences an expansion state. Subsequently, the tension loading creates a cooling effect. Therefore, when compared with the uncoupled case, relatively lower temperature change is observed in coupled case. Consequently, the lower temperature change gives rise to a smaller deformation response. This conclusion can also be applied loading condition 3, which results are presented in Fig. 6.
In conclusion, good agreements are obtained for three types of loadings. For both the thermal and mechanical fields, the results from ordinary state-based PD predictions agree well with those from BEM solutions obtained by Hosseini-Tehrani and Eslami [54]. Therefore, via these numerical simulations, the present nondimensional PD model is validated for 2D problems.  Figs. 7a and 8a, there is a remarkable agreement between PD and ANSYS solutions. When the plate is subjected to loading condition 1, temperature drop is observed due to tension loading as it can be seen from Fig. 7a. As the time progresses, temperature change increases to a final value of 6.5 K. When the plate is subjected to loading condition 2, temperature increases where there is local compression and temperature drops where there is local tension. Under both loading conditions displacement fields obtained from PD and ANSYS simulations match very well. Furthermore, it should be noted that even though the carbon steel has a relatively small coupling coefficient, i.e. = 0.002861, the generated temperature change due to mechanical shock loading is considerable. Therefore, if large strain rate exists, i.e. shock loading is applied, fully coupled thermomechanical analysis should be taken into consideration.

Block subjected to thermal loading
In order to validate the proposed PD model for 3D problems, a block subjected to temperature boundary condition is investigated. As shown in Fig. 9, the dimensionless length, width and height of the block are 5, 0.15 and 0.15, respectively. The Poisson's ratio is set as 0.33, and the coupling coefficient is 1.0. Regarding the PD model, the grid size is x = 0.0125 and the horizon is chosen as δ = 3.015 x. The integration time step size is π/4 × 10 −4 , and the total simulation time is π. On the other hand, directly coupled solid element-type SOLID 226 is adopted with mesh size of 0.05 and time step size of 0.02π in ANSYS model. The block is clamped at x =L. The block is gradually heated at x = 0, and all other surfaces are insulated. The temperature boundary condition is defined asT = sin(t), wheret is the dimensionless time. Fictitious layer is used to implement the boundary conditions [43].
The temperature distributions and horizontal displacements along the line of y =W /2 and z =H /2 are presented in Fig. 10a and b at dimensionless time oft = π/4, π/2, 3π/4, π, respectively. The results which are obtained from ANSYS solutions are also provided for comparison. It could be seen that both the temperature and displacement distributions match very well, indicating the capability of the derived PD formulations to accurately predict the thermal and mechanical responses for three-dimensional problems.
Derived formulations and explicit expressions of PD parameters including their dimensional and nondimensional forms are validated through these numerical simulations.

Numerical results for damage prediction
Peridynamics is a reformulation of classical equations that is better suited for modelling bodies with discontinuities. Therefore, analyses involving cracks are considered in this study by using the developed PD model. First, a three-point bending test is considered to predict the temperature distribution and deformed shape. Then a plate with a pre-existing central crack subjected to a pressure shock-loading condition is investigated. The temperature distributions and crack propagation paths are investigated. In addition, the crack propagation paths without considering coupling terms are also predicted for comparison. Finally, a numerical simulation based on Kalthoff experiment [55] is carried out. The crack path predicted by PD model is compared with the result of the experiment. The pressure shock loading is specified as where t 0 = 8 µs. The total node number in x or y direction for PD model is set to be 500 with a grid size of 0.2 mm. Thus, the critical stretch value s c is calculated as 0.0213 with G c being 42320 J/m 2 . The horizon is δ = 3.015 x. The time step size is set as 10 −9 s, and the total simulation time is 30 µs. On the other hand, directly coupled plane element-type PLANE223 is applied in the ANSYS model with the grid size of x = 0.00125 and time step size of 0.6 µs.
In order to better understand the existence of crack surface on the temperature and deformation field, initially failure is not allowed. The horizontal displacement predictions at different time steps are shown in Figs. 15, 16 and 17. It is observed that the peridynamic results coincide very well with ANSYS solutions. The displacement distribution along x-axis propagates uniformly in the vertical direction before the elastic wave reaches the crack, as shown in Fig. 15. After the elastic wave reaches the crack surface, displacements become non-uniform due to the discontinuity of the displacements along the crack surface as it can be seen from Figs. 16 and 17.   19 and 20 present the induced temperature distributions due to applied loading. It is observed that the peridynamic results agree very well with ANSYS solutions. The temperature distributions along x-axis propagate almost uniformly in the vertical direction before the thermal wave reaches the crack, as shown in Fig. 18. After the thermal wave reaches or passes the crack, thermal waves get disturbed by the existence of the crack. Higher temperature drop is observed near the crack tip as can be seen from Figs. 19 and 20. This also indicates the stress concentration near the crack tips. Besides, the region near the crack tips is under tension, indicating a tendency of crack growing in the vertical direction. The crack surfaces experience local compression. Therefore, temperature rise is observed in these regions as it can be seen from Fig. 20.
After verification of temperature and deformation field for a plate with a pre-existing crack, as a next example, crack propagation is allowed. The crack configurations at different time steps are provided in Figs. 21, 22 and 23. Crack propagation patterns are compared with coupled and uncoupled cases. Temperature distributions at corresponding stages from coupled cases are also provided.
In the early stage, the cracks grow in similar patterns for both simulation cases. Crack starts to propagate earlier for the uncoupled case (Fig. 21). Cracks both begin to propagate at around 16 µs. Up to 28 µs, the cracks propagate in a self-similar manner for both coupled and uncoupled cases. Cracks start branching at around 28 µs (Fig. 22) and split into visible asymmetrical branches (Fig. 23). In addition, the branches for uncoupled case grow faster than the coupled case.
For coupled case, it is clear that before 16 µs the temperature distribution is the same with the one obtained from the simulation without crack propagation. However, temperature drops at the crack tips move as the crack propagates (Figs. 22c, 23c). The cooling region at the crack tips creates temperature-induced local compression at these regions. Furthermore, the temperature rise around the crack surfaces creates local tension against the compression created by the pressure shock loading. In conclusion, the induced temperature change due to deformation influences the crack growth in the opposite direction against the applied mechanical loading, leading to a reduced degree of crack propagation response. Hence, different crack pattern from uncoupled simulation is obtained.
In conclusion, if shock loadings are applied, large strain rates are created, and thus the coupling term should be considered for more accurate crack propagation predictions.

Kalthoff problem simulation
Kalthoff and Winkler [57] and Kalthoff [55,58] performed a series of experiments where pre-notched plates were subjected to dynamic shear loads. In the experiments, a cylindrical projectile impacted on the notched side of the plate with a constant velocity V 0 parallel to the axis of the notch. The pre-existing crack in the upper half steel plate was observed to grow in an angle of approximately 70 • counter clockwise with the notch axis. The failure is in a brittle fracture mode under lower strain rate. These experiments have been successfully simulated by numerical methods such as phase-field simulation [59], finite element method [60]. In addition, Silling [61] and Dipasquale et al. [62] used PD to numerically simulate the Kalthoff problem within the realm of mechanical field. In this section, a fully coupled thermomechanical simulation is conducted based on the Kalthoff experiment. The problem is symmetric so that only the upper half plate is modelled. As shown in Fig. 24, a square plate is modelled with L = W = 100 mm and its thickness is 1 mm. A pre-existing crack of length being 50 mm is located above the x-axis with the distance of 25 mm. Due to the symmetric conditions, the lower horizontal edge of the plate is fixed in the y direction, i.e. u y (x, y = 0, t) = 0. The other boundaries are free. All the boundaries are thermally insulated. The impact is simulated by imposing a constant velocity to the nodes on the left surface between the crack and the lower horizontal boundary in the PD model. The velocity is parallel to the x direction, and its magnitude is where |v 0 | represents the magnitude of the applied velocity with v 0 = 16.5 m/s and t 0 = 1 µs [59]. The properties of the elastic material are E = 190 GPa, ρ = 8000 kg/m 3 , ν = 0.3, c v = 477 J/ (kgK), α = 17.6 × 10 −6 K −1 , k T = 16.2 W/ (mK). The critical energy release rate is G c = 2.213 × 10 4 J/m 2 . The reference temperature is 0 = 285 K. As to the PD discretization, the mesh size is x = 0.0005 m and the horizon is chosen as δ = 3.015 x. The time step size is 0.01 µs, and the total simulation time is 90 µs. The critical stretch value is calculated as s c = 0.0103. The crack is observed to propagate at t = 20 µs. The crack evolution at different times is provided in Fig. 25. The angle between the crack path and the positive x direction is observed to be 68 • , which is close to the It can be observed temperature rises near the crack and temperature drops in the crack, which agrees with the conclusions drawn in the last two simulation cases.

Conclusion
In this study, fully coupled thermoelastic equations in ordinary state-based peridynamic theory are provided, including their non-dimensional forms. To verify the PD model, some benchmark problems are solved by using both peridynamics and FEM solutions. The good agreement between PD and other methods indicates the validity of the proposed PD model. Finally, crack propagation patterns are predicted for three-point bending test, Kalthoff problem and a plate with a pre-existing crack subjected to a pressure load. The following conclusions are drawn 1. The ordinary state-based fully coupled thermomechanical PD model introduced in this paper is capable of predicting temperature and displacement responses accurately both for dimensional and for nondimensional problems. 2. When shock loadings are applied, the coupling effect on displacements and temperature should be taken into consideration for more accurate results. 3. The coupling terms do have an effect on crack propagation when shock loadings are applied. Therefore, fully coupled analysis is necessary in such cases.

Acknowledgements
The authors gratefully acknowledge financial support from China Scholarship Council (CSC No. 201506230126) and University of Strathclyde.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Appendix A
A.1 Ordinary state-based thermal modulus for a material point The thermal modulus for a material point can be obtained by using the PD local thermal modulus formulation provided in Eq. (21) as follows: When small deformation approximation is adopted, provided in Eq. (9) results in = 1. Consequently, the integration term in Eq. (21) becomes with V H x representing the integration volume. For 2D problems, the integration domain is a circle disc with its radius, δ, and thickness, h. On the other hand, the integration domain for 3D problems is a sphere; thus, the volume in Eq. (A.1) can be calculated as with ξ = x − x (Fig. 1). In Eq. (A.3) ϑ represents the bond angle with respect to x-axis and φ represents the bond angle with respect to the x-y plane for 3D case Recalling the classical material constants As it can be seen from the above formulations, bond-based PD has a limitation on Poisson's ratio [48].

Equation of motion
In bond-based PD by applying the restriction a = 0 in Chapter 4 of [49], the PD auxiliary parameters A and B provided in Eqs. (4) and (5)

Heat equation
The form of heat equation in bond-based peridynamics remains the same as Eq. (12) [43] ρc vṪ ( with bond-based thermal modulus provided in Eq. (A.10)