Abbas Tcharkhtchi
email: abbas.tcharkhtchi@ensam.eu

Sofiane Abdallah-Elhirtsi
email: sofiane.abdallah-elhirtsi@ensam.eu

Kambiz Ebrahimi
email: m.ebrahimi@bradford.ac.uk

Joseph Fitoussi
email: joseph.fitoussi@ensam.eu

Mohammadali Shirinbayan
email: mohammadali.shirinbayan@ensam.eu

Sedigheh Farzaneh
email: sedigheh.farzaneh@ensam.eu

Some New Concepts of Shape Memory Effect of Polymers

Keywords: polymers, shape memory effect, properties memory effect, polyurethane

In this study some new concepts regarding certain aspects related to shape memory polymers are presented. A blend of polylactic acid (PLA) (80%) and polybutylene succinate (PBS) (20%) was prepared first by extrusion, then by injection molding to obtain the samples. Tensile, stress-relaxation and recovery tests were performed on these samples at 70 °C. The results indicated that the blend can only regain 24% of its initial shape. It was shown that, this partial shape memory effect could be improved by successive cycles of shape memory tests. After a fourth cycle, the blend is able to regain 82% of its shape. These original results indicated that a polymer without (or with partial) shape memory effect may be transformed into a shape memory polymer without any chemical modification. In this work, we have also shown the relationship between shape memory and property memory effect. Mono and multi-frequency DMA (dynamic mechanical analyzer) tests on virgin and 100% recovered samples of polyurethane (PU) revealed that the polymer at the end of the shape memory tests regains 100% of its initial form without regaining some of its physical properties like glass transition temperature, tensile modulus, heat expansion coefficient and free volume fraction. Shape memory (with and without stress-relaxation) tests were performed on the samples in order to show the role of residual stresses during recovery tests. On the basis of the results we have tried to show the origin of the driving force responsible for shape memory effect.

Introduction

The capability of the shape memory polymers (SMPs) to remember their initial form under the action of stimulants is attractive in engineering applications. SMPs have a wide range of functions in the design of micro-biomedical components, medical devices, smart textiles, sensors, aerospace deployable equipment, and so on. This has resulted in a growing research interest in this field [START_REF] Hu | Recent advances in shape-memory polymers: Structure, mechanism, functionality, modeling and applications[END_REF][START_REF] Monkman | Advances in shape memory polymer actuation[END_REF][START_REF] Behl | Actively moving polymers[END_REF][START_REF] Behl | Shape-memory polymers and shape-changing polymers[END_REF][START_REF] Baer | Shape-memory behavior of thermally stimulated polyurethane for medical applications[END_REF][START_REF] Liu | Review of progress in shape-memory polymers[END_REF][START_REF] Wei | Shape-memory materials and hybrid composites for smart systems-Part I: Shape-memory materials[END_REF][START_REF] Ratna | Recent advances in shape memory polymers and composites: A review[END_REF][START_REF] Golbang | Temporary shape development in shape memory nanocomposites using magnetic force[END_REF][START_REF] Lendlein | Shape-memory polymers for biomedical applications[END_REF][START_REF] El Feninat | Shape memory materials for biomedical applications[END_REF].

The ability of SMPs to change form could be programmed in various ways such as thermo mechanical; chemical, electrical and magnetic. Shape memory effects (SME) can be obtained by a process cycle. In the case of heating stimuli, the process can be summarized in successive steps: heating, deformation, relaxation, fixation and recovery.

A review of the literature gives a great deal of valuable information [START_REF] Gall | Thermomechanics of the shape memory effect in polymers for biomedical applications[END_REF][START_REF] Nguyen | A thermoviscoelastic model for amorphous shape memory polymers: Incorporating structural and stress relaxation[END_REF][START_REF] Li | Thermomechanical behavior of thermoset shape memory polymer programmed by cold-compression: Testing and constitutive modeling[END_REF][START_REF] Heuchel | Relaxation based modeling of tunable shape recovery kinetics observed under isothermal conditions for amorphous shape-memory polymers[END_REF][START_REF] Chen | Influence of thermoviscoelastic properties and loading conditions on the recovery performance of shape memory polymers[END_REF][START_REF] Ferry | Introduction to polymer viscoelasticity[END_REF][START_REF] Srivastava | Thermally actuated shape-memory polymers: Experiments, theory, and numerical simulations[END_REF][START_REF] Tcharkhtchi | Viscoelastic properties of epoxy networks in the glass transition region[END_REF][START_REF] Prisacariu | Perspectives: Novel crosslinked polyurethanes as shape-memory materials[END_REF] about the SMPs; however the answers to a number of serious and foremost questions are still shrouded in mystery. For example, some macroscopic phenomena still lack a detailed understanding on the microscopic level. Thus, there is a need for further fundamental research on SMP and related concepts in order to develop their applications.

The first question is about the definition of SME. Does this effect concern only the macroscopic scale or is it reflected also in the microscopic shape change? In fact, it is usually common to consider SME for macroscopic changes. However, if the capacity to regain the initial (permanent) shape after certain stimulation, even at microscopic scale, is considered as SME, several phenomena like elasticity, shrinkage, swelling, thermal expansion, etc. can also be considered as various expressions of shape memory [START_REF] Zhang | A novel type of shape memory polymer blend and the shape memory mechanism[END_REF][START_REF] Pittini | Cole-cole plot analysis of dielectric behavior of monoalkyl ethers of polyethylene glycol (C n E m )[END_REF][START_REF] Merdas | Factors governing water absorption by composite matrices[END_REF][START_REF] Bronnec | Water absorption characteristics of diglycidylether of butane diol-3,5-diethyl-2,4-diaminotoluene networks[END_REF][START_REF] Tcharkhtchi | Thermal expansion of epoxide-amine networks in the glassy state[END_REF]. The definition of SME, in this way, helps us to understand its mechanism more fully on the microscopic and molecular scale. One of the objectives of this work is to demonstrate the importance of this new concept of SME.

The second question concerns the chemical, physical origins of SME and its mechanisms. It is generally accepted that shape memory effect is related to the presence of two different, hard and soft segments in the structure of polymers. The hard morphology stabilizes the permanent shape of the structure while the soft morphology allows the passage from permanent to temporary state or vice versa. In the case of cross-linked polymers, the cross-linkage points and the segments of the chains directly linked to these points may be considered as the hard zones of the morphology. These zones are responsible for the permanent state of the polymer because of their high stability. Moreover, there is no significant molecular motion in these zones. In the case of a blend, the mixture of two themoplastic polymers, one of the polymers (with higher glass transition temperature) may play the role of hard segment and the other one that of soft segment. The compatibility of these two polymers has a high influence. Even though this explanation is correct, it cannot respond to all of the questions related to the SME mechanism; for example it cannot illustrate the chemical and physical origin of SME. So our objective is also to bring in some new elements to explain the mechanism of SME. For example, we try to show the role of driving force during recovery.

In a study of SMPs one of the questions is about the relationship between SME and properties memory effect (PME). In our previous work [START_REF] Farzaneh | Shape memory effect and properties memory effect of polyurethane[END_REF], we showed that after recovery a SMP does not necessarily exhibit the same properties; that means even after regaining 100% of initial shape, the polymer may not always show the same properties. In the case of thermal stimuli sensible polymers, one of the ways to determine this relationship is to perform cyclic-thermo-mechanical tests. This investigation can be carried out by means of a tensile machine equipped with a thermo chamber.

Figure 1 shows a three-dimensional demonstration of each cycle. As can be seen, each cycle consists of three consecutive steps, stress-strain curve, stress-relaxation time evolution and recovery. This paper also introduces the concept of partial shape memory effect (PSME). It is worth mentioning that between the SMP and those without SME, there are polymers which have an intermediate behavior. This concept is important because it helps us to understand the difference between an SMP and a non-SMP via the behavior of the polymers with PSME.

Materials

In this study, two polymers were used: polyurethane (PU) and a blend of polylactic acid (PLA) and polybutylene succinate (PBS), which is one of the most promising polymers in the family of synthetic biodegradable polyesters. The PU exhibits 100% SME and the PLA/PBS shows PSME.

The blend PLA (80%)/PBS (20%) was made with the following polymers:

Polylactic acid (PLA): PLA, trade name "Ingeo™ biopolymer 2003D", is a NatureWorks LLC product (Blacksburg, VA, USA). It is a semi-crystalline thermoplastic, used naturally or as part of a formulated blend. This is a high molecular weight biopolymer grade that processes easily on conventional extrusion equipment.

Polybutylene succinate (PBS), trade name PBE002 commercialized by Natureplast (Caen, France) is a kind of fully biobased macromolecular polymer which is synthesized from succinic acid and 1,4-butanediol through a direct process of condensation polymerization.

The PLA that we used for this study is brittle, however different experiments showed that [START_REF] Bhatia | Compatibility of biodegradable poly(lactic acid) (PLA) and poly(butylene succinate) (PBS) blends for packaging application[END_REF][START_REF] Harada | Increased impact strength of biodegradable poly(lactic acid)/poly(butylene succinate) blend composites by using isocyanate as a reactive processing agent[END_REF], the blend of up to 20% (by weight) PBS content in PLA is expected to overcome the deficiencies of PLA. Polyurethane (PU), supplied by the RAIGI Company (Rouvray-Saint-Denis, France), is the result of the reaction between an isocyanate (diphenyl-methane di-isocyanate, (MDI) with a functionality of more than two) and a mixture of polyol (polyester with molecular weight of 800 g/mol) and chain extender (another polyol with a molecular weight of 400 g/mol). The catalysts of the types aluminosilicate and diamine were used for the reactions.

Processing Methods and Tests

To prepare the blend the homogenous mixture of PLA (80%) and PBS (20%) was extruded using a twin screw Brabender extruder machine (Brabender, Duisburg, Germany). The transformation temperature ranged from 180 to 200 °C. Extruded films were obtained by pulling the product through a set of cooling rolls.

In order to prepare the samples, extruded films are cut by punch using a manual press. NF ISO 6239 standard [START_REF]Plastics-Determination of Tensile Properties by Use of Small Specimens[END_REF] defines the test conditions and the shape of the samples. Two types of samples were obtained: parallel and perpendicular to the extrusion axis.

In the characterization methods the differential scanning calorimetry (DSC) measurements were carried out with the DSC Q10 V9.0 Build 275 (TA Instruments, Guyancourt, France). The samples of mass between 10 and 20 mg were placed in hermetic aluminum capsules. The sample was heated first up to 200 °C with a temperature rate of 5 °C/min, in order to eliminate the effect of the thermal history of the material. It was cooled down to -60 °C in order to be reheated up to 200 °C with the same temperature rate.

For the characterization of the viscoelastic behavior of PLA/PBS blends, the tests were carried out in tension using a dynamic mechanical analyzer (DMA) Q800 instrument, from the TA Company (Guyancourt, France). These tests were performed under the following condition: frequency 1 Hz, temperature ranging from -80 to 120 °C, temperature rate 2 °C/min, and dynamic force 0.1 N. The sample had a rectangular shape with 25 mm × 4 mm × 2 mm dimensions. Multi-frequency tests were performed on virgin and 100% recovered polyurethane (PU) in order to use the Williams-Landel-Ferry (WLF) equation for the determination of free volume fraction and thermal expansion coefficient.

A shape memory cycle test consists of the following successive steps: tensile test, stress-relaxation test, fixing and recovery test. For reproducibility of the results, the shape memory tests were carried out on at least three samples, each time.

(i) Tensile test: Tensile tests were achieved with the Instron 5881 (Elancourt, France), loading cell of 1 kN, according to the standard NF ISO 6239 [START_REF]Plastics-Determination of Tensile Properties by Use of Small Specimens[END_REF]. The strain rate is 5 mm/min. This machine is equipped with a temperature controlled thermo chamber for heating the samples and performing the tests at different temperatures.

(ii) Stress-Relaxation test: This test was performed at the end of the tensile test with the same machine. The strain was taken as constant (100%, 200% and 250%) and the decrease of stress was measured versus time. This test was performed to eliminate the residual stress induced during the tensile test in the samples. The stress-relaxation tests were performed at 70, 75 and 80 °C.

(iii) Fixing: By cooling the sample with a ventilator at the end of the stress-relaxation test (or tensile test) the dimension of the sample, of the 100%, 200% and/or 250% strain was conserved in the sample.

(iv) Recovery test: The recovery tests are carried out in a thermal chamber at 70, 75 and 80 °C and must follow the changes in size (length) of the samples as a function of time.

Because of the shape memory effect, the specimens regain a percentage of their initial shape. In the multi-cycle method, successive shape memory tests were conducted on the same sample. After the first cycle (at the end of the recovery test) the sample was used for a second cycle in order to carry out the same tests successively; this means tensile, fixing (without relaxation test) and recovery test. The cycle was repeated several times.

Experimental Results

In terms of the thermal characteristics of the blend the results of DSC tests on PLA, PBS and PLA (80%)/PBS (20%) blend are presented in Table 1. It can be seen that the glass transition temperature (T g ) of PBS and PLA are respectively -38 and 58 °C and the melting point (T m ) of PLA and PBS are respectively 154 and 115 °C. These results were confirmed by the same tests on two other samples. The DSC test on the blend shows two distinct T g (-38 and 54 °C) and two distinct T m (115 and 154 °C). These values are nearly the same as the T g and T m of PBS and PLA. This means that these polymers are not miscible in the blend. The same results were obtained by thermo-mechanical analysis (Figure 2). In this case, the loss modulus (E′′) spectra also shows two distinct transition regions, at -30 and 51.5 °C, related respectively to the glass transitions of PBS and PLA. The immiscibility of these two polymers was observed previously by scanning electron microscopy (SEM) [START_REF] Harada | Increased impact strength of biodegradable poly(lactic acid)/poly(butylene succinate) blend composites by using isocyanate as a reactive processing agent[END_REF] with which the separation of two phases was shown clearly. However, it can be seen that by mixing 80% of PLA with 20% of PBS, the glass transition temperature of PLA in the blend decreases clearly from 58 to 54 °C. This decrease of T g value may be contributing to the plasticizing effect of PBS. The results of tensile tests on pure PLA and the PLA/PBS blends confirm this hypothesis (Table 2). The pure polymer PLA with a very low strain at break (5.2%) is brittle at 24 °C, while the PLA/PBS blend (parallel sample) with a strain at break more than 110% is ductile. In fact PBS has a melting point lower than PLA; during extrusion of the blend at a temperature higher than the melting point of PLA, PBS is at very low viscosity; it may partially penetrate between the molecular chains of PLA, like a plasticizer, even it is not miscible in PLA.

Figure 3 compares parallel and perpendicular samples. For parallel samples the value of strain at break is about 110% at room temperature; that means by blending 20% of PBS and 80% of PLA, the brittleness of PLA is improved and it becomes ductile. This is not the case for the perpendicular sample which always remains brittle. The difference between parallel and perpendicular samples may be explained as follows:

(i) During extrusion, the molecular chains of PLA and the discontinuous phase PBS became oriented in the direction of the flow (Figure 4). So, because of this orientation, the parallel sample (E = 925 MPa) is more rigid than the perpendicular sample (E = 855 MPa).

(ii) The standard deviation related to the results obtained from the tests on the parallel samples is more important than the perpendicular samples. In fact, the mechanical properties of the polymer in the center and at the edge of the film are not the same. The samples which are taken from these two zones have different properties. After obtaining these results, only parallel samples were used, taken from the center of the extruded film. This type of sample has high ductility (ε r = 110% at room temperature), necessary for shape memory study.

Shape memory tests were carried out on the parallel extruded samples. In general, each shape memory cycle consists of the following successive tests: tensile, stress-relaxation, fixing and recovery. However, certain cycles were limited to tensile, fixing and recovery tests, without stress-relaxation tests. Tensile tests were performed on parallel samples at 70 °C. The result is shown in Figure 5. The sample has two different linear deformation zones; elastic zone (8%-10% strain) and plastic zone (>10% strain). At 100% strain the value of stress is about 4.3 MPa.

The stress-relaxation isothermal test was performed at the end of the tensile test with the same machine. The strain was taken constantly and the decrease of stress was measured versus time. This test was performed to eliminate the residual stress induced during the tensile test in the samples. The stress-relaxation tests were performed at 70, 75 and 80 °C (Figure 6). The rate of relaxation and the value of stress at the end of the test depend on the temperature. The rate of decrease of stress is more important when the temperature is higher. The value of the stress at the end of the relaxation test at 70 °C is more important than 80 °C. This value is the same at 75 and 80 °C.

Fixing can be performed after tensile tests or after stress-relaxation tests. During this step, the sample obtains a temporary shape. By cooling the sample with a ventilator, the maximum deformation (100%, 200% and/or 250%) obtained at the end of tensile tests is conserved in the sample.

The recovery tests are carried out in a thermal chamber at 70 and 75 °C after fixing at 100% of deformation. The percentage of recovery was determined using Equation (1):

100 (%) 0 1 2 1 L L L L R (1)
where, L 0 is the initial length of sample, L 1 represents the length of sample after tensile tests and L 2 is the length of sample after recovery test.

The results (Figure 7) show that the percentage of recovery in two cases (70 and 75 °C) is about 24%. But the recovery rate at 75 °C is higher than the recovery rate at 70 °C. For example after 5 min, the sample regains only 7% of its initial shape when the temperature is 70 °C, but this value is more than 21% (three times more) when the temperature is 75 °C.

The viscoelastic behavior of the blend is determined by comparing the thermodynamic behavior of the SMPs at different stages of the process. These results help to find the relationship between shape memory effect and the viscoelastic properties of the polymer. In fact, by increasing the temperature, the molecular motion becomes more and more important. The test on the stretched sample is especially interesting because, in the glass transition zone, the polymer tries to regain its initial shape. The results show that the stretching, relaxation test and fixing increase the value of the glass transition temperature. However, the comparison of the DMA spectrum of the samples before and after the recovery test reveals that the T g of the blend decreases during the recovery test. In parallel samples, before the recovery test, the T g is about 63.7 °C while after the recovery test it shifts to 62.3 °C. Indeed before the recovery test, the sample is stretched and its molecular mobility is decreased because of this stretching.

Analysis of a Multi-Cycle Shape Memory Test

Multi-cycle shape memory tests were performed on the same sample according to the procedure, explained before. The results of the tensile tests (at 70 °C) of the different cycles are presented in Figure 8. In this figure, real stress is plotted versus real strain.

These tests were repeated at least with two other samples; the same results were obtained. The different cycles may be analyzed as follows:

(1) Cycle 1 can be divided into two zones with: a. Elastic stage up to 6%-8% of deformation and with a Young's modulus of 50 MPa; b. Plastic linear stage up to 69% of real strain. After recovery test, the sample regains 24% of deformation (6%-8% of elastic part and 16%-18% of plastic part). (2) Cycle 2 begins at the end of the first cycle and after the recovery test. This cycle is obviously divided into three linear stages. a. First stage up to 3%-4% deformation and with a modulus of 110 MPa; b. Second stage up to 12% of real strain and with a modulus of 41 MPa. The stress-strain curve at this stage is relatively linear; c. Third stage with a modulus of 11 MPa. In this stage, the stress-strain curve remains first linear up to 30% of deformation; then its slope increases. (3) Cycle 3 begins at the end of the second cycle and after recovery test. This cycle is also divided into three linear stages. a. First stage up to 3%-4% of deformation and with a modulus of 115 MPa;

b. Second stage up to 20% of real strain and with a modulus of 42 MPa. The stress-strain curve at this stage is relatively linear; c. Third stage with a modulus of 20 MPa. In this stage, the stress-strain curve is relatively linear up to 45% of deformation. Then its slope increases of the molecular orientation. (4) Cycle 4 begins at the end of the third cycle and after recovery test. As in previous cycles, it is separated into three different stages. a. First stage up to 3%-4% of deformation and with a modulus of 130 MPa; b. Second stage up to 30% of deformation and with a modulus of 44 MPa. In this stage, the stress-strain curve is not linear; the increase of slope may be related to the molecular orientation; c. Third stage up to 50% of real strain and with a modulus of 37 MPa. In this stage, the stress-strain curve is relatively linear.

These results can be summarized by the following comments:

(1) It can be seen that during the first cycle, the stress-strain curve is a classical one. From the second cycle, the obtained tensile curves show three stages. (2) The slope measured during Stage 1, assumed to be a Young's modulus, E 1 increases with the number of cycles, from 50 to 130 MPa (Figure 9 and Table 3). This increase in the stiffness of the blend may be related to the reinforcement of the homogenous matrix, PLA phase (hard segment). The orientation of molecular chains in the PLA structure may be presumed as the cause of this reinforcement. However the change of E 1 between first and second cycles is more important (from 50 to 110 MPa). From the second cycle, E 1 changes less and the variation of its value is not very significant. (3) Regarding the formation and development of Stage 2, several remarks on the constitutive tensile curve should be made:

(a) The slope of the curves related to this stage remains in the same order of magnitude (41-44 MPa) and its value is nearly the same as that of the original stiffness of the material in cycle 1 (50 MPa); (b) The extent of this stage increases significantly with the number of cycles; (c) The slope of Stage 3 increases also with the number of cycles. In this stage, the last part of the curves presents a transition zone of kinetic change. After this zone, the material becomes stiff. This is classically attributed to the orientation of the macromolecular chains. One should draw attention to the fact that the level of stress measured at this point is always of the same order of magnitude as the level of stress reached at the transition ("yield point") between Stages 2 and 3 at the next cycle. For example, the transition zone of curve 1 (around 6 MPa) corresponds to the level of the stress observed at the "yield point" (between Stages 1 to 2) of the curve 2.

For a first approach, the Stage 2 of the curve 2 may be considered also as the expression of an elastic behavior. Indeed, in order to understand the origin of the appearance and development of Stage 2, we propose here an original phenomenological interpretation of the consecutive tensile curves. Let us first take a classical case. When the polymer is under plastic deformation (and more generally "anelastic deformation"), we can assume that the total strain (ε tot ) can be divided in two parts, elastic (ε e ) and plastic (ε p ) deformations:

tot e p ε = ε + ε (2)
In the case of linear elasticity:

e σ ε E (3)
At the end of the first cycle, during fixing, the temperature decreases under the glass transition zone and the sample passes from a rubbery state (70 °C) to a glassy state (room temperature) leading to an increase of the Young's modulus. Thus the elastic part decreases (Equation ( 3)) and tends to zero, leading to an increase of plastic (anelastic) strain which tends towards the imposed total strain (Equation ( 2)). So the yield stress on the curve related to cycle 2 corresponds to the stress at the end of fixing during cycle 1. On the other hand, from the phenomenological point of view, fixing leads to a hardening of the polymer.

When performing recovery, the specimen is submitted to a driving force. Figure 10 shows the elevation of temperature up to the glass transition which leads to a decrease of the stiffness simultaneously, due to decrease in the plastic strain (equilibrium between plastic and elastic strain, see Equation ( 2)), where the elastic part of the strain is released. Moreover the extent of Stage 2 increases significantly with the number of the cycle. This phenomenon may be understood by the hardening effect due to the fixing as it has been explained by the fact that the level of stress at this point is always of the same order of magnitude as the "yield point" of the previous cycle. This expression of a second elastic phase can be explained by different physical phenomena as, for example, the local effect due to the multi-phased configuration. Indeed, the localization of the deformation around the heterogeneity coupled with chain orientation can produce this kind of effect. This will be discussed in a further paper through the results of a micromechanical model using the concept of a coated heterogeneity.

(4) Additionally, the presence of a transition zone on Stage 2 for curve 4 (cycle 4) can be noticed. It seems after several cycles the orientation of the chain has been greatly facilitated and can be activated sooner even during cycle 2.

Shape Memory Effect

Some new concepts of shape memory effect are discussed by looking at the origin of the recovery and driving force. For example, partial shape memory effect, relationship between the shape and properties of SMP, and shape memory in micro-scale domains.

Origin of the Recovery and Driving Force

Because of the temperature increase during the recovery test, the stiffness of the polymer decreases leading to the decrease of the elastic (or viscoelastic) part of the strain. Consequently, the partition of the strain of Equation (1) leads to the decrease of the plastic strain. The release of the elastic part of the strain contributes partially to the recovery.

The driving force responsible for the shape memory effect is in fact the residual stress introduced in the sample during the tensile test.

In the case of polymers without shape memory, this deformation in the plastic zone is irreversible. The sample will recover only the elastic portion of deformation when it is released in the plastic zone. In this case if the test temperature is below the T β and the polymer is in its glassy state, due to a low mobility, the elastic deformation will be small and the sample remains relatively far from its original form. However, if the test temperature is above the T β and especially when the polymer is in a rubbery state, with a high molecular mobility, the elastic deformation may be very important and when the sample is released at the end of the tensile test, it can regain an important part of its original shape.

In the case of shape memory polymers, after the tensile test is at 100% of deformation, a stress will be introduced in the sample. By the relaxation test (Figure 11) after the tensile test, two portions of residual stress can be distinguished.

Figure 11. Schematic relaxation curve

The first portion, σ r1 , will be eliminated at the end of the relaxation test; but the second portion of this residual stress σ r2 , remains in the structure of the sample. This stress is the origin of shape memory effect. If the shape memory cycle test is performed without the relaxation test, as the value of residual stress is higher, the SMP regains its initial form faster. In any case the SMP regains its initial shape with or without the relaxation test. The residual stress then can be considered as the driving force of a SMP for regaining the initial shape.

In order to show the role of this driving force, the residual stress is recapitulated and plotted versus the percentage of recovery respectively in Table 4 and Figure 12. These results were obtained from the shape memory tests with and without stress-relaxation experiments.

Partial Shape Memory Effect

As we have seen before, the blend under study (PLA (80%)/PBS (20%)) is not a 100% shape memory polymer. At the beginning the blend has a partial shape memory effect with a recovery rate of 24%. Figure 13 shows that after each cycle, the blend has been modified and transformed to a new blend with SME. This effect increases after each cycle; 42% after the second cycle, 71% after the third cycle and reaches 82% after the fourth cycle. The evolution of the shape memory capacity of a polymer which is not originally a SMP, can reinforce the following ideas: As it was explained before, logically a polymer with partial shape memory effect may be transformed into a 100% SMP. However, this idea must be confirmed by experiments. In any case the improvement of shape memory properties of a PSMP is very attractive for industry.

Shape Memory in Micro-Scale Domains

In the previous paper [START_REF] Farzaneh | Shape memory effect and properties memory effect of polyurethane[END_REF], it was shown by dynamic mechanical analyzer (DMA) tests on polyurethane (with shape memory effect) that after recovery test and even after regaining 100% of the initial shape, the sample does not have the same properties. It was seen that, G′ (storage modulus) and G′′ (loss modulus) of recovered and virgin samples were different with different alpha transition temperatures. These results were followed by the successive cycles tests, performed on the same sample. The results indicated that the Young's modulus of the virgin sample of SMP and the same sample after recovery test was not the same.

Another physical property which is related to the microstructure is free volume, shown by free volume fraction f g . This parameter may be determined by the results of multi-frequency DMA tests on polyurethane (Figure 14), using the Williams-Landel-Ferry (WLF) equation:

1 α αr r2α αr () log () C T T F F C T T (4)
where F = frequency, T α = α transition temperature, F r = reference frequency, T αr = reference temperature, 1g 2.303

C B f , 2g α Cf , 1 B
, α = thermal expansion coefficient and g f = free volume fraction.

In order to verify the validity of this equation, the linear regression method was used and 5. These results show clearly that T g , ∆α and f g of virgin and recovered (100%) samples are not the same. That means that the PU under study shows SME but not memory effect for certain properties. After recovery tests and regaining initial shape, it is possible that the sample does not regain some of its initial properties because of microstructure irreversible changes. These microstructure changes may induce the micro-change of volume which cannot be regained after the recovery test. In a study of shape memory effect, generally, the authors measure the macro-change of volume. They do not give enough attention to the micro-change of volume. In fact on the macroscopic scale, the measurement of micro-changes is practically impossible.

The difference has been explained by the presence of two types of properties:

The properties like density or Poisson's ratio which are related directly to the macro-structure (shape) of the sample. These properties may be regained after recovery tests. The properties like modulus or glass transition temperature which is related to the micro-structure (in molecular or macromolecular scale). During shape memory test, the mechanism of their evolution is not the same as the mechanism of the shape memory effect. So these properties may not be regained totally after the recovery test.

Conclusions

This study proposes certain new concepts for shape memory polymers. For properties memory effect (PME) from the results of DSC tests, tensile tests, DMA test, and by using the WLF equation on virgin and 100% recovered polyurethane, we determined the value of: modulus, strength, glass transition temperature, heat coefficient expansion and free volume fraction. The results show that they do not have the same properties indicating that shape memory effect and properties memory effect may have two different mechanisms.

For the partial shape memory effect (PSME), the multi-cycle shape memory tests show that blend PLA (80%)/PBS (20%), regains only 24% of the initial shape at the end of the first cycle, indicating that certain polymers have a partial shape memory effect. This test also showed the shape capacity of polymer increases by the number of cycles. The same polymer after the second, third, and fourth cycles regains respectively, 42%, 71% and 82% of its initial shape. This new concept of shape memory effect is of very important value. We can in a legitimate way suppose that a polymer without shape memory effect can be transformed into a shape memory polymer.

Finally, for the shape memory effect in micro, as certain properties of polymers can be affected by the micro change of the shape (elasticity, free volume, etc.), it is very important to study the shape memory effect on the micro-scale. In fact, the shape memory effect has always been studied macroscopically. This study may not in itself be sufficient to fully explain the relationship between shape memory effect and properties memory effect.

Figure 1 .

 1 Figure 1. A tridimensional presentation (stress-time-strain) of one cycle of shape memory tests on polyurethane [26].

Figure 2 .

 2 Figure 2. Storage (E′) and loss modulus (E′′) spectra of polylactic acid (PLA)/polybutylene succinate (PBS) blend.

Figure 3 .

 3 Figure 3. Stress-strain curves of PLA (80%)/PBS (20%) at 24 °C: perpendicular (1) and parallel (2) samples.

Figure 4 .

 4 Figure 4. Molten film of PLA/PBS blend during extrusion (1); perpendicular and parallel samples (2).

Figure 5 .

 5 Figure 5. Nominal stress-strain curves of parallel sample at 70 °C.

Figure 6 .

 6 Figure 6. Relaxation tests at 70 °C (1); 75 °C (2) and 80 °C (3).

Figure 7 .

 7 Figure 7. Recovery tests at 70 °C (•) and 75 °C (♦).

Figure 8 .

 8 Figure 8. Real stress versus real strain (tensile tests) at the beginning of each cycle of shape memory tests at 70 °C. (1) First cycle; (2) Second cycle; (3) Third cycle; and (4) Fourth cycle.

  Stage 1 seems to correspond to an elastic behavior. Stage 2 is an intermediate behavior stage. Stage 3 follows Stage 2 until the end of the test;

Figure 9 .

 9 Figure 9. Relationship between tensile modulus and recovery percentage related to the first zone (E 1 ) and the second and third zones (E 2 and E 3 ) of stress-strain curves (at 70 °C) of different successive cycles.

Figure 10 .

 10 Figure 10. Variation of modulus versus temperature (fixing and recovery).

Figure 12 .

 12 Figure 12. Percentage of recovery versus residual real stress.

Figure 13 .

 13 Figure 13. Evolution of shape memory effect (SME) of PLA/PBS by succesive shape memory tests.

  If this plot is linear having 21 S C C as slope then the validity of the WLF equation is verified and f g may be determined by the following equation: this method, we performed multi-frequency 3 point bending DMA tests on PU. The test conditions are as follows: frequency range 0.1-100 Hz, temperature range 25-80 °C. The results are shown in Table

Figure 14 .Table 5 .

 145 Figure 14. Multi-frequency 3 point bending dynamic mechanical analyzer (DMA) test on virgin polyurethane (PU): (1) 0.1 Hz; (2) 1 Hz; (3) 10 Hz; and (4) 100 Hz.

Table 1 .

 1 Values of glass transition temperature (T g ) and melting point (T m ) of polylactic acid (PLA), polybutylene succinate (PBS) and the blend PLA/PBS.

	Polymer

T g (°C) T m (°C)

	Pure polylactic acid (PLA)	58	154
	PLA 80% + PBS 20%	PLA PBS	54 -38	154 115
	Pure polybutylene succinate (PBS)	-38	115

Table 2 .

 2 Results of tensile tests at 24 °C, on virgin PLA and blend of PLA/PBS (parallel samples).

	Polymer	E (

MPa) σ y (MPa) ε y (%) σ r (MPa) σ r (%)

	PLA	Average Standard deviation	1340 165	--	--	53 7.5	5.2 0.40
	PLA/PBS blend	Average Standard deviation	925 99	45 6.0	4.8 0.50	31.5 0.96	110 8.5

Table 3 .

 3 Relationship between percentage of recovery and modulus related to the different stages of different cycles of multi-cycle tests.

	Cycle No.	Stage	Residual stress (MPa) Recovery (%) E 1 (MPa) E 2 (MPa) E 3 (MPa)
	0	Initial	0	-	-	-	-
	1	before recovery after recovery	7.2 6.5	24%	50	-	5
	2	before recovery after recovery	13.1 9.5	42%	110	41	11
	3	before recovery after recovery	21 11	71%	115	42	20
	4	before recovery after recovery	22.9 6.5	82%	130	44	37

Table 4 .

 4 Relationship between residual stress and percentage of recovery.

	No

. Residual real stress (MPa) Temperature of recovery tests (°C) Recovery (%)

	1	1.22	80	4.72
	2	1.25	80	5.17
	3	3.81	80	17.14
	4	3.83	80	17.20
	5	5.37	80	18.20
	6	5.51	80	19.00
	7	7.10	70	19.12
	8	7.41	70	19.21
	9	7.73	70	20.51
	10	8.22	75	20.86
	11	8.47	70	23.50
	12	8.73	80	22.30
	13	7.78	70	22.35
	14	9.39	70	23.00
	15	11.02	75	24.00
	16	11.69	75	24.34
	17	29.32	70	42.00
	18	56.61	70	70.90
	19	90.17	70	81.80

Acknowledgments

The polyurethane (PU) samples used for this study have been supplied by RAIGI (Rouvray-Saint-Denis, France). The authors thank the company for its precious help.

Conflicts of Interest

The authors declare no conflict of interest.


A Belmiloudi

F Mahé

On nonlinear inverse problems of heat transfer with radiation boundary conditions. Application to dehydratation of gypsum plasterboards exposed to fire

Keywords: Optimal control, model calibration, numerical approximation, control constraint, adjoint model, dehydratation of gypsum, numerical stability analysis

The paper investigates boundary optimal controls and parameter estimates to the well-posedness nonlinear model of dehydratation of thermic problems. We summarize the general formulations for the boundary control for initial-boundary value problem for nonlinear partial differential equations modeling the heat transfer and derive necessary optimality conditions, including the adjoint equation, for the optimal set of parameters minimizing objective functions J. Numerical simulations illustrate several numerical optimization methods and several examples and realistic cases, in which several interesting phenomena are observed. A large amount of computational effort is required to solve the coupled state equation and the adjoint equation (which is backwards in time), and the algebraic gradient equation (which implements the coupling between the adjoint and control variables). The state and adjoint equations are solved using the finite element method.

Introduction and mathematical setting of the problem 1.Motivation

Since several years, a considerable effort has been made to develop materials having a good fire resistance. Such materials must provide a sufficient mechanical resistance to avoid the premature collapse of a building structure undergoing a fire. Consequently this type of material must withstand significant heating without burning and keep its mechanical resistance sufficient. Criteria which permit to appreciate the fire resistance of materials are given by several norms which define the minimum fire exposure duration that must support the building structure.

One of the building materials presenting the best fire resistance is gypsum plasterboard, which in turn is due to the dehydration phenomenon. This material presents the particularity to undergoing two chemical reactions of dehydration during its heating. These two endothermic reactions considerably slow down the heating of the material (since the dehydratation process consumes large amount of heat), and provides to the plasterboard an excellent fire resistance.

The main particularity of gypsum plasterboard is that it contains 21% of chemically bound water by weight. When the temperature reaches 100 o C in a point of the plasterboard, a reaction of dehydration occurs in the material. This chemical reaction dissociates a certain quantity of water which is combined to the crystal lattice. In terms of fire safety, the reaction of dehydration and the vaporization of free water absorb a certain amount of energy which significantly slows down the heating of the material and in particular the temperature rise on the unexposed side of plasterboard. (e-mail: Fabrice.Mahe@univ-rennes1.fr)

1
Although necessary, experimental testing is not a convenient way to estimate the fire resistance of a material. Indeed, full-scale testing poses the problem of the high cost of the experimental setup and the difficulty to implementing the experiment. In addition pilot-scale testing does not allow to accurately reproducing the real conditions of a fire exposure. Consequently, the development of a mathematical model, and the numerical simulation of the heating of gypsum plasterboard exposed to fire appears as a suitable means to study the thermal behaviour of the material during a fire exposure.

Thermo-chemistry of gypsum

Gypsum plasterboard is commonly used as construction material to improve fire resistance of building structures. The pure Gypsum, existing at the natural state as a more or less compact rock, is composed of calcium sulphate with free water at equilibrium moisture (approximately 3%) and approximately 20% per weight of chemically combined water of crystallization (see, e.g., [START_REF] Chang | Application of thermo-Raman spectroscopy to study dehydration of CaSO 4 .2H 2 O and CaSO 4 .0.5H 2 O[END_REF][START_REF] Freyer | Crystallization and phase stability of CaSO 4 and CaSO 4 -based salts[END_REF][START_REF] Gerlich | Buchanan Design of light steel-framed walls for fire resistance[END_REF]). Its chemical formula is CaSO 4 .2H 2 O (calcium sulphate di-hydrate). The industrials add various chemical elements (in small quantities) in order to increase their performance when exposed to elevated temperatures. The chemical reaction which consists in removing chemically combined water of crystallization is called calcination. During heating, gypsum plaster undergoes two endothermic decomposition reactions. The first dehydration reaction occurs at approximately 100 -120 o C when the calcium sulphate di-hydrate is converted to calcium sulphate hemihydrate (the reaction is always complete by 160 o C) as shown by the following reaction

CaSO 4 .2H 2 O -→ CaSO 4 . 1 2 H 2 O + 3 2 H 2 O (1.1)
The amount of energy required by this first dehydration is about 500kJ per kg of gypsum (see [START_REF] Mehaffey | A model for predicting heat transfer though gypsum-board/wood-stud walls exposed to fire[END_REF]).

The second dehydration reaction occurs when calcium sulphate hemihydrate is converted to calcium sulphate anhydrate as shown by the following reaction

CaSO 4 . 1 2 H 2 O -→ CaSO 4 + 1 2 H 2 O (1.2)
The amount of energy corresponding to this second reaction is about 169kJ per kg (see [START_REF] Axenenko | The modelling of dehydration and stress analysis of gypsum plasterboards exposed to fire[END_REF]).

Remark 1.1 Other amounts of energy can be found in the literature for these reactions, see [START_REF] Thomas | Thermal Properties of Gypsum Plasterboard at High Temperatures[END_REF] for a review.

Both reactions are endothermic, produce water vapour and absorb a large amount of energy. The effect of the endothermic reactions on the heating of the wall of plasterboard is taken into account by including the latent heats of reactions (1.1) and (1.2) in the specific heat evolution. The first dehydration reaction occurs at approximately 100 -160 o C, on the other hand there is some controversy to when the second dehydration reaction occurs. Andersson et al. [START_REF] Andersson | Analytical fire design with gypsum-A theoretical and experimental study[END_REF] (for example) estimate that the second reaction occurs between 210 o C and 300 o C. That consists in introducing two peaks in the evolution of the specific heat according to the temperature, corresponding to the temperatures to which the reactions occur. The areas under the two peaks are equal to the latent heats of the two chemical reactions. Other experiences show that this second reaction occurs immediately after the first one [START_REF] Ghazi Wakili | Gypsum board in fire modeling and experimental validation[END_REF]. In the numerical examples, we will choose this model with only one peak between 100 and 170 o C in the evolution of the specific heat. The information on the thermophysical properties of gypsum plasterboard, at high temperatures, are difficult to measure and then are limited, because the derived results are always complicated by the dynamic nature of the (fire resistive) materials, and vary considerably with the used method of measurement (a wide variety of experimental techniques exists for measuring these properties), and the rate of temperature change (for more details see, e.g., [START_REF] Alfawakhiri | Fire Resistance of Loadbearing Steel-Stud Walls Protected with Gypsum Board: A Review[END_REF][START_REF] Cooper | The Thermal Response of Gypsum-Panel/Steel Stud Wall Systems Exposed to Fire Environments : A Simulation for use in Zone Type Fire Models[END_REF][START_REF] Sultan | A model for predicting the heat transfer though noninsulated unloaded steelstud gypsum board wall assemblies exposed to fire[END_REF][START_REF] Thomas | Thermal Properties of Gypsum Plasterboard at High Temperatures[END_REF]).

Remark 1.2

The model of dehydration can be completed by other reactions at high temperatures. For example, a third endothermic reaction, for decomposition of CaCO 3 , is introduced at 650 o C in [START_REF] Ghazi Wakili | Gypsum board in fire modeling and experimental validation[END_REF] or 800 o C in [START_REF] De Korte | Calculation of thermal conductivity of gypsum plasterboards at ambient and elevated temperature[END_REF]. In [START_REF] Manzello | Influence of gypsum board type (X or C) on real fire performance of partition assemblies[END_REF], the authors note an exothermic reaction around 400 o C.

Outline of the paper

The paper is organized as follows. In the next subsection, we give a sketch of the modeling leading to problem and we establish the governing equations. The model is presented for dehydratation of gypsum plasterboards exposed to fire but the model is general and the method is valid for other applications. In Section 2 we give a description of the parameter estimates (identification problems) as nonlinear optimal control problems with boundary control. This includes results concerning the existence of the optimal solutions, necessary optimality conditions (necessary to develop numerical optimization methods), the optimization problem and adjoint model. Section 3 contains details of the computational algorithm and numerical simulation-optimizations of the optimal control problems. Numerical results for several examples are presented and a realistic situation is analyzed. Section 4 contains a summary and a discussion of future work.

Modeling of the wall of plasterboard heating and the direct forward model

This section is devoted to an introduction of the derivation of dehydration of gypsum plasterboards (exposed to fire) model we study. It is well known that the problem of heat and mass transfer in plasterboard exposed to fire is essentially one-dimensional, so the model is derived in one-dimensional formulation.

Let us consider a wall of plasterboard exposed to fire, which is located vertically on a retaining structural frame. The left hand side is exposed to a heat source, as may occur in furnace in which fire tests are conducted. We suppose that the depth and the width of the wall are much bigger than the thickness (L). Therefore heat fluxes in lateral (y) direction and vertical (z) direction can be neglected in front of the heat flux in direction x. On the other hand, the heat source is distributed uniformly on the heated side of the plasterboard. Applying these physical considerations, the heat transfer can be treated as 1-dimensional process in xdirection. The one-way heat transfer through a plane wall is described by the heat equation in its one-dimensional form

λ(u) ∂u ∂t (x, t) - ∂ ∂x (κ(u) ∂u ∂x )(x, t) = 0, in Q = Ω × (0, T ), u(x, 0) = u 0 (x), x ∈ Ω, (1.3)
where Ω is the boundary subset (0, L), u is the temperature, λ = ρc, ρ is the density of body material, c is the specific heat and κ is the thermal conductivity. The functions λ and κ are variables, positives and bounded.

The external surface of the wall of plasterboard exposed to fire receives a heat flux which consists of convective and radiative components. Consequently, the boundary condition on this side is written, for t in (0, T ), as:

-κ(u(0, t)) ∂u ∂x (0, t) = h 0 (u f (t) -u(0, t)) + ϕ 0 (|u f (t)| 3 u f (t) -|u(0, t)| 3 u(0, t)), (1.4)
where, for all t ∈ [0, T ], h 0 (t) > 0 is the convective heat transfer coefficient between the furnace and the plasterboard surface, u f (t) > 0 is the furnace temperature, ϕ 0 (t) = σǫ 0 (t) > 0, σ = 5.67.10 -8 W/m 2 K 4 is Stefan-Bolzmann's constant and ǫ 0 is effective emissivity of the surface.

The external surface of the plasterboard, which is not exposed to fire, transfers heat to the surroundings by means of convection and radiation. As the previous case, the boundary condition on this side is written, for t in (0, T ), as:

κ(u(L, t)) ∂u ∂x (L, t) = h 1 (u a (t) -u(L, t)) + ϕ 1 (|u a (t)| 3 u a (t) -|u(L, t)| 3 u(L, t)), (1.5)
where , for all t ∈ [0, T ], h 1 (t) > 0 is convective heat transfer coefficient between the surroundings and the plasterboard surface, u a (t) > 0 is the surroundings temperature , ϕ 1 (t) = σǫ 1 (t) > 0 and ǫ 1 is effective emissivity of the surface.

Remark 1.3 In the sequel, to simplify the notations, we denote u for u(x, t), u(0) for u(0, t), u(L) for u(L, t), u(T ) for u(x, T ), u f for u f (t) and u a for u a (t).

We assume that: (H1) functions f = (h 0 , h 1 ) ∈ (L 2 (0, T )) 2 and φ = (ϕ 0 , ϕ 1 ) ∈ (L 2 (0, T )) 2 and satisfy the following pointwise constraint

0 ≤ a i ≤ ϕ i ≤ b i a.e. t ∈ (0, T ), f or i = 0, 1 0 ≤ c i ≤ h i ≤ d i a.e. t ∈ (0, T ), f or i = 0, 1, (1.6)
for some positive constants a i , b i , c i , d i , i = 0, 1, (H2) operators κ(v), λ(v) are sufficiently regular with 0 < K 0 ≤ κ(v) ≤ K 1 and with 0 < C 0 ≤ λ(v) ≤ C 1 for any v > 0.

Remark 1.4 1. Emissivity of a material is defined as the ratio of energy radiated by a particular material to energy radiated by a black body at the same temperature. It is a dimensionless quantity (i.e. a quantity without a physical unit).

2. In the physical case there is not absolute values under the boundary conditions (since the temperature is non negative). For real physical data and operators λ and κ, we can prove by using the maximum principle that the temperature is positive and then we can remove the absolute values. It is useful for the mathematical study of the problem.

3. It is clear that we can obtain in the same way the model in N-dimensional for N ≤ 3 as follows

λ(u) ∂u ∂t -div(κ(u)∇u) = 0, in Q, u(., 0) = u 0 , in Ω, κ(u)∇u.n = h 0 (u f -u) + ϕ 0 (|u f | 3 u f -|u| 3 u), in Σ f = Γ f × (0, T ), κ(u)∇u.n = h 1 (u a -u) + ϕ 1 (|u a | 3 u a -|u| 3 u), in Σ c = Γ c × (0, T ), u = u a , in Σ s = Γ s × (0, T ), (1.7)
where

Ω ⊂ IR N is an open bounded domain with boundary Γ = ∂Ω, Γ = Γ c ∪ Γ f ∪ Γ s such that Γ c ∩ Γ f ∩ Γ s = ∅,
and n is the outward unit normal vector on Γ. Boundaries Γ c and Γ f denote the cold side and the fire side, respectively, and Γ s denotes the other surface of the plasterboard. The well-posedness of problem (1.7) in 3-dimensional, can be obtained in similar way as in [START_REF] Belmiloudi | Analysis of the Impact of Nonlinear Heat Transfer Laws on Temperature Distribution in Irradiated Biological Tissues: Mathematical Models and Optimal Controls[END_REF].

A priori, most researchers who have worked on the modeling of the behavior of gypsum board (literature in the public domain in this field is sparse, see, e.g., [START_REF] Andersson | Analytical fire design with gypsum-A theoretical and experimental study[END_REF][START_REF] Fuller | Temperature distribution in a nailed gypsumstud joint exposed to fire[END_REF][START_REF] Mehaffey | A model for predicting heat transfer though gypsum-board/wood-stud walls exposed to fire[END_REF][START_REF] Axenenko | The modelling of dehydration and stress analysis of gypsum plasterboards exposed to fire[END_REF][START_REF] Gerlich | Buchanan Design of light steel-framed walls for fire resistance[END_REF][START_REF] Sultan | A model for predicting the heat transfer though noninsulated unloaded steelstud gypsum board wall assemblies exposed to fire[END_REF]), have assumed the convective heat transfer coefficients h 0 and h 1 and the relative emissivity ǫ 0 as constants or/and neglected the relative emissivity ǫ 1 on cold surface. The choice of a constant for these coefficients is not a good physical representation of plasterboard exposed to fire. Moreover, the convective heat transfer coefficients and emissivities depend, among other things, on the state of the external surfaces. During the exposition to fire, mechanical resistance of the external surfaces decreases and that causes appearance and growing of crazes (degradation). These modifications of surface states of the external surfaces modify convective and radiative heat transfers. Consequently, the resultant emissivities and convective heat transfer coefficients depends on temperature and large uncertainties exist in regard to the quality of the data reported. Moreover the work of Belmiloudi and Le Meur [START_REF] Belmiloudi | Mathematical and numerical analysis of dehydratation of gypsum plasterboards exposed to fire[END_REF] shows very clearly that the radiative heat transfer between the unexposed surface and the surrounding can not be neglected. Then it is necessary to estimate the convective heat transfers and the emissivity coefficients.

To satisfy this requirement, we estimate these parameters, by using inverse problem techniques as optimal control methods. It is clear that the accuracy of the parameter estimate from furnace, fire test data and target observations (or measurement results) depend significantly on the thermal characteristics of a furnace, on the geometry of the studied element and on the input thermal properties of the material. So, it is important to have a consistent set of values for these data.

The inverse problem formulation 2.1 Problem formulation

We assume that there exists a unique solution u of problem (1.3), (1.6) with boundary conditions (1.4), (1.5), under some hypotheses for the data and some regularity of the operators λ and κ, satisfying the following regularity (by using [START_REF] Belmiloudi | Mathematical and numerical analysis of dehydratation of gypsum plasterboards exposed to fire[END_REF] and [START_REF] Ladyzhenkaya | Linear and quasilinear equations of parabolic type[END_REF]):

u ∈ L 2 (0, T ; H 1 (Ω)) ∩ L ∞ (0, T ; L 2 (Ω)), ∂u ∂t ∈ L 2 (0, T ; (H 1 (Ω)) ′ ), |u(0)| 3 u(0) ∈ L 5 4 (0, T ) and |u(L)| 3 u(L) ∈ L 5 4 (0, T ).
(2.1)

Introduce now the mapping F which maps the source term h of (1.3)-(1.5) into the corresponding solution u = F(h), where h = (h 0 , h 1 ). In this section we formulate the inverse problem as an optimal control problem. The control procedure consists of finding the optimal controls h opt = (h opt 0 , h opt 1 ) and the corresponding optimal temperature u opt = F(h opt ) which minimize a cost criterion J(h). The cost functional J measures the distance between a measured temperature (the observations u obs , u 0 obs , u L obs and u T obs ) and the corresponding predicted temperature obtained from the primal (or direct) model (1.3)-(1.5). Precisely we will study the following optimal control problem (P c ): find (u, h) such that the following objective function

J(h) = 1 2 Q |C(u -u obs )| 2 dxdt + 1 2 Ω |D(u(T ) -u T obs )| 2 dx + 1 2 T 0 γ 0 |u(0) -u 0 obs | 2 + γ 1 |u(L) -u L obs | 2 dt + 1 2 T 0 |Mh| 2 dt is minimized with respect to h ∈ Y ad subject to (1.3)-(1.5), (2.2)
where Y ad is the set of admissible controls, u = F(h), γ 0 and γ 1 are predefined nonnegative weights, M = diag( √ α 0 , √ α 1 ) are predefined nonnegative weights such that α 0 + α 1 = 0. The operators C and D are unbounded operators on L 2 (Ω) satisfying (∀v ∈ H 1 (Ω)) We consider here, an optimal control problem with Tichonov regularization terms [START_REF] Tikhonov | Nonlinear Ill-Posed Problems, V. 1, V. 2[END_REF]. The nonnegative weights α 0 and α 1 play the role of regularization parameters and take small values in the numerical simulations.

||Cv|| 2 L 2 (Ω) ≤ δ 1 ||v|| 2 L 2 (Ω) + δ 2 ||v|| 2 H 1 (Ω) and ||Dv|| 2 L 2 (Ω) ≤ δ 1 ||v|| 2 L 2 (Ω) + δ 2 ||v|| 2 H 1 (Ω) , with δ 1 + δ 2 > 0 and δ i ≥ 0, for i = 1, 2.
According to (1.6) the set of admissible controls describing the constraint is

Y ad = {h = (h 0 , h 1 ) ∈ (L 2 (0, T )) 2 such that 0 ≤ c i ≤ h i ≤ d i a.
e. in (0, T ), f or i = 0, 1}.

First-order necessary conditions

Assume that the nonlinear control problem (2.2) admits an optimal solution (for similar result see, e.g., [START_REF] Belmiloudi | Stabilization, optimal and robust control. Theory and applications in biological and physical sciences[END_REF] and [START_REF] Lions | Optimal control of systems governed by partial differential equations[END_REF]), the necessary conditions for this optimum is given by the following theorems (see [START_REF] Belmiloudi | Stabilization, optimal and robust control. Theory and applications in biological and physical sciences[END_REF])

Theorem 2.1 If J attains a (local) minimum at a point h * ∈ Y ad , then the following first optimality conditions hold.

J ′ (h * ).(h -h * ) ≥ 0, ∀h ∈ Y ad ,
where J ′ is the directional derivative of J.

In order to solve numerically the optimal control problem it is necessary to derive the gradient of the cost functional J with respect to the control h. For this we suppose that the operator solution F is continuously differentiable on Y ad and its derivative w = F ′ (h).g = lim ǫ-→0

F(h + ǫg) -F(h)
ǫ is the unique solution of the following system (for h = (h 0 , h 1 ) and g = (g 0 , g 1 ))

λ(u) ∂w ∂t - ∂ ∂x (κ(u) ∂w ∂x ) + F (u, w) = 0, in Q, -κ(u(0)) ∂w ∂x (0) = -h 0 -G(u)(0) + 4ϕ 0 |u(0)| 3 w(0) + g 0 (u f -u(0)), in (0, T ), κ(u(L)) ∂w ∂x (L) = -h 1 + G(u)(L) + 4ϕ 1 |u(L)| 3 w(L) + g 1 (u a -u(L)), in (0, T ), w(x, 0) = 0, x ∈ Ω, (2.3)
where operators F and G are defined as follows,

F (u, w) = ∂ ∂t (λ(u))w - ∂ ∂x (w ∂ ∂x (κ(u))), G(u) = ∂ ∂x (κ(u)).
Remark 2.2 We verify easily that

λ(u) ∂w ∂t - ∂ ∂x (κ(u) ∂w ∂x ) + F (u, w) = ∂(λ(u)w) ∂t - ∂ 2 ∂x 2 (κ(u)w), κ(u(0)) ∂w ∂x (0) + G(u)(0)w(0) = ∂(κ(u)w) ∂x (0), κ(u(L)) ∂w ∂x (L) + G(u)(L)w(L) = ∂(κ(u)w) ∂x (L).
In order to derive (2.3), we write the two systems satisfied by

u ǫ = F(h + ǫg) and u = F(h) as λ(u ǫ ) ∂u ǫ ∂t - ∂ ∂x (κ(u ǫ ) ∂u ǫ ∂x ) = 0, in Q, -κ(u ǫ (0)) ∂u ǫ ∂x (0) = (h 0 + ǫg 0 )(u f -u ǫ (0)) +ϕ 0 (|u f | 3 u f -|u ǫ (0)| 3 u ǫ (0)), in (0, T ), κ(u ǫ (L)) ∂u ǫ ∂x (L) = (h 1 + ǫg 1 )(u a -u ǫ (L)) +ϕ 1 (|u a | 3 u a -|u ǫ (L)| 3 u ǫ (L)), in (0, T ), u ǫ (x, 0) = u 0 (x), x ∈ Ω (2.4) and λ(u) ∂u ∂t - ∂ ∂x (κ(u) ∂u ∂x ) = 0, in Q, -κ(u(0)) ∂u ∂x (0) = h 0 (u f -u(0)) + ϕ 0 (|u f | 3 u f -|u(0)| 3 u(0)), in (0, T ), κ(u(L)) ∂u ∂x (L) = h 1 (u a -u(L)) + ϕ 1 (|u a | 3 u a -|u(L)| 3 u(L)), in (0, T ), u(x, 0) = u 0 (x), x ∈ Ω.
(2.5)

By subtracting (2.4) from (2.5), dividing the obtained system by ǫ and letting ǫ tend to zero, the system satisfied by w, called the tangent linear model (TLM) is given by the system (2.3).

We can now show the first-order necessary conditions (optimality conditions) and calculate the gradient of J, by using TLM and by introducing an intermediate co-state model. We introduce the following projection:

Π [a,b] (Ψ) = max(a, min(b, Ψ))
where Ψ is an arbitrary function and a, b are given real constants.

Theorem 2.2 If J attains a (local) minimum at a point h * ∈ Y ad , then

h * 0 = Π [c 0 ,d 0 ] (-u f + u * (0))ũ * (0) α 0 , h * 1 = Π [c 1 ,d 1 ] (-u a + u * (L))ũ * (L) α 1 ,
or in the variational inequality formulation (for all h ∈ Y ad )

T 0 (h 1 -h * 1 ) [(u a -u * (L))ũ * (L) + α 1 h * 1 ] dt + T 0 (h 0 -h * 0 ) [(u f -u * (0))ũ * (0) + α 0 h * 0 ] dt ≥ 0,
where u * = F(h * ) and ũ * = F * (h * ) with the function ũ = F * (h) which is the unique solution of the adjoint (co-state) problem (with initial value given at final time T ) given by

-λ(u) ∂ ũ ∂t - ∂ ∂x (κ(u) ∂ ũ ∂x ) + G(u) ∂ ũ ∂x = C * C(u -u obs ), in Q, -κ(u(0)) ∂ ũ ∂x (0) = -(h 0 + 4ϕ 0 |u(0)| 3 )ũ + γ 0 (u(0) -u 0 obs ), in (0, T ), κ(u(L)) ∂ ũ ∂x (L) = -(h 1 + 4ϕ 1 |u(L)| 3 )ũ + γ 1 (u(L) -u L obs ), in (0, T ), λ(u(., T ))ũ(., T ) = D * D(u(., T ) -u T obs ), in Ω. (2.6)
Moreover, the gradient of J at any element h of Y ad can be given by

∂J ∂h (h) = ((u f -u(0))ũ(0) + α 0 h 0 , (u a -u(L))ũ(L) + α 1 h 1 ).
Proof: see appendix A. We point out that the adjoint problem (2.6), which is backward in time, can be transformed into an initial-boundary value problem by the time transformation t := Tt, which allows to employ [START_REF] Belmiloudi | Mathematical and numerical analysis of dehydratation of gypsum plasterboards exposed to fire[END_REF] for the existence of a unique solution of (2.6) for a sufficiently regular data.

Optimization procedure

By solving successively the direct problem and the adjoint problem, we can therefore calculate the gradient of the objective function relative to the control parameters h = (h 0 , h 1 ). Once the gradient of the objective function J, ∇J, is known, we can seek a minimum of J. For a given observation (u obs , v obs ), the optimization algorithm is summarized in Table .1.

Control parameters

h ----→ Direct Problem u = F(h) ----→ u -u obs , u(., T ) -v obs h:=h-ρ∇J(h)     Convergence test ← ---- Gradient of J ∇J(h) ← ---- Adjoint Problem ũ = F * (h)   convergence Optimal solution (h * , u * = F(h * ))
Table 1: Optimization algorithm: J is minimized until some convergence criteria are attained.

Numerical analysis and simulations

In this section, we seek to estimate separately the parameter h 0 and h 1 . So, in order to facilitate the presentation, we denote by η the control function which plays the role of h i , for i = 0 or 1. Moreover, we consider two kinds of control problem:

• distributed observations control problems where C is the identity function, γ 0 = γ 1 = 0, D = 0 and

J(η) = 1 2 ||u -u obs || 2 L 2 (Q) , (3.1)
• boundary observations control problems with C = 0, γ 0 = γ 1 = 1, D = 0 and

J(η) = 1 2 ||u(0) -u 0 obs || 2 L 2 (0,T ) + 1 2 ||u(L) -u L obs || 2 L 2 (0,T ) , (3.2)
where u is the solution of the direct problem corresponding to the control function η which will be denoted in the sequel by F(η). Then, the expression of the corresponding gradient can be given (for both kinds of cost functions) by

g = ∇J(η) = [u f -u(0, .)] ũ(0, .) if η = h 0 , g = ∇J(η) = [u a -u(L, .)] ũ(L, .) if η = h 1 , (3.3)
where ũ is the solution of adjoint problem (2.6), corresponding to the direct model.

Numerical implementation and outline

As noted in the previous subsection, solving the nonlinear boundary control problem (2.2) by a gradient method requires, at each iteration of the optimization algorithm, to solve the direct problem and its corresponding adjoint problem. In order to solve numerically these two problems, we use first the discretization in space with Lagrange finite elements of order 2 and the derived differential systems are integrated by using variable-order, variable-step-size backward differentiation formulas [START_REF] Brown | Using Krylov methods in the solution of large-scale differential-algebraic systems[END_REF]. Finally the obtained non linear algebraic systems are solved with a Newton's Method [START_REF] Deuflhard | A modified Newton method for the solution of ill-conditioned systems of nonlinear equations with application to multiple shooting[END_REF] and at each iteration, a direct method is used to solve the considered linear system. In order to solve the optimization problems, we have used the Gradient Method (GM) and the Conjugate Gradient Method (CGM). For the CGM method we have considered the following well-known descent direction methods: Fletcher-Reeves, Polak-Ribière, Hestenes-Stiefel and the recent method of Dai and Yuan [START_REF] Dai | A nonlinear conjugate gradient method with a strong global convergence property[END_REF]. More precisely, for the iteration index k, we denote by g (k) the numerical approximation of the gradient function g (given by (3.3)), d (k) the descent direction, ρ k the descent step and η k the numerical approximation of the control function η, at the kth iteration of the algorithm, the considered gradient schemes are as follows:

1. Initialization: η 0 (given) (a) Compute u 0 = F(η 0 ) by solving the direct problem.

(b) Compute ũ0 by solving the adjoint problem.

(c) Gradient of J at η 0 , g (0) is given by (3.3).

(d) Determine the direction:

d (0) = -g (0)
.

(e) Determine η 1 = η 0 + ρ 0 d (0) and initialize k = 1.

2. Compute u k = F(η k ) by solving the direct problem. -------

Fletcher-Reeves CGM CGF R d (k) = -g (k) + β k-1 d (k-1) β k-1 = ||g (k) || 2 ||g (k-1) || 2 Polak-Ribieres CGM CGP R d (k) = -g (k) + β k-1 d (k-1) β k-1 = (g (k) -g (k-1) ; g (k) ) ||g (k-1) || 2 Hestenes-Stiefel CGM CGHS d (k) = -g (k) + β k-1 d (k-1) β k-1 = (g (k) ; g (k) -g (k-1) ) (d (k-1) ; g (k) -g (k-1) ) Dai-Yuan CGM CGDY d (k) = -g (k) + β k-1 d (k-1) β k-1 = ||g (k) || 2 (d (k-1) ; g (k) -g (k-1) )
6. Determine the descent step ρ k by the following methods:

Method Notation Description Constant step cstρ ρ k = ρ Decreasing step dimρ ρ 0 = ρ, ρ k = ρ k-1 2 if J(η k ) > J(η k-1 )(not changed else) Minimum step min ρ k = min(1, 1 ||g (k) || ∞ ) Optimal step opt ρ k is an approximate solution of min ρ≥0 J(η k + ρd (k) ) 7. Determine η k+1 = η k + ρ k d (k)
8. IF the gradient is sufficiently small THEN end; ELSE set k := k + 1 and GOTO 2. The approximation of the optimal solution η * is η k .

Remark 3.1 An approximate value of the optimal step is obtained using the linear approximation of F. For example, in the case of the distributed observations control problem, we have

ρ k ≃ - (u k -u obs , w k ) L 2 (Q) ||w k || L 2 (Q)
where w k is solution of (2.3) for the iteration k.

Remark 3.2 If the exact control function η ex is known, we can measure the efficiency of the method with the following relative error on η k :

errel(η k ) = ||η ex -η k || L 2 (Q) ||η ex || L 2 (Q) . (3.4)
Then, to solve the optimal control problem (2.2), we have developed a specific software, based on Comsol and Matlab tools, taking into account on: first, the nature of the nonlinearity in the operators λ and κ, and the nonlinear radiative term on the boundary of the domain, second, the nature of the adjoint problem which is backward in time and coupled with the direct problem.

In order to validate our approach, we have studied several examples in different situations and by using different descent direction methods for the optimization algorithm. In this paper we present only two application examples and a realistic simulation model. For the first application, we consider the model without the radiative boundary conditions and we assume that the observation u obs is an analytical given function. In the second one, the observation is computed by solving the model with the radiative boundary conditions, corresponding to given convective heat transfer parameters. For the realistic situation, we consider the real parameters and data given in [START_REF] Belmiloudi | Mathematical and numerical analysis of dehydratation of gypsum plasterboards exposed to fire[END_REF][START_REF] Mehaffey | A model for predicting heat transfer though gypsum-board/wood-stud walls exposed to fire[END_REF] and we construct the observation model as in the second application. With these examples, we validate the method and discuss two numerical questions: the improvement of the computed solution at final time and the stability of the algorithm for noisy data.

The numerical simulations are performed on a computer with a processor 2.8 GHz Core 2 Duo and 4 Gb memory and take between 1400 and 2200s for 50 optimization iterations (depending on the descent direction method and on the treated example).

Numerical examples and validation

In this section, we denote by η ex the exact value of the control function η that we want approximate by η k , at the kth iteration of the gradient algorithm. Moreover, we assume T = L = 1, we fix the time step δt and space step δx to 10 -3 .

To simplify the presentation, we detail only the results for the control on the convective heat transfer coefficient h 1 and give some comments for the control on the convective heat transfer coefficient h 0 .

Example without radiation conditions

In this example, the control function η plays the role of h 1 , the observation is the given function u obs = e t+L-x and η ex (= h 1 ) = 1 + cos(t). The other operators and data of the model are given by:

h 0 = 1, ϕ 0 = ϕ 1 = 0, λ(u) = 1 + u, κ(u) = 1 + u 2 , u 0 (x) = e L-x , u f (t) = 2 + e 2(t+L) e t+L , u a (t) = 1 -1+e 2t
1+cos t e t and we have added the following function F (x, t) = e 2(t+L-x) (1 -3e t+L-x ) such that (1.3) becomes

λ(u) ∂u ∂t - ∂ ∂x (κ(u) ∂u ∂x ) = F
to ensure that u obs is the analytical solution of the direct problem (corresponding to η ex ).

Remark 3.3

We want to emphasize on the fact that the addition of the artificial right-hand member F does not change the formulation of the adjoint problem, the expression of the gradient of the objective function, and then the optimization algorithm.

To validate our approach, we have tested all the gradient methods described above. With all these gradient methods we get computed values η k which converge towards the exact value η ex : the relative error (3.4) decreases when k increases (see Figures 1,2,3). But the convergence speed depends on the method and on the initial value. So, we are going to discuss these two points for the distributed control problem (3.1) and then compare with the results of the boundary control problems (3.2).

First, to use the optimization algorithm we have chosen an initial value η 0 for the control function η. We present on Figure 1 the convergence curves, log 10 (errel(η k )) versus k, for different initializations and on Figure 2 the function η k (t) for some k (computed with the CGDYOPT method: Dai-Yuan Conjugate Gradient with optimal step) and for these different initializations: (a) η 0 = 2.2, (b) η 0 = η ex (0), (c) η 0 = 1.7, (d) η 0 = η ex (T ), (e) η 0 = 1, (f) η 0 = [η ex (T )η ex (0)] t T + η ex (0). As expected, the method is convergent for each initialization but the convergence speed depend on the proximity of the initial value η 0 from the exact value η ex but the best accuracy that can be obtained with the method depends on the gap between η 0 (T ) and η ex (T ). Indeed, for each numerical simulation, we note that η k (T ) = η 0 for each k (see Figure 2 for the function η k (t) with (a) η 0 = η ex (0) = 2 and (b) η 0 = η ex (T ) = 1 + cos(T )): the last value of η k is never modified. For this kind of optimization problem, there is not final observation and the gradient (like the descent direction) is always zero for t = T . Then the parameter can not be well approximated near the final time T . Two methods will be proposed in the next subsection to solve this problem. Of course, this has no effect when the initialization is a continuous function with η 0 (T ) = η ex (T ) (see figure 2-(b)). From a physical point of view, the good choice for the initial value can be η 0 such that η 0 (0) = η ex (0) because the parameters are assumed known at ambient temperature (corresponding to t = 0). It is why, we use this initialization for the next examples. To compare the different gradient methods, we fix the initialization to η 0 = η ex (T ) and we present on Figure 3-(a) the convergence curves, log 10 (errel(η k )) versus k, for the CGPR Method with different optimization steps, and on Figure 3-(b) the convergence curves for the different gradient methods with optimal step. We conclude that the choice of the optimal step gives a significant improvment compared to other possibilities. As expected the Conjugate Gradient Methods have a better convergence speed than the gradient method. Among the Conjugate Gradient Methods, the Dai-Yuan method seems the better choice because it gives the better result for almost all iterations. This is in agreement with the literature results in optimization. It is why we present mainly the results for this method for the next examples. Now, we want to compare the distributed observations control problem and the boundary observations control problem. We present on Figure 4 versus k, for the CGDY Method and η 0 = η ex (0). In this case, the convergence is a little bit better for the boundary control problem but it depends on the initialization. The important point is that the results are close and it is sufficient to measure the values of u only in x = 0 and x = L.

Remark 3.4

To explore the difference between the control on the cold side of the domain (x = L) and the control on the hot side of the domain (x = 0), we have tested the control on the convective heat transfer coefficient h 0 (the function η plays here the role of h 0 ). We have considered the same data as in the previous example except for the values of h 0 and h 1 that are exchanged η ex (= h 0 ) = 1 + cos(t), h 1 = 1 and the related data

u f (t) = 1 + 1+e 2(t+L)
1+cos t e t+L , u a (t) = -e 3t , and we obtained similar results to previous. We propose here two methods to overcome the difficulty to compute the solution of the inverse problem at the final time, the gradient of the functional being null for each step of the algorithm. The first possibility is to add a term with final observation to the functional and the second to deduce the value of the solution at final time by a linear approximation of previous values as presented in [START_REF] Yang | The identification of a Robin coefficient by a conjugate gradient method[END_REF]. We present these methods for the distributed control problem.

The first method consists of replacing the expression of J in (3.1) by

J(η) = 1 2 ||u -u obs || 2 L 2 (Q) + β 2 ||u(T ) -u T obs || 2 L 2 (Ω) .
where β is a positive constant. More precisely, D is the identity function multiplied by √ β in 2.2. The expression of the gradient (3.3) is unchanged, but the definition of the adjoint problem (2.6) is modified because of the new value of D.

The choice of the value of β is important. If it is too large the algorithm does not converge to the expected solution for the same initialization because, with only a final observation (problem of exact controllability [START_REF] Lions | Contrlabilit exacte, perturbations et stabilisation de systmes distribus[END_REF][START_REF] Lions | Exact Controllability, Stabilization and Perturbations for Distributed Systems[END_REF]), different values of η give solutions u of the direct problem with the same final value. If it is too small, there is no change with the new functional. As we can see in Figure 5, a good choice consists to give to the final observation the same importance than each of the others observations with for β a value close to the time step δt = 10 -3 . We use here the data of the example 3.2.1 defined in the previous section.

The second method consists of a linear approximation of the solution at final time by previous values of the solution. In [START_REF] Yang | The identification of a Robin coefficient by a conjugate gradient method[END_REF], is proposed the following formula

η k (T ) = 2η k (T -δt) -η k (T -2δt)
which is efficient if the time step is not too small. For δt = 10 -3 , the value of η k (T ) converges very slowly to the exact value, see Figure 6. This convergence can be accelerated with the following generalization of the formula

η k (T ) = 2η k (T -n lc δt) -η k (T -2n lc δt)

Example with radiation conditions

For this second kind of tests, there is no more analytical expression for the observation u obs which is now computed as the solution of the direct problem (1.3)-(1.4)-(1.5) with the exact value η ex of the control function η. Then, solving the control problem from an initialization η 0 of η, we compute approximations η k of η ex using the CGDYOPT method. As in the previous section we detail the results for the control on the convective heat transfer coefficient h 1 and give a remark for the control on the convective heat transfer coefficient h 0 .

With the data h 0 = 1,

h 1 = 1 + cos(t), ϕ 0 = ϕ 1 = 10 -4 , λ(u) = 1 + u, κ(u) = 1 + u 2 , u a = 10, u f (t) = 20 + ln(t), u 0 = u a ,
H 0 = h 0 + ϕ 0 (u f + u obs )(u 2 f + u 2 obs ), H 1 = h 1 + ϕ 1 (u a + u obs )(u 2
a + u 2 obs ). The functions h 0 and H 0 , h 1 and H 1 are compared on Figure 8.

In this example, the control function η plays the role of H 1 with η ex = H 1 and we solve the same optimization problem as in Example 3.2.1 with F = 0 and (h 0 , h 1 ) replaced by (H 0 , H 1 ). The functions u obs and u a being known, we can deduce ϕ 1 from H 1 and h 1 , or h 1 from H 1 and ϕ 1 .

For the initialization η 0 = η ex (0), we present on Figure 9: (a) the convergence curves, log 10 (errel(η k )) versus k, for both kinds of problems (with distributed observations and boundary observations), (b) the function η k (t) for some k (corresponding to the boundary observation). The results are good but contrary to what is expected, we note again a better result with the boundary observations (it depends on the initialization, for example with η 0 = η ex (T ) the method with distributed observations is a little bit better). Remark 3.5 We observed that the result is worse if we initialize our algorithm with η 0 = η ex (T ) (instead of η 0 = η ex (0)). It can be explained by the fact that for η 0 = η ex (T ), η k is not changed for t small because without a source term in the equation (1.3), the value of H 1 , used in the boundary condition at the point x = L, has not a significant effect immediately on u. This comes from the fact that the source term u f is only on the boundary at the point x = 0, and there is a delay depending on the parameters of the equation. We checked that there is no more delay when the source term and the control are at the same point x = 0.

Numerical stability of the algorithm

In real problems, the observation u obs comes from experimental data with measurement errors and it is interesting to study the influence of these errors on the result given by the algorithm.

As in the stability study presented in [START_REF] Yang | The identification of a Robin coefficient by a conjugate gradient method[END_REF] for a linear problem, the noisy data are generated at each point number i of the mesh by the formula ũobs,i = u obs,i (1

+ δζ i ),
where δ gives the relative noise level and ζ i is a random number obtained by a normal distribution of mean zero and unit standard deviation (Matlab function randn).

On the Figure 10, we present the results for the example 3.2.1 with boundary observation and final observation (β = 10 -2 ) (the results are similar with n lc = 2 instead of final observation), initialization η 0 (t) = η ex (0), δt = 10 -3 , δx = 10 -2 . For the 3 noise levels (1%, 3% and 5%), the computed solution remains close to the exact solution with a gap comparable to the noise level. The results are completely different for the example 3.2.3 (see Figure 11). For noise level δ ≥ 0.0005, the computed solution does not remain close to the exact solution. The radiation condition on the boundary, with term (|u a | 3 u a -|u| 3 u), is very sensitive to the perturbations of the data. For the space step δx = 0.1 cm, the time step δt = 1 s and the initialization η 0 = η ex (0), we present on Figure 13: (a) the convergence curves, log 10 (errel(η k )) versus k, for both kinds of problems (with distributed observation and boundary observation), (b) the function η k (t) for some k (corresponding to the boundary observation). Despite the fact that the non linear parameters λ(u) and κ(u) are not regular, we get a good convergence of the method.

Commentary and conclusion

Optimal boundary heating control strategies for time-dependent thermal problems in spatially 1-dimensional domains, related to dehydratation of gypsum plasterboards exposed to fire, are developed. Optimal heating strategies and calibration of process models (parameter identification) are obtained as solutions of certain minimization problems and are computed from conjugate gradient method by considering the following well-known descent direction methods: Fletcher-Reeves, Polak-Ribière, Hestenes-Stiefel and Dai and Yuan. Other choices of control variables can be envisaged. The numerical results show the efficiency of the developed method. This method can be applied without additional argument to estimate emissivities. Indeed, we can introduce φ = (ϕ 0 , ϕ 1 ) and the following optimal control problem: find (u, h, φ) such that the following objective function

J(h, φ) = 1 2 Q |C(u -u obs )| 2 dxdt + 1 2 Ω |D(u(T ) -u T obs )| 2 dx + 1 2 T 0 γ 0 |u(0) -u 0 obs | 2 + γ 1 |u(L) -u L obs | 2 dt + 1 2 T 0 |Mh| 2 dt + 1 2 T 0 |N φ| 2 dt is minimized with respect to (h, φ) ∈ Y ad = V ad × W ad subject to (1.3)-(1.5), (4.1)
where Y ad is the set of admissible controls with

W ad = {φ = (ϕ 0 , ϕ 1 ) ∈ (L 2 (0, T )) 2 such that 0 ≤ a i ≤ ϕ i ≤ b i a.e. in (0, T ), f or i = 0, 1}, V ad = {h = (h 0 , h 1 ) ∈ (L 2 (0, T )) 2 such that 0 ≤ c i ≤ h i ≤ d i a.e. in (0, T ), f or i = 0, 1}, u = F(h, φ), γ 0 and γ 1 are predefined nonnegative weights, M = diag( √ α 0 , √ α 1 ) and N = diag( √ β 0 , √ β
(α i + β i ) = 0. The operators C and D are unbounded operators on L 2 (Ω) satisfying (∀v ∈ H 1 (Ω)) ||Cv|| 2 L ∞ (Ω) ≤ δ 1 ||v|| 2 L 2 (Ω) + δ 2 ||v|| 2 H 1 (Ω) and ||Dv|| 2 L 2 (Ω) ≤ δ 1 ||v|| 2 L 2 (Ω) + δ 2 ||v|| 2 H 1 (Ω) , with δ 1 + δ 2 > 0 and δ i ≥ 0, for i = 1, 2.
With same kind of calculus that previously, we obtain a more general form of theorem 2.2.

Theorem 4.1 If J attains a (local) minimum at a point (h * , φ * ) ∈ Y ad , then

h * 0 = Π [c 0 ,d 0 ] (-u f + u * (0))ũ * (0) α 0 , h * 1 = Π [c 1 ,d 1 ]
(-u a + u * (L))ũ * (L) It would be interesting to use the here developed method with observations coming from experimental data and a more complete description of the dehydration process including reactions for high temperatures . It is clear that, to get even closer to a more realistic calculation, it is necessary to study in the future the mathematical analysis and numerical validation of the controllability of the system, such as presented, for example, for linear parabolic problems in [START_REF] Respondek | Numerical analysis of controllability of a parabolic system with the delayed controls and non-zero boundary conditions[END_REF] and the references therein. In order to take into account the influence of noises and fluctuations, we can also study robust control problems by using the approach developed in Belmiloudi's Book [START_REF] Belmiloudi | Stabilization, optimal and robust control. Theory and applications in biological and physical sciences[END_REF].

A Proof of the theorem 2.2 By using the same technique as in [START_REF] Belmiloudi | Stabilization, optimal and robust control. Theory and applications in biological and physical sciences[END_REF], we start by calculating the variation of J.

According to the regularity of the operator solution F and the nature of the cost function J (J is composition of differentiable mappings), we have that J is differentiable and the directional G-derivative of J at point h along the direction g can be given by J ′ (h).g =

Remark 2 . 1

 21 The general form (2.2) of the objective function allows to consider different kinds of problem (we denote I the identity and d a positive scalar): 1. distributed observations control problems with C = dI, γ 0 = γ 1 = 0 and D = 0, 2. boundary observations control problems with C = 0, γ 0 = γ 1 = d and D = 0, 3. final observations control problems with C = 0, γ 0 = γ 1 = 0 and D = dI, 4. single control function, for example h 1 with α 1 = 0 and α 0 = 0.

3 . 4 .

 34 Compute ũk by solving the adjoint problem. Gradient of J at η k , g (k) is given by (3.3). 5. Determine the direction d (k) by one of the following expressions (where (.; .) and ||.|| are the scalar product and its associated norm): ) = -g(k)

Figure 1 :

 1 Figure 1: Example 3.2.1 -log 10 (errel(η k )) versus k for different initializations with the CGDY-OPT method

Figure 2 :

 2 Figure 2: Example 3.2.1 -η k (t) for some k with different initializations: (a) η 0 = η ex (0) = 2, (b) η 0 = η ex (T ) = 1 + cos(T ).

Figure 3 :

 3 Figure 3: Example 3.2.1 -log 10 (errel(η k )) versus k with (a) CGPR method and different optimization steps, (b) an optimal step and different gradient methods.

Figure 4 :

 4 Figure 4: Example 3.2.1 -log 10 (errel(η k )) versus k for the distributed and boundary control problems

Figure 5 :Figure 6 :

 56 Figure 5: Example 3.2.1 with Final observation: (a) log 10 (errel(η k )) versus k for different β, (b) η k (t) for some k with β = 0.001.

Figure 7 :

 7 Figure 7: Example 3.2.1 -Linear approximation (a) log 10 (errel(η k )) versus k for different n lc (b) η k (t) for some k with n lc = 2.

Figure 9 :

 9 Figure 9: Example 3.2.3 with n lc = 2 -(a) log 10 (errel(η k )) versus k, (b) η k (t) for some k with boundary observation.

Figure 10 :

 10 Figure 10: Example 3.2.1 with β = 0.01 and boundary observation -(a) log 10 (errel(η k )) versus k, (b) η 20 (t), for different values of noise level δ.

Figure 13 :

 13 Figure 13: Simulation on real data with n lc = 2: (a) log 10 (errel(η k )) versus k, (b) η k (t) for some k with boundary observation.

1 , 0 (h 1 0 (ϕ 1 0 (h 0

 1010100 [a 0 ,b 0 ] (-|u f | 3 u f + |u * (0)| 3 u * (0))ũ * (0) β 0 , ϕ * 1 = Π [a 1 ,b 1 ] (-|u a | 3 u a + |u * (L)| 3 u * (L))ũ * (L) β or in the variational inequality formulation (for all (h, φ) ∈ Y ad ) T h * 1 )((u au * (L))ũ * (L) + α 1 h * 1 )dt + T ϕ * 1 )((|u a | 3 u a -|u * (L)| 3 u * (L))ũ * (L) + β 1 ϕ * 1 )dt + T h * 0 )((u fu * (0))ũ * (0) + α 0 h * 0 )dt + T 0 (ϕ 0ϕ * 0 )((|u f | 3 u f -|u * (0)| 3 u * (0))ũ * (0) + β 0 ϕ * 0 )dt ≥ 0,where u * = F(h * , φ * ) and ũ * = F * (h * , φ * ) with the function ũ = F * (h, φ) which is the unique solution of the adjoint (co-state) problem (with initial value given at final time T ) given by C(u -u obs ), in Q, -κ(u(0)) ∂ ũ ∂x (0) = -(h 0 + 4ϕ 0 |u(0)| 3 )ũ + γ 0 (u(0)u 0 obs ), in (0, T ), κ(u(L)) ∂ ũ ∂x (L) = -(h 1 + 4ϕ 1 |u(L)| 3 )ũ + γ 1 (u(L)u L obs ), in (0, T ), λ(u(., T ))ũ(., T ) = D * D(u(., T )u T obs ), in Ω. Moreover, the gradient of J at any element (h, φ) of Y ad can be given by∂J ∂h (h, φ) = ((u fu(0))ũ(0) + α 0 h 0 , (u au(L))ũ(L) + α 1 h 1 ) ∂J ∂φ(h, φ) = ((|u f | 3 u f -|u(0)| 3 u(0))ũ(0) + β 0 ϕ 0 , (|u a | 3 u a -|u(L)| 3 u(L))ũ(L) + β 1 ϕ 1 ).

QC 0 γ 0 0 (M 2 h 1 + 0 [g 1 0 h 0 0 γ 1 0 [g 1 0 [g 0 (A. 3 ) 0 [g 1 0 g 1 [

 00021010001010030101 * C(uu obs )w dxdt + Ω D * D(u(T )u T obs )w(T ) dx + T (u(0)u 0 obs )w(0) + γ 1 (u(L)u L obs )w(L) dt + T h)g dt.Now, we simplify the directional derivative of J. For this we multiply the first part of (2.3) by some regular function ũ, using Green's formula and integrating by parts in times, we obtain(since w(., 0) = 0) (T ))ũ(T )w(T ) dx = 0, (A.1)According to the second and third parts of (2.3) (boundary conditions), we can deduce thatQ 4ϕ 1 |u(L)| 3 ũ(L) + κ(u(L)) ∂ ũ ∂x (L) w(L) dt -T (u au(L))] ũ(L) dt + T + 4ϕ 0 |u(0)| 3 ũ(0)κ(u((u fu(0))] ũ(0)dt + Ω λ(u(T ))ũ(T )w(T ) dx = 0. (A.2)In order to simplify the problem (A.2), we assume that ũ satisfies the "adjoint" (or dual) problem (2.6) (with initial value given at final time T ).Then the problem (A.2) becomesQ C * C(uu obs )w dxdt + Ω D * D u(T )u T obs w(T ) dx + T (u(L)u L obs )w(L) + γ 0 (u(0)u 0 obs )w(0) dt -T (u au(L))] ũ(L) dt -T (u fu(0))] ũ(0) dt = 0.According to the previous result, we can deduce thatJ ′ (h).g = T (u au(L))] ũ(L) dt + T 0 [g 0 (u fu(0))] ũ(0) dt + T 0 (M 2 h) g dtand then J ′ (h).g = T (u au(L))ũ(L) + α 1 h 1 ] dt + T 0 g 0 [(u fu(0))ũ(0) + α 0 h 0 ] dt.

  the convergence curves, log 10 (errel(η k ))

		2.1												2.2										η 0
		2												2.1										η ex η 1
																								η 2
		1.9												2										η 5
																								η 10
	η k (t)	1.7 1.8			η 1 η 0 η ex								η k (t)	1.8 1.9
					η 2
		1.6			η 5 η 10									1.7
		1.5			η 20 η 50									1.6
		1.4	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1	1.5	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1
								t												t

  Example 3.2.3 -Comparing data (a) h 0 and H 0 , (b) h 1 and H 1 .

				3													2.5
																								H 1
			2.8													2.4						h 1
			2.6									H 0				2.3
			2.4									h 0				2.2
	and H 0	2 2.2											and H 1	2 2.1
	h 0	1.8											h 1	1.9
			1.6													1.8
			1.4
																	1.7
			1.2
																	1.6
				1
				0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1			1.5	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1
									t															t
									(a)														(b)
	log 10 (errel(η k ))	-2.8 -2.6 -2.2 -2 -1.8 -1.6 -1.4 -1.2 -2.4	5 Figure 8: 0 10 15 CGDYOPT -n lc 20 25 =2 -η 0 30 = η ex (0) 35 distributed observation 40 45 boundary observation	50	η k (t)	2.05 2.1 2.2 2.25 2.3 2.35 2.4 2.45 2.15	0		0.1	0.2 η ex η 0 η 1 η 2 η 5 η 10	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1
									k															t

  (2.2) (with the cost function and its gradient given by (3.1) or (3.2) and (3.3)) by the optimization algorithm (using the CGDYOPT method) and compute error estimations on the control function η at each optimization iteration.

								η 0 =η ex (0)
		2.45
		2.4
		2.35
		2.3
	η 7 (t)	2.25
		2.2			η ex
		2.15			δ = 0 δ = 0.0001
					δ = 0.0005
		2.1			δ = 0.001
					δ = 0.005
		2.05	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1
								t

Figure 11: Example 3.2.3 with β = 0.01 and boundary observation -η 7 (t) for different values of noise level δ.

problem

Simulation on real data

We now apply our method to a realistic problem. To that purpose, we take similary data as in [START_REF] Belmiloudi | Mathematical and numerical analysis of dehydratation of gypsum plasterboards exposed to fire[END_REF][START_REF] Mehaffey | A model for predicting heat transfer though gypsum-board/wood-stud walls exposed to fire[END_REF]. We consider a gypsum plasterboard, with thickness L = 4 cm, exposed to fire in x = 0, during T = 1200s. The surroundings temperature is u a = 20 • C and the furnace temperature is

) (see Figure 12). The initial temperature is u 0 = u a . The thermal conductivity is given by We have λ(u) = ρ(u)C p (u) where the density is defined by

and the specific heat is defined with a peak between 100 and 170 • C, and a value corresponding to the total amount of energy of the two dehydratation reactions, 669kJ:

else .

We supposed that the convective heat transfer coefficient h 0 = 25 and the effective emissivity ǫ 0 = 0.82 are known in x = L. We deduce ϕ 0 = σǫ 0 with σ = 5.6697 10 -8 W/m 2 K 4 the Stefan-Bolzmann's constant. We want to estimate h 1 , supposing ǫ 1 = 0. The control function η plays here the role of h 1 . We choose the following exact expression for h 1 : η ex (= h 1 ) = 5 + cos(t/T ). Then we can compute the observation u obs as the solution of the direct problem (1.3), solve the control Consequently the gradient of J can be given by

Since h * is an optimal solution of J, then according to Theorem 2.1, we have that

By using a standard control argument concerning the sign of the variations (h 0 , h 1 ) (depending, respectively, on the size of (h * 0 , h * 1 )), we obtain that

(-u a + u * (L))ũ * (L) α 1 .

This completes the proof.