Identification and design of source term in a two-region heat conduction problem

Abstract

This work deals with the use of the conjugate gradient method with adjoint problem for the estimation of an unknown source term in a two-region problem. The physical problem consists of heat conduction in two contacting rectangular regions, where the source term is known to exist in only one of them. The source term is supposed to vary in time as well as within the region. Results are presented for the inverse identification problem by using simulated temperature measurements containing random errors. In addition, results are presented for an inverse design problem with application in hot embossing microreplication microfabrication technology employed for the fabrication of polymer-based Micro-electromechanical Systems (MEMS).

Keywords:

• Function estimation
• Conjugate gradient method
• Source term
• Hot embossing
• Micro-electromechanical Systems
• Design

1. Introduction

Inverse problems dealing with the identification of source terms have been drawing the attention of different research groups for more than a decade, because of its importance in many practical applications, such as predicting and controlling chemical reactions and phase-change, detecting inclusions or non-homogeneities, etc. For the solution of this kind of inverse problem, a variety of methods have been examined 1–20.

This study is concerned with the estimation of the source term in a two-region heat conduction problem. The physical problem of interest is typical of the manufacturing of Micro-electromechanical Systems (MEMS) and is called Hot Embossing Microfabrication Microreplication (HEMM). Such technique consists of pressing a mold on a substrate under a prescribed thermal history and can also be used for the fabrication of micro- and macro-scale polymer parts. In general, the polymer substrate and the mold are heated to a temperature just above the polymer glass transition temperature. At this temperature, the polymer undergoes significant variations in its physical properties, notably in the Young's modulus, and changes from hard and brittle to soft and pliable. The mold is then pressed on the polymer until features get replicated and then both the mold and the substrate are cooled down to a temperature below the glass transition temperature for de-embossing 21.

The estimation of the source term within the mold by using the temperature variation at selected points within the substrate is under observation in this article. Such an inverse problem can be aimed at the identification of the unknown source term by using temperature measurements within the substrate, or, alternatively, at the design of the source term that will result in a prescribed temperature history required for the HEMM technique. Results are presented for both the identification and the design types of inverse problems. For the estimation of the source term, no information is a priori assumed available regarding its functional form, except for the functional space that it belongs to. It is assumed that the unknown belongs to the Hilbert space of square integrable functions in the spatial and time domains of interest. The unknown source term is estimated by using the conjugate gradient method with adjoint problem formulation 22–24, as described below.

2. Physical problem and mathematical formulation

The physical problem considered in this work consists of the two-dimensional heat conduction in two contacting rectangular regions, as depicted in figure 1. The widths of both regions are a, while the thickness of region A is b and the thickness of region B is (c − b). Regions A and B are supposed to initially be at the uniform temperatures T_0A and T_0B, respectively. For t > 0, the two regions are put in contact and heat is generated in region B at a volumetric rate g(x,y,t). The contact conductance at the interface between the two regions, h_c, is supposed to be uniform and the other boundaries of regions A and B are supposed to be insulated. By assuming that the thermal properties of regions A and B are constant, the mathematical formulation for this problem is given in dimensionless form as:

Figure 1. Geometry and coordinates.

Region A (1) (2) (3) (4) (5) (6) Region B (7) (8) (9) (10) (11) (12) In order to write equations (1a–f) and (2a–f), the following dimensionless groups were defined: (13) (14) where k is the thermal conductivity, α is the thermal diffusivity, and the subscripts A and B refer to regions A and B, respectively.

3. Direct problem and inverse problem

The problem defined by equations (1a–f) and (2a–f), with known initial and boundary conditions, contact conductance, thermophysical properties and source term constitutes a direct problem that is concerned with the determination of the transient temperature fields θ_A(X,Y,τ) and θ_B(X,Y,τ).

For the inverse problem of interest here, the source term G(X,Y,τ) is regarded as unknown, while the other quantities appearing in the formulation of the direct problem are considered to be known with a high degree of accuracy. For the solution of the inverse problem, we consider available the transient temperature histories μ_Am(τ) at positions (X_m,Y_m), m = 1, … , M, in region A, as well as μ_Bn(τ), n = 1, … , N, at positions (X_n,Y_n) in region B. For the design inverse problem, μ_Am(τ) and μ_Bn(τ) represent desired temperatures for the fabrication procedure, while for the identification inverse problem, μ_Am(τ) and μ_Bn(τ) represent temperature measurements. The measurements contain errors, which are supposed to be additive, uncorrelated, normally distributed, with known and constant standard deviation and zero mean.

For the estimation of the unknown function G(X,Y,τ), we make no a priori assumption regarding its functional form, except that it belongs to the Hilbert space of square-integrable functions 22–24 in the domain 0 < X < 1, B < Y < C and 0 < τ< τ_f, where τ_f is the duration of the time interval of concern for the inverse analysis. For the solution of the present inverse problem, we consider the minimization of the following functional: (15) The minimization of the objective functional (4) is performed by the conjugate gradient method with adjoint problem formulation 22–24. The iterative procedure of such a method is presented below.

4. Conjugate gradient method

The iterative procedure of the conjugate gradient method, as applied to the estimation of the function G(X,Y,τ), is given by: (16) where the superscript k denotes the number of iterations and β^k is the search step size. The direction of descent, d^k, is obtained as a linear combination of the gradient direction at iteration k with directions of descent at previous iterations. It is given as 22–24: (17) Different expressions are available in the literature for the conjugation coefficient, γ^k. In this work, we use the so-called Fletcher–Reeves version of the conjugate gradient method, where the conjugation coefficient is written as 22–24: (18) For the implementation of the iterative procedure of the conjugate gradient method given by equations (5a–c), expressions are required for the gradient direction ∇ S [G ^k (X, Y, τ)], as well as for the search step size β^k. Two auxiliary problems are used to derive these expressions, namely the sensitivity and adjoint problems, as described next.

5. Sensitivity problem and search step size

The sensitivity problem is used to determine the variations Δθ_A(X,Y,τ) and Δθ_B(X,Y,τ) that the temperatures θ_A(X,Y,τ) and θ_B(X,Y,τ) undergo, respectively, when the unknown function is perturbed by ΔG(X,Y,τ) 22–24. The sensitivity problem is derived by substituting into the direct problem given by equations (1a–f) and (2a–f), θ_A(X,Y,τ) by [θ_A(X,Y,τ) + Δθ_A(X,Y,τ)], θ_B(X,Y,τ) by [θ_B(X,Y,τ) + Δθ_B(X,Y,τ)] and G(X,Y,τ) by [G(X,Y,τ) + ΔG(X,Y,τ)]. The original direct problem is then subtracted from the resulting equations in order to obtain the following sensitivity problem:

Region A (19) (20) (21) (22) (23) (24) Region B (25) (26) (27) (28) (29) (30)

The search step size is obtained by minimizing the objective functional with respect to β^k at each iteration 22–24. The following expression results: (31) where Δθ_A(X,Y,τ,d^k) and Δθ_B(X,Y,τ,d^k) are the solutions of the sensitivity problem given by equations (6a–f) and (7a–f), obtained by setting ΔG(X,Y,τ) = d^k(X,Y,τ).

6. Adjoint problem and gradient equation

The adjoint problem is derived by multiplying the governing equations of the direct problem, equations (1a) and (2a), by the Lagrange multipliers λ_A(X,Y,τ) and λ_B(X,Y,τ), respectively. The equations are then integrated in the spatial and time domains that they are valid and added to the original functional (4). The directional derivative of the extended functional in the direction of the perturbation of the unknown function is then obtained and the resultant expression, after some lengthy but straightforward manipulations, is allowed to go to zero 22–24. The following adjoint problem for the determination of the Lagrange multipliers λ_A(X,Y,τ) and λ_B(X,Y,τ) results:

Region A (32) (33) (34) (35) (36) (37) Region B (38) (39) (40) (41) (42) (43)

By applying the limiting process used to obtain the adjoint problem, the directional derivative of the objective functional along the direction of the perturbation ΔG(X,Y,τ) reduces to: (44)

We now invoke the hypothesis that the unknown function belongs to the Hilbert space of square integrable functions in the domain 0 < X < 1, B < Y < C, and 0 < τ < τ_f, so that we can write such directional derivative as: (45)

Therefore, by comparing equations (11a) and (11b) we obtain the gradient direction as: (46)

After developing the expressions for the search step size and for the gradient direction, the iterative procedure of the conjugate gradient method given by equations (5a–c) can be applied until a suitable convergence criterion is satisfied. In this article, the iterative procedure was stopped when the objective functional given by equation (4) became sufficiently small. The discrepancy principle 22–24 was used to specify the tolerance for the stopping criterion when simulated temperature measurements were used for the inverse analysis of source term identification. For the inverse design problem, the tolerance was specified as a small number.

7. Results and discussions

7.1. Identification problem

Let us consider first, the solution of the inverse problem for the identification of the source term function G(X,Y,τ) by using simulated temperature measurements containing random errors. For the results presented below, the dimensions of regions A and B were taken as a = 0.3 m, b = 0.015 m, and (b − c) = 0.1 m, so that B = 0.05 and C = 0.38. In order to examine strict cases involving materials with thermal conductivity and thermal diffusivity of different orders of magnitudes, we assumed regions A and B to be made of Teflon and steel, respectively, with the following thermophysical properties: k_A = 0.23 Wm⁻¹ K⁻¹, α_A = 1.005 × 10⁻⁷ m²s⁻¹, k_B = 45 Wm⁻¹ K⁻¹ and α_B = 1.240 × 10⁻⁵ m²s⁻¹, that is, K = 195.65 and Λ = 123.35. Both regions A and B were supposed to initially be at the same temperature T_0A = T_0B = 27°C and in perfect thermal contact. The duration of the experiment was taken as 2340 s, with a measurement frequency of 1.1 Hz. In dimensionless terms, the duration of the experiment was τ_f = 2.613 × 10⁻³. However, in order to overcome the difficulties encountered in the inverse problem solution in the neighborhood of the final time, which result from the null gradient at τ = τ_f, see equations (9f), (10f) and (12), the results are presented below up to τ = 2 × 10⁻³. The sensors were supposed to be located only within region A, which constitutes a difficult situation for the estimation of the source term in region B. For the cases examined below, the sensors were evenly distributed at the vertical position Y = 0.03, which corresponds to 0.006 m below the interface between regions A and B. The simulated measurements contained random errors with standard deviation of 0.003, which is equivalent to 1°C. Such simulated measurements were generated by using functions containing discontinuities and sharp-corners for the source term. These functions were supposed to vary from zero to a maximum value of 10⁶ Wm⁻³. The initial guess for the iterative procedure of the conjugate gradient method was taken as a null source term.

The direct, sensitivity, and adjoint problems were solved numerically by using the finite-volume method. Regions A and B were discretized with 30 volumes along the X direction and 25 and 50 volumes along their thickness, respectively. A step of Δτ = 10⁻⁶ was used to march the solution in time. Such numbers of volumes and time step were selected though a grid convergence analysis and by comparing the numerical and analytical solutions for the direct problem.

For the test-cases examined below in the inverse identification problem, the source-function was assumed to be separable in the form: (47)

Table 1 summarizes the test-cases examined for the inverse identification problem. Test-cases 1–3 deal with the identification of a spatially uniform function that only varies in time. These test-cases differ by the number of sensors used in the inverse analysis. For test-case 4, the unknown function is uniform along the Y direction and constant in time, thus varying only along the longitudinal X direction. Test-case 5 involves a function that varies in time and along the longitudinal X direction, but is uniform along the Y direction.

Table 1. Test-cases examined for the identification inverse problem

Test-case	Unknown function	Number of sensors
1	G(X, Y, τ) = Gτ (τ)	5
2	G(X, Y, τ) = Gτ (τ)	10
3	G(X, Y, τ) = Gτ (τ)	1740
4	G(X, Y, τ) = GX (X)	10
5	G(X, Y, τ) = GX (X)Gτ (τ)	10

Figure 2(a)–(d) presents the results for test-case 1, involving the estimation of the time variation of a spatially uniform source term, obtained by using five sensors. Figure 2(a)–(d) present the results obtained for the positions Y = 0.05, 0.16, 0.27, and 0.38, respectively, and for each of these figures the results for X = 0.02, 0.22, 0.42, and 0.98 are given. We notice in figure 2(a)–(d) that the accuracies of the estimated functions deteriorate for vertical positions closer to the interface at Y = 0.05, or closer to the upper boundary at Y = 0.38. On the other hand, the estimated functions are affected very little by the longitudinal position X, except for Y = 0.05. At Y = 0.05, the estimated functions are less accurate in regions near the surfaces at X = 0 and X = 1. It is interesting to notice in figure 2(a)–(d) that the differences between estimated and exact functions, in the regions where the exact function is underestimated, are practically identical to those differences in the regions where the exact function is overestimated. Therefore, the total energy input is accurately estimated, despite the fact that the source term functional form is not exactly recovered for some cases.

Figure 2. Test-case 1–inverse problem solution obtained with five sensors for a uniform exact function containing discontinuities in time: (a) Y = 0.05, (b) Y = 0.16, (c) Y = 0.27, and (d) Y = 0.38.

Figure 3(a)–(d) presents the results obtained for test-case 2, involving the same exact function of figure 2(a)–(d), but the measurements of 10 sensors instead of 5. A comparison of figures 2 and 3 shows that, as expected, the solution improves as the measurements of more sensors are used, since more information is available for the inverse analysis. Ideally, if a large number of sensors could be used, including sensors in the mold, the agreement between the exact and the estimated functions would be excellent even for regions distant from the interface, as illustrated in figure 4 (test-case 3). This figure shows the results for Y = 0.38, obtained with one sensor in each control-volume in the substrate and one sensor in two-thirds of the control-volumes of the mold.

Figure 3. Test-case 2–inverse problem solution obtained with 10 sensors for a uniform exact function containing discontinuities in time: (a) Y = 0.05, (b) Y = 0.16, (c) Y = 0.27, and (d) Y = 0.38.

Figure 4. Test-case 3–inverse problem solution at Y = 0.38, obtained with 1740 sensors for a uniform exact function containing discontinuities in time.

We now examine a test case involving the identification of a source term that varies along the X direction, but is constant in time and uniform along the Y direction. Figure 5(a)–(c) presents the results obtained for test-case 4, involving the estimation of a source term containing sharp corners for τ = 0.0001, τ = 0.001, and τ = 0.002, respectively, obtained by using 10 sensors in the inverse analysis. In each of these figures, the results for Y = 0.05, 0.16, 0.27, and 0.38 are given. These figures show that the accuracy of the inverse problem solution deteriorates for the positions located farther from the interface, because the measurements are less affected by the source term in such locations. However, as for the results shown in figures 2–4, the total energy input was accurately estimated. The inverse problem solution was not affected by the time variable in this case, even for τ = 0.002, which is near the final time considered for the problem.

Figure 5. Test-case 4–inverse problem solution obtained with 10 sensors for an exact function varying along the longitudinal direction: (a) τ = 0.0001, (b) τ = 0.001, and (c) τ = 0.002.

The results for test-case 5, dealing with a source term that varies in time and along the longitudinal X direction, but is uniform along the Y direction, are presented in figure 6(a)–(c) for τ = 0.0001, τ = 0.001 and τ = 0.002, respectively. Such results were obtained with the measurements of 10 sensors. As for the case involving only the variation of the source term along the X direction (figure 5(a)–(c), we notice in figure 6(a)–(c) that the results for vertical positions close to the interface are more accurate. Furthermore, the results for this case become less accurate as the final time is approached (figure 6(c), for τ = 0.002), when the estimated function is affected by the null gradient. We also examined test-cases involving the identification of source terms varying along the vertical Y direction. However, the exact function could not be recovered unless measurements were considered available within region B as well, which was an expected result from the analysis of Savateev 25 regarding the solvability of an inverse identification problem similar to the present one.

Figure 6. Test-case 5–inverse problem solution obtained with 10 sensors for an exact function varying in time and along the longitudinal direction: (a) τ = 0.0001, (b) τ = 0.001, and (c) τ = 0.002.

7.2. Design problem

We now consider the inverse problem of designing the heat source term in the mold so that the temperature at the mold/substrate interface follows a desired history imposed by the manufacturing process. Based on data of reference 21, we assumed the substrate and the mold to have the following properties: k_A = 0.19 Wm⁻¹ K⁻¹, α_A = 1.144 × 10⁻⁷ m² s⁻¹, k_B = 14.8 W m K⁻¹ and α_B = 9.056 × 10⁻⁶ m² s⁻¹, that is, K = 77.89 and Λ = 79.16. Both regions A and B were supposed to be initially at the same temperature T_0A = T_0B = 27°C. The mold and substrate dimensions (figure 1) were supposed to be a = 0.075 m, b = 0.001 m, and (b − c) = 0.065 m and the duration of the hot embossing process was taken as 1400 s, that is, B = 0.013, C = 1, and τ_f = 0.028. The regions were discretized with the same number of volumes used for the inverse identification problem examined above, but the time increment was taken as Δτ = 10⁻⁵. Uniform interface temperature and perfect contact between regions A and B were assumed for the results presented hereafter.

Figure 7 presents the desired temperature variation at the interface between regions A and B, as well the calculated temperatures obtained with the heat source term that resulted from the solution of the present inverse design problem. The agreement between desired and calculated temperatures is perfect within the graph scale, despite the fact that the desired temperature variation contains sharp corners. Figure 8(a)–(d) presents the resultant heat source term at positions Y = 0.013, 0.34, 0.67, and 1.0, respectively. Each of these figures shows the results for X = 0.02, 0.22, 0.42, and 0.98. We note in figure 8(a)–(d) that the source term resultant from the design procedure does not vary along the X direction for such test case. In addition, periods of heating (positive source) and cooling (negative source) can be clearly identified, as required for the sharp increase and decrease in temperature at times τ = 0 and τ = 0.02, respectively (figure 7). The obtained results reveal that the conjugate gradient method of function estimation is capable of providing smooth and stable solutions for the source term. The magnitude of the source term decreases for positions farther from the interface. This is due to the low sensitivity of the interface temperature with respect to the source term at these locations, similarly to the inverse identification problem examined above. We note that the present inverse design analysis assumes that the source term is continuous within mold; but for practical implementation, discrete heat sources and cooling channels will be required. The results obtained with the conjugate gradient method of function estimation permit the selection of the number and position of the heat sources and cooling channels, as well as the heating power, cooling fluid temperature, and convective heat transfer coefficient, which will result in the desired interface temperature variation.

Figure 7. Comparison between desired and calculated interface temperatures for the inverse design problem.

Figure 8. Inverse design solution: (a) Y = 0.013, (b) Y = 0.34, (c) Y = 0.67, and (d) Y = 1.

8. Conclusions

This work dealt with the solution of the inverse problem of estimating the source term function in a rectangular region, by using as input data the temperature history at points in another contacting region. The conjugate gradient method of function estimation was applied for the solution of the inverse problem. Several test cases were examined in the study, including the identification of functions varying only in time, only in space and varying simultaneously in space and time, by using simulated temperature measurements in the inverse analysis. In addition, results were also presented for the design of the source term for a prescribed interface temperature required for a manufacturing process.

Generally, for the inverse identification problem the agreement between estimated and exact functions improves for those locations closer to the interface between the two regions. For the inverse design problem, the conjugate gradient method is capable of recovering smooth functions that result in perfect agreement between desired and calculated interface temperatures.

Nomenclature

a,b,c	=	Dimensions of the rectangular region
d	=	Direction of descent
G(X,Y,τ)	=	Dimensionless source term
S	=	Objective functional
X,Y	=	Dimensionless Cartesian coordinates

Greeks

Δθ(X,Y,τ)	=	Sensitivity function, solution of the sensitivity problem given by equations (6) and (7)
β	=	Search step size
γ	=	Conjugation coefficient
λ(X,Y,τ)	=	Lagrange multiplier, solution of the adjoint problem given by equations (9) and (10)
μ(τ)	=	Measured or desired temperature
θ(X,Y,τ)	=	Dimensionless calculated temperature, solution of the direct problem given by equations (1) and (2)
τ	=	Dimensionless time
τf	=	Dimensionless final time

Subscripts

A,B

=

Regions A and B, respectively

Superscripts

k

=

Number of iterations

Acknowledgments

The Brazilian authors are thankful for the support provided by CNPq and FAPERJ, agencies for the fostering of science from the Brazilian and State of Rio de Janeiro governments, respectively. H. R. B. O and M. J. C are thankful to the hospitality of Prof. G. Dulikravich and of the Mechanical & Aerospace Engineering Department during their visits to UT-Arlington.