Toward the Finite-Time Blowup of the 3D Axisymmetric Euler Equations: A Numerical Investigation

©

by SIAM. Unauthorized reproduction

of this article is prohibited.

ULTISCALE

ODEL.

IMUL

2014 Society for I

ndustrial and Applied Mathematics

Vol. 12, No. 4, pp. 1722–1776

TOWARD THE FINITE-TIME BLOWUP OF THE 3D

AXISYMMETRIC EULER EQUATIONS: A NUMERICAL

INVESTIGATION

∗

GUO LUO

†

AND

THOMAS Y. HOU

‡

Abstract.

Whether the three-dimensional incompressible Euler equations can develop a singu-

larity in finite time from smooth initial data is one of the most challenging problems in mathematical

fluid dynamics. This work attempts to provide an affirmative answer to this long-standing open ques-

tion from a numerical point of view by presenting a class of potentially singular solutions to the Euler

equations computed in axisymmetric geometries. The solutions satisfy a periodic boundary condi-

tion along the axial direction and a no-flow boundary condition on the solid wall. The equations

are discretized in space using a hybrid 6th-order Galerkin and 6th-order finite difference method on

specially designed adaptive (moving) meshes that are dynamically adjusted to the evolving solutions.

With a maximum effective resolution of over (3

)

near the point of the singularity, we are

able to advance the solution up to

003505 and predict a singularity time of

≈

0035056,

while achieving a

pointwise

relative error of

(10

−

) in the vorticity vector

and observing a

)-fold increase in the maximum vorticity

∞

. The numerical data are checked against

all major blowup/non-blowup criteria, including Beale–Kato–Majda, Constantin–Fefferman–Majda,

and Deng–Hou–Yu, to confirm the validity of the singularity. A local analysis near the point of the

singularity also suggests the existence of a self-similar blowup in the meridian plane.

Key words.

3D axisymmetric Euler equations, finite-time blowup

AMS subject classifications.

35Q31, 76B03, 65M60, 65M06, 65M20

DOI.

10.1137/140966411

1. Introduction.

The celebrated three-dimensional (3D) incompressible Euler

equations in fluid dynamics describe the motion of ideal incompressible flows in the

absence of external forcing. First written down by Leonhard Euler in 1757, these

equations have the form

(1.1)

·∇

−∇

∇·

where

)

is the 3D velocity vector of the fluid and

is the scalar pres-

sure. The 3D Euler equations have a rich mathematical theory, for which the inter-

ested readers may consult the excellent surv

eys [2, 18, 24] and the references therein.

This paper primarily concerns

the existence or nonexistence of globally regular solu-

tions to the 3D Euler equations, which is regarded as one of the most fundamental

yet most challenging problems in mathematical fluid dynamics.

The interest in the global regularity or finite-time blowup of (1.1) comes from

several directions. Mathematically, the question has remained open for over 250 years

and has a close connection to the Clay Millennium Prize Problem on the Navier–

Stokes equations. Physically, the formation of a singularity in inviscid (Euler) flows

may signify the onset of turbulence in viscous (Navier–Stokes) flows, and it may

∗

Received by the editors April 24, 2014; accepted for publication (in revised form) August 15,

2014; published electronically November 18, 2014.

http://www.siam.org/journals/mms/12-4/96641.html

†

Applied & Computational Mathematics, California Institute of Technology, MC 9-94, Pasadena,

CA 91125. Current address: Department of Mathematics, City University of Hong Kong, Kowloon

Tong, Hong Kong (gluo@caltech.edu).

‡

Applied & Computational Mathematics, California Institute of Technology, MC 9-94, Pasadena,

CA 91125 (hou@acm.caltech.edu).

1722

Downloaded 01/29/15 to 131.215.70.231. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php

©

by SIAM. Unauthorized reproduction

of this article is prohibited.

FINITE-TIME SINGULARITY OF 3D EULER

1723

provide a mechanism for energy transfer to small scales. Numerically, the resolution

of nearly singular flows requires special numerical treatment, which presents a great

challenge to computational fluid dynamicists.

Considerable efforts have been devoted to the study of the regularity properties of

the 3D Euler equations. The main difficulty i

n the analysis lies in the presence of the

nonlinear vortex stretching term and the lack of a regularization mechanism, which

implies that even the local well-posedness of the equations can only be established

for sufficiently smooth initial data (see, for e

xample, [37]). Despite these difficulties,

a few important partial results [3, 44, 22, 46, 19, 20, 25] have been obtained over the

years which have led to improved understanding of the regularity properties of the 3D

Euler. More specifically, the celebrated theorem of Beale, Kato, and Majda [3] and

its variants [22, 46] characterize the regularity of the 3D Euler equations in terms of

the maximum vorticity, asserting that a smooth solution

of (1.1) blows up at

if and only if

(

)

∞

where

∇×

is the

vorticity vector

of the fluid. The non-blowup criterion of

Constantin, Fefferman, and Majda [19] focu

ses on the geometric aspects of Euler flows

instead and asserts that there can be no blowup if the velocity field

is uniformly

bounded and the vorticity direction

ω/

is sufficiently “well behaved” near the

point of the maximum vorticity. The theorem of Deng, Hou, and Yu [20] is similar

in spirit to the Constantin–Fefferman–Majda criterion but confines the analysis to

localized vortex line segments.

Besides the analytical results mentioned

above, there also exists a sizable litera-

ture focusing on the (numerical) search of a finite-time singularity for the 3D Euler

equations. Representative work in this direction include [27, 45], which studied Euler

flows with swirls in axisymmetric geometries, the famous computation of Kerr and his

collaborators [38, 8, 39], which studied Eu

ler flows generated by a pair of perturbed

antiparallel vortex tubes, and the viscous simulations of [5], which studied the 3D

Navier–Stokes equations using Kida’s high-symmetry initial data. Other interesting

pieces of work are [10, 47], which studied axisymmetric Euler flows with

complex

ini-

tial data and reported singularities in the complex plane. A more comprehensive list

of interesting numerical results can be found in the review article [24].

Although finite-time singularities were frequently reported in numerical simula-

tions of the Euler equations, most such singularities turned out to be either numerical

artifacts or false predictions, as a result of either insufficient resolution or inadver-

tent data analysis procedure (more to follow on this topic in section 4.4). Indeed, by

exploiting the analogy between the two-dimensional (2D) Boussinesq equations and

the 3D axisymmetric Euler equations away from the symmetry axis, E and Shu [21]

studied the potential development of finite-time singularities in the 2D Boussinesq

equations, with initial data completely analogous to those of [27, 45]. They found no

evidence for singular solutions, indicating that the “blowups” reported by [27, 45],

which were located away from the axis, are likely to be numerical artifacts. Likewise,

Hou and Li [32] repeated the computation of [38] with higher resolutions, in an at-

tempt to reproduce the singularity observed in that study. Despite some ambiguity in

interpreting the initial data used by [38], they managed to advance the solution up to

= 19, which is beyond the singularity time

=18

7 alleged by [38]. By using newly

developedanalytictoolsbasedonrescaledvorticitymoments,Kerralsoconfirmedin

Downloaded 01/29/15 to 131.215.70.231. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php

©

by SIAM. Unauthorized reproduction

of this article is prohibited.

1724

GUO LUO AND THOMAS Y. HOU

a very recent study [39] that solutions com

puted from initial data analogous to that

used in [38] eventually converge to supere

xponential growth and hence are unlikely

to lead to a singularity. In a later work, Hou and Li [33] also repeated the computa-

tion of [5] and found that the singularity reported by [5] is likely an artifact due to

insufficient resolution.. . . To summarize, the existing numerical studies do not seem

to provide convincing evidence to support the existence of a finite-time singularity,

and the question of whether initially smooth solutions to (1.1) can blow up in finite

time remains open.

By focusing on solutions with axial symmetry and other special properties, we

have discovered, through careful numerical studies, a class of potentially singular

solutions to the 3D

axisymmetric

Euler equations in a radially bounded, axially pe-

riodic cylinder (see (2.1)–(2.2) below). The reduced computational complexity in the

cylindrical geometry greatly facilitates the computation of the singularity. With a

specially designed adaptive mesh, we are able to achieve a maximum mesh resolution

of over (3

)

near the point of the singularity. This allows us to compute the

vorticity vector with four digits of accuracy throughout the simulations and to ob-

serve a (3

)-fold amplification in maximum vorticity. The numerical data are

checked against all major blowup/non-blowup criteria, including Beale–Kato–Majda,

Constantin–Fefferman–Majda, and Deng–Hou–Yu, to confirm the validity of the sin-

gularity. A careful local analysis also suggests the existence of a self-similar blowup

in the meridian plane. Our numerical method makes explicit use of the special sym-

metries built in the blowing-up solutions, which eliminates symmetry-breaking per-

turbations and facilitates a stable computation of the singularity.

The main features of the potentially singular solutions are summarized as follows.

The point of the potential singularity, which is also the point of the maximum vorticity,

is always located at the intersection of the solid boundary

= 1 and the symmetry

plane

=0. Itisa

stagnation point

of the flow, as a result of the special odd-even

symmetries along the axial direction and the no-flow boundary condition (see (2.3)).

The vanishing velocity field at this point could have positively contributed to the

formation of the singularity, given the potential regularizing effect of convection as

observed by [33, 31]. When viewed in the meridian plane, the point of the potential

singularity is a

hyperbolic saddle

of the flow, where the axial flow along the solid

boundary marches toward the symmetry plane

= 0 and the radial flow marches

toward the symmetry axis

= 0 (see Figure 16(a)). The axial flow brings together

vortex lines near the solid boundary

= 1 and destroys the geometric regularity

of the vorticity vector near the symmetry plane

= 0, violating the Constantin–

Fefferman-Majda and Deng–Hou–Yu geometric non-blowup criteria and leading to

the breakdown of the smooth vorticity field.

The asymptotic scalings of the various quantities involved in the potential finite-

time blowup are summarized as follows. Near the predicted singularity time

,the

scalar pressure and the velocity field remain uniformly bounded while the maximum

vorticity blows up like

(

−

)

−

,where

roughly equals

. Near the point of the

potential singularity, namely the point of the maximum vorticity, the radial and axial

components of the vorticity vector grow roughly like

(

−

)

−

while the angular

vorticity grows like

(

−

)

−

. The nearly singular solution has a locally self-similar

structure in the meridian plane near the point of blowup, with a rapidly collapsing

support scaling roughly like

(

−

)

along both the radial and the axial directions.

When viewed in

, this corresponds to a thin tube on the symmetry plane

evolved around the ring

= 1, where the radius of the tube shrinks to zero as the

singularity forms.

Downloaded 01/29/15 to 131.215.70.231. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php

©

by SIAM. Unauthorized reproduction

of this article is prohibited.

FINITE-TIME SINGULARITY OF 3D EULER

1725

We emphasize that the 3D axisymmetric Euler equations (2.1) are

different

from

their free-space counterpart (1.1) in that they have a constant of motion that is not

present in the nonsymmetric case [41]. In addition, it is well known that the choice

of the boundary conditions (periodic vs. no-flow) has a nontrivial impact on the

qualitative behavior of the solutions of the Euler equations, especially near the solid

boundaries [2, 18]. In view of these differences and the fact that the singularity we

discover lies right on the boundary, we stress that the work described in this paper

not

directly relevant to the Clay Millennium Prize Problem on the Navier–Stokes

equations, which is posed either in free space or on periodic domains.

Rather, it

should be viewed as an attempt at the understanding of the effect of solid bound-

aries in the creation of small scales and, in

the case of zero viscosity, the creation of

singularities in incompressible flows.

The rest of this paper is devoted to the study of the potential finite-time singu-

larity and is organized as follows. Sectio

n 2 contains a brief review of the 3D Euler

equations in axisymmetric form and defines the problem to be studied. Section 3

gives a brief description of the numerical method that is used to track and resolve the

nearly singular solutions. Section 4 examines the numerical data in great detail and

presents evidence supporting the existence of

a finite-time singularity. Finally section

5 concludes the paper with a brief discussion on future research directions.

2. Description of the problem.

The 3D Euler equations (1.1) with axial sym-

metry can be conveniently described in the so-called stream-vorticity form. To derive

these equations, recall first that in cylindrical coordinates (

r, θ, z

), an axisymmetric

flow

can be described by

the decomposition

(

r, z

(

r, z

)

(

r, z

)

(

r, z

)

where

=(cos

θ,

sin

θ,

−

sin

θ,

cos

θ,

,and

=(0

are coordi-

nate axes. The vorticity vector

∇×

has a similar representation,

(

r, z

(

r, z

)

(

r, z

)

(

r, z

)

−

,ω

−

,ω

(

)

where for simplicity we have used subscripts to denote partial differentiations. The

incompressibility condition

∇·

= 0 implies the existence of a stream function

(

r, z

(

r, z

)

(

r, z

)

(

r, z

)

for which

∇×

and

−

.Takingthe

-components of the velocity equation

(1.1), the vorticity equation

·∇

and the Poisson equation

−

gives an alternative formulation of the 3D Euler

equations

(2.1a)

)

(2.1b)

−

∂

+(3

)

∂

(2.1c)

Indeed, according to the partial regularity result of Caffarelli, Kohn, and Nirenberg [9], any

finite-time singularity of the 3D axisymmetric Navier–Stokes equations, if it exists, must lie on the

symmetry axis.

Downloaded 01/29/15 to 131.215.70.231. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php

©

by SIAM. Unauthorized reproduction

of this article is prohibited.

1726

GUO LUO AND THOMAS Y. HOU

where

/r, ω

/r,

and

are transformed angular velocity,

vorticity, and stream functions, respectively.

The radial and axial components of the

velocity can be recovered from

(2.1d)

−

rψ

for which the incompressibility condition

(

)

is satisfied automatically. As shown by [40], (

,ω

,ψ

) must all vanish at

=0if

a smooth velocity field. Thus (

,ω

,ψ

) are well defined as long as the corresponding

solution to (1.1) remains smooth. The reason we choose to work with the transformed

variables (

,ω

,ψ

) instead of the original variables (

,ω

,ψ

) is that the equations

satisfied by the latter,

−

)

have a formal singularity at

= 0, which is inconvenient to work with numerically.

We shall numerically solve the transformed equations (2.1) on the cylinder

(

r, z

): 0

≤

with the initial condition

(2.2a)

(

r, z

) = 100 e

−

30(1

−

)

sin

,ω

(

r, z

(

r, z

)=0

The solution is subject to a periodic boundary condition in

(2.2b)

(

r, L, t

)

,ω

(

r, L, t

)

,ψ

(

r, L, t

)

and a no-flow boundary condition on the solid boundary

=1:

(2.2c)

,z,t

)=0

The pole condition

(2.2d)

,z,t

)=0

is also enforced at the symmetry axis

= 0 to ensure the smoothness of the solution.

The initial condition (2.2a) describes a purely rotating eddy in a periodic cylinder,

and it satisfies special odd-even symmetries at the planes

L, i

Specifically,

is even at

,itisoddat

,and

,ψ

are both odd at all

’s. These symmetry properties are preserved by the equations (2.1), so instead of

solving the problem (2.1)–(2.2) on the entire cylinder

), it suffices to consider

the problem on the quarter cylinder

), with the periodic boundary condition

(2.2b) replaced by appropriate symmetry boundary conditions. It is also interesting

to notice that the boundaries of

) behave like “impermeable walls”:

(2.3)

−

rψ

=0 on

rψ

=0 on

which is a consequence of the no-flow boundary condition (2.2c) and the odd symmetry

These variables should not be confused with the components of the velocity, vorticity, and

stream function vectors.

Downloaded 01/29/15 to 131.215.70.231. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php

©

by SIAM. Unauthorized reproduction

of this article is prohibited.

FINITE-TIME SINGULARITY OF 3D EULER

1727

3. Outline of the numerical method.

The potential formation of a finite-time

singularity from the initial condition (2.2a) makes the numerical solution of the initial-

boundary value problem (2.1)–(2.2) a challenging and difficult task. In this section,

we describe a special mesh adaptation strategy (section 3.1) and a B-spline based

Galerkin Poisson solver (section 3.2), which are essential to the accurate computation

of the nearly singular solutions. The overall algorithm is outlined in section 3.3.

3.1. The adaptive (moving) mesh algorithm.

Singularities (blowups) are

abundant in mathematical models of physical systems. Examples include the semi-

linear parabolic equations describing the

blowup of the temperature of a reacting

medium, such as a burning gas [23]; the nonlinear Schr ̈

odinger equations describing

the self-focusing of electro

magnetic beams in a nonlinear medium [42]; and the ag-

gregation equations describing the concentration of interacting particles [36]. Often,

singularities occur on increasingly small length and time scales, which necessarily re-

quires some form of mesh adaptation. Further, finite-time singularities usually evolve

in a “self-similar” manner when singularity time is approached. An adaptive mesh de-

signed for singularity detection must corr

ectly capture this behav

ior in the numerical

solution.

Several methods have been proposed to compute (self-similar) singularities. In

[42], a dynamic rescaling algorithm is used to solve the cubic Schr ̈

odinger equation.

The main advantage of the method is that the rescaled equation is nonsingular and

the rescaled variable is uniformly bounded in appropriate norms. The disadvantage

is that the fixed-sized mesh is spread apart by rescaling, so accuracy is inevitably lost

far from the singularity.

In [4], a rescaling algorithm is proposed for the numerical solution of the semilinear

heat equation, based on the idea of adaptive mesh refinement. The method repeatedly

refines the mesh in the “inner” region of the singularity and rescales the inner solution

so that it remains uniformly bounded. The main advantage of the method is that it

achieves uniform accuracy across the entire computational domain and is applicable

to more general problems. The disadvantage is that it requires a priori knowledge of

the singularity and is not easily adaptable to elliptic equations (especially in multiple

space dimensions) due to the use of irregular mesh.

The moving mesh method [35] provides a very general framework for mesh adap-

tation and has been applied in various contexts, for example, the semilinear heat

equation [7] and the nonlinear Schr ̈

odinger equation [6]. The main idea of the method

is to construct the mesh based on a certain equidistribution principle, for example,

the equipartition of the arc length function. In one dimension this completely de-

termines the mesh, while in higher dimensions additional constraints are needed to

specify mesh shapes and orientations. The m

eshes are automatically evolved with the

underlying solution, typically by solving a moving mesh partial differential equation

(MMPDE).

While being very general, the “conventional” moving mesh method has the follow-

ing issues when applied to singularity detection. First, it requires explicit knowledge

of the singularity, for example, its scaling exponent, in order to correctly capture the

singularity [34]. Second, it tends to place too many mesh points near the singularity

while leaving too few elsewhere, which can cause instability. Third, mesh smoothing,

an operation necessary for maintaining stability, can significantly limit the maximum

resolution power of the mesh. Finally, the moving mesh method computes only a

dis-

crete approximation

of the mesh mapping function, which can result in catastrophic

loss of accuracy in the computation of a singularity (see section 3.3).

Downloaded 01/29/15 to 131.215.70.231. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php

©

by SIAM. Unauthorized reproduction

of this article is prohibited.

1728

GUO LUO AND THOMAS Y. HOU

For the particular blowup candidate considered in this paper, preliminary uniform

mesh computations suggest that the vorticity function tends to concentrate at a single

point. In addition, the solution appears to remain slowly varying and smooth outside

a small neighborhood of the singularity. These observations motivate the following

special mesh adaptation strategy.

The adaptive mesh covering the computational domain

) is constructed

from a pair of analytic mesh mapping functions

(

)

(

)

which are defined on [0

1], are infinitely differentiable, and have a density that is even

at both 0 and 1. The even symmetries of the mesh density ensure that the resulting

mesh can be extended

smoothly

to the full cylinder

). The mesh mapping

functions contain a small number of parameters, which are dynamically adjusted so

that a certain fraction of the mesh points (e.g., 50% along each dimension) is placed

in a small neighborhood of the singularity. Once the mesh mapping functions are

constructed, the computational domain

) is covered with a tensor-product

mesh:

(

): 0

≤

where

(

)

(

)

/N, h

/M.

The precise definition and construction of the mesh mapping functions are detailed

in Appendix A.

The mesh is evolved using the following

procedure. Starting from a reference

time

, the “singularity region”

is identified as the smallest rectangle in the

-plane that encloses the set

(

):=

(

r, z

)

∈

(

r, z, t

)

|≥

(

)

∞

,δ

∈

Once

is determined, an adaptive mesh

is fit to

and the solution is advanced

in the

ρη

-space by one time step to

. The singularity region

is then

computed and compared with

. If the ratios between the sides of

and

(in either

dimension) drop below a certain threshold (e.g., 80%), which indicates the support of

the maximum vorticity has shrunk by a sufficient amount, or if the maximum vorticity

is “too close” to the boundaries of

(3.1)

max

(

r,z

)

∈

∂S

(

r, z, t

)

|≥

(

)

∞

,δ

∈

(

which indicates the maximum vorticity is about to leave

, then a new mesh

is computed and adapted to

. In the event of a mesh update, the solution is

interpolated from

in the

ρη

-space using an 8th-order piecewise polynomial

interpolation in

and a spectral interpolation in

. The whole procedure is then

repeated with

replaced by

and

replaced by

We remark that the mesh update criterion (3.1) is designed to prevent the peak

vorticity from escaping the singularity region, as is the case in one of our earlier

computations where the singularity keeps moving toward the symmetry axis. Since

Downloaded 01/29/15 to 131.215.70.231. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php

©

by SIAM. Unauthorized reproduction

of this article is prohibited.

FINITE-TIME SINGULARITY OF 3D EULER

1729

in the current computation the singularity is fixed at the corner ̃

=(1

,the

criterion (3.1) has p

ractically no effect.

The mesh adaptation strategy described above has several advantages compared

with the conventional moving mesh method. First, it can automatically resolve a

self-similar singularity regardless of its scalings, provided that the singularity has a

bell-shaped similarity profile, which is what we observe in our case (see Figure 1(b)).

This is crucial to the success of our comput

ations, because the (a

xisymmetric) Euler

equations allow for infinitely many self-similar scalings (see section 4.7), which means

that the scaling exponent of the singularity cannot be determined a priori. Second,

the method always places enough mesh points (roughly 50% along each dimension)

outside the singularity region, ensuring a well-behaved and stable mesh (see section

4.1). Third, the explicit control of the mesh mapping functions eliminates the need

of mesh smoothing, which allows the mesh to achieve arbitrarily high resolutions.

Finally, the analytic representation of the mesh mapping functions ensures accurate

approximations of space derivatives, hen

ce greatly improving the quality of the com-

puted solutions (see section 3.3).

3.2. The B-spline based Galerkin Poisson solver.

A key observation we

made from our computations is that the overall accuracy of the numerical solution of

the initial-boundary value problem (2.1)–(2.2) depends crucially on the accuracy of the

Poisson solver. Among the methods commonly used for solving Poisson equations,

namely finite difference, finite element Galer

kin, and finite element collocation, we

have chosen the Galerkin method for both its high accuracy and for its rigorous

theoretic framework, which makes the error analysis much easier.

We have designed and implemented a B-spline based Galerkin method for the

Poisson equation (2.1c). Compared with the “conventional” Galerkin methods based

on piecewise polynomials, the B-spline based

method requires no mesh generation and

hence is much easier to implement. More importantly, the method can achieve

arbi-

trary global

smoothness and approximation order

with relative ease and few degrees

of freedom, in contrast to the conventional

piecewise polynomial based methods. This

makes the method a natural choice for our problem.

The Poisson equation (2.1c) is solved in the

ρη

-space using the following proce-

dure. First, the equation is recast in the

ρη

-coordinates:

−

ω,

(

ρ, η

)

∈

where for clarity we have written

for

and

for

. Next, the equation is multi-

plied by

and is integrated over the domain [0

where

∈

(to be defined

below) is a suitable test function. After a routine integration by parts, this yields the

desired weak formulation of (2.1c), which reads as follows: find

∈

such that

(

ψ,φ

):=

dρ dη

ωφr

dρ dη

(

)

∀

∈

(3.2a)

where (recall th

e odd symmetry of

=span

∈

(

−

ρ, η

(

ρ, η

)

,η

)=0

,φ

(

ρ,

−

(

ρ,

)

∀

∈

Downloaded 01/29/15 to 131.215.70.231. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php

©

by SIAM. Unauthorized reproduction

of this article is prohibited.

1730

GUO LUO AND THOMAS Y. HOU

To introduce Galerkin approximation, we define the finite-dimensional subspace

weighted uniform B-splines

[30] of even order

w,h

=span

(

)

j,h

(

)

i,h

(

)

∩

where

(

) is a nonnegative weight function of order 1 vanishing on

=1,

(

)

∼

−

)

,ρ

→

−

and

(

((

s/h

)

−

(

−

2)) are shifted and rescaled uniform B-splines of order

. The Galerkin formulation then reads as follows: find

∈

such that

(3.2b)

(

,φ

(

)

∀

∈

With suitably chosen basis functions of

, this gives rise to a symmetric, positive

definite linear system

which can be solved to yield the Galerkin solution

The detailed construction of the linear system is given in Appendix B.

The parameters used in our computations are

=6and

(

)=1

−

Using the theory of quasi interpolants, it can be shown that

(3.3)

|∇

−∇

dr dz

≤

(

)

−

|≤

−

∂

∇

dr dz,

where

∇

∂

,∂

)

∂

are differential operators in

-and

ρη

-planes,

respectively,

is a mesh-mapping dependent constant, and

is an absolute con-

stant. In our computations, the constant

is observed to be very close to 1 for all

times, which confirms the stability of the Galerkin solver.

The detailed error analysis of the Poisson solver will be reported in a separate

paper.

3.3. The overall algorithm.

Givenanadaptivemesh

and the data (

,ω

)

defined on

, the solution is advanced using the following procedure. First, the Pois-

son equation (2.1c) is solved for

in the

ρη

-space using a 6th-order B-spline based

Galerkin method (section 3.2). Second, the 2D velocity ̃

)

is evaluated at

the grid points using (2.1d). Third, an adaptive time step

is computed on

that the CFL condition is satisfied with a suitably small CFL number

(e.g., 0.5),

and the relative growth of the solution in o

ne step does not exceed a small threshold

(e.g., 5%). Finally, the solution (

,ω

) is advanced according to (2.1a)–(2.1b) by

using an explicit 4th-order Runge–Kutta method, and the mesh

is adapted to

the new solution if necessary (section 3.1).

In the last step of the algorithm, the evolution equations for

and

are

semidiscretized in the

ρη

-space, where the space derivatives are expressed in the

ρη

-coordinates and are approximated using 6

th-order centered difference formulas,

e.g.,

(

)=:(

)

(

)

(

)

≈

(

)

(

ρ,

)

Here, as usual,

ρ,

−

ρ,

−

ρ,

−

Downloaded 01/29/15 to 131.215.70.231. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php

©

by SIAM. Unauthorized reproduction

of this article is prohibited.

FINITE-TIME SINGULARITY OF 3D EULER

1731

denotes the standard 6th-order centered approximation to

∂

,and

(

ρ,

)

(

i,j

−

i,j

)

(

ρ,

)

(

i,j

−

i,j

−

)

denote the standard forward, backward, a

nd centered difference operators, respec-

tively. Note that the derivative

of the mesh mapping function is computed directly

from the analytic representation of

without any difference approximation. This is

crucial for the accurate evaluation of

, especially in “singularity regions” where the

inverse mesh density

is close to 0 and is nearly constant (Appendix A; in par-

ticular, see (A.3)). When

is small and nearly constant, a high-order difference

approximation of

tends to be contaminated by catastrophic cancellation, and the

discretely approximated values of

can have large relative errors or even become

negative

, causing failures of the entire computation. By computing

directly from

the analytic representation of

, this problem is avoided and the solution is ensured to

be accurately approximated even in regions where the singularity is about to form and

where

≈

1. This also explains why the conventional moving mesh method is

not suitable for singularity computations where high accuracy is demanded, because

the method computes only a

discrete approximation

of the mesh mapping function,

which necessarily requires a difference approximation of

in the evaluation of a space

derivative

. Without mesh smoothing, this can cause instability, while with mesh

smoothing the mesh resolution will be inevitably limited, which is undesired.

The centered difference formulas described above need to be supplemented by

numerical boundary conditions near

ρ, η

1. Along the

-dimension, the symmetry

condition

−

i,j

−

i,j

−

i,j

≤

is used near

=0and

= 1, where the + sign applies to

and the

−

sign applies

. Along the

-dimension, the symmetry condition

−

i,j

≤

is used near the axis

= 0 and the extrapolation condition

(

ρ,

−

)

≤

is applied near the solid boundary

=1.

The extrapolation condition is known to be

GKS stable for linear hyperbolic problems [28, Theorem 13.1.3], and it is expected to

remain stable when applied to the Euler equations as long as the underlying solution

is sufficiently smooth.

Using the superconvergence properties of the Poisson solver at the grid points (to

be proved elsewhere), it can be shown that the overall algorithm is formally 6th-order

accurate in space and 4th-order accurate in time. The details of this error analysis

will be reported in a separate paper.

4. Numerical results.

We have numerically solved the initial-boundary value

problem (2.1)–(2.2) on the quarter cylinder

)(with

). The results

suggest that the solution develops a singularity in finite time, and we shall provide,

While a 6th-order extrapolation condition (

ρ,

−

)

= 0 suffices to maintain a formal 6th-

order accuracy for the overall scheme, we choose the higher-order extrapolation condition for better

accuracy.

Downloaded 01/29/15 to 131.215.70.231. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php