A nonlinear constitutive theory for heat conduction in lagrangian description based on integrity

Abstract

Keywords: Nonlinear Heat Conduction; Solid Continua; Lagrangian Description; Generators and Invariants; Entropy Inequality; Integrity; Temperature Gra

Introduction

Fourier heat conduction law [1–3] stating that heat is proportional to the temperature gradient is the result of experimental observation from heat conduction experiments by Fourier with metal tubes and rods. Based on this, when generalized in R3 it states that the heat vector is equal to the product of a material coefficient called thermal conductivity and the temperature gradients. Thermal conductivity can be any desired function of temperature so long as it is continuous and differentiable. For many simple materials and simple applications this constitutive theory for heat vector may be sufficient to describe the physics of heat conduction. First we note that this constitutive theory is an analytical form to experimental observation. If the experimental observations are limited then the scope of the resulting theory will be limited as well. Secondly, if we need a constitutive theory for heat vector for solid media in Lagrangian description, then it is perhaps more fitting to derive such a theory based on laws and principles of continuum mechanics in conjunction with the physics of heat conduction of the solid matter and then use experiments to calibrate the theory in establishing the material coefficients. It is well known [4–6] that if the disturbed continuous matter is in thermodynamic equilibrium, then the laws of thermodynamcis can be used to describe the response of the continua. Conservation of mass, balance of linear momenta, balance of angular momenta, and the first law of thermodynamics utilize stress tensor and heat vector without regard to how these are arrived at, hence these conservation and balance laws contain no mechanism for deriving constitutive theories. This leaves us with the second law of thermodynamics, the entropy inequality, that must be explored for possible conditions or mechanisms that may be used in the derivations of the constitutive theories. In this approach we must ensure that other conservation and balance laws are not violated.

This paper was recommended for publication in revised for by Editor Mohammed Sajjad Mayeed

1. Department of Mechanical Engineering, University of Kansas, Lawrence, Kansas, USA

2 Department of Mechanical Engineering, Texas A & M University, College Station, Texas, USA

3. Department of Mechanical Engineering, Yildiz Technical University, Istanbul, TURKEY

*E-mail address: kssurana@ku.edu Manuscript Received 19 Septemper 2016, Accepted 21 November 2016

Journal of Thermal Engineering, Vol. 3, No. 6, Special Issue 6, pp. 1615-1631, December, 2017

Scope of Work In the work presented here we consider two possible approaches for deriving constitutive theory for heat vector. The first approach is strictly based on the inequality (or the condition) resulting from the second law of thermodynamics. In the second approach we use the conditions resulting from the entropy inequality to determine conjugate pairs, which establish the argument tensors of the heat vector. Now the theory of generators and invariants can be utilized to derive a constitutive theory for heat vector. This theory is based on integrity [4–6], hence is complete. Material coefficients in the constitutive theories resulting from the two approaches are discussed. The constitutive theory for heat vector resulting from the second approach is up to a cubic function of the temperature gradients while the constitutive theory from the first approach is a linear function of the temperature gradients. While the linear and the nonlinear constitutive theories for heat vector have also been presented by Surana [4], their usefulness in applications and their limitations have not been explored. One of the specific goals in the present paper is to present numerical studies using the constitutive theories presented here to demonstrate when such theories are necessitated and when the nonlinear theories are meritorious. This paper is the first presentation of this kind to our knowledge on the applications of the nonlinear constitutive theory for q in the solution of initial value problems. Two model problems consisting of initial value problems in one spatial dimension (R1 ) and time are used to present numerical studies to illustrate the important features of the two constitutive theories and to demonstrate significant merits of the constitutive theory derived using theory of generators and invariants. Both constitutive theories for heat vector are incorporated in the energy equation (in Lagrangian description) to construct mathematical model describing evolution of temperature in the two initial value problems considered here. The mathematical model is a nonlinear partial differential equation (PDE) in temperature, space coordinate, and time. Numerical solutions of the nonlinear PDE are obtained using space-time finite element method based on space-time residual functional [7,8] in which the local space-time approximations are in higher order scalar product Hilbert spaces H k,p (Ω̄ext ) that permit higher order global differentiability in space and time. Solutions of the nonlinear algebraic equations resulting from the finite element process are obtained using Newton’s linear method with line search. CONSTITUTIVE THEORIES FOR HEAT VECTOR: LAGRANGIAN DESCRIPTION In this section we present details of the constitutive theories for heat vector q . The conditions resulting from the entropy inequality [4] require that q ·gg ≤ 0 (1) be satisfied by the constitutive theories for q regardless of how they are derived. In (1), g is temperature gradient, i.e. g = ∇ θ , θ being absolute temperature. We can take two approaches to derive constitutive theories for q . In the first approach [6, 9], we strictly use (1) to derive the constitutive theory for q . Such a constitutive theory for q will naturally satisfy the entropy inequality as it is derived using the conditions resulting from it. In the second approach we use g as the argument tensor of q (since q and g are conjugate due to (1)) and then use the theory of generators and invariants [Reiner(1945), Todd(1948), 10–23]. The constitutive theories derived using this approach must satisfy (1) so that the disturbed matter will be in thermodynamic equilibrium during evolution. We present the derivation of the constitutive theories for q using both approaches and present comparisons of the resulting constitutive theories, discuss assumptions, and make some remarks regarding their merits and shortcomings. Constitutive theory for q using entropy inequality This derivation based on (1) is fundamental and can be found in any textbook on continuum mechanics [4–6, 9]. We present details in the following to point out the assumptions used in the derivation, as they play a significant role when comparing this constitutive theory with the theories resulting from the theory of generators and invariants. Following references [6, 9], we begin with (1). Equation (1) implies that q ·gg = β ≤ 0 Using equality, we obtain 1616

Journal of Thermal Engineering, Vol. 3, No. 6, Special Issue 6, pp. 1615-1631, December, 2017

That is, heat flux vanishes in the absence of temperature gradient. Thus, the constitutive theory for q must be a function of g . At this stage, many possibilities exist; the simplest of course is assuming that q is proportional to −gg , i.e. q is a linear function of −gg .

∂qq ∂ 2β = = −kk (θ ) ≤ 0 2 ∂ g ∂gg ∂ qi ∂ 2β = = −ki j (θ ) ≤ 0 ∂ g j ∂ gi ∂gj

From (7), we conclude that [k] is positive-semidefinite and all its eigenvalues are non-negative. Equation (5) is the Fourier heat conduction law in Lagrangian description. The thermal conductivity matrix [k] does not have to be symmetric but is often assumed to be. In general, in this constitutive theory for q , the coefficients of [k] can be functions of temperature θ . This constitutive theory is based on the assumption that q is a linear function of g . While mathematically (5) is justified, we note that for isotropic, homogeneous matter, q cannot exhibit directional dependence, hence (5) must be reduced to q = −k(θ )II ·gg

This constitutive theory requires only one material coefficient, thermal conductivity that can be temperature dependent. Constitutive theories for q using theory of generators and invariants In this approach, the heat vector q , a tensor of rank one, is expressed as a linear combination of the combined generators (only tensors of rank one) of its argument tensors. The material coefficients in the linear combination are assumed to be functions of the combined invariants of the argument tensors and temperature θ . The material coefficients are derived by expanding each coefficient in the linear combination in Taylor series about a known configuration. In this approach it is obvious that the explicit form of the constitutive theory for q depends on the argument tensors of q and the terms retained in the Taylor series expansion of the coefficients in the linear combination. In this derivation, we assume that q = q (gg , θ ) (9) This is justified as q and g are conjugate (from (1)). Tensors q and g are tensors of rank one and θ is a tensor of rank zero. The only combined generator of rank one of the argument tensors g and θ is g , hence based on the theory of generators and invariants, we can write q = −qαgg

The coefficient qα is a function of the combined invariants of g , θ , i.e. {g}T {g} and temperature θ . Let us define qI = {g}T {g} to simplify the details of further derivation. We note that (10) holds in the current e 1617

Journal of Thermal Engineering, Vol. 3, No. 6, Special Issue 6, pp. 1615-1631, December, 2017

configuration in which the deformation is not known. Hence in (10), qα = qα(qI , θ ) is not yet deterministic. To e about a known configuration determine material coefficients from (10), we expand qα(qI , θ ) in Taylor series e q q Ω in I and θ and retain only up to linear terms in I and θ . e e q
∂ qα ∂ α q q α = qα Ω + q I − (qI )Ω + (θ − θΩ ) (11) ∂ I Ω e ∂θ Ω e e Substituting from (11) into (10)

 
∂ qα ∂ qα q I − (qI )Ω + (θ − θΩ ) g q = − qα Ω + q ∂ I Ω e ∂θ Ω e e

We note that qα|Ω , ∂∂ qαI Ω and ∂∂ θα Ω are functions of (qI )Ω and θΩ , whereas qα in (10) is a function of qI e e e and θ in the current configuration. From (12) we can write the following, noting that qI = {g}T {g} e qα qα qα ∂ ∂ ∂ q = −qα|Ω g − q {g}T {g}gg + q ({g}T {g})Ω g − (θ − θΩ )gg (13) ∂ I Ω ∂ I Ω ∂θ Ω e e or   ∂ qα ∂ qα ∂ qα (θ − θΩ )gg (14) q = − qα|Ω + q ({g}T {g})Ω g − q {g}T {g}gg − ∂ I Ω ∂ I Ω ∂θ Ω e e Let ∂ qα k(θΩ , (qI )Ω ) = qα|Ω + q ({g}T {g})Ω ∂ I Ω e e q ∂ α k1 (θΩ , (qI )Ω ) = q ∂ I Ω e e ∂ qα k2 (θΩ , (qI )Ω ) = ∂θ Ω e

then If we neglect the θ term then we can write
q = − k + k1 (gg · g ) g

This is the simplest possible constitutive theory based on the theory of generators and invariants using (9). The only assumption in this theory beyond (9) is the truncation of the Taylor series in (11) beyond linear terms in qI and θ . The constitutive theory for q given by (8) is recoverable from (17) if the second term e side of (17) is neglected and we assume that k = k(θ ), i.e. k only depends on temperature θ . on the right Clearly (8) is a subset of (17). We note that using q from (17) we obtain q · g = −k(gg · g ) − k1 (gg · g )(gg · g )

Since g · g > 0 for g 6= 0 and k, k1 > 0, we obtain the following from (18). q·g <0 That is, the constitutive theory for q in (17) satisfies the condition (1) resulting from entropy inequality.

Journal of Thermal Engineering, Vol. 3, No. 6, Special Issue 6, pp. 1615-1631, December, 2017

COMPLETE MATHEMATICAL MODEL: ENERGY EQUATION When the velocity field and stress field are zero, for isotropic, homogeneous solid matter the energy equation in Lagrangian description constitutes the complete mathematical model. In the absence of sources and sinks we can write the following for the energy equation [4] in R1 (one dimensional case). ∂θ ∂q + = 0 ∀(x,t) ∈ Ωxt = Ωx × Ωt = (0, L) × (0, τ) (20) ∂t ∂x ρ0 is density and cv is specific heat at constant volume, both in reference configuration. θ is absolute temperature. q in (20) is defined by (17). For standard Fourier heat conduction we set k1 = 0, k = k(θ ) in (17). Substituting q from (17) in (20) and assuming k, k1 to be constant, we can derive another form of the energy equation, a single nonlinear PDE in temperature θ , spatial coordinate x, and time t. ρ0 cv

Dimensionless form of the energy equation First, we express all quantities with hat ( ˆ ) in (21), implying that they all have their usual dimensions. ∂ θ̂ ∂ 2 θ̂ ρ̂0 ĉv − k̂ 2 − 3k̂1 ∂ tˆ ∂ x̂

We choose the following reference quantities (subscript 0 or re f ) and define the dimensionless variables (without hat ( ˆ )). x=

Using (23) in (22), the following form of the mathematical model can be obtained. (ρ0 )re f cv0 θ0 t0 or

 2    2 2 θ0 ∂ θ k0 θ0 ∂ 2 θ θ0 ∂θ ∂ θ − 2 k 2 − =0 (3 k̂ ) 1 ∂t ∂ x ∂ x2 L0 ∂ x L02 L02

 2 2  2  θ0 ∂θ ∂ θ ∂ 2θ ∂θ (3k̂1 ) ρcv −k 2 − =0 ∂t ∂x ∂ x ∂ x2 L02 k0

 2   2 2 θ0 ∂θ ∂ 2θ ∂θ ∂ θ −k 2 − (k ) (3k ) =0 1 re f 1 ∂t ∂x ∂ x ∂ x2 L02 k0

Equation (28) is the final form of 1D energy equation with constant material coefficients that constitutes the mathematical model for model problems considered in this paper.

Journal of Thermal Engineering, Vol. 3, No. 6, Special Issue 6, pp. 1615-1631, December, 2017

Space-Time Finite Element Method And Numerical Studies

In this section, two numerical studies are presented using the energy equation (28). Computations are performed using space-time least squares finite element formulation based on space-time residual functional [7, 8] for nonlinear space-time differential operators in which a space-time strip for an increment of time ∆t is discretized using a row of p-version hierarchical space-time higher order global differentiability finite elements. The nonlinear algebraic equations resulting from the space-time integral form are solved using Newton’s linear method with line search. The evolution is computed for an increment of time and when converged, is then time-marched to compute the entire evolution for desired final value of time. The computational processes in this approach are unconditionally stable throughout the evolution due to varationally consistent space-time integral form [7, 8]. A space-time finite element method using a space-time strip for an increment of time with time marching is used to compute the evolution. Let Aθ − f = 0

Ω̄xt (as in Figure 3) in which Ω̄ext is a space-time finite element. Let θh (x,t) be the approximation of θ over Ω̄Txt such that θh (x,t) = θhe (x,t) S

is the residual function over discretization Ω̄Txt of the space-time strip Ω̄xt and E e = Aθhe − f

Let I and I e be the residual functionals over Ω̄Txt and Ω̄ext , then we can write (scalar product of E with itself) I(θh ) = (E, E)Ω̄Txt = ∑ I e = ∑(E e , E e )Ω̄ext e

If I(θh ) is differentiable in its arguments, then its first variation set to zero is a necessary condition for an extremum of I(θh ). δ I(θh ) = 2(E, δ E)Ω̄Txt = 2 ∑(E e , δ E e )Ω̄ext = 2 ∑{ge } = 2{g} = 0

δ 2 I(θh ) ∼ = 2(δ E, δ E)Ω̄Txt = 2 ∑(δ E e , δ E e )Ω̄ext > 0

is a unique extremum principle or sufficient condition, hence the integral form in (34) is variationally consistent implying that this computational process is unconditionally stable. Since {g} in (34) is a nonlinear function of θh as the differential operator A is nonlinear, we use Newton’s linear method with line search to find a solution θh that satisfies the necessary condition (34). Let θh0 be an assumed starting solution (generally initial conditions can be used for this), then the improved solution θh is given by (see reference [7] for details) θh = θh0 + α∆θh 1 {g}θ 0 ∆θh = − [δ 2 I(θh )]−1 θh0 h 2 0 < α ≤ 2 such that I(θh ) ≤ I(θh0 )

If max|gi (θh )| ≤ ∆ where ∆ is a preset tolerance for zero, then θh is the converged solution in the iterative i

process, otherwise we set θh0 to θh and repeat computation of ∆θh , θh , and {g(θh )}. 1620

Journal of Thermal Engineering, Vol. 3, No. 6, Special Issue 6, pp. 1615-1631, December, 2017

Remarks (1) The converged solution θh is calculated for the first space-time strip corresponding to time increment ∆t by using mesh refinement and/or p-level increase. (2) Solution θh at t = ∆t is used to determine ICs for the second space-time strip (∆t ≤ t ≤ 2∆t) and then a converged solution is calculated for the second space-time strip. This procedure is continued for subsequent space-time strips until the desired time t = τ is reached. (3) Accuracy of the solution obtained from space-time finite element method is important to discuss. The choices of the orders of the approximation space in x and t are such that the space-time integrals over the discretization Ω̄Txt of each space-time strip are Riemann. With these choices, when I(θh ) over Ω̄Txt approaches zero (i.e. O(10−15 ) or lower) due to mesh refinement and p-level increase in space and time, we have I(θh ) → 0 =⇒ E → 0 in pointwise sense due to the fact that I(θh ) = (E, E)Ω̄Txt is Riemann. This of course means that Aθh − f = 0 holds in the pointwise sense (within the tolerance of O(10−15 ) or lower for I(θh )). The solutions reported in the paper satisfy the conditions discussed above, confirming the numerical solutions presented in the paper have extremely high accuracy and can in fact be treated as theoretical solutions. The first numerical study consists of 1D heat conduction in a rod of length L (0 ≤ x ≤ L). The left end (x = 0) is insulated and the right end (x = L) is subjected to a change in temperature over time ∆t from θ1 to θ2 (θ1 being initial condition). The right end is held at temperature θ = θ2 for all values of time beyond ∆t. The second model problem is similar to “model problem I” except that both left and right boundaries are insulated and the initial condition consists of a Gaussian distribution centered at x = 0.2 (L = 1.0). We consider evolution of temperature θ using space-time coupled finite element method for an increment of time ∆t. Upon convergence, the solution is time-marched until the entire evolution is obtained. Model problem I: 1D transient heat conduction in a rod Consider a rod of uniform cross-section. The left end of the rod and its surface are insulated and the right end of the bar is subjected to a continuous and differentiable temperature change from θ1 to θ2 (θ2 > θ1 ) over an increment of time ∆t and θ (L,t) = θ2 for t ≥ ∆t. For 0 ≤ t ≤ ∆t, θ is a cubic distribution in time such that θ = θ1 and ∂∂tθ = 0 at t = 0 and θ = θ2 and ∂∂tθ = 0 at t = ∆t. Figure 1 shows a schematic with boundary conditions and Figure 2 shows θ at the right boundary as a function of time for 0 ≤ t ≤ τ, where τ is the final value of time. BC: ∂θ =0 ∂x

   θ1 + θ2 − θ1 (3t 2 ∆t − 2t 3 ) ; t ≤ ∆t ∆t 3 θ (L,t) =   θ ; t ≥ ∆t 2

Figure 1: Model problem I: schematic of 1D heat conduction with boundary conditions

Journal of Thermal Engineering, Vol. 3, No. 6, Special Issue 6, pp. 1615-1631, December, 2017

We consider the following material properties and reference quantities. L0 = L̂ = 0.1 m ; cv0 = ĉv = 903 J/kg-K ;

(ρ0 )re f = ρ̂0 = 2700 kg/m3 which yield the following dimensionless quantities. L=1;

k1 = 37974.68k̂1 We consider purely 1D heat conduction in the rod of Figure 1, i.e. we only consider centerline of the rod. At t = 0 (commencement of evolution) initial temperature θ̂ (x̂, 0) = θ̂1 defines the initial condition. For an increment in time ∆t the space-time strip [0, L] × [0, ∆t] = [0, 1] × [0, 0.1] is discretized using a twenty element uniform space-time finite element discretization with element length of 0.05 in space and 0.1 in time (Figure 3). Initial p-convergence studies with this discretization suggest that p = 11 with k = 2, local approximations of class C1,1 (Ωext ), are sufficient for good accuracy of results. The residual or least squares functional values for the space-time strip remain O(10−16 ) – O(10−14 ), indicating that the PDEs are satisfied accurately. Newton’s linear method used for solving the nolinear algebraic equations converges in less that 7 iterations for all increments of time. In the numerical studies, we consider the following dimensionless values of temperature θ2 : t a nine-node p-version space-time element t = ∆t = 0.1

Figure 3: Model problem I: twenty element uniform discretization of the first space-time strip

which correspond to 360 K, 390 K, 510 K, and 750 K. We choose the following values of the coefficient k̂1 : k̂1 = 0 W/m-K3 (k1 = 0), k̂1 = 0.0001 W/m-K3 (k1 = 3.7975), k̂1 = 0.0005 W/m-K3 (k1 = 18.987). When k̂ = 0, (17) reduces to standard Fourier heat conduction law. The values of k̂1 are chosen so that with progressively increasing temperature gradient, we can clearly observe the influence of the nonlinear terms and hence, deviation from the Fourier heat conduction law. When the temperature gradients are high enough such that their squares and higher powers cannot be neglected compared to them, then perhaps consideration of nonlinearity in temperature gradients is meritorious in the constitutive theory for heat vector q . First, we examine temperature gradients for linear heat conduction. Figure 4(a)–(d) shows plots of temperature gradient versus x for first four increments of time for four pairs of (θ1 , θ2 ): (1, 1.2), (1, 1.3), (1, 1.7), and (1, 2.5) values. ∂∂θx is most pronounced for the first time increment as expected. ∂∂θx increases with progressively increasing value of θ2 . Even at θ2 = 1.2,  2 peak value of ∂∂θx ≃ 0.8 suggests that ∂∂θx ≃ 0.64 is not negligible. At θ2 = 2.5 the peak value of ∂∂θx ≃ 4  2 (Figure 4(d)) gives rise to ∂∂θx ≃ 16 which is certainly not negligible compared to ∂∂θx ≃ 4. All these graphs  2 in Figure 4 show that in the initial stages of the evolution ∂∂θx are quite significant compared to ∂∂θx . Next, we examine ∂∂θx when evolution is computed using nonlinear constitutive thory for heat vector. These studies are considered for two values of k1 (given earlier). Figure 5(a)–(d) shows graphs of ∂∂θx versus x for first four increments of time for the same values of (θ1 , θ2 ) as used in case of linear heat conduction 1622

Journal of Thermal Engineering, Vol. 3, No. 6, Special Issue 6, pp. 1615-1631, December, 2017

Figure 4: Model problem I: temperature gradient ∂∂θx versus x for first four time steps, linear heat conduction studies. Compared to linear heat conduction, ∂∂θx values are reduced due to enhanced heat conduction but are still significant enough to warrant use of nonlinear constitutive theory for q . With progressively increasing θ2 , ∂∂θx continues to increase and its highest value is for θ2 = 2.5. In these studies also we note that the ∂∂θx values are most pronounced for the first increment of time as expected. A significant point to note is that in figure 5(a)–(d) (nonlinear conduction) the values of ∂∂θx remain far below those in case of linear heat conduction due to enhanced heat conduction. As expected, higher value of k1 results in increased heat conduction hence lower values of ∂∂θx as can be seen in Figure 5(a)–(d). These studies confirm that significantly higher temperature gradients can result due to applied temperature boundary condition during  2 and higher initial stages of the evolution. In such situations linear heat conduction that neglects ∂∂θx power terms may not be a good constitutive theory for the heat vector. Figures 6 – 9 show graphs of the evolution of temperature θ versus x for four different values of θ2 with the three choices of k1 . We observe that with progressively increasing values of θ2 , the numerical solutions obtained for nonzero k1 begin to deviate significantly from Fourier heat conduction law. For θ2 = 2.5 (Figure 9), the temperature distribution along the length of the rod differs significantly compared with Fourier heat conduction law. Larger value of k1 obviously results in greater deviation from Fourier heat conduction law. This study demonstrates that when temperature gradients are high, the nonlinear constitutive theory for the heat vector may be a more realistic representation of the physics as opposed to standard Fourier heat conduction law.

Journal of Thermal Engineering, Vol. 3, No. 6, Special Issue 6, pp. 1615-1631, December, 2017

Figure 5: Model problem I: temperature gradient ∂∂θx versus x for first four time steps, nonlinear heat conduction

Figure 6: Model problem I: temperature θ versus x for θ2 = 1.2

Journal of Thermal Engineering, Vol. 3, No. 6, Special Issue 6, pp. 1615-1631, December, 2017

Figure 7: Model problem I: temperature θ versus x for θ2 = 1.3

Figure 8: Model problem I: temperature θ versus x for θ2 = 1.7

Figure 9: Model problem I: temperature θ versus x for θ2 = 2.5

Journal of Thermal Engineering, Vol. 3, No. 6, Special Issue 6, pp. 1615-1631, December, 2017

Model problem II: 1D transient heat conduction with Gaussian distribution as initial condition In this case we consider 1D transient heat conduction with different boundary conditions and initial conditions. Figure 10 shows a schematic of the problem with boundary conditions (both ends are insulated and the surface of the rod is insulated as well). For this model problem, a 20 element uniform discretization is used with element length of 0.05 and ∆t = 0.0005 BC:

Figure 10: Model problem II: schematic of 1D heat conduction with boundary conditions Figure 11 shows initial condition at t = 0, a Gaussian distribution of temperature θ centered at x = 0.2 and with a support of one unit and a peak value of θ = 1 at x = 0.2. The distribution ranges over 0.1 ≤ x ≤ 0.3 in continuous and differentiable manner. Equation (37) describes the Gaussian distribution with mean x0 = 0.2 and standard deviation σ0 = 0.3. θ (x) 1

Figure 11: Model problem II: initial condition θ (x, 0), Gaussian distribution   (x − x0 )2 θ (x, 0) = exp − 2(σ0 )2

Initial p-convergence studies suggest that p = 7 with k = 2, local approximations of class C1,1 (Ω2xt ), are sufficient for accurate evolution. The residual or least squares functional values remain O(10−18 ) – O(10−14 ) indicating that the PDEs are satisfied accurately. Newton’s linear method used for solving the nonlinear algebraic equations converges in less than 7 iterations for all space-time strips. The following values of the coefficient k̂1 are used. k̂1 = 0.0 W/m-K3 ; k̂1 = 10

For this model problem also we first examine the temperature gradient ∂∂θx versus x for linear and nonlinear heat conduction during initial stages of evolution. Figure 12(a)–(d) shows plots of ∂∂θx versus x for t = ∆t, t = 3∆t, t = 5∆t, and t = 10∆t for k̂1 = 0 (linear heat conduction) and k̂1 = 10−6 , 5 × 10−5 , and 10−5 (nonlinear heat conduction). In case of linear heat conduction at t = ∆t peak value of ∂∂θx is close to ten. At t = 5∆t the 1626

Journal of Thermal Engineering, Vol. 3, No. 6, Special Issue 6, pp. 1615-1631, December, 2017

peak value of ∂∂θx is still greater than one. The values of ∂∂θx indicate that

significant compared to ∂∂θx and are worthy of consideration in the constitutive theory for heat conduction. Values of ∂∂θx versus x shown in Figure 12(b)–(d) for the same values of time as in Figure 12(a) for nonlinear heat conduction show lower values of ∂∂θx compared to linear heat conduction but still significant enough so  2 that ∂∂θx is not negligible compared to ∂∂θx . As expected, higher values of k̂1 result in lower values of ∂∂θx due to enhanced heat conduction. These studies confirm that higher temperature gradients in the initial conditions result in higher temperature gradients during initial stage of evolution. In such situations linear  2 heat conduction law that neglects ∂∂θx may be a significant approximation to the actual physics of heat conduction. 25

Figure 12: Model problem II: temperature gradient ∂∂θx versus x

Journal of Thermal Engineering, Vol. 3, No. 6, Special Issue 6, pp. 1615-1631, December, 2017

Graphs of temperature θ versus x for t = ∆t, 3∆t, 5∆t, 10∆t, and 50∆t for k̂1 = 0 (linear heat conduction) and k̂1 = 10−6 , 5 × 10−5 , and 10−5 (nonlinear heat conduction) are shown in Figures 13 – 16. In Figure 13 we observe that as time elapses peak of the Gaussian distribution progressively reduces followed by progressive elongation of the base of the Gaussian distribution. Similar behavior is also observed for nonlinear heat conduction. When the temperature evolution in Figures 14 – 16 is compared with linear heat conduction (Figure 13) we observe that (i) the peak values of θ are lower in case of nonlinear heat conduction, (ii) peak values of θ are progressively reduced with increasing values of k̂1 due to enhanced heat conduction, and (iii) the base of the Gaussian distribution for nonlinear heat conduction is larger compared to linear heat conduction. This is quite obvious from the values of θ at x = 0 in Figures 14 – 16 when these are compared with the corresponding values in Figure 13. IC t=1∆t t=3∆t t=5∆t t=10∆t t=50∆t

Figure 13: Model problem II: temperature θ versus x for k̂1 = 0.0

Figure 14: Model problem II: temperature θ versus x for k̂1 = 10−6

Journal of Thermal Engineering, Vol. 3, No. 6, Special Issue 6, pp. 1615-1631, December, 2017

Figure 15: Model problem II: temperature θ versus x for k̂1 = 5 × 10−6 IC t=1∆t t=3∆t t=5∆t t=10∆t t=50∆t

Figure 16: Model problem II: temperature θ versus x for k̂1 = 10−5 SUMMARY AND CONCLUSIONS Constitutive theories for heat conduction in solid matter are derived in Lagrangian description using the conditions resulting from entropy inequality. In the first derivation the condition resulting from the entropy inequality is used directly. This results in standard Fourier heat conduction law in which the thermal conductivity can be a continuous and differentiable function of temperature. In the second approach the condition resulting from the entropy inequality is used to infer the conjugate pairs q and g , and since g can be derived by simple differentiation of θ , q must be the dependent variable in the constitutive theory. g and θ serve as argument tensors of q . Theory of generators and invariants is used to derive constitutive theory for q that is based on integrity, hence complete. This constitutive theory contains up to third power of the temperature gradients as well as a linear term in temperature gradient and requires material coefficients k and k1 . k is standard thermal conductivity associated with the temperature gradient while k1 is an additional material coefficient associated with the cubic terms in temperature gradient. Both k and k1 can be functions of θ and {g}T {g} (invariant of {g}) in a known configuration Ω . When k1 is set to zero and k = k(θ ), the standard Fourier heat conduction law is recoverable from this nonlinear constitutive theory for heat vector q . Evolutions of temperature for two model problems are presented. Both initial value problems contain high temperature gradients (dimensionless values >> 1) in the initial stages of the evolution, hence their squares are quite significant in comparison with the temperature gradients. Linear heat conduction studies

Journal of Thermal Engineering, Vol. 3, No. 6, Special Issue 6, pp. 1615-1631, December, 2017

show 1 ≤ ∂∂θx ≤ 4.5 for up to 3 to 4 initial time steps. In such cases

From the studies and the constitutive theory we observe that nonlinear heat conduction enhances heat conduction, hence results in significantly lower temperatures, especially during the initial stages of the evolution when ∂∂θx is high. As the evolution proceeds, enchanced heat conduction effectively lowers ∂∂θx (more than linear heat conduction), hence the influence of the nonlinear terms in the constitutive theory for q becomes progressively less significant during later stages of the evolution. This can be observed clearly by comparison of the linear and nonlinear heat conduction studies during the later stages of the evolution. In high temperature physics that invariably is accompanied by isolated high temperature gradients, the nonlinear constitutive theory for q may be of significant benefit in describing more realistic physics of heat conduction.

Acknowledgements

The first and fifth authors are grateful for the support provided by their endowed professorships during the course of this research. The computational infrastructure provided by the Computational Mechanics Laboratory (CML) of the Mechanical Engineering department of the University of Kansas is gratefully acknowledged. The financial support provided to the second author by the Naval Air Warfare Center is greatly appreciated.

References

1. Bennett, C. O.; Myers, J. E. Momentum, Heat, and Mass Transfer ; McGraw-Hill, 1962.
2. Anderson, D. A.; Tannehill, J. C.; Pletcher, R. H. Computational Fluid Mechanics and Heat Transfer ; CRC Press, 1984.
3. Incropera, F. P.; DeWitt, D. P. Introduction to Heat Transfer ; John Wiley & Sons, 1985.
4. Surana, K. S. Advanced Mechanics of Continua; CRC/Taylor & Francis: Boca Raton, 2015.
5. Reddy, J. N. An Introduction to Continuum Mechanics; Cambridge University Press, 2013.
6. Eringen, A. C. Mechanics of Continua; Robert E. Krieger Publishing Co., 1980.
7. Surana, K. S.; Reddy, J. N. The Finite Element Method for Boundary Value Problems: Mathematics and Computations; CRC/Taylor & Francis: Boca Raton, 2017.
8. Surana, K. S.; Reddy, J. N. The Finite Element Method for Initial Value Problems: Mathematics and Computations; CRC/Taylor & Francis: Boca Raton, 2018.
9. Eringen, A. C. Nonlinear Theory of Continuous Media; McGraw-Hill, 1962.
10. Reiner, M. A Mathematical Theory of Dilatancy. American Journal of Mathematics, 1945, 67, 350–362.
11. Todd, J. A. Ternary Quadratic Types. Philosophical Transactions of the Royal Society of London. Series A: Mathematical and Physical Sciences, 1948, 241, 399–456.
12. Rivlin, R. S.; Ericksen, J. L. Stress-Deformation Relations for Isotropic Materials. Journal of Rational Mechanics and Analysis, 1955, 4, 323–425.
13. Rivlin, R. S. Further Remarks on the Stress-Deformation Relations for Isotropic Materials. Journal of Rational Mechanics and Analysis, 1955, 4, 681–702.
14. Wang, C. C. On Representations for Isotropic Functions, Part I. Archive for Rational Mechanics and Analysis, 1969, 33, 249.
15. Wang, C. C. On Representations for Isotropic Functions, Part II. Archive for Rational Mechanics and Analysis, 1969, 33, 268. 1630 Journal of Thermal Engineering, Vol. 3, No. 6, Special Issue 6, pp. 1615-1631, December, 2017
16. Wang, C. C. A New Representation Theorem for Isotropic Functions, Part I and Part II. Archive for Rational Mechanics and Analysis, 1970, 36, 166–223.
17. Wang, C. C. Corrigendum to ‘Representations for Isotropic Functions’. Archive for Rational Mechanics and Analysis, 1971, 43, 392–395.
18. Smith, G. F. On a Fundamental Error in two Papers of C.C. Wang, ‘On Representations for Isotropic Functions, Part I and Part II’. Archive for Rational Mechanics and Analysis, 1970, 36, 161–165.
19. Smith, G. F. On Isotropic Functions of Symmetric Tensors, Skew-Symmetric Tensors and Vectors. International Journal of Engineering Science, 1971, 9, 899–916.
20. Spencer, A. J. M.; Rivlin, R. S. The Theory of Matrix Polynomials and its Application to the Mechanics of Isotropic Continua. Archive for Rational Mechanics and Analysis, 1959, 2, 309–336.
21. Spencer, A. J. M.; Rivlin, R. S. Further Results in the Theory of Matrix Polynomials. Archive for Rational Mechanics and Analysis, 1960, 4, 214–230.
22. Spencer, A. J. M. Theory of Invariants In: Treatise on Continuum Physics, I A. C. Eringen, Ed.; Aca- demic Press, 1971.
23. Boehler, J. P. On Irreducible Representations for Isotropic Scalar Functions. Journal of Applied Math- ematics and Mechanics / Zeitschrift für Angewandte Mathematik und Mechanik, 1977, 57, 323–327.
24. Zheng, Q. S. On the Representations for Isotropic Vector-Valued, Symmetric Tensor-Valued and Skew- Symmetric Tensor-Valued Functions. International Journal of Engineering Science, 1993, 31, 1013– 1024.
25. Zheng, Q. S. On Transversely Isotropic, Orthotropic and Relatively Isotropic Functions of Symmetric Tensors, Skew-Symmetric Tensors, and Vectors. International Journal of Engineering Science, 1993, 31, 1399–1453. 1631

Share and Cite