Isogeometric analysis of coupled thermo-elastodynamic problems under cyclic thermal shock

Abstract

The isogeometric analysis (IGA) method was extended for the solution of the coupled thermo-elastodynamic equations. The dimensionless formulation was accepted in discretization of the uncoupled and coupled thermoelasticity equations and the Generalized Newmark method was used in the time integration procedure. First, the performance of the proposed method was verified against a two-dimensional benchmark example subjected to constant thermal shock with available exact analytical solutions. Then a two-dimensional half-space benchmark example under thermal shock was solved. Finally, cyclic thermal shock (CTS) loading was applied on the half-space problem. The results dedicated that IGA can be used as a suitable approach in the analysis of the general thermomechanical problems.

Keywords

isogeometric analysis / coupled thermo-elastodynamic / dynamic analysis / generalized newmark / cyclic thermal shock

Introduction

The solution of the thermomechanical problems has been widely used in engineering. Thermomechanic means considering effect of the thermal loading on a mechanical body. This problem was defined in coupled and uncoupled states. The studies of initial boundary value problems generated by linear coupled thermoelasticity equations began in 1957 [¹]. Danilovskaya [²] worked on an analytical solution for thermal stresses in an elastic half-space due to a sudden heating of its boundary. Sternberg and Chakravorty [³] studied inertia effects in a transient thermoelastic problem. Boley and Tolins [⁴] obtained the analytical solution for the transient coupled thermoelastic boundary value problems in the half-space. Nickell and Sackman [⁵] presented an approximate solution in linear coupled thermoelasticity. In the numerical field, boundary element method (BEM) has been reported for the coupled thermoelasticity problem [⁶–⁸]. Tehrani and Eslami [⁹,¹⁰] employed the BEM to investigate the effect of the coupling and inertia terms in dynamical thermal loading problems. In finite element method, KC and Kim [¹¹] worked on the evaluation of the non-singular T-stress and mixed-mode stress intensity factors in functionally graded materials under steady-state thermal loads via interaction integral.

The isogeometric analysis (IGA) is a newly numerical method introduced by Hughes et al. [¹²] based on the converting suitable meshes from developments in Computer Aided Design (CAD) technology for numerical analysis to able to define complex geometries. For this purpose, splines and their new versions, i.e., non-uniform rational B-splines (NURBS) and T-splines, have been usually employed. Very promising results were reported in current researches [¹³–²⁴] where the idea was employed in several advanced problems in the area of computational mechanics. In the IGA geometry construction and boundary definition are flexible; hence, the problem discretization and the boundary satisfaction are accurate [¹²,²²]. In Ref. [¹⁶], IGA is extended for analysis of thin shells using polynomial splines over hierarchical T-meshes (PHT-splines) to satisfy C1 continuity for the Kirchhoff-Love theory and to avoid rotational degrees of freedom. Nguyen-Thanh et al. [¹⁷] used an alternative method to NURBS-based isogeometric analysis based on polynomial splines to improve efficiency of T-meshes for local refinement in the IGA.

In this paper, IGA is proposed for solution of the uncoupled and coupled planar dynamic thermoelasto-dynamic equations. The coupled dynamical system of equations obtained from the IGA discretization was solved by the Generalized Newmark method in time domain. The NURBS functions does not satisfy the Kronecker-delta property. Hence, it is require some post processing to get true values for the unknowns. Here, least squares minimization procedure [²⁵] was used in IGA formulation to apply the essential boundary conditions directly. Performance of the method was verified by solving and comparison of a benchmark example with available exact analytical solutions under constant shock. Then a thermal shock with step function was applied and the results of the uncoupled and coupled systems were compared. When cyclic thermal shock was applied on a solid domain, after a certain number of thermal shock damage will appear by a decrease in Young’s modulus of elastic materials. Hence, in this paper, the CTS loading was applied in thermoelasticity problem to better show the capability and efficiency of the proposed IGA to analysis of the thermomechanical problems.

This paper is organized as follows. The dimensionless governing equations of coupled thermoelasticity problem are discretized in section 2. In section 3, some examples are solved to verify the proposed method. Then, some concluding remarks are addressed in section 4.

Governing equations and discretization

In this section, first, the governing coupled thermo-elastodynamic equations are described and discretized. The least square minimization approach is introduced to directly impose the essential boundary conditions in IGA. Finally, the Generalized Newmark (GN) time integration procedure is applied to discretize the dynamic coupled equations in time domain.

Governing equations

In absence of body forces and heat flux, the governing equations for the coupled thermo-elastodynamic in the time domain can be written as follows [¹⁰]:

(1)(\lambda +\mu ){{\displaystyle u}}_{j,ij}+{{\displaystyle \mu u}}_{i,jj}-{{\displaystyle \rho \ddot{u}}}_i-{{\displaystyle \gamma T}}_{,i}=\mathrm{0},

(2){{\displaystyle k}}_{ij}{{\displaystyle T}}_{,ii}-{{\displaystyle \rho c\dot{T}}}_e-{{\displaystyle \gamma {{\displaystyle T}}_{\mathrm{0}}\dot{u}}}_{j,j}=\mathrm{0},

where {{\displaystyle u}}_i is the displacement vector components, \rho is solid density, T is current absolute temperature, {{\displaystyle T}}_{\mathrm{0}} is the initial reference temperature, {{\displaystyle c}}_e is the specific heat capacity at constant deformation expresses as (energy).(unit mass .{}^{\circ }C)^-1, {{\displaystyle k}}_{ij} is the coefficient of thermal conduction for a general anisotropic material, \gamma is stress–temperature modulus. The stress tensor components was defined by the Neumann-Duhamel relation as,

(3){{\displaystyle \sigma }}_{ij}=\mu ({{\displaystyle u}}_{i,j}+{{\displaystyle u}}_{j,i})+[{{\displaystyle \lambda u}}_{k,k}-\beta \theta ]{{\displaystyle \delta }}_{ij},

where \beta =\alpha (\mathrm{3}\lambda +\mathrm{2}\mu ) in which \alpha is the coefficient of thermal expansion.

It is convenient to introduce dimensionless variables [²⁶] as follows,

(4)\widehat{x}={\displaystyle \frac{x}{a}},\widehat{t}={\displaystyle \frac{{{\displaystyle tc}}_{\mathrm{1}}}{a}},\widehat{\sigma }={\displaystyle \frac{{{\displaystyle \sigma }}_{ij}}{{{\displaystyle \gamma T}}_{\mathrm{0}}}},{{\displaystyle \widehat{u}}}_i={\displaystyle \frac{(\lambda +\mathrm{2}\mu ){{\displaystyle u}}_i}{{{\displaystyle a\gamma T}}_{\mathrm{0}}}},\widehat{T}={\displaystyle \frac{T-{{\displaystyle T}}_{\mathrm{0}}}{{{\displaystyle T}}_{\mathrm{0}}}},

where \alpha =k/{{\displaystyle \rho c}}_e{{\displaystyle C}}_{\mathrm{1}} is the dimensionless unit length and {{\displaystyle C}}_{\mathrm{1}}=\sqrt{(\lambda +\mathrm{2}\mu )/\rho } is the velocity of propagation of the longitudinal wave. With above dimensionless definitions Eqs. (1) and (2) take the form (with dropping the hat for convenience)

(5){\displaystyle \frac{\mu }{\lambda +\mathrm{2}\mu }}{{\displaystyle u}}_{j,ij}+{\displaystyle \frac{\lambda +\mu }{\lambda +\mathrm{2}\mu }}{{\displaystyle u}}_{i,jj}-{{\displaystyle \ddot{u}}}_i-{{\displaystyle T}}_{,i}=\mathrm{0},

(6){{\displaystyle T}}_{,ii}-\dot{T}-{\displaystyle \frac{{\gamma }^{\mathrm{2}}{{\displaystyle T}}_{\mathrm{0}}}{{{\displaystyle \rho c}}_e(\lambda +\mathrm{2}\mu )}}{{\displaystyle \dot{u}}}_{j,j}=\mathrm{0.}

The imposed boundary conditions are

{{\displaystyle u}}_i={{\displaystyle \overline{u}}}_i,\hspace{0.17em}\mathrm{o}\mathrm{n}\hspace{0.17em}\Gamma ={{\displaystyle \Gamma }}_u,

{\sigma }_{ij}{n}_j={\overline{t}}_i,\hspace{0.17em}\mathrm{o}\mathrm{n}\hspace{0.17em}\Gamma ={{\displaystyle \Gamma }}_t,

q_i{n}_i={{\displaystyle \overline{q}}}_i,\hspace{0.17em}\mathrm{o}\mathrm{n}\hspace{0.17em}\Gamma ={{\displaystyle \Gamma }}_q,

(7){{\displaystyle u}}_i={{\displaystyle \bar{u}}}_i,\hspace{0.17em}on\hspace{0.17em}\Gamma ={{\displaystyle \Gamma }}_u,

where \Gamma ={{\displaystyle \Gamma }}_t\cup {{\displaystyle \Gamma }}_u={{\displaystyle \Gamma }}_q\cup {{\displaystyle \Gamma }}_{\theta } and \bar{u},\bar{t},\bar{T} and \bar{q} are the prescribed displacement, traction, temperature and heat flux vector on the boundaries, respectively.

Discretization

Applying the weighted residual integral to the coupled dimensionless thermo-elastodynamic Eqs. (5) and (6) with respect to the weighting functions W, the formal Galerkin approximations reduce to,

(8){\displaystyle \underset{V(e)}{\int }\left({\displaystyle \frac{\mu }{\lambda +\mathrm{2}\mu }}{{\displaystyle u}}_{j,ij}+{\displaystyle \frac{\lambda +\mu }{\lambda +\mathrm{2}\mu }}{{\displaystyle u}}_{i,jj}-{{\displaystyle \ddot{u}}}_i-{{\displaystyle T}}_{,i}\right)}{{\displaystyle \mathit{W}}}_i\mathrm{d}V=\mathrm{0},i=\mathrm{1},\mathrm{2},\mathrm{...},ns,

(9){\displaystyle \underset{V(e)}{\int }\left(T_{ii}-\dot{T}-{\displaystyle \frac{{\gamma }^2{T}_0}{\rho c_e(\lambda +\mathrm{2}\mu )}}{\dot{{\displaystyle u}}}_{j,j}\right)}{{\displaystyle W}}_i\mathrm{d}V=\mathrm{0},i=\mathrm{1},\mathrm{2},\mathrm{...},ns,

where ns is the number of the NURBS shape functions of the element (e) and {{\displaystyle \mathit{W}}}_i is the component of the vector \mathit{W}.

Thus the coupled u-T formulation can be written in tensor index form:

(10)\left[\begin{array}{ll}{{\displaystyle M}}_{KijL}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}\end{array}\right]\left[\begin{array}{l}{{\displaystyle \ddot{u}}}_{Lj}\\ {{\displaystyle \ddot{T}}}_N\end{array}\right]+\left[\begin{array}{ll}\mathrm{0}& \mathrm{0}\\ {{\displaystyle G}}_{LjM}& {{\displaystyle S}}_{MN}\end{array}\right]\left[\begin{array}{l}{{\displaystyle \dot{u}}}_{Lj}\\ {{\displaystyle \dot{T}}}_N\end{array}\right]+\left[\begin{array}{ll}{{\displaystyle K}}_{KijL}& -{{\displaystyle Q}}_{KiN}\\ \mathrm{0}& {{\displaystyle H}}_{MN}\end{array}\right]\left[\begin{array}{l}{{\displaystyle u}}_{Lj}\\ {{\displaystyle T}}_N\end{array}\right]=\left[\begin{array}{l}\mathrm{0}\\ \mathrm{0}\end{array}\right],

where,

(11)\begin{array}{l}{{\displaystyle M}}_{KijL}={{\displaystyle \delta }}_{ij}\left({{\displaystyle \int }}_{\Omega }{{\displaystyle N}}_K^u{{\displaystyle N}}_L^{u}d\Omega \right),\\ {{\displaystyle G}}_{LjM}=C{{\displaystyle \int }}_{\Omega }{{\displaystyle N}}_{L,j}^u{{\displaystyle N}}_M^{T}d\Omega ,{{\displaystyle S}}_{MN}={{\displaystyle \int }}_{\Omega }{{\displaystyle N}}_M^T{{\displaystyle N}}_N^{T}d\Omega ,\\ {{\displaystyle K}}_{KijL}={\displaystyle \frac{\mathrm{1}}{\lambda +\mathrm{2}\mu }}{{\displaystyle \int }}_{\Omega }{B}^{T}DBd\Omega ,{{\displaystyle H}}_{MN}={{\displaystyle \int }}_{\Omega }{{\displaystyle N}}_{M,i}^T{{\displaystyle N}}_{N.i}^{T}d\Omega ,{{\displaystyle Q}}_{KiN}={{\displaystyle \int }}_{\Omega }{{\displaystyle N}}_{L,j}^u{{\displaystyle N}}_M^{T}d\Omega ,\end{array}

where \mathrm{0}\le C\le \mathrm{1} is the coupled parameters, the solid displacement components {{\displaystyle u}}_i and the temperature T can be approximated using NURBS shape functions and nodal values.

Least squares minimization

The least-squares minimization procedure was applied for imposing the essential boundary condition directly. In this method, it is not required to find the penalty coefficient or Lagrange multiplier for imposing the essential boundary conditions. The least squares error functional J for a set of interpolating points was picked on the boundary of the domain x^p,p\in P where P is the index of influenced nodes with the boundary nodes:

(12)\mathit{J}={\displaystyle \frac{\mathrm{1}}{\mathrm{1}}}\underset{p\in P}{\Sigma }{\Vert {\displaystyle \sum _{A\in A^p}{{\displaystyle N}}_A(S^p){{\displaystyle q}}_A-{{\displaystyle u}}_d(x^p)}\Vert }^{\mathrm{2}},

where A^p represents the set of basis functions for the p^th interpolation point. Minimization of above functional leads to the linear system as follows:

(13){\displaystyle \frac{\partial J}{\partial {{\displaystyle q}}_B}}=\underset{p\in P}{\Sigma }\left({\displaystyle \sum _{A\in A^p}{{\displaystyle N}}_A(s^p){{\displaystyle q}}_A-{{\displaystyle u}}_d(x^p)}\right){{\displaystyle N}}_B(s^p)=\mathrm{0}\forall i\in A^B=\{\cup A^p,p\in P\}.

Leads to

(14)\sum _{p\in P}{\displaystyle \sum _{A\in A^p}{{\displaystyle N}}_B(s^p)N_A(s^p)q_A}=\sum _{p\in P}{u}_d(x^p){{\displaystyle N}}_B(s^p)\forall B\in A^B\mathrm{.}

Rearrange N and q in matrix form

(15)\underset{p\in P}{\Sigma }\left[\begin{array}{lll}{{\displaystyle N}}_{\mathrm{1}}^p{{\displaystyle N}}_{\mathrm{1}}^p& \cdots & {{\displaystyle N}}_{\mathrm{1}}^p{{\displaystyle N}}_n^p\\ ⋮& \ddots & ⋮\\ {{\displaystyle N}}_n^p{{\displaystyle N}}_{\mathrm{1}}^p& \dots & {{\displaystyle N}}_n^p{{\displaystyle N}}_n^p\end{array}\right]\left[\begin{array}{lll}{{\displaystyle q}}_x^{\mathrm{1}}& {{\displaystyle q}}_y^{\mathrm{1}}& {{\displaystyle q}}_z^{\mathrm{1}}\\ ⋮& ⋮& ⋮\\ {{\displaystyle q}}_x^n& {{\displaystyle q}}_y^n& {{\displaystyle q}}_z^n\end{array}\right]=\underset{p\in P}{\Sigma }\left[\begin{array}{lll}{{\displaystyle N}}_x^p{{\displaystyle N}}_{\mathrm{1}}^p& {{\displaystyle N}}_y^p{{\displaystyle N}}_{\mathrm{1}}^p& {{\displaystyle N}}_z^p{{\displaystyle N}}_{\mathrm{1}}^p\\ ⋮& ⋮& ⋮\\ {{\displaystyle u}}_x^p{{\displaystyle N}}_n^p& {{\displaystyle u}}_y^p{{\displaystyle N}}_n^p& {{\displaystyle u}}_z^p{{\displaystyle N}}_n^p\end{array}\right],

where n equals the total number of base functions involved in the boundary displacements. Enough points have to be taken for the system to be invertible. A direct solver for symmetric definite matrices was used to compute the displacements of boundary nodes {{\displaystyle q}}_A.

Generalized Newmark time integration procedure

Generalized Newmark (GN) is a well-known method from the Newmark family which is used for time discretization of the coupled equations [²⁷]. According to GN, problem variables and their temporal derivatives in the time interval of [t_n, t_n+1] is given as:

(16)\{\begin{array}{l}{{\displaystyle \ddot{u}}}_{n+\mathrm{1}}^{(k)}={{\displaystyle \ddot{u}}}_n+{{\displaystyle \Delta \ddot{u}}}_n\Rightarrow {{\displaystyle \ddot{u}}}_{n+\mathrm{1}}^{(k-\mathrm{1})}+{\Delta \ddot{u}}^{(k)}={{\displaystyle \ddot{u}}}_n+{{\displaystyle \Delta \ddot{u}}}_n\Rightarrow {{\displaystyle \Delta \ddot{u}}}_n={{\displaystyle \ddot{u}}}_{n+\mathrm{1}}^{(k-\mathrm{1})}+{\Delta \ddot{u}}^{(k)}-{{\displaystyle \ddot{u}}}_n\\ {{\displaystyle \dot{u}}}_{n+\mathrm{1}}={{\displaystyle \dot{u}}}_{n+\mathrm{1}}+\Delta t{{\displaystyle \ddot{u}}}_n+{{\displaystyle \beta }}_{\mathrm{1}}\Delta t\left({{\displaystyle \ddot{u}}}_{n+\mathrm{1}}^{(k-\mathrm{1})}+{\Delta \ddot{u}}^{(k)}-{{\displaystyle \ddot{u}}}_n\right)\\ {{\displaystyle u}}_{n+\mathrm{1}}={{\displaystyle u}}_n+\Delta t{{\displaystyle \dot{u}}}_n+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\Delta t}^{\mathrm{2}}{{\displaystyle \ddot{u}}}_n+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{{\displaystyle \beta }}_{\mathrm{2}}{\Delta t}^{\mathrm{2}}\left({{\displaystyle \ddot{u}}}_{n+\mathrm{1}}^{(k-\mathrm{1})}+{\Delta \ddot{u}}^{(k)}-{{\displaystyle \ddot{u}}}_n\right)\end{array},

and

(17)\{\begin{array}{l}{{\displaystyle \dot{T}}}_{n+\mathrm{1}}^{(k)}={{\displaystyle \dot{T}}}_{n+\mathrm{1}}^{(k-\mathrm{1})}+{\Delta \dot{T}}^{(k)}={{\displaystyle \dot{T}}}_n+{{\displaystyle \Delta \dot{T}}}_n\Rightarrow {{\displaystyle \Delta \dot{T}}}_n={{\displaystyle \dot{T}}}_{n+\mathrm{1}}^{(k-\mathrm{1})}+{\Delta \dot{T}}^{(k)}-{{\displaystyle \dot{T}}}_n\\ {{\displaystyle T}}_{n+\mathrm{1}}={{\displaystyle T}}_n+\Delta t{{\displaystyle \dot{T}}}_n+{{\displaystyle \bar{\beta }}}_{\mathrm{1}}\Delta t\left({{\displaystyle \dot{T}}}_{n+\mathrm{1}}^{(k-\mathrm{1})}+{\Delta T}^{(k)}-{{\displaystyle \dot{T}}}_n\right)\end{array},

where k is the iteration number in the Newton-Raphson iteration and Dt is the length of time step.

The approximation is stable unconditionally when [²⁸]:

(18){{\displaystyle \beta }}_{\mathrm{2}}\ge {{\displaystyle \beta }}_{\mathrm{1}}\ge {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathrm{a}\mathrm{n}\mathrm{d}{{\displaystyle \overline{\beta }}}_{\mathrm{1}}\ge {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}.

By considering the governing equation in t_n+1 and implementing the variable at t_n+1 in terms of the variable at t_n, the governing equation lead to:

(19)\left[\begin{array}{ll}M+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{{\displaystyle \beta }}_{\mathrm{2}}{\Delta t}^{\mathrm{2}}K& -{{\displaystyle \bar{\beta }}}_{\mathrm{1}}\Delta tQ\\ -{{\displaystyle \bar{\beta }}}_{\mathrm{1}}\Delta tG& -S-{{\displaystyle \bar{\beta }}}_{\mathrm{1}}\Delta tH\end{array}\right]\left[\begin{array}{l}{\Delta \ddot{u}}^{(k)}\\ {\Delta \dot{T}}^{(k)}\end{array}\right]=\left[\begin{array}{l}{F}^u\\ -F^T\end{array}\right],

where the equivalent forces are

(20)F^T=-G\left({{\displaystyle \dot{u}}}_{n+\mathrm{1}}+\Delta t{{\displaystyle \ddot{u}}}_n+{{\displaystyle \beta }}_{\mathrm{1}}\Delta t\left({{\displaystyle \ddot{u}}}_{n+\mathrm{1}}^{(k-\mathrm{1})}-{{\displaystyle \ddot{u}}}_n\right)\right)-{{\displaystyle S\dot{T}}}_{n+\mathrm{1}}^{(k-\mathrm{1})}-H\left({{\displaystyle T}}_n+\Delta t{{\displaystyle \dot{T}}}_n+{{\displaystyle \bar{\beta }}}_{\mathrm{1}}\Delta t\left({{\displaystyle \dot{\theta }}}_{n+\mathrm{1}}^{(k-\mathrm{1})}-{{\displaystyle \dot{T}}}_n\right)\right).

Results and discussion

A half-space subjected to constant thermal shock

A benchmark problem of a square plate with available exact analytical solution subjected to constant thermal shock [²⁹] was considered as a first example for verification of the implemented Matlab code for the solution of uncoupled and coupled thermo-elastodynamic problems as shown in Fig. 1. The dimensionless square plate length is l=10 and the plane stress condition was considered. Left edge (x=0) was subjected to constant temperature T= 1 and the other edges were thermally isolated and are traction-free. Table 1 shows the material properties defined in the solution. For discretization of time domain of coupled and uncoupled systems by generalized Newmark procedure, the equal time step \Delta t=\mathrm{0.01} for the whole time simulation (\mathrm{0}\le t\le \mathrm{2}) was used.

The results were obtained for a point located at dimensionless length x=1 (which is the location of the elastic wave front at the dimensionless time t=1) along the axis of symmetry of the plate (x-axis) and compared with exact analytical solution. Figure 2 shows a comparison of the dimensionless temperature at dimensionless length x=1 for the coupled (C=1) and uncoupled (C=0) systems versus dimensionless time. The analytical solution is only available for the uncoupled system obtained from [²⁶]. It is indicated that the temperature results have good agreement with the analytical results. It is observed that for initial times the heat wave has small effects on temperature gradient at x=1 but at time about t=0.1 the heat wave front reach to this point and the temperature starts to increase considerably. The temperature in uncoupled system was increased from t=0.1 to reach the dimensionless temperature value about 0.62 at the end of simulation time t=2, but considering the coupled system, the temperature of point at x=1 was increased until time t=0.7 to value about 0.35 and then it was reduced slightly to reach value about 0.26 at time t=1. Then, with the same behavior as uncoupled system, the dimensionless temperature at point x=1 was increased to reach the maximum value about 0.42 at end of time simulation t=2. It is clear that considering damping effect parameter C in the coupled equations, the dimensionless temperature at the end of simulation, t=2, has smaller value in comparison with final temperature obtained for the uncoupled system.

The contours of the temperature distribution on problem domain for coupled and uncoupled systems versus time were compared in Fig. 3. The contours were plotted for first (t=0.01), middle (t=1) and end (t=2) of the simulation times. It is clear that the results of the temperature distributions are symmetric with respect to the horizontal center line of plate. It can be observed the temperature distribution in the problem domain in the uncoupled system is higher than the coupled system in which verify the same conclusion obtained from Fig. 2. This can be inferred from wider distribution of the highest value (red color) for the temperature in the uncoupled system at all of the simulation times. Another phenomenon that can be seen from the figures is the difference in the shape of heat front wave propagated in the coupled and uncoupled systems. The heat wave propagates uniformly in uncoupled system whereas the heat wave propagates in a non-uniform manner for the uncoupled system.

The results of the dimensionless axial displacement at point x=1 of the uncoupled and coupled systems with analytical solution versus time were compared in Fig. 4. As it can be seen, the thermal loading has not much effect on axial displacement of point at x=1 until time about t=0.2; Then it was increased due to temperature effects to reach to the maximum value of axial displacement at time t=1. Then the axial displacement at this point was decreased. At time about t=1.7 the dimensionless axial displacement of point at x=1 reached zero value. It means at this time the point does not possess displacements due to thermal loading. Then the direction of the axial displacement at x=1 was changed. It should be noted that at time t=1 the axial displacement has a shocked variation. As a conclusion from Figs. 2‒4, the temperature and the axial displacement of the uncoupled system has a higher value than those for the coupled system. It should be noted that the coarse elements and large time steps is not capable to represent exactly the sharp variations of displacement. The error in temperature and axial displacement distribution increase with the enhancement of the distance from the free edge x = 0; where the temperature shock was applied. Using smaller time step and finer elements will increase the accuracy in detection of any shock in the temperature and displacement.

A half-space subjected to thermal shocks

A square half-space plate of isotropic and homogeneous material was considered as second example according to Fig. 5. The dimensionless size of the plate is l=10 with unit thickness that leads to plane stress condition. The simulations for the coupled and uncoupled systems were done for 8 seconds with equal time steps \Delta t=\mathrm{0.01}. The material properties and constants are similar to previous example. The mechanical and thermal boundary conditions are

(21)\begin{array}{l}{{\displaystyle u}}_{\mathrm{1}}={{\displaystyle u}}_{\mathrm{1}}=\mathrm{0}\hspace{0.17em}\mathrm{a}\mathrm{t}\hspace{0.17em}x=\mathrm{10};{{\displaystyle \tau }}_i=\mathrm{0}\hspace{0.17em}\mathrm{a}\mathrm{t}\hspace{0.17em}y=\pm \mathrm{5},\\ {{\displaystyle T}}_{,\mathrm{1}}=\mathrm{0}\hspace{0.17em}\mathrm{a}\mathrm{t}\hspace{0.17em}x=\mathrm{10};{{\displaystyle T}}_{,\mathrm{2}}=\mathrm{0}\hspace{0.17em}\mathrm{a}\mathrm{t}\hspace{0.17em}y=\pm \mathrm{5},\end{array}

where, {{\displaystyle \tau }}_i are the traction components. The edge x=\mathrm{0} is kinematically free and thermally exposed to a thermal shock with the known equation as

(22)T(t)=\mathrm{5}t\exp (-\mathrm{2}t),

where, t is the dimensionless time. Then thermal shock T(t) was applied on the left edge (x=\mathrm{0})\tau (t)=\mathrm{0} on x=\mathrm{0}. Temperature shock pattern is in the form of heat input was applied in the positive x-direction (See Fig. 6).

The uncoupled and coupled equations were considered and solved. Figures 7 and 8 show the comparison of the temperature and axial displacement at the center of the plate along the x-axis versus dimensionless time t=2, t= 4, t= 6, and t=8.

The maximum value of temperature was decreased versus time at the center of the plate under thermal shock for the uncoupled and coupled systems as indicated in Fig. 7. Also, the axial displacement of the coupled system is small compared to the uncoupled system according to Fig. 8. The range of the changes between maximum and minimum values are almost the same for each time in the uncoupled system. However, the difference between maximum and minimum value decreases versus time in the coupled system.

A half-space subjected to cyclic thermal shock

The previous example, subjected to the cyclic thermal shock, was considered as third example as shown in Fig. 9.

(23)T(t)=\mathrm{5}t{e}^{-\mathrm{2}t}\sin (\mathrm{2}\pi t).

The dimensionless temperature at the center of the plate under thermal cyclic shock is shown in Fig. 10. The dimensionless axial displacement at the center of the plate under thermal cyclic shock was shown in Fig. 11. As it can be seen that the displacement close to the left edge was increased to its maximum value at t=2 and the right side of the plate has zero displacement. The axial shock displacement also progressed in problem domain in x direction with progression in time. The fluctuation of the left side of the wave after it crosses over that time is an important observation made from the comparison of the figures. In other words, the left side of the plate oscillated afterwards for a large time. Also, it can be seen that the maximum value of axial displacement (displacement in positive direction) progressed in time domain. But the minimum value of displacement (displacement in negative direction) was decreased. The range of the oscillation is approximately the same for all simulation times.

Conclusions

Isogeometric analysis was proposed for solution of uncoupled and coupled thermo-elastodynamic problems under constant step function and cyclic thermal shocks. The dynamical system was discretized and solved and the Generalized Newmark method was used in time integration procedure. The performance of the proposed method was verified by solving and comparison of the results of a benchmark example with available exact analytical solutions. Also, examples for step function thermal shock and cyclic thermal shock were solved and the results for the uncoupled and coupled equations were compared. The results showed the fidelity and efficiency of the proposed method in analysis of thermomechanical problems. The conclusions can be categorized as follows:

1) The coupled system showed smaller temperature change and axial displacement with respect to the uncoupled system. This is because of the damping effect considered on the coupled system. Temperature and the axial displacement of the uncoupled system were higher than those for the coupled system.

2) The shapes of the propagated heat wave for the coupled and uncoupled systems acquire different patterns.

References

[1] Eringen A C. Mechanics of Continua (2nd ed). Huntington, New York: Robert E. Krieger Publishing Co., 1980, 1–606

[2] Danilovskaya V I. Thermal stresses in an elastic half-space due to a sudden heating of its boundary. Prikl Mat Mekh, 1950, 14: 316–318 (in Russian)

[3] Sternberg E, Chakravorty J G. On inertia effects in a transient thermoelastic problem. J Appl Mech, 1959, 26: 503–509

[4] Boley B A, Tolins I S. Transient coupled thermoelastic boundary value problems in the half-space. Journal of Applied Mechanics, 1962, 29(4): 637–646

[5] Nickell R E, Sackman J J. Approximate solutions in linear, coupled thermoelasticity. J Appl Mech, 1968, 35(2): 255–266

[6] Suh I G, Tosaka N. Application of the boundary element method to 3-D linear coupled thermoelasticity problems. Theor Appl Mech, 1989, 38: 169–175

[7] Ishiguro S,Tanaka M.Analysis of two-D transient thermoelasticity by the time stepping BEM (in consideration of coupling effect). Trans Jpn Soc Mech Engng, 1993, 59(A): 2593–2598

[8] Dargush G F, Banerjee P K. A new boundary element method for three-dimensional coupled problems of conduction and thermoelasticity. J Appl Mech, 1991, 58(1): 28–36

[9] Tehrani P H, Eslami M R. Two-dimensional time harmonic dynamic coupled thermoelasticity analysis by boundary element method formulation. Engineering Analysis with Boundary Elements, 1998, 22: 245–250

[10] Tehrani P H, Eslami M R. Boundary element analysis of coupled thermoelastic problem with relaxation times in a finite domain. AIAA Journal, 2000, 38(3): 534–541

[11] KC A, Kim J H. Interaction integrals for thermal fracture of functionally graded materials. Engineering Fracture Mechanics, 2008, 75(8): 2542–2565

[12] Hughes T J R, Cottrell J A, Bazilevs Y. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. Computer Methods in Applied Mechanics and Engineering, 2005, 194(39-41): 4135–4195

[13] Hughes T J R, Reali A, Sangalli G. Efficient quadrature for NURBS-based isogeometric analysis. Computer Methods in Applied Mechanics and Engineering, 2010, 199(5–8): 301–313

[14] Nguyen-Thanh N, Zhou K, Zhuang X, Areias P, Nguyen-Xuan H, Bazilevs Y, Rabczuk T. Isogeometric analysis of large-deformation thin shells using RHT-splines for multiple-patch coupling. Computer Methods in Applied Mechanics and Engineering, 2017, 316: 1157–1178

[15] Nguyen-Thanh N, Muthu J, Zhuang X, Rabczuk T. An adaptive three-dimensional RHT-splines formulation in linear elasto-statics and elasto-dynamics. Computational Mechanics, 2014, 53(2): 369–385

[16] Nguyen-Thanh N, Kiendl J, Nguyen-Xuan H, Wüchner R, Bletzinger K U, Bazilevs Y, Rabczuk T. Rotation free isogeometric thin shell analysis using PHT-splines. Computer Methods in Applied Mechanics and Engineering, 2011, 200(47‒48): 3410–3424

[17] Nguyen-Thanh N, Nguyen-Xuan H, Bordas S P A, Rabczuk T. Isogeometric analysis using polynomial splines over hierarchical T-meshes for two dimensional elastic solids. Computer Methods in Applied Mechanics and Engineering, 2011, 200(21‒22): 1892–1908

[18] Ghasemi H, Park H S, Rabczuk T. A level-set based IGA formulation for topology optimization of flexoelectric materials. Computer Methods in Applied Mechanics and Engineering, 2017, 313: 239–258

[19] Ghorashi S Sh, Valizadeh N, Mohammadi S, Rabczuk T. T-spline based XIGA for fracture analysis of orthotropic media. Computers & Structures, 2015, 147: 138–146

[20] Nguyen-Thanh N, Valizadeh N, Nguyen M N, Nguyen-Xuan H, Zhuang X, Areias P, Zi G, Bazilevs Y, De Lorenzis L, Rabczuk T. An extended isogeometric thin shell analysis based on Kirchhoff–Love theory. Computer Methods in Applied Mechanics and Engineering, 2015, 284: 265–291

[21] Nguyen B H, Tran H D, Anitescu C, Zhuang X, Rabczuk T. An isogeometric symmetric Galerkin boundary element method for two-dimensional crack problems. Computer Methods in Applied Mechanics and Engineering, 2016, 306: 252–275

[22] Nguyen V P, Anitescu C, Bordas S P A, Rabczuk T. Isogeometric analysis: an overview and computer implementation aspects. Mathematics and Computers in Simulation, 2015, 117: 89–116

[23] Bazilevs Y, Beirao Da Veiga L, Cottrell J, Hughes T J R, Sangalli G. Isogeometric analysis: approximation, stability and error estimates for h-refined meshes. Mathematical Models and Methods in Applied Sciences, 2006a, 16(7): 1031–1090

[24] Bazilevs Y, Calo V M, Zhang Y, Hughes T J R. Isogeometric fluid structure interaction analysis with applications to arterial blood flow. Computational Mechanics, 2006b, 38(4–5): 310–322

[25] De Luycker E, Benson D J, Belytschko T, Bazilevs Y, Hsu M C. X-FEM in isogeometric analysis for linear fracture mechanics. International Journal for Numerical Methods in Engineering, 2011, 87(6): 541–565

[26] Hosseini-Tehran P Hi M R, Eslami . BEM analysis of thermal and mechanical shock in a two-dimensional finite domain considering coupled thermoelasticity. Engineering Analysis with Boundary Elements, 2000, 24(3): 249–257

[27] Zienkiewicz O C, Chan A H C, Pastor M, Schrefler B A, Shiomi T. Computational Geomechanics with Special Reference to Earthquake Engineering. New York: John Wiley, 1999

[28] Bathe K J. Finite Element Procedures. Englewood Cliffs, NJ: Prentice Hall, 1996

[29] Sternberg E, Chakravorty J G. On inertia effects in a transient thermoelastic problem. J Appl Mech, 1959, 26: 503–509

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag GmbH Germany, part of Springer Nature