# PDE Modeling of a Microfluidic Thermal Process for Genetic Analysis Application

Authors: Reza Banaei Khosroushahi; Horacio J. Marquez; Jose Martinez-Quijada; Christopher J. Backhouse
Journal: Journal of Applied Mathematics
 (2013)

---

## Abstract

This paper details the infinite dimensional dynamics of a prototype microfluidic thermal process that is used for genetic analysis purposes. Highly effective infinite dimensional dynamics, in addition to collocated sensor and actuator architecture, require the development of a precise control framework to meet the very tight performance requirements of this system, which are not fully attainable through conventional lumped modeling and controller design approaches. The general partial differential equations describing the dynamics of the system are separated into steady-state and transient parts which are derived for a carefully chosen three-dimensional axisymmetric model. These equations are solved analytically, and the results are verified using an experimentally verified precise finite element method (FEM) model. The final combined result is a framework for designing a precise tracking controller applicable to the selected lab-on-a-chip device.


---

## Body

## 1. Introduction

Biomedical microdevices usually require the very tight performance specification to be satisfied so that they become properly functional. Carefully designed controllers can guarantee certain performance specifications such as robustness and stability. Successful controller design process requires a precise mathematical model. A perfect mathematical model is the one that captures the highest amount of the dynamics of the system. Likewise, numerical tools can provide precise insights into a system’s dynamics. Finite element method (FEM) is a common numerical tool that can provide precise simulation of the steady-state and transient behaviour of a microdevice. However, FEM models are not suitable for controller design because they are not in a compact mathematical form. In this work, we present a mathematical model that captures infinite dimensional dynamics of the heat distribution inside a prototype system. The selected system is a lab-on-a-chip (LOC) device that is designed and built to perform the Polymerase Chain Reaction (PCR) and capillary electrophoresis (CE) processes. The three-dimensional complete FEM model of the system is used as the reference to derive the mathematical model. The FEM model is experimentally verified with the actual setup.Lab-on-a-chip (LOC) devices are finding growing attention in various applications ranging from chemical analysis and water contamination sensors to genetic analysis [1]. A LOC device can be thought of as a miniaturized instrument capable of performing single or multiple laboratory processes that integrate various functionalities such as sample pretreatment, sample transportation, genetic amplification, mixing, reaction, separation, and detection. LOC devices offer advantages beyond conventional biochemical and clinical laboratory settings. Dramatic reduction of sample volume and reduced consumption of costly reagents are two distinct benefits of miniaturized laboratory processes. Most significantly, miniaturization improves analysis characteristics, such as shorter analysis time and higher separation performance, and can even enable innovative applications that are not otherwise attainable.Genetic and molecular analyses are among the major applications of LOC devices and have become feasible by the adoption of the Polymerase Chain Reaction (PCR) technique. Using PCR, specific regions of even a single DNA molecule can be replicated up to about 107 to 108 times. This process is also called selective amplification of DNA. The amplified DNA segments can be later analyzed and compared to a standard pattern to find out if they correspond to some disease condition. PCR has been traditionally a lengthy molecular biology operation, but the LOC concept, proposed in the early 1990s [2, 3], has pushed this technology to attain ever shorter reaction times. DNA amplification using PCR involves repeatedly cycling the temperature of the reaction chamber in three predefined temperatures. Optimum reaction temperature and exposure times of the sample to each temperature are key factors for successful DNA amplification [4]. The PCR-LOC research community emphasizes the importance of ±1°C temperature precision and rapid transition time for high amplification efficiency. At the same time, overshoots and undershoots, which are likely to happen at higher speeds, have a detrimental impact on the amplification yield. For example, overshoots at denaturation temperature will reduce the lifetime of replicating enzyme (Taq polymerase). Undershoots at annealing temperature may cause the amplification of unspecified products. Nevertheless, rapid thermal cycling of reaction chamber in DNA amplification not only reduces the amplification time, but it is also important to ensure the stability of biological reagents during amplification as the Taq polymerase has a half lifetime that decreases rapidly at temperatures close to 100°C [5]. It is important to note that about half of the time in a standard PCR process is spent in transitions which are undesirable. As a result, researchers have focused their attention to reduce transition time between temperature steps in PCR process. To achieve this goal, both structural design of the device and thermal control of the process need to be optimized. Consequently, various platforms for PCR-LOC devices have been reported in the literature [6] which can be roughly divided into two major categories: continuous-flow PCR and static PCR [7]. Furthermore, different heating techniques have been chosen, where thermoelectric modules (TEMs) and thin-film resistive heaters [8–11] are the most commonly used techniques. Other techniques, such as infrared radiation [12] and microwaves [13, 14], have also been reported. These different structural designs are discussed in detail in the excellent review papers by Zhang et al. [7], Roper and coworkers [6], Verpoorte [15], and DeMello [16]. The latest achievements in PCR-LOC devices design can be found in the recently published survey paper by Zhang and Ozdemir [1].Most of the published work on PCR-LOC devices has been devoted to the structural design aspects, and less attention has been paid to optimize the thermal control of the process. As a result, there is very little work in the literature on the dynamics and control of these microsystems. Although neural network based controllers [17] and hybrid controllers [18] are used for thermal control of PCR-LOC devices, PI/PID controllers are still the most commonly preferred controllers to deal with the thermal control problem of PCR-LOCs [5, 19, 20]. What makes the thermal control of most of the PCR-LOC devices a challenging task is the fact that the temperature inside the PCR-LOC chips not only depends on time but is also spatially distributed. Thus lumped modeling and controller design cannot exactly address the problem, fundamentally because in this approach the distributed parameter nature of the system is neglected. We observed this critical problem to some extent in the formerly designed PCR-LOC chip [21] where switching controllers were employed [8]. The PCR-LOC system studied in this paper is the next generation of the one described in our previous work [22]. The new design has undergone significant improvement in structural design with respect to its predecessor.While lumped modeling can precisely describe the dynamics of lumped systems, it usually cannot provide a fair description of the dynamics of the distributed parameter systems. Those systems that can be described by partial differential equations are often referred to as distributed parameter systems or PDE systems which feature distributed dynamics and are categorized under a general set of the systems known as complex dynamical systems [23]. Distributed parameter systems are subsets of a more general set called infinite dimensional systems because they can be mathematically formulated as differential equations on an abstract linear vector space of infinite dimension [24]. In contrast, lumped systems can be formulated mathematically as ordinary differential equations on a finite dimensional linear vector space. Designing controllers for distributed parameter systems is called infinite dimensional controller design or PDE controller design [24–26].This paper characterizes the heat distribution dynamics of a microfluidic thermal process. The work begins with a description of the process and experimental setup. Then results of FEM simulations on the system were analyzed to choose a suitable structure for the model. Next the PDE equations describing the model in steady-state and transient parts were derived. These equations were solved analytically from which the dynamic temperature distribution in the whole domain can be expressed by an infinite dimensional model. Next, we optimize the calculated model to minimize the error due to the simplification of the actual multidomain system to an abstract single-domain system. FEM simulation and laboratory experiments are used to verify the model. The result of this work is a framework for boundary control design for thermal control of the aforementioned PCR-LOC device.


## 2. Materials and Methods

### 2.1. Chip Configuration

Genetic analysis usually consists of three steps: sample preparation, amplification, and detection. In the first step the DNA is extracted from a raw sample such as blood. Then, a selected fragment of this DNA is amplified through the PCR process to obtain sufficient DNA copies to perform detection. Next, the amplified region of DNA is labelled with a fluorescent dye, and finally it is detected through various analysis techniques [27].The integrated genetic analysis platform that has been developed in the Applied Miniaturization Laboratory (AML) at the University of Alberta is capable of performing complete genetic analysis on a single microfluidic chip [28]. This system is depicted in Figure 1. The microchip is composed of two layers of Borofloat glass with a polydimethylsiloxane (PDMS) layer in the middle. Reaction chamber and fluid channels are carved on one of the glass layers, and a Pt heater/sensor is patterned on the bottom glass plate using standard microfabrication techniques. A close-up view of the microfluidic chip, depicted in Figure 2, clearly shows the configuration of the platinum heater-sensor ring. A schematic representation of the components within the chip is shown in Figure 2(a). Perspective view of the glass/PDMS/glass chip with horizontal laser beam is shown in Figure 2(b) and the chip atop a heatsink and light detector is depicted in Figure 2(c). To improve the cooling rate during thermal cycling, the microchip is mounted on a copper heat sink. The heat sink has a circular opening aligned to the center of the PCR chamber and heater/sensor, ring axis. The radius of this opening is optimized to increase the cooling rate without having negative effect on the heating rate [21]. We are mostly interested in heat transfer dynamics around the reaction chamber and heater-sensor ring. The chamber, heater/sensor and heat-sink opening are circular and share the same axis. This suggests that we should look for an axisymmetric cylindrical structure for our model. The half cross-section view of the chip along a vertical plane parallel to the smaller side and passing through the center of the chip is shown in Figure 3. The irregular and multilayer multimaterial structure of the PCR-LOC system is not suitable to form a set of solvable analytic equations; additional simplifications are needed. To find a valid simplified model describing the dynamics of the chip where the heater and chamber are located, we first have to analyze the temperature profile around the heater and chamber area. Figure 4 depicts the temperature profile along the chip thickness for some selected distances from the center of the chip, ranging from 2.5 mm to 5.0 mm. Referring to Figure 4, the temperature profile along the thickness of the chip at the upper glass layer (1.354 mm to 2.454 mm from the bottom of the chip) shows a semilinear trend. By choosing the radius of 3.5 mm from the center of the chip, the slope is negligible, as the temperature varies for less than 0.8°C, which is even less than the acceptable error for the chamber temperature. Therefore, with reasonable approximation, we can consider a constant temperature wall as the boundary, where this boundary is positioned at a distance of 3.5 mm from center of the chip. These calculations, in addition to several other FEM simulations, have been performed using COMSOL Multiphysics 3.5a on a comprehensive FEM model of the chip that describes its layered geometry, material properties, and boundary conditions with a high degree of accuracy.  More FEM simulations disclosed some other important facts about this chip which can be summarized as follows.(1)
The amount of heat going to the upper portion of the chip through the PDMS and top glass layer has an approximately constant ratio to the total heat being dissipated by the heater.(2)
If we consider a disk-shape bounded volume at the center of the chip, the amount of heat lost by natural convection in air through the upper surface of the disk is negligible compared to the amount of heat lost by conduction through the wall.(3)
If the PDMS layer was replaced by a glass layer, only a small percentage of the dynamic temperature distribution would be affected, due to the alteration of thermal properties (thermal conductivity and specific heat capacity). This small effect on the dynamics can be compensated later by using a correction factor.(4)
The vertical temperature profile in the top glass layer at a distance of 3.5 mm from the center is approximately constant.Figure 1
Integrated genetic analysis platform comprising the PCR-LOC device and the control electronics. The platform is operated through a graphical user interface [29].The modelled PCR-LOC device. On (a) a schematic of the chip, showing the PCR chamber as a yellow circle. The chamber is surrounded by the heater ring. The squares in the middle are the heater contact pads. The crossing channels form the CE section. On (b) a 3D view of the chip showing the three layers that comprise the device. On (c) a cross-section view of the chip showing the chip placed on the heatsink that provides passive cooling [28].
(a)
(b)
(c)Figure 3
Cross-section view and dimensions of PCR LOC microchip. Due to the axisymmetric geometry of the chip, only half the chip is shown. The axis of symmetry is along the left vertical side. A circular hole is cut in the middle of the heatsink to provide thermal insulation under the heater and chamber. Heat is lost by free convection on the top and by conduction through the heatsink [21].Figure 4
Steady-state vertical temperature profile predicted by FEM. Each curve shows the temperature along a vertical line that crosses through the chip. Such a line is placed at a defined distance from the centre. The legend shows the radial distance in millimetres. The horizontal axis represents the vertical distance from the bottom to the top of the chip, which is 2.454 mm thick. As the distance from the center increases from 2.5 to 5 mm, the curves flattened out towards 22°C. The temperatures in the top glass plate (right side of the curves) are almost constant.


### 2.2. Model Structure

Using these observations we can conclude that, by picking a single layer glass disk, we can model the temperature distribution at a closed disk in the upper side of the chip with a good level of approximation relative to the temperature accuracy of ±1°C demanded by PCR-LOC systems. This model structure is shown in Figure5 where the heater is defined by the circular strip at the bottom surface with k2 width and k1 distance to the vertical axis. It is worth noting that we only modelled the area of interest in the chip, and we have put aside the rest of the chip. Furthermore, we defined a comparison aperture that can be characterized by a fixed size rectangular area in cross-section view of the actual PCR-LOC microchip and our model. The comparison aperture contains the reaction chamber, and it is discussed in detail in the simulation results section. We chose this scenario to arrive at a model structure which is simple enough to find its analytical solution meanwhile it is accurate enough to characterize the actual thermal process. A simple estimate of the thermal capacity and radial thermal resistance of the PDMS layer shows that they are small compared to those of the glass layers. It is expected that treating the PDMS layer as if it has the same thermal properties of glass will make the systems dynamic behavior 9% too fast (because the radial resistance is predicted to be 10% too low and the thermal capacity is estimated to be about 1% too high). It is possible to correct for this by scaling the thermal parameters of the glass. The assumptions made will result in underestimating the vertical gradient, but that gradient cannot be controlled by control mechanism and must instead be dealt with at the design and simulation level (i.e., dynamic simulations of the present design ensured that the PCR chamber radius and height were sufficiently small that the vertical gradient could not lead to temperature variations greater than ±1°C). The FEM model by itself was verified by extensive experiments, wherein thermochromic liquid crystal (TLC) sensitive to different temperatures was deposited in the chamber to read the temperature within its volume in steady-state conditions.Figure 5
Geometrical representation of the PDE model. Taking advantage of the axial symmetry of the chip, the entire system is represented as a two-dimensional rectangular region. Similarly the heater is represented as a line instead of being a surface. This reduction of one dimension allows the model to be solved for the temperature in a three-dimensional space but at a small fraction of the time and computational resources required by a FEM simulation.


## 2.1. Chip Configuration

Genetic analysis usually consists of three steps: sample preparation, amplification, and detection. In the first step the DNA is extracted from a raw sample such as blood. Then, a selected fragment of this DNA is amplified through the PCR process to obtain sufficient DNA copies to perform detection. Next, the amplified region of DNA is labelled with a fluorescent dye, and finally it is detected through various analysis techniques [27].The integrated genetic analysis platform that has been developed in the Applied Miniaturization Laboratory (AML) at the University of Alberta is capable of performing complete genetic analysis on a single microfluidic chip [28]. This system is depicted in Figure 1. The microchip is composed of two layers of Borofloat glass with a polydimethylsiloxane (PDMS) layer in the middle. Reaction chamber and fluid channels are carved on one of the glass layers, and a Pt heater/sensor is patterned on the bottom glass plate using standard microfabrication techniques. A close-up view of the microfluidic chip, depicted in Figure 2, clearly shows the configuration of the platinum heater-sensor ring. A schematic representation of the components within the chip is shown in Figure 2(a). Perspective view of the glass/PDMS/glass chip with horizontal laser beam is shown in Figure 2(b) and the chip atop a heatsink and light detector is depicted in Figure 2(c). To improve the cooling rate during thermal cycling, the microchip is mounted on a copper heat sink. The heat sink has a circular opening aligned to the center of the PCR chamber and heater/sensor, ring axis. The radius of this opening is optimized to increase the cooling rate without having negative effect on the heating rate [21]. We are mostly interested in heat transfer dynamics around the reaction chamber and heater-sensor ring. The chamber, heater/sensor and heat-sink opening are circular and share the same axis. This suggests that we should look for an axisymmetric cylindrical structure for our model. The half cross-section view of the chip along a vertical plane parallel to the smaller side and passing through the center of the chip is shown in Figure 3. The irregular and multilayer multimaterial structure of the PCR-LOC system is not suitable to form a set of solvable analytic equations; additional simplifications are needed. To find a valid simplified model describing the dynamics of the chip where the heater and chamber are located, we first have to analyze the temperature profile around the heater and chamber area. Figure 4 depicts the temperature profile along the chip thickness for some selected distances from the center of the chip, ranging from 2.5 mm to 5.0 mm. Referring to Figure 4, the temperature profile along the thickness of the chip at the upper glass layer (1.354 mm to 2.454 mm from the bottom of the chip) shows a semilinear trend. By choosing the radius of 3.5 mm from the center of the chip, the slope is negligible, as the temperature varies for less than 0.8°C, which is even less than the acceptable error for the chamber temperature. Therefore, with reasonable approximation, we can consider a constant temperature wall as the boundary, where this boundary is positioned at a distance of 3.5 mm from center of the chip. These calculations, in addition to several other FEM simulations, have been performed using COMSOL Multiphysics 3.5a on a comprehensive FEM model of the chip that describes its layered geometry, material properties, and boundary conditions with a high degree of accuracy.  More FEM simulations disclosed some other important facts about this chip which can be summarized as follows.(1)
The amount of heat going to the upper portion of the chip through the PDMS and top glass layer has an approximately constant ratio to the total heat being dissipated by the heater.(2)
If we consider a disk-shape bounded volume at the center of the chip, the amount of heat lost by natural convection in air through the upper surface of the disk is negligible compared to the amount of heat lost by conduction through the wall.(3)
If the PDMS layer was replaced by a glass layer, only a small percentage of the dynamic temperature distribution would be affected, due to the alteration of thermal properties (thermal conductivity and specific heat capacity). This small effect on the dynamics can be compensated later by using a correction factor.(4)
The vertical temperature profile in the top glass layer at a distance of 3.5 mm from the center is approximately constant.Figure 1
Integrated genetic analysis platform comprising the PCR-LOC device and the control electronics. The platform is operated through a graphical user interface [29].The modelled PCR-LOC device. On (a) a schematic of the chip, showing the PCR chamber as a yellow circle. The chamber is surrounded by the heater ring. The squares in the middle are the heater contact pads. The crossing channels form the CE section. On (b) a 3D view of the chip showing the three layers that comprise the device. On (c) a cross-section view of the chip showing the chip placed on the heatsink that provides passive cooling [28].
(a)
(b)
(c)Figure 3
Cross-section view and dimensions of PCR LOC microchip. Due to the axisymmetric geometry of the chip, only half the chip is shown. The axis of symmetry is along the left vertical side. A circular hole is cut in the middle of the heatsink to provide thermal insulation under the heater and chamber. Heat is lost by free convection on the top and by conduction through the heatsink [21].Figure 4
Steady-state vertical temperature profile predicted by FEM. Each curve shows the temperature along a vertical line that crosses through the chip. Such a line is placed at a defined distance from the centre. The legend shows the radial distance in millimetres. The horizontal axis represents the vertical distance from the bottom to the top of the chip, which is 2.454 mm thick. As the distance from the center increases from 2.5 to 5 mm, the curves flattened out towards 22°C. The temperatures in the top glass plate (right side of the curves) are almost constant.


## 2.2. Model Structure

Using these observations we can conclude that, by picking a single layer glass disk, we can model the temperature distribution at a closed disk in the upper side of the chip with a good level of approximation relative to the temperature accuracy of ±1°C demanded by PCR-LOC systems. This model structure is shown in Figure5 where the heater is defined by the circular strip at the bottom surface with k2 width and k1 distance to the vertical axis. It is worth noting that we only modelled the area of interest in the chip, and we have put aside the rest of the chip. Furthermore, we defined a comparison aperture that can be characterized by a fixed size rectangular area in cross-section view of the actual PCR-LOC microchip and our model. The comparison aperture contains the reaction chamber, and it is discussed in detail in the simulation results section. We chose this scenario to arrive at a model structure which is simple enough to find its analytical solution meanwhile it is accurate enough to characterize the actual thermal process. A simple estimate of the thermal capacity and radial thermal resistance of the PDMS layer shows that they are small compared to those of the glass layers. It is expected that treating the PDMS layer as if it has the same thermal properties of glass will make the systems dynamic behavior 9% too fast (because the radial resistance is predicted to be 10% too low and the thermal capacity is estimated to be about 1% too high). It is possible to correct for this by scaling the thermal parameters of the glass. The assumptions made will result in underestimating the vertical gradient, but that gradient cannot be controlled by control mechanism and must instead be dealt with at the design and simulation level (i.e., dynamic simulations of the present design ensured that the PCR chamber radius and height were sufficiently small that the vertical gradient could not lead to temperature variations greater than ±1°C). The FEM model by itself was verified by extensive experiments, wherein thermochromic liquid crystal (TLC) sensitive to different temperatures was deposited in the chamber to read the temperature within its volume in steady-state conditions.Figure 5
Geometrical representation of the PDE model. Taking advantage of the axial symmetry of the chip, the entire system is represented as a two-dimensional rectangular region. Similarly the heater is represented as a line instead of being a surface. This reduction of one dimension allows the model to be solved for the temperature in a three-dimensional space but at a small fraction of the time and computational resources required by a FEM simulation.


## 3. Heat Distribution Model

### 3.1. Problem Formulation

The partial differential equation describing temperature distribution in the chosen domain at any time is given by the heat conduction equation as follows [30]:
(1)K∇2u=ρ¯𝒞∂u∂t,
where u(ρ,φ,z,t) is the temperature of every point of the domain at any given time. K is thermal conductivity, ρ¯ is density, and 𝒞 is heat capacity. Since the chosen model features axial symmetry, the temperature in the cylinder does not depend on φ. This means that ∂u/∂φ=0=∂2u/∂φ2. Therefore φ can be dropped and the Laplace operator can be written in a 2D cross-section domain using ρ and z as follows:
(2)∇2u=∂2u∂ρ2+1ρ∂u∂ρ+∂u2∂z2.To avoid any confusion betweenρ and the notation used for density, we introduce here notation κ as follows:
(3)κ=Kρ¯𝒞
which is thermal diffusivity of the material in m2/s. Finally the partial differential equation giving the temperature u(ρ,z,t) at the half part cross-section of our cylindrical domain with ρ and z as its horizontal and vertical axes can be written as
(4)κ∇2u=∂u∂t.The definition of our domain implies thatρ and z are limited to 0≤ρ≤a and 0≤z≤b. The initial and boundary conditions are given by
(5)u(ρ,z,0)=u0(ρ,z),(6)u(a,z,t)=Tw,(7)uρ(0,z,t)=0,(8)uz(ρ,b,t)=0,(9)uz(ρ,0,t)=-f(ρ),
where u0(ρ,z) is the initial distribution of temperature. Tw is the temperature of the outside wall of the cylinder. The boundary condition (7) comes from the axisymmetrical cylindrical structure of the model. The boundary condition (8) is based on the assumption that the upper surface of the chip is isolated. f(ρ) describes the geometric distribution of the control action signal at the bottom surface, which is defined by the heater location. A general solution of (4) can be found by adding the solution of homogeneous and nonhomogeneous equations describing the steady-state and a transient temperature distribution, respectively:
(10)u(ρ,z,t)=u^(ρ,z,t)+u-(ρ,z).Substituting (10) in (4) we get
(11)1κu^t=u^ρρ+u-ρρ+1ρu^ρ+1ρu-ρ+u^zz+u-zz=(u^ρρ+1ρu^ρ+u^zz)+(u-ρρ+1ρu-ρ+u-zz).Without the loss of generality, we can assumeu-ρρ+(1/ρ)u-ρ+u-zz=0 as we defined u-(ρ,z) to describe the steady-state temperature. By applying (10) to the initial condition (5) and boundary conditions (6), (7), (8), and (9), we arrive at two PDEs:
(12)1κu^t=u^ρρ+1ρu^ρ+u^zzwith{u^(ρ,z,0)=u^0(ρ,z)u^(a,z,t)=0u^ρ(0,z,t)=0u^z(ρ,b,t)=0u^z(ρ,0,t)=0,(13)u-ρρ+1ρu-ρ+u-zz=0with{u-(a,z)=Twu-ρ(0,z)=0u-z(ρ,b)=0u-z(ρ,0)=-f(ρ).The first one is an initial value problem with homogeneous boundaries, while the second one is a nonhomogeneous boundary value problem.


### 3.2. Steady-State Part

To solve (13), we employ the method of separation of variables [31]:
(14)u-(ρ,z)=P-(ρ)Z-(z)+Tw,
where P- and  Z- are functions of ρ and z, respectively. The constant part is used to normalize the boundary conditions and to achieve the homogeneous boundary condition. Using (14) in (13) and dividing P- Z- we arrive to
(15)P-′′(ρ)P-(ρ)+1ρP-′(ρ)P-(ρ)=-Z-′′(z)Z-(z).The left side contains functions ofρ alone, while the right side contains functions of z alone. Since this equality must hold for all ρ and z in the given interval, the common value of the two sides must be a constant; say -λ-2 varying neither with ρ nor with z. The constant here was chosen as a negative value to avoid a trivial solution [31]. Now we have two ODEs for the two factor functions:
(16)ρP-′′+P-′+λ-2ρP-=0,Z-′′-λ-2Z-=0.The boundary conditions of (13) can also be stated in the product form (14) which results in the following:
(17)P-(a)=0,(18)P-′(0)=0,(19)Z-′(b)=0,(20)P-(ρ)Z-′(0)=-f(ρ).The solutions of (16) are given by [31]
(21)P-=c-1J0(λ-ρ)+c-2Y0(λ-ρ),(22)Z-=c-3coshλ-(b-z)+c-4sinhλ-(b-z),
where J0 is the Bessel function of the first kind of order 0 and Y0 is the Bessel function of the second kind of order 0. From the boundedness condition at ρ=0 stated by (18), we must have c-2=0 because Y0 is unbounded at zero. Thus the solution (21) becomes
(23)P-=c-1J0(λ-ρ).From the boundary condition (19) and the fact that
(24)Z-′(z)=-λ-(c-3sinhλ-(b-z)+c-4coshλ-(b-z))
we see that
(25)Z-′(b)=-λ-c-4=0
which implies c-4=0, so the solution (22) becomes
(26)Z-(z)=c-3coshλ-(b-z).The boundary condition (17) can be employed to determine λ-:
(27)P-(a)=c-1J0(λ-a)=0
which can be satisfied only if J0(λ-a)=0 so that
(28)aλ-=r1,r2,r3,…,rm,…,
where rm(m=1,2,3,…) is the mth positive root of J0(x)=0. Thus λ-’s can be expressed as
(29)λ-m=rma,m=1,2,3,….So a solution satisfying the boundary conditions is(30)u-(ρ,z)-Tw=P-m(ρ)Z-m(z)=A-J0(rmaρ)coshrma(b-z),
where m=1,2,3,… and A-=c-1c-3. The general solution of (13) can be obtained by employing the superposition principle [31]. By replacing A- with Am and summing up all eigenfunctions, it follows that
(31)u-(ρ,z)=Tw+∑m=1∞‍AmJ0(rmaρ)coshrma(b-z).From(32)∂u-∂z=∑m=1∞‍-AmrmaJ0(rmaρ)sinhrma(b-z)
and the last boundary condition (20) we can write
(33)∑m=1∞‍(Amrmasinhrmba)J0(rmaρ)=f(ρ).The right side is the Bessel series expansion off(ρ) with coefficients equal to (Am(rm/a)sinh(rmb/a)) [31], so we get
(34)Am=2armsinh(rmb/a)J12(rm)∫0a‍ρJ0(rmaρ)f(ρ)dρ.


### 3.3. Transient Part

We employ the separation of variables again to solve (12)
(35)u^(ρ,z,t)=P^(ρ)Z^(z)T^(t),
where P^, Z^, and T^ are functions of ρ, z, and t, respectively. Following the same procedure as the steady-state part and defining μ, ν, and λ^ as eigenvalues for P^, Z^, and T^, respectively, we arrive at the following ODEs:
(36)ρP^′′+P^′+μ2ρP^=0,(37)Z^′′-ν2Z^=0,(38)T^′+κλ^2T^=0,
where ν2=μ2-λ^2. By applying (35) into (12), the boundary conditions become
(39)P^(a)=0,(40)P^′(0)=0,(41)Z^′(b)=0,(42)Z^′(0)=0.The solutions of (36), (37), and (38) are given by
(43)P^=c^1J0(μρ)+c^2Y0(μρ),(44)Z^=c^3coshν(b-z)+c^4sinhν(b-z),(45)T^=c^5e-κλ^2t.From the boundedness condition atρ=0 expressed by (40), we must have c^2=0. Thus the solution of (43) becomes
(46)P^(ρ)=c^1J0(μρ).From the boundary condition (41), we have
(47)Z^′(b)=-νc^4=0
which implies c^4=0, so the solution (44) becomes
(48)Z^(z)=c^3coshν(b-z).From boundary condition (39) we can determine μ. Consider
(49)P^(a)=c^1J0(μa)=0
so we have
(50)μm=rma,
where rm(m=1,2,3,…) is the mth positive root of J0(x)=0. A solution of (36) satisfying the boundary conditions is
(51)P^m(ρ)=c^1J0(rmaρ),
where m=1,2,3,…. From the boundary condition (42) we can determine ν. We have
(52)Z^′(0)=-νc^3sinhνb=0
which can be satisfied only if sinhνb=0 or
(53)ν=kπibk=0,1,2,….
Using (53) in (48) and the fact that
(54)coshix=cosx
the solution (44) becomes
(55)Z^k(z)=c^3coskπb(b-z).Now from (53) and (50) it follows that
(56)λ^2=μ2-ν2=(rma)2-(kπib)2=rm2a2+k2π2b2.So a solution satisfying all boundary conditions is given by(57)u^(ρ,z,t)=A^e-κ(rm2/a2+k2π2/b2)tJ0(rmaρ)coskπb(b-z),
where A^=c^1c^3c^5, k=0,1,2,3,…, and m=1,2,3,…. Replacing A^ with Akm and summing up k and m we obtain the solution by the superposition principle
(58)u^(ρ,z,t)=∑k=0∞‍∑m=1∞‍Akme-κ(rm2/a2+k2π2/b2)t×J0(rmaρ)coskπb(b-z).From the initial conditionu^(ρ,z,0)=u^0(ρ,z) we have
(59){∑k=0∞‍∑m=1∞‍AkmJ0(rmaρ)coskπb(b-z)}=u0(ρ,z)-u-(ρ,z).This can be written as(60)∑k=0∞‍∑m=1∞‍AkmJ0(rmaρ)coskπb(b-z)=(u0(ρ,z)-Tw)+∑m=1∞‍-AmJ0(rmaρ)coshrma(b-z),
where Am is defined by (34). To simplify the equation, without the loss of accuracy, we consider that u0(ρ,z)=Tw which is reasonable because we start the experiment when the system is in a uniform temperature distribution. We can write (60) in this form
(61)∑m=1∞‍{∑k=0∞‍Akmcoskπb(b-z)}J0(rmaρ)=∑m=1∞‍{-Amcoshrma(b-z)}J0(rmaρ).Clearly because of orthogonality of the Bessel functions, these coefficients have to be equal; that is,(62)∑k=0∞‍Akmcoskπb(b-z)=-Amcoshrma(b-z).Written in standard form, we have(63)A0m+∑k=1∞‍Akmcoskπb(b-z)=-Amcoshrma(b-z).Recalling theFourier cosine series expansion [31], the left side of (63) is the Fourier cosine series expansion of the right side with respect to b-z. Thus substituting x=b-z and replacing A0m and Akm with a0/2 and ak, respectively, we have
(64)a02+∑k=1∞‍akcoskπbx=-Amcoshrmax.Soak for k=0,1,2,… are given by
(65)ak=2b∫0b‍(-Amcoshrmax)coskπbxdz.Equation (65) can be written as
(66)ak=-2bAm∫0b‍coshrmaxcoskπbxdz
and solving the integral we have
(67)ak=-2bAmaba2k2π2+b2rm2×(+brmcos(kπ)sinh(rmab)akπcosh(rmab)sin(kπ)+brmcos(kπ)sinh(rmab)).Sincek takes only positive integers, we can simplify (67) and arrive to
(68)ak=2Amrmab(-1)k+1rm2/a2+k2π2/b2sinh(rmab).Thus we have(69)A0m=a02=-Amabrmsinh(rmab)
and Akm is defined by (68) as
(70)Akm=ak=2Amrmab(-1)k+1rm2/a2+k2π2/b2sinh(rmab),
where k=1,2,…. Hence, it follows that the complete solution to the problem is given by (10), (31), and (58), where their coefficients are defined by (34), (69), and (70).


### 3.4. Stability Analysis

Stability analysis of the calculated model can be easily performed by studying the transient part of the general solution. A quick visit of transient equation (58) reveals that the infinite poles or eigenvalues of the system are given by the coefficients of the exponential term in infinite series as follows:
(71)SystemPoles:-κ(rm2a2+k2π2b2),
where k=0,1,2,… and m=1,2,3,…. Because thermal diffusivity, κ, takes only positive real values for different materials, all the poles are negative real poles, and therefore the system is stable. Calculating the slowest poles of the system using (71) and the parameters given in Table 1, we arrive to the following results:
(72)SystemPoles=-0.14,-0.76,-1.5,-1.86,-2.11,….Table 1
Parameter values of PCR-LOC chip model.
Parameter
Symbol
Value
Unit
Inner radius of heater
k
1
2.3
mm
Heater trace width
k
2
200
μm
Chip radius
a
4.94
mm
Chip height
b
2.10
mm
Heater power
P
0.77
W
Thermal conductivity of Borofloat glass
K
1.11
W/(m·K)
Density of Borofloat glass
ρ
-
2200
Kg/m3
Thermal capacity of Borofloat glass
𝒞
830
J/K


## 3.1. Problem Formulation

The partial differential equation describing temperature distribution in the chosen domain at any time is given by the heat conduction equation as follows [30]:
(1)K∇2u=ρ¯𝒞∂u∂t,
where u(ρ,φ,z,t) is the temperature of every point of the domain at any given time. K is thermal conductivity, ρ¯ is density, and 𝒞 is heat capacity. Since the chosen model features axial symmetry, the temperature in the cylinder does not depend on φ. This means that ∂u/∂φ=0=∂2u/∂φ2. Therefore φ can be dropped and the Laplace operator can be written in a 2D cross-section domain using ρ and z as follows:
(2)∇2u=∂2u∂ρ2+1ρ∂u∂ρ+∂u2∂z2.To avoid any confusion betweenρ and the notation used for density, we introduce here notation κ as follows:
(3)κ=Kρ¯𝒞
which is thermal diffusivity of the material in m2/s. Finally the partial differential equation giving the temperature u(ρ,z,t) at the half part cross-section of our cylindrical domain with ρ and z as its horizontal and vertical axes can be written as
(4)κ∇2u=∂u∂t.The definition of our domain implies thatρ and z are limited to 0≤ρ≤a and 0≤z≤b. The initial and boundary conditions are given by
(5)u(ρ,z,0)=u0(ρ,z),(6)u(a,z,t)=Tw,(7)uρ(0,z,t)=0,(8)uz(ρ,b,t)=0,(9)uz(ρ,0,t)=-f(ρ),
where u0(ρ,z) is the initial distribution of temperature. Tw is the temperature of the outside wall of the cylinder. The boundary condition (7) comes from the axisymmetrical cylindrical structure of the model. The boundary condition (8) is based on the assumption that the upper surface of the chip is isolated. f(ρ) describes the geometric distribution of the control action signal at the bottom surface, which is defined by the heater location. A general solution of (4) can be found by adding the solution of homogeneous and nonhomogeneous equations describing the steady-state and a transient temperature distribution, respectively:
(10)u(ρ,z,t)=u^(ρ,z,t)+u-(ρ,z).Substituting (10) in (4) we get
(11)1κu^t=u^ρρ+u-ρρ+1ρu^ρ+1ρu-ρ+u^zz+u-zz=(u^ρρ+1ρu^ρ+u^zz)+(u-ρρ+1ρu-ρ+u-zz).Without the loss of generality, we can assumeu-ρρ+(1/ρ)u-ρ+u-zz=0 as we defined u-(ρ,z) to describe the steady-state temperature. By applying (10) to the initial condition (5) and boundary conditions (6), (7), (8), and (9), we arrive at two PDEs:
(12)1κu^t=u^ρρ+1ρu^ρ+u^zzwith{u^(ρ,z,0)=u^0(ρ,z)u^(a,z,t)=0u^ρ(0,z,t)=0u^z(ρ,b,t)=0u^z(ρ,0,t)=0,(13)u-ρρ+1ρu-ρ+u-zz=0with{u-(a,z)=Twu-ρ(0,z)=0u-z(ρ,b)=0u-z(ρ,0)=-f(ρ).The first one is an initial value problem with homogeneous boundaries, while the second one is a nonhomogeneous boundary value problem.


## 3.2. Steady-State Part

To solve (13), we employ the method of separation of variables [31]:
(14)u-(ρ,z)=P-(ρ)Z-(z)+Tw,
where P- and  Z- are functions of ρ and z, respectively. The constant part is used to normalize the boundary conditions and to achieve the homogeneous boundary condition. Using (14) in (13) and dividing P- Z- we arrive to
(15)P-′′(ρ)P-(ρ)+1ρP-′(ρ)P-(ρ)=-Z-′′(z)Z-(z).The left side contains functions ofρ alone, while the right side contains functions of z alone. Since this equality must hold for all ρ and z in the given interval, the common value of the two sides must be a constant; say -λ-2 varying neither with ρ nor with z. The constant here was chosen as a negative value to avoid a trivial solution [31]. Now we have two ODEs for the two factor functions:
(16)ρP-′′+P-′+λ-2ρP-=0,Z-′′-λ-2Z-=0.The boundary conditions of (13) can also be stated in the product form (14) which results in the following:
(17)P-(a)=0,(18)P-′(0)=0,(19)Z-′(b)=0,(20)P-(ρ)Z-′(0)=-f(ρ).The solutions of (16) are given by [31]
(21)P-=c-1J0(λ-ρ)+c-2Y0(λ-ρ),(22)Z-=c-3coshλ-(b-z)+c-4sinhλ-(b-z),
where J0 is the Bessel function of the first kind of order 0 and Y0 is the Bessel function of the second kind of order 0. From the boundedness condition at ρ=0 stated by (18), we must have c-2=0 because Y0 is unbounded at zero. Thus the solution (21) becomes
(23)P-=c-1J0(λ-ρ).From the boundary condition (19) and the fact that
(24)Z-′(z)=-λ-(c-3sinhλ-(b-z)+c-4coshλ-(b-z))
we see that
(25)Z-′(b)=-λ-c-4=0
which implies c-4=0, so the solution (22) becomes
(26)Z-(z)=c-3coshλ-(b-z).The boundary condition (17) can be employed to determine λ-:
(27)P-(a)=c-1J0(λ-a)=0
which can be satisfied only if J0(λ-a)=0 so that
(28)aλ-=r1,r2,r3,…,rm,…,
where rm(m=1,2,3,…) is the mth positive root of J0(x)=0. Thus λ-’s can be expressed as
(29)λ-m=rma,m=1,2,3,….So a solution satisfying the boundary conditions is(30)u-(ρ,z)-Tw=P-m(ρ)Z-m(z)=A-J0(rmaρ)coshrma(b-z),
where m=1,2,3,… and A-=c-1c-3. The general solution of (13) can be obtained by employing the superposition principle [31]. By replacing A- with Am and summing up all eigenfunctions, it follows that
(31)u-(ρ,z)=Tw+∑m=1∞‍AmJ0(rmaρ)coshrma(b-z).From(32)∂u-∂z=∑m=1∞‍-AmrmaJ0(rmaρ)sinhrma(b-z)
and the last boundary condition (20) we can write
(33)∑m=1∞‍(Amrmasinhrmba)J0(rmaρ)=f(ρ).The right side is the Bessel series expansion off(ρ) with coefficients equal to (Am(rm/a)sinh(rmb/a)) [31], so we get
(34)Am=2armsinh(rmb/a)J12(rm)∫0a‍ρJ0(rmaρ)f(ρ)dρ.


## 3.3. Transient Part

We employ the separation of variables again to solve (12)
(35)u^(ρ,z,t)=P^(ρ)Z^(z)T^(t),
where P^, Z^, and T^ are functions of ρ, z, and t, respectively. Following the same procedure as the steady-state part and defining μ, ν, and λ^ as eigenvalues for P^, Z^, and T^, respectively, we arrive at the following ODEs:
(36)ρP^′′+P^′+μ2ρP^=0,(37)Z^′′-ν2Z^=0,(38)T^′+κλ^2T^=0,
where ν2=μ2-λ^2. By applying (35) into (12), the boundary conditions become
(39)P^(a)=0,(40)P^′(0)=0,(41)Z^′(b)=0,(42)Z^′(0)=0.The solutions of (36), (37), and (38) are given by
(43)P^=c^1J0(μρ)+c^2Y0(μρ),(44)Z^=c^3coshν(b-z)+c^4sinhν(b-z),(45)T^=c^5e-κλ^2t.From the boundedness condition atρ=0 expressed by (40), we must have c^2=0. Thus the solution of (43) becomes
(46)P^(ρ)=c^1J0(μρ).From the boundary condition (41), we have
(47)Z^′(b)=-νc^4=0
which implies c^4=0, so the solution (44) becomes
(48)Z^(z)=c^3coshν(b-z).From boundary condition (39) we can determine μ. Consider
(49)P^(a)=c^1J0(μa)=0
so we have
(50)μm=rma,
where rm(m=1,2,3,…) is the mth positive root of J0(x)=0. A solution of (36) satisfying the boundary conditions is
(51)P^m(ρ)=c^1J0(rmaρ),
where m=1,2,3,…. From the boundary condition (42) we can determine ν. We have
(52)Z^′(0)=-νc^3sinhνb=0
which can be satisfied only if sinhνb=0 or
(53)ν=kπibk=0,1,2,….
Using (53) in (48) and the fact that
(54)coshix=cosx
the solution (44) becomes
(55)Z^k(z)=c^3coskπb(b-z).Now from (53) and (50) it follows that
(56)λ^2=μ2-ν2=(rma)2-(kπib)2=rm2a2+k2π2b2.So a solution satisfying all boundary conditions is given by(57)u^(ρ,z,t)=A^e-κ(rm2/a2+k2π2/b2)tJ0(rmaρ)coskπb(b-z),
where A^=c^1c^3c^5, k=0,1,2,3,…, and m=1,2,3,…. Replacing A^ with Akm and summing up k and m we obtain the solution by the superposition principle
(58)u^(ρ,z,t)=∑k=0∞‍∑m=1∞‍Akme-κ(rm2/a2+k2π2/b2)t×J0(rmaρ)coskπb(b-z).From the initial conditionu^(ρ,z,0)=u^0(ρ,z) we have
(59){∑k=0∞‍∑m=1∞‍AkmJ0(rmaρ)coskπb(b-z)}=u0(ρ,z)-u-(ρ,z).This can be written as(60)∑k=0∞‍∑m=1∞‍AkmJ0(rmaρ)coskπb(b-z)=(u0(ρ,z)-Tw)+∑m=1∞‍-AmJ0(rmaρ)coshrma(b-z),
where Am is defined by (34). To simplify the equation, without the loss of accuracy, we consider that u0(ρ,z)=Tw which is reasonable because we start the experiment when the system is in a uniform temperature distribution. We can write (60) in this form
(61)∑m=1∞‍{∑k=0∞‍Akmcoskπb(b-z)}J0(rmaρ)=∑m=1∞‍{-Amcoshrma(b-z)}J0(rmaρ).Clearly because of orthogonality of the Bessel functions, these coefficients have to be equal; that is,(62)∑k=0∞‍Akmcoskπb(b-z)=-Amcoshrma(b-z).Written in standard form, we have(63)A0m+∑k=1∞‍Akmcoskπb(b-z)=-Amcoshrma(b-z).Recalling theFourier cosine series expansion [31], the left side of (63) is the Fourier cosine series expansion of the right side with respect to b-z. Thus substituting x=b-z and replacing A0m and Akm with a0/2 and ak, respectively, we have
(64)a02+∑k=1∞‍akcoskπbx=-Amcoshrmax.Soak for k=0,1,2,… are given by
(65)ak=2b∫0b‍(-Amcoshrmax)coskπbxdz.Equation (65) can be written as
(66)ak=-2bAm∫0b‍coshrmaxcoskπbxdz
and solving the integral we have
(67)ak=-2bAmaba2k2π2+b2rm2×(+brmcos(kπ)sinh(rmab)akπcosh(rmab)sin(kπ)+brmcos(kπ)sinh(rmab)).Sincek takes only positive integers, we can simplify (67) and arrive to
(68)ak=2Amrmab(-1)k+1rm2/a2+k2π2/b2sinh(rmab).Thus we have(69)A0m=a02=-Amabrmsinh(rmab)
and Akm is defined by (68) as
(70)Akm=ak=2Amrmab(-1)k+1rm2/a2+k2π2/b2sinh(rmab),
where k=1,2,…. Hence, it follows that the complete solution to the problem is given by (10), (31), and (58), where their coefficients are defined by (34), (69), and (70).


## 3.4. Stability Analysis

Stability analysis of the calculated model can be easily performed by studying the transient part of the general solution. A quick visit of transient equation (58) reveals that the infinite poles or eigenvalues of the system are given by the coefficients of the exponential term in infinite series as follows:
(71)SystemPoles:-κ(rm2a2+k2π2b2),
where k=0,1,2,… and m=1,2,3,…. Because thermal diffusivity, κ, takes only positive real values for different materials, all the poles are negative real poles, and therefore the system is stable. Calculating the slowest poles of the system using (71) and the parameters given in Table 1, we arrive to the following results:
(72)SystemPoles=-0.14,-0.76,-1.5,-1.86,-2.11,….Table 1
Parameter values of PCR-LOC chip model.
Parameter
Symbol
Value
Unit
Inner radius of heater
k
1
2.3
mm
Heater trace width
k
2
200
μm
Chip radius
a
4.94
mm
Chip height
b
2.10
mm
Heater power
P
0.77
W
Thermal conductivity of Borofloat glass
K
1.11
W/(m·K)
Density of Borofloat glass
ρ
-
2200
Kg/m3
Thermal capacity of Borofloat glass
𝒞
830
J/K


## 4. Simulation Results

The parameters with the values listed in Table1 are used in simulation. We assume that Tw=22∘C and that an input power of one Watt is being applied to the heater. The input to the system is given by the heat flux transferred from the heater to the system. This heat flux is considered in definition of f(ρ) which by itself determines the Am coefficients. Definition of f(ρ) is given based on the chosen model configuration by the following equation:
(73)f(ρ)={qKifk1<ρ<(k1+k2)0otherwise,
where 0<k1,k2≤a, 0<k1+k2≤a, q is heat flux from the heater to chip, and K is the thermal conductivity of the chip. For discontinuities we have f(k1)=q/2K=f(k1+k2). Now for Am we can write
(74)Am=2qaKrmsinh(rmb/a)J12(rm)×∫k1k1+k2‍ρJ0(rmaρ)dρ.For heat fluxq, we have
(75)q=PA,
where P is the amount of the power transferred from the heater to the chip and A is the area of the heater which can be calculated as follows:
(76)A=π(k1+k2)2-πk12=πk2(2k1+k2).Recalling the recurrence relations for the Bessel function [31], we have
(77)ddx[xnJn(x)]=xnJn-1(x)(n=1,2,…).Forn=1, (77) can be written as follows:
(78)∫‍xJ0(x)dx=xJ1(x)+c.Using a change of variable,x=(rm/a)ρ, for the integral part of Am we can write
(79)∫k1k1+k2‍ρJ0(rmaρ)dρ=(arm)2∫(rm/a)k1(rm/a)(k1+k2)‍xJ0(x)dx.Using (78) in (79) we get
(80)∫k1k1+k2‍ρJ0(rmaρ)dρ=(arm)2[xJ1(x)]|(arm)2(rm/a)k1(rm/a)(k1+k2)=arm((k1+k2)J1(rma(k1+k2))-k1J1(rmak1)).Taking into account (75) and (80) in (74), we have
(81)Am=2P((k1+k2)J1(rma(k1+k2))-k1J1(rmak1))×(πk2(2k1+k2)Krm2sinh(rmba)J12(rm))-1.Now we have all the required equations to calculate the temperature distribution using our analytic solution.The next important step before start running simulations is how to compensate the error caused due to the simplification of the actual system into a single-domain model with Borofloat glass as the only material involved. Substituting the PDMS layer with a glass layer could result in a significant amount of error and mismatch between model and the actual system. To overcome this issue, we constructed an optimization problem to optimize our model parameters arriving to a best fit between single domain analytic model and multidomain actual system. Chip radius, chip thickness, and input power in the mathematical model are selected as the optimization degrees of freedom. We focus our attention on the two regions that have the highest importance. The first region is a rectangular area which starts from the left bottom side of the chamber and continues up to the top of the constant temperature wall and is all made from glass in the analytical model. We called this region comparison aperture. The second region which is a subset of comparison aperture and has higher importance is reaction chamber. To get quantitative measure of the error in both the chamber and the selected aperture, we define the following parameters in both regions: absolute average error, maximum absolute error, and uniformity of the error. The definitions of the first two measures are obvious from their naming, and the last measure is defined as follows:(82)Uap≐max(eap(ρ,z))-min(eap(ρ,z)),ρ,z∈𝒟1,Uch≐max(ech(ρ,z))-min(ech(ρ,z)),ρ,z∈𝒟2,
where eap is error at selected aperture indicated by 𝒟1 and ech is error at chamber region indicated by 𝒟2. Next we define our weighted-sum objective function including 𝕃2 norm of the error at chamber location, uniformity measure at chamber location, and uniformity measure at selected aperture as follows:
(83)J=w11A2∫A2‍ech2dA2+w2Uch+w3Uap,
where w1, w2, and w3 are weight coefficients and A2 is area of chamber region. By selecting proper weights, optimization results in a maximum absolute error of 0.14°C in chamber temperature, which is less than 0.03 of the error when compared to the results of a nonoptimized model (a=3.5 mm, b=1.354 mm, and P=1 W), satisfy the acceptable error criteria for reliable PCR operation. Error distribution in whole aperture area is depicted in Figure 8. It is clear that, although the maximum absolute error reaches 3.6°C in the area close to the heater, the error in the chamber area is kept lower than 0.14°C.The resulted steady-state temperature distributions based on optimized analytical solution and the FEM simulation are shown in Figures7 and 6, respectively. Comparison aperture is indicated by dash lines. The solutions of the two results show an excellent degree of correlation, verifying accuracy of the steady-state solution.Figure 6
Steady-state temperature distribution of the cross-section of the chip calculated by the FEM. The axis of symmetry is along the left vertical side. The lines bound four domains assigned to the thermals conductivities of glass, PDMS, and water. The glass plates and PDMS layer are 1.1 mm thick and 254 m thick, respectively. The thin rectangle on the left defines the PCR chamber. The heater is located at the point of the highest temperature in red.Figure 7
Steady-state temperature distribution calculated by the PDE model. The geometry of the model is a small portion of the geometry of real system and the FEM model. The area shown corresponds to the upper left part of the graph in Figure6. The PDE model represents the system with a single material. The dashed lines indicate the location of the chamber and PDMS layer but do not define domains of different materials.Figure 8
Error distribution in selected aperture after optimization. Obviously, the maximum absolute error within the selected aperture is kept under 3.6°C. It reaches to its maximum in the area close to the heater. The error in the chamber area is well kept under 0.14°C. The chamber area is a rectangular area with 0.9 mm thickness and 1.5 mm diameter just on the left bottom corner of the selected aperture.Dynamic simulation of temperature evolution inside the reaction chamber is the next important result that we would like to study about our analytic model. Precise FEM simulations disclosed that lower center point of the chamber and upper edge of the chamber have the highest and the lowest temperatures inside the chamber, respectively (SW and NE edges of the chamber in Figure5). The results follow the fact that upper center point of the chamber (NW edge of the chamber in Figure 5) has a temperature close to mean temperature of the chamber. We selected upper center point of the chamber to study the dynamic evolution of temperature in our model and in actual system. An input power which consisted of three predefined values that are required to take the chamber temperature to denaturation, annealing, and extension temperatures is applied to both models. Both systems were powered up after being settled in ambient temperature and kept in each power level for 60 seconds. The experiment ran over a whole PCR cycle in a 240-second interval in open-loop condition. The results are depicted in Figure 9. The optimization was performed for the lowest temperature in PCR cycle. However the error at the highest temperature is still under 1°C. Model dynamics nicely follows the FEM results, and as we expected, the transient response is about 1% faster than that of the actual system, due to the difference between the thermal properties of PDMS and Borofloat glass.Figure 9
Transient response of the PDE model to a predefined input power trajectory in comparison to the FEM results. The input power represents a complete PCR cycle. There is less than 1% error in dynamic response of the derived PDE model and the complete FEM simulation.


## 5. Conclusion

This paper presents the derivation and the solution of the heat distribution equation of a microfluidic chip designed and manufactured at the Applied Miniaturization Laboratory at the University of Alberta, Canada. The partial differential equation featuring three-dimensional domain was used to describe the heat distribution of the chip during a DNA amplification process. The axisymmetric architecture of the chip helped reduce the domain to two-dimensional cylindrical form. The PDE equation describing the heat distribution of the chip was divided into steady-state and transient equations and then solved explicitly. Precise FEM simulations in addition to laboratory measurements were used to verify the developed model. The modeling approach is considered to be appropriate for similar thermal processes with collocated sensors and actuators. Future work is being carried out toward developing closed-loop tracking control based on boundary control techniques.A dynamic simulation of the system was made using the developed model to confirm that the scaling of the glass thermal properties could be adjusted to match the temperature predicted by the FEM simulations at all points in the PCR chamber with less than 1°C error at all times during a PCR temperature cycle.



---

*Source: 767853-2013-10-09.xml*