SQU Journal for Science, 2015, 20(2), 60-77 ©2015 Sultan Qaboos University 60 Numerical Solution to a One-Dimensional, Nonlinear Problem of Thermoelasticity with Volume Force and Heat Supply in a Slab Wael M. Mohamed, Ahmed F. Ghaleb*, Enaam K. Rawy, Hassan A.Z. Hassan and Adel A. Mosharafa Department of Mathematics, Faculty of Science, Cairo University, Giza 12613, Egypt. *Email: afGhaleb@gmail.com. ABSTRACT: 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 half- space [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. KEYWORDS: Finite difference method; Heat supply; Nonlinear thermoelasticity; Nonlinear wave propagation; Volume force. الحل العذدي لمسألت غيز خطيت أحاديت البُعذ في المزونت الحزاريت بقوة حجميت وتغذيت حزاريت في شزيحت مسطحت فت وائل محمود وأحمذ فؤاد غالب وإنعام خليفت راوي وحسن أحمذ سكي حسن وعادل عطيت مشزَّ ًجٌد لٌة حجمْت مخغْزة ًحغذّت حزارّت معنمذَ حال عذدّا ٌمسأٌت شزًط حذّت أحادّت اٌبُعذ فِ اٌمزًنت اٌحزارّت اٌالخطْت فِ شزّحت النيائْت، ملخص: ىنان ذٌهو ًجٌد شزط رًبن اٌحزارُ، بْنما ُّثبَّج اٌسطح آخز ححج درجت حزارة صفزّت. معاٌشزّحت إساحت دًرّت ِ. ُّعطَ أحذ سطح فِ اٌشزّحت حُحاوِ اٌشزًط اٌمسخخذ مت حأثْز شعاع من اٌجسْماث اٌساخنت اٌخِ حسمظ عٍَ سطح اٌشزّحت ًحخخزليا ٌمسافت ما.ُّعخبز ًإمىانْت ٌفزض شزًط أخزٍ. ار اٌمٌجاث اٌمزنت اٌحزارّت ، ُدرسج فْو ىذه اٌمسأٌت فِ نصف اٌفزاغ. حُبْن األشىاي اٌمعزًضت انخش [1]اٌبحث اٌحاٌِ امخذادا ٌبحث سابك ٌنفس اٌمؤٌفْن اٌشزّحت. وما حبْن أّضا حأثْز اٌمٌة اٌحجمْت ًاٌخغذّت اٌحزارّت عٍَ حٌسّع وً من اإلساحت اٌمْىانْىْت ًاٌحزارة ِاٌالخطْت اٌمخشاًجت ًانعىاسيا من سطح .داخً اٌٌسظ. ُّبْن اٌحً اٌممذََّ ًجٌد ظاىزة اٌضزباث فِ األسمنت اٌىبْزة نسبْا اٌمزًنت اٌحزارّت اٌالخطْت، االنخشار اٌمٌجِ اٌالخطِ، اٌمٌة اٌحجمْت، اٌمنبع اٌحزارُ، طزّمت اٌفزًق اٌمنخيْت. :مفتاحيتكلماث 1. 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-17]. The formation, the development of discontinuities and blow-up of the solutions are treated in [8] and[18-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-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- 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 NUMERICAL SOLUTION TO A ONE-DIMENSIONAL, NONLINEAR PROBLEM 61 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- 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 one- dimensional 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-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. WAEL M. MOHAMED ET AL 62 2. 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: 2 1 2 (1 2 3 ) = ( ) ( , ), tt xx x x x x x u u u u u f x t         (1) 21 ( ) ((1 ) )) = ( , ), 0 < < , 2 x x t x x x au bu u g x t x L      (2) under the initial conditions ( , 0) = ( , 0) = 0, t u x u x (3) ( , 0) = ( , 0) = 0, t x x  (4) 0 x L  and the boundary conditions 0 (0, ) = (1 cos ), (1 (0, )) (0, ) = , 0, x x u t u t u t t H t     (5) and ( , ) = 0, ( , ) = 0.u L t L t (6) 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: 1 3 1 = (1), = (1 10 ), = (10 ),O O to O     3 1 1 2 = (10 ), = (10 ), = (10 ),O a O b O    = (1), = (1).O H O 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 ) , x x Q u   (7) 21 = , 2 x x s au bu   (8) 2 3 1 2 1 = ( ). 2 x x x x u u u u          (9) 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]. 3. 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-52]. Details of the method may be found in these references, and also in [1]. Let ( )nU 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 ( )nu and ( )n  to the solution of the original nonlinear problem. For the computational work, consider a finite interval on the x -axis. The domain in the x t plane is discretized by a grid with step length =x h and time step =t h . The exact values of u and  at the grid point ( , ) ( , )x t rh sk are denoted by ,r sU and ,r s . One takes 2 , 1, 1 1, 1 1, 1 1, 1 1 | = [ ] ( ), 4 x r s r s r s r s r s u u u u u O h h             2 2 2 2 , 1, 1 1, 1 1, 1 1, 12 1 | = [( ) ( ) ] ( ), 8 x r s r s r s r s r s u u u u u O h h             , 1, 1 , 1 1, 1 1, 1 , 12 1 | = [ 2 2 2 xx r s r s r s r s r s r s u u u u u u h              2 1, 1 ] ( ), r s u O h    NUMERICAL SOLUTION TO A ONE-DIMENSIONAL, NONLINEAR PROBLEM 63 2 , , 1 , , 12 1 | = [ 2 ] ( ), tt r s r s r s r s u u u u O k k      2 2 , 1, 1 1, 1 1, 1 1, 1 1 | = [ ] ( ), 4 xt r s r s r s r s r s u u u u u O h k hk              , , 1 , 1 1 | = [ ], 2 r s r s r s       2 , 1, 1 1, 1 1, 1 1, 1 1 | = [ ] ( ), 4 x r s r s r s r s r s O h h                  2 , , 1 , 1 1 | = [ ] ( ), 2 t r s r s r s O k k        2 , 1, 1 , 1 1, 1 1, 1 , 1 1, 12 1 | = [ 2 2 ] ( ), 2 xx r s r s r s r s r s r s r s O h h                        Substituting this into (1) and (2) and neglecting the truncation error, the resulting equations take the form 1 1 1 1 = 4 r r r r qU q U qU W d       2 2 1 3 1 4 5 5 2 4 9 3 4 2 [ ( ) ( ( ) 4 ) ( )] r r r q V d V d d q d q R d q q d         2 2 1 23 1 24 25 5 22 4 21 3 4 22 [ ( ) ( ( ) 4 ) ( )] r r r q Y d Y d d q d q Z d q q d         2 2 1 3 1 4 25 5 2 4 21 4 2 [ ( ) ( ( ) 4 ) ( )] r r r q Y d Y d d q d q R d q d        2 2 2 1 23 1 24 5 5 22 4 9 4 22 [ ( ) ( ( ) 4 ) ( )] 2 , r r r r q V d V d d q d q Z d q d k          (10) 2 2 1 1 1 8 4 2 3 2 4 28 3 28 = ( ) ( ) r r r p p p d p d p d p d p d             5 1 6 1 7 2 5 1 26 1 27 28 [ ] [ ] r r r r r r p R d R d R d p Z d Z d Z d           5 1 6 1 7 2 [ ] r r r p Z d Z d Z d      5 1 26 1 27 28 [ ] 2 , r r r r p R d R d R d k        (11) = 1, 2, , 1, = 0,1, 2, , = 0,1, 2, ,r N s n ( ) ( ) ( ) ( 1)1 , 1 , , 1 , 1 , , , , n n n ns s r s r r s r r s r r s r R S Z              ( ) ( ) ( ) ( 1)1 , 1 , , 1 , 1 , , , , n n n ns s r s r r s r r s r r s r U V U W U Y U U         2 22 2 2 1 2 1 2 3 4 5 ,2 3 2 4 3 = , = 2(1 ), = , = , = , = , , 22 8 8 r s r k kk k k q q q q q q q f hh h h h       1 2 3 4 5 ,2 2 3 = , = 2 1, = 2 1, = , = , = , , 28 2 r s r k b a k p p p p p p p p g hh h h      1 1 1 2 1 1 3 1 = 1 , = , = 2 , r r r r r r r d qV q V qV d V V d V V          4 1 5 1 1 6 1 = 2 , = 2 , = , r r r r r r r d V V d V V V d V V         7 1 8 1 2 1 9 1 1 = , = , = , r r r r r r r d V V d pR p R pR d R R          21 1 1 22 1 1 23 1 = , = , = 2 , r r r r r r d Z Z d Y Y d Y Y         24 1 25 1 1 26 1 = 2 , = 2 , = , r r r r r r r d Y Y d Y Y Y d U U         27 1 28 1 1 = , = , r r r r d U U d U U      and the superscript n denotes the final number of iterations required to obtain an acceptable approximation to the values of ,r sU and ,r s at the grid points on the line =t sk subject to ( 1) ( ) 8 , , | |< 10 , 1 < < 1max n n r s r s r U U r N     (12) and ( 1) ( ) 8 , , | |< 10 , 1 < < 1.max n n r s r s r r N      (13) WAEL M. MOHAMED ET AL 64 The local truncation error of schemes (10) and (11) is of the order 2 2 ( )O h k . 4. 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 1 1 1 = , r r r qU q U qU F     (14) 1 1 1 = , r r r p p p G        (15) = 1, 2, , 1, = 0,1, 2, , = 0,1, 2, ,r N s n where 2 1 2 1 3 1 4 5 5 2 4 9 3 4 2 = 4 [ ( ) ( ( ) 4 ) ( )] r r r r F W d q V d V d d q d q R d q q d           2 2 1 23 1 24 25 5 22 4 21 3 4 22 [ ( ) ( ( ) 4 ) ( )] r r r q Y d Y d d q d q Z d q q d         2 2 1 3 1 4 25 5 2 4 21 4 2 [ ( ) ( ( ) 4 ) ( )] r r r q Y d Y d d q d q R d q d        2 2 2 1 23 1 24 5 5 22 4 9 4 22 [ ( ) ( ( ) 4 ) ( )] 2 , r r r r q V d V d d q d q Z d q d k          (16) 2 2 8 4 2 3 2 4 28 3 28 5 1 6 1 7 2 = ( ) ( ) [ ] r r r G d p d p d p d p d p R d R d R d           5 1 26 1 27 28 5 1 6 1 7 2 [ ] [ ] r r r r r r p Z d Z d Z d p Z d Z d Z d           5 1 26 1 27 28 [ ] 2 , r r r r p R d R d R d k        (17) and 1 2 3 4 5 1 2 3 , , , , , , , , ,q q q q q q p p p p and 4 p are defined above. For = 0,s one has 1,1 1 ,1 1,1 ,0 1, 1 1 , 1 1, 1= 4 ( )r r r r r r rqU q U qU U qU q U qU            2 1,1 1, 1 1,1 ,1 1, 1 , 1[( )( 2 2 )r r r r r rq U U U U U U           1,1 1, 1 ,1 1,1 , 1 1, 1( )(2 2 )]r r r r r rU U U U U U           1,1 ,1 1,1 , 1 1, 1[ 2 2 ]r r r r rU U U U U         2 2 5 1,1 1,1 1, 1 1, 1 [ (( ) ( ) ) r r r r q U U U U           4 ,1 , 1 1,1 1,1 1, 1 1, 14 ( )] [ ]r r r r r rq               3 4 1,1 1,1 1, 1 1, 1[ ( )]r r r rq q U U U U          2 ,0 2 , r k f (18) 1,1 1 ,1 1,1 1, 1 2 , 1 1, 1=r r r r r rp p p p p p                 4 1,1 1,1 1, 1 1, 1( )r r r rp U U U U         2 2 3 1,1 1,1 1, 1 1, 1 [( ) ( ) ] r r r r p U U U U           5 1,1 1, 1 1,1 1, 1 ,1 , 1[( )( )r r r r r rp U U U U            1,1 1, 1 ,1 , 1 1,1 1, 1( )( )r r r r r rU U U U            ,1 , 1 1,1 1, 1 1,1 1, 1( )( )]r r r r r rU U U U            ,02 .rkg (19) From initial conditions (3), (4) , 1 , 1 , 1 , 1 , , | = , | = , 2 2 r s r s r s r s t r s t r s U U U k k         then NUMERICAL SOLUTION TO A ONE-DIMENSIONAL, NONLINEAR PROBLEM 65 ,1 , 1 ,1 , 1 ,0 ,0 | = = 0, | = = 0. 2 2 r r r r t r t r U U U k k       Hence ,1 , 1 ,1 , 1= , = .r r r rU U    Therefore, (14) and (15) finally assume the form 1 1 1 1 = , r r r qU q U qU F     (20) 1 1 1 2 = , r r r p p p G        (21) = 1, 2, , 1, = 0,1, ,r N n where 2 2 1 2 1 23 1 24 25 5 22 4 21 3 4 22 = 2 2 ( ) 2 ( ( 4 ) ( 2 )) 2 , r r r r r F V q Y d Y d d q d q Z d q q d k            1 5 1 26 1 27 28 = 2 ( ) 2 . r r r r G p Z d Z d Z d k        Schemes (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 ,r sU and , ,r s = 0,1, 2, ,r N 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 ( ) ,0 ,1 , , , n r r r r r r V U W V Y U   ( ) ,0 ,1 , , , n r r r r r r R S R Z     for all r . 2. Calculate ( 1) ,1 ( ) n r r U U   from scheme (20) by the backward sweep method using the given boundary conditions and the initial conditions. 3. Calculate the values of ( 1) ,1 ( ) n r r     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 (set ( 1) ( 1) ,1 ,1 , n n r r r r Y U Z      ) to calculate improved values of ( 1) ,1 n r U  from scheme (20). 5. From scheme (21), calculate improved values of r  using the improved values of r U calculated in the previous step. 6. Repeat the iterative procedure (4) and (5) until conditions (12) and (13). (II) Repeat the same at further times, where the values obtained for r U , ( = 1, 2,..., )r N I and r  ( = 0,1...., )r N from step (I) and initial conditions (3), (4) are the initial conditions for schemes (14), (15), while their boundary conditions are given by (5). 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 1 = , = 1, 2, ,1. r r r r U A U B r N n     (22) Substitute the values of 1r U  from (22) into (20) and compare the coefficients of the resulting equations with (22). This gives the relations: 1 1 = , r r E q qA   = , r r q A E (23) 1 1 = , = 1, 2, , 1.r r r qB F B r N E    (24) Taking = 0r in (22) and using the boundary condition 0 = (1 cos ),u u t we get 0 0 1 0 (1 cos ) = ,u t A U B  then 0 0 0 = (1 cos ), = 0,B u t A (25) 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 r A and , r B the solution r U is obtained by a backward sweep subject to the known boundary condition N U . System (14) can be solved in a similar manner. WAEL M. MOHAMED ET AL 66 2. To solve (20), take the solution in the form 1 = , = 1, 2, ,1, r r r r C D r N n       (26) then substitute the values of 1r   from (26) into (21) and compare the coefficients of the resulting equations with (26). This gives the relations: 1 = 2 , r r K C   1 = , r r C K (27) 1 1 = , = 1, 2, , 1.r r r pD G D r N pK    (28) In this case, one uses the boundary condition (1 ) = . x x u H   One has 1 1 1= , = , 2 r r r r x x U U U h h        hence, 1 0 1 1 =0 =0 | = , | = . 2 x r x r U U U h h      Therefore, 1 0 1 1 0 (1 ) = ( )( ) = , 2 x x h U U u H h h            so we may write 1  in the from 2 1 0 1 1 0 2 = . h H h U U         (29) Equation. (21) at = 0r is 1 0 1 =0 2 = 1| . r p p p G       (30) Substitute the value of 1  from (29) into (30) to get 2 =0 1 0 1 0 1 | (1 ) = , 2 r Gh H h U U p       (31) which can be written in the form 1 0 1 =0 1 0 0 2 2 1 0 1 0 ( ) ( 1 | )( ) = . 2 ( ) r h U U G h U U h U U h H p h U U h H                      (32) Taking = 0r in (26) and comparing with (32), we get 2 1 0 1 0 = , h U U h H KK h U U          0 1 = ,C KK (33) =0 0 1 | = , 2 r G D pKK  (34) where 1 U is obtained from (1), =0 1 | r G can be determined from (31), knowing the initial condition. Similarly as in (1), Cr and Dr are computed from (27) in a forward sweep manner ( = 1,..., 1)r N  subject to the initial values given by (33). Using the stored values of r C and , r D the solution r  is obtained by a backward sweep subject to the known boundary condition = 0. N  System (15) can be solved in a similar manner. 5. 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 , 1 , = , r s r s A B     where ,r s denotes the vector of unknowns at level , .r s The non-homogeneous terms are obviously included in the matrix term ,B hence they do not affect the eigenvalues of .A Thus, the stability is the same as for the associated NUMERICAL SOLUTION TO A ONE-DIMENSIONAL, NONLINEAR PROBLEM 67 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 310L  [32, 47]. 6. 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 0u of the prescribed harmonic displacement at the boundary: 0 = 600, = 0.008.L u The volume force and heat supply functions are both taken as smooth functions decreasing exponentially with x and t in the form: 2 ( , ) = ( , ) = (1 ) . t t x f x t g x t t e e e     (35) 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 = = 0f g . Labels (b) and (c) correspond respectively to the cases when 0, = 0f g and = 0, 0.f g  Label (d) corresponds to the case when 0, 0.f g  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 (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. Figures 3-4 represent the distributions of displacement and temperature as functions of time at the location = 400x 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 = 2100t 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 . u t t u t t       (36) The corresponding results are shown on Figures 9-14. In particular, Figure 10-c shows the thermoelastic wave propagating on the background of the heat supply, while Figures 13-14 show the propagation of beats. WAEL M. MOHAMED ET AL 68 Figure 1. Displacement U. Figure 2. Temperature . NUMERICAL SOLUTION TO A ONE-DIMENSIONAL, NONLINEAR PROBLEM 69 Figure 3. Displacement U as a function of time at x = 400. Figure 4. Temperature as a function of time at x = 400. WAEL M. MOHAMED ET AL 70 Figure 5. Displacement U for large times. Figure 6. Temperature for large times. NUMERICAL SOLUTION TO A ONE-DIMENSIONAL, NONLINEAR PROBLEM 71 Figure 7. Displacement U as a function of time at x = 400 for two values of . Figure 8. Temperature as a function of time at x = 400 for two values of . WAEL M. MOHAMED ET AL 72 Figure 9. Displacement U for one pulse. Figure 10. Temperature for one pulse. NUMERICAL SOLUTION TO A ONE-DIMENSIONAL, NONLINEAR PROBLEM 73 Figure 11. Displacement U as a function of time at x = 400 for one pulse. Figure 12. Temperature as a function of time at x = 400 for one pulse. WAEL M. MOHAMED ET AL 74 Figure 13. Displacement U for one pulse for large times. Figure 14. Temperature for one pulse for large times. NUMERICAL SOLUTION TO A ONE-DIMENSIONAL, NONLINEAR PROBLEM 75 7. Conclusions The following conclusions are due: 1. A finite difference iterative method has been used to solve a one-dimensional, nonlinear problem of thermoelasticity involving bulk force and heat supply in a slab. A periodic boundary displacement and Robin thermal condition pertain at the left boundary of the slab, while complete rest prevails at the other boundary. In a previous work, the used numerical method was shown to exhibit unconditional stability. An iteration process at each time level and the use of Thomas' algorithm to avoid round-off error growth and minimize the storage machine computation guarantees the efficiency of the method and good accuracy. The present calculations confirm these facts up to sufficiently large values of time to allow multiple reflections at both boundaries. The (dimensionless) time period, however, could not exceed the dimensionless breaking distance, so that the calculations did not capture the shock wave formation. 2. The method was successfully applied to other boundary displacements, not necessarily periodic. In particular, the case of a boundary displacement consisting of only one time period was also considered. Plots are provided for this case (c.f. Figures 9-14). Other boundary conditions on the right face of the slab are also possible, without any further complications. 3. As concerns the used set of equations, the present work is closely related to [32, 39, 47] and is a continuation of [1]. 4. The wave nature of the solution is clearly represented, while the diffusive character of heat propagation can also be seen on the plots (c.f. Figures 2-c and 10-c). Thus the mixed parabolic-hyperbolic character of the solution is evidenced. Therefore, one may conclude that the proposed numerical method has correctly reproduced the main characteristics of the solution. 5. The successive reflections at the boundaries as viewed by an observer located inside the medium are illustrated on Figures 3,4 and 11,12. 6. The small increase in the velocity of propagation of the thermoelastic wave due to heat is illustrated on some plots by comparing the results for two values of the thermoelastic coefficient (c.f. Figure 7). 7. The study of the wave propagation at sufficiently large time values has put in evidence a phenomenon of propagation of beats. This is a result of the interaction of different wave components with slightly different frequencies. 8. The solutions presented may be used to study the propagation of nonlinear coupled thermoelastic waves at sufficiently small times, before shock formation, so that the results of [1] for the half-space are fully retrieved. In addition, the results may be used to detect the existence of a force field or a heat supply distribution in the medium. References 1. Mahmoud, W., Ghaleb, A.F., Rawy, E.K., Hassan, H.A.Z. and Mosharafa, A. Numerical solution to a nonlinear, one-dimensional problem of thermoelasticity with volume force and heat supply in a half-space. Arch. Appl. Mech., 2014, 84, 9-11, 1501-1515. Arch. Appl. Mech., DOI 2014, 10.1007/s00419-014-0853-y. 2. Maugin, G.A. Physical and mathematical models of nonlinear waves in solids. In: Jeffrey, A. and Engelbrecht, J. (Eds.), Nonlinear Waves in Solids, 109-233, Springer-Verlag, Vienna [CISM Lecture Notes 1993], 1994. 3. Maugin G.A. Nonlinear Waves in Elastic Crystals. Oxford Mathematical Monographs, Oxford University Press, USA, 2000. 4. Slemrod, M. Global existence, uniqueness, and asymptotic stability of classical smooth solutions in one- dimensional nonlinear thermoelasticity. Arch. Rat. Mech. Anal., 1981, 76, 97-133. 5. Chirita, S. Continuous data dependence in the dynamical theory of nonlinear thermoelasticity on unbounded domains. J. Thermal Stresses, 1988, 11, 57-72. 6. Racke, R. Initial boundary-value problems in one-dimensional nonlinear thermoelasticity. Math. Meth. Appl. Sci., 1988, 10, 517-529. 7. Jiang, S. and Racke, R. On some quasilinear hyperbolic-parabolic initial boundary value problems. Math. Meth. Appl. Sci., 1990, 12, 315-319. 8. Racke, R. Blow up in nonlinear three-dimensional thermo-elasticity. Math. Meth. Appl. Sci, 1990, 12, 267-273. 9. Ponce, G. and Racke, R. Global existence of small solutions to the initial value problems for nonlinear thermoelasticity. J. Differential Equations, 1990, 87, 70-83. 10. Racke, R. and Shibata, Y. Global smooth solutions and asymptotic stability in one-dimensional nonlinear thermoelasticity. Arch. Rat. Mech. Anal, 1991, 116, 1-34. 11. Shibata, Y. On one-dimensional nonlinear thermoelasticity. In: Murthy, M.V., Spagnolo, S., (Eds.), Nonlinear Hyperbolic Equations and Field Theory, 178-184, Longman Sci. and Tech., Harlow, Essex, England; John Wiley and Sons, Inc., New York, 1992. 12. Cui, Xia, Yue, J.Y. and Yuan, G.W. Nonlinear scheme with high accuracy for nonlinear coupled parabolic- hyperbolic system. J. Comp. Appl. Math., 2011, 235, 3527-3540. 13. Munoz Rivera, J.E. and Racke, R. Smoothing properties, decay, and global existence of solutions to nonlinear coupled systems of thermoelasticity type. SIAM Journ. Math. Anal., 1995, 26, 1547-1563. 14. Munoz Rivera, J.E. and Barreto, R.K. Existence and exponential decay in nonlinear thermoelasticity. Nonlinear Analysis, 1998, 31, 149-162. WAEL M. MOHAMED ET AL 76 15. Kalpakides, V.K. On the symmetries and similarity solutions of one-dimensional, nonlinear thermo elasticity. Int. J. Eng. Sci., 2001, 39, 1863-1879. 16. Munoz Rivera, J.E. and Qin, Y. Global existence and exponential stability in one-dimensional nonlinear thermo elasticity with thermal memory. Nonlinear Analysis, 2002, 51, 11 - 32. 17. Khalifa, M.E. Existence of almost everywhere solution for nonlinear hyperbolic-parabolic system. Appl. Math. Comp., 2003, 145, 569-577. 18. Dafermos, C.M. and Hsiao, L. Development of singularities in solutions of the equations on nonlinear thermoelasticity. Q. Appl. Math., 1986, 44, 463-474. 19. Racke, R. On the time-asymptotic behaviour of solutions in thermoelasticity. Proc. Roy. Soc. Edinburgh, 1987, 107A, 289-298. 20. Racke, R. Initial boundary-value problems in one-dimensional nonlinear thermoelasticity. Lecture Notes in Mathematics, 1988, 1357, 341-358, Springer, Berlin. 21. Hrusa, J.W. and Messaoudi, S.A. On formation of singularities in one-dimensional nonlinear thermo elasticity. Arch. Rat. Mech. Anal., 1990, 3, 135-151. 22. Babaoglu, C., Erbay, H.A. and Erkip, A. 2012. Global existence and blow-up of solutions for a general class of doubly dispersive nonlocal nonlinear wave equations. arXiv: 1208.1003v1 [math.AP], 5 Aug. 23. Messaoudi, S.A. and Said-Houari, B. Exponential stability in one-dimensional nonlinear thermo elasticity with second sound. Math. Meth. Appl. Sci., 2005, 28(2), 205-232. 24. Ghaleb, A.F. Coupled thermo electro elasticity in extended thermodynamics. In: Hetnarski R.B. (Ed.) Encyclopedia of Thermal Stresses 2014, 1, 767-774, Springer Dordrecht, Heidelberg, New York, London. 25. Didenko, A.N., Ligachev, A.E. and Kurakin, I.B. The interaction of charged particle beams with the metals and alloys surfaces, Energoatomizdat, 1987. 26. Andreaus, U and Dellisola, F. On thermo kinematic analysis of pipe shaping in cast ingots: A numerical simulation via FDM. Int. J. Engng Sci., 1996, 34(12), 1349-1367. 27. Amirkhanov, I.V., Zemlyanaya, E.V., Puzynina, I.V., Puzynina, T.P. and Sarhadov, I. Numerical simulation of the thermoelastic effects in metals irradiated by pulsed ion beam. JCSME, 2002, 2(N1s-2s), 213-224. 28. Jiang, S. Far field behavior of solutions to the equations of nonlinear 1-D thermoelasticity. Appl. Anal., 1990, 36, 25-35. 29. Jiang, S. Numerical solution for the Cauchy problem in nonlinear 1-D thermoelasticity. Computing, 1990, 44, 147-158. 30. Jiang, S. An uncoupled numerical scheme for the equations of nonlinear one-dimensional thermoelasticity. J. Comp. Appl. Math., 1991, 34, 135-144. 31. Jiang, S. On global smooth solutions to the one-dimensional equations of nonlinear inhomogeneous thermoelasticity. Nonlinear Anal.,1993, 20, 1245-1256. 32. Abd-Alla, A.N., Ghaleb, A.F. and Maugin, G.A. Harmonic wave generation in nonlinear thermoelasticity. Int. J. Eng. Sci., 1994, 32, 1103-1116. 33. Sweilam, N.H. Harmonic wave generation in nonlinear thermoelasticity by variational iteration method and Adomian's method. J. Comp. Appl. Math., 2007, 207, 64-72. 34. Sweilam, N.H. and Khader, M.M., Variational iteration method for one-dimensional nonlinear thermoelasticity. Chaos, Solitons and Fractals, 2007, 32, 145-149. 35. Sadighi, A. and Ganji, D.D., A study on one-dimensional nonlinear thermoelasticity by Adomian decomposition method. World Journal of Modelling and Simulation, 2008, 4, 19-25. 36. Cui, Xia and Yue, J.Y. A nonlinear iteration method for solving a two-dimensional nonlinear coupled system of parabolic and hyperbolic equations. J. Comp. Appl. Math., 2010, 234, 343-364. 37. Mohyud-Din, S.T., Yildirim, A. and Gülkanat, Y. Analytical solution of nonlinear thermoelasticity Cauchy problem. World Applied Sciences Journal, 2011, 12, 2184-2188. 38. Akbari Alshati, R. and Khorsand, M. Three-dimensional nonlinear thermoelastic analysis of functionally graded cylindrical shells with piezoelectric layers by differential quadrature method. Acta Mech., 2012, 232, 2565-2590. 39. Amirkhanov, I.V., Sarhadov, I., Ghaleb, A.F. and Sweilam, N.H. Numerical simulation of thermoelastic waves arising in materials under the action of different physical factors. Bulletin of PFUR. Series Mathematics. Information Sciences. Physics, 2013, N2, 64-76. 40. Wei, P., Wang, M.Y. and Xing, X. A study on X-FEM in continuum structural optimization using a level set model. Computer-Aided Desig, 2010, 42, 708-719. 41. Bellman, R.E. A new method for the identification of systems. Math. Biosci., 1969, 5(1-2), 201-204. 42. Bert, C.W. and Malik, M. Differential quadrature method in computational mechanics: A Review. Appl. Mech. Rev.,1996, 49(1), 1-28. 43. Wu, T.Y. and Liu, G.R. The generalized differential quadrature rule for initial-value differential equations. Journal of Sound and Vibration, 2000, 233(2),195-213. 44. Ghaleb, A.F. and Ayad, M.M. Nonlinear waves in thermo-magnetoelasticity. (I) Basic equations. Int. J. Appl. Electromagn. Mat. Mech., 1998, 9(4), 339-357. NUMERICAL SOLUTION TO A ONE-DIMENSIONAL, NONLINEAR PROBLEM 77 45. Eringen, A.C. and Maugin, G.A. Electrodynamics of Continua: Foundations and solid media. Springer-Verlag, New York, 1990. 46. Ghaleb, A.F. and Ayad, M.M. Nonlinear waves in thermo-magnetoelasticity. (II) Wave generation in a perfect electric conductor, Int. J. Appl. Electromagn. Mat. Mech., 1998, 9(4), 359 - 379. 47. Rawy, E.K., Iskandar, L. and Ghaleb, A.F. Numerical solution for a nonlinear, one-dimensional problem of thermoelasticity. Journal of Comp. Appl. Math., 1998, 100, 53-76. 48. Mahmoud, W., Ghaleb, A.F., Rawy, E.K., Hassan, H.A.Z. and Mosharafa, A. Numerical solution to a nonlinear, one-dimensional problem of anisotropic thermoelasticity with volume force and heat supply in a half-space. Interaction of displacements. Arch. Appl. Mech., 2015, 85(4), 433-454. Arch. Appl. Mech.DOI 2014, 10.1007/s00419-014- 0921-3. 49. Jain, P.C. and Iskandar, L. Numerical solutions of the regularized long-wave equation. Comp. Meth. Appl. Mech. Eng., 1979, 20, 195-201. 50. Jain, P.C., Iskandar, L. and Kadalbajoo, M.K. Iterative techniques for non-linear boundary control problems. Proc. Int. Conf. on Optimization in Statistics, Academic Press, New York. 1979, 289-300 51. Iskandar, L. New numerical solution of the Korteweg-deVries equation. Appl. Numer. Math., 1989, 5, 215- 221. 52. Elzoheiry, H., Iskandar, L. and Mohamedein, M. Sh.El-Deen. Iterative implicit schemes for the two-and three- dimensional Sine-Gordon equation. J. Comp. Appl. Math., 1991, 34, 161-170. 53. Luongo, A., Paolone, A. and Di Egidio, A. Multiple Timescales Analysis for 1:2 and 1:3 Resonant Hopf Bifurcations. Nonlinear Dynamics, 2003, 34(3), 269-291. 54. Bellman, R.E. and Kalaba, R.E. Quasilinearization and Nonlinear Boundary-Value Problems. American Elsevier, New-York. 1965. 55. Mitchell, A.R. and Griffiths, D.F. The Finite Difference Method in Partial Differential Equations, Wiley, New York, 1990. 56. Thomas, J.W. Numerical Partial Differential Equations: Finite Difference Methods. Texts in Applied Mathematics, Springer, 1995. Received 1 August 2014 Accepted 9 November 2014