Numerical Solution to a One-Dimensional , Nonlinear Problem of Thermoelasticity with Volume Force and Heat Supply in a Slab

A numerical solution is presented for a one-dimensional, nonlinear boundary-value problem of thermoelasticity with variable volume force and heat supply in a slab. One surface of the body is subjected to a given periodic displacement and Robin thermal condition, while the other is kept fixed and at zero temperature. Other conditions may be equally treated as well. The volume force and bulk heating simulate the effect of a beam of hot particles infiltrating the medium. The present study is a continuation of previous work by the same authors for the halfspace [1]. The presented Figures display the process of propagation and reflection of the coupled nonlinear thermoelastic waves in the slab. They also show the effects of volume force and heat supply on the distributions of the mechanical displacements and temperature inside the medium. The propagation of beats provides evidence for sufficiently large time values.


Introduction
he propagation of nonlinear waves in thermoelastic solids is one of the main topics of Continuum Mechanics.Mathematical models describing this phenomenon have been treated by many authors [2,3].In classical thermodynamics, the basic equations of thermoelasticity yield systems of nonlinear partial differential equations of mixed type for which few exact solutions exist.Of theoretical importance are investigations dedicated to the existence, uniqueness and stability of the solutions for such systems as in [4][5][6][7][8][9][10][11][12][13][14][15][16][17].The formation, the development of discontinuities and blow-up of the solutions are treated in [8] and [18][19][20][21][22]. Nonlinear thermoelasticity with finite velocities of propagation of thermal disturbances was considered in [23], and in relation to an electric field in [24].Problems including moving boundaries and multiphase systems with interfaces related to melting, solidification and evaporation processes were treated in [25][26][27].The methods of solution of the various systems of equations of nonlinear thermoelasticity are numerous.They are mainly based on finite elements or finite differences [6,12,20,[28][29][30][31][32][33][34][35][36][37][38].In [30] the author uses an uncoupled numerical scheme to investigate wave propagation driven by initial conditions only.It relies on a finite element technique for the spatial variable, in combination with an uncoupled difference scheme for the time variable.The method of finite volumes is used in [27,39].An extended finite element method to treat crack propagation, stress concentration or flows with interfaces is used in [40].In [32], the authors treat T a one-dimensional problem for a half-space with prescribed harmonic displacement at the boundary.A Poincaré expansion in a small parameter was used to obtain the near-field solution, while the multiple-scale technique took care of the far-field solution.In both methods, the authors obtained a particular solution for the thermoelastic problem, a solution that does not satisfy any thermal boundary conditions at the bounding surface.Other techniques for tackling nonlinear initial/boundary-value problems of continuum mechanics are also available in the literature.Examples are the Method of Finite Volumes, which relies on the possibility of writing the basic field equations in the form of conservation laws [39], and the Differential Quadrature Method introduced by Richard Bellman and collaborators [41][42][43], which is an efficient computational tool for finding solutions of nonlinear initial/boundary-value problems and was successfully applied for solving various engineering problems.This last method makes it possible to avoid difficulties that may arise from quasi-linearization.
In the present work, we use a restriction to thermoelasticity of a fully nonlinear model introduced in [44] to solve a concrete problem of nonlinear wave propagation in a slab, under volume force and bulk heating.One face of the slab is under periodic displacement and Robin thermal condition, while the other one is fixed and kept at zero temperature.This is a continuation of previous work by the same authors for the half-space [1].The main differences here reside in the fact that the far boundary is now fixed, and the time values are allowed to become large enough so as to capture the reflected waves from both surfaces of the slab.Other boundary conditions are also possible without any further complications.In this model, the governing field equations, constitutive relations and boundary conditions are derived in material form within the frame of rigorous thermodynamics.For details concerning the material form of the equations of Continuum Mechanics, we refer to [45].An illustrative application of the present model on onedimensional wave propagation was treated in [46] using a multiple scale technique.Another application was considered in [47] for a one-dimensional, nonlinear thermoelastic wave propagation in a half-space using a finite difference scheme.Amirkhanov [39] investigated a similar problem for a slab using finite volumes, but with no bulk force and for a different boundary condition and for a different heat supply.A problem for a half-space was treated in [1].It was shown that the proposed numerical method correctly represents the phenomenon of thermoelastic wave propagation.The case with three displacement components was investigated in [48].The present model of thermoelasticity relies on the hypothesis that all material coefficients of the medium are functions of strain.This is equivalent to using a free energy function which can be expanded in Taylor's series in its arguments: strain and temperature.Such an approach provided a consistent way to derive the governing system of equations for thermoelasticity for such media, uniformly approximated up to the desired degree of accuracy in terms of the small parameter of strain.Details may be found in [44].
The motivation for introducing body force and heat supply is as follows: Let a beam of heated particles be directed towards the slab.A fraction of those particles penetrates the medium, as a result of which momentum and heat supply are transmitted to the particles of the medium.Thus a field of volume force is created and a bulk heat source is established inside the medium.The authors are not aware of any available experimental evidence that can be used to infer a mathematical form of these fields.For the present purposes, this phenomenon has been modeled using sufficiently smooth functions that decrease exponentially with time and with depth for both the volume force and the heat supply.Other forms could also be used for a concrete form of the heat supply function [39].The left bounding surface of the slab is subject to a given displacement and to a heat radiation condition, while the right boundary is kept fixed and at zero temperature.These rest conditions allow the use of the same numerical technique as in [1].
The study of systems of nonlinear parabolic, or mixed parabolic-hyperbolic equations presents many technical difficulties and may be carried out using different numerical techniques, as is apparent from the given list of references.The merits and disadvantages of each method may be found in specialized textbooks or monographs.In view of the simplicity of the domain geometry, we have opted for a finite-difference scheme which has been tested for effectiveness and reported in previous publications [47,[49][50][51][52].In each case, the numerical results have rendered correctly all the expected characteristics of the solution.
The present system of equations was investigated in [1,47].In the second of these references, it was shown that one of the characteristic curves of this system clearly expresses parabolicity, meaning that a temperature field installs instantaneously in the half-space.The other characteristic corresponds to the propagation of a coupled thermoelastic wave.In dimensionless variables, the velocity of propagation of this latter wave is close to unity, being slightly affected by both strain and temperature.This is clear from the form of the equation of motion.Thus there is only one time scale to be considered, as distinct to other studies where two time scales are necessary for the description of the model [53].The model does not involve any phase transitions as in other publications that deal with the problems of melting, solidification and evaporation [26,27] and the extension to include such effects requires more elaboration on the method used and the computational techniques.Moreover, the calculations have avoided the regions where shock formation occurs because the method does not allow consideration of allow one to effectively consider this phenomenon, especially in the presence of reflected waves.Attention has been mainly focused on the process of reflection of waves at the boundaries.Beats are shown to develop due to the interaction of different wave components with slightly different frequencies.

Formulation of the Problem
Let L be the thickness of the slab.The basic field equations and boundary conditions for the in-depth mechanical displacement and temperature as given in [1,47] read: ) (( 1) )) = ( , ), 0 < < , 2 under the initial conditions ( ,0) = ( ,0) = 0, ( , ) = 0, ( , The symbols u and  denote, respectively, the dimensionless elastic displacement and the dimensionless temperature, x and t are the spatial (depth) coordinate and the time.The suffixes suffices denote differentiation.The constants involved in these equations have obvious physical significance as expressing the dependence of the different material coefficients on strain [32].They will be assumed to have the following orders of magnitude: .


The above governing equations were derived in [32] on the basis of rigorous thermodynamics and in material configuration.Therefore, the boundary conditions are set on a well-defined boundary and do not involve any approximations.
The dimensionless heat flow vector component , Q specific entropy , s and stress component  are determined from the constitutive relations by the expressions: = (1 ) , A finite difference method is used to find a numerical solution of the previous coupled equations under the prescribed initial and boundary conditions.The numerical method and the results are briefly discussed in the following sections.More details may be found in [1,47].

Finite Difference Scheme
A combined approach of quasilinearization and the finite difference iterative method is used to solve the problem.This method has been tested for accuracy and efficiency in previous work [49][50][51][52].Details of the method may be found in these references, and also in [1].
Let () n U and () n  denote the numerical values of u and  at the n -th iteration and let (0) U and (0)  be the initial guess.We consider the Picard approximation [54] for equations ( 1) and ( 2) under which the functions in the sequence () n u and () n  satisfy the boundary conditions specified for u and .
 There is linear convergence of the sequence () n u and () n  to the solution of the original nonlinear problem.For the computational work, consider a finite interval on the x -axis.
Substituting this into (1) and ( 2) and neglecting the truncation error, the resulting equations take the form , The local truncation error of schemes ( 10) and ( 11) is of the order

Numerical Method
The difference scheme (10) and (11) as presented in [47] is a three-level iterative scheme.It may be written in the form , , , , , , , , , q q q q q q p p p p and 4 p are defined above.For = 0, s one has ) From initial conditions (3), ( 4) Therefore, ( 14) and ( 15) finally assume the form  14), ( 15) and ( 20), ( 21) form two tridiagonal linear systems of equations that are used, together with conditions ( 5) and ( 6), to find the values of , rs U and , , rs  = 0,1, 2, , rN and = 1, 2, s .The numerical method for solving schemes ( 14), ( 15) and ( 20), ( 21) subject to the given initial and boundary conditions is as follows: (I) For the first time level 1.Put for all r .
In order to solve (20), (21) one uses Thomas' algorithm [55] to avoid round-off error growth in machine computation, as follows: 1.To solve (20), take the solution of (20) in the form Substitute the values of 1 r U  from ( 22) into (20) and compare the coefficients of the resulting equations with (22).This gives the relations: Taking =0 r in (22) and using the boundary condition where 0 u is given.The scalars r A and r B are computed from 23 in a forward sweep-manner subject to the initial values given by ( 25) and initial conditions (3), (4).Using the stored values of  14) can be solved in a similar manner.
2. To solve (20), take the solution in the form then substitute the values of 1 r   from ( 26) into (21) and compare the coefficients of the resulting equations with (26).This gives the relations: In this case, one uses the boundary condition (1 ) = .
Equation.( 21) at =0 r is Substitute the value of 1   from ( 29) into (30) to get 2 =0 10 10 which can be written in the form Taking =0 r in (26) and comparing with (32), we get where 1 U is obtained from (1), =0 1| r G can be determined from (31), knowing the initial condition.Similarly as in (1), C r and D r are computed from (27) in a forward sweep manner ( = 1,..., 1) rN  subject to the initial values given by (33).

Stability
The study of stability of the difference scheme used for the system of equations introduced earlier relies on the investigation of the eigenvalues of the amplification matrix A defined by A Thus, the stability is the same as for the associated homogeneous difference scheme [56].This has been investigated in [47].According to this reference, the homogeneous schemes associated with the nonhomogeneous difference schemes (10)(11) are unconditionally stable and the local truncation error is of the second-order in both temporal and spatial dimensions.The numerical results show good stability up to sufficiently large values of time.However, the present scheme does not capture the formation of discontinuities which must occur at the (dimensionless) breaking distance

Numerical Results and Discussion
The numerical calculations for the mechanical displacement and temperature distributions in the slab were carried out for the following values of the slab thickness L and the amplitude 0 u of the prescribed harmonic displacement at the boundary: 0 = 600, = 0.008.Lu The volume force and heat supply functions are both taken as smooth functions decreasing exponentially with x and t in the form: ) This is different from the expression for the heat supply function used in [39].
The Figures show the effect of three factors: the boundary displacement, the volume force and the heat supply.Each Figure has four components, labeled from (a) to (d): Label (a) corresponds to the case when the motion is caused only by the boundary displacement, i.e. when = = 0 fg .Labels (b) and (c) correspond respectively to the cases when 0, = 0 fg  and = 0, 0. fg  Label (d) corresponds to the case when 0, 0. fg  Relatively large values of time were needed to evidence the reflection of the waves from both boundaries.One notices the remarkable stability of the solution on these large time intervals.
Figures 1-2 show in perspective the distributions of the mechanical displacement and of the temperature as functions of the distance, for the different time values.The reflection of the wave is obvious and is accompanied for both the displacement and the temperature by a reverse of the amplitude.This is due to the nature of the imposed boundary conditions at the far boundary of the slab.For the chosen values of the different physical parameters, it is seen from Figure 1 that the mechanical displacement produced only by the heat supply is always negative due to dilatation, and is much smaller than that driven only by the volume force, hence this force shadows the effect of a heat supply.Thus the Figures corresponding to labels Figures 3-4 represent the distributions of displacement and temperature as functions of time at the location = 400 x in the slab.One clearly sees on these Figures the instants of time when the waves reflected from the left and from the right boundaries reach the considered location in the slab.One can notice from Figure 3 that the amplitude of the temperature decays with time, a fact that is compatible with the chosen function of heat supply.
Figures 5-6 show the distributions of displacement and temperature in the slab at two time levels = 2100 t and = 2200.t Here, the wave front is moving to the left after the second reflection on the right boundary and before reaching the left one.One can clearly see a phenomenon of beats ahead of the wave.The region occupied by the beats disappears gradually as the wave front moves further.
Figures 7-8 show the effect of variation of the thermoelastic coefficient ( 1  ) on wave propagation, as viewed by an observer fixed inside the slab at the location = 400.
x As this coefficient gets larger, one clearly sees on Figure 7-a the shift of the wave front conditioned by the small increase of the velocity of propagation of the thermoelastic wave.
We have also treated the case with only one period of the driving mechanical displacement at the left boundary.The mechanical boundary condition in this case is: 0 (1 cos ), 0 < 2 ; (0, ) = 0, > 2 .
20) by the backward sweep method using the given boundary conditions and the initial conditions.3.Calculate the values of from (21) by the backward sweep method using the calculated values of r U (calculated from (2)), the boundary conditions and the initial conditions.4. Use the calculated values of r U (

U
is obtained by a backward sweep subject to the known boundary condition N U .System (


is obtained by a backward sweep subject to the known boundary condition = 0.N  System (15) can be solved in a similar manner.


denotes the vector of unknowns at level ,. rsThe non-homogeneous terms are obviously included in the matrix term , B hence they do not affect the eigenvalues of .
(b) and (d) (the two bottom Figures) are almost identical.One cannot draw the same conclusion about the distribution of temperature in the slab.In fact, Figures 2-c and 2-d show marked differences near the left boundary where the heat supply is concentrated.A small thermal disturbance associated with the thermoelastic wave is seen on Figure 2-c progressing in the positive direction and then reflecting from the far boundary.Figures 2-c and 2-d clearly show the mixed parabolic-hyperbolic character of this heat wave.
are shown on Figures 9-14.In particular, Figure10-cshows the thermoelastic wave propagating on the background of the heat supply, while show the propagation of beats.

Figure 3 .
Figure 3. Displacement U as a function of time at x = 400.

Figure 4 .
Figure 4. Temperature as a function of time at x = 400.

Figure 5 .
Figure 5. Displacement U for large times.

Figure 6 .
Figure 6.Temperature for large times.

Figure 7 .
Figure 7. Displacement U as a function of time at x = 400 for two values of .

Figure 8 .
Figure 8. Temperature as a function of time at x = 400 for two values of .

Figure 9 .
Figure 9. Displacement U for one pulse.

Figure 10 .
Figure 10.Temperature for one pulse.

Figure 11 .
Figure 11.Displacement U as a function of time at x = 400 for one pulse.

Figure 12 .
Figure 12.Temperature as a function of time at x = 400 for one pulse.

Figure 13 .
Figure 13.Displacement U for one pulse for large times.

Figure 14 .
Figure 14.Temperature for one pulse for large times.
The domain in the xt  plane is discretized .