Research Article | Open Access
On Nonlinear Inverse Problems of Heat Transfer with Radiation Boundary Conditions: Application to Dehydration of Gypsum Plasterboards Exposed to Fire
The paper investigates boundary optimal controls and parameter estimates to the well-posedness nonlinear model of dehydration 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 . Numerical simulations illustrate several numerical optimization methods, 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.
1. Introduction and Mathematical Setting of the Problem
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 appreciating 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 dehydration process consumes large amount of heat) and provide the plasterboard excellent fire resistance.
The main particularity of gypsum plasterboard is that it contains of chemically bound water by weight. When the temperature reaches 100°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.
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 implement the experiment. In addition pilot-scale testing does not allow to accurately reproduce 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 appear as a suitable means to study the thermal behaviour of the material during a fire exposure.
1.2. Thermochemistry 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 ) and approximately per weight of chemically combined water of crystallization (see, e.g., [1–3]). Its chemical formula is CaSO4·2H2O (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°C when the calcium sulphate di-hydrate is converted to calcium sulphate hemihydrate (the reaction is always complete by 160°C) as shown by the following reaction: The amount of energy required by this first dehydration is about 500 kJ per kg of gypsum (see ).
The second dehydration reaction occurs when calcium sulphate hemihydrate is converted to calcium sulphate anhydrate as shown by the following reaction: The amount of energy corresponding to this second reaction is about 169 kJ per kg (see ).
Remark 1. Other amounts of energy can be found in the literature for these reactions; see  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) and (2) in the specific heat evolution. The first dehydration reaction occurs at approximately 100–160°C; on the other hand, there is some controversy to when the second dehydration reaction occurs. Andersson et al.  (e.g.,) estimate that the second reaction occurs between 210°C and 300°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 . In the numerical examples, we will choose this model with only one peak between 100 and 170°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., [6, 9–11]).
Remark 2. The model of dehydration can be completed by other reactions at high temperatures. For example, a third endothermic reaction, for decomposition of CaCO3, is introduced at 650°C in  or 800°C in . In , the authors note an exothermic reaction around 400°C.
1.3. 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 dehydration 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.
1.4. 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 that 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 (). Therefore, heat fluxes in lateral () direction and vertical () direction can be neglected in front of the heat flux in direction . 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 -direction. The one-way heat transfer through a plane wall is described by the heat equation in its one-dimensional form where is the boundary subset , is the temperature, , is the density of body material, is the specific heat, and is the thermal conductivity. The functions and are variables, positive, 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 in , as where, for all , is the convective heat transfer coefficient between the furnace and the plasterboard surface, is the furnace temperature, , is Stefan-Boltzmann's constant, and 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 in , as where, for all , is convective heat transfer coefficient between the surroundings and the plasterboard surface, is the surroundings temperature, , and is effective emissivity of the surface.
Remark 3. In the sequel, to simplify the notations, we denote for , for , for , for , for , and for .
We assume that(H1)functions and and satisfy the following pointwise constraint:
for some positive constants ,
(H2)operators , are sufficiently regular with and with for any .
Remark 4. 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).
In the physical case there are not absolute values under the boundary conditions (since the temperature is nonnegative). 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.
It is clear that we can obtain in the same way the model in -dimensional for as follows: where is an open bounded domain with boundary , such that , and is the outward unit normal vector on . Boundaries and denote the cold side and the fire side, respectively, and denotes the other surface of the plasterboard. The well-posedness of problem (7) in 3-dimensional can be obtained in similar way as in .
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., [3–5, 7, 11, 15]), have assumed the convective heat transfer coefficients and and the relative emissivity as constants or/and neglected the relative emissivity 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 depend on temperature and large uncertainties exist in regard to the quality of the data reported. Moreover, the work of Belmiloudi and Le Meur  shows very clearly that the radiative heat transfer between the unexposed surface and the surrounding cannot 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) depends 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.
2. The Inverse Problem Formulation
2.1. Problem Formulation
We assume that there exists a unique solution of problem (3), (6) with boundary conditions (4), (5) under some hypotheses for the data and some regularity of the operators and , satisfying the following regularity (by using [16, 17]):
Introduce now the mapping which maps the source term of (3)–(5) into the corresponding solution , where . In this section we formulate the inverse problem as an optimal control problem. The control procedure consists of finding the optimal controls and the corresponding optimal temperature which minimize a cost criterion . The cost functional measures the distance between a measured temperature (the observations , , , and ) and the corresponding predicted temperature obtained from the primal (or direct) model (3)–(5). Precisely we will study the following optimal control problem .
Find such that the following objective function:
is minimized with respect to subject to (3)–(5), where is the set of admissible controls, , and are predefined nonnegative weights, and are predefined nonnegative weights such that . The operators and are unbounded operators on satisfying (): with and , for .
Remark 5. The general form (9) of the objective function allows considering different kinds of problems (we denote as the identity and as a positive scalar): (1)distributed observations control problems with , , and ,(2)boundary observations control problems with , , and ,(3)final observations control problems with , , and ,(4)single control function, for example, , with and .
We consider here an optimal control problem with Tichonov regularization terms . The nonnegative weights and play the role of regularization parameters and take small values in the numerical simulations.
According to (6) the set of admissible controls describing the constraint is
2.2. First-Order Necessary Conditions
Assume that the nonlinear control problem (9) admits an optimal solution (for similar result see, e.g., [19, 20]); the necessary conditions for this optimum is given by the following theorems (see ).
Theorem 6. If attains a (local) minimum at a point , then the following first optimality conditions hold: where is the directional derivative of .
In order to solve numerically the optimal control problem it is necessary to derive the gradient of the cost functional with respect to the control . For this we suppose that the operator solution is continuously differentiable on and its derivative is the unique solution of the following system (for and ):
where operators and are defined as follows:
Remark 7. We verify easily that
In order to derive (13), we write the two systems satisfied by and as
We can now show the first-order necessary conditions (optimality conditions) and calculate the gradient of by using TLM and by introducing an intermediate costate model. We introduce the following projection: where is an arbitrary function and are given real constants.
Theorem 8. If J attains a (local) minimum at a point , then
or in the variational inequality formulation (for all )
where and with the function which is the unique solution of the adjoint (costate) problem (with initial value given at final time ) given by
Moreover, the gradient of at any element of can be given by
Proof. See the Appendix.
We point out that the adjoint problem (21), which is backward in time, can be transformed into an initial-boundary value problem by the time transformation , which allows to employ  for the existence of a unique solution of (21) for a sufficiently regular data.
2.3. 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 . Once the gradient of the objective function , , is known, we can seek a minimum of . For a given observation , the optimization algorithm is summarized in Figure 1.
3. Numerical Analysis and Simulations
In this section, we seek to estimate separately the parameters and . So, in order to facilitate the presentation, we denote by the control function which plays the role of , for or . Moreover, we consider two kinds of control problems:(i)distributed observations control problems, where is the identity function, , , and (ii)boundary observations control problems with , , , and where is the solution of the direct problem corresponding to the control function which will be denoted in the sequel by . Then, the expression of the corresponding gradient can be given (for both kinds of cost functions) by where is the solution of adjoint problem (21), corresponding to the direct model.
3.1. Numerical Implementation and Outline
As noted in the previous subsection, solving the nonlinear boundary control problem (9) by a gradient method requires, at each iteration of the optimization algorithm, solving 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 and the derived differential systems are integrated by using variable-order, variable-step-size backward differentiation formulas . Finally, the obtained nonlinear algebraic systems are solved with Newton's method  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 . More precisely, for the iteration index , we denote by the numerical approximation of the gradient function (given by (25)), the descent direction, the descent step, and the numerical approximation of the control function ; at the th iteration of the algorithm, the considered gradient schemes are as follows.(1)Initialization: (given):(a)compute by solving the direct problem;(b)compute by solving the adjoint problem;(c)gradient of at , is given by (25);(d)determine the direction: ;(e)determine and initialize .(2)Compute by solving the direct problem.(3)Compute by solving the adjoint problem.(4)Gradient of at , is given by (25).(5)Determine the direction by one of the expressions shown in Table 1 (where and are the scalar product and its associated norm).(6)Determine the descent step by the methods shown in Table 2. (7)Determine .(8)If the gradient is sufficiently small then end; else set and go to . The approximation of the optimal solution is .
Remark 9. An approximate value of the optimal step is obtained using the linear approximation of . For example, in the case of the distributed observations control problem, we have where is solution of (13) for the iteration .
Remark 10. If the exact control function is known, we can measure the efficiency of the method with the following relative error on :
Then, to solve the optimal control problem (9), we have developed a specific software, based on Comsol and Matlab tools, taking into account the following: first the nature of the nonlinearity in the operators and and the nonlinear radiative term on the boundary of the domain and 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 used 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 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 [4, 16] 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 taken between 1400 and 2200 s for 50 optimization iterations (depending on the descent direction method and on the treated example).
3.2. Numerical Examples and Validation
In this section, we denote by the exact value of the control function that we want to approximate by , at the th iteration of the gradient algorithm. Moreover, we assume that and we fix the time step and space step to .
To simplify the presentation, we detail only the results for the control on the convective heat transfer coefficient and give some comments for the control on the convective heat transfer coefficient .
3.2.1. Example without Radiation Conditions
In this example, the control function plays the role of ; the observation is the given functions and . The other operators and data of the model are given by , , and we have added the following function such that (3) becomes to ensure that is the analytical solution of the direct problem (corresponding to ).
Remark 11. We want to emphasize on the fact that the addition of the artificial right-hand member does not change the formulation of the adjoint problem, the expression of the gradient of the objective function, and 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 which converge towards the exact value : the relative error (27) decreases when increases (see Figures 2, 3, and 4). 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 (23) and then compare with the results of the boundary control problems (24).
First, to use the optimization algorithm we have chosen an initial value for the control function . We present on Figure 2 the convergence curves, versus , for different initializations and on Figure 3 the function for some (computed with the CGDYOPT method: Dai-Yuan conjugate gradient with optimal step) and for these different initializations: (a) , (b) , (c) , (d) , (e) , and (f) . As expected, the method is convergent for each initialization but the convergence speed depends on the proximity of the initial value from the exact value but the best accuracy that can be obtained with the method depends on the gap between and .
Indeed, for each numerical simulation, we note that for each (see Figure 3 for the function with (a) and (b) ): the last value of is never modified. For this kind of optimization problem, there is no final observation and the gradient (like the descent direction) is always zero for . Then, the parameter cannot be well approximated near the final time . 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 (see Figure 3(b)). From a physical point of view, the good choice for the initial value can be such that because the parameters are assumed known at ambient temperature (corresponding to ). It is why we use this initialization for the next examples.
To compare the different gradient methods, we fix the initialization to and we present on Figure 4(a) the convergence curves, versus , for the CGPR method with different optimization steps and on Figure 4(b) the convergence curves for the different gradient methods with optimal step. We conclude that the choice of the optimal step gives a significant improvement 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 to be 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 5 the convergence curves, versus , for the CGDY method and . 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 only in and .
Remark 12. To explore the difference between the control on the cold side of the domain () and the control on the hot side of the domain (), we have tested the control on the convective heat transfer coefficient (the function plays here the role of ). We have considered the same data as in the previous example except for the values of and that are exchanged , and the related data , and we obtained similar results to the previous ones.
3.2.2. How to Improve the Result at Final Time?
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 . We present these methods for the distributed control problem.
The first method consists of replacing the expression of in (23) by where is a positive constant. More precisely, is the identity function multiplied by in (9). The expression of the gradient (25) is unchanged, but the definition of the adjoint problem (21) is modified because of the new value of .
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 [25, 26]), different values of give solutions 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 6, a good choice consists to give the final observation the same importance than each of the other observations with being a value close to the time step . We use here the data of the example of Section 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 , is proposed the following formula: which is efficient if the time step is not too small. For , the value of converges very slowly to the exact value; see Figure 7. This convergence can be accelerated with the following generalization of the formula: with a small integer which should not be too large to avoid instability. We present on Figure 8(a) the convergence curves for different values of and on Figure 8(b) computed with . We will use this value in the sequel.
3.2.3. Example with Radiation Conditions
For this second kind of tests, there is no more analytical expression for the observation which is now computed as the solution of the direct problem (3)-(4)-(5) with the exact value of the control function . Then, solving the control problem from an initialization of , we compute approximations of using the CGDYOPT method.
As in the previous section we detail the results for the control on the convective heat transfer coefficient and give a remark for the control on the convective heat transfer coefficient .
The observation being known on the boundary and because of the time dependence of and , we can consider an optimization problem without radiation condition on the boundary. Indeed, the boundary conditions of the direct problem become with
The functions and , and are compared in Figure 9.
In this example, the control function plays the role of with and we solve the same optimization problem as in Example of Section 3.2.1 with and replaced by . The functions and being known, we can deduce from and or from and .
For the initialization , we present on Figure 10(a) the convergence curves, versus , for both kinds of problems (with distributed observations and boundary observations) and (b) the function for some (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, e.g., with the method with distributed observations is a little bit better).
Remark 13. We observed that the result is worse if we initialize our algorithm with (instead of ). It can be explained by the fact that for , is not changed for small because without a source term in (3), the value of , used in the boundary condition at the point , has no significant effect immediately on . This comes from the fact that the source term is only on the boundary at the point , 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 .
3.2.4. Numerical Stability of the Algorithm
In real problems, the observation 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  for a linear problem, the noisy data are generated at each point number of the mesh by the formula where gives the relative noise level and is a random number obtained by a normal distribution of mean zero and unit standard deviation (Matlab function randn).
On Figure 11, we present the results for the example of Section 3.2.1 with boundary observation and final observation () (the results are similar with instead of final observation) and initialization , . For the 3 noise levels (, , and ), 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 of Section 3.2.3 (see Figure 12). For noise level , the computed solution does not remain close to the exact solution. The radiation condition on the boundary, with term , is very sensitive to the perturbations of the data.
3.3. Simulation on Real Data
We now apply our method to a realistic problem. To that purpose, we take similarly data as in [4, 16]. We consider a gypsum plasterboard, with thickness cm, exposed to fire in , during s. The surroundings temperature is C and the furnace temperature is (see Figure 13). The initial temperature is . The thermal conductivity is given by
We have , where the density is defined by and the specific heat is defined with a peak between and and a value corresponding to the total amount of energy of the two dehydration reactions, 669 kJ:
We supposed that the convective heat transfer coefficient and the effective emissivity are known in . We deduce that with , the Stefan-Boltzmann's constant.
We want to estimate , supposing that . The control function plays here the role of . We choose the following exact expression for : . Then, we can compute the observation as the solution of the direct problem (3), solve the control problem (9) (with the cost function and its gradient given by (23) or (24) and (25)) by the optimization algorithm (using the CGDYOPT method), and compute error estimations on the control function at each optimization iteration.
For the space step cm, the time step s, and the initialization , we present in Figure 14(a) the convergence curves, versus , for both kinds of problems (with distributed observation and boundary observation) and (b) the function for some (corresponding to the boundary observation). Despite the fact that the nonlinear parameters and are not regular, we get a good convergence of the method.
4. Commentary and Conclusion
Optimal boundary heating control strategies for time-dependent thermal problems in spatially 1-dimensional domains, related to dehydration 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 and the following optimal control problem:
Find such that the following objective function:
is minimized with respect to subject to (3)–(5), where is the set of admissible controls with where , and are predefined nonnegative weights, , and are predefined nonnegative weights such that . The operators and are unbounded operators on satisfying ():