A new adaptive local mesh refinement algorithm and its application on fourth order thin film flow 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 formStarting 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 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,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 , a symmetric positive definite matrix H ∈ Rn×n is majorizing the Hessian of u iffor 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)
- et al.
Higher order nonlinear degenerate parabolic equations
J. Differ. Equat.
(1990) - et al.
Computing three-dimensional thin film flows including contact lines
J. Comp. Phys.
(2002) - et al.
Diffusional phase transitions in multicomponent systems with a concentration dependent mobility matrix
Physica D
(1997) Variational mesh adaptation: isotropy and equidistribution
J. Comput. Phys.
(2001)- et al.
Variational mesh adaptation II: error estimates and monitor functions
J. Comput. Phys.
(2003) - et al.
On adaptive grid refinement in the presence of internal or boundary layers
IMPACT Comput. Sci. Eng.
(1990) Software for C1 surface interpolation
Spatial coupling of plant and herbivore dynamics: the contribution of herbivore dispersal to transient and persistent waves of damage
Theor. Popul. Biol.
(1994)- et al.
ADI schemes for higher-order nonlinear diffusion equations
Appl. Numer. Math.
(2003) - et al.
Stable and unstable singularities in the unforced hele-shaw cell
Phys. Fluids
(1996)
Error estimates for adaptive finite element computations
SIAM J. Numer. Anal.
The finite element method and its reliability
Refinement algorithms and data structures for regular local mesh refinement
Mesh smoothing using a posteriori error estimates
SIAM J. Numer. Anal.
Asymptotically exact a posteriori error estimators, part I: grids with superconvergence
SIAM J. Numer. Anal.
Asymptotically exact a posteriori error estimators, part II: general unstructured grids
SIAM J. Numer. Anal.
Finite element approximation of a fourth order degenerate parabolic equation
Numer. Math.
The thin film equation: recent advances and some new perspectives
J. Phys. Condens. Mat.
Nonnegative solutions of a fourth order nonlinear degenerate parabolic equation
Arch. Rational Mech. Anal.
Loss and gain of regularity in a lubrication equation for thin viscous films, in Free boundary problems: theory and applications
Stability of compressive and undercompressive thin film travelling waves
Eur. J. Appl. Math.
The lubrication approximation for thin viscous films: regularity and long time behaviour of weak solutions
Commun. Pure Appl. Math.
Symmetric singularity formation in lubrication-type equations for interface motion
SIAM J. Appl. Math.
Linear stability and transient growth in driven contact lines
Phys. Fluids
Singularities and similarities in interface flow
Anisotropic adaptive mesh generation in two dimensions for cfd
Splines (with optimal knots) are better
Appl. Anal.
The Cahn–Hilliard equation with a concentration dependent mobility: motion by minus the laplacian of the mean curvature
Eur. J. Appl. Math.
On the origin of the bump in the profile of surface-tension-gradient-driven spreading films
Mat. Res. Soc. Symp. Proc.
Finger instability of this spreading films driven by temperature gradients
Nature
High Accuracy Theory of Finite Element Methods
Optimal anisotropic simplicial meshes for minimizing interpolation errors in Lp-norm
Math. Comp.
Droplet breakup in a model of the Hele–Shaw cell
Phys. Rev. E
Optimal triangular mesh generation by coordinate transforma tion
SIAM J. Sci. Stat. Comp.
On optimal interpolation triangle incidences
SIAM J. Sci. Stat. Comp.
Optimal triangular mesh generation by coordinate transformation
SIAM J. Sci. Stat. Comp.
Cited by (16)
An adaptive moving mesh method for two-dimensional thin film flow equations with surface tension
2019, Journal of Computational and Applied MathematicsSimulation of thin film flows with a moving mesh mixed finite element method
2018, Applied Mathematics and ComputationCitation 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 MathematicsNumerical investigations on self-similar solutions of the nonlinear diffusion equation
2013, European Journal of Mechanics, B/FluidsCitation 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 PhysicsCitation 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.
A numerical scheme for particle-laden thin film flow in two dimensions
2011, Journal of Computational Physics