Numerical solutions to two-dimensional fourth order parabolic thin film equations using the Parabolic Monge-Ampere method

1.   Introduction

Thin-film flows caused by outside forces, like surface tension or gravity gradients, are used in business ^[1,2], geophysics, the environment ^[3,4], and biology ^[4,5]. Finite-time singularities that cause instability in fingers and film breakage, as well as wave propagation and regression, are some of the fascinating dynamics of thin-film liquid flows ^[6]. When both physical and mathematical factors are taken into account simultaneously, a certain class of thin-film flow issues, such as moving waves with fast oscillations and the effects of surface tension, become very important ^[7]. In these situations, the apparent changes usually happen in areas with very small scales, like surface tension, especially near the connection lines, where the free surface of the film has significant internal changes in how it curves. The effect diminishes and becomes more insignificant the farther the interior layers are from the surface tension. It is essential to keep in mind that these effects and how they change show up as shaky spatial transversal fingers in the inner layers ^[8,9,10,11]. Thus, it is crucial to research and appreciate the root of finger instability. Because the answer to these problems changes over time, it is best to use a moving mesh to find the numerical answer. Because of this, it has more precise grids in some areas and more dense grids in others, creating many variations in space. Often, an even grid with equal steps is employed. But this kind of mesh is not good enough because the size of the diffusion flows gets much bigger and it takes a long time for the fingers of instability to show up. Using the Parabolic-Monge-Ampere (PMA) method, the fourth-order parabolic thin film equations are numerically solved in two dimensions with constant flux boundary conditions. Instead of a mesh solution with the same step size, these inner layers will be accurately resolved by computers. An extended lubrication approximation or wavelength ^[6] is usually used to get the given equations. Surface tension changes the boundary conditions and controlling variables of thin-film liquid flows into a non-linear PDE that describes how the free surface of the film changes ^[12]. For example, a parabolic PDE, usually second-order, could be used to describe the surfactant concentration. For thin film liquid flow equations on a uniform or non-uniform grids, the finite difference method in two-dimensional PDEs has been applied in various numerical investigations ^{[9,10,11,13,14,15,16,17,18,19,20]}. Finite element and spectral methods have both been used to investigate numerical solutions ^[21,22,23].

Many numerical analyses with an adaptable mesh have shown varying outcomes. By applying open-source solvers, Warner et al. ^[15,24] are able to solve PDEs based on adaptive mesh techniques, such as PDECOL ^[25], and TOMS731 ^[26,27,28]. Sun et al. ^[29] use the h-adaptive method to solve a two-dimensional flow problem caused by gravity. Li et al. ^[30] solved a two-dimensional flow equation that was driven by gravity by combining an r-adaptive moving mesh technique with a fully discrete and nonlinear multi-grid scheme. They used an h-adaptive based on a finite-difference approach. All of the above numerical techniques were able to find the traveling wave and the fingering instabilities and fix them because they used fewer mesh points and took less time to calculate than uniform solutions. The proposed problem is solved numerically over a plane using a complete approximation repository multi-grid technique with running automated grid adaptivity ^[31].

We use the Parabolic-Monge-Ampere (PMA) method ^[32,33,34], which uses a monitor function with a smoothing monitor function and a two-dimensional PMA equation to adjust and control the movement of the mesh connected to the physical PDE that represents a two-dimensional fourth-order parabolic thin film problem with constant flux boundary conditions.

1.1.   Fourth-order $2D$ parabolic thin films problem for gravity-driven downward diffusion on an inclined, pre-wetted plane

The two-dimensional fourth-order parabolic thin film problem ^[16,35] is given by

where, ${\boldsymbol{\nabla} = \left(\frac{\partial}{\partial x}, \frac{\partial}{\partial y}\right)^T}$ and ${\boldsymbol{\nabla}^2 = \left(\frac{\partial^2}{\partial x^2}, \frac{\partial^2}{\partial y^2}\right)}$ are the gradient and the Laplacian operators, respectively. A derivation of this equation can be found in Kondic ^[36].

We emphasize that equation was created by balancing the viscous forces with the horizontal component of gravity (1.1). It is a fourth-order nonlinear degraded parabolic PDE for the development of $U(x, y, t)$. The expression degenerates because the high-derivative coefficient equals $0$ when $U$ equals $0$. The horizontal gravity component controls the fourth term in Eq (1.1). In contrast, the third term in Eq (1.1), the vertical component of gravity, frames second-order diffusion, while the second term in Eq (1.1), surface tension, frames fourth-order diffusion. The three parameters are: $D\left(\theta\right)$, which represents the equilibrium between the horizontal and vertical components of gravity, which indicates the angle of inclination of the substrate, and the capillary number, $\beta$, which is related to the idea of surface tension. $U = U(x, y, t),$ where $t > 0$ is the temporal variable, and $x$ and $y$ are space variables that depend on the direction of flow, which stands for the film height.

The boundary conditions (BCs) for Eq (1.1) are provided by

where, $\alpha\ll 1$ refers to the precursor film's thickness, $a,$ $b$ are the endpoints of $x$-direction and $c,$ $d$ are the endpoints of $y$-direction.

The initial condition is given by

where $H(x)$ is the Heaviside function, $j$ denotes the wavenumber of every sinusoidal mode with period $j/2,$ $A_j$ denotes the amplitude of the traveling wave solution, and $K_i$ denotes the breadth of the localized concerns at $x_i$. Figure 1 shows the initial condition (1.3) with the following values: $m = 2,$ $A_1 = 0,$ $A_2 = 0,$ $B_1 = 20,$ $B_2 = 100,$ precursor film thickness $0.1,$ $a = 0,$ $b = 20,$ $c = 1,$ and $d = 1.$ The capillary ridge and the point where the front connects to the precursor film, respectively, are located at $x_1$ and $x_2,$, which govern a single transverse perturbation $(j = n = 1).$

To obtain the evolution of $U(x, y, t)$ using finite difference techniques, we search the numerical results of Eq (1.1) with the constant flux border Eq (1.2) and the beginning conditions Eq (1.3). Both the Parabolic Monge-Ampere method and the uniform mesh with equal step sizes are used to carry out the numerical operations.

2.   Numerical solution of Eq (1.1) on a uniform mesh

The rectangular spatial domain $[a, \, b] \times [c, \, d]$ is uniformly split into $(M+1) \times (N+1)$ mesh points, respectively, as follows: $(x_j, y_i) = ((j-1) h_x+a, (i-1) h_y+c), \; \; j = 1, 2, \ldots, M+1, i = 1, 2, \ldots, N+1$, where ${(h_x, h_y) = \left(\dfrac{b-a}{M}, \dfrac{d-c}{N}\right) }$ is the width per sub-rectangular domain. The fundamental idea is to discretize the spatial derivatives in Eq (1.1) utilizing finite difference notations while maintaining the continuity of the temporal derivative.

Hence, a semi-discretisation of Eq (1.1) is given by

The second-term in Eq (2.1) with writing $\boldsymbol{\nabla}^2 U(x, y, t) = \Gamma(x, y, t)$ is discretised as

The third term in Eq (2.1) is discretised as

The fourth term is discretised as

where

with the same manner for $U^3_{j-\frac{1}{2}, i}$ and $U^3_{j, i-\frac{1}{2}}$. Estimating the interior points at $j = 1, M+1$ and $i = 1, N+1$ must have fictitious points $U_{0, i}$, $U_{M+2, i}$, $U_{j, 0}$, $U_{j, -1}$, $U_{j, N+2}$ and $U_{j, N+3}$ which are evaluated via discretising Eq (1.2) utilizing finite differences. Consequently, $U_{j, 0} = U_{j, 2}$, $U_{j, -1} = U_{j, 3}$, $U_{j, N+2} = U_{j, N}$, $U_{j, N+3} = U_{j, N-1}$, $U_{0, i} = U_{2, i}$ and $U_{M+2, i} = h_{M, i}$. Their ODE form returns the boundary conditions $U(x = 0, y, t) = 1$, $U(x = L_x, y, t) = \alpha$:

With a $13-$point stencil and second-order accuracy, the numerical scheme (Eqs (2.1–2.4)) is applied. As a result, when employing a centered finite difference approach on the $x$ and $y$ axes, it utilizes less bandwidth than a $26-$point stencil. One can estimate the system of ODEs that is produced by using the ODE solver DASPK ^[37]. The linearized system uses an iterative strategy based on Krylov sub-space techniques, including preconditioning via $LU-$factorizing the Jacobian matrix. It significantly affected the implementation, in contrast to other solvers like DASSL. We must select sufficient large fill-ins to decompose the Jacobian matrix using $LU.$ In this situation, the convergence of the iterative solver would be pretty slow.

The preceding numerical scheme converges as long as $(h_x, h_y)\rightarrow 0$, which we have confirmed (but have not presented here). As a result, the numerical findings for $h_x = 0.01,$ and $h_y = 0.01$ are accurate enough to be regarded as good approximation results. The surface plots of $U(x, y, t)$ at a fixed time $t = 30$ are shown in Figure 2 (a, b). It is obvious that, given the parameter values used here, the construction of a single finger agrees with the linear stability analysis. Figure 3 (a, c, f) show the surface plots of $U(x, y, t)$ exhibiting its progression over time at the indicated times of $t = 10, 20,$ and $30.$ Early $t = 10$ slow growth of a finger is observed. According to the linear stability study, the finger evolves with a preferred breadth as time t increases. In the contour plots shown in Figure 4 (b, d, e), where the finger's width is around one, and it lengthens as time $t$ grows, this is more clearly shown.

3.   Numerical results of Eq (1.1) using PMA equation

Consider the Eq (1.1) is supplemented by boundary conditions (1.2) and the initial condition (1.3). We seek the numerical results of Eq (1.1) using PMA equation ^[34]. First, we require to Express Eq (1.1) as

Here, the spatial variables $\boldsymbol x(\boldsymbol\xi, t) = (x(\xi, \eta, t), y(\xi, \eta, t))$ are depend on time and the transformation are achieved using $\boldsymbol{x} = \boldsymbol{\nabla}_{\boldsymbol\xi} P,$ by

where $\boldsymbol\xi = (\xi, \eta)$, $\nabla_{\xi} = \left(\frac{\partial}{\partial \xi}, \frac{\partial}{\partial \eta}\right)^t$ and $\boldsymbol{H}$ is the Hessian matrix given by ^[33]

Thus, the spatial derivative of Eq (3.1) is given by

The high-derivatives terms can be given by

Now, we present the PMA equation described in ^[34]:

where $\boldsymbol{x} = \left(x, y\right)$ and $\boldsymbol{\xi} = \left(\xi, \eta\right)$.

Subject to the boundary conditions:

The initial condition is given by

Here, We choose the mesh density function to be ^[32]:

Figure 5 (a, c, f) shows the surface plots for the temporal development of $U$ obtained by the PMA approach employing the curvature mesh density function (3.8). With $M = 200$ and $N = 40,$ the parameter values are provided by $\beta = 0.001$, $D(\theta) = 0$, $\alpha = 0.1$, $a = 0, \, b = 20$, $c = -1, \, d = 1$ (thus initial $h_x = 0.1$ and $h_y = 0.05$). At time $t = 10,$ a finger begins to form gradually. The finger develops over time with a specific width that is consistent with the linear stability investigation. The best representation of it may be found in Figure 6 (b, d, e), where it is approximately one finger's breadth wide and grows with time. Using the Parabolic-Monge-Ampere scheme with $M = 200,$ $N = 40$ (thus, starting $h_x = 0.1$ and $h_y = 0.05$), Figure 7 show the moving mesh at times $(a)$ $t = 10,$ $(b)$ $t = 20$, and $(c)$ $t = 30$, respectively. When the propagating finger gradually takes shape, we can plainly see the mesh adaption in both $x$ and $y$.

We use the curvature mesh density function and the following values: $a = 0, \, b = 20, \, c = 1, \, d = 1$ and $\beta = 0.001,$ $D(\theta) = 0$, $\alpha = 0.1$, (3.8). We observe the mesh transforming in both the $x$ and $y$ directions as the propagating finger gradually expands. The surface plots of $U$ at time $t = 30$ are shown in Figure 8 (a, b) using the PMA moving mesh equation with the curvature mesh density function (3.8) and $M = 200, \, N = 40.$ It is evident that the single finger visually resembles the uniform solution shown in Figure 2 (a, b) (note $M = 2000$ and $N = 200$ for the uniform solution). We can see the clustering of mesh points in the direction of y along the finger in Figure 8 (c), which shows the matching PMA moving mesh.

It can be noticed that a distinct grouping of mesh points in the $y$-direction along the finger in Figure 8 (c), which depicts the matching PMA moving mesh. Figure 9 (a, b) illustrates this as well by plotting the adaptive meshes $y$ and $x$, respectively, at time $t = 30.$ Figure 9 (b, c) shows that when $y$ is plotted as a function of $\eta$ for fixed $\xi\in[0.6, 0.8],$ it deviates from a straight line (indicating a uniform mesh) with more points clustered in the interval $[0.5, 0.5],$ which is roughly the width of a finger. The variable $y$ seems to remain uniform in the direction of $\xi$. Like $y$ changes non-uniformly as a function of, $x$ also exhibits non-linear variation, with more points concentrated around $\xi\in[0.6, 0.8],$ which is roughly the finger length scale. The variable $x$ is uniform in the direction of $\eta$. The curvature-based mesh density function at $t = 30$ is depicted in Figure 9 (c), where there are significant fluctuations in the curvature close to the finger region. As can be seen in Figure 9, this causes the mesh to be refined $(a, b).$

In Table 1, we give a number of metrics for accuracy, with a focus on the most important finger properties, so that you can check how accurate the numerical results are when the PMA moving mesh approach is used instead of the uniform mesh solution. As follows: Finger-length is the length of the finger, and finger-width is the width of the finger at a certain point (chosen as $\xi = 8.5$). % $L_2$ "Relative-Error" (using the relative error, we evaluated and compared the uniform mesh solution with $h_x = h_y = 0.01$ as the best solution.) At a certain time, $t = 30,$ we guess the amounts based on how the finger was made. They are calculated for a number of PMA meshes (designated by "PMA" in Table 1) with different values for $M$ and $N.$ We contrast the PMA solution with the equivalent values obtained using two uniform solutions (designated as "uniform" in Table 1). $M = 2000$ and $N = 200$ $(h_x = h_y = 10^2)$ are used to calculate the error as compared to the uniform solution. We use linear interpolation to figure out the error while we interpolate the solution on this uniform grid. It is noted that for PMA moving mesh solutions using between $(15-41)\times 1063$ elements, the relative inaccuracy to the uniform solution produced with $M = 2000$ and $N = 200$ is within $4.6\%$ accuracy to that of the uniform solution using $4.1\times10^5$ mesh points. The relative error, on the other hand, is about $6.5\%$ for the numerical solution using an insensitive uniform mesh $(M = 400$ and $N = 40)$ with the same number of mesh points as the PMA meshes. The PMA mesh makes an adjustment in this region in the spatial directions $x$ and $y$ (see Figure 7) to reduce the inaccuracy, which is primarily caused by the finger.

4.   Conclusions

I have successfully used the PMA moving mesh method to achieve the numerical solutions for the two-dimensional fourth-order parabolic thin film equations. These solutions are confirmed with numerical results obtained using the uniform mesh scheme with $N = 200$ and $M = 2000$. The main highlights of the solutions are presented in Figures 5, 7, 8 and Table 1, which permits a direct comparison of the uniform mesh solution obtained with $N = 200$, $M = 2000$, and the numerical solutions obtained using the PMA moving mesh method with various mesh numbers in $x$ and $y$ directions.

Through these, the numerical results obtained using PMA adaptive method are similar to a considerable extent. The mean error goes to zero whenever $h_x, h_y\to 0$. The numerical schemes for both uniform and PMA adaptive moving mesh are always stable for setting the parameter values $\beta = 0.001$, $D(\theta) = 0$, $\alpha = 0.001$, $a = 0, \, b = 20$, $c = -1, \, d = 1$.

This study effectively addressed a prototype two-dimensional thin film flow problem with surface tension and a moving contact line using the PMA moving mesh method based on mesh density functions. In order to align with the propagating finger with local refinement along the finger width, our results demonstrate how the mesh adapts in both coordinate directions (see Figure 7). It is demonstrated that, while having a significantly less number of mesh points, the PMA moving mesh method is 4.5% more accurate than one that uses a fixed uniform grid. As the non-uniform mesh is improved, we notice the convergence to a more accurate solution than a fixed uniform mesh solution by observing critical aspects of the solution behavior. The PMA moving mesh solution may be made more accurate by increasing the starting mesh points, but doing so would require more CPU processing time because extra mesh PDEs would need to be solved in addition to the underlying PDE. If the required inaccuracy is not very small, the CPU time difference could be better, and it is also significantly faster than using a fixed uniform mesh. When the initial number of points in the $y$ direction is increased from $80$ to $100,$ for instance (initial $y = 0.025$ and $0.02,$ respectively), Table 1 indicates a doubling of CPU time even if the solution appears to have converged. In this instance, the solution corresponding to the adaptive mesh with $M = 400$ and $N = 80$ may be reasonable in accuracy and computing time. Moreover, the locally revised mesh structure is similar to results previously reported for the identical test problem utilizing $h$ or $hp$-adaptive mesh refinement algorithms for unstructured meshes ^[29,30]. While we have not conducted quantitative comparisons with these studies on the accuracy of the solution, a visual examination of their solution and the meshes produced reveals very similar qualitative traits as ours. The PMA employed here offers a straightforward foundation for dynamic mesh adjustment and improvement. This contrasts with the approaches employed by Li et al. ^[30] and Sun et al. ^[29], which construct a hierarchy of multigrid meshes based on the gradient of the solution and mesh refinement based on minimizing an a posteriori error estimate and an estimate of the mesh refinement error, respectively. Using conventional methods for solving a system of parabolic partial differential equations by either the finite difference method or the finite element method, the MMPDEs may be easily implemented and solved concurrently with the thin film problem rather effectively. The thin film flow research community would significantly benefit from the PMA moving mesh approach shown here in this regard. As a result of the need to solve two additional equations for the mesh, the CPU time required by this method may be higher than that of the two studies mentioned above. This can be avoided by using mesh generation techniques inside the PMA moving mesh framework, including the r-adaptive moving mesh (MMPDEs) equations (see Budd et al. ^[38,39]), in which there is just one equation for the mesh. Future reports on these are currently being investigated.

The efficiency of PMA moving adaptive mesh approaches (compared to fixed uniform mesh schemes and maybe with h or hp-adaptive methods) and their ease of implementation are both very promising, according to our data, and they can be used often in thin-film flow problems. In this test scenario, the PMA moving mesh was not used and thus prevented mesh tangling and loss of regularity. However, this does not sound confident that other thin film flow issues generally will not also exhibit these harmful properties. We would need to investigate different mesh generation techniques inside the PMA framework, such as the so-called optimal transport equations, such as the Monge-Ampère and Parabolic Monge-Ampere (PMA) equations (Budd et al. ^[38,39]), which produce regular meshes without mesh tangling. Before their viability can be assured, these must be tested on more difficult two-dimensional thin film issues prone to more dramatic dendritic fingering instabilities, such as the fingering instabilities seen when a drop or sheet laden with surfactant spreads on a horizontal or inclined plane ^[16,17,35].

In addition to being realistic, the choice of these parameters, especially $\alpha, \, \beta \ll 1,$ makes the problem stiff and difficult to resolve numerically for the interior layers linked to the capillary ridge region. Hence, the PMA moving mesh approach, based on moving mesh PDEs, has been effectively used to solve the two-dimensional thin film equations under consideration in this study. Our findings show that this method holds tremendous promise for frequent application to thin film equations because of its ease of use and effectiveness. Before its success can be guaranteed, there are complex challenges in the thin film literature where this technology needs to be further tested.

Conflict of interest

The author declares that he has no potential conflict of interest in this article.