5.1.2 Modeling Droplet Transport

Forces acting on the droplet can be used in order to calculate the location where the droplet should reach the surface. This is a challenge because the simulation environment must now be divided into several segments. For an in-depth explanation of the forces involved during the droplet transport, refer to Section 4.1.3. In this section, the modeling of those forces and the droplet movements will be described. The first segment we must treat separately is the thermal zone which is within 10mm of the wafer surface. In this area, the temperature gradient shown in Figure 4.5 is high and the thermal forces play a significant role in the droplet speed as well as size, due to evaporation. In addition, when the electrical force is included, the complexity of the problem is significantly increased.

This section will first examine the droplet's motion with no influence from the electrical force outside of the thermal zone. Subsequently, the electrical force is introduced and its influence on the droplet motion is described. Finally, the thermal zone is covered by the introduction of the thermophoretic force and the idea of droplet evaporation. Figure 5.1 shows how the simulation space is divided in order to accommodate the thermal zone for droplets.

Figure 5.1: The droplet transport in the space above the substrate surface and the accelerations which are considered in the transport model. T$_{th}$ is the height of the thermal zone ($\sim$10mm for ESD, $\sim$5mm for PSD), and H is the distance between the substrate and atomizer.

For a more detailed explanation regarding the derivation of the motion equations used in this section, refer to app:transport.

5.1.2.1 Droplet transport using gravity and Stokes forces

In order to follow the trajectory of a droplet after leaving the atomizer and under the influence of gravity and the Stokes force, the distance required to reach the thermal zone $\left(H-t_{th}\right)$, the initial velocity $v_0$ are the droplet radius $r_d$ are known. The equations for the gravitational force and the Stokes force can be reformulated to equations for acceleration $a\left(t\right)$ which govern droplet movement. It is important to note that the vertical motion must be calculated first in order to know how much time $t$ is required to reach the thermal zone. The time $t$ is then used in order to calculate how much radial distance the droplet underwent during its vertical ``fall''. For the sake of simplicity, the vertical acceleration, velocity, and distance will be marked with a subscript $v$, while radial acceleration, velocity, and distance will be marked with a subscript $r$ in the equations that follow.

Vertical trajectory:
Given the droplet mass (4.2), the force of gravity (4.3) and the Stokes force (4.7), the vertical acceleration experienced by the particle is given by:

$$a_{v}\left(t\right)=g-\cfrac{9\,\eta_{a}}{2\,\rho_{d}\,r_d^{2}}\,v_{v}\left(t\right),$$ (162)

where the drag component of the acceleration depends only on velocity and is direction-independent. It is defined as

$$s=\cfrac{9\,\eta_{a}}{2\,\rho_{d}\,r_d^{2}}.$$ (163)

(5.3) resembles a first order differential equation which can be solved for the vertical velocity $v_{v}\left(t\right)$ using characteristic equations

$$v_{v}\left(t\right)=\cfrac{g}{s}\left(1-e^{-s\,t}\right)+v_{v0}e^{-s\,t},$$ (164)

where $v_{v0}$ represents the initial vertical velocity. The velocity can now be integrated in order to find an explicit equation for the distance which is the known value $\left(H-t_{th}\right)$.

$$d_{v}\left(t\right)=\cfrac{1}{s}\left(\cfrac{g}{s}-v_{v0}\right)\left(1-e^{-s\,t}\right)+\cfrac{g}{s}t,$$ (165)

where it is assumed that $d_{v0}=0$ and $d_v\left(t\right)=\left(H-t_{th}\right)$, although it is also possible to set up the problem such that $d_{v0}=H$ and $d_v\left(t_{final}\right)=t_{th}$. In order to know how far the droplet moved in the radial direction, (5.6) must be inverted, giving the time $t$

$$t=\cfrac{gW\left[\cfrac{\left(g-s\,v_{v0}\right)\cdot e^{1-\cfrac... ...\right)\right)}{g}}}{g}\right]- g+s\,v_{v0}+s^{2}\left(H-t_{th}\right)}{g\,s},$$ (166)

where $W\left[\cdot \right]$ is the Lambert-W function defined by

$$z=W\left(z\right)e^{W\left(z\right)}.$$ (167)

The Lambert-$W$ function is solved iteratively, as described in [215]. Now that the time required to reach the heat zone is known, the radial trajectory can be calculated.

Radial trajectory:
In order to calculate the radial trajectory, the time $t$, radius $r_d$, and initial velocity $v_{r0}$ must be known. Those parameters are derived using the vertical trajectory discussion. The force of gravity does not influence the radial movement of the droplet; therefore, only the Stokes force must be considered. The acceleration of the droplet is defined as

$$a_{r}\left(t\right)=-\cfrac{9\,\eta_{a}}{2\,\rho_{d}\,r_d^{2}}\,v... ...ft(t\right)\qquad\textrm{or}\qquad a_{r}\left(t\right)=-s\,v_{r}\left(t\right).$$ (168)

Similar methods used for the vertical trajectory can be used to solve the radial velocity

$$v_{r}\left(t\right)=v_{r0}e^{-s\,t},$$ (169)

and the radial distance

$$d_{r}\left(t\right)=\cfrac{v_{r0}\left(1-e^{-s\,t}\right)}{s}.$$ (170)

The vertical and radial location as well as the velocity of the droplet is now known as it contacts the thermal zone. The radius remains unchanged, as it was experimentally shown that the droplet radius is largely unaffected outside of the heat zone [171]. The droplet carries information regarding its position, velocity, and radius as it passes to the next step in the model.

5.1.2.2 Droplet transport inside the heat zone

As the droplet enters the heat zone, it experiences a large temperature gradient resulting from a rapid increase in temperature. This causes the droplet to be exposed to an additional retardant force component, which pushes it away from the hot surface (4.9). This force is assumed to have a uniform behavior throughout the heat zone, making it straightforwardly implemented in the previously-discussed trajectory calculations. The thermophoretic force is independent of the droplet velocity, therefore it causes a reduction of the gravitational force effect. This is performed by replacing $g$ in (5.3) to (5.7) with a value $\left(g-l\right)$, where $l$ is the addition of the negative acceleration caused by the thermal gradient

$$\left(g-l\right)=g-\cfrac{27\,\eta_{a}^{2}\,\kappa_{a}\,\nabla T}{4\,\rho_{a}\,\rho_{d}\,T\left(2\kappa_{a}+\kappa_{d}\right)r_d^{2}},$$ (171)

where the temperature $T$ is defined in Kelvin and $\nabla T$ in Kelvin$\slash$m.

An additional effect which occurs in the heat zone is the significant increase in the mean droplet radius noticed in measurements [171]. The reason behind the increased mean radius is that droplets with a small radius evaporate before reaching the surface, while larger droplets, although also evaporating slightly, stay relatively complete until fully in contact with the surface [171]. Tracking of the detailed changing droplet size during its travel through the heat zone does not modify the model enough to merit the additional computational expense. Therefore, the model automatically excludes droplets which bounce away from the surface due to the thermophoretic force, while other droplets which reach the surface have their radius reduced as they enter the heat zone. The approximate relationship which governs the small droplet's lifetime is given by [81]

$$t_{life}=\cfrac{4\,r_{init}^2}{q_0\,\triangle T}\,,$$ (172)

where $r_{init}$ is the initial droplet radius, $q_0$ is an evaporation rate constant, and $\triangle T$ is the temperature difference within the droplet. After calculating the time required for a droplet of initial radius $r_{init}$ to traverse through the heat zone, if this time is larger than $t_{life}$, then the droplet is removed from the simulation, since it will not contribute to the film growth.

The exact solution for the decrease of the radius of a droplet requires the solution of a diffusion equation, since the evaporation of a droplet is given by [81]

$$\cfrac{dr_d}{dt}=\cfrac{-M_{W}D_{v,f}}{r_d\,\rho_{d}\,R\,T_{f}}\triangle p\left(1+0.276\,Re^{1/2}Sc^{1/3}\right),$$ (173)

where $M_W$ is the molecular weight of the evaporating liquid, $D_{v,f}$ is the average diffusion coefficient for vapor molecules in the saturated film around the droplet at the final temperature $T_f$, $R$ is the gas constant, $\triangle \rho$ is the pressure difference between the vapor pressure near the drop and the ambient pressure, and $Re$ and $Sc$ are the dimensionless Reynold's and Schmidt's numbers, respectively, given by

$$Re=\cfrac{2\rho_{a}v\left(t\right)r_d}{\eta_{a}},\quad\textrm{and}\quad Sc=\cfrac{\eta_{a}}{\rho_{a}D_{v,f}}.$$ (174)

The average diffusion coefficient $D_{v,f}$ is estimated using the Stokes-Einstein relationship

$$D_{v,f}=\cfrac{R\,T_f}{N_A}\,\cfrac{1}{6\pi\,\eta_{d}\,r_d}\:,$$ (175)

where $N_A$ is Avogadro's number. The droplet radius $r_d$ and velocity $v\left(t\right)$ are time dependent; therefore, the change in radius through the heat zone can only be solved numerically using a time discretization technique. However, large droplets do not experience significant size reduction through a zone with a non-zero temperature gradient $\nabla T$. Also, it can be assumed that a sedimentation velocity is reached relatively quickly and does not change with time, the diffusion of the droplet can be approximated linearly by

$$K=q_{0}\triangle T\left(1+2q_{1}r_d\right),$$ (176)

where $K$ represents the surface evaporation rate in $\left(m^2 \slash s \right)$ with $K=-\cfrac{dr_d^2}{dt}$, while $q_0$ and $q_1$ are given by

$$q_{0}=\cfrac{2a}{\triangle T}\left(1+b\,s_{0}\right),\qquad q_{1}=\cfrac{b\,r_{0}}{1+bs_{0}},$$ (177)

where $r_0$ and $s_0$ are constants given by $r_0=64.65$s$^{-0.5}$ and $s_0=-1.117\times 10^{-3}$ms$^{-0.5}$ and the variables $a$ and $b$ are given by

$$a=\cfrac{4\gamma M_{W}D_{v,f}}{\rho_{d}R}\cfrac{\triangle T}{T_{f}},\qquad b=0.276\left(\cfrac{\rho_{a}}{\eta_{a}D_{v,f}^{2}}\right)^{1\slash6},$$ (178)

where $\gamma$ is a constant [81]. With the goal of a topographical simulation in mind, a full detailed analysis for the droplet size, as it changes in the heat zone is not merited. Therefore, in order for (5.13)-(5.19) to be included in the model, some assumptions are made:

1. Droplets with a radius which is too small, giving a small $t_{life}$ will not be taken into account since those droplets will never reach the surface.
2. The droplet radius does not change during its transport through the heat zone for the calculations of the forces acting on the droplet. Rather, it is represented by a single drop in radius, which is calculated following the above discussion.
3. The velocity of the droplet through the heat zone is assumed to be constant in order for the above analysis to be valid and to calculate the size reduction due to thermal effects. The change in velocity from the time when the droplet enters the heat zone until it reaches the surface will, nevertheless, be calculated using the forces at play and the modified droplet size from (5.13)-(5.19).

5.1.2.3 Droplet transport including the electrical force

When an ESD system is used for spray pyrolysis deposition, an electric field is introduced between the atomizing needle and the metal plate on which the substrate is placed. This electric field provides additional acceleration for the droplets so that they do not need to rely only on the gravitational force, as is the case for PSD systems. The addition of the electrical force to the equations which govern droplet motion is not as straight-forward as was the case with the thermophoretic force, because the electric field which is generated in the experimental region is not uniform [228]. The external electric field is calculated by representing the atomizer as a semi-infinite line of charge and the substrate as an infinite conducting plane, separated by a distance $H$. The equation governing the external electric field is given by [96]

$$\vec{E}_{ext}=\cfrac{\Phi_{0}}{H}\nabla\Phi^{*},$$ (179)

where $\Phi_0$ is the applied electrical potential between the nozzle and the substrate and $\Phi^*$ is a normalized potential ( $\Phi\slash\Phi_0$), given by

$$\Phi^{*}\left(r^*,z^*\right)=\cfrac{K_{V}}{\log\left(4H\slash R\r... ...(1-z^*\right)}{\sqrt{r^{*2}+\left(1+z^*\right)^{2}}+\left(1+z^*\right)}\right],$$ (180)

where $R$ is the outer radius of the nozzle, $K_V$ is a value which ranges from 0 to 1 depending on the $H\slash R$ ratio [58], and $r^*$ and $z^*$ are normalized radial and vertical distances, respectively, given by

$$r^{*}=\cfrac{d_{r}\left(t\right)}{H}\quad\textrm{and}\quad z^{*}=\cfrac{d_{v}\left(t\right)}{H}.$$ (181)

For the purposes of spray deposition, it is often assumed that the value of $K_V$ is 1, because the ratio of $H\slash R$ is on the order of several hundreds, which gives a value close to 1 for $K_V$. This value is adjusted in the model using the relationship

$$K_V=1-e^{-0.021\frac{H}{R}}.$$ (182)

In fact, assuming that $K_V=1$ can cause erroneous results for the electric field. The negative exponential dependence on $R$ from (5.23) is in the numerator of (5.21), which shows an additional inverse logarithmic dependence in the denominator. A plot of the fraction $K_{V}\slash\left[\log\left(4H/R\right)\right]$ for various $R$ values is shown in Figure 5.2, when the variation in $K_V$ is taken into account and when it is assumed that $K_V=1$. It is clear that the effects of $K_V$ should be included in the droplet transport model.

Figure 5.2: The effects of varying the atomizing nozzle's outer radius on the strength of the electric field with and without the inclusion of the $K_V$ effects.

The individual vertical and radial components of the electric field when taking the gradient of $\Phi*$ from (5.21) can be found

$$\cfrac{\delta\Phi^*}{\delta z^*}=-\cfrac{K_{V}}{\log\left(4H\slas... ...t(1+z^{*}\right)^{2}}}+\cfrac{1}{\sqrt{r^{*2}+\left(1-z^{*}\right)^{2}}}\right]$$ (183)

and

$$\cfrac{\delta\Phi^*}{\delta r^{*}}=\cfrac{K_{V}}{r^*\log\left(4H\... ...*}\right)^{2}}}-\cfrac{1-z^{*}}{\sqrt{r^{*2}+\left(1-z^{*}\right)^{2}}}\right].$$ (184)

Figure 5.3 shows the value for the normalized potential $\Phi^*=\Phi\slash\Phi_0$ and its distribution in an ESD deposition setup. The atomizing nozzle is located at $\left(r^{*},z^{*}\right)=\left(0,1\right)$. The inset shows the electric field distribution in the same simulation space. It is evident that the strength of the electric field is not uniform or linear, but that the field causes charged droplets to spread radially.

Figure: Magnitude of the normalized electric potential $\Phi\slash\Phi_0$ during ESD processing. The distance between needle and deposition plate as well as the radial distance from the center are normalized to the distance between the atomizer and the substrate. The inset is the normalized electric field distribution.

Given the electric field distribution shown in Figure 5.3 (inset), which is dependent on the droplet position, the electrical force (4.4) cannot be implemented to find the vertical displacement using an explicit function in the form (5.6) from the analysis where the droplet is affected only by the gravitational and Stokes forces. The vertical acceleration, when the electric field is considered is no longer (5.3), but it changes to

$$a_{v}\left(t\right)=(g-l)-\cfrac{9\,\eta_{a}}{2\,\rho_{d}\,r_d^{2... ... \sqrt{\cfrac{\gamma\,\epsilon_{0}}{r_d^{3}}}E\left(d_{v}\left(t\right)\right).$$ (185)

In addition, the radial acceleration changes from (5.9) to

$$a_{r}\left(t\right)=-\cfrac{9\,\eta_{a}}{2\,\rho_{d\,}r_d^{2}}\,v... ... \sqrt{\cfrac{\gamma\,\epsilon_{0}}{r_d^{3}}}E\left(d_{r}\left(t\right)\right).$$ (186)

Simplifying (5.26)-(5.27) with the assumption of a linear dependence of the electrical force on displacement $d\left(t\right)$ and dividing by the droplet mass gives the electrical component of the displacement-dependent acceleration

$$c_v=\cfrac{6}{\rho_{d}}\sqrt{\cfrac{\gamma\,\epsilon_{0}}{r_{d}^{... ...cdot\cfrac{\Phi_{0}}{H}\cdot\cfrac{K_{V}}{\log\left(4H\slash R\right)}\,c_{ev},$$ (187)

where $c_{ev}$ is a vertical linearization constant which is meant to estimate the $\left[\cdot\right]$ term in (5.24) with

$$c_{ev}\cdot z^{*}\approx\left[\cfrac{1}{r^{*2}+\left(1+z^{*2}\right)}+\cfrac{1}{r^{*2}+\left(1-z^{*2}\right)}\right].$$ (188)

A similar estimate is made to find the radial linearization constant to estimate the electrical force effect on the droplet's acceleration for (5.25). Realizing that the location of the droplet depends on the time it travels, the $z^*$ from (5.29) can be viewed as the time-dependent displacement. With this linearization assumption, a system for the vertical acceleration from (5.26) can be re-stated as a combination of the displacement-dependent electrical acceleration term $c_v$, velocity-dependent Stokes acceleration term $s$, and the independent gravitational $g$ and thermal $l$ acceleration terms

$$a_{v}\left(t\right)=(g-l)-s\, v_{v}\left(t\right)+c_{v}\, d_{v}\left(t\right).$$ (189)

Similarly, for the radial acceleration, (5.27) is re-stated as

$$a_{r}\left(t\right)=-s\, v_{r}\left(t\right)+c_{r}\, d_{r}\left(t\right).$$ (190)

The derivation of a displacement equation, stemming from the acceleration (5.30) is discussed in Section A.2, and the resulting displacement is given by

$$d_v\left(t\right)=C_{1}e^{-r_{1}\,t}+C_{2}e^{-r_{2}\,t}+C_{3},$$ (191)

where

$$r_{1}=\cfrac{-s+\sqrt{s^{2}-4\,c_{v}}}{2},\quad r_{2}=\cfrac{-s-\sqrt{s^{2}-4\,c_{v}}}{2},$$ (192)

and

$$\begin{array}{l} C_{1}=r_{1}v_{v0}-\cfrac{r_{1}^{2}\, r_{2}}... ...\\ \\ C_{3}=\cfrac{c_{v}}{1+\left(g-l\right)+s} \end{array},$$ (193)

where a similar expression can be found for the radial displacement using $v_{r0}$ and $c_r$ as parameters in (5.33) and (5.34). The electric field is strongest at the nozzle tip, as shown in Figure 5.3. Therefore, the field will be calculated using the given equations once the droplets are generated in order to find the droplet's initial direction and to give it its vertical and radial acceleration. However, a linear approximation for the vertical and radial electrical force will be used afterwards.

An explicit equation for the time required to achieve the displacement (5.32) cannot be found. Therefore, the droplet motion must be solved by time discretization, MC methods, or iteratively in order to obtain the droplet trajectory through the electric field. For the presented topography simulator, an iterative method is used to solve for $t$ with the initial guess $t_{init}$ given by (5.7). Since the introduction of the electric field causes the required time for a droplet to reach the surface $t$ to decrease, the time $t$ with the electrical force included must be $0<t<t_{init}$. For the purposes of a topography simulation, an assumption for a constant electric field is sufficient [171]. It is also important to note that the full extent of the radial distance $r^*$ need not be incorporated in the linearization process. When considering the deposition on a chip which extends 1mm radially from a spray source located at a vertical distance of 100mm, this means the radial distance $r^*\leq0.01$ is of importance. Additional complexity may be introduced unnecessarily into the model leading to increased computational costs. The equations governing the electric field are used in the simplified model in order to determine the initial direction of the droplet. The main effect from the electrical force is seen directly as the droplet exits the nozzle. At this location, the droplet's initial speed and direction is decided.