Computational fluid dynamics simulation of the roll-to-roll coating process for the production of thin film composite membranes including validation

On the way toward a carbon-neutral economy, a rise in the demand for separation technologies can be expected. This holds especially for membranes, which are energy efficient and thus very promising technologies. However, this challenges membrane researchers to consider the sustainability and scalability of their membrane fabrication processes. At our institute, we employ a roll-to-roll coating process for the production of thin film composite membranes. This procedure is relatively easy to scale up and can be adapted for different polymers/solvents and thus applications. To increase the production efficiency and optimize the process for any polymer of interest, it is necessary to develop a solid understanding of the physics of this production system. Therefore, we would like to present a numerical model based on computational fluid dynamics that can predict the film thickness of a polydimethylsiloxane-based polymer coated in a roll-to-roll setup. In the future, this model should improve the production efficiency and fine-tuning of process parameters. We verify the numerical procedure with a mesh refinement study and validate the predicted film thicknesses with experimental results for different roll speeds and polymer concentrations. The predicted variability of the thickness is assessed by a design of experiments study and compares relatively well to the measured variations of the coated membrane thickness.


| INTRODUCTION
In the following study, we would like to present a numerical model for the prediction of coating film thickness of thin film composite (TFC) membranes for gas separation, as produced in a roll-to-roll coating device designed in our group at Helmholtz-Zentrum Geesthacht (HZG). This method allows for the large scale production of membranes with an effective membrane area of 100 m 2 or more. The production of TFC membranes in our group goes back many years, [5] with one prominent application being CO 2 selective membranes using PolyActive as the selective material. [4,6] More recently, membranes of this type have been used in pilot scale setups to remove CO 2 from industrial sources. [7,8] For a basic understanding of the roll-to-roll coating process the reader is referred to the cited literature and Section 3.1, where we discuss the membrane preparation setup in more detail.

| Motivation for membrane based gas separation processes
Since the first industrial revolution, the economic model of the world is a linear one. Materials are extracted from the earth, converted to useful products via manifold value chains across the globe, distributed and then, after consumption, finally disposed as waste. Thermodynamically, these processes create complex structures and thereby decrease entropy; hence, they need energy input and release waste heat. Until today, the prevalent energy sources are fossil fuels, but the extraction of materials and energy sources, as well as the disposal of waste, puts an unsustainable burden on our eco-systems, drives climate change and causes biodiversity loss. A new industrial revolution is needed, which limits the extraction of materials and replaces it with a circular economy that is fueled by non-fossil energy sources. This vision implies a great demand in efficient separation technologies of all sorts. Independent of implementation, they are a necessity to create circular material flows. For instance, we need separation technologies for water desalination, for biotechnological processes that seek to replace the current petrochemical industry or we need processes for gas separation to close the loop between exhaust gasses and power-to-gas or power-to-liquid technologies.
This can be seen in recent publications elucidating the way towards a net-zero CO 2 energy system fueling the economy. It seems inevitable, that many elements of a new post-carbon industry require gas separation and purification technologies at different points. [1] Carbon capture and storage or utilization might be the most prominent applications for CO 2 separation technologies, for which several methods are possible. These include processes based on special liquids or solids for ab-or adsorption, respectively, mineralization or physical processes like cryogenic cooling and membranes.
Because of the lack of toxic chemicals in the process, the lower energy input and less waste, membrane processes appear to be the most environmentally friendly overall. [2] As for now, polymeric TFC membranes have the advantages of much easier large scale production, flexibility of geometry with respect to installation in membrane modules, price and operational ruggedness compared to ceramic membranes which in turn have advantages in terms of temperature and pressure stability as well as constancy of permeation performance. [3] Moreover, further scale up of existing membrane systems is relatively straightforward due to the modularity.
Nevertheless, as big the promises of polymeric membranes might be, to be suitable for industrial production, it is necessary that the fabrication procedure is very easily scalable and the adoption to different separation tasks is feasible. To produce large areas of flat sheet membrane with very thin selective polymer films, we employ a rollto-roll coating process where dense polymer layers are applied on top of a porous support. Our in-house produced TFC membranes consist of a highly permeable and easy to apply gutter, a selective and a protective layer. [4] With this production setup, it is possible to tailor the membrane for a specific separation tasks by choosing the right polymer and thickness.
However, this freedom does not come without a challenge, as a change in polymer and solvent system also changes the physical properties of the polymer solution, especially viscosity and surface tension. This, in turn, changes the coating behavior from case to case. Currently, the fine-tuning of process parameters relies primarily on empirical knowledge and personal experience. To be ready for an efficient and agile production on industrial scale, a more structured approach guided by a model is desirable.

| Challenges of modelling meniscus coating flows
Some of the first mathematical analysis of roll-to-roll processing go back to the 70s and were done by Greener and Middleman. [9,10] Among other aspects, the authors were interested in the understanding of ribbing behavior leading to wave-like surface irregularities during roller coating. For the lack of computational power, they used a lubrication approximation, that is, ignored inertial terms in the Navier-Stokes equations. However, the lubrication approximation is only useful, when the flow field is quasi-unidirectional, that is, the horizontal velocity components parallel to the walls in the gap are much larger than the vertical ones. For the roller coating device this limits the approximation to the "flooded" regime, where the applicator roll delivers more liquid to the gap than what could go through. When the arriving film thickness is much smaller than the gap size, the regime is called "starved" and the flow structure in the bead, which contains a primary transfer jet and two sub-eddies becomes important, besides the pressure gradients caused by the surface tension. The amount of liquid transferred to the membrane essentially depends on the position of the stagnation point at the exit of the gap. For a schematic image of the starved flow structure, see Figure S5 in the Supporting Information. For experimental images of the flows in reverse roll coating the reader is referred to the publications of Colye, Gaskell and Savage, among others, who made further numerical and experimental contributions to this research field in the 90s. [11][12][13][14] With the rise of information technology and faster computational resources, the phenomenon was also tackled with CFD, often solved with finite element methods. [15][16][17] A large part of the work was motivated by problems that arise when applying this coating technique for coatings of paint, organic photovoltaic or semi-conductor materials.
Today, despite the highly increased computational power, the problem still is a difficult one to solve. The main challenges arise from the free surface flow, that is, the fact that the boundaries of the flow are not known a priori and have to be solved simultaneously with the flow field. In former studies, the computation of the free surface was done with a Volume of Fluid method, [18] an "interface capturing method" that tracks the surface with a scalar fraction function. Here instead, we will use an "interface tracking method" that tracks the interface position directly, as it coincides with the boundaries of a moving mesh.

| Novelty and outline of this study
While dip and slot die coating are commonly investigated membrane fabrication methods in terms of experiments and modeling, the application of a roll-to-roll coating process to produce large-scale TFC membranes is not tackled yet in the literature in terms of modeling and simulation. As has been shown in Section 1.2, simulations of roll-to-roll processes have been done, but it is not clear in how far they can be adapted in a membrane fabrication context, especially for the relevant polymers and the very thin films arising from the solutions of low viscosity applied here. Since the solvent is constantly evaporating in such a setup, the resulting concentration gradients of the polymer might alter the local density and viscosity in a way that the models discussed in Section 1.2 cannot describe, as they usually consider the liquid to be homogeneous. In this study, we investigate a computational fluid dynamics (CFD) model that accounts for these effects and verify it experimentally for a range of roller speeds and polymer concentrations. The presented study should also answer the question if a 2-dimensional simplification can nevertheless be enough to improve the understanding of the coating process and to enhance the optimization potential. We envisage that this model may also serve as a valuable tool in the upscaling of the production process.

| RESULTS
In the results Section we show measurements of the toplayer thickness of a TFC membrane produced in a roll-toroll coating device (see Section 3.2). In all cases, PDMS was used as the coating polymer, as this standard silicone is often utilized in polymeric membranes. We compare the measured thicknesses and variations thereof for several process parameters to a two-dimensional CFD model, which was numerically validated in a mesh refinement study. Moreover, we show experimental results that estimate the evaporation rate and show the temperature profiles during the coating procedure.

| Mesh refinement study
First, a mesh refinement study demonstrated that the simulation results will not change substantially anymore once the mesh reaches a certain resolution. The results can be seen in Figure 1. While the global output parameters total volume of liquid V liquid in the domain and total evaporation rate _ M e do not vary with mesh resolution, the local output parameters like polymer thickness δ (ie, after evaporation of the solvent) and fluid film thicknesses H 2 and H 3 deviate strongly at low resolution meshes. When approaching a mesh size of 360 000 triangular elements none of the parameters yielded a F I G U R E 1 Mesh refinement study of the computational fluid dynamics (CFD) model. The number of elements corresponds to the part of the domain containing the liquid. Output parameters H 2 and H 3 refer to the liquid film thickness on the membrane and applicator roll after the gap, respectively. The parameter _ M e denotes the integrated evaporation rate, V liquid the total volume of liquid in the domain and the parameter δ means the thickness of the polymer on the membrane after evaporation of solvent. All parameters are given as a relative value with respect to the value at the highest resolution substantially different result from the highest resolution. Therefore, all simulations presented here were done with meshing parameters that correspond to the mesh with 360 000 triangular elements.

| Evaporation and temperature measurements
By measuring the mass loss of polymer solution on a rotating applicator roll with no membrane roll present in the coating machine, we estimated the average evaporation rate and its deviation. The evaporation rate was determined for different applicator roll speeds and concentrations, but since there was no apparent correlation of evaporative mass loss and roll speed, the mean value of mass loss for all roll speeds is considered and amounts to 0.025 and 0.023 g s for PDMS concentrations of 1 and 3wt%, respectively. A bar chart with corresponding SD can be found in Figure S1 in the Support Information.
Since no general trend for evaporation rate dependence on concentration or roll speed could be concluded, further simulations were carried out with a mean evaporation rate of 0.0244 g s . With an estimated surface area for the evaporating film of 0.03 m 2 this yields average evaporation rates of approximately 0.81 g m 2 s . The SD based on this sample is 0.0037 g m 2 s . To assess the validity of the assumption of isothermal evaporation, we carried out temperature profile measurements for three different applicator roll speeds. The measurements were done during the coating process and, as can be seen in Figure S2 in the Supporting Information, the temperature on the applicator roll is approximately at room temperature and does not go down substantially before the coated membrane is leaving the roll gap. This is the same for all roll speeds.

| Model verification
For the verification of the numerical model we measured the mass loss during several coating experiments. By subtracting the mean evaporation rate of 0.0244 g s (see Section 3.2) we determined the mass transferred to the membrane roll. By using the initial weight fraction of polymer in the solution, we determined the mass of polymer per unit area on the membrane. Division by the PDMS density yields the estimated average polymer thickness that is compared to the respective simulation result in Figure 2. The empirical error bars are based on the SD of the evaporation rate creating an uncertainty in the evaporative correction to be subtracted from the total mass loss rate. The comparison shows a reasonably good agreement of the simulation results with the experimentally estimated values for all parameters except for the highest PDMS concentration, where the experimental values are consistently around 20% lower than the simulation predictions.
F I G U R E 2 Comparison of simulation results of the final PDMS layer thickness for different applicator roll speeds and polymer concentrations (solid lines) versus experimentally estimated thicknesses based on mass loss measurements during coating (symbols). Errors bars represent maximum and minimum due to uncertainties in the evaporation rate. The colors indicate the different concentrations, as denoted on the respective right side F I G U R E 3 Comparison of simulation results (white dots) of the final PDMS layer thickness for different applicator roll speeds and polymer concentrations versus experimentally estimated thicknesses based on cross-sectional SEM images (boxplots). Small black dots show data points that are outliers, that is, deviate from the median by more than the 1.5-fold interquartile range. Error bars on the white dots denote the estimated SD of 25%, respectively Additionally, the thickness for selected samples was determined from cross sectional scanning electron microscope (SEM) images. Figure 3 depicts median values as well as SD for the top layer thickness in the cross-sectional images. Representative cross-sectional images are given in Figure 4. The thickness of the high concentration solution (3 wt%) shows the biggest variation across different samples, while the spectrum of thicknesses for the 1 w% solution is most narrow. As can be seen, the mean values for both samples with 1 and 2 wt% PDMS are a bit smaller than the simulation results. The median values in the SEM measurements for 1 wt% are around 25-28 nm and around 80-90 nm for 2 wt%, while the simulations predict 40-50 nm, and 120-150 nm, respectively. For the higher concentration of 3 wt%, the median of 330 nm is close to the simulation results. One interesting observation is that the measured variation is not symmetric around the median but shows higher deviations above it, with maximum values being more than twice of the median in some cases. The simulation results have error bars of +/− 25%, which is the expected deviation based on the results in Section 2.3.2.

| Variability of membrane thickness
For the parameters, which are fluctuating or not exactly known, we estimated the maximum variability in PDMS thickness using a design of experiments (DoE) study combined with a response surface method. The study was done for 1 and 2 wt% polymer concentration and in both cases the evaporation rate had the lowest effect on the resulting thickness with a higher evaporation rate yielding a slightly increased thickness by few nanometers (corresponding to not more than 1-2% of the total thickness). The strongest impact comes from the variation of membrane roll velocity, with an increase in thickness that is almost directly proportional to the increase in velocity. Thus, an increase of membrane roll speed by 20% will also increase the coated thickness by roughly 20%. When the gap size is in the range of 135 to 195 μm, it changes the thickness by 5 to 10 nm and hence is also not very substantial. Taken together, the DoE study suggests a deviation of approximately +/− 25% in the chosen parameter range. Another result of the DoE study is the lack of any strong interaction effect between parameters. In the chosen parameter range, to give an example, the impact of membrane roll speed does not depend on the evaporation and gap size and vice versa. This is also expressed by the statistical information on the response surface results, which show a non-significant parameter fit for all interaction terms. This can also be seen by the almost straight contour lines of the response surface plots (see Figures S3 and S4 in the Supporting Information). However, the characteristics of the system change rapidly when a certain minimum gap size is reached, as can be seen in Figure 5. In that case, the coated thickness may suddenly increase by a factor of two or more. The threshold for this minimum gap size depends on the polymer concentration, for a concentration of 0.5 wt% the effect cannot be seen (data not shown).

| Membrane preparation
All TFC membranes were produced in-house at HZG, as already described elsewhere. [4] The fabrication is a multistep process. First, on top of a non-woven fabric made of polyester, a support layer of polyacrylonitrile (PAN) was casted by means of phase inversion. A more detailed description can be found in the literature. [19] The resulting PAN support layer has a thickness of around 50 μm, while the total membrane thickness (support and non-woven) is usually around 200 μm. Because of fabrication tolerances of the non-woven base material, a few micrometers deviation is possible.
The subsequent thin films consist of a polydimethylsiloxane (PDMS) based gutter-and protective layer, as well as the selective layer, for example, PolyActive, in between. However, for the verification of the model, we only consider the PDMS-based protective layer. The coating step is done with a roller-coating device, as depicted schematically in Figure 6 and further described in the cited literature. [4,6] The polymer solution in the basin is transferred to the membrane via an applicator roll. A meniscus is formed in the gap between roller and membrane and the forces in this meniscus region primarily determine the resulting transferred film thickness on the membrane. As in the gap region the rotation of the rolls is in the opposite directions, this is also called reverse roll-coating. During the roller coating process, the solvent will start to evaporate already but for the most part the remaining solvent is removed by drying the wetted membrane in an oven that immediately follows the coating step. As the solvent evaporates, the remaining polymer forms a dense film, where the polymer chains are also bound covalently by a catalyzed crosslinking reaction. Typically, the flat-sheet membrane rolls that are produced with this setup have a width of 0.3 m and a length of 200 m.

| Experimental verification
During membrane coating, the mass loss of the polymer feed solution is measured online as a function of time. A data acquisition and control software written in LabVIEW is employed for this task. Since the overflow is fed back to a reservoir on the same balance as the feed, the measured mass loss represents the sum of solvent evaporation in the coating apparatus and the transferred fluid.
To determine the evaporative flux, a series of measurements were conducted with only the applicator roll spinning in the solution basin while no membrane roll was present. This mass loss then accounts for the evaporation only and the calculated average evaporative mass loss from those measurements was subtracted from the total mass loss measurements during membrane coating to obtain the transferred mass flow. A scheme of the experimental setup is depicted in Figure 6.
In some cases, additional cross-sectional images were taken and the thickness of the PDMS layer was determined with a SEM (Merlin, Carl Zeiss Microscopy GmbH, Oberkochen, Germany). For the sample preparation, the membrane surface was coated with a thin layer of platinum (approximately 2 nm). This step was followed by fixing a self-adhesive copper tape on top to protect the surface. The cross sections of the membranes were prepared by argon ion milling in a PECS II (AMETEK GmbH, München, Germany). Finally, to eliminate charging effects during SEM investigations, the cross sections were coated with a 4 nm thick carbon layer. The section of the membrane sample was randomly chosen. For each sample, the top layer thickness was determined nine times in subsequent 100 μm long frames.
Additionally, the temperature distribution in the coating device was measured with an infrared camera (FLIR i60, InfrarotTec Systems, Ranstadt, Germany).
The process parameters of all conducted experiments are summarized in the Supporting Information (Table S1).
F I G U R E 6 Scheme of the mass loss measurement setup. Both feed and waste reservoir are placed on a balance and the mass measurement is transferred digitally to a computer with a custom data acquisition program saving the values in certain time intervals (≈20 s). A detailed description of the production process can be found in the literature. [4,6] 3.3 | Numerical analysis

| Simplifying the problem to two dimensions
In reality, the basin containing the polymer solution has a central feed inlet and a free overflow at the sides. Since we measured no big change in viscosity and density when comparing the stock solution and the liquid leaving the basin at the overflow, we expect no substantial concentration gradient along the axis of the applicator roll. Hence, we simplified to the problem to a 2-dimensional model including the overflow as shown in Figure 7.
The applicator roll speed is v 1 while the membrane roll speed is v 2 . Thick lines in the geometry represent walls of the rolls and basin, while the thin lines are the open boundaries of the system. The dashed line represents the liquid-gas interface.
The incoming feed flow is _ M in , the evaporation mass flow is _ M e , the mass transferred to the membrane roll is _ M m and the flow of the remaining mass on the applicator roll is _ M a . Subtracting the evaporation, transferred and remaining mass flows from the feed flow yields the outflow _ M out , which was estimated to be 0.1 g s . The exact value is not known but based on visual observation it can be guessed to be somewhere between 0.05 and 0.5 g s . It should be noted that a more precise estimation of the overflow rate is unfeasible in the current setup, as it depends on the volumetric flow of the feed pump, which is adjusted manually and will also fluctuate slightly between experiments and during the coating process. However, preliminary simulations yielded no substantial difference in the results when the estimated minimum or maximum overflow rate was used, thus an overflow of 0.1 g s was arbitrarily chosen for all further studies. Note also, that the inflow was not set to a fixed value but calculated during the simulation. It is a variable that, at each time step, equals the sum of the fixed overflow, the evaporation _ M e , and mass flow of liquid leaving the domain on the applicator roll ( _ M a ) and on the membrane ( _ M m ). Effectively, this means that the inflow is constantly adjusted during the computation to ensure that the total amount of liquid in the domain stays constant.
For clarification, the initial film height H 1 in Figure 7 refers to the thickness of the film when it leaves the dynamic meniscus region. The denoted angle α of 15.8 is used later on for the estimation of the liquid surface area in Section 2.2. The value for α is known as it follows from the geometry of the roll and basin, which were constructed in-house. Additionally, the scaling of the model is represented by the arrow in the lower left corner.

| CFD model
For computing the fluid dynamics during the reverse-roll coating, we employ the finite element method that is implemented in COMSOL Multiphysics. [20] Besides the F I G U R E 7 Scheme of the roll coating device and 2-dimensional model of the meniscus coating between the applicator and membrane roll operating in reverse mode fluid mechanics, the system features transport of chemical species, with the local concentration affecting the properties like density and viscosity. Moreover, the gas and liquid phases need to be considered in a free surface flow, in which the location of the free interface is not known a priori but must be solved simultaneously. Here, the interface of the two phases is modeled by an Arbitrary-Lagrange-Eulerian moving mesh method. [21] This is a method belonging to the class of interface tracking procedures, where the interphase coincides with the corresponding mesh boundary.
As for any model, we need to make some simplifying assumptions. Firstly, the system is assumed to be isothermal. The enthalpy loss due to evaporation as well as heat generated by dissipation is not considered in this first version of the model. In consequence, a constant evaporation flux is assumed. Secondly, we do not consider a third spatial dimension in the model. Thus, gradients along the direction of the roller axis are not taken into account.
Thirdly, the local polymer concentration has no effect on the surface tension and the concentration is low enough to assume Newtonian viscosity. Note that this, in combination with isothermal conditions, means that we do not consider any secondary flows like the Marangoni effect.
Lastly, we assume a perfectly smooth and homogenous surface, on which the thin films are coated. The porous structure of the support material will have an effect primarily on the gutter layer, therefore we focused on the top layer for now.
The model contains two phases, the gas and the liquid phase. However, for this model the gas phase is only important in terms of the atmospheric pressure working against the surface tension. The flow of gas is very slow and irrelevant for now, thus the following description will focus on the liquid phase only. The description of fluid motion is given by the well-known Navier-Stokes equations: where u ! is the velocity vector, ρ is the density and the divergence of their product forms the mass balance (Equation (1)). The source term S e equals zero everywhere, except for the free surface where it equals the evaporation rate. For the evaporation rate, we used experimental data from evaporative mass loss measurements (see Section 2.2). However, for the model we need to compute the evaporation per surface area. Therefore, we divided the evaporation rate by the estimated wetted surface of the applicator roll. The wetted surface is approximated by the following formula where the angle α (see Figure 7) equals 15.8 . For the estimation we use the applicator roll radius r a (0.0271 m) and the roll width r w (0.3 m).
180 + 2*α 360 Á π Á 2 Á r a Á r w = 0:03m 2 In the momentum balance (2), P denotes gauge pressure, η denotes viscosity, g ! the gravitational acceleration and F ! st represents the surface tension in vector notation, which is also only present at the interface and is given by the following: In Equation (3), n ! is the interface normal vector and r s is the surface gradient, that is, the conventional gradient with the absence of the component normal to the interface. The surface tension σ is assumed to be constant and does not depend on the polymer concentration, so the last term becomes zero. As discussed earlier, the viscosity will change with polymer concentration, but the shear rate is assumed to be low enough to be in the Newtonian regime. Analogous to former finite element analyses [17] the free surface boundary condition is: where the index 1 refers to the evaporating liquid and the index 2 to the gas phase and A i denotes the interfacial area. Note that this implies equal velocities of both phases in the tangential direction of the surface. It also means that the normal velocity in both phases equals the interface velocity. As for the boundary condition at the outlets, the normal stress τ vanishes that is, .
On all walls of the system a no-slip condition is employed, hence u ! = u ! wall . The contact angle between liquid, gas and wall in the meniscus region is a complicated issue, since here the surface forms a dynamic wetting line, the description of which is still under scientific debate. [22] The implementation of a conventional no-slip condition would cause a stress singularity, as two contradicting boundary conditions arise at the dynamic contact point. The conventional no-slip condition of fluid mechanics would assign the velocity of the wall directly adjacent to the fluid element while reaching a steady state position at the contact point requires the local velocity to be zero. The implementation in COMSOL takes the approach of Gerbeau et al. [23] by applying a force to the fluid which is tangential to the wall surface and moves the fluid towards the equilibrium contact angle (here it is assumed to be 90 ).
In addition to the fundamental fluid mechanical equations, the dissolved PDMS is modelled as a diluted chemical species with concentration c and its distribution follows from the convection-diffusion equation For a short-chained, pre-crosslinked PDMS, an estimated polymer radius of 3.3 nm was used based on information of the PDMS supplier that cannot be disclosed for confidentiality reasons. The diffusion coefficient D is approximated by the Einstein formula: where k B is the Boltzmann constant, T is temperature (293.15 K) and r polymer is the polymer 'radius', that is half of the polymer chain length. The source term S vanishes in the bulk fluid but it is positive at the free surface, with V m being the local mesh element volume: The polymer concentration influences the local density and viscosity as described by the empirical Equations (9) and (10). Both equations are quadratic functions, which were fitted to experimental values of density and viscosity of the silicone solution with 0.5, 1, 2 and 3 wt% using a quadratic function in which w P refers to the weight percentage also given in percent.
For computing the molar polymer concentration, we use the values found for PDMS and assumed a molecular weight MW PDMS of 70000 g mol . [23] and the respective PDMS density ρ PDMS of 0.97 g/cm 3 at room temperature. The silicone is dissolved in isooctane with a density ρ Solvent of 0.692 g cm3 at room temperature. For a 3 wt% solution this yields approximately: Therefore, we obtain the following relation between weight percentage and molar concentration: Finally, the resulting polymer thickness δ after evaporation has to be deduced. For this, we convert the weight percentage (which depends on the local concentration of polymer) to volumetric percentage by multiplying with ρ Solvent ρ PDMS . Then, we integrate the volumetric percentage over the film thickness H 2 at the membrane roll outlet: So, to paraphrase, we compute the percentage of the liquid film thickness at the outlet that is made up of polymer molecules only and consider this to be the final thickness after evaporation.

| Moving mesh
To solve the problem of free surface flows, COMSOL provides several algorithms. As we do not need to model topological changes but need to consider mass transfer, the most convenient algorithm is the Arbitrary Lagrange-Eulerian (ALE) method. Here we will only outline the core ideas of the method, for a detailed explanation the interested reader is referred to the literature. [21] Usually in fluid dynamics, the governing equations are formulated with respect to a fixed point in space, that is in the spatial reference frame. This viewpoint is called the Eulerian Frame and is used in cases, when it is not useful nor feasible to follow the movement of individual material particles. If, on the other hand, the deformation of a solid body is simulated, it is more useful to formulate the PDE's with respect to the respective, moving material elements. This is called the Lagrange viewpoint. The latter method is not feasible for a vast number of fluid particles, but has the advantage of being able to compute the movement of boundaries. The ALE method represents a mixture of both approaches, which is used for free surface flows, where the fluid movement is computed in a Eulerian framework, but the surface boundary moves according to a Lagrangian viewpoint in the material framework. The mesh is then deformed to follow the deformed boundary. To prevent strong deformation of the adjacent mesh elements, while not altering the rest, the mesh movement is distributed more or less equally across several mesh nodes. Thus, the movement of the individual mesh elements is somewhat arbitrary, hence the name Arbitrary Lagrange-Eulerian method.
COMSOL provides four different algorithms to govern the mesh deformation across the grid. In this study we only use the "Yeoh" smoothing method, which is based on a hyperelastic material model. This algorithm requires an additional parameter called "stiffening factor" that we have set to 10.
If the initial guess for the interface position is far away from the equilibrium value, the mesh deformation will cause numerical problems at some point. To circumvent this, COMSOL allows automatic re-meshing, which stops the solver if the quality any mesh element (ratio of shortest to longest edge) goes beneath a threshold of 0.2. The geometry is then re-meshed and the solution from the last time-step is interpolated to the new mesh. For a visual impression of the mesh used in the presented simulations, see Figure 8.
Even though we are only interested in the steady-state solution for a given parameter set, the moving mesh feature only allows transient simulations. Consequently, all simulations were conducted as transient studies until the relevant output variables would not change anymore. That way, it is possible to start with initial starting conditions far away from steady state and then slowly converge to the steady state solution. For an example of how the simulations are typically converging see Figure S6 in the Supporting Information.

| Parameter study
The parameter study is divided into three parts. The first part determines the mesh resolution that is needed for the simulation results to be numerically reliable. The second part considers the polymer concentration and applicator roll speed as the main process parameters and compares the simulation results to the estimated mean thickness based on mass loss measurements and crosssectional images taken by a SEM. The third part aims at estimating the variability of the thickness across the membrane for constant concentration and applicator roll speed. This variability arises from the fluctuation of other influencing parameters. Here, we consider gap size, evaporation rate and membrane roll speed to be most important. We use the established CFD model and vary these three stated parameters, while keeping concentration and applicator speed constant. With the use of a statistical tool called DoE, we thereby predict the thickness variation for a given setup, which can then be compared to the experimentally observed variations in cross-sectional SEM images.

| Mesh independence
As a first step, the numerical validity was checked in a mesh-independence study. Here, the resolution parameters in the liquid domain were increased in several steps, yielding a series of meshes with a total of approximately 74 000, 124 000, 173 000, 366 000, and 710 000 triangular elements. For the mesh study, we only looked at the smallest gap size analyzed in our work, as this should yield highest velocity gradients, thus requiring the smallest element size. Hence, we chose an applicator roll speed of 0.045 m/s and a gap size of 65 μm. The polymer concentration was set to 1 wt%. For a visual representation of the mesh, see Figure 8. We then ensured that this resolution is also sufficient for a simulation with polymer concentration of 0.5 wt%. Additionally, a simulation with 3 wt% for a gap size of 100 μm was done, using the same mesh resolution parameters (the smallest gap size is not feasible here, as it causes complete flooding in this case) and again the result did not substantially change when choosing a higher mesh resolution.
Moreover, we confirmed that two different solver types available in COMSOL (called MUMPS and PARDISO) as well as decreasing the tolerance of residual error did not change the results substantially for the stated concentrations and gap size (data not shown).
Building on this mesh study, the best compromise of mesh resolution parameters was used for all subsequent simulations (see Section 2.1).

| Model verification
To check the predictive quality of the model, simulations and coating experiments were done for a combination of polymer concentrations and applicator roll velocities. For concentrations of 0.5, 1, 2 and 3 wt% the coating was carried out with applicator roll velocities of 0.027, 0.036 and 0.045 m s , respectively. For a conclusive overview of the conducted experiments, see Table S1 in the Supporting Information.
In the corresponding simulations, the average gap size during coating was determined to be 165 μm by measuring the average total gap on several points of the applicator roll (approximately 375 μm with a SD of 25 μm) and subtracting the average thickness of the membrane support layer (approximately 200 μm with a SD of 5 μm). The overflow was 0.1 g s and the membrane roll speed v 2 = 0.011 m s . The evaporation rate of 0.8 g m 2 s was also constant.

| Estimating the variability of the predicted thickness
Moreover, a statistical analysis using DoE was done to deduce a typical variation of polymer thickness given a certain fluctuation or uncertainty of three process parameters, namely: gap size, evaporation rate and membrane roll speed.
For all statistical analysis we used the open source software R. [24] A Central Composite Design (CCD) plan with three parameters was employed and the results for predicted thicknesses were fitted to a full quadratic response surface model. For an introduction to the methodology see the R documentation online. [25] The design plan was set up with intervals ranging from 135 to 195 μm for the gap size, which was based on the expected SD of 30 μm for the gap size, the membrane roll speed varied between 0.0083 and 0.0133 m s , based on the typical fluctuations of membrane roll speed measured during the coating procedure. For the evaporation, the flux can only be estimated but based on the experiments presented in Section 2.2., it will be somewhere between 0.5 and 1 g m 2 s , which thus became the interval for the evaporation parameter. Table S2 in the Supporting Information shows the full CCD plan. For this study, the polymer concentration was set to be 1 wt% and the applicator roll speed was v 1 =0.045 m s . The complete study was repeated for a polymer concentration of 2 wt% to see if there is an increase in the relative deviations of polymer thickness.
Additionally, Gaskell et al. [14] also reported the occurrence of a secondary transfer jet in the meniscus region as an effect of very small gap sizes. Therefore, the effect of very small gaps going down to a minimum of 65 μm was analyzed additionally.

| DISCUSSION
The goal of this study was to set up a starting point for a conclusive mathematical model that can predict the coating thickness of polymers like PDMS in the roll-to-roll fabrication of TFC membranes. To demonstrate its value, the model was validated and then verified. Here, "validation" refers to the test of the numerical procedure, that is to say: showing that the implemented equations are solved reasonably well. By "verification" we mean the comparison with experiments, that is, checking if the implemented equations are sufficient to yield a good approximation of the physical system.
For the validation part, several mesh resolutions were compared and it could be shown that the polymer thickness, that is, the main parameter of interest, as well as other characteristic parameters, do not change significantly, when a certain resolution of the mesh is reached. Very coarse mesh resolutions cause strong deviations in the film thickness at the outlets of the domain. This seems to stem from problems of the solver to fulfill the boundary conditions for the free surfaces at the outlet, which are more pronounced when the resolution is too small. Further numerical validation might also include robustness tests of other parameters in the solver configuration, for example, the residual error threshold.
To verify the model, the coating thickness was estimated by means of cross-sectional SEM images and mass loss measurements during the coating procedure. The SEM measurements are on the same order of magnitude and indicate that our simplified initial model is already a relatively good approximation of the real system. The variation of thickness appearing on the SEM images goes up to 100% and more, which can only be explained in parts by the statistical DoE analysis. However, if we consider the effect of a more reduced gap size, the strong upwards deviations might well be explained. When the gap becomes so small that the incoming liquid film cannot go through the gap completely, the remaining liquid is essentially directly transferred to the membrane, thus increasing the coated thickness drastically (see Figure S5 in the Supporting Information for a schematic explanation of the phenomenon). Such small gap sizes are not implausible, since 30 μm variation of the radius of the applicator roll was only the SD, bigger deviations are possible.
Additionally, further instabilities in the process might realistically come into play. The most probable is a variation of the applicator roll speed. Further examination of this source of error should be tackled by future studies. Further irregularities in the coating process might be due to inhomogeneous polymer concentration, dust grains or small air bubbles in the gap region or increasing crosslinking effects of the polymer over time. It should be noted that the liquid film on the membrane roll, after leaving the roll gap, soon enters the drying oven, where the remaining solvent evaporates quickly and the PDMS undergoes cross-linking. Factors like oven temperature might also influence the final polymer thickness, but this question was not addressed so far.
As a secondary verification, the results in Figure 2 show the comparison of the mean polymer thickness from simulations with the polymer thickness estimated from mass loss measurements. In contrast to the SEM images, which show the thickness in a small cross -sectional area on arbitrary sample spots of the membrane, the measured mass loss yields an integral measure for the total amount of liquid (and thus polymer) that is transferred to the membrane. As the mass loss rate is a linear regression fit over several minutes of coating, it can be considered as an average value. However, these measurements have a relatively high uncertainty as the evaporation was fluctuating in experiments and cannot be measured independently during the coating procedure. The fluctuation of the evaporation rate may be caused by changes of air convection, temperature or humidity during the fabrication process. More probable though, most of the fluctuation might be explained by inaccuracies of the measurement method, for example due to air bubbles in the tubing that connect the overflow back to the reservoir and hinder the outflow. This might cause larger volumes of the liquid to stay in the overflow basin before going back into the reservoir, causing fluctuations in the outflow rate and thus an uncertainty in the deduced regression slope.
Nevertheless, the estimated transferred mass and calculated remaining polymer are generally in good agreement with the SEM images and the CFD model. In Figure 2, the deduced thickness for the 3 wt% solution is about 20% smaller than expected given the simulation and SEM results. One possible explanation could be, that within the chosen time interval, that is, coating area, the average gap size might have been increased, which would have strong enough effect on the thickness (see Figure 5). Since SEM images only represent a very small area corresponding to a coating interval of 10 −3 m, while mass loss measurements were averaged over several meters, this possibility cannot be ruled out and so these deviating observations are not necessarily contradicting.
In summary, the presented model delivers a goodenough prediction of the coated thickness to be used as a guide for the experimentalist. As one interesting phenomenon, the model shows an increase of coated thickness when the gap size is reduced beyond a certain threshold. Notably, a decreased applicator roll speed can also increase the thickness slightly. From a fluid mechanical standpoint, this can be understood, as the decreased velocity changes the force balance in the gap region. When the applicator roll speed decreases, the relative membrane roll speed, which is dragging the liquid to the opposite direction, gets bigger. In consequence, a larger part of the liquid is transferred to the top. The accuracy of the experiments presented in this study does not allow this effect to be validated. In principle, both phenomena described above can also be seen in former studies on reverse roll flows, as presented by Gaskell et al. [14] Finally, we would like to discuss the applicability of the presented model for interested researchers in this field. The coating process in the roll-to-roll device is a relatively complex interplay of several factors that need to be accounted for. Primary factors are the geometry and constructive details of the apparatus, the physical properties of the dissolved polymer and its solvent and the roller speeds, all of which interact and must be usually balanced out by the experimentalist. All experiments presented here were done solely with PDMS dissolved in Isooctane where the typical viscosity is in the range of one to a few mPa Á s. Changing the solvent or polymer (or both) can drastically change the viscosity. To sustain a similar liquid film on the roll (which is necessary for a stable meniscus), a much lower viscosity would require higher applicator roll speeds and vice versa. While for high velocities one might encounter a shear-rate dependent effect on viscosity (depending on the polymer), lower velocities increase the residence time of the liquid on the roll and in the meniscus gap, thereby increasing the impact of evaporation. In turn, this could mean that the process can no longer be considered to be isothermal, as temperature gradients of a few degrees can also have an effect on the local viscosity and also on surface tension. The same caveats may apply for solvents with a much lower saturation pressure than Isooctane or with higher enthalpy of vaporization, especially when paired with lower heat conductivity. Another factor to be aware of in this context is the possibility of a gel-like layer of polymer at the phase interphase during evaporation, which is usually referred to as skin-formation. In this study, our model suggests a concentration increase of approximately 10% at the interphase, for which we do not anticipate skin formation to occur. However, if the initial polymer concentration or evaporation rate is higher, this might cause skin-formation to occur at some point, which will impact the fluid mechanics, especially in the meniscus region. Additionally, the mentioned effects also depend on the material and geometry of the rolls and the gap size, because this will also affect heat conduction and residence time of the liquid.
All of the above points have to be kept in mind when applying this model to any other membrane coating setup where some of these parameters may deviate strongly. Any other coating setup will likely require some amount of adjustment to account for the stated effects. Moreover, from other experiments not reported in this study, we know that chemically more complicated block-copolymers cause inaccurate model predictions for reasons still unclear. In such cases, the presented model should thus be regarded as a starting point that provides a reasonable initial guess for more sophisticated models. We estimate that the model in its current formulation is only applicable when dealing with homogeneous polymers with similar concentration and viscosity in the solvent system. When the geometry is of comparable size, velocities are in the order of cm/s and the evaporation rate of the solvent is similar, we would guess that isothermal assumptions can be made. In those cases, the presented model might also be valid.

| CONCLUSION
In this study, we presented a multiphase CFD model implemented in COMSOL to describe the fluid flows in the roll-to-roll coating apparatus that is used in our institute for the large-scale production of TFC membranes. The comparison with experiments shows, that this relatively simple fluid mechanical model can still successfully predict the coating thickness for a dense PDMS layer with a reasonable accuracy. The next research steps should include the test of the same model for higher membrane roll speeds or other process parameters and different chemical systems, especially those solvent and polymer systems used for the coating of selective layers.
As a concluding remark, it should be noted that the presented model might appear more complex than necessary given the quality of predictions. The evaporation and concentration gradients do not play a big role in the given system and could be discarded without losing much predictive quality. We nevertheless decided to present the model in this way, as we want to build upon it when tackling more complex polymers (eg, block copolymers like PolyActive). As mentioned above, such polymers might show non-Newtonian viscosity or surface-tension gradient effects, that depend on the local flow field as well as on the concentration and other factors.