Modelling and simulation of micro-well formation

Abstract

Physico-chemical processes on the micro-scale require new modelling concepts because some effects become dominating that are negligible for macroscopic systems. This is illustrated by a new method for the production of micro-wells based on the placement of a small drop of toluene on a plate of polystyrene. After droplet evaporation, a micro-well is left. A mathematical model has been developed to understand the elementary processes of the micro-well formation. The model accounts for: (1) growth of the drop on the substrate, (2) evaporation process of the solvent, (3) dissolution of the substrate, (4) flow rate in the evaporating drop caused by the pinning effect, including the vertical velocity profile, and (5) increase in the concentration of dissolved material followed by precipitation. In the modelling and simulation process, it could be shown that the method of drop production also has a significant influence on the shape of the micro-wells.

Keywords:

• Micro-wells
• Polymer
• Droplet
• Micro-chemistry
• Evaporation
• Dissolution
• Precipitation
• Pinning effect

1 Introduction

In many fields of modern engineering, miniaturization has become increasingly important. The same holds for different fields in chemistry. Methods for micro-structuring surfaces are required, such as ‘High-Throughput Screening’ in the pharmaceutical industry. Even though there are already several methods for surface engineering, not all demands can be satisfied by any one of these procedures. Thus, new methods in this field are still of high interest for research and industrial production. With the method presented here, it is possible for instance to provide the pharmaceutical industry with ‘test boards' with several hundreds of test facilities on a few square centimetres. In particular, expensive experiments can become much cheaper by miniaturization. Owing to improved understanding of the processes in production of the micro-structured surfaces, the processes in the cavities themselves (e.g. crystallization) can also be better understood.

As shown in figure 1, it is possible to form micro-wells by placing solvent drops on to a polymer substrate. In contrast to the procedures described in the literature [1 –  4], here the bottom is dissolved by the liquid. The solute remains after the drying process as a ring stain, and a micro-well is left. The process of evaporation is the subject of many publications [3, 4], as is the precipitation of dissolved materials [1, 2, 5]. An understanding of the influence of physical and chemical parameters such as concentration of the solvent in the vicinity, composition of the liquid phase, etc. is of high interest for the optimization of the process.

Figure 1. Atomic force microscope (AFM) image of a micro-well produced with the ‘sessile drop method’ (note the different scale on the z-axis).

The micro-wells can be used as small cavities for chemical experiments or for example as masks for the production of micro-lenses. Thus, it is also possible to think about forming micro-landscapes on demand for chemical processes and micro-reactors.

Though the process of diffusion of the solute molecules in the solvent is very well known, there is surprisingly no model for the solving step itself, i.e. the change from the solid phase into the liquid phase (see Section 3.2). Here, a model analogue to the heat transfer between two materials and the subsequent heat transfer in the materials was used. The constant describing this process was subjected to a fitting procedure.

2 Experimental results and methods

Several methods have been studied to investigate the major influences on the shape and the depth of the micro-wells produced, one of which will be presented here, namely, the ‘sessile drop method’. This will be explained in detail in the following, as well as the simulations concerned with this method.

With a syringe, a small amount of solvent (toluene) is placed onto the polymer substrate (polystyrene). The needle is almost in contact with the substrate. The toluene is brought to the surface with some pressure. The photographs in figure 2 were taken from a video sequence of about 2 s. The diameter of the micro-well was about 300 μm. A smaller micro-well (diameter = 88 μm) was imaged with an atomic force microscope (see figure 1).

Figure 2. Video sequence of a toluene drop evaporation on polystyrene (∼2.5 s, drop diameter ∼300 μm).

3 Theoretical background

The mathematical model developed to understand the process of micro-well formation recognizes the following physical and chemical effects of the system (figure 3):

1.	creation of the drop;
2.	evaporation process of the solvent;
3.	dissolution process of the substrate;
4.	flow in the evaporating drop caused by the pinning effect, including the vertical velocity profile;
5.	increase in the concentration of dissolved material followed by precipitation.

With the exception of the drop creation, the theoretical background for the processes will be described below. The technical details of the drop creation will be described in Section 5.1.

Figure 3. Physical and chemical processes in the evaporating drop.

3.1 Evaporation process

The process of evaporation starts when molecules of the solvent pass the liquid – gas interface. After this, the vapour molecules disappear in a diffusion process. Directly above the drop surface, the maximum possible concentration in air is reached. From here, the concentration decreases with increasing distance from the drop until the concentration of the environment is reached (normally 0 for toluene) [6]. If the concentration in the solvent in air is denoted by c, the diffusion constant with D and the vapour flow rate with J, the following equations are obtained:

The shape of the drop can be considered as a sphere since the diameters of the drops are very small ( < 1 mm). For such small radii, the spherical assumption is justified since the influence of the gravity on the shape of the surface can be neglected [3]. The maximum contact angle ϑ (≈4°) in our case was determined by analysing the measured data (figure 4).

Figure 4. Nomenclature of the variables.

Hu and Larson [3] have developed the following formula for the calculation of mass loss rate

due to evaporation through the surface S (normal to the surface) of the drop.

The factor

is the evaporation rate at the centre.

reflects the singular behaviour at the border of the drop.

accounts for the slope of the surface at a distance

from the centre. These two terms are transferred into one expression by an approximation with the phenomenological constant Λ.

To determine α₁ and α₂ in the expression for Λ, numerical values from FEM calculations were used in a previous study [3], which found values of α₁ = 0.2239 and α₂ = 0.3619, respectively. By a fitting procedure, the angular dependency of J ₀ (ϑ) can be expressed by J ₀ (π/2).

where

for the evaporation over the full surface, this produces:

for a fixed angle, finally an expression for the integrated mass loss rate is found:

Comparison of time for evaporation of the full drop with measured data gave a very good agreement for pure water. In the experiment with toluene and polystyrene, large discrepancies of the order of one magnitude are found. This can be explained by the increased number of solute molecules at the surface of the drop, since here the evaporation takes place. Thus, the ‘effective’ surface is reduced, and the total evaporation time increases. This phenomenon is described by Raoult's law [7] that scales the mass loss rate for the evaporating solvent with the factor:

Here, the concentration of the solute at the surface is denoted by C _surf, and C _max is the maximum possible solute concentration (see also Section 3.2). The concentration C _surf was extrapolated from the two highest cells of the discretization grid (see Section 5), so the fraction can be greater than 1. The factor γ was determined by a fitting procedure to be γ = 0.75.

3.2 Dissolution process of the substrate

Molecules of the substrate change into the liquid phase. Quantitatively, this is described by a desorption constant D _esorp. Up to the factor AΔt (area and time interval), the amount of material that is dissolved is proportional to the difference between the maximum possible concentration C _max of the solute in the solvent and the actual concentration C

The molecules then drift into regions with a lower concentration (i.e. the liquid – gas border) according to Fick's diffusion law. The change in concentration over time is thus described by the diffusion equation:

3.3 Flow in the evaporating drop

In the process of evaporation, the volume of the drop is decreased. At the same time, the three-phase border line remains fixed (pinning effect). Thus, the diameter of the drop stays (nearly) constant, and liquid has to be transported to the border. This causes the flow inside the drop. If the black line in figure 5 below represents the shape of the drop at a time t, the dash-dotted line is its shape—if it were not pinned—after some of the solute has evaporated. In contrast, if the drop is pinned, the grey line will emerge.

Figure 5. Pinning effect transporting liquid with the solute into the outer region of the drop.

Owing to the assumption that the flow rate near the liquid – substrate interface of the drop is almost 0, the profile of the flow rate is basically a function of the position in the drop. Hu and Larson [3] have given the following result for the profile of the radial component.

Here, the vertical position is denoted by z, h is again the height of the drop at a given radius r (figure 4), and v is the height-averaged flow velocity.

3.4 Increase in concentration of dissolved material followed by precipitation

Beginning at the border of the drop, the missing liquid has to be transported from the inner part of the drop. This flow mixes liquid with different concentrations of the solute. Owing to the pinning effect, solvent and dissolved material are thus transported to the outer region of the drop. The amount of dissolved material in the drop and the concentration of the solute, respectively, are increased due to

1.	the dissolution process;
2.	the fact that the amount of liquid is decreased by the evaporation process and;
3.	especially in the outer part of the drop, the transport of dissolved material.

if the concentration of solute in the solvent is above the maximum possible concentration, the dissolved material starts to precipitate. This process cannot be described in the same way as the process of dissolution. Since there is no detailed model for the precipitation available, the amount of dissolved material is simply reduced such that the averaged concentration at that distance is equal to the maximum allowed concentration C _max. As soon as the concentration of the solvent exceeds its maximum value, it starts to precipitate. In the implemented program (see Section 5) in one time step, the amount of material m _sed is calculated, which is necessary to decrease the concentration to the value. If the height of the liquid is h (see figure 6), the following expression is obtained

from which follows

Figure 6. Calculation of the amount of material that precipitates.

4 Data processing

Since the radial symmetry of the problem is exploited in the simulations, the measured data are averaged as described below to be able to compare the measured and the simulated data. For the preparation of the chart at the right of figure 7, the centre of the micro-well has been determined by a fitting procedure. After this, the height of the different pixels at a given distance from that point has been averaged.

Figure 7. Measured (left) and averaged (right) profile of a micro-well (‘sessile drop method’).

To determine whether the pinning of the drops has taken place directly or after some solvent has already evaporated, and whether dissolved material has its origin outside the micro-well, the mass conservation was checked. This procedure also checks if the density of the precipitated material agrees with that of the substrate material. To do this, the following integral is calculated as a function of its upper border

The zeros of this function are the radii at which mass conservation is obtained. With these values, the function

was also calculated (shown on the right-hand side in figure 8). In the example, the radius of the micro-well was determined to be 44 μm.

Figure 8. Integral function (left) and function (6) for the zero at 44 μm (right).

5 Simulation steps and parameter

The implementation of the simulation algorithm was done in MATLAB. The description of the simulation steps follows the order in Section 3. Since no dependency in the azimuth was expected, the calculations were performed in a cylindrical coordinate system. To be able to simulate the diffusion process of the dissolved material and to take into account the profile of the flow rate, the drop was discretized not only in the r-direction but also in the z-direction (figure 9).

Figure 9. Coordinate system used for the simulation.

The height of the drop at a given radius r from the centre is due to the spherical assumption (see Section 3):

for a given volume V, radius and contact angle, the height in the centre can be calculated by

Since, in all the processes described below, the height of at least one cell is changed, the cells have to be rescaled after each step, and thus the contents of different cells are mixed.

5.1 Growth of the drop

Since the flow rate of the solvent from the syringe is not known precisely at present, this was a parameter in the program which was successively adapted. For the initialization, the drop was assumed to have a volume that was below 0.001 of the maximum volume, after which the process of drop growing was stopped. To take into account the fact that the solvent is injected with some pressure, the added solvent is put into the lowest cell of the computational grid in the centre of the drop. Thus, it is possible to use the flow model that is used for the description of the pinning effect for the spreading of the liquid. Since the drop now is considered not to be pinned, the radius of the drop is increased with increasing volume such that the contact angle ϑ is kept constant. The other processes (evaporation, etc.) are, of course, simulated, too. The flow rate was kept constant while the drop was filled.

5.2 Evaporation process of the solvent

From the formula (1) for the mass loss, the evaporation volume can be calculated, since the density of the liquid is known. In the program, the mass loss for the full surface was compared with equation (2). Discrepancies owing to numerical discretization errors could be avoided by rescaling the evaporated amount of liquid after each step. This is especially important for small radii, because the number of intervals in radial direction is low. The contact angle changes after the growing drop process has ended due to evaporation and the pinning effect. The angle is calculated from the drop radius and the height in the centre of the drop

While filling the drop, the length of a time step is calculated such that the height of the cell in which the solvent is filled is not increased by a factor greater than 2. When the filling of the drop has ended, the length of a time step is determined such that in the interval with the highest evaporation, half of the solvent is left. Of course, there is a lower limit for the length of the time step so that the solvent in the cell can be completely removed.

5.3 Dissolution process of the substrate

The amount of material from the substrate that is dissolved is determined by equation (3). In the liquid phase, the solute drifts to regions with a lower concentration by the diffusion process. In the program, this diffusion process was implemented with the implicit Euler method for equation (4).

5.4 Flow in the evaporating drop

The height of the drop at a distance r from the centre is changed by the injection of liquid, the evaporation, the dissolution of the substrate and the deposition of material that precipitates. It should be taken into account that some volume of the liquid is in the hole formed by the dissolution process, and the deposited material also forms a socket (see also figure 3). The shape of the drop is determined for the remaining volume as a sphere. On the other hand, the remaining height in the calculation (that reflects all described processes) at a given radius is different (see figure 10 left). Thus, the amount of liquid that has to flow to this interval is determined (see also figure 5). The resulting flow is shown on the right-hand side in figure 10.

Figure 10. Height difference (left) and resulting flow rate (right).

5.5 Increase in concentration of dissolved material followed by precipitation

The concentration of the solute in the outer region of the drop is drastically increased (figure 11).

Figure 11. Increased concentration in the outer part of the drop (with time as parameter).

The amount of material that is deposited in the outer part of the drop is higher than the liquid after some of the volume has evaporated. In the program, the diameter of the drop is reduced at that moment. For the calculation of the sphere, the volume in the remaining socket and also in the ‘hole’ in the substrate has to be taken into account.

5.6 Parameters used in the program

The simulation was terminated as soon as the remaining volume of the drop was below 10 ^– 4 of the maximum volume of the drop (see table 1).

Table 1. Physical parameters.

Meaning	Value
Parameters of the drop:
Radius of the drop
Maximum contact angle
Diffusion constant of polystyrene in toluene
Desorption constant
(Concentrations are calculated in volume per cent)	R = 44 μm
ϑ = 0.0675 rad
D solv = 2.7 × 10^– 7 cm^2 s^– 1
D esorp = 0.013 cm s^– 1
Parameter of the solvent:
Density of toluene	ρ = 0.86 g cm^– 3 (from [8])
Parameters of the gas phase:
Saturation concentration for toluene
Diffusion constant	C sat = 0.00016 g cm^– 3
D = 0.0806 cm^2 s^– 1, both values from [8]
Parameter of the substrate:
Density of polystyrene	ρ = 1 g cm^– 3

6 Comparison of the experimental and simulated results

To compare the simulated results with the experimental results, the parameters for flow rate of the solvent V _flow, the desorption constant D _esorp (figure 13), and the maximum possible concentration C _max (figure 12) are changed.

Figure 12. Measured data (left) and simulated (right) data (parameter C _max).

Figure 13. Variation of the parameters D _esorp (left) and V _flow (right).

The results from the simulations are in reasonable agreement with the experiments. With a more sophisticated and detailed model of the flow in the drop, especially during the formation process of the drop, the accuracy can certainly be increased, and the results will even fit the experimental data better. The mismatch between simulated and measured evaporation time indicates that the concentration change at the drop surface should be considered in a refined model.

7 Conclusion

It turns out from the simulation that the drop growth is crucial for the shape of the micro-well. This is because in the centre of the drop near the liquid – solid interface, the concentration of the dissolved material in the solvent is always low. This results in a high dissolution rate of the polymer into the solvent drop and thus increases the amount of polymer deposited at the drop rim. Thus, it is possible to ‘drill’ deep holes into the substrate.

By combining simulation and experiment, it is possible to determine the dissolution constant of a polymer into a solvent. This constant is experimentally difficult to measure because the dissolution in general is accompanied and determined by a diffusion of the polymer. Thus, the dissolution constant and the diffusion constant are strongly coupled in generally.

8 Outlook

The work will be continued as a part of research group founded by the Deutsche Forschungsgemeinschaft (DFG) FOR 516/1 in the projekt WI 1705/7-1. Subject of the further analysis will be a detailed parameter study to determine the experimental sensitivity of the whole system with respect to the chosen parameters values and related processes. These theoretical efforts will be accompanied by experimental investigations. Hence, the intention is to prove the valid range of the models presented before and especially to focus on arising discrepancies. From the viewpoint of simulation, it would seem worthwhile focusing on implementing more sophisticated numerical methods.

Acknowledgements

This work was partially funded by the Deutsche Forschungsgemeinschaft within the DFG research group 516 ‘Physical and chemical foundations, components and systems for lab-on-chip technology’. We gratefully thank Guangfen Li and Elmar Bonaccurso for providing us with the experimental data.