Flow of viscoelastic fluids past a cylinder at
high Weissenberg number: stabilized
simulations using matrix logarithms
Martien A. Hulsen
∗
,
Materials Technology, Eindhoven University of Technology, P.O. Box 513,
5600MB, Eindhoven, The Netherlands
Raanan Fattal
School of Computer Science and Engineering, The Hebrew University, Jerusalem,
91904, Israel
Raz Kupferman
Institute of Mathematics, The Hebrew University, Jerusalem, 91904, Israel
Abstract
The log conformation representation proposed in [1] has been implemented in
a fem context using the DEVSS/DG formulation for viscoelastic fluid flow. We
present a stability analysis in 1D and attribute the high Weissenberg problem to
the failure of the numerical scheme to balance exponential growth. A slightly dif-
ferent derivation of the log based evolution equation than in [1] is also presented.
We show numerical results for the flow around a cylinder for an Oldroyd-B and
a Giesekus model. We provide evidence that the numerical instability identified in
the 1D problem is also the actual re ason for the failure of the standard fem im-
plementation of the problem. With the log conformation representation we are able
to obtain solutions far beyond the limiting Weissenberg numbers in the standard
scheme. However it turns out that, although in large parts of the flow the solution
converges, we have not been able to obtain convergence in localized regions of the
flow. Possible reasons for this are: artefacts of the model (Oldroyd-B) or unresolved
small scales (Giesekus). However, more work is necessary, including the use of more
refined meshes and/or higher-order schemes, before any conclusion can be made on
the local convergence problems.
Key words: finite element method, matrix logarithm, flow around a cylinder,
viscoelastic flow simulation
Preprint submitted to J. Non-Newtonian Fluid Mech. 7 September 2004
1 Introduction
The high Weissenberg number problem (hwnp) has been the maj or stumbling
block in computational rheology for the last three decades. The term “hwnp”
coins the empirical observation that all numerical methods break down when
the Weissenberg number exceeds a critical range. The precise critical value at
which computations break down varies with the problem, the method and the
mesh used. It is now widely accepted that the hwnp arises when computations
are not resolved well enough to capture the sharp stress gradient that form
within the flow, yet, the precise mechanism that leads to the breakdown has
remained somewhat of a mystery (see [2] and the references therein).
The mechanism responsible for the hwnp has been recently explained in [3].
In essence, it is a numerical instability caused by the failure to balance the
exponential growth of the stress (due to deformation) with convection. Such
a failure is common to all methods that approximate the stress by polynomial
base functions (the choice of polynomial base function is explicit in finite
elements and implicit in finite differences). The proposed remedy was a change
a variables into new variables that scale logarithmically with the stress tensor.
Specifically, the constitutive relations were reformulated as equations for the
matrix logarithm of the conformation tensor, exploiting the fact that the latter
is symmetric positive definite. Numerical experiments for a two-dimensional
lid-driven cavity using a second-order finite difference scheme indicate that
schemes based on the new logarithmic formulation are immune to the high
Weissenberg breakdown. The “old” hwnp is rather replaced with a “new”
problem: accurate computations at high Weissenberg require high resolution.
To some extent, the situation may now be compared to classical CFD, whe re
stable calculations can be performed at arbitrarily large Reynolds numbers,
but accuracy is lost when resolution is insufficient.
In this paper we implement the new logarithmic formulation with a finite
element method (fem), and test it for a classical benchmark problem—flow
past a cylinder. More specifically, we have in mind the following goals:
(1) Develop a fem scheme immune to the high Weissenberg instability.
(2) Get more insight into the hwnp by analyzing the numerical breakdown
within a fem point of view.
(3) Present benchmark results that confirm the validity of the new method
at low and moderate Weissenberg numbers, where comparison to existing
data is possible.
∗
Corresponding author. Tel: +31-40-247-5081; Fax: +31-40-244-7355.
Email addresses: m.a.hulsen@tue.nl (Martien A. Hulsen),
raananf@cs.huji.ac.il (Raanan Fattal), raz@math.huji.ac.il (Raz
Kupferman).
2
(4) Establish new benchmark results at higher Weissenberg numbers.
(5) Investigate the limitations of the method, and in particular, understand
how accuracy is lost at high Weissenberg numbers. Such an understanding
is important for the future design of higher-order methods.
The paper is structured as follows: In Section 2 we describe the governing
equations and introduce the notations used henceforth. In Section 3 we de-
scribe the standard fem discretization based on DEVSS and discontinuous
Galerkin. In Section 4 we analyze the stability of the standard discretization,
and in particular obtain a simple stability criterion whose violation causes the
high Weissenberg breakdown. In Section 5 we derive the constitutive equation
for the matrix logarithm of the conformation tensor; the derivation is based
on a slightly different approach than in [1]. In Section 6 we present numerical
results for the flow past a cylinder for an Oldroyd-B and a Giesekus model. In
particular we will verify the criterion for breakdown in the standard discretiza-
tion. We also show the much improved stability of the matrix log formulation.
A discussion follows in Section 7.
2 Governing equations
We will consider the flow of viscoelastic fluids at creeping flow conditions
(inertia can be neglected) in a spatial region Ω, having a boundary denoted
by Γ. For this, we need the following set of equations: the momentum balance,
the mass balance for incompressible flows and a constitutive equation decribing
the stress.
The momentum balance and mass balance are given by
∇p − ∇ · (2η
s
D) − ∇ · τ = 0, (1)
∇ · u = 0, (2)
where p is the pressure, u is the velocity vector, the rate-of-deformation tensor
is given by D =
1
2
(L + L
T
), with the velocity gradient tensor L = (∇u)
T
and τ is the extra-stress (or polymer stress). The coefficient η
s
is the solvent
viscosity. The polymer stress τ is given by
τ =
η
p
λ
(c − I), (3)
where c is the conformation tensor, η
p
is the zero-shear-rate viscosity of the
polymer part of the stress and λ is the relaxation time. We will use two different
models for the evolution of the conformation tensor: the Oldroyd-B and the
Giesekus model, which can be written as follows
˙
c = L · c + c · L
T
+ f(c), (4)
3
with
f(c) =
−
1
λ
c − I
Oldroyd-B,
−
1
λ
c − I + α(c − I)
2
Giesekus.
(5)
The Oldroyd-B model is identical to a Giesekus model with α = 0. We will
solve Eq. (4) in a Eule rian frame and theref ore write
˙
c =
∂c
∂t
+ u · ∇c. (6)
Boundary conditions will be dis cussed for the flow around a cylinder problem
that will be analysed in Sec. 6.
3 Numerical discretization
For the spatial discretisation of the system of equations we will use the finite
element method (fem). We will basically use the same implementation as
described in [4] and give a brief summary here.
In order to obtain a proper mixed velocity-stress formulation we use the Dis-
crete Elastic-Viscous Split Stress (DEVSS) formulation of Gu´enette & Fortin
[5] for the discretisation of the linear momentum balance and the continu-
ity equation. For this we introduce an extra variable e = 2η
p
D. Note, that
Gu´enette & Fortin [5] introduce D as an extra variable, however using e in-
stead leads to a symmetric system matrix. We rewrite the momentum balance
Eq. (1) and the continuity equation Eq. (2) as follows
∇p − ∇ · (2η
s
D(u) + τ ) − ∇ · (2η
p
D(u) − e) = 0, (7)
−∇ · u = 0, (8)
−D(u) +
1
2η
p
e = 0. (9)
For the weak f ormulation of Eqs. (7)–(9) we introduce separate functional
spaces for u, p and e, which we denote by U, P and E, respectively. The weak
formulation can be found by multiplying with testfunctions and integration
by parts: find (u, p, e) ∈ U ×P ×E such that for all (v, q, g) ∈ U ×P ×E we
have
−(∇ · v, p) + (∇v, 2ηD(u) − e + τ ) = (v, σ)
Γ
, (10)
−(q, ∇ · u) = 0, (11)
−(g, ∇u) +
1
2η
p
(g, e) = 0, (12)
4
where (·, ·) and (·, ·)
Γ
are proper L
2
inner products on the domain Ω and on
the boundary Γ, respectively. The viscosity η is the zero-shear-rate viscosity
η
s
+ η
p
and σ is the traction vector on the boundary. The system for (u, p, e)
is symmetrical.
The discrete form of the equations is obtained by requiring that the weak form
is valid on approximating subspaces U
h
× P
h
× E
h
which consist of piecewise
polynomial spaces. The discrete solutions and the discrete testfunctions are
denoted with subindex
h
: (u
h
, p
h
, e
h
) and (v
h
, q
h
, g
h
). We will us e quadrilat-
eral elements with continuous biquadratic polynomials (Q
2
) for the velocity
space U
h
, discontinuous linear polynomials (P
1
) for the pressure space P
h
and
continuous bilinear polynomials (Q
1
) for the viscous polymer stress space E
h
.
For the discretisation of the constitutive equation we use the discontinuous
Galerkin method (DG) [6]. This leads to the following weak formulation of
the constitutive equation Eq. (4): find c ∈ T on all elements e
i
such that for
all w ∈ T we have
w,
∂c
∂t
+ u · ∇c − L · c − c · L
T
− f(c)
e
i
+
Z
γ
in
i
(n · u)w : (c
+
− c) dγ = 0, (13)
where (·, ·)
e
i
denotes an L
2
inner product on element e
i
only, γ
in
i
is the part of
the element boundary where u ·n < 0 (inflow boundary), n is the outwardly
directed unit normal on γ
in
i
. Furthermore, c
+
is the value of c in the upstream
neighbour element or the imposed value at the inflow boundary part of Γ.
The functional space for c is denoted by T . We will use discontinuous bilinear
polynomials (Q
1
) for the space T
h
. The combination of the DEVSS formulation
with DG has been introduced by B aaijens et al. [7] and it has been proved by
Fortin et al. [8] that it leads to a proper mixed velocity-stress scheme.
For the time discretisation of the constitutive equation Eq. (13) we use an
explicit Euler scheme, where the time derivative is discretised by ∂c/∂t ≈
(c
n+1
− c
n
))/∆t and all other terms are evaluated at t
n
. For time-accurate
solutions we can use a higher-order scheme, such as Adams-Bashforth, but for
obtaining steady state solutions this is not necessary. Note, that the equations
can be solved at element level. Next we substitute τ
n+1
= η
p
/λ(c
n+1
−I) into
Eq. (10). The system matrix for solving (u
n+1
, p
n+1
, e
n+1
) is symmetrical and
LU decomposition is performed at the first time step. Since this matrix is con-
stant in time, solutions at later time steps can be found by back substitution
only. This results in a significant reduction of the CPU time.
5
