Modeling and characterization of surface roughness effect on fluid flow in a polydimethylsiloxane microchannel using a fractal based lattice Boltzmann method

A finite composite of patches and roughness exist on almost all surfaces. Mostly they are unnoticed or ignored. With the evolution of miniaturization, the surface to volume ratio (S/V) in almost all devices is increasing rapidly and surface properties show effect on physical parameters. In this work, this is being demonstrated for microfluidics application. The surface related properties have now emerged as a crucial phenomenon in determining characteristics and performance of the device especially in nano and some micro devices. Thus, it has become necessary to model the intrinsic surface property such as surface roughness to understand the microscopic phenomenon and to develop a model based on the surface roughness information.^1,2 The surface roughness can affect the physical, chemical, hydrodynamic and electrical properties of any device.³ Recently, microfluidics devices are increasingly getting popular due to their versatile applications in micro and nano domain. Materials used in these devices are crucial. In the recent past, the microfluidics researchers have moved from the traditional MEMS fabrication materials like silicon, quartz to the new materials such as polydimethylsiloxane(PDMS), polymethylmethacralate(PMMA) and other polymers in order to minimize the cost of the microfluidic device and various reasons including biocompatibility.⁴ Though the surface fabrication and interfacing skills of such devices are well developed, the investigation on fundamental phenomenon such as surface roughness has been least addressed in the past. Usually, in any microfluidics device, the size of the microchannel is very small. So, the fluid flow rate is generally expected to be low. Moreover, the turbulence inside a microchannel is usually unnoticed. Thus, in general the flow is assumed to be fully laminar⁵ due to the reduction in dimensionless Reynolds number up to the order of 10⁻³. However, the oscillations can be present atleast locally on the surface due to the intrinsic variations like roughness. At sub-micron and nano scale, these variations cannot be ignored especially for multi scale varying microfluidics problem. To model such a multiphysics, multicomponent critical microfluidic problem, the traditional Navier-Strokes(NS) approach is not generally convenient due to the possible failure in satisfying continuum hypothesis throughout the structure. Moreover, obtaining convergence to the solution in a reasonable computational time is a challenging task. Recently a mesoscopic Lattice Boltzmann method (LBM) has been effectively proven as an alternative fluid flow simulations due to its tremendous advantages like massive parallelism and easy incorporation of multiphysics.^6,7

In this work, the effect of surface roughness on a fluid flow in a PDMS microfluidic device using the fractal based Weierstrass-Mandelbrot (WM) model is studied. At first, we have fabricated a PDMS microchannel using relatively a simple procedure. Then, Atomic force microscroscope (AFM) images of the microchannel were captured. Using the roughness detail, the statistical fractal dimension (D) is measured. Finally, the corresponding rough surface is modeled using WM function for LBM simulations. Generally, to fabricate a microfluidic device, a soft lithography process is used.⁸ In such soft lithography process, a photoresist (SU-8) mold is usually created using UV lithography. But in contrary to that, we have not used the expensive UV-lithography set up involved, by using a simple copper(Cu) printed circuit board (PCB) as a mold as shown in Fig. 1. Though this method has limitations in terms of the achievement of minimum feature size, it still can be used to fabricate the microfluidic devices with dimensions in the range of few tens or hundreds of microns. In this study, the minimum width of the fabricated microchannel is 300μm and depth is 70μm. The devices with these features can be used for many low-cost Bio-MEMS applications. For fabrication, an elastomer precursor (sylgard 184) and its curing agent are mixed in 10:1 ratio and degassed using a desiccator. The solution is then poured over PCB-mold on a flat surface in a glass petri dish and cured at 80°c for 40 min. After hardening, the replica is peeled off from the mold, and the reservoir holes were punched. The PDMS channels are then bonded to a glass slide using an adhesive bonding technique.⁹ For bonding, a fresh PDMS solution is spin coated on a glass slide at 8000rpm for 8min. As a result, an uncured PDMS thickness of 1 - 1.5μm is obtained on a glass substrate. Above the fresh PDMS coated glass substrate, the fabricated PDMS micro channel is gently placed without allowing air bubbles in the assembly. Finally the assembled device is kept flat at 90°c for 15min to cure the adhesive bonding layer. All the processes are done in a class 10000 semi clean room environment. In general, the roughness of a surface is represented by its average roughness or root mean square (RMS) roughness value.¹⁰ Additionally, it can also be represented by its frequency properties such as power spectrum. RMS roughness information provides the vertical magnitude of roughness, but it does not give spatial information. Various approaches have been discussed in the past to model the surface roughness such as aspect ratio, the height of protuberances, distribution, shape, and density of surface features.¹⁰ In all these measurements, the parameters such as power spectrum, RMS values, and autocovariance are used to simulate the surfaces. However, these methods are not adequate for the surface roughness modeling of the typical microelectronic films and polymers. Moreover, all these parameters vary greatly across different instruments, scan size and scan location.¹¹ On the other hand, fractal based models are proving to be a better approach to the roughness measurement.¹² In this method, the roughness is represented as a function of fractal dimension(D), which is independent of scan size and also describes the spatial information adequately. The fractal dimension is described by

where N is the number of parts scaled down and r is the ratio of a number of parts scale down to the total scale. In this case, the fractal dimension of the fabricated PDMS coated glass substrate of the microchannel is calculated using the simple flooding method of the WSxM SPM software tool¹³ as shown in Fig. 2, where S is the scan area and P is the perimeter of the filled holes and hills of the AFM image. From the slope of the plotted curve, it is found that the PDMS coated bottom surface posses fractal dimension of D=1.4. In order to model the roughness of the fabricated PDMS microchannel in LBM, the popular Weierstrass-Mandelbrot(WM)¹⁴ function is used. The WM function is described by

In this equation, the parameter D is the fractal dimension, ϕ is randomly generated phase, b is the frequency multiplier value varies typically between 1.1 to 3.0. Usually, for the fairly regular surfaces,¹ k is kept as 1 and b is chosen to be 2. So, by adjusting D and k, any rough profile can be generated as shown in Fig. 3. In fact, from this generated topography one can also find the RMS values and other statistical parameters. To calculate the RMS value, for each coordinate (x, y) in the plane, the square root of(x² + y²) should be found and substituted for x in Eq. (2). So, given the fractal dimension, one can realize and tune any roughness profile using the WM function. Thus, to model the surface roughness of the PDMS microchannel in the LBM, the roughness profile corresponds to D=1.4 is generated using WM function. For the experimental characterization of the fabricated device, a microfluidics characterization setup is used.¹⁵ The complete setup consists of a Universal serial Bus (USB) digital camera for capturing images/video in real time, XYZ stage for holding the sample, two syringe pumps, Arm7 processor based electronic system to handle image processing, along with High definition multimedia interface (HDMI) display. The two different fluids used are paper ink of blue and red color. The fluids enter the mixing chamber through inlet1 and inlet2 of the Y-mixer simultaneously with the controlled flow. The flow rates of both the fluids are fixed to be 0.1ml/hour using a syringe pump. The fluid studied has the kinematic viscosity(ν) similar to water as 0.9e⁻6 m²/s at room temperature. In this case, the dimensionless Reynolds number(Re) is obtained as 0.02231. The LBM simulation parameters are chosen for this Reynolds number. To calculate the velocity of the flow experimentally in the microchannel, the chain of images is captured from the top of the device using a digital camera at 30frames/second. The velocity of the fluid inside the channel was measured with this set-up. To capture the hydrodynamics by simulation, we have used a popular three dimensional nineteen velocities (D3Q19) LBM model. In this approach, the evolution of the particle distribution function (f_i) at the point x at time t is given by the lattice Boltzmann(LB) equation in the predefined 19 directions(e_i) of the model as

Here, the left-hand side of the above equation is responsible for streaming and the right-hand side for collision with i=0, 1,..18 and τ is the dimensionless single relaxation time in the Bhatnagar-Groos-Krook (BGK)¹⁶ collision operator approximation. Here $fieq$ is the equilibrium distribution function which will be defined initially as

where c is defined as$13.$ and $w$_i are pre defined standard lattice weights. The macroscopic density and velocity can be computed by adding all f_i using Eqs. (5) and (6) respectively.

In simulation, the optimum lattice domain of 300 X 300 X 70 is computed for 2000 lattice time steps. We have done the relative lattice grid refinement study on three cases on 270 x 270 x 63, 300 x 300 x 70, 330 x 330 x 77 and relative error is found to be in the order of e-8. Hence, 300 x 300 x 70 is optimum with δ_x = 3.33e⁻³, δ_t = 0.5e⁻³. The top, bottom, front and back walls are defined as no slip boundary walls. The left and right walls are defined to be periodic. The dimensionless relaxation time(τ = ν/c²δ_t + (1/2)) is set to be 0.6 and applied flow rate(0.1ml/hr) is set as 2.7778e⁻11.

In order to study the effect of fractal dimension on the velocity, D is varied between 1 to 2. The velocity is measured for all realized surfaces across the height of microchannel from the bottom most lattice point to the top most lattice point at steady state and plotted as shown in Fig. 4. From the plot, it is clear that increase in D, results in an increase in the roughness and hence offered more shear resistance to fluid flow. This agrees with other work also.¹⁷ This leads to the corresponding decrease in the velocity of the fluid flow near surface. For the most rough surface (D=1.9), the roughness can even block the flow partially in the channel. In this case, this partial channel blockage is observed upto 28.5% of the total channel height. Note that in simulation, only bottom surface is modeled with measured roughness. If we include roughness on all other side, then this factor will be muliplied accordingly.

Fig. 5 shows the comparison of LBM simulated result with the experimental result. Experimentally, velocity is measured using a kymograph image processing technique using characterization setup described above. The values are taken after averaging it to the five iterations. Simulation results are plotted across the width of microchannel at the middle from the left corner to the right corner at time 2ms. The simulation is performed without adding roughness and then with the roughness. For the simulation with roughness, the experimentally measured D, (D=1.4) is substituted. From the plot, it is clear that the simulation result, without considering roughness is flat and ideal, indicating velocity is uniform across the channel. However, experimental results show that the velocity is not constant and has oscillatory behaviour. The simulation without roughness failed to capture the oscillatory behaviour in the channel. These oscillations are occurred due to the presence of roughness in the surface of the channel. On the other hand, the inclusion of a WM surface roughness model significantly improved the accuracy by capturing the oscillatory behaviour in the middle of the channel. Thus, this approach is proved to be effective in capturing the real phenomenon in the numerical simulations. The simulation result also provided the insight on the vorticity behaviour near the surface. From the Fig. 6, it can be observed that the presence of roughness on the surface created the local oscillations. This oscillatory behaviour can well be utilized for the microfluidic activities like mixing. In fact, one can even tune the mixing by tuning the required roughness by fabrication process optimization.

In conclusion, a fluid flow in the microchannel for a PDMS based microfluidic application has been studied. The experimental results show that the velocity is not constant near the surfaces but shows oscillatory behavior. This behavior can be explained by accounting details of roughness on the surfaces of the channel. The local oscillations in the microchannel due to roughness in the wall has been captured. The Weierstrass-Mandelbrot function has been used to model fractal based surface roughness. The fractal dimension of the fabricated device is measured experimentally using AFM and substituted in WM function of the model. The fluid flow was simulated using 3D LBM method. Obtained results show the oscillatory behaviour of the fluid flow. The result also shows an increase in the fractal dimension increases the roughness and hence reduces the fluid velocity. The result suggests that by synthesizing rough surfaces, one can tune the turbulence in the microchannel and improve the mixing efficiency.