Predicting layer thicknesses by numerical simulation for meniscus-guided coating of organic photovoltaics

To achieve maximum efficiency in organic photovoltaics (OPV), functional layers with uniform and exactly predefined thickness are required. An in-depth understanding of the coating process is therefore crucial for an accurate process control. In this paper, the meniscus-guided blade coating process, which is the most commonly used process for the manufacturing of organic electronics, is investigated by experimental and numerical methods. A computational fluid dynamics (CFD) model is created to simulate the coating behaviour of P3HT:O IDTBR, an industrial state-of-the-art active material system used in OPV, and its results' independence of numerical parameters is ensured. In particular, the influence of the coating velocity and the initially injected fluid volume on the resulting wet film thickness is studied. The developed CFD analysis is able to reproduce the experimental results with very high accuracy. It is found that the film thickness follows a power law dependence on the velocity (˜v2/3) and a linear dependence on the ink volume (˜V). Accordingly, an analytical expression based on our theoretical considerations is presented, which predicts the wet film thickness as a function of the coating velocity and the ink volume only based on easily accessible ink properties. Consequently, this CFD model can effectively substitute time-consuming and expensive experiments, which currently have to be performed manually in the laboratory for a multitude of novel material systems, and thus supports highly accelerated material research. Moreover, the results of this work can be used to achieve homogeneous large-area coatings by utilising accelerated blade coating.


Introduction
Organic photovoltaics (OPV) provide a huge potential to unlock novel applications for photovoltaics, because of their unique characteristics, for instance lightweight, mechanical flexibility, and semi-transparency (Hu et al., 2022).In lab experiments, OPV maximum efficiencies above 19% have been reported (Zhan et al., 2022;Zhu et al., 2022).The current world record efficiency of an OPV module with an area > 200 cm 2 is 11.7% (Distler et al., 2021).To further improve the device efficiency, an immense research effort is being invested, not only concerning practical innovations and optimisations, but also the utilisation of numerical simulations to gain a better understanding of the optoelectronic properties of printed thin-film photovoltaics and the underlying working mechanisms (Aeberhard et al., 2019;Chen et al., 2019;De Falco et al., 2012;Zandi & Razaghi, 2019), to optimise the device layout (Lo Sciuto et al., 2016;Lucera CONTACT Fabian Gumpert fabian.gumpert@th-nuernberg.deFabian Gumpert, Keßlerplatz 12, 90489 Nuremberg, Bavaria, Germany et al., 2015), or to explore the potential of performanceenhancing features (Lo Sciuto et al., 2020;Mann & Rastogi, 2020).
A thin-film solar cell consists of different layers with different functionalities, which in combination generate electrical energy from the absorbed light.In the case of printed photovoltaics, such as OPV, these layers are deposited from solution, i.e. an ink that comprises the functional material dissolved in a solvent.Upon upscaling from small-area lab cells to large-area modules, a drop in power conversion efficiency of OPV devices is observed.One major reason for this is related to the production process and the homogeneity of the active layer, which becomes harder to maintain with increasing device area (Ibrahim & Hassan, 2019).For solution-processed organic electronics, like OPV, meniscus-guided coating techniques are most commonly used to deposit the functional layers (Diao et al., 2014).For industrial scale applications with high throughputs, roll-to-roll (R2R) slot-die coating is the most relevant approach among these (Gu et al., 2017).For early-stage research, a very similar technique called blade coating (aka 'horizontal dipping', see Figure 1) is used because of its minimal material consumption.
For blade coating, an applicator (red) is positioned over a substrate (e.g.glass) at a predefined distance, and a defined volume of ink is injected into the gap between the applicator and the substrate.The applicator (or the substrate) is then moved at a defined speed in one direction.While most of the fluid remains inside the gap between the applicator and the substrate, some of the fluid is deposited onto the substrate to leave a thin liquid film behind.A predefined uniform layer thickness is required to achieve maximum efficiency in an OPV device.However, the complex interplay of many different parameters (coating speed, fluid volume, gap height, temperature, fluid kinematic viscosity, fluid surface tension. . . ) influences the process and makes it difficult to predict the resulting coating liquid film thickness for novel inks based on new material combinations.
On one hand, there are some practical studies reported, which investigate the meniscus-guided coating of organic electronic materials and the influence of coating parameters on the properties of the deposited thin film (Chen et al., 2019;Lu et al., 2020).On the other hand, there are also a few simulative works, which numerically model the film formation process itself (Iliopoulos & Scriven, 2005;Singh & Ormiston, 2021) or the morphology formation upon drying of the solvent (Michels et al., 2020;Ronsin & Harting, 2022;Zhang & Wang, 2021), in order to contribute to a better process understanding of meniscus-guided coating techniques.
Time-consuming and expensive experiments need to be performed to experimentally establish correlations between process parameters and the resulting film thickness for a specific material system.In order to overcome this difficulty, we developed a computational fluid dynamics (CFD) model based on finite element analyses (FEM), which simulates the blade coating process.The model acts as a 'digital twin' of the process that can predict OPV active layer thicknesses as a function of fluid properties and process parameters.In this work, we describe the development of this first-of-its-kind model and validate it by comparing the results of simulation and experiment.An analytical expression is derived that allows for the prediction of the wet film thickness as a function of the initially injected fluid volume and the coating velocity of the applicator.

CFD model
For the numerical modelling of the coating process, a CFD analysis is conducted.The simulations are performed in the COMSOL Multiphysics R platform, which utilises FEM to solve differential equations.For the CFD analysis of the coating process, a moving mesh modelling approach is chosen, since it is superior compared to field-based methods in terms of computation time.However, the moving mesh method cannot handle topological changes.This limitation is overcome by the introduction of an artificial area in the simulation domain (see chapter 3.2.1).Two basic equations are solved in a coupled manner.The first equation is the Navier-Stokes equation: Here, u is the velocity vector (SI unit: m/s), ρ is the density (SI unit: kg/m 3 ), T is the stress tensor of the fluid (SI unit: Pa), and F is the volume force vector (SI unit: N/m 3 ).
The stress tensor is defined as T = [−pI + K], where p refers to the pressure (SI unit: Pa), I to the 3x3 identity matrix, and K to the viscous stress tensor (SI unit: Pa).
The second equation is the continuity equation: Since the Mach number is very small (MA < < 0.3), an incompressible flow can be assumed.Thus, the continuity equation can be rewritten as The fluid at the applicator surface is assumed to be stationary in the model.A no slip boundary condition is implemented for boundaries 7, 8, and 9 (see Figure 3) To simulate the movement of the applicator during the coating process, boundaries 2, 3, and 4 in Figure 3 are set as Navier slip boundary conditions.This condition prescribes a no-penetration at the boundary with n being the unit normal vector of the boundary.An additional friction term F fr at the fluid-solid interface can be defined as Here, γ is a numerical parameter, which depends on the size of the mesh elements, μ is the dynamic viscosity (SI unit: Pa s), and u slip is the tangential velocity of the fluid at the boundary, which is defined as where u tang regards to the tangential velocity of the fluid which can be written as u tang = u − (u • n)n, u wall to the magnitude of the wall movement which equals the coating velocity, and e x to a unit vector, pointing parallel to the boundary in the positive x direction.In Figure 3, the boundaries 1, 5, and 6 represent the system boundaries and are set up with open boundary conditions with f 0 being the magnitude of the stress in normal direction (SI unit: N/m 2 ) at the boundary.For the CFD analysis, zero normal stress is assumed.The surface of the coating fluid is in contact with the surrounding air, making a fluid-fluid interface boundary condition for the boundaries of the menisci (10, 11, and 13 in Figure 3) necessary.
On the one hand, the fluid-fluid interface boundary condition connects the velocities in both phases, air and coating ink.Since no mass transfer across the interface is assumed, the continuity of both velocity fields across the interface can be expressed by where u air and u ink represent the velocity fields in air and ink, respectively.On the other hand, the stress difference across the fluid-fluid interface can be written as with n air being the unit normal vector at the interface pointing into the air domain, ∇ s being the surface gradient operator (∇ s = ∇ − n(n • ∇), and σ being the surface tension coefficient (SI unit: N/m) (COMSOL Multiphysics R , n.d.).The fluid-fluid interface boundary condition assumes a contact angle for the intersection points of the interface with solid surfaces.
The contact angles of the air-ink interface with the applicator (α) and the glass substrate (β) are experimentally determined (see Table 1).The contact angle at the intersection point of boundaries 5, 6, and 13 is assumed to be 90 • .

Coating liquid
For all investigations, an OPV ink based on the photoactive material combination of the donor polymer poly(3hexylthiophene) (P3HT) and the non-fullerene acceptor (5Z,5'Z)-5,5'- ((7,7'-(4,4,9,9-tetraoctyl-4 was selected, a material system that is of very high relevance for large-area OPV and its industrial application due to processing and stability reasons (Armin et al., 2021;Strohm et al., 2018).The ink consists of the solute P3HT:O-IDTBR in a 1:1.25 (w/w) ratio, which is dissolved in a solvent mixture of o-xylene and 1methylnaphthalene in a 19:1 (v/v) ratio.The solid content of the ink is 36 mg/ml and the resulting ink density is 0.92 g/cm 3 .After depositing a thin film of the OPV ink on a glass substrate by blade coating (using a Proceq/Zehntner ZAA 2300 thin film coater), the formed wet film dries quickly and leaves a P3HT:O-IDTBR dry film behind (i.e. the active layer of an organic solar cell).After measuring the dry film thickness (in our case with a confocal microscope), the respective original wet film thickness can be calculated.This method is used in this work to determine the experimentally obtained wet film thicknesses.
In addition, the dynamic viscosity μ and the surface tension σ of the ink, as well as its contact angles with the blade applicator (α), and the glass substrate (β) are measured experimentally.All fluid-related properties that are used as input data for the simulations are summarised in Table 1.

Experimental motivation
For blade coating of low-viscosity inks (as being used in OPV manufacturing), our experimental results show that the deposited wet film thickness is a function of coating speed and injected ink volume (see Figure 2).On one hand, the film thickness strongly increases with increasing velocity.On the other hand, the film thickness also significantly increases with increasing ink volume, which has an important implication: Upon coating, the ink volume inside the gap between substrate and applicator steadily decreases (because part of the fluid is deposited on the substrate), which means that coating over longer distances at constant speed will inevitably lead to a thickness gradient, i.e. a decreasing film thickness along the coating direction.However, if there were a mathematical expression that exactly described the film thickness as a function of coating speed and ink volume, it would potentially be possible to compensate this thickness gradient caused by the ink loss with a defined acceleration of the blade movement, in order to obtain films with constant thickness over large distances.In the following, such a theoretical correlation will be established based on simulations and verified by experiments.

Design of the modelling domain
Figure 3 shows the modelling domain of the CFD analysis.Regions A and C are air domains, whereas regions B and D represent the domains of the coating liquid.The trapezoid shape on the top stems from the applicator.Boundaries 7, 8, and 9 depict the interface of the applicator.The surface of the substrate on which the coating liquid is applied is represented by the boundaries 2, 3, and 4.
For the investigation, the geometry of the model is developed in a two-step procedure.This improves the solvability of the problem and reduces its computing time.The initial shape of the air-fluid interface is assumed to be circular.Different centres and radii of the underlying circles correspond to different injected volumes (50-90 μl).The exact shape of the menisci (boundaries 10 and 11 in Figure 3) is then modelled by the quasi-stationary study in the first step.In the second step, the geometry is supplemented with an artificial region (D in Figure 3) to avoid topology changes of the model and maintain a fine enough discretisation.

Numerical discretisation
During the transient simulation, the geometry will experience a deformation.A 'roller' feature is applied to all boundaries where the air-liquid interface intersects with the ambient boundaries of the simulation domains (boundaries 2-9) to allow this modification.This feature allows the discretisation elements to move tangentially along the boundary and prevents mesh elements from moving perpendicular to the boundary.For the discretisation of the simulation domain, particular emphasis was put on meshing domain D since film formation takes place in this area.The mesh that is used for the investigation of the wet film thickness is shown in Figure 4.
A conventional triangular discretisation type is chosen for the air domains (A and C) and for the domain of the injected fluid volume (domain B).The geometry of domain D will undergo significant changes during the transient study focussing on the wet film thickness.The height of domain D will decrease during the study, along with the liquid film thickness.Using triangular elements, the shrinking of this domain in y-direction would cause deformed mesh elements.Resulting from this suboptimal discretisation, only simulation results with poor spatial resolution and/or challenging numerical convergence behaviour would be possible.As a better approach, distribution features are set up for boundaries 5 and 12 to define a fixed number of mesh elements at these boundaries.Afterwards, a mapped mesh is applied to domain D, resulting in a rectangular discretisation (inset in Figure 4).To quantify the impact of the numerical discretisation on the simulation results, a mesh refinement study is conducted.This means that multiple studies were performed with different numbers of mesh elements along boundaries 5 and 12.The considered injected fluid volume was 90 μl and the coating velocity 90 mm/s, since these parameters result in the largest film thicknesses.Thus, dependencies of the simulation results on the numerical discretisation should be most distinct.In Figure 5, the result of the mesh refinement study is shown.
Since the size of the mesh elements in domain D has been varied, the number of elements has increased, which corresponds to more degrees of freedom (DOF).The DOF is a typical measure of the magnitude of a numerical problem; generally, the smaller the DOF, the less computing time is necessary.The linear dependency of the computing time on the DOF can be seen in Figure 5.
In our specific case, the simulation runs on an Intel(R) Xeon(R) Gold 6142 CPU @ 2.60GHz machine with 128 GB of RAM.For a specific discretisation, the liquid film thickness depends highly on the DOF until ˜150k DOF.At larger values, the film thickness is relatively constant: it varies by only 0.1 μm, which is a 0.3% relative variation of the total film thickness.
In order to efficiently investigate all further wet film thicknesses, a discretisation of ˜150k DOF is chosen, as it results in the same wet film thickness as simulations with finer discretisation, while having a reduced computational time.

Comparison of simulation and experimental results
In order to validate the developed CFD model, the simulated results are compared to the experimental results of the blade coating process.To this end, firstly, crosssectional photographs of the OPV ink in the gap between the applicator and the glass substrate were taken at rest and at the maximal coating speed of 90 mm/s (see Figure 6(a ,b)).The reddish trapezoid in all images is the lower part of the applicator, which holds the coating liquid (black).At rest (a)), the coating liquid is distributed symmetrically with respect to the centre of the applicator, forming equal menisci on both sides.In  the case of a moving applicator (b)), the fluid is no longer symmetrically centred underneath the applicator, but is rather dragged behind it due to the friction at the liquid-substrate interface.Consequently, the higher the coating speed, the greater the difference in fluid geometry between the idle state and the steady state of movement.
The respective simulated counterparts of the images are shown in Figure 6(c ,d), in which the liquid domains are indicated in black.In both cases, i.e. the idle (c)) and the moving (d)) state, the simulated menisci of the coating fluid match the experimental observations very well.In Figure 6(c), where the applicator is at rest, the supplemented artificial domain D from the geometry initialisation procedure, as described in Section 3.2.1, can be seen (indicated by the blue arrow).It is important to note that, in this state (i.e.before the coating started), the height of this area does not yet represent the experimentally observed wet film thickness at this position, because there is no deposited wet film yet.However, for the simulation of a moving blade, as depicted in 6(d), this area will change and adapt to the steady-state condition for the respective velocity, here 90 mm/s.In this case, the height of the artificial domain equals the deposited wet film thickness that is predicted by the simulation (indicated by the orange arrow).
The wet film thickness, generated by the blade coating process, is investigated in a transient, two-dimensional CFD study.In the beginning of the transient study, the liquid film thickness oscillates due to the rapid height reduction of domain D. Therefore, the time domain of the simulation is chosen in such a way that a steady state and a stable coating process are established, which is the case after a coating distance of 20 mm.Thus, the time domain depends on the coating velocity.Interestingly, also in case of the experiment, stable coating with homogeneous film thicknesses only starts after 10-20 mm coating distance, depending on the coating speed.
To compare the wet film thickness values obtained by simulation and experiment, a variation matrix of different coating speeds ranging from 1 to 90 mm/s and different volumes of injected ink ranging from 50 to 90 μl is chosen and investigated.The experimental wet film An excellent agreement between simulation and experimental results is observed for injected volumes of 50 to 80 μl.For 90 μl, the experimental and simulation data match well for small coating velocities, while they slightly deviate at the higher velocities.This mismatch is due to experimental values that are lower than the simulated ones, which is most probably caused by the applicator's inability to properly drag the injected fluid at the start of the coating process for high fluid volumes in combination with high accelerations.Consequently, an inordinately large amount of the coating fluid remains at the starting position, which has the same effect on the resulting layer thickness as an initially smaller injected volume V.Both simulation and experiment show the same trend of increasing wet film thickness with increasing coating speed.Furthermore, the increase in wet film thickness with higher injected ink volumes predicted by the simulations matches very well the observed experimental results.

Comparison of simulation and experiment with theory
In 1940s, Landau and Levich (Landau & Levich, 1988) and Derjaguin (Derjaguin, 1993) developed the theory to describe the deposition of a thin liquid film on a solid.The theoretical considerations hold only true for Newtonian fluids with a capillary number C a << 1 (so-called Landau-Levich-regime).In this case, the wet film thickness h as a function of the coating velocity v can be described by the power-law, derived by Landau and Levich In previous literature reports, the predicted powerlaw dependency of the wet film thickness versus the coating speed has already been experimentally verified for blade coating of OPV inks (Nickel et al., 2012;Park & Han, 2009).However, the dependency of Equation ( 11) on the injected volume V has not been studied directly yet.Since the deposited wet film thickness is determined by the radius of the formed downstream meniscus, its dependence on the injected ink volume depends itself on the applicator geometry.For example, Park et al. derived an analytical expression for the wet film thickness for a cylindrical applicator (i.e. a circular cross-section) (Park & Han, 2009).In our setup, which is a commonly used blade coater for OPV devices, the applicator has a trapezoidal cross-section.In order to determine the wet film thickness as a function of ink volume, Equation ( 11) can be written as: where a(V) is a constant (SI unit: m 1/3 s 2/3 ), which depends on the injected fluid volume V. Accordingly, Equation ( 12) is used to fit our simulated and experimental data to the theoretical predictions.The results are shown in Figure 8(a,b), respectively.Remarkably, both the simulation and the experiment match very well with the theoretical correlation.In order to further determine the dependency on the injected ink volume, the values of the coefficient a, which correspond to the fits of the simulation and experimental data, are plotted versus the injected volume V in Figure 9. Again, the experimentally and simulatively determined coefficients a are in great agreement, with only a mismatch at 90 μl (experimental data point too low), because of the previously described insufficient ability of the applicator to drag the high amount of injected fluid at the start of the process Furthermore, the values of a suggest a linear dependency on the injected volume V, which is confirmed by a linear fit to the simulation data (dashed blue line in Figure 9).Consequently, a can be written as a(V) = b + cV, where b refers to the intersection and c to the slope,  respectively.The values corresponding to the linear fit in Figure 9 are listed in Table 2. Note: The fact that for V = 0, a(V) is > 0, as implied by the positive intersection value b, suggests that this correlation cannot be valid for very small volumes.Due to the trapezoidal applicator geometry, there is also a practical minimum volume that always has to be present between the applicator and the substrate (filling the gap underneath boundary 8 in Figure 3) in order to cover the whole substrate width upon coating, which is ˜23 μl in our case (gap height (400 μm) * substrate width (25 mm) * applicator bar length (2.3 mm)).Therefore, the validity range of the correlation is also limited to V > 23 μl.
To the best of our knowledge, the expression of the film thickness h as a function of V is a novel establishment.This correlation is certainly determined by the applicator geometry as well as the contact angle at the fluid-vapourapplicator interface, since both have an effect on the shape of the down-stream meniscus.In our case, the contact angle does not vary for the different injected volumes, due to the trapezoidal geometry of the applicator.Thus, the determination of the coefficients b and c enables an analytical expression to predict the wet film thickness h for an ink, as a function of the coating velocity v and the injected volume V: With our developed CFD simulation model, we have shown to be able to precisely determine these coefficients only based on the ink properties density, viscosity, surface tension, and contact angles.Consequently, this approach is able to predict wet film thicknesses and to suggest coating parameters for desired wet film thicknesses, which saves a lot of temporal and economical expenditure, when it comes to the qualification and process development of new ink formulations, and thus enables highly accelerated material research.Moreover, based on Equation ( 13) and the corresponding parameters in Table 2, the deposited wet film thickness can be calculated for every velocity-volume combination.This allows for the determination of an acceleration profile that is able to compensate for the loss of ink volume during coating, in order to prevent a thickness gradient as described in Section 3.1.Consequently, this approach enables the fabrication of large-area coatings with constant film thickness.

Conclusions
The CFD analysis of the blade coating process developed and presented in this work shows an excellent agreement with the experimental data and is able to precisely predict the wet film thickness as a function of coating velocity and ink volume only based on easily accessible ink properties.This approach has been demonstrated and verified for one OPV active material ink and will be extended to further material systems with different fluid properties in the future, in order to prove the general applicability of our simulations to model the coating behaviour of novel inks.By this, our CFD model can effectively substitute time-consuming and expensive experiments to establish such correlations for new material systems and thus constitutes a resource-efficient approach that supports the necessary high-throughput optimisations in a rapidly evolving technology field such as OPV.Moreover, the development of an accelerated blade coating process based on the results of this work to achieve homogeneous coatings over large distances, is currently ongoing and will be published in due course.Finally, our model is not limited to OPV or other types of printed electronics, but can be utilised for blade coated films of all kind, which require a precise and homogeneous layer thickness.

Disclosure statement
No potential conflict of interest was reported by the author(s).

Figure 1 .
Figure 1.3D illustration of the blade coating process.The 2D plane for the simplified CFD model is highlighted.

Figure 2 .
Figure 2. Contour plot of experimentally determined wet film thickness of blade coated P3HT:O-IDTBR layers as a function of coating velocity v and injected ink volume V.The plot comprises a total of 40 data points (a variation matrix of eight different coating speeds and five different ink volumes) and each data point is the average of two measured samples coated with identical parameters.

Figure 3 .
Figure 3. Computation region for the CFD simulation for an injected volume of 90 μl.For the sake of clarity, individual areas and boundaries are labelled by capital letters (A-D) and numbers (1-13), respectively.

Figure 4 .
Figure 4. Numerical discretisation of the simulation domain with an inset to highlight domain D.

Figure 5 .
Figure 5. Influence of the degrees of freedom on the liquid film thickness (indicated by blue squares) and the computing time (plotted as orange circles).

Figure 6 .
Figure 6.Cross-sectional view of the blade coating process with 90 μl injected fluid volume.Figures (a) and (b) show photographs taken from the experiments.The red trapezoid is the applicator, the black area is the coating ink, and the remaining two areas (left and right) are air.The substrate surface equals the bottom edge of the photographs.In each image, the coating direction is indicated by an arrow labelled with the corresponding velocity.Images (c) and (d) show the simulated counterparts of the experiments.In (a) and (c), the applicator is at rest (0 mm/s); (b) and (d) show the steady state for an applicator movement of 90 mm/s.In images (c) and (d), the artificial simulation domain (domain D) is indicated by coloured arrows.

Figure 7 .
Figure 7. Simulated and experimentally determined wet film thickness as a function of coating speed for different injected fluid volumes.Simulation results are plotted as solid lines, whereas experimental data is represented by squares with corresponding error bars.

Figure 8 .
Figure 8.In (a), simulated wet film thicknesses for different injected volumes are displayed by solid lines, the regarding fit of Equation (12) is shown by dashed lines.In (b), the experimentally determined wet film thicknesses for different injected volumes are indicated by square data points with associated error bars, while the regarding power-law fits are plotted as dashed lines.

Figure 9 .
Figure 9.The power-law coefficient a(V) from Equation (12) as extracted from the fits to the simulated (blue) and experimental data (orange) data with corresponding error bars as a function of the injected volume V.The blue dashed line represents a linear fit to the data obtained from simulation.

Table 2 .
Parameters of the linear fit to the simulatively determined coefficient a.