Inverse heat transfer analysis for design and control of a micro-heater array

Abstract

A non-linear inverse heat source identification problem is described and solved. The inverse problem analysis is used in the design of an embedded micro-heater array and to estimate the required control settings, which are the input currents to each heating element, to generate as close as possible to a prescribed temperature profile on the surface of a thin copper film. The purpose of the micro-heater array is to control the local copper microstructure through control of the local temperature field. A finite element model of the micro-heater system is used to define a discrete set of non-linear equations used as a basis for the inverse problem solution. Two methods are explored to solve the inverse problem, a direct minimization method with Tikhonov regularization and a passivity-based feedback control algorithm. A uniform and a linear temperature distribution could be attained in the central region above the micro-heater array, but the temperatures near the edges of the domain could not be controlled due to heat loss at the edges. Thus, to control the temperature field over the full width of the domain, the heater array must extend beyond the domain of interest. Both methods to solve the inverse problem are found to perform well. The regularization method allows for a smoother solution, while the feedback control method is simpler as the coefficient matrix for which the update remains unchanged for each iteration.

Keywords:

• Inverse heat source problem
• micro-heater array
• regularization
• feedback control

AMS Subject Classifications:

• 80A23

1. Introduction

Defining the thermal–mechanical processing conditions to obtain desired microstructural characteristics and resulting material properties is a major priority for materials manufacturers.[1,2] In thermal processing of metals, energy, in the form of dislocations and grain boundaries, is released through recovery, recrystallization and grain growth leading to a change in the microstructure.[3] The resulting grain size distribution is one of the most important characteristics of the microstructure because it directly affects the yield strength, ductility and fatigue resistance.[4–6] Thus, improvements in temperature control during thermal processing can have a large impact on material performance. The temperature field is normally controlled through manipulating the heat sources, and the heat source parameters are typically obtained from an experience-based trial and error approach, leading to a long material development cycle. In order to better understand the relationship between the local temperature field and microstructure and to provide a means for locally controlling the microstructure through controlling the local temperature field, a micro-heater array has been proposed.[7] This paper presents and compares two methods to define the input parameters to the micro-heater array to generate a prescribed temperature profile on the surface of a thin copper film, which represents a first step in ultimately controlling the microstructure during thermal processing. Specifically, this can be thought of as a feed-forward control problem, where the goal is to find the required current history that should be applied to each heater line in a heater array to generate a desired temperature profile history that will drive the copper microstructure to, as close as possible, some desired microstructure. The simulation tools developed are also used for refining the design of the micro-heater array.

The problem of finding the heat source in each heating element in an embedded micro-heater array to obtain as close as possible to a prescribed surface temperature distribution represents an inverse problem. The problem is non-linear because of the temperature-dependent material properties and the presence of a radiation boundary condition. Several researchers have also considered the inverse heat source identification problem for a variety of applications. Much of that work focuses on identifying the location and strength of point sources based on measurements of surface temperatures.[8,9] Other researchers have considered spatially dependent and time-dependent sources. A boundary element formulation for solving the linear, two-dimensional steady-state problem with both unknown boundary conditions and heat sources combined with a singular value decomposition to handle the resulting ill-conditioned system was proposed by Martin and Dulikravich [10]. In Yan et al. [11], a spatially dependent heat source is found given interior temperature measurements using a finite difference approximation of the governing equations and applying a singular value decomposition with regularization to solve the inverse problem. The problem of determining a time-dependent heat source that leads to a prescribed temperature field at a given time has been investigated by Jiang et al.[12] using a modified conjugate gradient method. Chen et al. [13] use an extended Kalman filter and weighted recursive least-squares estimator to solve for a time-varying heat source. In Oden et al. [14], an adjoint finite element method is used to determine the best combination of heating parameters to generate a prescribed temperature field. Frackowiak et al. [15] studied the inverse heat source problem to find interior heat sources from measured temperatures on the surface of a turbine blade. They used an iterative optimization process based on polyharmonic functions to minimize the error functional on the temperature at the outer boundary with respect to the heat source control. An overview of the more general class of inverse source problems that discusses mathematical issues of uniqueness, existence and stability of solutions can be found in Isakov [16].

In this work, two approaches to solve the optimal heater settings in an embedded micro-heater array to drive the surface temperature profile towards a prescribed target temperature field are described and compared. The two approaches investigated are: a direct inverse method using Zeroth-order Tikhonov regularization and a passivity-based feedback control algorithm. The micro-heater array is a linear array that generates a two-dimensional temperature field through the volume of the device with a one-dimensional temperature profile at the surface. Furthermore, as will be shown, a steady-state temperature profile is obtained in about 10 ms, as compared to the grain growth experimental time frame of minutes and hours. Thus, the time-dependent feed-forward control problem of finding the current history that should be applied to each heater to achieve a temperature profile history can be reduced to a sequence of steady-state inverse problems at discrete times. The primary innovation in the work presented here is the application of the methods to the design of a novel micro-heater array, and the comparison of the direct inverse method to the feedback control approach.

The paper is organized as follows. A description of the proposed micro-heater array design is presented first together with the forward transient model and finite element formulation. The inverse problem is then defined. The two methods for solving the inverse problem are formulated next. Then, results from two examples are presented and the performance of the two methods is discussed. Implications on the heater array design are also considered. Finally, conclusions are made.

Figure 1. Schematic of the heater array model (a) cross section and (b) top view showing titanium resistive heater lines.

2. Methods

2.1. Model description

A micro-heater array that can be operated in a scanning electron microscope (SEM) is proposed for local temperature control in order to study and control local grain growth.[7] Micro-fabricated heater arrays to generate local temperature gradients, similar to that described here, have been demonstrated previously.[17,18] A schematic of the initial design is shown in Figure 1. The device is to be fabricated on a standard $500\mathrm{\mu }\mathrm{m}$ thick silicon (Si) substrate. Ten titanium (Ti) resistive heating strips will be used to generate a temperature profile along x over a distance of 1 mm in a $1\mathrm{\mu }\mathrm{m}$ thick copper (Cu) film. The heating strips will be relatively long, 2.7 mm, compared to the width, and thus, the temperature field in the central section of the heater array will be essentially two-dimensional in the x–y plane with little variation in the z direction except near the ends of the lines. The titanium heater strips will be embedded in silicon dioxide (${\text{SiO}}_2$) to insulate the heaters from the Si substrate and Cu film. The width w and thickness h of the Ti heating elements are design parameters, and the thickness of the ${\text{SiO}}_2$, d, and position of the Ti heaters in the ${\text{SiO}}_2$, denoted b, must also be prescribed. When current is delivered to the Ti heaters through aluminium (Al) lines, the volumetric resistive heat generation rate Q in the Ti is(1) $$\begin{array}{c}\hfill Q=\mathit{\rho }{j}^2\mathrm{i}n{\mathrm{\Omega }}_4\end{array}$$(1)

where $\mathit{\rho }(T)$ is the temperature-dependent electrical resistivity of the Ti, $j(\mathit{x})$ is the current density, which can be different in each heater strip, and ${\mathrm{\Omega }}_4$ is the region occupied by the heater strips.

The usual forward problem is to find the resulting temperature field $T(\mathit{x},t)$, which depends on position $\mathit{x}\in \mathrm{\Omega }$ where $\mathrm{\Omega }$ is the union of all regions in the domain, ${\mathrm{\Omega }}_1$–${\mathrm{\Omega }}_4$, at time $t\in \mathit{\tau }$ given prescribed current densities j in each of the Ti heater lines that satisfies the following system,(2) $$\begin{array}{cc}& \mathrm{\nabla }·(k\mathrm{\nabla }T)+Q=C\frac{\mathit{\partial }T}{\mathit{\partial }t}\mathrm{i}n\mathrm{\Omega }\times \mathit{\tau }\hfill \end{array}$$(2) (3) $$\begin{array}{cc}& T(\mathit{x},0)=T_0(\mathit{x})\text{in}\mathrm{\Omega }\hfill \end{array}$$(3) (4) $$\begin{array}{cc}& T(\mathit{x},t)=\check{T}(\mathit{x},t)\text{on}{\mathrm{\Gamma }}_1\hfill \end{array}$$(4) (5) $$\begin{array}{cc}& \mathit{q}(\mathit{x},t)·\mathit{n}(\mathit{x})=\mathit{\epsilon }\mathit{\sigma }\left(T{(\mathit{x},t)}^4-T_{\infty }^4\right)\text{on}{\mathrm{\Gamma }}_2\hfill \end{array}$$(5) (6) $$\begin{array}{cc}& \mathit{q}(\mathit{x},t)·\mathit{n}(\mathit{x})=0\text{on}{\mathrm{\Gamma }}_3\hfill \end{array}$$(6)

where $k(\mathit{x},T)$ and $C(\mathit{x},T)$ are the material and temperature-dependent thermal conductivity and volumetric heat capacity, respectively, $\mathit{q}=-k\mathrm{\nabla }T$ is the heat flux assuming Fourier’s Law for heat conduction, $\mathit{\epsilon }(T)$ is the emissivity on the Cu surface, $\mathit{\sigma }$ is the Stefan–Boltzmann constant, $\mathit{n}$ is the outward unit normal on the boundary $\mathrm{\Gamma }$, and the temperature-dependent heat source Q(T) is as defined in Equation (1(1) $$\begin{array}{c}\hfill Q=\mathit{\rho }{j}^2\mathrm{i}n{\mathrm{\Omega }}_4\end{array}$$(1) ). Equation (2(2) $$\begin{array}{cc}& \mathrm{\nabla }·(k\mathrm{\nabla }T)+Q=C\frac{\mathit{\partial }T}{\mathit{\partial }t}\mathrm{i}n\mathrm{\Omega }\times \mathit{\tau }\hfill \end{array}$$(2) ) represents conservation of energy. The initial condition on the temperature field is represented by Equation (3(3) $$\begin{array}{cc}& T(\mathit{x},0)=T_0(\mathit{x})\text{in}\mathrm{\Omega }\hfill \end{array}$$(3) ), and the boundary conditions are given in Equations (4(4) $$\begin{array}{cc}& T(\mathit{x},t)=\check{T}(\mathit{x},t)\text{on}{\mathrm{\Gamma }}_1\hfill \end{array}$$(4) )–(6(6) $$\begin{array}{cc}& \mathit{q}(\mathit{x},t)·\mathit{n}(\mathit{x})=0\text{on}{\mathrm{\Gamma }}_3\hfill \end{array}$$(6) ). The temperature is assumed prescribed on the bottom boundary ${\mathrm{\Gamma }}_1$ and a radiation boundary condition on the top boundary ${\mathrm{\Gamma }}_2$, while the side boundaries, ${\mathrm{\Gamma }}_3$ are treated as insulated. We note that since the micro-heater array is designed to operate in the vacuum environment of an SEM, convection is not considered. The simulation domain was made wide enough in the x direction, so that the boundary conditions on ${\mathrm{\Gamma }}_3$ did not have a significant effect.

The temperature dependency of the material parameters k, C, $\mathit{\epsilon }$ and $\mathit{\rho }$ are fit to quadratic functions from experimental data in the literature.[19,20] The temperature dependency of the parameters for each relevant material are as follows. The thermal conductivity k ($\frac{\text{W}}{\text{m}·\text{K}}$) is(7a) $$\begin{array}{cc}\hfill k(T)& =405.7-0.05699T-(9.085\times {10}^{-6})T^2\text{in}{\mathrm{\Omega }}_1\text{(Cu)}\hfill \\ \hfill k(T)& =24.24-0.06094T+(4.800\times {10}^{-5})T^2\text{in}{\mathrm{\Omega }}_2\text{(}{\text{SiO}}_2\text{)}\hfill \\ \hfill k(T)& =295.9-0.5435T+(2.723\times {10}^{-4})T^2\text{in}{\mathrm{\Omega }}_3\text{(Si)}\hfill \\ \hfill k(T)& =22.72-0.01653T+(1.375\times {10}^{-5})T^2\text{in}{\mathrm{\Omega }}_4\text{(Ti)}\hfill \end{array}$$(7a) The volumetric heat capacity C ($\frac{\text{J}}{{\text{m}}^3·\text{K}}$) is(8a) $$\begin{array}{cc}\hfill C(T)& =3.741\times {10}^6-1544T+1.615T^2\text{in}{\mathrm{\Omega }}_1\text{(Cu)}\hfill \\ \hfill C(T)& =1.034\times {10}^6+3721T-1.159T^2\text{in}{\mathrm{\Omega }}_2\text{(}{\text{SiO}}_2\text{)}\hfill \\ \hfill C(T)& =1.209\times {10}^6+2318T-1.497T^2\text{in}{\mathrm{\Omega }}_3\text{(Si)}\hfill \\ \hfill C(T)& =3.212\times {10}^6-2809T+3.218T^2\text{in}{\mathrm{\Omega }}_4\text{(Ti)}\hfill \end{array}$$(8a) Copper emissivity $\mathit{\epsilon }$ is(9) $$\begin{array}{c}\hfill \mathit{\epsilon }(T)=0.07681+0.0003696T-(1.932\times {10}^{-7})T^2\text{on}{\mathrm{\Gamma }}_2\end{array}$$(9)

Titanium electrical resistivity $\mathit{\rho }$ ($\mathrm{\Omega }·$m) is(10) $$\begin{array}{c}\hfill \mathit{\rho }(T)=-1.893\times {10}^{-7}+(2.462\times {10}^{-9})T-(6.674\times {10}^{-13})T^2\text{in}{\mathrm{\Omega }}_4\end{array}$$(10)

The temperature T is in K in each of the relations given above.

2.2. Forward problem formulation – finite element model

A standard finite element formulation is used to solve the non-linear forward problem described in Equations (1(1) $$\begin{array}{c}\hfill Q=\mathit{\rho }{j}^2\mathrm{i}n{\mathrm{\Omega }}_4\end{array}$$(1) )–(6(6) $$\begin{array}{cc}& \mathit{q}(\mathit{x},t)·\mathit{n}(\mathit{x})=0\text{on}{\mathrm{\Gamma }}_3\hfill \end{array}$$(6) ), which is summarized here for completeness. The temperature field is approximated in terms of finite element basis functions at a given time t such that(11) $$\begin{array}{c}\hfill T(\mathit{x},t)\approx \widehat{T}(t,\mathit{x})={\overline{T}}_{\mathit{\alpha }}(t){\mathit{\psi }}_{\mathit{\alpha }}(\mathit{x}),\mathit{\alpha }=1,N_T\end{array}$$(11)

where $\widehat{T}$ is a finite dimensional approximation of T, ${\mathit{\psi }}_{\mathit{\alpha }}$ are finite element basis functions, and $N_T$ is the number of temperature degrees of freedom (DOF). Greek subscripts refer to finite element basis function numbers, and overbar denotes coefficients of interpolating functions, which are the values of the solution at the finite element nodes. Summation is implied on the repeated index $\mathit{\alpha }$. Substituting into the standard Galerkin formulation yields the following equation for the temperature field $\widehat{T}$(12) $$\begin{array}{c}\hfill \underset{\mathrm{\Omega }}{\int }\left(C\frac{\mathit{\partial }\widehat{T}}{\mathit{\partial }t}\widehat{w}+k\mathrm{\nabla }\widehat{T}·\mathrm{\nabla }\widehat{w}-Q\widehat{w}\right)\mathrm{d}\mathrm{\Omega }+\underset{{\mathrm{\Gamma }}_2}{\int }\mathit{\epsilon }\mathit{\sigma }\left({\widehat{T}}^4-T_{\infty }^4\right)\widehat{w}\mathrm{d}\mathrm{\Gamma }=0\end{array}$$(12)

where $\widehat{w}$ is a finite dimensional weighting function, which is interpolated with the same functions as the temperature. It should also be emphasized that k, C, Q and $\mathit{\epsilon }$ depend on the solution $\widehat{T}$. After substituting in the finite element interpolations (11(11) $$\begin{array}{c}\hfill T(\mathit{x},t)\approx \widehat{T}(t,\mathit{x})={\overline{T}}_{\mathit{\alpha }}(t){\mathit{\psi }}_{\mathit{\alpha }}(\mathit{x}),\mathit{\alpha }=1,N_T\end{array}$$(11) ), a non-linear system of equations results of the form(13) $$\begin{array}{c}\hfill \mathit{a}(\dot{\overline{\mathit{T}}},\overline{\mathit{T}})=\mathbf{0}\end{array}$$(13)

where $\overline{\mathit{T}}$ represents the vector of nodal temperatures ${\overline{T}}_{\mathit{\alpha }}$, $\mathit{\alpha }=1,N_T$, and $\dot{\overline{\mathit{T}}}$ represents its time derivative. This system is linear in $\dot{\overline{\mathit{T}}}$ and non-linear in $\overline{\mathit{T}}$. The system is integrated in time using a backward difference approximation for the time derivative, and then the resulting non-linear system is solved with a standard Newton–Raphson method at each time. The procedure continues in time until a steady-state temperature is reached.

In simulations of typical micro-heater array designs, it is found that the temperature reaches steady-state very fast relative to the rate at which the microstructure evolves because the domain is very small and thus has a very low thermal inertia. Specifically, steady-state is reached in about 10 ms while grain growth occurs on the scale of minutes and hours. Thus, designing for steady-state (quasi-static) is valid.

2.3. Inverse problem formulation

The inverse problem of interest here is to find the heat source Q that should be applied to each heater strip to generate as close as possible to a desired temperature distribution $T^d$ on the top surface ${\mathrm{\Gamma }}_2$. Because of the low thermal inertia, we treat the system as quasi-static. The equations to solve for T and Q are then(14) $$\begin{array}{cc}& \mathrm{\nabla }·(k\mathrm{\nabla }T)+Q=0\mathrm{i}n\mathrm{\Omega }\hfill \end{array}$$(14) (15) $$\begin{array}{cc}& T(\mathit{x})=\check{T}(\mathit{x})\text{on}{\mathrm{\Gamma }}_1\hfill \end{array}$$(15) (16) $$\begin{array}{cc}& \mathit{q}(\mathit{x})·\mathit{n}(\mathit{x})=0\text{on}{\mathrm{\Gamma }}_3\hfill \end{array}$$(16) (17) $$\begin{array}{cc}& \mathit{q}(\mathit{x})·\mathit{n}(\mathit{x})=\mathit{\epsilon }\mathit{\sigma }\left[{(T(\mathit{x}))}^4-T_{\infty }^4\right]\text{on}{\mathrm{\Gamma }}_2\hfill \end{array}$$(17) (18) $$\begin{array}{cc}& T(\mathit{x})\approx T^d(\mathit{x})\text{on}{\mathrm{\Gamma }}_2\hfill \end{array}$$(18)

We can see from this that we now have over-specified boundary conditions on ${\mathrm{\Gamma }}_2$ with both the temperature $T^d$ and effectively the flux $\check{q}=\mathit{q}·\mathit{n}\approx \mathit{\epsilon }\mathit{\sigma }({T^d}^4-T_{\infty }^4)$ prescribed approximately. We say approximately because it may not be possible to find a set of heat sources Q that leads to exactly the desired temperature $T=T^d$ on ${\mathrm{\Gamma }}_2$, and thus, we do not strictly enforce that condition. Also note that once the optimal heat source in each heater strip Q is determined, since the temperature will also be known from the solution of the governing equations, the current density j that should be applied in each heater strip may then be backed out from Equation (1(1) $$\begin{array}{c}\hfill Q=\mathit{\rho }{j}^2\mathrm{i}n{\mathrm{\Omega }}_4\end{array}$$(1) ).

Since the heat source is only applied in each of the 10 heater strips and is uniform in each strip, we can write(19) $$\begin{array}{cc}& Q(\mathit{x})={\overline{Q}}_i{\mathit{\varphi }}_i(\mathit{x}),i=1,N_H\hfill \end{array}$$(19) (20) $$\begin{array}{cc}& {\mathit{\varphi }}_i(\mathit{x})=\left\{\begin{array}{cc}1\hfill & \text{if}\mathit{x}\in {\mathrm{\Omega }}_{4i},\hfill \\ 0\hfill & \text{otherwise}.\hfill \end{array}\right.\hfill \end{array}$$(20)

where ${\overline{Q}}_i$ is the heat source in heater strip i, $N_H=10$ is the number of heater strips, and ${\mathrm{\Omega }}_{4i}$ is the region occupied by heater strip i. Equation (12(12) $$\begin{array}{c}\hfill \underset{\mathrm{\Omega }}{\int }\left(C\frac{\mathit{\partial }\widehat{T}}{\mathit{\partial }t}\widehat{w}+k\mathrm{\nabla }\widehat{T}·\mathrm{\nabla }\widehat{w}-Q\widehat{w}\right)\mathrm{d}\mathrm{\Omega }+\underset{{\mathrm{\Gamma }}_2}{\int }\mathit{\epsilon }\mathit{\sigma }\left({\widehat{T}}^4-T_{\infty }^4\right)\widehat{w}\mathrm{d}\mathrm{\Gamma }=0\end{array}$$(12) ) can now be written as(21) $$\begin{array}{c}\hfill \underset{\mathrm{\Omega }}{\int }\left(k\mathrm{\nabla }\widehat{T}·\mathrm{\nabla }\widehat{w}-{\overline{Q}}_i{\mathit{\varphi }}_i\widehat{w}\right)\mathrm{d}\mathrm{\Omega }+\underset{{\mathrm{\Gamma }}_2}{\int }\check{q}\widehat{w}\mathrm{d}\mathrm{\Gamma }=0\end{array}$$(21)

or symbolically, after substituting in finite element interpolations (11(11) $$\begin{array}{c}\hfill T(\mathit{x},t)\approx \widehat{T}(t,\mathit{x})={\overline{T}}_{\mathit{\alpha }}(t){\mathit{\psi }}_{\mathit{\alpha }}(\mathit{x}),\mathit{\alpha }=1,N_T\end{array}$$(11) ), a system of the following form results(22) $$\begin{array}{c}\hfill \mathit{b}(\overline{\mathit{T}},\overline{\mathit{Q}})=\mathbf{0}\end{array}$$(22)

which is linear in $\overline{\mathit{Q}}$ and non-linear in $\overline{\mathit{T}}$ because k and $\check{q}$ depends on T.

The inverse problem can be expressed as the following minimization problem(23) $$\begin{array}{c}\hfill \underset{\overline{\mathit{Q}}}{min}\frac{1}{2}{\left(\mathit{R}\overline{\mathit{T}}(\overline{\mathit{Q}})-{\overline{\mathit{T}}}^d\right)}^T\left(\mathit{R}\overline{\mathit{T}}(\overline{\mathit{Q}})-{\overline{\mathit{T}}}^d\right)\end{array}$$(23)

where $\mathit{R}$ is a matrix that extracts the components of $\overline{\mathit{T}}$ that are at nodes on ${\mathrm{\Gamma }}_2$, ${\overline{\mathit{T}}}^d$ are the prescribed values of $T^d$ at the nodes on ${\mathrm{\Gamma }}_2$, and the dependency of $\overline{\mathit{T}}$ on $\overline{\mathit{Q}}$ is defined through Equations (21(21) $$\begin{array}{c}\hfill \underset{\mathrm{\Omega }}{\int }\left(k\mathrm{\nabla }\widehat{T}·\mathrm{\nabla }\widehat{w}-{\overline{Q}}_i{\mathit{\varphi }}_i\widehat{w}\right)\mathrm{d}\mathrm{\Omega }+\underset{{\mathrm{\Gamma }}_2}{\int }\check{q}\widehat{w}\mathrm{d}\mathrm{\Gamma }=0\end{array}$$(21) ) and (22(22) $$\begin{array}{c}\hfill \mathit{b}(\overline{\mathit{T}},\overline{\mathit{Q}})=\mathbf{0}\end{array}$$(22) ).

Because of the non-linear relationship between $\overline{\mathit{T}}$ and $\overline{\mathit{Q}}$, an iterative solution is required. We linearize Equation (21(21) $$\begin{array}{c}\hfill \underset{\mathrm{\Omega }}{\int }\left(k\mathrm{\nabla }\widehat{T}·\mathrm{\nabla }\widehat{w}-{\overline{Q}}_i{\mathit{\varphi }}_i\widehat{w}\right)\mathrm{d}\mathrm{\Omega }+\underset{{\mathrm{\Gamma }}_2}{\int }\check{q}\widehat{w}\mathrm{d}\mathrm{\Gamma }=0\end{array}$$(21) ) approximately by holding $k^{(r)}$ and ${\check{q}}^{(r)}$ fixed in terms of the temperature from the rth inverse iteration leading to the following system after substituting in the finite element interpolations(24) $$\begin{array}{c}\hfill \underset{\mathrm{\Omega }}{\int }\left(k^{(r)}\mathrm{\Delta }{\overline{T}}_{\mathit{\alpha }}\mathrm{\nabla }{\mathit{\psi }}_{\mathit{\alpha }}·\mathrm{\nabla }{\mathit{\psi }}_{\mathit{\beta }}-\mathrm{\Delta }{\overline{Q}}_i{\mathit{\varphi }}_i{\mathit{\psi }}_{\mathit{\beta }}\right)\mathrm{d}\mathrm{\Omega }=0\end{array}$$(24)

or in matrix form(25) $$\begin{array}{c}\hfill {\mathit{K}}^{(r)}\mathrm{\Delta }\overline{\mathit{T}}=\mathit{H}\mathrm{\Delta }\overline{\mathit{Q}}\end{array}$$(25)

where(26a) $$\begin{array}{cc}\hfill K_{\mathit{\beta }\mathit{\alpha }}^{(r)}& =\underset{\mathrm{\Omega }}{\int }{k}^{(r)}\mathrm{\nabla }{\mathit{\psi }}_{\mathit{\beta }}·\mathrm{\nabla }{\mathit{\psi }}_{\mathit{\alpha }}\mathrm{d}\mathrm{\Omega }\hfill \\ \hfill H_{\mathit{\beta }i}& =\underset{\mathrm{\Omega }}{\int }{\mathit{\psi }}_{\mathit{\beta }}{\mathit{\varphi }}_i\mathrm{d}\mathrm{\Omega }\hfill \end{array}$$(26a) $\mathrm{\Delta }\overline{\mathit{T}}={\overline{\mathit{T}}}^{(r+1)}-{\overline{\mathit{T}}}^{(r)}$, and $\mathrm{\Delta }\overline{\mathit{Q}}={\overline{\mathit{Q}}}^{(r+1)}-{\overline{\mathit{Q}}}^{(r)}$. Note that ${\mathit{K}}^{(r)}$ is a symmetric, positive definite matrix.

The function to be minimized in Equation (23(23) $$\begin{array}{c}\hfill \underset{\overline{\mathit{Q}}}{min}\frac{1}{2}{\left(\mathit{R}\overline{\mathit{T}}(\overline{\mathit{Q}})-{\overline{\mathit{T}}}^d\right)}^T\left(\mathit{R}\overline{\mathit{T}}(\overline{\mathit{Q}})-{\overline{\mathit{T}}}^d\right)\end{array}$$(23) ) to update to iteration $(r+1)$ can now be rewritten as(27) $$\begin{array}{cc}\hfill m& =\frac{1}{2}{\left(\mathit{R}{\overline{\mathit{T}}}^{(r+1)}-{\overline{\mathit{T}}}^d\right)}^T\left(\mathit{R}{\overline{\mathit{T}}}^{(r+1)}-{\overline{\mathit{T}}}^d\right)\hfill \\ \hfill & =\frac{1}{2}{\left[\mathit{R}{({\mathit{K}}^{(r)})}^{-1}\mathit{H}\mathrm{\Delta }\overline{\mathit{Q}}-\left({\overline{\mathit{T}}}^d-\mathit{R}{\overline{\mathit{T}}}^{(r)}\right)\right]}^T\left[\mathit{R}{({\mathit{K}}^{(r)})}^{-1}\mathit{H}\mathrm{\Delta }\overline{\mathit{Q}}-\left({\overline{\mathit{T}}}^d-\mathit{R}{\overline{\mathit{T}}}^{(r)}\right)\right]\hfill \\ \hfill & =\frac{1}{2}{\left[{\mathit{S}}^{(r)}\mathrm{\Delta }\overline{\mathit{Q}}-\mathit{\delta }{\overline{\mathit{T}}}^{o(r)}\right]}^T\left[{\mathit{S}}^{(r)}\mathrm{\Delta }\overline{\mathit{Q}}-\mathit{\delta }{\overline{\mathit{T}}}^{o(r)}\right]\hfill \end{array}$$(27)

where ${\mathit{S}}^{(r)}=\mathit{R}{({\mathit{K}}^{(r)})}^{-1}\mathit{H}$ is the sensitivity matrix evaluated at iteration (r) and $\mathit{\delta }{\overline{\mathit{T}}}^{o(r)}={\overline{\mathit{T}}}^d-\mathit{R}{\overline{\mathit{T}}}^{(r)}$ is the distance the nodal temperatures are from the target value on ${\mathrm{\Gamma }}_2$ at iteration (r).

Inverse problems are typically ill-posed, which is manifested by instabilities in the solution. Overfitting of the prescribed temperature profile can lead to spatial oscillations in the heat source solution. Here, zeroth-order Tikhonov regularization [21] is used to limit the variability in the solution by augmenting the function to be minimized as follows(28) $$\begin{array}{c}\hfill m=\frac{1}{2}\left[{\left({\mathit{S}}^{(r)}\mathrm{\Delta }\overline{\mathit{Q}}-\mathit{\delta }{\overline{\mathit{T}}}^{o(r)}\right)}^T\left({\mathit{S}}^{(r)}\mathrm{\Delta }\overline{\mathit{Q}}-\mathit{\delta }{\overline{\mathit{T}}}^{o(r)}\right)+\mathit{\lambda }\mathrm{\Delta }{\overline{\mathit{Q}}}^T\mathrm{\Delta }\overline{\mathit{Q}}\right]\end{array}$$(28)

where $\mathit{\lambda }$ is the regularization parameter. Finally, minimizing m with respect to $\mathrm{\Delta }\overline{\mathit{Q}}$ yields our update formula(29) $$\begin{array}{c}\hfill \left[{({\mathit{S}}^{(r)})}^T{\mathit{S}}^{(r)}+\mathit{\lambda }\mathit{I}\right]\left({\overline{\mathit{Q}}}^{(r+1)}-{\overline{\mathit{Q}}}^{(r)}\right)={({\mathit{S}}^{(r)})}^T\mathit{\delta }{\overline{\mathit{T}}}^{o(r)}\end{array}$$(29)

that is used in our inverse analysis to update the heat supply vector ${\overline{\mathit{Q}}}^{(r+1)}$ at iteration $(r+1)$.

In summary, our inverse analysis algorithm is as follows:

1.	Set r = 0, specify an initial guess for ${\overline{\mathit{T}}}^{(0)}$ (for example, ${\overline{\mathit{T}}}^{(0)}={\overline{\mathit{T}}}_0$ the initial conditions in Equation (3(3) $$\begin{array}{cc}& T(\mathit{x},0)=T_0(\mathit{x})\text{in}\mathrm{\Omega }\hfill \end{array}$$(3) ));
2.	Calculate sensitivity matrix ${\mathit{S}}^{(r)}=\mathit{R}{({\mathit{K}}^{(r)})}^{-1}\mathit{H}$;
3.	Find ${\overline{\mathit{Q}}}^{(r+1)}$ using Equation (29(29) $$\begin{array}{c}\hfill \left[{({\mathit{S}}^{(r)})}^T{\mathit{S}}^{(r)}+\mathit{\lambda }\mathit{I}\right]\left({\overline{\mathit{Q}}}^{(r+1)}-{\overline{\mathit{Q}}}^{(r)}\right)={({\mathit{S}}^{(r)})}^T\mathit{\delta }{\overline{\mathit{T}}}^{o(r)}\end{array}$$(29) );
4.	Solve forward problem (13(13) $$\begin{array}{c}\hfill \mathit{a}(\dot{\overline{\mathit{T}}},\overline{\mathit{T}})=\mathbf{0}\end{array}$$(13) ) to steady-state for ${\overline{\mathit{T}}}^{(r+1)}$ with prescribed ${\overline{\mathit{Q}}}^{(r+1)}$;
5.	Calculate error norm $E^{(r+1)}=||{\overline{\mathit{T}}}^d-\mathit{R}{\overline{\mathit{T}}}^{(r+1)}||$;
6.	If $E^{(r+1)}<{10}^{-5}$ or $0.98<\frac{E^{p+1}}{E^p}<1.0$ ($p=r-4,\dots ,r$), stop and output ${\overline{\mathit{Q}}}^{(r+1)}$; else $r=r+1$ and go to step 2.

2.4. Passivity-based quasi-static feedback control

Considering the system as an input-output relationship, the forward steady-state finite element model can be represented as(30) $$\begin{array}{c}\hfill H_1:{\overline{\mathit{T}}}^o=\mathit{F}(\overline{\mathit{Q}})\end{array}$$(30)

where $\mathit{F}(·)$ is a memory-less function representing the mapping from applied heat sources $\overline{\mathit{Q}}$ to the output temperatures ${\overline{\mathit{T}}}^o$ related to $\overline{\mathit{T}}$. To retain the same dimension of input and output, ${\overline{\mathit{T}}}^o$ is chosen as the average temperature in each of 10 evenly separated regions on ${\mathrm{\Gamma }}_2$ above the 10 heater strips. In this way, the target temperature profile is reduced to a description with only 10 parameters.

Inspired by the energy dissipating nature of this thermal system, passivity, a powerful tool for stability analysis of non-linear systems, is first analysed. The memory-less system in Equation (30(30) $$\begin{array}{c}\hfill H_1:{\overline{\mathit{T}}}^o=\mathit{F}(\overline{\mathit{Q}})\end{array}$$(30) ) is defined to be passive (in this case, positive semi-definite) if ${({\overline{\mathit{T}}}^o)}^T\overline{\mathit{Q}}\ge 0$. It is input-strictly passive (positive definite) if ${({\overline{\mathit{T}}}^o)}^T\overline{\mathit{Q}}\ge \mathit{\beta }{\overline{\mathit{Q}}}^T\overline{\mathit{Q}}$. where $\mathit{\beta }$ is a positive scalar. Similar concepts are defined for dynamical systems, e.g. with state space description as(31) $$\begin{array}{c}\hfill H_2:\left\{\begin{array}{c}\dot{\mathit{x}}=\mathit{f}(\mathit{x},\mathit{u})\\ \mathit{y}=\mathit{h}(\mathit{x},\mathit{u})\end{array}\right.\end{array}$$(31)

where $\mathit{x}$ denotes the state, $\mathit{f}(·,·)$, $\mathit{h}(·,·)$ are non-linear functions, $\mathit{y}$ and $\mathit{u}$ are, respectively, output and input. The system is passive if there exists a continuously differentiable positive semidefinite function $V(\mathit{x})$ (called storage function) so that(32) $$\begin{array}{c}\hfill {\mathit{y}}^T\mathit{u}\ge \dot{V}=\frac{\mathit{\partial }V}{\mathit{\partial }\mathit{x}}\mathit{f}(\mathit{x},\mathit{u}).\end{array}$$(32)

The Passivity Theorem states that if two passive systems are connected as a negative feedback, ($\mathit{u}=\overline{\mathit{Q}}$, $\mathit{y}=-{\overline{\mathit{T}}}^o$), stability of the resulted system around the origin is guaranteed [22]. Moreover, if the memoryless non-linearity $H_1$ is input-strictly passive, then the connected system is asymptotically stable.

Suppose the target output is ${\overline{\mathit{T}}}^d$ and the corresponding input is ${\overline{\mathit{Q}}}^d$, then ${\overline{\mathit{T}}}^d=\mathit{F}({\overline{\mathit{Q}}}^d)$. Denoting $\mathit{\delta }{\overline{\mathit{T}}}^o={\overline{\mathit{T}}}^o-{\overline{\mathit{T}}}^d$, $\mathit{\delta }\overline{\mathit{Q}}=\overline{\mathit{Q}}-{\overline{\mathit{Q}}}^d$, then the system can be shifted so that ${\overline{\mathit{Q}}}^d$, ${\overline{\mathit{T}}}^d$ is the origin:(33) $$\begin{array}{c}\hfill \mathit{\delta }{\overline{\mathit{T}}}^o=\mathit{F}({\overline{\mathit{Q}}}^d+\mathit{\delta }\overline{\mathit{Q}})-\mathit{F}({\overline{\mathit{Q}}}^d)=\mathit{g}(\mathit{\delta }\overline{\mathit{Q}})\end{array}$$(33)

where $\mathit{g}(·)$ is also a memory-less function. Passivity of $\mathit{g}(·)$ is verified by randomly sampling $\mathit{\delta }\overline{\mathit{Q}}$ around numerous ${\overline{\mathit{T}}}^d$, input-strict passivity is found to hold for all the samples.

Now, we connect this system with another passive integral controller as negative feedback:(34) $$\begin{array}{c}\hfill \mathit{\delta }\dot{\overline{\mathit{Q}}}=-{\mathit{K}}^c\mathit{\delta }{\overline{\mathit{T}}}^o\end{array}$$(34)

where ${\mathit{K}}^c$ is a positive definite matrix. Then the origin is guaranteed to be asymptotically stable. This means that it’s stable and $\overline{\mathit{Q}}\to {\overline{\mathit{Q}}}^d$.

Applying this into a discrete control scheme, the integral controller can be transformed into:(35) $$\begin{array}{c}\hfill {\overline{\mathit{Q}}}^{(r+1)}={\overline{\mathit{Q}}}^{(r)}-{\mathit{K}}^c\mathrm{\Delta }{t}^c({\overline{\mathit{T}}}^o-{\overline{\mathit{T}}}^d)\end{array}$$(35)

where $\mathrm{\Delta }{t}^c$ is the time interval of the discrete system, which actually can be any value due to the quasi-static nature of the plant. Equation (35(35) $$\begin{array}{c}\hfill {\overline{\mathit{Q}}}^{(r+1)}={\overline{\mathit{Q}}}^{(r)}-{\mathit{K}}^c\mathrm{\Delta }{t}^c({\overline{\mathit{T}}}^o-{\overline{\mathit{T}}}^d)\end{array}$$(35) ) is the iterative formula used to update $\overline{\mathit{Q}}$, making the output ${\overline{\mathit{T}}}^o$ gradually approach the target ${\overline{\mathit{T}}}^d$.

Lastly, it remains to define the feedback control matrix ${\mathit{K}}^c$ and pseudo-time increment $\mathrm{\Delta }{t}^c$. We combine them into one matrix $\mathcal{K}={\mathit{K}}^c\mathrm{\Delta }{t}^c$. In practice, $\mathcal{K}$ can be specified as the following form(36) $$\begin{array}{c}\hfill \mathcal{K}={\mathcal{K}}_{diag}+\mathcal{L}\end{array}$$(36)

where ${\mathcal{K}}_{diag}$ is a diagonal with positive entries (known as decentralized feedback). $\mathcal{L}$ introduces an augmentation of consensus feedback [23]; it’s a weighted Laplacian of a graph we construct so that every heater is connected to its first, second, third and fourth neighbour. Specifically, $\mathcal{L}={\mathit{D}}^T\mathit{W}\mathit{D}$ where $\mathit{D}$ is incidence matrix of the graph stated above, and $\mathit{W}$ is a diagonal weighting matrix. The augmentation of $\mathcal{L}$ allows the feedback to penalize the difference between entries of $({\overline{\mathit{T}}}^o-{\overline{\mathit{T}}}^d)$, so that in the process of converging to ${\overline{\mathit{T}}}^d$, the shape of the intermediate profile is kept close to the desired one.

The iterative formula in Equation (35(35) $$\begin{array}{c}\hfill {\overline{\mathit{Q}}}^{(r+1)}={\overline{\mathit{Q}}}^{(r)}-{\mathit{K}}^c\mathrm{\Delta }{t}^c({\overline{\mathit{T}}}^o-{\overline{\mathit{T}}}^d)\end{array}$$(35) ) is different from the one in the inverse problem method, Equation (29(29) $$\begin{array}{c}\hfill \left[{({\mathit{S}}^{(r)})}^T{\mathit{S}}^{(r)}+\mathit{\lambda }\mathit{I}\right]\left({\overline{\mathit{Q}}}^{(r+1)}-{\overline{\mathit{Q}}}^{(r)}\right)={({\mathit{S}}^{(r)})}^T\mathit{\delta }{\overline{\mathit{T}}}^{o(r)}\end{array}$$(29) ), as the coefficient matrix $\mathcal{K}$ is fixed for all iterations while the analogous coefficient matrix $[{({\mathit{S}}^{(r)})}^T{\mathit{S}}^{(r)}+\mathit{\lambda }\mathit{I}]$ is updated at the beginning of each iteration in the inverse analysis. The feedback control algorithm follows the same basic steps as those in the inverse problem method listed at the end of Section 2.3, where Equation (35(35) $$\begin{array}{c}\hfill {\overline{\mathit{Q}}}^{(r+1)}={\overline{\mathit{Q}}}^{(r)}-{\mathit{K}}^c\mathrm{\Delta }{t}^c({\overline{\mathit{T}}}^o-{\overline{\mathit{T}}}^d)\end{array}$$(35) ) is used in place of Equation (29(29) $$\begin{array}{c}\hfill \left[{({\mathit{S}}^{(r)})}^T{\mathit{S}}^{(r)}+\mathit{\lambda }\mathit{I}\right]\left({\overline{\mathit{Q}}}^{(r+1)}-{\overline{\mathit{Q}}}^{(r)}\right)={({\mathit{S}}^{(r)})}^T\mathit{\delta }{\overline{\mathit{T}}}^{o(r)}\end{array}$$(29) ).

3. Results and discussion

3.1. Forward model analysis and design

Before solving the inverse problem to determine the required heat input into each heater strip in the micro-heater array, the geometric design parameters w, h, d and b must be specified (see Figure 1(a)). Of these, the most critical are w and h defining the cross-sectional dimensions of the heater strips. The cross section of the heater strips must be such that they can carry enough power per unit length to generate sufficient heat to cause grain growth in the Cu film while keeping the current density sufficiently small to avoid electromigration-induced failure in the Al and Ti lines. The power per unit length P is(37) $$\begin{array}{c}\hfill P=\frac{I^2\mathit{\rho }}{A}=Ij\mathit{\rho }\end{array}$$(37)

where I is the current in a line and $A=wh$ is a heater strip cross-sectional area. In addition, the lines should be sufficiently wide, or the gap between the lines $s=p-w$ (pitch $p=100\mathrm{\mu }$m), sufficiently small, so that the temperature profile is smooth.

Grain growth in copper films occurs on a time scale of tens of minutes in the range 300–${400}^{\circ }\mathrm{C}$, which is on par with the time scale required for experimental observations. Thus, the heater array should be able to bring the temperature in the copper to this range and should allow for a gradient of at least ${50}^{\circ }\mathrm{C}{\mathrm{mm}}^{-1}$ over the 1 mm domain so that a gradient in the microstructure may be created. To limit the required heat generation from the micro-heater array, the device can be placed on a heating stage that can bring the initial temperature to $T_0=200{}^{\circ }\mathrm{C}$ and maintain this temperature on ${\mathrm{\Gamma }}_1$, i.e. $\check{T}=200{}^{\circ }\mathrm{C}$. Then the micro-heater array only needs to raise the temperature $100{}^{\circ }\mathrm{C}$ to cause grain growth. From forward simulations, it was determined that the minimum required power per unit length to achieve this is $P>20\mathrm{W}{\mathrm{cm}}^{-1}$. In the temperature range of interest, $\mathit{\rho }\approx 100\times {10}^{-6}\mathrm{\Omega }\text{cm}$, and thus from Equation (37(37) $$\begin{array}{c}\hfill P=\frac{I^2\mathit{\rho }}{A}=Ij\mathit{\rho }\end{array}$$(37) ) $Ij>2\times {10}^5{\mathrm{A}}^2{\text{cm}}^{-2}$.

To avoid electromigration failure, the current density j must be limited and the surrounding ${\text{SiO}}_2$ should be sufficient to contain the line without cracking under stress due to electromigration-induced diffusion of either the Al or Ti. While there have been many studies of electromigration in Al conductor lines,[24,25] electromigration in Ti lines has not been studied. Based on available Al data and estimating the electromigration behavior in Ti, we estimate that electromigration will not be significant if the current density is kept at $j<{10}^6\text{A}{\mathrm{cm}}^{-2}$ provided the ${\text{SiO}}_2$ has a thickness of at least $1.5\mathrm{\mu }\mathrm{m}$ above the Ti lines. To satisfy the power requirement $P>20\mathrm{W}{\mathrm{cm}}^{-1}$ and the current density requirement $j<{10}^6\text{A}{\mathrm{cm}}^{-2}$,the current must be limited such that $I>0.2$ A and the cross-sectional area must be $A=wh>20\mathrm{\mu }{\mathrm{m}}^2$. A design that satisfies all of these requirements is $w=70\mathrm{\mu }\mathrm{m}$ ($s=30\mathrm{\mu }\mathrm{m}$), $h=0.4\mathrm{\mu }\mathrm{m}$, $d=5.9\mathrm{\mu }\mathrm{m}$, and $b=2.5\mathrm{\mu }\mathrm{m}$.

Now that the micro-heater array geometry is completely defined, it must be discretized into a finite element mesh for both the forward and inverse finite element analysis. The mesh must be sufficiently refined to maintain a satisfactory accuracy in the forward solution. A mesh convergence study was performed, and a mesh with 1540 four-node quadrilateral elements was found to be adequate to generate a solution accurate to within $3{}^{\circ }\mathrm{C}$. A graded mesh was used that is highly refined near the top where the heat is generated and coarse near the bottom of the Si substrate. The uncertainty due to the finite element approximation is relatively small compared to uncertainty in the physical model. For example, there is uncertainty in the layer thicknesses because of the fabrication process, and the resistivity in the Ti lines shows variability due to microstructure and imperfections. Model calibration to experiments will be done as the experiments come on-line, and the model uncertainty will be more precisely quantified.

Results from a transient forward simulation are shown in Figure 2 where a uniform current of $I=0.24$ A is prescribed in each heater strip for a current density of $j=8.57\times {10}^5\text{A}{\mathrm{cm}}^{-2}$. From this plot, we see that steady-state is achieved in about 10 ms, the temperature is elevated by $100{}^{\circ }\mathrm{C}$ to about $300{}^{\circ }\mathrm{C}$ in the central part of the region above the heater array, but drops off near the edges due to the heat conducting away on the sides. We also see that the temperature profile is not completely smooth with small bumps in the temperature above each heater and dips in between.

Figure 2. Evolution of the temperature distribution on the top surface ${\mathrm{\Gamma }}_2$ for a uniform current of $I=0.24$ A in each heater strip.

Figure 3. Temperature distribution on the top of the copper thin film from inverse analysis with different regularization parameters. The units for the regularization parameter are (${\text{K-mm}}^3{\text{W}}^{-1}{)}^2$.

3.2. Inverse problem solution

In this section, we compare the two methods for solving the inverse heat transfer source problem described in Sections 2.3 and 2.4, direct minimization and feedback control, solving problems relevant to the micro-heater array design. Specifically, we consider the cases of finding the required heater settings for generating a prescribed uniform temperature of $300{}^{\circ }\mathrm{C}$ and a linear gradient of $50{}^{\circ }\mathrm{C}{\mathrm{mm}}^{-1}$, where part of the distribution is above $300{}^{\circ }\mathrm{C}$.

The goal of the inverse problem is to find the heat source, and thus the current, that must be applied to each micro-heater in the array to achieve a desired temperature distribution on the top surface. In the inverse analysis approach with regularization, the regularization parameter $\mathit{\lambda }$ must be selected large enough to generate a smooth, stable solution, but not so large as to drive the solution away from satisfying the objective function. To show the effect of the regularization parameter, we first solve the inverse problem for a prescribed uniform temperature distribution of $T^d=300{}^{\circ }\mathrm{C}$ (573 K) on the top surface with varying values for $\mathit{\lambda }$. Since the regularization parameter is added on the diagonal to ${({\mathit{S}}^{(r)})}^T{\mathit{S}}^{(r)}$ in Equation (29(29) $$\begin{array}{c}\hfill \left[{({\mathit{S}}^{(r)})}^T{\mathit{S}}^{(r)}+\mathit{\lambda }\mathit{I}\right]\left({\overline{\mathit{Q}}}^{(r+1)}-{\overline{\mathit{Q}}}^{(r)}\right)={({\mathit{S}}^{(r)})}^T\mathit{\delta }{\overline{\mathit{T}}}^{o(r)}\end{array}$$(29) ), it should not be very large compared to those diagonal elements, which in this case are on the order of ${10}^{-5}$ (${\text{K-mm}}^3{\text{W}}^{-1}{)}^2$. Thus, regularization parameters considered were near this order of magnitude. The results are shown in Figures 3 and 4. Figure 3 shows the resulting temperature profile on the surface. As the regularization parameter $\mathit{\lambda }$ increases, the temperature at both ends decreases, and the overall shape of the temperature profile becomes smoother, though the small peaks with respect to the micro-heaters remain. While it is possible to obtain a solution close to the prescribed temperature near the centre of the domain, it is not possible near the edges. When the regularization parameter is too high, the solution starts to move further away from the target temperature $T^d$. There is almost no difference between the result without regularization and when the regularization parameter is set at $\mathit{\lambda }={10}^{-6}$ (${\text{K-mm}}^3{\text{W}}^{-1}{)}^2$.

In addition, we calculated the condition number of $\mathit{S}$ at convergence and found it to be 9.1, which indicates very weak instability of solution [26] and is consistent with the relatively small oscillations observed even with no regularization. In Figure 4, the current input I to each heater is plotted for the cases with $\mathit{\lambda }=0$ and $\mathit{\lambda }={10}^{-5}{({\text{K-mm}}^3{\text{W}}^{-1})}^2$, where I is computed from the heat source Q solved for in each heater using Equation (1(1) $$\begin{array}{c}\hfill Q=\mathit{\rho }{j}^2\mathrm{i}n{\mathrm{\Omega }}_4\end{array}$$(1) ) and $j=I/A$. Without regularization, oscillations in the required current occur at the ends with a large spike in the current predicted for the heaters at the ends. Such spikes are undesirable from a practical point of view because high current densities may lead to rapid failure of the lines. The current distribution near both ends becomes smoother as $\mathit{\lambda }$ changes from 0 to ${10}^{-5}{({\text{K-mm}}^3{\text{W}}^{-1})}^2$.

We also wanted to see if we could get a good solution ignoring the non-linearity, and instead, solve a linear system setting the material parameters at the target temperature of $300{}^{\circ }\mathrm{C}$. When we did this, we found the solution for the optimal current input to each line deviated by as much as 30% from the value found including the non-linearity. Thus, for an accurate solution, the non-linear material behaviour is very important.

Figure 4. Currents in the 10 heaters from the solution of inverse analysis for $\mathit{\lambda }=0$ and $\mathit{\lambda }={10}^{-5}$. The units for the regularization parameter are (${\text{K-mm}}^3{\text{W}}^{-1}{)}^2$.

To quantitatively evaluate the deviation of the solution ${\overline{\mathit{T}}}^s=\mathit{R}\overline{\mathit{T}}$ from the prescribed temperature ${\overline{\mathit{T}}}^d$ on the copper thin film surface, we set a metric ${\mathit{\sigma }}^{\ast }$ by calculating the standard deviation as(38) $$\begin{array}{c}\hfill {\mathit{\sigma }}^{\ast }=\sqrt{\frac{\sum \Vert {\overline{\mathit{T}}}^s-{\overline{\mathit{T}}}^d{\Vert }^2}{N}}\end{array}$$(38)

Figure 5 shows the metric ${\mathit{\sigma }}^{\ast }$ for varying regularization parameters. In addition, ${\mathit{\sigma }}^{\ast }=9.322$ K for the case with the regularization parameter $\mathit{\lambda }$ = 0. It is shown that ${\mathit{\sigma }}^{\ast }$ increases with $\mathit{\lambda }$, which means the solution increasingly deviates from the prescribed temperature as $\mathit{\lambda }$ increases. Therefore, there is a trade-off between the smoothness of the profile and the minimization of the metric ${\mathit{\sigma }}^{\ast }$. However, we also see that ${\mathit{\sigma }}^{\ast }$ does not increase significantly until $\mathit{\lambda }>{10}^{-5}{({\text{K-mm}}^3{\text{W}}^{-1})}^2$, and thus, this appears to be an optimal value. Solving the inverse problem with $\mathit{\lambda }={10}^{-5}{({\text{K-mm}}^3{\text{W}}^{-1})}^2$ leads to a relatively small metric ${\mathit{\sigma }}^{\ast }$ ($\approx$10 K) and a relatively uniform temperature profile on the surface of the area of interest.

Figure 5. Comparison of standard deviation of ${\mathit{\sigma }}^{\ast }$ (K) as regularization parameter $\mathit{\lambda }$ (${\text{K-mm}}^3{\text{W}}^{-1}{)}^2$ varies.

To solve the inverse problem with the feedback control strategy, it remains only to define the diagonal matrices ${\mathcal{K}}_{diag}$ and $\mathit{W}$. Both of these matrices are set with constant values on the diagonals, with the diagonal elements of ${\mathcal{K}}_{diag}$ set at $200\text{W}/({\mathrm{mm}}^3$-K) and those of $\mathit{W}$ set at $30\text{W}/({\text{mm}}^3$-K).

We now compare the results and performance for the inverse problem and feedback control approaches for two cases, one with the objective of creating a uniform temperature profile, as described earlier when investigating the effect of the regularization parameter, and one with the objective of creating a linear temperature profile. Although the inverse analysis and the feedback control approach both use an iterative scheme to find the optimal solution for the heat sources, the difference in the update equations leads to dissimilarity in their performance.

Figures 6 and 7 present the results for the case with the target of a uniform temperature profile. In Figure 6, the resulting temperature profiles are shown. The results in the central region for both methods are nearly identical, and the results only deviate near the ends, above the last two heaters at each end. The feedback control generates a higher peak at the end, which is similar to the inverse analysis result without regularization shown in Figure 3. In Figure 7, the current input I to each heater is plotted. We see that the two micro-heaters near both ends of the area of interest have relatively large values, especially for the solution from the feedback control. As mentioned in Section 3.1, we would like to keep the current density $j<{10}^6\text{A}{\mathrm{cm}}^{-2}$ to reduce the likelihood for failure by electromigration for the given cross-section geometries, that corresponds to $I<280$ mA. We can see in Figure 7 that the currents exceed the recommended limit in the last heaters on the ends because these heaters have to provide more heat to compensate for loses due to lateral conduction. This may not be a major issue as the resistivity $\mathit{\rho }$ in the Ti lines assumed here is a reported value for bulk Ti, and the resistivity in the 400 nm thick heater strips is likely to be higher, thus requiring less current to generate the same amount of heating. In fact, preliminary experimental measurements of the resistance in the Ti lines for fabricated micro-heaters shows a resistance of approximate twice the bulk value. In addition, a relatively thick layer of ${\text{SiO}}_2$ is used, which also will act to hinder electromigration.

Figure 6. Comparison of temperature distribution of nodes on the top of the copper thin film between the inverse analysis and feedback control approachs.

Figure 7. Comparison of currents in the 10 heaters between the inverse analysis and the feedback control approachs.

As for the computational efficiency, the inverse analysis with regularization converges to the solution faster than the feedback control method. Figure 8 shows that the inverse analysis converges after five iterations, and for the feedback control it takes about eight iterations. Most likely this is due to their difference in solving the iterative formula. In Equation (29(29) $$\begin{array}{c}\hfill \left[{({\mathit{S}}^{(r)})}^T{\mathit{S}}^{(r)}+\mathit{\lambda }\mathit{I}\right]\left({\overline{\mathit{Q}}}^{(r+1)}-{\overline{\mathit{Q}}}^{(r)}\right)={({\mathit{S}}^{(r)})}^T\mathit{\delta }{\overline{\mathit{T}}}^{o(r)}\end{array}$$(29) ), e.g. the iterative formula of the inverse analysis, the coefficient matrix ${({\mathit{S}}^{(r)})}^T{\mathit{S}}^{(r)}+\mathit{\lambda }\mathit{I}$ is updated in each iteration. However, for feedback control, Equation (35(35) $$\begin{array}{c}\hfill {\overline{\mathit{Q}}}^{(r+1)}={\overline{\mathit{Q}}}^{(r)}-{\mathit{K}}^c\mathrm{\Delta }{t}^c({\overline{\mathit{T}}}^o-{\overline{\mathit{T}}}^d)\end{array}$$(35) ) uses a predefined feedback matrix $\mathcal{K}$, which is fixed in all iterations. Therefore, the inverse analysis captures the temperature dependence of the model better than the feedback control, speeding up the convergence. On the other hand, fewer computations are required in each iteration for the feedback control strategy because $\mathcal{K}$ only needs to be evaluated once.

Figure 8. Comparison of the standard deviation ${\mathit{\sigma }}^{\ast }$ at different inverse analysis or feedback control steps.

Figure 9. Comparison between the inverse analysis results using $\mathit{\lambda }={10}^{-5}$ (${\text{K-mm}}^3{\text{W}}^{-1}{)}^2$ and feedback control results for reaching a linear temperature profile.

In addition, we compared the performance between the inverse analysis with regularization and feedback control methods for reaching a linear temperature profile as shown in Figure 9. Here, the prescribed temperature distribution is a linear distribution that varies over the range 275–$325{}^{\circ }\mathrm{C}$ (547–598 K), i.e. a temperature gradient of $50{}^{\circ }\mathrm{C}$ across the 1 mm area of interest. The solution of the inverse analysis is regularized by the same regularization parameter ($\mathit{\lambda }={10}^{-5}$ (${\text{K-mm}}^3{\text{W}}^{-1}{)}^2$) as before. From Figure 9, both the inverse analysis and feedback control are able to drive the temperature at the centre region close to the prescribed temperature profile $T^d$ in the central region with peaks and a rapid drop-off near the ends, similar to the results for the constant profile case. As in the previous case, the feedback control method has larger spikes at the ends, which are smoothed with regularization.

4. Conclusion

In this paper, we present two methods to solve the inverse heat source problem to find the optimal current inputs to elements in a micro-heater array to achieve a prescribed temperature field that may be used to control microstructure evolution. First, we show that the micro-heater array reaches a steady-state temperature field in about 10 ms, which is negligible compared to the grain growth experiment time frame of minutes and hours. Thus, the feed-forward control problem of finding the current as a function of time that should be applied to each heater element to achieve a prescribed temperature profile history, can be reduced to solving the steady-state problem at discrete times. A direct inverse method with zeroth-order Tikhonov regularization and a passivity-based feedback control method were developed, used to solve two cases of the inverse problem, and compared. The methods are based on a non-linear finite element model of the micro-heater array. Both methods performed well.

The accuracy of the finite element model in representing the actual physical system must be tested. For example, early tests show that the resistivity in the Ti heater lines is roughly twice that reported for bulk Ti used here. Thus, it is expected that the forward finite element model will require refinement as data becomes available. The experimental data can also be used to quantify the model uncertainty which may be used to estimate uncertainty in the solution for the required heater settings using Bayesian methods.

Nevertheless, the model has proven useful in generating an initial design for a micro-heater array, and the methods to solve the inverse problem can be applied to the refined model to generate improved estimates for the required currents for different temperature profiles for feed-forward control. The results show that both a uniform and a linear temperature distribution can be attained near the central part of the region of interest, but it is not possible to match the temperatures well at the edges of the region. If we want to accurately control the temperature field over the full 1 mm region, the heater array must extend beyond this region with either more or wider heater elements. Comparing the inverse analysis with the feedback control approach, we found that the inverse analysis generated a smoother distribution of current inputs and resulting temperature field due to the regularization. The feedback control algorithm was relatively simpler because the coefficient matrix used to solve for the heat source update was fixed for all iterations, but required more iterations to converge.

5. Nomenclature

• Q   volumetric heat generation rate

• T   temperature

• $\overline{\mathit{Q}}$   vector of heat sources in each heater strip

• $\mathit{\rho }$   electrical resistivity

• j   current density in a heater strip

• ${\mathrm{\Omega }}_1$   Cu region

• ${\mathrm{\Omega }}_2$   ${\text{SiO}}_2$ region

• ${\mathrm{\Omega }}_3$   Si region

• ${\mathrm{\Omega }}_4$   ${\text{TiO}}_2$ region

• $\mathrm{\Omega }$   union of all regions

• ${\mathrm{\Gamma }}_1$   bottom boundary

• ${\mathrm{\Gamma }}_2$   top boundary

• ${\mathrm{\Gamma }}_3$   side boundaries

• $\mathit{q}$   heat flux

• k   thermal conductivity

• C   volumetric heat capacity

• $\mathit{\epsilon }$   emissivity

• $\mathit{\sigma }$   Stefan–Boltzmann constant

• $\mathit{n}$   outward unit normal on the boundary

• w   width of the Ti heating elements

• h   thickness of the Ti heating elements

• d   thickness of the ${\text{SiO}}_2$ layer

• b   position of the Ti heaters in the ${\text{SiO}}_2$ layer

• $\widehat{T}$   finite dimensional approximation of temperature field

• ${\mathit{\psi }}_{\mathit{\alpha }}$   finite element basis functions

• $N_T$   number of temperature degrees of freedom

• $\widehat{w}$   finite dimensional weighting function

• $\overline{\mathit{T}}$   vector of nodal temperatures

• $\dot{\overline{\mathit{T}}}$   time derivative of $\overline{\mathit{T}}$

• ${\overline{\mathit{T}}}^d$   vector of prescribed nodal temperatures on ${\mathrm{\Gamma }}_2$

• $\check{q}$   heat flux due to radiation on ${\mathrm{\Gamma }}_2$

• $N_H$   number of heater strips

• $\mathit{R}$   matrix that extracts nodal temperature on ${\mathrm{\Gamma }}_2$ from $\overline{\mathit{T}}$

• $\mathit{S}$   sensitivity matrix

• $\mathit{\lambda }$ regularization parameter

• $T_{\infty }$   ambient temperature

• $\mathit{I}$   identity matrix

• ${\overline{\mathit{T}}}^o$   vector of average temperature on ${\mathrm{\Gamma }}_2$ above each heater strip

• ${\mathit{K}}^c$   feedback control matrix

• P   power per unit length of a heater strip

• A   heater strip cross-sectional area

• p   pitch of heater strips

• $T_0$   initial temperature

• I   current in a heater strip

• s   gap between heater strips

• r   iteration number

• $\mathcal{K}$   coefficient matrix in feedback control

• $\mathit{y}$   output variable in feedback control

• $\mathit{x}$   state variable in feedback control

• $\mathit{u}$   input variable in feedback control

• ${\overline{\mathit{T}}}^{\mathit{s}}$   nodal temperature on ${\mathrm{\Gamma }}_2$ calculated from finite element model

• ${\mathit{\sigma }}^{\ast }$   metric to evaluate the deviation of ${\overline{\mathit{T}}}^s$ from ${\overline{\mathit{T}}}^d$

 

Additional information

Funding

This research is sponsored through a grant from the National Science Foundation, Award CMMI-1334283, DMREF: Real Time Control of Grain Growth in Metals. The authors acknowledge our co-workers on this grant: Professors Robert Hull and Dan Lewis, postdoctoral fellow Prabhu Balasubramanian, and graduate students David Christ, Genevieve Kane and Christina Gilliland.

Notes

No potential conflict of interest was reported by the authors.