A Test Problem for Flow Codes
Henry von Wahl L. Ridgway Scott






A Test Problem for Flow Codes
    
Henry von Wahl, Friedrich-Schiller-Universität Jena 

L. Ridgway Scott, University of Chicago


    April 29, 2024
====================================================================================================

arabic



We propose a test problem for Navier–Stokes solvers based on the flow around
a cylinder. We choose a range of Reynolds numbers for which the flow is
time-dependent but can be characterized as essentially two-dimensional. The
test problem requires accurate resolution of chaotic dynamics over a long time
interval. It also requires the use of a relatively large computational domain,
part of which is curved, and it requires evaluation of derivatives of the
solution and pressure on the curved boundary. We review the performance of
different finite element methods for the proposed range of Reynolds numbers.
These tests indicate that some of the most established methods do not capture
the correct behavior.


We present a classic fluids problem, flow around a cylinder, as a challenge
for numerical schemes. As reviewed in <cit.>, there is ample
experimental data with which to compare simulations. To make the challenge more
tractable, we focus on 2D simulations and lower Reynolds numbers. We list some
references where the computations are at variance with existing information.

Flow around a cylinder is not only a fundamental problem in fluid dynamics but
also a practical problem of interest in energy generation. So-called bladeless
turbines <cit.> have been
proposed as a viable energy generation method. The design of such systems
requires accurate simulations of the flow around cylinders.

We review some published examples where erroneous results were obtained. We
also show that some commonly used techniques fail to give acceptable results
at Reynolds numbers of interest. As computational simulation has emerged as a
new form of experimentation in physics, it is essential to have some benchmarks
that give guidance, which we use here. For many flow problems, there is no
dispute regarding the Navier–Stokes equations as a suitable model, as we do
here. However, this requires reliable methods to decide fundamental physics
questions, such as whether the resulting flow is chaotic at a given Reynolds number.

Flow around a cylinder has been proposed as a test problem
before <cit.>. What is
different here is that we focus on a range of Reynolds numbers for which the
flow appears to be chaotic <cit.>. Thus, different metrics must be
used to compare different numerical methods. Moreover, the time interval for
the simulation is quite long, so this simulation evaluates the ability of
different discretizations with regard to accumulation of numerical error and
challenges the implementations with regard to efficiency.

We focus on simple metrics so that different methods can be compared easily.
The lift and drag on a cylinder are two scalar quantities computed by many
methods and measured in many experiments. Another metric is the Strouhal number
or period of the time-dependent flow. In <cit.>, a simple approach was
used to estimate the periodicity of chaotic flow. We present a slightly
improved method compared to <cit.> for evaluating the Strouhal period
that provides an indication of the transition to chaos. Additional metrics,
such as the Lyapunov exponent and the fractal dimension of the attractor,
provide a finer level of detail regarding the flow dynamics <cit.>.

When the Reynolds number is more than 10^4, the flow becomes three-dimensional
<cit.>. Thus, we restrict ourselves to Reynolds numbers
in the range [10^2,10^4] for simplicity. Three-dimensional effects of a
different kind are seen <cit.> at lower
Reynolds numbers due to the finite length of cylinders used in experiments.
However, these are potentially like so-called blocking effects in
two-dimensional flow <cit.>, which become less
significant at higher Reynolds numbers.

In this work, We consider several methods that are well-studied and have a
solid theoretical base. We compare their performance to a fifth method used
in <cit.>, and we find that the simulations in <cit.> were
under-resolved for higher Reynolds numbers. Our main conclusion is that
high-order methods are necessary to obtain reliable results at higher Reynolds
numbers. Surprisingly, the commonly used, lowest-order Taylor–Hood method is
not accurate enough to be used in such settings. It appears that methods that
result in exactly divergence-free velocity solutions have a significant
advantage at higher Reynolds numbers, even though the one we investigate is
non-conforming.

One of the main conclusions from <cit.> was that the simulation drag
values differ substantially from experiments, as shown in Figure <ref>.
We confirm this conclusion here, suggesting that cylinder vibration likely
plays a significant role in practice <cit.>.






§ SETTING THE PROBLEM AND MODEL EQUATIONS


Suppose that (,p) is a solution of the time-dependent Navier–Stokes
equations in a domain Ω⊂^d containing an obstacle with
boundary Γ⊂. This fulfills the equations

    ∂_t-νΔ    + ·∇+∇ p =  in Ω,
        =0 in Ω,

with the kinematic viscosity ν, and together with boundary conditions

    = on ∂Ω\Γ,  = on Γ.

For the well-posedness of these equations, see <cit.>.



 §.§ Problem Set-up

For our simulations, we consider the domain

    Ω=(x,y) : -30<x<300, |y|<30,   x^2+y^2>1.

The boundary condition is set as = (1, 0)^T and the cylinder diameter
as the reference length, the Reynolds number is given by Re = 2 / ν. We
consider the time interval [0, 500].





 §.§ Weak formulation of the Navier-Stokes equations

While this test problem is not limited to finite-element
simulations, we expect many will be done with finite elements. To this end,
we require a variational formulation of (<ref>).
The Navier-Stokes equations can be written in a weak (or variational)
form as follows:

Find ∈ (H^1_g(Ω))^d and p ∈ L^2_0(Ω) such that

    (∂_t,)_L^2(Ω)+ a(, ) + b(, p) +c(,,)   = 0
      ∀∈ (H^1_0(Ω))^d, 
    
    b(, q)    = 0   ∀ q ∈ L^2_0(Ω).

The space H^1_g(Ω) is the space of H^1 functions with trace g on the
(outer) boundary, H^1_0(Ω) is the space of H^1 functions with trace
0 on the boundary and L^2_0(Ω) is the space of L^2 functions with
mean zero. The bilinear forms a(·, ·) and b(·, ·) are
defined as

    a(, ) = ∫_Ων/2𝒟():𝒟(),
         b(, q) = -∫_Ω q (∇·) ,

respectively, where =∇+∇^t and the colon (:) indicates the
Frobenious inner-product of matrices. The convective (nonlinear) term takes various forms.
By default, we take it to be

    c(,,)=c_conv(, , ) = ∫_Ω (·∇)· .

For divergence-free functions, the various forms of the convective term
are equivalent. Hence, we may take the above so-called convective form
of the trilinear form for exactly divergence-free methods. However, for
non-divergence-free methods, other forms are beneficial, see below.



 §.§ Drag and lift evaluation


The flow of fluid around an obstacle generates a force called drag, which is a
fundamental concept in fluid dynamics <cit.>.
It plays a critical role in determining the behavior of objects in flight and
has been studied since the time of d'Alembert <cit.>.
Drag is composed of two components: pressure drag β_p and
viscous drag β_v, which can be calculated by evaluating the
following functions using =(1,0):

    β_v()    =∮_((ν))· s,     β_p()  =∮_  -p· s.

The full drag β is defined by β=β_v+β_p.
Similarly, lift is computed using =(0,1).




  §.§.§ Alternate drag evaluation


Another way <cit.> to evaluate β is
to test the weak formulation of (<ref>) (see (<ref>)
below) with a non-conforming test function. This results in the functional

    ω()   =∫_Ω∂_t +
       ∫_Ων/2: +(·∇)·-p,

Then ω()=β() for all ∈ H^1(Ω)^d
<cit.>.
This form is potentially more accurate, but it requires an accurate
approximation of the time derivative ∂_t.

This approach may be traced back to <cit.>. It is also known,
that for the case of strongly imposed Dirichlet boundary conditions, this
approach of testing the residual with a non-conforming test-function doubles
the rate of convergence for the drag and lift, see, for
example <cit.> for the proof in the steady case.

Having two ways to compute drag and lift provides a valuable internal
check of the accuracy of the simulation.



§ CHALLENGE COMPUTATION


The data in Figure <ref> suggests a significant disagreement
between computation and experiment.
This could be for a variety of reasons.
The most obvious is some flaw in the computational scheme.
Thus, we propose this problem as a computational challenge to see
what other methods predict.

Another reason for the disagreement could be that the model is wrong.
Experimentalists have raised the issue that the cylinders used could vibrate.
In particular, the Relf data is based on flows past thin wires,
used in musical instruments for their vibrational qualities.
Thus, one response to the challenge is to allow the cylinder to vibrate
back and forth in response to the oscillating lift and drag.
The question would be to evaluate how oscillations of the cylinder change
parameters describing the flow, such as drag, lift, and Strouhal number,
representing the frequency of the generated vortices.
One study with forced vibrations <cit.>
suggests that the change could be substantial.




§ DISCRETIZATION METHODS



The simulations in this paper are done using finite element methods. We,
therefore, begin py presenting some details of the finite element approach.

To discretize the weak formulation (<ref>) with finite element
methods, we take a mesh of the domain Ω with characteristic length
h, denoted as . Finite element spaces then look to discretize
(<ref>) by approximating the spaces H^1(Ω) and
L^2(Ω) with spaces of piecewise polynomials on each element of the
mesh, with varying degree of discontinuity of these polynomials over
element edges.

We utilize two well-established simulation methods for which extensive
numerical analysis is available, establishing not only convergence but
also detailed simulation properties.



 §.§ H(div) conforming finite element discretization



We consider an -conforming finite element method based on
<cit.>. This discretization
is based on mass-conserving mixed formulation with symmetric stresses, also
known as the MCS formulation. This approach discretizes the velocity in 
using Brezzi–Douglas–Marini (BDM) elements <cit.>. That is, the
normal component of the velocity is continuous across element boundaries, while
the tangential component is discontinuous. The pressure is discretized using
discontinuous elements. This method is exactly mass conserving and
pressure-robust.
From a numerical point of view, pressure-robust methods have several advantages.
First, the velocity remains unchanged by irrotational changes in any external
forcing term, a property that most methods do
not preserve on the discrete level <cit.>.
More importantly, pressure-robust methods have the property that the velocity
error is independent of the pressure error <cit.>.
This contrasts classical Galerkin discretizations of the
(Navier)-Stokes equations, where the velocity error depends on
1/ν times the pressure error. From a physical point of view,
exactly divergence-free methods have the advantage that, in addition to mass,
kinetic energy, linear momentum, and angular momentum are conserved on the
discrete level <cit.>.

Furthermore, this method is based on a formulation where an additional variable
σ = ν∇ which is discretised in the function space .
Thus the energy form is modified to

    a'(σ, τ) = 1/ν∫_Ωσ:τ,

and a new form is defined:

    b'(τ,) = ∑_T∈∫_T(τ)·
    +∑ _e∈∫_e [τ_nn]· s,

where  is the triangular mesh, T denotes triangles in , and
e denotes edges in .
The expression [τ_nn] denotes the jump in the normal-normal component
of the tensor τ across the edge e, and  is the normal to the edge.
Note that the orientation of the normal is not material.
Thus the formulation is

    a'(σ_h, τ_h) + b'(τ_h,_h)   = 0    ∀τ_h ∈Σ_h, 
    
    (∂_t_h,_h)_L^2(Ω)+ b'(σ_h,_h) + c(_h,_h,_h)
         + b(_h, p_h)    = 0   ∀_h ∈ V_h,
    
    b(_h, q_h)    = 0   ∀ q_h ∈ Q_h.

Here, Σ_h is a space of discontinuous piecewise polynomials of degree k,
whose values are trace-zero tensors satisfying continuity of the
`normal-tangential component' across edges and that the normal-tangential
component is a polynomial of order k-1 on each edge,
V_h denotes the BDMk finite-element space of order k,
and Q_h consists of discontinuous piecewise polynomials of degree k-1.

As  spaces are characterized by normal continuity across
facets, boundary condition in the normal components are imposed strongly.
The implementation used uses hybridisation for the tangential component of the
velocity. That is, a second finite element space is used to couple the
tangential component across element facets and the tangential part of
the boundary condition is enforced strongly on the facet space.

For time-stepping we use a second order IMEX scheme known as SBDF2 <cit.>.
This scheme uses the BDF2 scheme to discretise the time-derivative, the Stokes
term is fully implicit, and the non-linear convective term is treated explicitly
with a second order extrapolation:


    a'(σ_h^n, τ_h) + b'(τ_h,_h^n)   = 0    ∀τ_h ∈Σ_h, 
    [b]
    1/Δ t(_h^n - 2 _h^n-1 + _h^n-2,_h)_L^2(Ω)+ b'(σ_h^n,_h) + b(_h, p_h^n)
    
    + 2 c(_h^n-1,_h^n-1,_h) - c(_h^n-2,_h^n-2,_h)
       = 0   ∀_h ∈ V_h,
    
    b(_h^n, q_h)    = 0   ∀ q_h ∈ Q_h.


To regularize the saddle-point system and to fix the pressure constant in the
case of pure Dirichlet boundary conditions, we may add a small perturbation to
the lower right pressure block of the block-system resulting from
(<ref>). That is, we modify
(<ref>) to

    -∫_Ω q_h (∇·_h^n) + ϵ∫_Ω p_h^n q_v
       = -∫_Ω q_h ((∇·_h^n) - ϵ p_h^n)
       =0,   ∀ q_h∈ Q_h,

with some ϵ > 0. Since our discrete spaces V_h, Q_h satisfy the
discrete de-Rahm complex <cit.>, we have the property that
∇· V_h = Q_h. Consequently, we can replace q_h with
∇·_h in (<ref>). This gives us an
equation for the pressure, which we may plug into the second equation in
(<ref>). This leads to the discrete system,

    a'(σ_h^n, τ_h) + b'(τ_h,_h^n)
       + 1/Δ t(_h^n,_h)_L^2(Ω)
       + b'(σ_h^n,_h) - b(_h, 1/ϵ∇·_h).

This is a smaller linear system to solve, and the system is no longer a saddle
point problem, making it an easier system to solve numerically. We want to
choose ϵ small to keep the perturbation from the divergence constraint
small. However, this means that the 1/ϵ-scaled term causes the
linear system (<ref>) to be very ill-conditioned.
Consequently, we need to use direct linear solvers to handle these
ill-conditioned linear systems, which in turn limits us with respect to the
maximal size of the system we can solve and the parallel scalability.
With the small perturbation, the divergence of _h is no longer exactly zero,
but _h^n_L^2≤ϵp_h^n_L^2.
In the computations, we chose ϵ=10^-12 / ν.


We will refer to approach above as _k, where k is the order of the velocity
finite element space. This approach is implemented using
Netgen/NGSolve <cit.> and based on a template
available on github[<https://github.com/NGSolve/modeltemplates>].



 §.§ Taylor-Hood

To compare with more established methods, we also consider Taylor-Hood elements
consisting of H^1-conforming finite elements of order k≥ 2 for the
velocity and order k-1 for the pressure.

To stabilize this scheme, we optionally add grad-div stabilisation to improve
the conservation of mass<cit.>. That is, we add the
term

    j_h(_h,_h) γ_gd∫_Ω∇·_h ∇·_h,

to our bilinear form. In our computations, we take γ_gd=10^3, and we
refer to this method as _k, where k is the order of the velocity space.
In case no grad-div stabilization is used, we refer to the method as _k.
Furthermore, we note that it is known that for γ_gd→∞,
the grad-div stabilized Taylor-Hood solution converges to the Scott-Vogelius
solution <cit.>.

For the convective term, there are multiple choices possible, we take the
divergence form <cit.>,

    c_div(,,) = c_conv(,,) 
          + 1/2∫_Ω(∇·)·,

which is equivalent to the standard convective form on the continuous level,
but vanishes under less restrictive conditions for = on the discrete
level.

Thus the variational formulation takes the form

    (∂_t_h,_h)_L^2(Ω)+ a(_h, _h) +j_h(_h,_h)
    + b(_h,p_h) +c_div(_h,_h,_h)   = 0   ∀_h ∈ V_h, 
    
    b(_h, q_h)    = 0   ∀ q_h ∈ Q_h,

where V_h is the standard Lagrange space of continuous,
piecewise polynomials of degree k≥ 2, and
Q_h is the standard Lagrange space of continuous,
piecewise polynomials of degree k-1.

To discretise in time-we use the fully implicit Crank-Nicolson scheme, that is

    1/Δ t(_h^n-_h^n-1,_h)_L^2(Ω)
       +(a(^n_h,_h) + c_div(^n_h,u^n_h,_h))
       +j_h(^n_h,_h) + b(_h,p^n_h)
    
          + (a(^n-1_h, _h) +c_div(^n-1_h,^n-1_h,_h)) 
                = 0    ∀_h∈ V_h,
    
          b(^n_h, q_h)    = 0    ∀ q_h ∈ Q_h,

Note that we do not take into account the pressure from the previous time step,
nor do we add the grad-div stabilization term from the last time step. This is
because the method needs to correctly reflect that the is Langrange-multiplier
to weakly enforce the divergence constraint <cit.>.
Similarly, the grad-div form is also there to penalize deviations from the
divergence constraint.

We use a quasi-Newton method with line search to solve the nonlinear systems
arising in every time step. More precisely, we only recompute the factorization
of the inverse Jacobi-matrix if the residual decreases by less than a
factor of 0.3. The line search consists of a check that the Newton update
has decreased the residual. Otherwise, we reduce the size of the update by half
recursively until the residual has decreased compared to the previous iteration.
We start each quasi-Newton iteration with (u_h,p_h) from the previous time step.
As the scheme is fully implicit,
we can use time steps independent of the mesh size, and as the Jacobi-matrix
rarely has to be refactorized, the method is essentially as efficient as
IMEX schemes.

The Taylor–Hood methods are also implemented using Netgen/NGSolve. The meshes
used are constructed as those for the -conforming method in
Section <ref>.




 §.§ Scott–Vogelius


As a final method, we also consider the Scott–Vogeluis finite element pair.
This finite element pair is very similar to the Taylor–Hood pair, using
H^1-conforming finite elements for the velocity but discontinuous elements
of order k-1 for the pressure. Consequently, we use the same variational
formulation time-stepping scheme used for the unstabilized Taylor-Hood method.
For k≥ 4 in two dimensions, this method is inf-sup stable on meshes
without nearly singular vertices <cit.> and for k≥ d=2 on
barycentric refined meshes <cit.>. This method is also exactly
divergence-free and pressure-robust.

There are different ways that the Scott–Vogelius method <cit.>
can be implemented.
In <cit.>, the iterated penalty method was used to avoid explicit
knowledge of the pressure space. However, this approach was plagued by nearly
singular vertices <cit.> in the meshes used.
Here, we utilize the Netgen mesh generator, which, for the cylinder problem,
does not have any singular or nearly singular vertices.
Thus, the mixed-method implementation of Scott–Vogelius can be used, and in
this case, there is only one change from the (unstabilized) Taylor–Hood method,
namely, the pressure elements are discontinuous instead of continuous.
One drawback of this approach is the linear systems are larger, and it
appears that for Reynolds number 8000 (Table <ref>),
it is no longer possible to keep the divergence very small.
Upon closer examination, it was found that the Newton step was not converging
at a few time steps, and the divergence also spiked there.
In general, the divergence was very small except at these few times.
Furthermore, it is worth mentioning that, unlike the unstabilized Taylor-Hood
method, the lack of convergence of Newton's method did not lead to a blow-up
of the solution.




 §.§ Meshes


We used several meshes for our experiments. For the methods described in
Section <ref>, the meshes are unstructured simplicial
meshes constructed using Netgen <cit.> using a bulk mesh
parameter h_max and a local mesh size h_min=h_max/250
on the boundary of the cylinder. The part of the coarsest mesh
(h_max=8) around the cylinder is visualized in Figure <ref>.

For h_max=8,4,2, the Taylor-Hood finite element space has a total of
6994, 163107, and 437439 degrees of freedom (dofs), respectively, for the
saddle point system. The Scott-Vogelius saddle point systems have 87858, 205176,
and 551434 dofs, respectively, due to the larger pressure space. On the same
meshes, the MSC method has 95856, 224487, and 604047 dofs for the stress and
velocity space in which we solve the problem. However, we emphasize that a
direct comparison of dofs is not appropriate, as different types of systems
must be solved for each method, and the methods have significant differences
regarding their approximation power.



 §.§ Computing the period


Following <cit.>, we estimate the period of the drag-lift trajectory 
as follows. For selected points along the trajectory (about 100 per period), we 
consider points from roughly one-half period away to three-halfs away and 
determine the point of closest approach. We start with an estimated average 
period, and for each initial point we, obtain an estimated period. Histograms 
of such estimates are given in Figure <ref>. Taking the average 
period from this data gives an estimate of the Strouhal period that can be used 
to refine the process, providing a new estimated average period. To allow for 
estimates of the period at a finer grain than the individual time steps, we 
used a quadratic interpolation of estimated periods to allow for prediction of 
the period to second order. That is, we found the trajectory error at the 
optimal point and then considered the error at the time steps on either side of 
this. Fitting this data with a quadratic and finding the minimum of the 
quadratic provided an estimated period to within a fraction of Δ t.




§ COMPUTATIONAL EXPERIMENTS


We consider Reynolds numbers between 120 and 8000. The results for Reynolds
numbers 120, 250, 500, 1000, 2000, 4000, and 8000 are shown in
Table <ref>, Table <ref>,
Table <ref>, Table <ref>,
Table <ref>, Table <ref>, and
Table <ref>, respectively. Furthermore, results for
Reynolds numbers between 1100 and 1800 are shown in Table <ref>.
Finally, the vorticity of the resulting flow computed on the finest mesh
with the MSC method is shown for Reynolds numbers 500, 1000, and 2000
in Figure <ref>, Figure <ref>,
and Figure <ref>, respectively.



 §.§ Reynolds 120–500


We began with a small Reynolds number (120) to provide a low bar for diverse 
methods. Nevertheless, Table <ref> shows that unstabilized, 
low-order (k=2) Taylor–Hood, which is commonly thought of as the workhorse 
for fluid flow, struggled to get a viable result. By contrast, all of the 
high-order methods agree to several digits in the prediction of the Strouhal 
period and drag. The stabilized Taylor–Hood for k=2 is between the 
high-order predictions and the unstabilized k=2 results. Note that at this 
Reynolds number, stabilization of high-order (k=4) Taylor–Hood does not make 
a difference, except on the size of the divergence error.

Table <ref> tells a similar story with Reynolds number 250.
The high-order methods are in agreement and the low-order methods are
substandard, in particular suggesting chaotic behavior.

For Reynolds numbers 500 and higher, we decided to drop the lowest-order
Taylor–Hood simulation data. This is due to the method's poor performance
at the lower Reynolds numbers. In some cases, the unstabilized did not even
continue to time 500, with the Newton iteration failing and the solution
blowing up. For Reynolds number 500, c.f. Table <ref>, even
the stabilized lowest-order Taylor–Hood simulation data is not viable.
For higher Reynolds numbers, we no longer report results using this method.
Finally, the vorticity of the velocity solution from the MCS method on the
finest mesh is shown in Figure <ref>. Here, we see the
periodic flow in the wake of the cylinder for about 70 spatial units, after
which the flow transitions away into more chaotic behavior. However, since our
mesh is very coarse in this part of the domain, compared to the resolution
around the cylinder, we do not claim that this is the fully resolved behavior
so far downstream in the channel.













 §.§ Reynolds 1000 and higher


For Reynolds numbers 1000 and higher, only three methods remain viable: the
pressure-robust MSC and Scott–Vogelius methods, and the grad-div stabilized
high-order Taylor–Hood method. The results are presented in
Table <ref> for Reynolds number 1000. Additionally, the
vorticity of the MSC method on the finest mesh at time t=300 is shown in
Figure <ref>.

At Reynolds number 1000, the drag-lift dynamics are very close to periodic but
have a very small standard deviation in the computed period. In particular, we
see in Figure <ref> that the solution remains periodic for
about 40 units downstream of the cylinder before the solution transitions into
chaos. However, we again cannot claim high accuracy of the solution in this
part of the domain, compared to close to the cylinder. Furthermore, the
_4 and _4 drag-lift dynamics are only periodic on the finest
mesh considered. Nevertheless, the Scott–Vogelius method on the second mesh
presents a significantly smaller standard deviation of the Strouhal period
compared to the solution of the stabilized Taylor–Hood method. Overall, it is
harder to get a close agreement for the average drag and even harder for the
Strouhal period at this and higher Reynolds numbers.

In Table <ref>, we see for Reynolds number 2000 that the
drag-lift dynamics is no longer periodic, as also observed in
<cit.>. Similar
to <cit.>, we observe an oscillation of the
wake of the flow behind the cylinder in Figure <ref>.

As the drag-lift dynamics is no longer periodic at Reynolds number 2000, we
computed the flow for a sequence of Reynolds number values between 1000 and
2000 to see how the chaos emerges, as shown in Table <ref>.
We obtain a plausibly periodic solution at Reynolds number 1125, but with an
ever-increasing standard deviation of the Strouhal period as the Reynolds
number increases, even on the finest meshes.

For Reynolds numbers 4000 and higher, we no longer attempt to provide consensus
values for the Strouhal period or drag; see Table <ref> and
Table <ref>. Furthermore, Figure <ref>
illustrates how the spread in the computed Strouhal period increases with
increasing Reynolds number.



  §.§.§ Better to be div-zero than conforming?


As far as we are aware, there are no direct comparisons between conforming and
non-conforming pressure-robust methods, so our results give some data on this
issue. A direct comparison between fully discontinuous Galerkin and a
-conforming method was performed in <cit.>. Our results
indicate that the -conforming MCS method is slightly better than stabilized
Taylor–Hood in many cases. This suggests that there is an advantage to being
exactly divergence-zero, with little penalty due to the nonconformity of the
velocity. This may be because the MCS method, like Scott–Vogelius, is pressure
robust, i.e., the velocity error is independent of the pressure error. On the
other hand, for Reynolds number 4000 (Table <ref>), we see
that the MCS predictions have a significant variation with regard to mesh
variation. The Scott–Vogelius predictions tend to be smaller than
Taylor–Hood. For Reynolds number 8000 (Table <ref>), we
see that the MCS prediction for the average drag coefficient differs
significantly from the Scott–Vogelius and Taylor–Hood predictions. However,
as noted above, the Newton iteration failed to converge in a few time steps in
the Scott–Vogelius computation, although it did converge in the stabilized
Taylor–Hood simulation.

























§ OTHER METHODS


We survey briefly other methods that have been used for simulating flow around
a cylinder.



 §.§ Vortex blob


The publication <cit.> suggests that drag
and lift for flow around a cylinder is nearly periodic at R=1000, in agreement
with our findings, with chaotic flow for R≥ 2000.



 §.§ OpenFoam

OpenFoam may not be the ideal discretization approach for low Reynolds numbers.
In <cit.>, Strouhal numbers and frequencies
computed with OpenFoam were found to be in error by greater than a factor of two
for R∈[55,161].




§ CONCLUSIONS


We considered several finite element methods for solving flow past a cylinder.
All of these methods are well studied both computationally and theoretically.
We found that high-order and pressure-robust methods are necessary to obtain
reliable results at higher Reynolds numbers. The widely-used, lowest-order
Taylor–Hood method was not accurate enough to give appropriate results. We
were able to confirm some earlier results by other researchers. Exact
incompressibility appears to be important at higher Reynolds numbers, even
though the one such method that we investigated is non-conforming. We studied
in detail how chaos appears to emerge for Reynolds numbers just above 1000.



§ DATA AVAILABILITY STATEMENT


The code used to realize the presented results, the raw data generated by
this code and videos of key simulations are available on zenodo at
<https://doi.org/10.5281/zenodo.11068455>.



§ ACKNOWLEDGMENTS

This material is based upon work supported by the National Science Foundation
under Grant No. DMS-1929284 while the authors were in residence at the Institute
for Computational and Experimental Research in Mathematics in Providence, RI,
during the Numerical PDEs: Analysis, Algorithms, and Data Challenges program.

Part of this research was conducted using computational resources and services
at the Center for Computation and Visualization, Brown University.

40


urlstyle

[Ascher et al.(1995)Ascher, Ruuth, and Wetton]ARW95
U. M. Ascher, S. J. Ruuth, and B. T. R. Wetton.
Implicit-explicit methods for time-dependent partial differential
  equations.
SIAM J. Numer. Anal., 320 (3):0 797–823,
  1995.
10.1137/0732037.

[Babuška and Miller(1984)]babuvska1984post
I. Babuška and A. Miller.
The post-processing approach in the finite element method—part 1:
  calculation of displacements, stresses and other higher derivatives of the
  displacements.
Internat. J. Numer. Methods Engrg., 200 (6):0
  1085–1109, 1984.
10.1002/nme.1620200610.

[Bahadur(2022)]bahadur2022dynamic
I. Bahadur.
Dynamic modeling and investigation of a tunable vortex bladeless wind
  turbine.
Energies, 150 (18):0 6773, 2022.
10.3390/en15186773.

[Braack and Richter(2006)]braack2006solutions
M. Braack and T. Richter.
Solutions of 3D Navier–Stokes benchmark problems with adaptive
  finite elements.
Comput. Fluids, 350 (4):0 372–392, 2006.
10.1016/j.compfluid.2005.02.001.

[Case et al.(2011)Case, Ervin, Linke, and
  Rebholz]ref:LinkeRebholzSVconnection
M. A. Case, V. J. Ervin, A. Linke, and L. G. Rebholz.
A connection between Scott–Vogelius and grad-div stabilized
  Taylor–Hood FE approximations of the Navier–Stokes equations.
SIAM J. Numer. Anal., 490 (4):0 1461–1481,
  2011.

[Charnyi et al.(2017)Charnyi, Heister, Olshanskii, and Rebholz]CHOR17
S. Charnyi, T. Heister, M. A. Olshanskii, and L. G. Rebholz.
On conservation laws of Navier-Stokes Galerkin discretizations.
J. Comput. Phys., 337:0 289–308, May 2017.
10.1016/j.jcp.2017.02.039.

[de Frutos et al.(2017)de Frutos, García-Archilla, John, and
  Novo]dFGJN17
J. de Frutos, B. García-Archilla, V. John, and J. Novo.
Analysis of the grad-div stabilization for the time-dependent
  navier–stokes equations with inf-sup stable finite elements.
Adv. Comput. Math., 440 (1):0 195–225, 2017.
10.1007/s10444-017-9540-1.

[Durante et al.(2022)Durante, Pilloton, and
  Colagrossi]ref:wolfchaocylinderflow
D. Durante, C. Pilloton, and A. Colagrossi.
Intermittency patterns in the chaotic transition of the planar flow
  past a circular cylinder.
Phys. Rev. Fluids, 70 (5):0 054701, 2022.
10.1103/physrevfluids.7.054701.

[El-Shahat et al.(2018)El-Shahat, Hasan, and Wu]el2018vortex
A. El-Shahat, M.-M. Hasan, and Y. Wu.
Vortex bladeless wind generator for nano-grids.
In 2018 IEEE Global Humanitarian Technology Conference (GHTC),
  pages 1–2. IEEE, 2018.
10.1109/GHTC.2018.8601572.

[Fehn et al.(2019)Fehn, Kronbichler, Lehrenfeld, Lube, and
  Schroeder]FKLLS_IJNMF_2019
N. Fehn, M. Kronbichler, C. Lehrenfeld, G. Lube, and P. W. Schroeder.
High-order DG solvers for under-resolved turbulent incompressible
  flows: A comparison of L2 and H(div) methods.
Internat. J. Numer. Methods Fluids, 910 (11):0
  533–556, 2019.
10.1002/fld.4763.

[Francis et al.(2021)Francis, Umesh, and Shivakumar]francis2021design
S. Francis, V. Umesh, and S. Shivakumar.
Design and analysis of vortex bladeless wind turbine.
Materials Today: Proceedings, 47:0 5584–5588, 2021.
10.1016/j.matpr.2021.03.469.

[Girault and Raviart(1986)]giraultraviart
V. Girault and P.-A. Raviart.
Finite Element Methods for Navier–Stokes Equations.
Springer Series in Computational Mathematics, volume 5.
  Springer–Verlag, Berlin, 1986.
10.1007/978-3-642-61623-5.

[Gjerde and Scott(2023a)]lrsBIBjn
I. G. Gjerde and L. R. Scott.
Resolution of d'Alembert's Paradox using Navier's slip boundary
  conditions.
2023a.
10.48550/arXiv.2204.12240.

[Gjerde and Scott(2023b)]lrsBIBkj
I. G. Gjerde and L. R. Scott.
Verification and validation of cylinder drag: Pressure and stress
  approximations on curved boundaries.
2023b.
10.48550/arXiv.2306.12362.

[Gopalakrishnan et al.(2020a)Gopalakrishnan, Lederer, and
  Schöberl]gopalakrishnan2020mass
J. Gopalakrishnan, P. L. Lederer, and J. Schöberl.
A mass conserving mixed stress formulation for the Stokes
  equations.
IMA J. Numer. Anal., 400 (3):0 1838–1874,
  2020a.
10.1093/imanum/drz022.

[Gopalakrishnan et al.(2020b)Gopalakrishnan, Lederer, and
  Schöberl]gopalakrishnan2020weak
J. Gopalakrishnan, P. L. Lederer, and J. Schöberl.
A mass conserving mixed stress formulation for Stokes flow with
  weakly imposed strss symmetry.
SIAM J. Numer. Anal., 580 (1):0 706–732,
  2020b.
10.1137/19m1248960.

[Hayakawa and Hussain(1989)]ref:FazleHayakawatreedee
M. Hayakawa and F. Hussain.
Three-dimensionality of organized structures in a plane turbulent
  wake.
J. Fluid Mech., 206:0 375–404, 1989.
10.1017/S0022112089002338.

[John(2004)]ref:refvalcyliftdragVolkerJohn
V. John.
Reference values for drag and lift of a two-dimensional
  time-dependent flow around a cylinder.
Internat. J. Numer. Methods Fluids, 440 (7):0
  777–788, 2004.
10.1002/fld.679.

[John(2016)]JohnBook16
V. John.
Finite element methods for incompressible flow problems.
Springer, Cham, 2016.
10.1007/978-3-319-45750-5.

[John et al.(2017)John, Linke, Merdon, Neilan, and Rebholz]JLM17
V. John, A. Linke, C. Merdon, M. Neilan, and L. G. Rebholz.
On the divergence constraint in mixed finite element methods for
  incompressible flows.
SIAM Rev., 590 (3):0 492–544, Jan. 2017.
10.1137/15m1047696.

[Kornbleuth(2016)]ref:KornbleuthFoamCylinder
M. Kornbleuth.
Studying the Viscous Flow Around a Cylinder Using OpenFoam.
PhD thesis, Astronomy Department, Boston University, Boston, MA, USA,
  2016.

[Lamb(1932)]lamb1993hydrodynamics
S. H. Lamb.
Hydrodynamics.
Cambridge University Press, sixth edition, 1932.

[Linke(2014)]Link2_helmholtz_2014
A. Linke.
On the role of the helmholtz decomposition in mixed methods for
  incompressible flows and a new variational crime.
Comput. Methods Appl. Mech. Engrg., 268:0 782–800,
  2014.
10.1016/j.cma.2013.10.011.

[Linke et al.(2011)Linke, Rebholz, and Wilson]LINKE2011612
A. Linke, L. G. Rebholz, and N. E. Wilson.
On the convergence rate of grad-div stabilized taylor–hood to
  scott–vogelius solutions for incompressible flow problems.
J. Math. Anal. Appl., 3810 (2):0 612–626,
  2011.
https://doi.org/10.1016/j.jmaa.2011.03.019.

[Meneghini(1993)]ref:Meneghini1994Thesis
J. R. Meneghini.
Numerical simulation of bluff body flow control using a
  discrete vortex method.
PhD thesis, Imperial College London, 1993.

[Olshanskii et al.(2009)Olshanskii, Lube, Heister, and
  Löwe]OLSHANSKII20093975
M. Olshanskii, G. Lube, T. Heister, and J. Löwe.
Grad–div stabilization and subgrid pressure models for the
  incompressible navier–stokes equations.
Comput. Methods Appl. Mech. Engrg., 1980
  (49):0 3975–3988, 2009.
10.1016/j.cma.2009.09.005.

[Qin(1994)]Qin94
J. Qin.
On the convergence of some low order mixed finite elements for
  incompressible fluids.
PhD thesis, Pennsylvania State University, 1994.
URL
  <https://www-users.math.umn.edu/ arnold/papers/QinThesis.pdf>.

[Richter(2017)]Ric17
T. Richter.
Fluid-structure Interactions, volume 118 of Lecture
  Notes in Computational Science and Engineering.
Springer, Cham, 2017.
10.1007/978-3-319-63970-3.

[Sahin and Owens(2004)]ref:blockagecylinderatio
M. Sahin and R. G. Owens.
A numerical investigation of wall effects up to high blockage ratios
  on two-dimensional flow past a confined circular cylinder.
Phys. Fluids, 160 (5):0 1305–1320, 2004.
10.1063/1.1668285.

[Schäfer et al.(1996)Schäfer, Turek, Durst, Krause, and
  Rannacher]schafer1996benchmark
M. Schäfer, S. Turek, F. Durst, E. Krause, and R. Rannacher.
Benchmark computations of laminar flow around a cylinder.
In E. H. Hirschel, editor, Flow Simulation with
  High-Performance Computers II: DFG Priority Research Programme Results
  1993–1995, pages 547–566. Vieweg+Teubner Verlag, Wiesbaden, 1996.
ISBN 978-3-322-89849-4.
10.1007/978-3-322-89849-4_39.

[Schöberl(1997)]schoeberl_netgen
J. Schöberl.
NETGEN an advancing front 2D/3D-mesh generator based on
  abstract rules.
Comput. Vis. Sci., 10 (1):0 41–52, 1997.
10.1007/s007910050004.

[Schöberl(2014)]schoeberl_ngsolve
J. Schöberl.
C++11 implementation of finite elements in NGSolve.
Technical report, 2014.

[Schroeder and Lube(2017)]SL17
P. W. Schroeder and G. Lube.
Pressure-robust analysis of divergence-free and conforming FEM for
  evolutionary incompressible Navier-Stokes flows.
J. Numer. Math., 250 (4), Dec. 2017.
10.1515/jnma-2016-1101.

[Scott(2018)]lrsBIBih
L. R. Scott.
Introduction to Automated Modeling with FEniCS.
Computational Modeling Initiative, 2018.

[Scott and Durst(2024)]lrsBIBki
L. R. Scott and R. Durst.
Chaotic dynamics of two-dimensional flows around a cylinder.
Phys. Fluids, 360 (2):0 024118, 2024.
10.1063/5.0186496.

[Scott and Vogelius(1985)]SV85
L. R. Scott and M. Vogelius.
Norm estimates for a maximal right inverse of the divergence operator
  in spaces of piecewise polynomials.
190 (1):0 111–143, 1985.
10.1051/m2an/1985190101111.

[Vogelius(1983)]Vog83
M. Vogelius.
A right-inverse for the divergence operator in spaces of piecewise
  polynomials: Application to thep-version of the finite element method.
Numer. Math., 410 (1):0 19–37, Feb. 1983.
10.1007/bf01396303.

[Williamson(1989)]ref:obliqueparallelmodesWilliamson
C. H. K. Williamson.
Oblique and parallel modes of vortex shedding in the wake of a
  circular cylinder at low Reynolds numbers.
J. Fluid Mech., 206:0 579–627, 1989.
10.1017/S0022112089002429.

[Williamson and Govardhan(2004)]ref:annurevfluidvortexvibrations
C. H. K. Williamson and R. Govardhan.
Vortex-induced vibrations.
Annu. Rev. Fluid Mech., 36:0 413–455, 2004.
10.1146/annurev.fluid.36.050802.122128.

[Zhang(2004)]Zha04
S. Zhang.
A new family of stable mixed finite elements for the 3d stokes
  equations.
Math. Comp., 740 (250):0 543–555, Aug. 2004.
10.1090/s0025-5718-04-01711-9.



