Temporal response of an initially deflected PDMS channel

Microfluidics is used in a wide array of applications including photonics, catalysis, optoelectronics, particle synthesis, bio-diagnosis and tissue engineering. Poly(dimethyl) siloxane (PDMS) is one of the most widely used materials for the fabrication of such devices. In many applications, applied flow gives rise to pressure variations along the length of the channel and a commensurate spatially varying channel height. We have previously performed a scaling analysis of the relaxation time of an initially deflected PDMS channel (Dendukuri et al 2007 Lab Chip 7 818–28). The prior study did not consider the dynamic evolution of the initially deflected channel and, being a scaling theory, required a fitting parameter for a priori predictions. In this paper, we analyze the dynamic behavior of a retracting PDMS wall after the removal of an external stress. For small deflections, the problem lends itself to a regular perturbation analysis that is analytically solvable at zeroth order. We find that the zeroth-order solution is very close to the numerically solved full solution for most situations of experimental interest. We further demonstrate that our model captures both qualitative and quantitative trends seen in the prior (Dendukuri et al 2007) and newly performed experiments.


Introduction
The field of microfluidics has grown exponentially over the past decade, driven by the realization that miniaturization of processes leads to distinct advantages in diverse fields such as biological separations, synthetic chemistry and combinatorial screening. While devices on the micron scale can be constructed using 'hard' materials like silicon, the rapid growth of microfluidics has been driven by replica molding technologies which enable the economical and fast fabrication of disposable plastic devices. Poly(dimethyl) siloxane (PDMS) has become one of the most widely used materials in such technologies because of its ease in molding, high fidelity, low surface energy, transparency, porosity and conductivity, among other qualities.
Among the properties of PDMS which may or may not be useful, depending on the intended application, is its low bulk modulus (in the range of 0.5-4 MPa) [2,3] which leads to deformation under external pressure. This deformation has been cleverly exploited to synthesize valves [4] and microlenses [5]. Flow in deforming channels has also been used to model cardiovascular flow behavior [6]- [8]. PDMS deformation is, however, undesirable in applications where wall rigidity is essential, causing sagging of channels [9], increasing Taylor dispersion and diverging from predictions of flow profiles and mass transfer that have been made for rectangular cross-sections. In pressure-driven flows, low aspect ratio channels will have a tendency to bulge and thus any microfluidic application that requires a precise knowledge of the flow profile, pressure and channel geometry [10] will be compromised if the actual channel deformation cannot be quantified and taken into account. Some common microfluidic applications which are dependent upon the channel deformation are shear-flow assays (cellular or protein) [11,12] and surface plasmon resonance (SPR) sensing technology [13,14]. Deformation of PDMS channels can also affect fluorescence intensity signals, which depend on the absorption path length across the channel, thus biasing measurements if deformation is not taken into account. The synthesis of monodisperse particles using microfluidic devices [15] also depends on an accurate prediction of the deformation profile of PDMS.
In a previous work [1], we have analyzed the dynamic behavior of a PDMS channel under external pressure by using a scaling argument. While this prior work gives insight into how the relaxation of the system scales with fluid properties and channel dimensions, it does not give a priori predictions for absolute values for the relaxation times nor does it provide information about the local temporal response of the system. For example, the scaling theory does not show that transient flow reversal occurs in a portion of the channel. Moreover, previous work on the deformation of PDMS channels has depended on simulations to make predictions of channel deformations, pressure distribution and flow profiles. The purpose of the current work is to develop a more comprehensive analysis of the dynamics of a deforming PDMS microfluidic channel and use a perturbation analysis to provide analytic expressions to describe the dynamic response. The theory will be compared to existing and new experimental data.

Theory
We begin with the equation for channel height as a function of time and space which has been derived earlier [1] In the above equation, x represents the direction of the primary flow, t is the time, h is the height of the channel and is a function of position and time, E is the Young's modulus of the channel walls (here PDMS), W is the width of the channel and µ is the viscosity of the Newtonian fluid flowing through the channel. As we might recall, the height represents an average height across the width of the channel at any cross-sectional interface. Further, for this creeping flow (Re 1) of fluid through a thin channel, the flow is mainly in the x-direction and the lubrication approximation is applied to the Navier-Stokes equation. The derivation also assumes the curvature of the PDMS top wall to be small in relation to the length of the device from which it follows that the velocity of the retracting membrane is only in the negative z-direction. Lastly, this equation is valid for small deformations in the range where the linear stress-strain relation of PDMS is valid, i.e. ε < 0.25 (see section 3.3). ε is the maximum deformation scaled by the relaxed channel height. An illustration of the initial deformation on application of the pressure pulse is shown in figure 1(a). We now rewrite this equation in terms of the deformation of the channel ( h(x, t)) using h(x, t) = H + h(x, t), where H is the undeformed height of the channel and h(x, t) is the deformation of the channel The dependent variable in the problem is now the extent of deformation of the channel, h(x, t). The initial condition for the problem is the steady state value of channel deformation and can be calculated from earlier work [3] h( Here P in is the gauge pressure (pressure above the atmospheric pressure) at the inlet of the device where the channel deformation is initially at its largest value (H + h max , where h max = P in W/E) and L is the length of the channel. The boundary conditions are given by Nondimensionalizing equations (2)-(4) using we arrive at with the boundary and initial conditions given by where ε is defined as and is the maximum deformation in the channel scaled by the channel height. For deformations which are small compared to the channel height, ε 1, we can solve equation (5) using a regular perturbation approach by expanding The zeroth-order problem governing θ 0 is obtained by setting ε = 0 in the original equations to give The solution to equation (11) with boundary conditions specified in equation (12) is obtained using finite Fourier transforms (FFT) [16] and is given by The higher order perturbation equations cannot be solved analytically and thus we stop here with the zeroth-order solution and will present a numerical solution to the full problem later in the paper. The nondimensionalized pressure (℘ = P/P in ) is related to the deformation of the channel by [1] ∂℘ ∂η = ∂θ ∂η .
The zeroth-order solution for the nondimensionalized pressure is given by From our previous work [1], the velocity in the x-direction is given by The velocity at z=H /2 is given by Nondimensionalizing equation (17) using 6 we obtain At small deformations, the second term in the bracket in equation (19) can be neglected and the nondimensionalized centerline velocity using equations (13) and (19) is given by Now, using equations (12) and (20), the initial steady state nondimensionalized velocity is 1 for the zeroth-order solution. For the velocity predictions from the full equation (equation (5)), we concentrate on the initial centerline velocity at the exit, which was derived in an earlier work [1] and can be nondimensionalized to obtain .
The above expression can also be derived using equations (7), (8) and (19). We have a squeeze flow problem wherein, at all η < 0.5, the velocity changes direction from initially flowing toward the exit (positive) to toward the inlet (negative). Therefore, the response time (κ) is defined as the time required for the velocity to reduce to 0.01 (1% of the initial steady state value), for η > 0.5 and −0.01, for η < 0.5. The nondimensionalized response time (λ) is given by For large values of η and τ , n = 1 is the dominating term in the FFT. Using the expression for V x (H/2) in equation (20), the time required for the deformation to reduce to c% of its initial value as η → 1 is given by Finally, the position of the stagnation point (η s ) can be found by solving ∂θ/∂η(η s ) = 0. The stagnation point is the position along the length in the channel where the velocity of the fluid in the x direction is zero. For experimental validation, the nondimensionalized time (α) for the stagnation point to reach any η < 0.5 is used. α at any η is defined as the nondimensionalized time at which ∂θ/∂η = 0 for predictions by the zeroth-order solution. For experimental observations, however, α at any η is defined as the nondimensionalized time at which ∂θ/∂η changes sign from negative to positive values.

Zeroth order and full solution
The zeroth-order problem (equation (11)) is solved analytically using FFT. The solution to the zeroth-order problem given in equation (13) and for the velocity given in equation (20) was truncated at five terms (the contribution of the sixth term is less than 1% of the sum of five terms). We solved the full problem (equation (5)) numerically using the PDE solver in MATLAB, to obtain θ as a function of η and τ . This solution for θ was used to obtain the velocity from the full solution.  table were used  whereas for TMPTA and TMPTA

Microfluidic devices
Microfluidic devices were made by pouring PDMS (Sylgard 184, Dow Corning) on a silicon wafer containing positive-relief channels patterned in SU-8 photoresist (Microchem) using polymer and curing agent in a ratio 10 : 1. The thickness of the PDMS devices was always maintained to be 5 mm. Straight channels with a rectangular cross-section and variable widths, lengths and heights were used as reported in table 1.

Materials
In all experiments involving PEG-DA, a 1% solution of 1.57 µm PMMA beads in poly(ethylene glycol) (400) diacrylate (PEG-DA, Sigma Aldrich) was used. In experiments involving TMPTA and TMPTA esters, a 1% solution of 1.5 µm Fluoresbrite carboxylate beads (Polysciences) in tri(methylol propane) triacrylate (TMPTA, Polysciences) or ethoxylated tri(methylol propane) triacryalte esters (TMPTA esters, Sartomer) was used. PEG-DA, TMPTA and TMPTA esters have been reported by the manufacturers to have viscosities of 56, 106 and 230 cP, respectively, at 25 • C. The linear stress-strain regime of PDMS was found to be valid for ε < 0.25 using dynamic mechanical analysis carried out using a DMA Q800 from TA Instruments. The elastic modulus of PDMS, in the linear stress-strain regime was found to be 2 MPa by the same experiment.

Bead tracking
A dilute solution of beads was used to track the flow in the device. A repeating square wave pressure profile comprising two parts-1 s of applied pressure and 4 s of turning the pressure to atmosphere was applied to the system. The beads in the mid-plane of the channel (z = H/2) were followed with a 20× or 40× microscope objective (Zeiss) with an optivar setting of 8 2.5× leading to effective magnifications of 50× and 100× respectively. Movies of the beads were recorded using a video tape recorder (DSR-25, Sony) using a CCD camera that captured images at the rate of 30 frames s −1 using an exposure of 1/1500 s. The frame-to-frame position of beads was measured as they moved across the screen using macros set up in NIH image. The central difference approximation was used to calculate the velocity. All bead tracking experiments were performed at a pressure of 3 psi. The bead tracking experiments at the exit of the channel were performed using PEG-DA, TMPTA and TMPTA esters for the channel dimensions reported in table 1. The bead tracking experiments at different positions along the channel were performed using PEG-DA for channels having a height of 10 or 20 µm and length and width of 1 cm and 200 µm, respectively. Wall effects were avoided by choosing beads close to the center of the channel (y = W/2).

Results and discussion
The PDMS channel deforms on the application of a pressure pulse as shown in figure 1(a). The deflected PDMS channel then relaxes back to its undeformed state after the pressure pulse has been turned off as shown in figure 1(b). In microfluidic devices fluid flows are commonly driven using syringe pumps. These devices compress the fluid in the tubing external to the microfluidic device, causing transients that can be several minutes or longer for micron-scale systems [17]. A rapid dynamic response is desired [17,18]. Therefore, we use air driven flows in our experiments using the SFL setup, which lead to finite transients associated with the deformation of the PDMS device. The channel dimensions used for this work lead to deformations up to ∼ 5 µm.
Measuring these deformations accurately using our SFL setup is difficult. We used bead tracking to measure local velocity in the channel and calculate the response times. The local velocity in the channel is dependent on the wall deformation and also is the more critical measurable in many applications. Therefore, we use the response time at a given position along the channel to compare the zeroth-order solution with experimentally observed values. We begin by numerically solving for the deformation profile using the full equation (equation (5)) presented in figure 2. It can be seen from figure 2 that at initial times, the relaxation is rapid followed by a more gradual change. Also, for values of τ > 0.2, the deformation is almost symmetric about the center of the channel. Past studies on flow in deforming microchannels have depended on solving complicated equations numerically [3], [6]- [8]. The ability to make accurate predictions using an analytical solution would be useful. Hence, we studied the applicability of the analytical zeroth-order solution by studying the effect of increasing ε on: (i) the response times at the exit of the channel as shown in figure 3(a) and (ii) the deformation profile of a PDMS channel for an early time (τ = 0.1) as shown in figure 3(b). It can be seen that the deviation of the full solution from the zerothorder solution increases systematically on increasing ε. The maximum relative deviation of the zeroth-order solution from the full solution for response times at the exit of the channel is 5.8% ( figure 3(a)). For profile predictions the maximum relative deviation of the zeroth-order solution from the full solution is as high as 40% for values of η 1 ( figure 3(b)). This high relative deviation is primarily because of small channel deformations as η → 1. The relative deviation is less than 15% for large channel deformations i.e. for 0.4 < η < 0.6. Further, we can see from the zeroth-order solution (equation (13)) that at large values of τ , the profile can be approximated by a sine curve with decreasing amplitude which agrees with the prediction from solving the full solution in figure 2. The above analysis indicates that the zeroth-order solution is  figure 2 (full solution) we would expect the velocity to decay from 1 to 0, reach a negative minimum, and then decay back to 0 for η < 0.5 and to decay gradually from 1 to 0 over time for η > 0.5. This behavior is also predicted by the zeroth-order solution as shown in the inset of figure 4(b). Further, velocity predictions by the zeroth-order solution match experimental observations as shown in the inset of figure 4(a). In the inset to figure 4(a), the velocity monotonically reduces from positive to negative values before reaching a negative minimum and gradually attaining a value of 0. This time for the velocity to crossover from positive to negative values at any η < 0.5 is the time required for the stagnation point to reach the corresponding η. The predictions by the zeroth-order solution for the time required for the stagnation point to reach a certain η match experimental observations as shown in figure 4(a). It can also be seen from figure 4(a) that the stagnation point approaches η = 0.5 asymptotically. We could not obtain times for the stagnation point to reach η < 0.3 due to the limitations in the frame rate of the camera. Further, the predictions by the zeroth-order solution for response times at all η are in good agreement with experimentally observed values as shown in figure 4(b). The V-shaped curve, almost symmetric about η = 0.5, is expected since the initial velocity is 1 for all η and the profile is almost symmetric about η = 0.5 for large values of τ (figure 2). We observe that the response time converges to a value ∼ 1.6 at the exit  (5) and (21). ε = 0 has been used for the zeroth-order solution. (b) Comparison of the deformation profiles for a PDMS channel after the application of a pressure pulse predicted by the analytical zeroth order and numerical full solution for (ε = 0.1 and 0.25) obtained using equation (5) for an early time (τ = 0.1).
of the channel (η = 1). This value of λ is also predicted from equation (23) by setting c = 1 (velocity reduces to 1% of initial velocity).
The measured value of the elastic modulus, 2 MPa, which is within previously measured limits (0.5 to 4 MPa) [2,3], was used to obtain the dimensional zeroth-order solution. The zeroth-order solution predicts a linear dependence of κ with (4 µWL 2 )/(H 3 E). The difference between this and our prior work is that the coefficient is explicitly calculated in the current work and hence can provide a priori predictions. The response times obtained experimentally are a good match to the response times predicted by the zeroth-order solution (solid black line in figure 5).

Conclusions
In conclusion, we have presented a comprehensive analysis of the temporally complex squeeze flow which occurs in a relaxing PDMS channel. Using a regular perturbation approach, we found an analytical solution to the problem which agrees with experiments for small channel deformations. The simple form of the resulting governing zeroth-order equation (a diffusion equation) both lends itself to closed-form analytic solutions and simple physical arguments to predict system trends such as the migration of the stagnation point toward the channel center or the dependence of the relaxation (response) time on channel dimensions and material properties. This analysis lends insight into many microfluidic processes which exploit channel deformation or, conversely, try to mitigate it.