University of Dundee
An asymptotic fitting finite element method with exponential mesh refinement for
accurate computation of corner eddies in viscous flows
Shapeev, Alexander V.; Lin, Ping
Published in:
SIAM Journal on Scientific Computing
DOI:
10.1137/080719145
Publication date:
2009
Document Version
Publisher's PDF, also known as Version of record
Link to publication in Discovery Research Portal
Citation for published version (APA):
Shapeev, A. V., & Lin, P. (2009). An asymptotic fitting finite element method with exponential mesh refinement
for accurate computation of corner eddies in viscous flows. SIAM Journal on Scientific Computing, 31(3), 1874-
1900. https://doi.org/10.1137/080719145
General rights
Copyright and moral rights for the publications made accessible in Discovery Research Portal are retained by the authors and/or other
copyright owners and it is a condition of accessing publications that users recognise and abide by the legal requirements associated with
these rights.
• Users may download and print one copy of any publication from Discovery Research Portal for the purpose of private study or research.
• You may not further distribute the material or use it for any profit-making activity or commercial gain.
• You may freely distribute the URL identifying the publication in the public portal.
Take down policy
If you believe that this document breaches copyright please contact us providing details, and we will remove access to the work immediately
and investigate your claim.
Download date: 09. Aug. 2022
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
SIAM J. SCI. COMPUT.
c
2009 Society for Industrial and Applied Mathematics
Vol. 31, No. 3, pp. 1874–1900
AN ASYMPTOTIC FITTING FINITE ELEMENT METHOD WITH
EXPONENTIAL MESH REFINEMENT FOR ACCURATE
COMPUTATION OF CORNER EDDIES IN VISCOUS FLOWS
∗
ALEXANDER V. SHAPEEV
†
AND PING LIN
‡
Abstract. It is well known that any viscous fluid flow near a corner consists of infinite series of
eddies with decreasing size and intensity, unless the angle is larger than a certain critical angle [H. K.
Moffat, J. Fluid Mech., 18 (1964), pp. 1–18]. The objective of the current work is to simulate such
infinite series of eddies occurring in steady flows in domains with corners. The problem is approached
by high-order finite element method with exponential mesh refinement near the corners, coupled with
analytical asymptotics of the flow near the corners. Such approach allows one to compute position
and intensity of the eddies near the corners in addition to the other main features of the flow. The
method was tested on the problem of the lid-driven cavity flow as well as on the problem of the
backward-facing step flow. The results of computations of the lid-driven cavity problem show that
the proposed method computes the central eddy with accuracy comparable to the best of existing
methods and is more accurate for computing the corner eddies than the existing methods. The
results also indicate that the relative error of finding the eddies’ intensity and position decreases
uniformly for all the eddies as the mesh is refined (i.e., the relative error in computation of different
eddies does not depend on their size).
Key words. finite element method, asymptotic expansion matching, Moffatt eddies near sharp
corners
AMS subject classifications. 76D05, 76M10, 65M60
DOI. 10.1137/080719145
1. Introduction. The two-dimensional flow of a viscous fluid near the corner
between two steady rigid planes was first examined by Moffatt [28]. He established
that when the angle between planes is less than a certain critical angle, any flow near
the corner consists of infinite series of eddies with decreasing size and intensity as the
corner point is approached.
One of the most famous examples of flow in domain with corners is a flow in the
lid-driven cavity. The lid-driven cavity problem has become a benchmark problem
for researchers to test the performance of numerical methods designed for computa-
tion of viscous fluid flow. Particularly, among other criteria, the researchers examine
the accuracy of their methods based on how accurately they can compute the corner
eddies. However, in the previous works only a few eddies were computed (maximum
four corner eddies [4, 18] for certain Reynolds numbers). In addition, the accuracy of
finding intensity and position of the smaller eddies was less than the accuracy for the
larger eddies.
The only attempt known to the authors to compute a large number of corner ed-
dies for the lid-driven cavity problem is the work of Gustafson and Leben [25]. They
∗
Received by the editors March 24, 2008; accepted for publication (in revised form) December 1,
2008; published electronically March 13, 2009. This research is partially supported by the Singapore
Academic Research Funds R-146-000-064-112 and R-146-000-099-112.
http://www.siam.org/journals/sisc/31-3/71914.html
†
Department of Mathematics, National University of Singapore, 2, Science Drive 2, Singapore
117543. Current address: Lavrentyev Institute of Hydrodynamics SB RAS, 15 Lavrentyev pr.,
Novosibirsk, Russia, 630090 (alexander@shapeev.com).
‡
Department of Mathematics, National University of Singapore, 2, Science Drive 2, Singapore
117543. Current address: Division of Mathematics, University of Dundee, 23 Perth Road, Dundee,
Scotland DD1 4HN, UK (plin@maths.dundee.ac.uk).
1874
Downloaded 10/27/16 to 134.36.50.219. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
AN ASYMPTOTIC FITTING METHOD FOR EDDY COMPUTATION 1875
computed a large number of eddies (up to ten) for the Stokes flow (Re = 0) on a
sequence of subregions contracting to a corner point, setting the boundary conditions
for the smaller subregion by interpolation of solution on the larger subregions. How-
ever, their method starts with a large error due to initial grid being coarse, and this
error does not decrease when interpolating the solution onto the finer grids. Gustafson
and Leben pointed out that “Global interaction with the coarser grids is needed to
improve the solutions on all levels.” However, no works implementing this are known
to the authors of the present work.
Flows near the corner between two steady rigid planes have a weak singularity
near the corner: Flows of such type decay at a rate proportional to some power of
distance to the corner point. Therefore, the derivatives of sufficiently high degree are
not bounded in the neighborhood of the corner point. Because of these properties,
special treatment of singularities might be required to solve numerically the problem
with corner singularities.
It has been noticed that the solutions in domains with corners for problems of fluid
mechanics as well as in other disciplines have singularities which cause a slow conver-
gence rate (or sometimes divergence) of numerical methods. It has been found out
that local mesh refinement near corners and use of analytical formulas of asymptotic
solution near corners produce better results for problems with corner singularities.
Some of the popular techniques to overcome slow convergence are as follows: singu-
lar function method [19, 35], singular complement method [2], dual singular function
method [6, 7, 10, 11], introduction of analytical constraints to finite element formula-
tion [33], truncation of corners and introduction of Dirichlet-to-Neumann boundary
conditions for domains with truncated corners [21], and other methods based on the
similar ideas [26, 34]. Also, various methods based solely on mesh refinement (with-
out using asymptotic expansion of the solution) were developed (see, for example,
[1, 3, 14, 30, 31]). Mesh refinement for biharmonic boundary-value problems is dis-
cussed in [5]. Most of the works devoted to solving problems with singularities at
corners, however, either used unrefined mesh [2, 6, 10, 11], or algebraically refined
mesh[1,3,14,21,30,31].
The aim of this paper is not simply to obtain better results, but to develop a sys-
tematic method that can accurately compute position and intensity of infinite series
of eddies in addition to computing the other main features of flow in domains with
corners. The proposed method is based on the techniques developed for problems
with corner singularities, namely: Local mesh refinement near the corners and use
of asymptotic solution. The proposed local mesh refinement is exponential in the
polar radius r and uniform in the polar angle θ. A standard C
1
-continuous finite
element discretization (namely, Argyris elements) was applied to the stream function
equation. Theoretical and numerical justification of the proposed method is provided.
The proposed method was applied to the lid-driven cavity problem as well as to the
backward-facing step problem. The computations indicate that the proposed method
allows one to accurately compute the infinite series of eddies, with the relative error
of finding intensity and position of different eddies being independent of their size.
The words “asymptotic fitting” in the name of our method are motivated by the long-
existing exponential fitting method which is designed to uniformly resolve exponential
layers in singular perturbation problems (see a collection of such methods in [32]).
In this paper, by computing an infinite series of eddies we mean producing an
approximate formula of computing eddies’ intensity and position depending on the
number of the eddy. However, strictly speaking, the number of eddies we can practi-
cally compute is limited by floating point arithmetic.
Downloaded 10/27/16 to 134.36.50.219. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
1876 ALEXANDER V. SHAPEEV AND PING LIN
The structure of the paper is as follows. In section 2 we give the problem for-
mulation and discuss the properties of flows with infinite series of eddies. In section
3 we describe the proposed method for computing the infinite series of eddies. In
section 4 we present and discuss the results of computation of two problems: The lid-
driven cavity problem and the backward-facing step problem. Finally, the concluding
remarks are given in section 5.
2. Problem formulation. The problem of viscous fluid flow in domain Ω is
governed by the Navier–Stokes equations, which in 2D can be written in the form of
a single equation for the stream function ϕ:
ΔΔϕ +Re
∂Δϕ
∂x
∂ϕ
∂y
−
∂Δϕ
∂y
∂ϕ
∂x
=0, (x, y) ∈ Ω,(2.1)
where Re is the Reynolds number. This equation will be referred to as the stream
function formulation of the Navier–Stokes equations. For simplicity, we consider only
the Dirichlet boundary conditions, which cover nonslip, moving wall, and inlet/outlet
boundary conditions:
ϕ|
∂Ω
= ϕ
0
,
∂ϕ
∂n
∂Ω
= ϕ
1
, (x, y) ∈ ∂Ω,(2.2)
where ∂Ω is the boundary of Ω, and
∂
∂n
is the outward normal derivative on ∂Ω. The
variational formulation of (2.1) and (2.2) is: Find ϕ ∈ H
2
(Ω) such that
⎧
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎨
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎩
Re
Ω
∂ϕ
∂x
∂ϕ
∂y
∂
2
ψ
∂x
2
−
∂
2
ψ
∂y
2
−
∂ϕ
∂x
2
−
∂ϕ
∂y
2
∂
2
ψ
∂x∂y
dxdy
+
Ω
ΔϕΔψdxdy =0
∀ψ ∈ H
2
0
(Ω)
,
ϕ|
∂Ω
= ϕ
0
,
∂ϕ
∂n
∂Ω
= ϕ
1
.
(2.3)
The structure of the flow depends on the problem under consideration. Our
particular interest is the structure of the flow in the vicinity of the corners. As was
found by Moffatt, any flow near the corner with angle smaller than the critical one
consists of a series of eddies with decreasing size and intensity as the corner point is
approached [28]. The first (i.e., largest) eddies can be affected by the flow far from the
corner as well as by the nonlinear forces. However, such impact on the smaller eddies
can be neglected, and therefore their behavior is expected to be close to the behavior
of the family of asymptotic solutions. To summarize, the flow domain consists of
1. the part without the corner eddies,
2. the part with the relatively large corner eddies that might not be well de-
scribed by the asymptotic solution due to the impact of the flow far from the
corner as well as the impact of the nonlinear forces, and finally,
3. the part with the relatively small eddies that are well described by the asymp-
totic solution.
To compute such structure of the flow, the computational method should have specific
properties. Namely, in the first part of the domain the mesh can be uniform (unless
there are other singular features of the solution that are of interest); in the second
part the mesh should be refined in such a way that all the eddies are represented
with approximately the same number of triangles in order to compute the eddies
Downloaded 10/27/16 to 134.36.50.219. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
AN ASYMPTOTIC FITTING METHOD FOR EDDY COMPUTATION 1877
uniformly accurately; in the third part the asymptotic solution itself can be used as
a discretization.
To derive the asymptotics for the solution near the corner, following the work of
Moffatt [28], we can neglect the nonlinear terms because the velocity near the corner
between two rigid planes tends to zero. The polar coordinates, with the corner point
as the origin, can be separated in (2.1), and hence the main term in the asymptotic
solution can be found as the real part of the following complex-valued function:
ϕ = Cr
λ
f
λ
(θ),(2.4)
where
f
λ
(θ)=d
1
cos(λθ)+d
2
sin(λθ)+d
3
cos((λ − 2)θ)+d
4
sin((λ − 2)θ).
Parameters d
1
, d
2
, d
3
,andd
4
are found from the nonslip boundary conditions and λ
is defined to satisfy the Stokes equation. Particularly, for the case of right angle
f
λ
(θ)=sinθ sin ((π/2 − θ)(λ − 1)) + sin (π/2 − θ)sin(θ (λ − 1)) ,(2.5)
and λ ≈ 3.74 + 1.12i. See [8, 27] for a rigorous mathematical theory on asymptotic
expansion of the biharmonic equation near a corner.
This asymptotic solution allows one to find the asymptotic ratio of eddies’ position
and intensity, which are defined as position and value of stream function ϕ at a local
extrema. However, absolute position and intensity of eddies depend on the complex-
valued constant C which depends on the particular problem. This constant can be
found numerically for each corner of the domain. By finding the constant C,wecan
compute position and intensity of the infinite series of eddies in each corner of the
domain in the following way.
We find position (θ
k
, r
k
) and intensity (ϕ
k
) of the eddies as local extrema of the
real part of ϕ in (2.4):
∂(ϕ)
∂r
=0,
∂(ϕ)
∂θ
=0,
or after substituting (2.4):
Cλr
λ−1
f
λ
(θ)
=0,(2.6)
Cr
λ
d
dθ
f
λ
(θ)
=0.(2.7)
Here denotes the real part of a complex number. Simple analysis shows that these
equations can be satisfied only on the bisector θ = π/4, in which case (2.7) is satisfied
automatically. Hence r can be found by substituting θ = π/4into(2.6):
r
λ−1
Cλf
λ
(π/4)
=0.
Taking into account that f
λ
(π/4) =0and
r
λ−1
= e
(λ−1) ln r
= e
((λ)−1) ln r
(cos ((λ)lnr)+i sin ((λ)lnr)),
where is the imaginary part of a complex number, we can write position (r
k
,θ
k
)of
the eddies as
θ
k
=
π
4
,r
k
= e
(λ)
(
−πk+arccot
(
arg
(
Cλf
λ
(
π
4
)
)))
,(2.8)
Downloaded 10/27/16 to 134.36.50.219. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
