Dynamics of thin liquid films on

Cooperation with: A. Münch (Humboldt-Universität zu Berlin, Heisenberg Fellow at WIAS), T. Rieckmann (Fachhochschule (FH) Köln, Fakultät für Anlagen, Energie- und Maschinensysteme)

Supported by: DFG: ``Mathematical modeling and analysis of spreading polymer films'', DFG Research Center MATHEON, project C10

In this new project we develop mathematical models and algorithms for the numerical simulation of the dynamics of thin liquid films on rotating discs, typically used in so-called PET (polyethylenterephthalat) reactors.

Since PET is a commodity product with a wide range of applications, ranging from the production of technical yarns to plastic bottles, it is of great practical interest to optimize reactor design and hence PET synthesis, through the development of mathematical models that are able to predict and control the shape and thickness of the thin film.

PET is produced in polycondensation reactors that typically consist of a horizontal cylinder which is partially filled with polymer melt and contains disks rotating about the horizontal axis of the cylinder, thus picking up and spreading the melt in form of a thin film over a large area of the disks (see Figure 1). At the same time, the melt, which is fed into the cylinder from one end, slowly flows to the other end where it is removed. Therefore, the film profiles and thickness will vary from disk to disk, since the viscosity increases with the degree of polymerization.

For current reactor designs and operating conditions, it is typically assumed that the liquid is Newtonian. While the thickness of the film influences the reaction and output rate, they in turn do not significantly affect the film flow and thickness. Hence one can assume the chemical kinetics to decouple from the fluid dynamical part. Potential mechanisms for a coupling of the chemical kinetics to the fluid dynamics are Marangoni effects due to chemical or temperature gradients, which can be induced by the thickness-dependent evaporation of by-products, or by local release/consumption of heat by exo/endothermic reactions. However, the time scales for these effects to influence the film flow are much longer than the time required for one revolution of the disk (for typical rotation rates).

**Fig. 1:** Schematic of PET reactor (left). Front and side view of rotating disk (right).
$\makeatletter \@ZweiProjektbilderNocap[h]{0.45\textwidth}{reactor-crop.eps}{PETreac.eps} \makeatother$

The primary problem is therefore essentially fluid-mechanical and leads to a mathematical model for the rotating-disk problem involving a coupled system of a three-dimensional mass and momentum balance (the Stokes or Navier-Stokes equations for incompressible fluids) with a free surface, where surface tension (i.e. curvature) is relevant. It is related to the classical problem of wetting a plate by pulling it out of a liquid reservoir, [1, 2]. This problem comprises two scales w.r.t. the direction normal to the plate. The small scale is the height of the thin film, merging into the reservoir via a meniscus, where the height becomes large. For such problems, surface tension plays a decisive role on how the meniscus fixes the height of the thin film. A direct three-dimensional numerical simulation of such a problem which takes into account the disparate scales is computationally very intensive. We could, however, achieve a model reduction by exploiting the two scaling regimes, making use of ideas by [3] and [4], who developed extended lubrication models from the full equations for Marangoni-driven thin film flow. We developed these methods further for the geometry of a flat, vertical, and partially submerged rotating disk and by including centrifugal forces into the model. As a result, the PDE for the film thickness, written in cylindrical coordinates, has the form:

$\displaystyle {\frac{{\partial h}}{{\partial t}}}$ = $\displaystyle {\frac{{1}}{{r}}}$ $\displaystyle {\frac{{\partial }}{{\partial r}}}$ $\displaystyle \left(\vphantom{r \frac{h^3}{3} (p_r+B sin \theta) }\right.$ r $\displaystyle {\frac{{h^3}}{{3}}}$ (p_r + Bsin $\displaystyle \theta$ ) $\displaystyle \left.\vphantom{r \frac{h^3}{3} (p_r+B sin \theta) }\right)$ + $\displaystyle {\frac{{1}}{{r}}}$ $\displaystyle {\frac{{\partial }}{{\partial \theta}}}$ $\displaystyle \left(\vphantom{ \frac{h^3}{3} (\frac{1}{r}p_{\theta} + B cos \theta) + r\Omega h }\right.$ $\displaystyle {\frac{{h^3}}{{3}}}$ ( $\displaystyle {\frac{{1}}{{r}}}$ p + Bcos $\displaystyle \theta$ ) + r $\displaystyle \Omega$ h $\displaystyle \left.\vphantom{ \frac{h^3}{3} (\frac{1}{r}p_{\theta} + B cos \theta) + r\Omega h }\right)$ ,

p = - 2 $\displaystyle \left[\vphantom{\frac{\partial}{\partial r}\left(\frac{h_r}{N}\ri... ...}{r}\frac{\partial}{\partial \theta}\left(\frac{h_{\theta}}{rN}\right) }\right.$ $\displaystyle {\frac{{\partial }}{{\partial r}}}$ $\displaystyle \left(\vphantom{\frac{h_r}{N}}\right.$ $\displaystyle {\frac{{h_r}}{{N}}}$ $\displaystyle \left.\vphantom{\frac{h_r}{N}}\right)$ + $\displaystyle {\frac{{1}}{{r}}}$ $\displaystyle {\frac{{\partial }}{{\partial \theta}}}$ $\displaystyle \left(\vphantom{\frac{h_{\theta}}{rN}}\right.$ $\displaystyle {\frac{{h_{\theta}}}{{rN}}}$ $\displaystyle \left.\vphantom{\frac{h_{\theta}}{rN}}\right)$ $\displaystyle \left.\vphantom{\frac{\partial}{\partial r}\left(\frac{h_r}{N}\ri... ...}{r}\frac{\partial}{\partial \theta}\left(\frac{h_{\theta}}{rN}\right) }\right]$ withN = $\displaystyle \left(\vphantom{1+h_r^2+\frac{h_{\theta}^2}{r^2}}\right.$ 1 + h_r² + $\displaystyle {\frac{{h_{\theta}^2}}{{r^2}}}$ $\displaystyle \left.\vphantom{1+h_r^2+\frac{h_{\theta}^2}{r^2}}\right)^{{1/2}}_{}$ .

p_r(r, $\displaystyle \theta$ , t) = 0forr, $\displaystyle \theta$ $\displaystyle \to$ $\displaystyle \Gamma_{{out}}^{}$ andr, $\displaystyle \theta$ $\displaystyle \to$ $\displaystyle \Gamma_{{pool}}^{}$ .

For a stable and efficient numerical solution of this boundary value problem for this nonlinear fourth-order PDE, we developed modifications of the Keller-box method. The algorithm is based on FEM using linear, quadratic, cubic, and mixed elements together with either direct or a numerical integration for computing the coefficients of the stiffness and/or mass matrices.

In order to test and calibrate the implementation of our numerical scheme, we first implemented our new scheme for cases of Marangoni- and gravity-driven thin film flows. For these geometrically simpler problems, our group already has accurate ADI/finite difference schemes which we used as a reference. We have also implemented the above problem for the simulation of the thin film dynamics on a rotating disk and for comparison to experimental results, performed by the group of Prof. Thomas Riekmann at the FH Cologne.

Dynamics of thin liquid films on rotating disks