A new adaptive local mesh refinement algorithm and its application on fourth order thin film flow problem

https://doi.org/10.1016/j.jcp.2006.11.005Get rights and content

Abstract

A new adaptive local mesh refinement method is presented for thin film flow problems containing moving contact lines. Based on adaptation on an optimal interpolation error estimate in the Lp norm (1 < p  ∞) [L. Chen, P. Sun, J. Xu, Multilevel homotopic adaptive finite element methods for convection dominated problems, in: Domain Decomposition Methods in Science and Engineering, Lecture Notes in Computational Science and Engineering 40 (2004) 459–468], we obtain the optimal anisotropic adaptive meshes in terms of the Hessian matrix of the numerical solution. Such an anisotropic mesh is optimal for anisotropic solutions like the solution of thin film equations on moving contact lines. Thin film flow is described by an important type of nonlinear degenerate fourth order parabolic PDE. In this paper, we address the algorithms and implementation of the new adaptive finite element method for solving such fourth order thin film equations. By means of the resulting algorithm, we are able to capture and resolve the moving contact lines very precisely and efficiently without using any regularization method, even for the extreme degenerate cases, but with fewer grid points and degrees of freedom in contrast to methods on a fixed mesh. As well, we compare the method theoretically and computationally to the positivity-preserving finite difference scheme on a fixed uniform mesh which has proven useful for solving the thin film problem.

Introduction

Adaptive procedures for the numerical solution of partial differential equations (PDEs), actively investigated since the late 1970s, are now standard tools in science and engineering – e.g. see [71] for references on adaptivity for elliptic PDEs. Adaptive finite element methods (FEMs) are a particularly meaningful approach for handling multiscale phenomena and making realistic computations feasible, especially on irregular domains and in higher spatial dimensions with complex boundary conditions, where a posteriori error estimators are available as an essential ingredient of adaptivity. Such estimators are computable quantities depending on the computed solution(s) and data which provide information about the quality of approximation and may thus be used to make judicious mesh modifications. The ultimate purpose is to construct a sequence of meshes which will eventually equidistribute the approximation errors and, as a consequence, the computational effort. To this end, the a posteriori error estimators are split into element indicators which are then employed to make local mesh modifications by refinement (and sometimes coarsening). This naturally leads to loops of the formSolveEstimateRefine.Starting from a coarse mesh, such an iteration has been widely successful in applications. Nevertheless, except for the rather complete description of the one-dimensional situation by Babuška and Rheinboldt [3], convergence of (1.1) in the multidimensional case is still largely an open issue. The fundamental paper [38] of Dörfler for the Poisson equation shows a linear error reduction rate for the energy norm towards a preassigned tolerance in finite steps. Recently Nochetto et al. [60] have constructed an adaptive FEM algorithm for elliptic PDEs and proved its linear convergence rate for the energy norm. Any prescribed error tolerance is thus achieved in a finite number of steps. It is in this context that there are some of the best current convergence results for adaptive finite element methods. All of the above a posteriori error estimators fall into a class of residual-type methods because they are all based on residual error on each element.

Besides the above developments for a posterior error estimators’ convergence analysis, there is an alternative approach to produce a similar effect, viz., constructing (nearly) optimal meshes for suitable order piecewise finite element interpolation of a given function through a priori interpolation error estimators. Specifically, let ΩRn be a bounded domain, T a simplicial finite element mesh of Ω¯ with a fixed number N of elements, and uI a piecewise finite element interpolation of a given function u defined on Ω¯. An optimal mesh could be obtained by minimizing the error |||u  uI||| in some sense, where the norm ||| · ||| is a classical Sobolev space norm.

This approach can be traced back to de Boor [34], [33] where the problem of the best approximation by free knots splines was studied in one spatial dimension. In this work, the equidistribution principle was introduced specifically for computing equidistributed meshes. Actually the concept of equidistribution was first used by Burchard [21], and then by a number of researchers for studying grid adaptation. Pioneering work for adaptive finite element methods was done in [3] where a finite element mesh was shown nearly optimal in the sense of minimizing the H1 norm error if the local errors are approximately equal for all elements. Thus, to get an optimal mesh, elements where the error is large are marked for refinement, while elements with a small error are left unchanged or coarsened.

Most of the adaptive finite element methods in the literature [4] are concerned with meshes that are shape-regular (which in two dimensions means that no element has a very small angle). This type of shape-regular finite elements is appropriate for physical problems that are fairly isotropic. But for many anisotropic problems (e.g. with sharp boundary layers or internal layers), the shape of elements can be further optimized, and an equidistribution of a scalar error density is not sufficient to ensure that a mesh is optimally efficient [32]. Nadler [62] studied the optimal triangulation for the discontinuous piecewise linear approximation for a quadratic function in the sense of minimizing error in the L2 norm. For this optimal mesh, each triangle is equilateral under the Hessian metric, and the error is equidistributed on each triangle. The L case was studied by D’Azevedo and Simpson [30], [31]. In [31] they developed a local linear interpolation error formula for quadratic functions. Based on this formula, D’Azevedo [30] obtained the same condition as Nadler’s. Thereafter anisotropic mesh adaptation which aims to generate equilateral triangles under the metric induced by Hessian matrix was developed in [20], [48], [37] and successfully applied to computational fluid dynamic problems in two spatial dimensions [48], [70].

Recently there have been some a priori interpolation error estimates for anisotropic finite elements [2], [56], [42]. Apel [2] obtained some estimates under a condition on the coordinate orientation and on the maximal allowable mesh angle. Formaggia and Perotto [42] exploited the spectral properties of the affine map from the reference triangle to the general triangle to get anisotropic estimates for the L2 and H1 interpolation error on linear finite elements in two dimensions. Kunert [56] introduced the matching function to measure the alignment of an anisotropic function and an anisotropic mesh and presented error estimates using the matching function. Yet the overall optimal convergence rate in terms of the number of degree of freedom is not easy to get from these approaches.

Motivated by the fundamental papers of Huang [49], [50], a general result which develops this approach of a priori interpolation error estimators further is given by Chen et al. [28]. In this paper, a new local edge-based error estimate and global interpolation error bound in the Lp (1  p  ∞) norm for any spatial dimension are given. Specifically,u-uILp(Ω)CN-2/ndet(H)nLpn2p+n(Ω),where N is the number of elements in the triangulation, the constant C is independent of u and N, and H is a majorizing Hessian matrix of u. This estimate is optimal in the sense that it is a lower bound if u is strictly convex or concave.

This estimate can be viewed as a modification and generalization of the special case p = ∞ in D’Azevedo [30] and the case p = 2 in Huang and Sun [50]. By requiring the mesh to be quasi-uniform under the new metric [det(H)]−1/(2p+n)H, the simplices become locally isotropic and their volumes are globally equidistributed under the new metric. On the other hand, the same simplices may simultaneously transform back to be locally anisotropic under the standard Euclidean metric if the majorizing Hessian matrix is anisotropic, or more specifically, if its eigenvalues differ markedly from each other.

The estimate (1.2) is the theoretical foundation of our method, as the new adaptive local mesh refinement algorithm aims at minimizing (or at least reducing) this interpolation error by iteratively modifying the grids. Then the corresponding generated adaptive mesh is able to anisotropically concentrate on the sharp interface or boundary layer so that the solution singularities are fully resolved. The thin film fluid flow problem is exactly such a case where the solution has anisotropic sharp internal layers.

While some numerical results are given in [27], [28], no algorithm discussing how to adequately implement the adaptive mesh method based on the new interpolation error estimator (1.2) is included. In this paper, we describe a fairly detailed implementation for producing anisotropic adaptive locally refined meshes and the techniques of local mesh improvement which aim at minimizing the interpolation error. We adopt the thin film fluid flow problem, being typical of many realistic problems, to test our new adaptive finite element method.

Thin film flow problems, which lead to fourth order nonlinear degenerate parabolic equations [45], [19], [12], [53], represent just one type of interesting fourth order problems which have been of considerable recent interest and for which adaptive numerical methods need to be developed and analyzed. One reason we study them for our application is that not only do their solutions include natural moving contact lines – one kind of internal layers which need to be fully resolved in order to get accurate numerical solutions – but also because the solution is relatively simply distributed, being almost constant except around the moving contact lines. We can then use a local refined mesh only near contact lines and coarse meshes elsewhere. In contrast, fixed uniform meshes can be unacceptable in higher spatial dimensions because of the relatively large spatial scale as compared to the slim contact line regions and the long time scale (starting from the initial data to rupture and afterwards). Our adaptive local anisotropic mesh refinement algorithm based on minimizing the interpolation error (1.2) is able to efficiently overcome these difficulties and accurately solve the problem.

The paper is organized as follows. In Section 2 we present an interpolation error estimate and show that it is also a lower bound for convex or concave functions. In Section 3 we describe in detail the numerical implementation of the adaptive mesh method based upon the above interpolation error estimate, and we develop local mesh improvement techniques such as refinement, coarsening and smoothing strategies which aim to minimize the interpolation error. In Section 4 we introduce the thin film fluid flow problem, give its background and solution properties, and review previous efforts for solving it numerically. An implementation of our new adaptive local refinement method is given in Section 5. Some concluding remarks are made in Section 6.

Section snippets

Theoretical foundation

In this section, we give an interpolation error estimate derived by Chen et al. [28] which provides the theoretical foundation for our algorithm based upon minimizing (or at least reducing) this interpolation error by iteratively modifying the grids.

The estimate. Let Ω be a bounded domain in Rn. Given a function uC2(Ω¯), a symmetric positive definite matrix H  Rn×n is majorizing the Hessian of u if|ξt(2u)(x)ξ|c0ξtH(x)ξ,ξRn,xΩ,for some positive constant c0.

One example of H can be constructed

General implementation issues

In this section, we elaborate on our new adaptive local mesh refinement method. Recall that our goal is to construct a quasi-uniform mesh under the new metric (2.3). Our strategy is to do the local refinement on a coarse grid by equidistributing the interpolation error in the whole triangulation, viz., we make the edge lengths of each triangle element under the new metric equal to each other. In grid generation, one often begins with a uniform or regular (in the Euclidean metric) grid which

Thin film fluid flows

Many physical phenomena exhibit strong anisotropic behaviors, e.g. boundary layer flows in porous media, currents and concentrations in fuel cells, characteristics of semiconductors, and stresses and strains in thin plates, shells and anisotropic materials. When an anisotropic phenomenon occurs and an accurate approximation is required, it is natural to use stretched anisotropic meshes according to the variation of the solution. This reduces the number of elements needed to partition the

Mixed finite element method for fourth order PDE

In principle, there are a number of natural finite element methods to numerically solve fourth order PDEs. One is the standard FEM with high order piecewise polynomials which include derivatives as degrees of freedom – say, the cubic Hermite elements; another is the mixed finite element method. Considering the difficulties of constructing a Hermite element in an arbitrary triangular element or a tetrahedron element in the 3D case, we prefer the mixed finite element method. For a fourth order

The consistency of FEM and PPS

In [78] the authors indicate that it is possible to generalize the positivity-preserving scheme (PPS) to finite element methods on arbitrary element spaces (including those involving nonuniform grids). They begin by showing that the 1D PPS finite difference scheme is equivalent to a specific finite element discretization in which a nonlinear function of the solution is represented in the element basis. This representation suggests a general abstract finite element approximation of the problem

Time stepping

Explicit time stepping procedures present a major challenge for fourth order equations, with the general time step restriction that Δt be on the order of h4. This restriction, far worse than the corresponding stability requirement for second order problems, is particularly severe when fine grids are needed to obtain highly resolved solutions accurately, as is the case for the thin film problem. Reduction of the thin film equation from fourth order form to a second order system (5.1) is critical

Numerical results

In this section, we consider the thin film fluid flow equation (4.4) where α = 1, β = 0, γ = 1, and consider the thin film dynamics in a rectangle domain whose lateral dimension is comparable to the wavelength of maximum growth. We take the domain size as 100 × 16 in the following simulations.

The experiments and numerical simulations done in [54] demonstrate that, when a thin film flows down an inclined plane, after some time the initially straight line where the liquid, gas and solid phase meet

Conclusions

In this paper, we have presented a detailed implementation of an adaptive finite element method (AFEM) and demonstrated that it can perform very well on a challenging problem like the thin film problem with a moving contact line. The results for our anisotropic AFEM are consistent with those obtained previously, most using the positivity-preserving finite difference method [19], [14], [1], [53], [35], [54]. It requires much fewer grid points for comparable accuracy than with the uniform grid

Acknowledgements

P.S. and R.D.R. are partially supported by NSERC Grant A8781, Canada. J.X. is partially supported by NSF Grant #DMS-0074299 and Center for Computational Mathematics and Applications, Penn State.

References (80)

  • T. Apel, Anisotropic finite elements: local estimates and applications, in: J.P. Pier (Ed.), Advances in Numerical...
  • I. Babuška et al.

    Error estimates for adaptive finite element computations

    SIAM J. Numer. Anal.

    (1978)
  • I. Babuška et al.

    The finite element method and its reliability

  • R.E. Bank et al.

    Refinement algorithms and data structures for regular local mesh refinement

  • R.E. Bank et al.

    Mesh smoothing using a posteriori error estimates

    SIAM J. Numer. Anal.

    (1997)
  • R.E. Bank et al.

    Asymptotically exact a posteriori error estimators, part I: grids with superconvergence

    SIAM J. Numer. Anal.

    (2003)
  • R.E. Bank et al.

    Asymptotically exact a posteriori error estimators, part II: general unstructured grids

    SIAM J. Numer. Anal.

    (2003)
  • J. Barrett et al.

    Finite element approximation of a fourth order degenerate parabolic equation

    Numer. Math.

    (1998)
  • J. Becker et al.

    The thin film equation: recent advances and some new perspectives

    J. Phys. Condens. Mat.

    (2005)
  • E. Beretta et al.

    Nonnegative solutions of a fourth order nonlinear degenerate parabolic equation

    Arch. Rational Mech. Anal.

    (1995)
  • F. Bernis, Viscous flows, fourth order nonlinear degenerate parabolic equations and singular elliptic problems, in free...
  • A. Bertozzi

    Loss and gain of regularity in a lubrication equation for thin viscous films, in Free boundary problems: theory and applications

  • A. Bertozzi et al.

    Stability of compressive and undercompressive thin film travelling waves

    Eur. J. Appl. Math.

    (2001)
  • A. Bertozzi et al.

    The lubrication approximation for thin viscous films: regularity and long time behaviour of weak solutions

    Commun. Pure Appl. Math.

    (1996)
  • A.L. Bertozzi

    Symmetric singularity formation in lubrication-type equations for interface motion

    SIAM J. Appl. Math.

    (1996)
  • A.L. Bertozzi et al.

    Linear stability and transient growth in driven contact lines

    Phys. Fluids

    (1997)
  • A.L. Bertozzi et al.

    Singularities and similarities in interface flow

  • H. Borouchaki et al.

    Anisotropic adaptive mesh generation in two dimensions for cfd

  • H.G. Burchard

    Splines (with optimal knots) are better

    Appl. Anal.

    (1974)
  • J. Cahn et al.

    The Cahn–Hilliard equation with a concentration dependent mobility: motion by minus the laplacian of the mean curvature

    Eur. J. Appl. Math.

    (1996)
  • P. Carles et al.

    On the origin of the bump in the profile of surface-tension-gradient-driven spreading films

    Mat. Res. Soc. Symp. Proc.

    (1990)
  • A.M. Cazabat et al.

    Finger instability of this spreading films driven by temperature gradients

    Nature

    (1990)
  • C.M. Chen et al.

    High Accuracy Theory of Finite Element Methods

    (1995)
  • L. Chen, Robust and accurate algorithms for solving anisotropic singularities, Ph.D. thesis, Pennsylvania State...
  • L. Chen, P. Sun, J. Xu, Multilevel homotopic adaptive finite element methods for convection dominated problems, in:...
  • L. Chen et al.

    Optimal anisotropic simplicial meshes for minimizing interpolation errors in Lp-norm

    Math. Comp.

    (2007)
  • P. Constantin et al.

    Droplet breakup in a model of the Hele–Shaw cell

    Phys. Rev. E

    (1993)
  • E. D’Azevedo

    Optimal triangular mesh generation by coordinate transforma tion

    SIAM J. Sci. Stat. Comp.

    (1991)
  • E. D’Azevedo et al.

    On optimal interpolation triangle incidences

    SIAM J. Sci. Stat. Comp.

    (1989)
  • E.F. D’Azevedo

    Optimal triangular mesh generation by coordinate transformation

    SIAM J. Sci. Stat. Comp.

    (1991)
  • Cited by (16)

    • Simulation of thin film flows with a moving mesh mixed finite element method

      2018, Applied Mathematics and Computation
      Citation Excerpt :

      One is the standard finite element method (FEM) with high order piecewise polynomials which includes derivatives as degrees of freedom. Another approach is to use a low order basis such as the mixed finite element method [13]. Now we discretize the system (8) in space using the standard piece-wise linear finite element method (FEM) and we use the same space for both variables.

    • An adaptive moving mesh method for thin film flow equations with surface tension

      2017, Journal of Computational and Applied Mathematics
    • Numerical investigations on self-similar solutions of the nonlinear diffusion equation

      2013, European Journal of Mechanics, B/Fluids
      Citation Excerpt :

      Adaptive mesh refinement for thin film equations is developed in [22]. A detailed implementation of an adaptive finite element method was presented in [23]. In [24], the authors numerically investigated the effect of the convection term treatment using the Godunov scheme, the WENO scheme, and an upwind-type scheme of a driven thin film equation.

    • The cutoff method for the numerical computation of nonnegative solutions of parabolic PDEs with application to anisotropic diffusion and Lubrication-type equations

      2013, Journal of Computational Physics
      Citation Excerpt :

      Russell et al. [44] solve the same regularized equation using a moving collocation method. Sun et al. [48] solve (24) in two dimensions using an adaptive finite element method and show that proper mesh adaptation can provide accurate resolution and there is no need to use regularization in the numerical simulation of the singularity development in (24). In the previous section we have studied the cutoff method for the numerical computation of nonnegative solutions of parabolic partial differential equations.

    View all citing articles on Scopus
    View full text