Two-layer modeling of thermally induced Bénard convection in thin liquid films: Volume of fluid approach vs thin-film model

Fabrication of micro- and nano-sized structures is an innovative solution for microelectronics, optoelectronics, and micro/nano-fluidic applications. A variety of techniques has been employed to fabricate micrometer- and submicrometer-sized patterns to the surface of thin-film (TF) layers, such as photolithography, electron-beam lithography, and nanoimprint lithography.^1,2 Although these methods are widely used to manufacture combined circuits in the electronics industry,³ these techniques’ multi-step fabrication process is both time- and cost-intensive. Furthermore, they consist of an etching step, which is an environmentally destructive and pollution-causing process. An impressive amount of work has been done to find an alternative method to address limitations associated with the conventional micro/nano-manufacturing.^4–8 These methods include hot embossing,⁹ drop-on-demand and printing technologies,^10–12 and lithography-induced self-assembly (LISA) of thin liquid films.^7,13–15 The latest method is typically driven by an external force such as electrical,^13,16–19 mechanical,⁷ intermolecular,^20–23 thermal,^24–28 or their combination.^29–32 One of the remarkable advantages of LISA is its ability to create an optically smooth surface due to the limited physical contact between the thin film and master mask, which is of interest in optoelectronic applications. Despite the simplicity in applying LISA to create micrometer- and nano-sized features compared to the conventional soft lithography techniques, the non-linearity and complexity of the governing mechanisms are not trivial as it relates to the growth of instabilities and the pattern formation. Hence, a significant amount of effort is required to understand the underlying mechanism of such instabilities and to harness them to create well-ordered features with consistent shape and sizes.^{14,19,31–34}

A typical LISA system consists of a molten polymer thin film, or polymer solution, on a substrate. In the films with micrometer- and submicrometer-sized thickness, the intermolecular interactions are relatively strong, and the film may undergo interfacial deformations, leading to the creation of holes or drops. Under such conditions, the material properties of the film, substrate, and bounding layer define the dominant intermolecular forces,³⁵ dewetting mechanisms, and consequently the shape and size of structures.^23,36 In the electrically induced patterning, the interface deformation is triggered by applying an electric field, which results in the application of net electrostatic force on the interface.³⁷ Similarly, there are instabilities caused by temperature gradients in the fluid layer at the interface, which is referred to as Bénard–Marangoni (BM) instability. If the temperature gradient at the interface results in non-uniformity of the surface tension at the interface, the BM instability is referred to as the thermocapillary (TC) effect.^38–47 The TC-induced instabilities are categorized into two fundamental modes: short-wavelength (SW) and long-wavelength (LW). The former leads to a cellular convection pattern without any nonlinear surface deformations, and it occurs commonly in relatively thicker or less viscous liquid films or a combination of both. In contrast, the LW mode is accompanied by a surface deformation appearing in thinner and more viscous liquid films.⁴⁸ Most of the previous studies that investigated the TC-induced instabilities of heated liquid films have focused on thermal convection⁴⁹ resulting from a density gradient across the film. As such, the SW mode of instabilities led to cellular convection structures with small height-to-width ratios.^50,51

Early studies on the TC-induced patterning of thin films are primarily devoted to exploring the corresponding mechanism for interface deformation. In the TC-induced patterning processes, the solid thin layer of film is initially heated to above its glass transition temperature (T_g). Then, the transverse temperature gradient is applied using the controlled cooling of the film from the top plate. Chou and Zhuang⁵² were the pioneers who reported the creation of micrometer-sized features on the ultrathin molten polymer films using high transverse thermal gradients. They introduced the surface charge (SC) model, based on the image charge-induced electrohydrodynamic instabilities, as the underlying mechanism for TC-induced instability. Soon later, Schäffer et al.^7,53 reported the formation of pillars after the grounding of substrates and liquid films. They introduced the acoustic phonon (AP) model as an alternative to SC. According to the AP model, the low-frequency acoustic phonons generate significant destabilizing radiation pressure. The AP model was well accepted in the scientific community until recent studies by Dietzel and Troian¹⁴ and McLeod et al.,¹⁵ who proposed the TC model to predict the spontaneous formation of periodic pillar arrays on the molten nanofilms exposed to a transverse thermal gradient. In contrast to SC and AP models, the TC model is applicable to any nanofilm and represents an extreme range of BM instabilities for LW fluctuations, where gravity forces are absent, and TC forces overcome the stabilization by capillary forces.¹⁵ While the previous efforts were limited to the AP and TC prediction of instabilities, based on the initial stages of pattern formation, we recently studied the non-linear stages of pattern formation for both AP and TC models under the long-wave approximation limit.⁵⁴ 

Nevertheless, the previous studies are mainly limited to simplified governing equations using long-wave approximation either in linearized or non-linear modes.^54–56 A majority of these studies utilized an analytical approach to evaluate the length and time scales in the TC-induced patterning process.^14,33 This non-linear model, called the “thin-film” (TF) equation, relies on the Stokes flow and the lubrication approximation. The linearized predictions can be obtained by replacing the interface height with the sinusoidal perturbations and ignoring non-linear terms.^14,32,54 However, neglecting the higher-order non-linear terms is a critical drawback, preventing effective prediction of the dynamic process of the thermal-induced patterning. A non-linear numerical model with high accuracy is required to address this challenge, which is the focus of this study.

The spatiotemporal evolution of the thin liquid film includes multiphase flow interface tracking, which is commonly modeled by computational fluid dynamics (CFD). These models track the interface between two phases and primarily include volume of fluid (VOF), level set, and phase field models.^57–61 In particular, the VOF method has attracted significant interest owing to its remarkable accuracy in the modeling of two-phase flow systems with a sharp interface. This method has been mostly implemented in examining the hydrodynamics of bubbles rising in a liquid.^62,63 Tomiyama et al.⁶⁴ proved the high accuracy of the VOF method in simulating the shapes and terminal velocities of a single rising bubble as compared to the published experimental data. Kawaji et al.⁶⁵ reported new findings for the two-dimensional simulation of a Taylor bubble rising in a stagnant liquid-filled tube. Given the high accuracy of the VOF method, it is selected in the present work to numerically investigate the thermally induced patterning of heated liquid films.

Here, we studied the dynamics and the interface evolution of nanofilms exposed to the transverse thermal gradient either using a flat or patterned top mask. The patterned top mask resembles the localized cooling of the thin film, which triggers the instability locally and induces the motion in the film. The characteristic wavelength was obtained by VOF, TF model, and linear stability (LS) analysis for many cases. First, we derived a correlation, which shows the relation between the characteristic wavelengths obtained by VOF and LS analysis. Moreover, we not only considered the effect of distance between the top and bottom substrates on the characteristic wavelength, but we also investigated the impact of the initial film thickness and temperature difference. This will advance the initial empirical model developed by Dietzel and Troian,¹⁴ which constitutes a novel contribution to the field. We studied the impact of the width of the top mask protrusions on the final shape of the structures, the number of pillars, and the processing time. It was demonstrated that, by smart adjustment of width and height of protrusions as well as the initial thickness of the polymer and the distance between the protrusions, well-organized pillars could be formed. We also investigated the formation of secondary pillars, originating from the residual thin film after the growth of primary pillars. The presented analyses aim to provide an insight into the mechanism of thermal-induced patterning and to save time and cost in designing related experiments.

A three-dimensional (3D) schematic of the ultra-thin liquid film with the mean initial thickness of h₀ is shown in Fig. 1. The liquid film and air bounding layer are placed between two plates that are heated from the bottom boundary (T_H) and cooled from their top boundary (T_C). It is assumed that the heat dissipation within the confining plates is negligible. The film is assumed to be flat and initially stable. The TC-induced instabilities are originated from imbalances at the temperature distribution along with the interface. This non-uniformity of temperature distribution either can be achieved using the patterned upper plate, inhomogeneity of the upper plate, or by imposing the random perturbations on the film interface. In the present study, the patterned upper plate is used to induce the TC instabilities. Using a patterned top plate, the non-uniform temperature distribution is achieved on the air–polymer interface, where the thickness of the air gap affects the cooling rate. Since the geometry has a long length in the y-direction, the variations in velocity and temperature in this direction are negligible; therefore, a 2D geometry is used to capture the physics. As shown in Fig. 1, d₁ and d₂ are the separation distances between the plates, where Δd = d₁ − d₂ is the protrusion height. The protrusion width is w, and the spacing distance between the neighbor protrusions is l_p. Although the viscosity of the molten polymer changes at temperatures above T_g, it is assumed that the ultra-thin film behaves as a Newtonian and incompressible fluid with constant viscosity evaluated at the mean working temperature: $μ=μTH+TC2$ . The variation in density due to temperature gradient within the film and during the interface deformation is included in our simulation using the Boussinesq approximation. However, it should be noted that the effect of microgravity condition and density variation is negligible compared to viscous forces in nanofilms.

The numerical analysis considered in this study is based on solving the fluid motion governing equations and continuity equation using the VOF method. The details on the VOF model and the computational setup are provided in Secs. II A 1 and II A 2.

This study numerically investigates the dynamics, instability, and morphological evolution of the air–liquid interface. The air–liquid system is considered a two-phase flow interface tracking problem, and the governing equations for two-phase flow are solved using the VOF method. ANSYS-Fluent is used as our base-solver, which has been extensively used for similar studies and it is known to provide accurate and stable results.^62,63 The calculations are accomplished using the pressure-based segregated and implicit solver. Pressure–velocity coupling is performed using the pressure implicit with the splitting of the operator (PISO), considering two corrections to account for neighbor and skewness. The momentum and energy equations are discretized using a second-order upwind differencing scheme, while the pressure staggering option (PRESTO) is used for pressure interpolation. The implicit scheme with a courant number of 0.2 and a time step of 10⁻⁶ is used to obtain convergence. No gravitational force is imposed on the simulation, and all cases are run in the double-precision mode. The VOF method is used to track the polymer–air interface utilizing a geometric reconstruction scheme based on the piecewise linear interface calculation (PLIC) method.⁶⁶ Moreover, the geo-reconstruction module is added to the existing VOF scheme to enhance the prediction accuracy of the polymer interface.^67,68 In the VOF model, a so-called fraction function of α is defined as the integral of fluid’s characteristic function on an Eulerian fixed grid cell, which is unity at any point occupied by the fluid and zero elsewhere. The cells with 0 < α < 1 state the partially filled case that indicate the interface.^69,70 The cutoff value for determining phases in the VOF model is 1⁻⁶. For a gas–liquid two phase flow system, the volume fraction of gas and liquid phases is α_g and α_l, respectively, and must satisfy the conservation of mass (with respect to the phase concentration) in each computational cell: α_g + α_l = 1.

In general, the VOF method is applied to track the interface movement between the immiscible liquid/gas and the liquid/liquid phases. A single momentum equation is solved throughout the domain to obtain a velocity field that is shared among the phases,

where p is the pressure and ρ and μ are the volume-averaged density and viscosity computed for each cell,

and u is treated as the mass-averaged variable, depending on the volume fraction of each phase α_l and α_g and their material properties,

The term F_V in Eq. (1) is the volumetric force at the air–polymer interface, resulting from the surface tension force. The volumetric surface tension force at cells containing both phases is evaluated using the continuum surface force (CFS) model,⁷¹ 

where σ is the coefficient of surface tension, n is the normal vector, which is calculated from the gradient of volume fraction, and k is the local surface curvature,

Equation (5) is added to the governing equations as a user-defined function (UDF). The surface tension varies by non-uniform temperature distribution at the gas–liquid interface, which leads to TC-induced motion.^23,26,56 Here, the focus is on the TC-induced instabilities, in which we only consider the variations in σ due to inhomogeneity in temperature distribution. We also considered a linear dependence of σ on the temperature T(x, z, t),^54,56

where σ_T is the temperature derivative of the surface tension, σ₀ is the surface tension evaluated at reference temperature T₀, and T is the liquid film temperature evaluated at the interface $T(x,z,t)z=h$⁠. The tracking of the air–film interface is accomplished by the solution of a volume fraction continuity equation for only one phase (the gas phase here),

In order to capture the TC effect, the temperature field is computed within the domain using the energy balance equation. The energy equation that is shared among the phases is given by

Here, E is the energy and can be computed as follows:

where H is the enthalpy. H can be calculated for incompressible flows as

Here, Y_j is the mass fraction of species j and

where T₀ = 298.15 K. The VOF treats energy (E) and temperature (T) as mass-averaged variables using the following definition:

The corresponding version of Eq. (13) is written for each phase based on the shared energy/temperature and the specific heat of that phase. The thermal conductivity is also shared by the phases similar to the density and viscosity defined in Eqs. (2) and (3).

A two-dimensional computational model is used to solve the governing equations [Eqs. (1), (8), and (9)] based on the finite volume method (FVM). The three-dimensional variations (gradients) along the span of the domain are expected to be minimal, and the system is expected to behave two-dimensional for the scales considered here.^56,61 The simulation domain in the Cartesian coordinate system is shown in Fig. 2, which consists of 2 subdomains (one for each phase), 15 boundaries, and 14 vertices. Domain-1 and domain-2 are initialized by the air and polymer film, respectively. The red and blue lines represent the applied hot and cold constant temperatures. The no-slip boundary condition is also applied to the top and bottom plates. It is necessary to have the periodic boundary conditions at edges 3, 4, 14, and 15. In all the simulations, a domain with a length that is greater than the characteristic wavelength (obtained from LS analysis) is chosen. This would be enough to capture the growth of instabilities when the periodic boundary conditions are employed; otherwise, fluctuations will dampen out over time.⁷² It is assumed that the liquid film, shown with the purple line at boundary 2, is initially flat. The boundary conditions are summarized in Table I. A non-homogenous grid is generated with a higher grid density near the walls (ratio of 1:1) to improve the tracking accuracy of the simulations close to the solid walls. This leads to a higher resolution for the dynamics and morphological evolution of the interface as the polymer touches the top plate or the bottom plate’s dewetting.

Table II shows the range of geometrical parameters and the material properties used in the TC simulation. The value is adopted from the experimental studies on the thermal-induced patterning in the literature.^14,55

The parameters, such as vertical and lateral coordinates, velocity field, temperature, and time, are normalized as follows:

The initial film thickness (h₀) is used to nondimensionalize the vertical dimensions of z, d₁, and d₂, whereas the characteristic wavelength (λ_c) is used to normalize the lateral dimensions of x, l_p, and w. Terms $w̃,l̃p,t̃,x̃,z̃,T̃$⁠, and $h̃$ denote the scaled parameters of the width of the protrusions, the distance between the protrusions, time, length, height, temperature, and interface height, respectively. The velocity and pressure fields are also scaled using the TC-induced characteristic velocity u_c defined as $uc=σTkrd2ΔTμ[1−kr+krd̃2]2$⁠.^14,54 The nondimensional term of $d̃2$ is referred to as the inverse filling ratio, and k_r is the ratio of the thermal conductivity of the polymer film to the air layer. Applying the above-mentioned scaling factors to the governing equations results in the following dimensionless ratios for Marangoni, capillary, and Prandtl numbers:

The governing equations (continuity, momentum, and energy equations) and boundary conditions (no slip at the top and bottom substrate and kinematic and stress balance at the interface) are normalized.^14,54 To find the characteristic wavelength relation, the shorter distance between two substrates (⁠$d̃2$⁠) is used (see Sec. S1 of the supplementary material). The normalized equations are then simplified using long-wavelength approximation (h₀/λ_c ≪ 1), leading to the following dimensional TF equation, which shows the spatiotemporal evolution of a thin liquid film subjected to the transverse thermal gradient (the derivations are available in the supplementary material, S2):

To derive the insight into the thin-film’s instabilities and dynamics, we solve the nonlinear TF equation that is the fourth order nonlinear partial differential equation (PDE). Using the finite difference (FD) discretization and applying the second order central finite difference scheme for spatial derivatives and a differential-algebraic system solver of DASSL with adaptive time step solver settings. The corresponding boundary conditions (BCs) for the thin-film equation are (a) no-slip BCs at the top and bottom substrates, (b) periodic BCs at the right and left walls, and (c) kinematic and stress balance BCs at the polymer/air interface.^14,54 More details about this numerical scheme are available in our previous studies.^32,54,73

Linear stability (LS) analysis is the conventional technique in predicting the characteristic wavelength for the growth of instabilities, λ_c. The corresponding wavelength is calculated by replacing the uniform interface height (h) in the so-called thin-film equation with a periodic perturbation of $h̃=1+ξ⁡exp(κ̄x̃i+Sκ̄t̃)$⁠, where $i=−1$⁠, $κ̄$ is the wavenumber, $Sκ̄$ is the growth coefficient, and ξ (≪1) is the infinitesimal amplitude coefficient. The details for the derivations of the dispersion relation and the characteristic wavelength are available in the supplementary material for the LS analysis section (S3). The following equation shows the characteristic wavelength obtained by LS analysis:

The characteristic wavelength obtained from LS analysis (λ_LS) will be used as a lateral normalizing factor. To show the dependencies of the characteristic wavelength on the distance between substrates, film thickness, and temperature gradient, the following equations can be derived based on Eq. (17):

Different spatial grid element sizes are examined to find the optimum grid size that accurately predicts the interface deformation. In this effort, the case is considered for the following conditions: d₁ = 200 nm, d₂ = 150 nm, w = 600 nm, l_p = 4000 nm, and ΔT = 100 K. Four cases (grid 1, grid 2, grid 3, and grid 4) are generated to verify the effect of the spatial grid on the results. Here, the number of elements is increased systematically from grid 1 (12 000) to grid 4 (24 000), while maintaining the same grid distribution and quality. The exact location of the interface for each mesh is compared between consecutive cases (e.g., grid 4 with grid 3 and grid 3 with grid 2), and the results are presented in Table III. The error in Table III was calculated as

where n is the number of extracted data points for interface location. This variation between the grids falls below 5% for grid 3, which indicates that it is sufficient to accurately track and predict the interface deformation for the geometry with the biggest size and most complex features in the present study.

The air bounding layer has a lower thermal conductivity than the liquid film. Therefore, the protrusion areas, which have a smaller air gap, are exposed to the higher temperature gradient and vice versa. The induced non-uniform temperature gradient results in a gradient of surface tension at the interface. The generated tangential stress transfer to the lower layers of fluid by viscosity triggers the fluid flow and subsequently leads to the film’s morphological evolution. Using a top plate with protrusions generates a temperature gradient along the polymer–air interface and adds more control to the pattern formation, resulting in well-ordered pillar arrays. The details of this observation will be discussed in Secs. IV A–IV C.

Many studies have explored the formation of hexagonal or random pillar arrays using a flat plate as a top substrate.^7,54 Here, we used a patterned cold plate to reduce the coarsening effect and generate well-ordered pillars. The temperature gradient at the interfaces and the resultant interfacial tension gradient are responsible for the TC-induced instabilities and pattern formation on the film.

The pattern formation steps and the film’s temperature distribution are presented over the evolution time in Figs. 3(a)–3(d). The deformation of the interface is divided into three stages. At the early stage of deformation, the air gap is large, and the temperature gradient and the resulting TC force at the interface are small, leading to a slow amplification of interface deformation. In the second stage, the interface has already deformed, and the air gap is smaller, which generates a higher thermal gradient at the interface. This leads at this stage to a faster deformation of the interface. The pattern evolution process continues until the liquid film touches the top plate in the last stage (⁠$t̃=403$⁠). The wetting mechanism of the top plate and dewetting of the bottom plate by the polymer film depend on the intermolecular interactions between the polymer material and the plates as well as the ratio of the initial film thickness to the plate separation distances.^23,74 These factors affect the transition in the shape of created pillars from cone to columnar and finally to the inverted cone-shaped structures. Tracking the interface height over time [shown in Fig. 3(d)] elucidates a non-linear increase in the rate of interface deformation as the interface moves toward the top plate. In this study, the contact angle that shows the tendency of the polymer film in the wetting of the plates is constant and set to π/2 for all cases.

The temperature distribution within the film and along the interface, including the interface evolution over time, is presented in Fig. 3(b). At the initial condition (⁠$t̃=0$⁠), the temperature in the film and the bounding layer is set to $T̃$ = 0.5, which is the average temperature of the top and bottom surfaces. The temperature profile along the interface shows that the temperature gradient on the interface increases as the film deforms and eventually reaches the top plate [Fig. 3(c)]. Thermal conductivity ratios of the film and the bounding layer (k_r) impact the rate of interface deformation and thus more significantly alter the pillar growth rate [shown in Fig. 3(e)]. The temperature distribution at the interface becomes uniform when the thermal conductivity of the polymer is higher. Therefore, the TC force becomes smaller, leading to the slower deformation of the interface. A similar case is numerically simulated using the TF equation, and the results are presented in Fig. 4. The snapshots of the interface height at different stages show a single pillar formation under the protrusion with the exponential growth rate in time.

The characteristic length for the growth of instabilities in the TC-induced patterning process is predicted using the LS analysis, which relies on the small deformation from the original state (base case). Experimental observations and monitoring of the interface deformation at its early stages showed that the LS predicted wavelength (λ_LS) matches well with the experimental results.^14,15,33,55 However, the interface deformation undergoes non-linear deformations soon after the initial undulation. This leads to uncertainties in the validity of λ_LS to predict the final form of the structures. Here, we used a numerical simulation to predict the TC-induced interface deformation at both early linear and late non-linear stages of the pattern formation process. In this section, the relationship between the VOF-predicted characteristic length scale and the one obtained using LS analysis is predicted.

The simulation results are first compared with the experiments of Dietzel and Troian¹⁴ to verify their accuracy. In the VOF model predictions, the characteristic wavelength (λ_c) is found after pattern formation via measuring the center to center distance of neighboring pillars. For the case considered by Dietzel and Troian¹⁴ (d₁ = 310 nm, d₂ = 210 nm, w = 6 μm, l_p = 500 nm, h₀ = 100 nm, and ΔT = 46 K), our VOF-predicted characteristic wavelength is 1.82 μm, which is in good agreement (∼10% error) with their experimental value of 2.04 μm. However, the value of λ_c obtained from the LS analysis [using Eq. (17)] is 4.52 μm, which expectedly deviates significantly (∼121% error) from the experimental observation. In addition, the characteristic wavelength obtained from TF is 4.2 μm, which is a better prediction compared to the LS analysis but far from the value obtained from the experimental data and VOF model. Table IV presents the characteristic wavelength obtained from VOF, TC, and TF. The VOF simulation showed the lowest deviation from experimental observation compared to the TF and LS predictions. Given the simplicity of LS theoretical analysis compared to other two numerical predictions, it is of interest to find a correlation between the LS predicted characteristic wavelength and the more accurate VOF model results.

To find the relationship between the characteristics wavelengths predicted by the VOF model (λ_VOF) and the LS analysis (λ_LS), 13 cases are simulated based on different initial film thicknesses (h₀), separation distances between the plates (d₁ and d₂), and temperature gradients (ΔT). The values of λ_VOF and λ_LS are then compared in Table V.

The VOF-predicted characteristic wavelengths in Table V are obtained by measuring the average center to center distances of pillars formed at the quasi-steady stage of pattern formation. The width of the protrusion is 6 μm (w = 6 μm), and the distance between two protrusions is 7 μm (l_p = 7 μm). Figure 5 shows the snapshot of the film spatiotemporal evolution for case 1 in Table V, showing different pattern formation stages. The interface deformation started under the edge of protrusion, where the interfacial tension (and temperature) has the highest gradient. Subsequently, the generated surface waves traveled toward the center. The initial undulations on the interface have led to the formation of the first set of pillars near the edges with the concurrent fluid depletion around these pillars [Figs. 5(a-ii) and 5(a-iii)]. The remnant fluid layer between the two pillars experienced the TC force due to its uneven surface, which resulted in the formation of second sets of pillars at t > 119 ms [Figs. 5(a-iv) and 5(a-v)]. The pillar formation process continued until the remained polymer was sufficient to form one pillar [Fig. 5(a-vi)]. After the pattern formation is completed, when all pillars touched the upper mask, as shown in Fig. 5(a-vi), the characteristic wavelength is measured by averaging the center of the created pillars.

The measured λ_VOF is plotted against the LS-predicted wavelength (λ_LS) in Fig. 6. The relationship between λ_VOF and λ_LS is determined using an implicit radical function dependence,

where c₁ = 1.5034, c₂ = −0.0728, and c₃ = −1.5959 are fitting coefficients.

Figure 7(a) shows the variation of the characteristic wavelength predicted by the VOF model (red symbols), TF model (blue symbols), and LS analysis (blue line). Figure 7(b) shows the Marangoni number (Ma) and the capillary inverse (Ca⁻¹) with inverse filling ratio (⁠$d̃2=d2/h0$⁠). The Marangoni number (Ma) represents the relative strength of the TC force to the viscous force, while the capillary inverse (Ca⁻¹) represents the relative strength of the surface tension force acting on the interface to the viscous force. At constant plate separation distances (d₁ and d₂), the initial film thickness is increased to study the effect of different filling ratios on the final structure of the pillars. As can be observed in Fig. 7(b), the Marangoni number decreased by increasing the filling ratio. An explanation is that increasing the filling ratio (thicker film) leads to a lower temperature gradient at the interface, decreasing the TC force and Marangoni number. Similarly, Ca⁻¹ decreases as the filling ratio increases.

The characteristic wavelength predicted by the LS analysis (λ_LS) increased sharply for small $d̃2$ up to a maximum at $d̃2=2.44$⁠. To calculate the maximum point, the derivative of λ_LS with respect to h₀ is set to zero (dλ/dh₀ = 0) in Eq. (17). Thus, a relation is obtained for $d̃2$ as a function of k_r: $d̃2=3(1−kr)/(−kr)$⁠. Since k_r is a material property, the extremum value for $d̃2$ can be adjusted based on the geometrical limitations with respect to the thermal properties of the film. In the current study, where the ratio of the thermal conductivity of polymer to that of air is k_r = 5.37, $d̃2$ is calculated as 2.44. The red and blue symbols in Fig. 7(a) show the variations in characteristic wavelength predicted by the VOF model and TF model, respectively. λ_VOF has a similar trend as λ_LS. It increased for small values of $d̃2$ and dropped rapidly after it reached a maximum point at $d̃2=2.5$⁠. It must be noted that although λ_VOF followed a similar parabolic trend as λ_LS, its magnitude is about three times smaller than λ_LS. On the other hand, λ_TF has a similar trend as λ_VOF; however, it reached a maximum point at $d̃2=2.44$ (like LS). The magnitude of λ_TF is less than λ_LS and higher than λ_VOF.

Next, the effect of the temperature gradient (ΔT in the range of 50–100 K) on the capillary inverse (Ca⁻¹) and the characteristic wavelength predicted by VOF, TF, and LS analysis is investigated, and the results are presented in Fig. 8. Ca⁻¹ decreased with an increase in the nondimensional temperature as it is inversely related to ΔT. As the temperature difference along the film interface increases, the resulting TC force applied to the interface increases, leading to a higher characteristic velocity within the film. The higher u_c (and shear rate) results in higher viscous resistance and thus lower Ca⁻¹.

All the characteristic wavelengths from the LS analysis (blue line), TF model (blue symbols), and VOF simulation (red symbols) decreased with an increase in the applied thermal gradient [presented in Fig. 8(a)]. From Eq. (17), λ_LS is proportional to the inverse of the square root of the temperature difference $(λLS∝ΔT−1/2)$⁠. With an increase in the temperature difference from 0.5 to 1, the reduction in characteristic wavelengths is about 27%, 33%, and 33% for VOF, TF, and LS methods, respectively. Given the large difference in lambda’s absolute values between these approaches, the predictions for the relative reductions are in a close match.

So far, the similarity and differences in predictions of the characteristic length using three methods, the LS analysis, the TF model, and the VOF model, have been discussed. In what follows, the effect of localized cooling using a patterned top plate on the pattern formation’s fidelity will be compared using both VOF and TF nonlinear simulations.

The size and periodicity of the protrusions need to be tailored to avoid secondary pillar formations when using plates with protrusions to generate the non-uniform thermal gradient and enable forming well-ordered patterns. In addition, the shape, size, and number of the formed pillars are affected by the geometry of the protrusions at the top cold substrate. Later in this section, we discuss these conditions and outline the range of geometrical parameters predicted by VOF numerical simulations and the TF models. To simplify the comparison, λ_LS is used for the lateral normalizing scale.

The size and number of created pillars in TC-induced patterning are affected by the protrusions’ width. The snapshot of the thin-film evolution for the case of $d̃1=3,d̃2=2,w̃=0.52,l̃p=0.26,andT̃=1$ is shown in Fig. 9. The polymer interface is assumed to be initially flat [Fig. 9(a-i)]. At $t̃=5644$ in Fig. 9(a-ii), two twin pillars formed under each protrusion. The separation distance between the formed pillars was long enough to prevent them from merging [Fig. 9(a-iii)] until they touch the top plate at $t̃≥14940$ in Fig. 9(a-iv). The interface deformation is strongly related to the interface temperature profile since it causes the non-uniformity in the interfacial tension. Figure 9(b) demonstrates the normalized temperature profile at the air–film interface over the deformation time. There exist two low-temperature regions close to $x̃∼0.4$ and $x̃∼1.5$ at the initial stages of interface deformation (⁠$t̃=5644$⁠). The surface tension at these regions is higher, pulling the polymer toward the neighboring areas with lower interfacial tensions. With the significant jump in the surface tension at the corners, the two twin pillars ascended at the initial stage of the pattern formation [Figs. 9(a-ii) and 9(a-iii)]. The temperature difference between the peak and trough at the interface increased as the pillar grew over time, which triggered the growth rate over time. This also prevented the merging of the neighboring pillars.

Figure 10 shows the snapshot of the thin-film evolution predicted by the TF model for the case of $d̃1=3,d̃2=2,w̃=1.47,l̃p=0.26,andT̃=1$⁠. The size and number of pillars in the TF model are similarly affected by the protrusions’ width. Similar to the VOF model, the thermal gradient subjected to the interface is higher at the corner of the protrusions, resulting in a significant gradient at the interfacial tension underneath protrusions. The width of the protrusion is large enough for the growth of instabilities and generating two pillars. However, the magnitude of $w̃$ for making two twin pillars is found to be higher than that in the VOF simulation.

An interesting observation in Figs. 9 and 10 is that by tuning the protrusion’s width, multiple pillars can be generated under each protrusion. This enhances the pattern miniaturization capability in TC patterning since the formed pillars’ width is less than half of the protrusion width. Depending on the width of the top plate protrusions, the number of created pillars varies.

Figure 11 shows the 2D snapshots of the interface over the evolution time for three different conditions. First, the case with $w̃<0.34$ leads to the formation of only one pillar under the protrusion, as shown in Figs. 11(a-i)–11(a-iii). This is followed by the transition stage, at which the protrusion width is $0.34<w̃<0.49$⁠. Here, two peaks are observed at the early stages. As the film evolved, however, the initially formed peaks merged, and only one pillar was formed at the end [Figs. 11(b-i)–11(b-iii)]. For $w̃=0.52$ in Figs. 11(c-i)–11(c-iii), the initially formed peaks are well separated from one another, and they did not merge over time to form two isolated pillars under one protrusion. Hence, the width of protrusions can be increased to form more pillars (e.g., five pillars in Fig. 5) under a single protrusion. However, using very large protrusions leads to the formation of secondary and tertiary pillars, whose size and shape will no longer be under control. Moreover, it should be noted that the threshold value of $w̃=0.52$ found here is for a specific filling ratio of 2. This threshold varies by changing the filling ratio, which will be discussed later in this section.

Figure 11 also highlights the streamlines for three different conditions. These conditions show that flow circulation takes place within both the polymer and the surrounding air. Viscous forces transfer the TC force at the interface to different fluid layers, which induces motion. The lateral flow of liquid at the surface takes place from the warmer regions to the cooler side under the protrusion, which forces the cooler liquid to descend. The downward movement of cooler liquid creates two convection cells, known as Bénard convection cells, for each interface undulation. The formation of two recirculation regions for each pillar matches well with experimental and numerical data reported by Vanhook et al.⁴³ 

A 2D map for both VOF and TF models is generated based on the normalized width of the protrusions (⁠$w̃$⁠) and the inverse filling ratios $(d̃2)$ to predict the threshold width value of $w̃*$⁠, in which the transition from single- to two-pillar formation occurred. The other parameters, such as temperature difference and thermal conductivities of the layers, are kept constant to generate this map. As shown in Fig. 12(a), there is a linear relationship between the inverse filling ratio and the threshold width. As the inverse filling ratio increased, the threshold width increased as well. The equation of the dashed line in Fig. 12(a) is $w̃VOF=0.2d̃ave−0.08$⁠, where $w̃$ is the scaled width (⁠$w̃VOF=w/λLS$⁠). The dimensional form of this equation [$wVOF=0.2d2h0λLS−0.08λLS$] can be used to obtain the width of patterns as a preliminary choice for fabricating patterned masks. On the other hand, a similar equation for the TF model is obtained [Fig. 12(b)], which is $w̃=0.32d̃ave−0.62$⁠, and its dimensional form is $wVOF=0.32d2h0λLS−0.62λLS$⁠.

The distance between two protrusions at the top mask is another geometrical property of the top plate that affects the fidelity of the pattern formation in TC-induced patterning. In all the cases discussed so far, no pillar formations are observed in the regions between each protrusion. Under certain conditions, secondary pillars can also form in the area between each protrusion. Here, we discuss the dynamics of interface evolution when the secondary pillars form. We will also look at the conditions required to avoid or mitigate the formation of such features.

The 2D time snapshots of the interface evolution for the case with $d̃1=2,d̃2=1.5,w̃=0.31,(l̃p)VOF=0.82,andT̃=1$ show the formation of the secondary pillar in the gap between two neighboring protrusions, as shown in Fig. 13. Like the previous cases, the thin film is initially flat [Fig. 13(a-i)]. The non-uniform temperature distribution along the interface destabilizes the liquid film under the protrusions [Fig. 13(a-ii)]. This leads to the formation of primary pillars beneath each protrusion at later stages [Fig. 13(a-iii)]. Once the primary pillars are formed, the residual polymer between these pillars, which no longer has a flat interface, moves toward the top plate due to the heterogeneous distribution of the interfacial tension [Fig. 13(a-iv)]. The residual polymer connected to the created structures through a concave meniscus and channel of evolution in this area starts with the film rupture at the meniscus. The TC force applied on the polymer/air interface near the meniscus results in the diffusion of fluid from thinner to thicker areas of the polymer, resulting in a film rupture at the meniscus. The thinning of the film between two crests continues until the interface touches the top plate, and the secondary pillar is formed [Fig. 13(a-v)]. The thickness of the remnant polymer after the formation of primary pillars, the magnitude of temperature gradient in those regions, and the wetting properties of both top and bottom plates affect the dynamics and the evolution of the secondary pillars.

Figure 14 shows the snapshots of the thin-film interface in the TC patterning obtained by the TF model. The formation of the primary and secondary pillars is similar to that of the VOF model. The primary pillar formed under each protrusion, and then, the secondary pillar grew at later stages until it touched the top substrates (longer distance between substrates).

Figure 15 reveals the effect of the distance between neighboring protrusions (⁠$l̃p$⁠) on the formation of secondary pillars. Figure 11(a-i) shows the case with $l̃p<0.8$⁠, where no secondary pillars are formed. For the distances in the range of $0.8<l̃p<1.4$⁠, only one secondary pillar was formed, while multiple secondary pillars were generated for the higher spacing distances (⁠$l̃p>1.4$⁠). Evidently, the secondary pillars should be avoided when the high fidelity of the pattern replication is of concern since there is no control on the shape and size of the features.

Similar to the map developed for the width of protrusions, we established a 2D map to find the threshold value for $l̃p$ that separates the patterns with primary and secondary pillars in both VOF and TF models. The effect of the distance between the protrusions on the number of pillars at different inverse filling ratios (⁠$d̃2$⁠) is presented in Fig. 16. In both models, by increasing the distance between two protrusions beyond a threshold (say $l̃p*$⁠), the secondary pillars appear. It is worth noting that even large spacing did not lead to the interface deformation for relatively high $d̃2$⁠. Moreover, only primary pillars are formed under the protrusions. In all these cases, the Ma number remained constant since it affects the predicted value of $l̃p*$⁠. In the VOF model, the equation of the dashed line in Fig. 16(a) is $(l̃p)VOF=0.39d̃2−0.04$⁠, where $(l̃p)VOF$ is the scaled distance between two protrusions [$(l̃p)VOF=(lp)VOF/λLS$]. By substituting the dimensionless parameters [$(lp)VOF=0.9d2h0λLS−0.36λLS$], the dimension size of the distance between two protrusions can be calculated and used as a preliminary choice for fabricating the patterned mask. On the other hand, a similar equation is obtained for the TF model, which is $(l̃p)VOF=0.85d̃2−0.02$⁠, and its dimensional form is $(lp)TF=0.85d2h0λLS−0.02λLS$⁠.

The volume of fluid (VOF) model, thin-film (TF) model, and linear stability (LS) analysis frameworks were used to simulate the thermocapillary (TC)-induced instability and morphological evolutions in thin liquid films. In the VOF model, a multiphase flow interface tracking approach and the full continuity, Navier–Stokes, energy, and volume fraction equations were solved. This also enables investigating more complex thin-film problems with multi-component films. The numerical approach presented in this study surpasses the long-wave limit inherent in the TF and LS models. It successfully predicts different interface deformation stages, including the initial linear and the late non-linear deformation. Compared to the experimental data,¹⁴ the VOF model predicted the characteristic wavelength (λ) with significantly higher accuracy compared to the TF model and LS analysis. Several cases with different polymer thicknesses, gap spacing, and temperature differences were investigated to correlate the characteristic wavelength predicted by the VOF model and TF model with that theoretically found by the LS formulation. The effect of localized cooling/heating was also investigated by using a patterned top cold plate. The results showed that the interface deformation and pillar formation are sensitive to the localized cooled area (w, under the protrusions). It leads to the transition from single to twin pillars. The 2D maps for both VOF and TF models are generated based on the normalized width of the protrusions (⁠$w̃$⁠) and their separation distance (⁠$l̃p$⁠), and the inverse filling ratio to predict the threshold values in which the transition from single- to two-pillar formation occurred and the secondary pillars formed.

The supplementary material includes the (i) vertical characteristic wavelength, (ii) thin-film equation, and (iii) linear stability (LS) analysis.

Financial support for this work through Canada’s Oil Sands Innovation Alliance (COSIA) and the Natural Sciences and Engineering Research Council of Canada is gratefully acknowledged (Grant No. CRDPJ 501857-16).

The data that support the findings of this study are available within the article.