CFD model of slot die coating for lithium-ion battery electrodes in 2D and 3D with load balanced dynamic mesh refinement enabled with a local-slip boundary condition in OpenFOAM

Slot die coating is a state-of-the-art process to manufacture lithium-ion battery electrodes with high accuracy and reproducibility, covering a wide range of process conditions and material systems. Common approaches to predict process windows are one-dimensional calculations with a limited expressiveness. A more detailed analysis can be performed using CFD simulations, which are often based on in-house code or closed-source software. In this study, a two-phase CFD model in two and three dimensions was created in OpenFOAM with the intent to provide a method for more detailed investigations of the slot die coating process with open access to source code and files. A custom boundary condition enables the proper description of the wetting behavior in the two-dimensional model. The combination of standard no-slip boundary conditions at the substrate boundary with the volume-of-fluid solution algorithm leads to a method-related air entrainment, which was prevented by allowing local slip at the dynamic wetting line at the upstream meniscus in the two-dimensional model. Additionally, a load-balancing dynamic refinement algorithm was implemented to minimize the computational effort and increase the ease of use of the simulation environment. The simulation was validated by comparing the simulated process limits to experimental observations, showing good agreement. As a result, this model enables detailed analyses regarding the influences of slot die geometries, material properties, and process parameters on the coating stability and wet-film profile.


Introduction
In state-of-the-art production lines, lithium-ion battery electrode slurries are coated onto a metal substrate with a slot die, which is embedded in a large-scale rollto-roll process chain. 1 Slot die coating is a proven method, which enables a reproducible application of thin films for a wide range of material systems and process conditions with high accuracy. As a premetered coating method, process parameters and the slot die dimensions determine the width and height of the evolving film, such as the flow rate of the coating fluid and web speed. The range of process parameters, within which the coated film is stable, is referred to as coating window. 2 When exceeding the range limits of this coating window, defects may disrupt the formation of a stable and homogeneous film. [3][4][5] Depending on the process parameters, different mechanisms cause the formation of defects. For excessive flow rates or insufficient coating speeds, the coating fluid may swell out of the coating gap against the coating direction. In contrast, when the flow rate falls below the lower limit, the fluid menisci, that bridge the substrate and die surfaces, may become unstable. Therefore, air is entrained into the film, or defects are caused in a downstream direction. 2 Another quality criterion is the homogeneity of the wet-film surface regarding the elevations of the film edges, which may lead to several complications in the downstream process 6 and can be interpreted as a defect as well.
Usually, simplified mathematical approaches estimate the coating window. These models are based on a one-dimensional flow analysis with semi-empirical correlations and assumptions, [7][8][9][10][11] leading to a restricted scope of validity with regard to the process conditions. One example is the ratio between the viscosity g, web speed u Web, and surface tension r, known as the capillary number.
Usual pressure boundary conditions at the downstream free surface assume Ca to be smaller than unity. 10,11 In the battery production area, these assumptions may not necessarily be justified, for example, due to higher coating speeds and characteristics of the material systems. In addition, the informational content is limited, as one-dimensional approaches do not provide full detail of velocity profiles, flow mechanical influences, fluid-fluid and fluid-solid surface interactions, interphase effects, interphase positions and shapes, and wetting behavior. Moreover, effects in cross-web direction are neglected and these approaches are limited to single-layer coatings and multilayer coatings with identical material properties. 12 Further optimizations of the slot die coating process require a more detailed analysis, for example, to investigate the influence of the inner die and die-lip geometry on flow characteristics, coating stability, and wet-film profile. Additionally, some defect mechanisms are not understood completely, such as the formation of edge elevations, representing a general challenge in the coating industry. Current research has shown that adjustments to the (inner) slot die geometry influence the formation of edge elevations. 13 This motivates the establishment of more advanced methods to investigate the complex interactions of the various influencing parameters, which have an effect on edge elevations. One way to meet these requirements is to implement a CFD-simulation model and deriving novel solution approaches.
Literature offers ample studies of CFD simulations for coating applications. In the case of two dimensions, several studies on the operability windows are available. [14][15][16][17] Further investigations on the formation of defects and the flow fields, in general, are provided by extending the simulation domain to three dimensions. [18][19][20] However, these models are not easy to implement as they are built with in-house code or commercial (closed-source) software (Fluent, FLOW3D, Comsol).
In this work, a multiphase CFD model of the slot die coating process in OpenFOAM is presented, which is entirely built on open-source software. The simulation includes two-and three-dimensional approaches with extended mesh functionality regarding a generation on CAD-basis with static and/or dynamic refinement methods, exceeding the scope of the standard Open-FOAM libraries. This allows local areas to be resolved at a higher level without a significant increase in computational effort, such as the interphase or the edges of the wet film. In addition, a custom boundary condition is included for proper handling of the method-related air entrainment at the dynamic wetting line. Conceivable applications for the two-dimensional model are investigations of flow phenomena in coating direction comprising an extended prediction of coating windows, investigations of the defect mechanisms, and the influence of material properties on the flow characteristics and wetting behavior. Furthermore, the three-dimensional model enables the extension of the study to the influences on the wet-film profile, edge effects, and other cross-web-directed mechanisms.

Experimental materials and methods
For validation purposes, several experiments were carried out using a development coater (TSE Troller AG). It consists of a steel roll with a high surface precision, which is coated using a slot die oriented in -30°position. On the opposite side, the coating is continuously removed with a doctor blade. In this way, influences on the coating process can be investigated on an infinite coating. The experimental setup and the dimensions of the slot die, which are used for the experiments, can be seen in Fig. 1. The coating stability is evaluated optically with a camera setup on the upand downstream side. The coating width was adjusted by changing the shim geometry.

Slurry formulation
In this study, material systems for lithium-ion battery anodes were used. The slurry consisted of the active material graphite, binders (carboxymethylcellulose, styrene-butadiene rubber), conductive carbon black, and water as a solvent. The associated concentrations The mixing of the slurry ingredients was performed with a dissolver Dispermat CN F2 (VMA-GETZ-MANN GmbH). After a dry mixing step to break agglomerates, solids were added stepwise and dispersed in a CMC solution. Subsequently, the mixture was dispersed for approximately 1 h and SBR-binder was added.

Rheology
In order to characterize the flow behavior of the coating fluid, the viscosity was measured in dependence on the shear rate between 0.01 and 3000 s -1 , using an Anton Paar MCR101 rheometer equipped with a plate-plate measurement head with a diameter of 25 mm. The measurements were repeated three times to calculate average values. For the CFD simulations, the Cross-Power-Law was used to approximate the rheology in dependence in the shear rate, as described in the following equation: The ratio between the viscosity m(c) at a specified shear rate c moves within the boundaries m 0 (zero shear rate) and m ¥ (infinite shear rate), approximated by a proportionality factor j and an exponent n. The results from the rheology measurements are visualized in Fig. 2 alongside the Power-Law by Ostwald/de Waele and Cross-Power-Law approximations. 21

Surface tension
The pendant-drop method was used to determine the surface tension with an EasyDrop measuring instrument (Krü ss GmbH) by analyzing an image of a drop under the influence of gravity. A syringe was filled with the coating fluid and a drop was created below the tube. The established drop shape was parametrized using the Krü ss Drop Shape Analysis Software to obtain the value of surface tension.

Contact angles
The static contact angles as a material property were measured with the sessile-drop method (EasyDrop, Krü ss GmbH), where a drop of liquid is placed onto a substrate and viewed from the side. In regard to slot die coating, the fluid bridge between the slot die and the substrate leads to a formation of four menisci, one at the up-and downstream lips and two lateral menisci (see Figs. 3 and 4). For the model buildup, the up-and downstream menisci and the associated contact angles are of great importance and are determined by using a camera-lens combination to obtain images of the coating gap in side view. The angles were determined by graphical analysis of the captured photos with an image processing software. While the contact angles at the die surface are static, the contact line at the substrate is in continuous motion, leading to a depen-

Numerical methods
The CFD simulation in the present work was performed with the OpenFOAM version v2012 using the two-phase volume-of-fluid algorithm interFoam. This solution method introduces the parameter a, which reflects the phase fraction in a grid cell, where a = 0 is pure gas and a = 1 is pure liquid. In each grid cell, the local viscosity and density are calculated according to the defined models from the locally present phase fraction. The resulting phase fraction is transported by convection, which results in the moving interphase at a = 0.5, 22 described by the following transport equation 3. In this model, the simulation is resolved in discrete time with pressure-velocity-coupling in form of a PIMPLE loop. 23 The simulation model is based on the development coater and the slot dies in the TFT laboratories. As simplification, the orientation of the gravitational vector is aligned with the slot and the substrate is approximated as a flat surface. An estimation of the curvature of the roll surface underneath the slot die shows that the influence is smaller than the experimental uncertainty for adjusting the gap height, which justifies this assumption. As a visualization, a rendered simulation result for exemplary process parameters (substrate speed u Web = 5 m min -1 , gap height h Gap = 180 lm, and a dimensionless gap height G* = 1.15) is shown in Fig. 3.

2D-model
First, a two-dimensional model was created to investigate the influence of the numerical method, wetting behavior, interphase/surface effects, and velocity profiles in the coating gap. Therefore, a section plane through the center of the film is considered, while cross-sectional effects are neglected.
For the mesh generation, the standard procedure of the blockMesh utility is used according to the slot die dimensions mentioned in the experimental section. Thereon, boundary conditions are assigned to the system boundaries, as illustrated in Fig. 4.
The boundary conditions for the fields of pressure (p), velocity (u), and the phase fraction (a) are presented in Table 2.

Dynamic wetting line
The contact line between the coating fluid and the substrate at the upstream meniscus is a key criterion for the coating stability. In general, the boundary condition at solid walls is the no-slip condition, which equates to the velocities of the fluid and the wall. In case of a moving contact line, this boundary condition is not the condition of choice, due to a method-related entrainment of air into the film. When applying a fixed velocity in coating direction, the average velocity in a boundary cell ū will also result in a positive contribution in coating direction, which then leads to a convective transport of phase fractions with a < 1 into the film. This mechanism is visualized in Fig. 5, where the gas phase (a = 0) is colored blue and the liquid phase (a = 1) is gray. The area in between (0 < a < 1) is marked by a color gradient.
The described entrainment of gas phase into the liquid film is non-physical and falsifies the calculations of material properties in the respective cells. Moreover, an influence on the surface properties is possible, which illustrates the relevance of the avoidance of the method-related air entrainment.
One possible way to bypass the method-related air entrainment is to allow local slip in proximity to the dynamic wetting line. This method was proposed by Kistler (1984), where the velocity at the moving wall follows a linear transition from the slip velocity to the web velocity over a defined distance. 18,24 This method is also known for the hydrodynamic theory of contact angles, where local slip is allowed to avoid a mathematical singularity at the contact line due to the abrupt change of material properties at the interphase. 2 In principle, the tolerance of local slip leads to a minimization of the wall shear gradient to a value of zero. With the linear transition up to substrate speed, diverging wall shear stress is avoided. Since the local slip tolerance at the boundary cell at contact line position dampens the velocity component in coating direction as well, a positive effect on the methodrelated air entrainment is expected.
The boundary condition that was developed in this work is based on a no-slip condition, which was extended by the localization of the contact line and the application of the linear interpolation of the wall velocity from slip velocity to substrate velocity. This initial version of the code will only work for flat surfaces, which need to be oriented along the origin planes. For example, the substrate wall normal vector The linear interpolation of the velocity field from a generic slip condition to the wall velocity at the substrate boundary layer is then applied in proximity to the dynamic wetting line, as shown in Fig. 6. On a basic level, the program flow can be described in four steps. A more detailed description of the flow is accessible in the online appendix material.
1. Pass through all the cells at the wall until the contact line is found 2. Move one cell further in coating direction 3. Calculate the distance from the current cell to the contact line 4. Apply the interpolation if the distance of the current cell from the interphase cell is smaller than the slip distance, otherwise, leave the loop 5. Go back to step 2

Contact angles
In a non-moving 3-phase-system, for example, a resting droplet on a solid surface surrounded by gas, the Young equation describes the ratio and directions of the surface energies and the resulting contact angle that is formed by the liquid surface with respect to the solid surface. If a force is applied or the surface shape is deflected from its equilibrium state, Young's equation loses validity due to an emerging (de-)wetting force. 25,26 In the case of slot die coating, this situation occurs at the moving contact line, where the liquid surface shape is continuously deflected by the moving substrate. Therefore, a boundary condition with a velocity-dependent dynamic contact angle is necessary. The standard OpenFOAM libraries offer the dynamicAlphaContactAngle boundary condition, which applies the following calculation:   The dynamic contact angle H dynamic is derived from the static contact angle H static , the advancing/receding contact angle H advancing / H receding, and the fit parameter u H as a function of the relative velocity u rel perpendicular to the wall normal vector. As mentioned in the experimental section, the dynamic contact angles were measured in dependence on the substrate speed and the static contact angle is obtained by the sessile-drop measurements. The parameters H advancing , H receding , u H, and u rel were fitted to the equation of the OpenFOAM boundary condition using a least-square-fit method.

3D-model
The three-dimensional mesh generation is performed by using CAD files and the mesh generation utility snappyHexMesh, which allows quick and straightforward changes in the slot die geometry and other process-related dimensions. The CAD model was created in OpenSCAD with comprehensive parametrization, which is exported in .stl-format and then imported into OpenFOAM. With the inclusion of the third dimension, the inner geometry of the slot die is adjustable by defining a chamfer at the outlet. In view of the minimization of the computational effort, a symmetry plane is placed in the center of the film, which leads to a halving of the cell number (Fig. 7).
The boundary conditions are based on the twodimensional model. Since the substrate area was divided into several sub-planes for mesh generation purposes, they are all assigned the same condition. Table 3 provides an overview of the used boundary conditions of the three-dimensional model for pressure (p), velocity (u) and phase fraction (a). It is noted here that the velocity boundary condition at the moving wall is chosen to be the conventional no-slip condition since the effect of the method-related air entrainment is generally lower in three dimensions and the usage of a local-slip boundary condition would require a presimulation in order to provide initial values to keep the simulation numerically stable.

Mesh generation
The simplest way to generate a mesh for CFD simulations is equidistant meshes, in which the smallest dimension, which is relevant for the simulation, determines the global resolution. In this work, the smallest dimension to be resolved is the edge elevation of the film with a magnitude of about 10 to 50 lm. When considering the dimensions of the slot die, for example with a width of 60 mm, the number of cells in equidistant meshes will be very large and may exceed the available computing power. A reduction of the number of cells while maintaining accuracy and mesh independence is possible by using local refinement methods, which can be applied statically and dynamically.
In case of static refinement, the snappyHexMesh utility provides a wide range of functionality. Refinement regions with generic shapes (boxes, circles, etc.) are placed inside the simulation domain with a specified refinement level. Equally, the system boundaries or other defined planes can be defined as refinement surfaces with refinement layers on top, which is for example suitable for walls. However, the static methods, especially in transient simulations, inevitably lead to an inflated number of cells, since the exact position of relevant areas, for example, the interphase position, is not accurately predictable (Fig. 8).
A more flexible method is the dynamic refinement, where the areas to be refined are specified by criteria. During simulation, the mesh is modified accordingly so that the specified areas have the assigned mesh resolution. This method is ideally suited for a multiphase simulation in which a given resolution is attached to an entire phase and/or the phase boundary. As a result, the total number of cells and the computational effort is reduced. In the OpenFOAM standard libraries, the dynamic refinement algorithm dynamicRefineFvMesh is available. One challenge of this standard method in combination with transient simulations is the uneven distribution of mesh cells to the processor cores. In unfortunate cases, this can lead to a processor receiving a multiple number of cells and  Fig. 7: Three-dimensional model with system boundaries and patch names the simulation time increases significantly. In order to restore an equal cell distribution, Rettenmaier et al. 27 developed a load balancing algorithm to extend the dynamic refinement toolbox. With this technique, a maximum imbalance parameter is set in the simulation controls. If this threshold is exceeded, a redistribution of the cells is performed. Another advantage of this algorithm is the ability to use the dynamic refinement in two-dimensional simulations, which is not possible with the standard libraries. To use this method in the current version of OpenFOAM, a port to v2012 by Scheufler 28 was used.

Dynamic wetting line
In this chapter, the results for the two-dimensional model are presented. First, the local-slip boundary condition, which was created in the scope of this study, is analyzed. The motivation for this development was the prevention of a method-related entrainment of gas phase into the liquid film and all resulting falsifications. As a performance test, it is suitable to investigate the phase fraction value in the immediate vicinity of the substrate wall, comparing the no-slip to the local-slip boundary conditions. In Fig. 9, the value of the phase fraction is visualized in dependence on the position in coating direction x.
On the upstream side, the phase fraction is zero in both cases, as expected for the gas phase. Below the upstream lip, the phase boundary occurs and its position is slightly shifted, due to the influence of the boundary conditions on the momentum equation. In case of the no-slip condition, the phase fraction increases up to a value of 0.95, followed by a settling to approximately 0.92. The profile also shows oscillations, which are areas of lower or higher phase fractions that move through the coating bead from the upstream to the downstream side. In contrast, the simulation with the local-slip boundary condition only has values of zero and one outside of the range of the phase boundary and the interphase itself is more sharply resolved by a factor of around 3. In this case, no oscillations are observed.

Contact angles
In order to select the combinations of boundary conditions concerning the static and dynamic contact angles, a systematic study of different combinations was carried out. For both velocity boundary conditions (no-slip and local-slip), the combination of a constant contact angle at the die surface and an empirical dynamic contact angle model at the substrate's contact  (2) substrate (1) wet film surface (2) edge area (2) The experimental data and the associated fitted correlations are presented in Fig. 10.

Mesh analysis
As a quantification method for discretization errors, the Grid Convergence Index (GCI) was used. 29 This method compares an integral functional value of simulations with different mesh resolutions to provide a measure for the discretization errors and their change for different mesh resolutions. For the two-dimensional model, the position of the upstream meniscus was selected as the integral value. Four grades of resolution were compared, coarse (5040 cells), medium (11580 cells), fine (25450 cells), and very fine mesh (38470 cells). The highest mesh resolution results in a GCI of 0.02 % (no-slip condition) and 2.8 % (local-slip condition), speaking for a sufficient grade of refinement.

Validation
The stability limits of the coating process from the simulations, analytical calculations, and experiments are compared. The investigation on the coating windows is performed in two dimensions since the defect mechanisms considered here are two-dimensional. In addition, the effort due to a three-dimensional simulation is not in relation to the effort involved. The simulations were performed by varying the dimensionless gap height G* in steps of 0.05 with a constant gap height of 180 lm. A coating is considered unstable if the coating fluid either swells in upstream-direction (upper limit), the meniscus shape fluctuates or air is entrained into the film (lower limit), as shown in Fig. 11.
Besides the simulative results from the no-slip and local-slip boundary conditions, the conventional onedimensional calculations are contrasted. For the experimental process limits, the coating speed was set to 1, 5, 10 and 20 m min -1 . Experience with this setup has shown that process limits can be obtained with an accuracy of about 10-15 lm. In Fig. 12, the coating windows are illustrated as the minimum (left) and maximum (right) wet-film height, that is necessary to ensure coating stability.
The deviation of the results for the upper process limit swelling is strongly noticeable. In general, the simulation with the local-slip boundary condition has the closest values to the experimental data, followed by the no-slip condition. The one-dimensional approach shows a significantly higher deviation, especially for low coating speeds. Moreover, the experimental determination of the swelling limit has an uncertainty since a marginally wetting of the upstream lip shoulder does not necessarily lead to an instable coating. A possible explanation for the underestimation of the upper limit is the neglect of the third dimension. In reality, with rising volume flow, the coating fluid has the ability to flow outwards to the lateral side of the die. Furthermore, the capillary pressure of the two lateral menisci is able to compensate for a rising bead pressure, before the upstream meniscus breaks and swelling occurs.
In regard to the lower process limits, the simulation with the local-slip boundary condition generates the lowest deviation of the minimal wet-film height from the experimental observations. Also, the value of the minimal wet-film height is generally higher than the results from the simulations with the no-slip boundary condition and the analytical calculations, resulting in an upwards shifted coating window.

Velocity profile in the coating gap
For the purpose of the mathematical prediction of coating windows, one-dimensional attempts are used to couple the pressure drops from the upstream to the downstream side. 3,4,7 In this case, velocity profiles are accessible by solving a transport equation iteratively with the assumption of a fully developed flow with Power-Law behavior and the neglect of inertial forces. To assess the results from the CFD simulations, the velocity profiles under the upstream (left) and downstream lip (right) are compared in Fig. 13 for process parameters u Web = 5 m min À1 , h Gap = 180 lm, and G* = 1.15. Both, the spatial coordinates and the velocity are presented as dimensionless values related to the gap height and the substrate speed, respectively.
In the case of the no-slip boundary condition, the velocity profiles match up very well with the simplified, one-dimensional approach. This is expected since the boundary conditions are identical. Besides, the difference in the viscosity models of the CFD simulations (Cross-Power-Law) and the Power-Law is not noticeable in the velocity profiles. In contrast, the local-slip boundary condition deviates, which may be associated with the influence of the boundary condition on the entire velocity field. In general, these results speak for a sufficient representation of the state of knowledge by the developed CFD model.

3D-model
For an estimation of the discretization errors, four different mesh resolutions were used to simulate the wet-film profile. Additionally, the static and dynamic methods of local refinement were put into proportion. In Fig. 14, the results for the coarse, medium, fine, and very fine equidistant meshes are shown on the left, while on the right the static and two levels of dynamic refinement are compared to the highest resolution of the equidistant meshes. This data was taken at a distance of 7 mm after the downstream die outlet. The coating width was set to 12 mm, as the equidistant mesh studies result in a large number of cells. In the case of the present dimensionless gap G* = 1.15 and the wet-film height of 156 lm, specified by the volume flow rate, the width of the coating edges only extend over 2-3 mm on each side. The plateau in the center of the film is independent of the coating width. Thus, the simulation is still relevant for productionscale coating widths. Due to the symmetrical structure of the simulation model, the half wet-film profile is shown.
In the case of equidistant meshes, it is obvious that the coarse and medium resolutions result in significant fluctuations of the wet-film surface. These fluctuations decrease for increasing resolutions so that the simulation with the very fine mesh results in a smooth surface. When applying local refinement methods, the fluctuations are generally lower and the wet-film profiles are in good agreement with the equidistant results using high-resolution meshes. For the quantification of the discretization errors, the GCI was calculated using the lateral position of the liquid phase as the integral function value, since the wet-film height is not suitable due to the high fluctuations at lower mesh resolutions. In the case of the highest resolution, the GCI amounts to 0.02 % and, thus, stands for sufficient accuracy. It is noted at this point, that the lateral liquid boundary line shows certain fluctuations that emerge within millimeters after the die outlet. The fluctuations are reduced in magnitude to < 10 lm with an increasing mesh resolution but remain existent for the time step ranges included in this study. In relation to the coating of 12-60 mm, the fluctuations amount to < 0.1 %.
In order to evaluate the local refinement methods, a more detailed consideration of the simulation effort is required. From an application perspective, important parameters are the computational effort in terms of the number of required processor cores and the simulation time, which is also influenced by the time step. The number of processor cores was set so that each core receives approximately 20000 cells. For each simulation condition, these parameters are comparably listed in Table 4.
As a result of the dynamic time step calculation, a higher mesh resolution leads to smaller time steps. This fact is attributable to the Courant number, which is defined as the product of velocity and time step divided by the length of a cell edge.
During simulation, the time step is adjusted, so that a maximum of Co = 0.9 is not exceeded in the entire mesh. The reduction of the time steps leads to an overall non-linear increase in simulation time, which becomes apparent from the fine mesh on to higher resolutions. With regard to the ease of use of this model, the mesh resolution should only be selected as high as necessary to ensure sufficient accuracy while still obtaining acceptable simulation times. As shown in the table below, the time steps decrease significantly for the higher mesh resolutions or higher levels of refinement. In the case of the static refinement, it has to be considered that the preparation of the mesh before actually starting the simulation involves several iterations to ensure the mesh orthogonality, quality, and the position of the refined mesh areas, which can be disadvantageous for practicability.
An insightful illustration is provided by the comparison of the simulations of the equidistant mesh medium and the dynamic refinement level 1. While the computational effort and the simulation times are similar, the quality of the results for the refined mesh is significantly higher. Moreover, the dynamic refinement level 2 only uses about half of the processor cores compared to the equidistant mesh while maintaining a decent quality of the results.

Conclusions and outlook
In this study, a two-phase CFD model of the slot dies coating process of lithium-ion battery electrodes was developed in OpenFOAM with the use of the volumeof-fluid method. The focus was set on the optimization of different mesh generation and refinement methods, the mapping of the wetting behavior of the battery slurry on the substrate, and the transfer to the threedimensional simulation. For validation, experiments were carried out to compare the empirically determined coating windows to the simulation results.
The wetting behavior was firstly analyzed in the twodimensional model, which includes the physical phenomena oriented in coating direction and along the slot axis. At the dynamic wetting line at the upstream meniscus, gas is entrained into the liquid film as a result of the combination of the used simulation methods, which leads to a falsification of the fluid properties in the coating bead and limits the general accuracy of the simulation. As a countermeasure, a boundary condition was developed, which allows local slip at the contact line to dampen the shear gradient and prevent the gas phase from entraining into the liquid film. The dependence of the dynamic contact angle on the coating speed was implemented using an empirical correlation based on results from coating trials with an extended camera setup. The functionality of the boundary condition was tested by analyzing the homogeneity of the phase fraction in the simulated liquid phase along the substrate wall boundary layer. As a result, the influence on the method-related air entrainment was shown to be eliminated. The comparison with experimental observations of upper and lower process limits showed that the simulation results were noticeably closer to the conventional, one-dimensional approaches for mathematical estimations of coating stability. In addition, the developed boundary condition was even closer to the experimental data than the simulations with the ordinary no-slip boundary condition. Subsequently, the findings of the twodimensional model were transferred to a simulation model with three dimensions. For this application, a fully parametrized CAD-based model was created in OpenSCAD to simplify modifications of the slot die dimensions. Due to the unbalanced spatial extent of the coated film, conventional mesh generation methods were found not to be sufficient in terms of computational effort and simulation times. Thus, different static and dynamic methods of local refinement techniques were compared and integrated into the model. For the latter, it was also necessary to implement an algorithm to ensure an equal distribution of the computation effort to the processor cells to minimize simulation times. The advantage of the dynamic refinement method in regard to the simulation time, result quality, and computational effort was shown. The developed simulation model enables in-depth reviews of the influence of process conditions or material properties on the process limits, flow characteristics and, in the case of the three-dimensional model, the wet-film profile. In addition, the influence of shaped slot die geometries can be investigated. Furthermore, this model can be used for a detailed analysis of the formation of edge elevations, since these are challenges in the coating industry and the underlying physical mechanism is not fully understood. As a long-term objective, this two-phase model can be transferred to three phases to simulate multilayer coating.

A1 Gitlab repository
Additional files and the annotated source code of the local-slip boundary condition is provided in the following GitLab repository.
https://git.scc.kit.edu/kw3888/cfd-model-of-slot-die-c oating Acknowledgments The authors would like to acknowledge financial support of the Federal Ministry of Education and Research (BMBF) via the ProZell cluster projects ''Sim4Pro'' (Grant No.: 03XP0242C) and ''HiStructures'' (Grant No.: 03XP0243C). The authors would like to thank all involved mechanics at the institute, Konrad Dubiel and Christian Wachsmann for particular contributions to this work.
Funding Open Access funding enabled and organized by Projekt DEAL.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/b y/4.0/.