arXiv:1508.01041v2 [cs.CE] 20 Apr 2016

Implementation of the Log-Conformation Formulation for
Two-Dimensional Viscoelastic Flow
Kristian Ejlebjerg Jensen
kristianejlebjerg+cpc@gmail.com
Department of Micro- and Nanotechnology
Technical University of Denmark

Fridolin Okkels
Department of Micro- and Nanotechnology
Technical University of Denmark

Peter Szabo
ps@kt.dtu.dk
Department of Chemical and Biochemical Engineering
Technical University of Denmark
April 21, 2016

Abstract
We have implemented the log-conformation method for two-dimensional viscoelastic flow in COMSOL, a
commercial high-level finite element package. The code is verified for an Oldroyd-B fluid flowing past a
confined cylinder. We are also able to describe the well-known bistability of the viscoelastic flow in a crossslot geometry for a FENE-CR fluid, and we describe the changes required for performing simulations with
the Phan-Thien-Tanner (PTT), Giesekus and FENE-P models. Finally, we calculate the flow of a FENECR fluid in a geometry with three in- and outlets. The implementation is included in the supplementary
material, and we hope that it can inspire new as well as experienced researchers in the field of differential
constitutive equations for viscoelastic flow.

1

Introduction

The flow of viscoelastic fluids in complex geometries has many practical applications due to the fact that
any dissolved elastic component, be it polymers or biological molecules, will give rise to viscoelastic effects.
Such effects require the use of models that take the fluid’s memory of past deformations into account.
Differential constitutive equations is able to do just that, and they provide a good quantitative agreement
with experiments [1]. Complex geometries call for the use of numerical methods, and therefore significant
effort has been devoted towards various discretization techniques as well as model reformulations. The
appearance of singularities, that caused codes to break down in smooth geometries at moderate elasticity,
led to the definition of a ”High Weissenberg Number Problem” (HWNP). This was effectively solved, when
loss of positive definiteness for the conformation tensor was recognized as the cause of the singularity and
the log-conformation method was introduced as a remedy [2].

1

Commercial numerical tools provide a simple workflow due to the integration of geometry description,
unstructured mesh generation, discretization, solvers, and post-processing, which makes the treatment of
complex geometries effortless compared to research grade code, but to our knowledge there does not exist
any commercial tool that implements the log-conformation method. This is perhaps due to the fact that
the method involves the calculation of eigenvectors and eigenvalues for the conformation tensor, which
complicates the formulation of the constitutive equation. In two (and three) dimensions explicit expressions
exists for these quantities, and it is thus possible to formulate the governing equations solely in terms of
functions recognized by commercial simulation packages. We have included a COMSOL implementation of
the Oldroyd-B model with log-conformation reformulation in the supplementary material, and we hope this
can help other researchers working with differential constitutive equations.
The article is structured in four parts: First we discuss implementation details of the log-conformation
method for the most simple constitutive equations. We then calculate the flow of an Oldroyd-B fluid past a
confined cylinder using COMSOL for comparison with other works [2]. As a third point we demonstrate the
ability to predict bistable behavior of a FENE-CR fluid in a cross-slot geometry [3]. Finally we consider a
geometry with three in- and outlets, which has not previously been described in the context of viscoelastic
flow, at least to our knowledge.

2

Log-Conformation Implementation

In the so-called dumbbell models, the elastic part of a viscoelastic fluid is approximated by two springconnected point masses – an elastic dumbbell. The orientation and elongation of this dumbbell is described
using the end-to-end vector a with the equilibrium length aeq , see figure 1. The conformation tensor, A, is

Figure 1: An elastic dumbbell is illustrated with its two spring-connected point masses, and the end-to-end
vector describing orientation and extension.
the statistical average, h...i of the dyadic product between the end-to-end vector and itself normalized with
the squared equilibrium length,
A =

ha ⊗ ai
.
a2eq

(1)

The equation system for the Oldroyd-B model is written below. We neglect the effect of temperature
and compressibility, but keep inertia for the sake of generality.


τ
s
z
}|
{


Dv


ρ
= ∇ · −pI + ηs ∇ v + (∇ v)T +τ e 
(2)
Dt
∇·v
τe
0

=

0

ηp
=
A−I
λ
1
= − (A − I)
λ

DA 
+
− A · ∇ v + (∇ v)T · A
Dt

(3)
(4)

(5)

where D is the material derivative, ρ is the density, v is the velocity vector, t is time, p is the pressure ηs is
the solvent viscosity, ηp is the elastic viscosity, λ is the relaxation time, τ s is the solvent stress tensor and
2

τ e is the elastic stress tensor. The Navier-Stokes equation (2) guarantees that Newton’s 2nd law is obeyed,
while the continuity equation (3) insures mass conservation. The expression for the elastic stress tensor
in equation (4) is best understood in light of the conformation tensor definition (1) with ηp /λ acting as a
spring constant. The evolution of the conformation tensor is described in equation (5) with relaxation on
the first line balanced by the upper convected time derivative, which consists of convection in terms of the
total derivative together with rotation and extension of the conformation tensor as imposed by the velocity
gradient in square brackets.
When the velocity vector is approximated by C 0 continuous polynomials, the velocity gradient, ∇ v
becomes discontinous. This can cause numerical difficulties [4], and therefore aC 0 continuous approximation

G is usually constructed for use in equation (5). Furthermore the zero, ηp ∇ v + (∇ v)T − G − GT

is

added on the right-hand side of equation (2) to preserve the elliptic nature of the equation for small solvent
viscosities in what is called the discrete elastic viscous stress splitting (DEVSS). Similarly, we use streamline
upwind diffusion (SUPG) to introduce an elliptic component in the constitutive equation (5), see section 2.2.
The conformation tensor is symmetric, and the point of the log-conformation method is to exploit this
property to rewrite equation (5) to the form
Ds
− Π(s, G) = 0,
Dt
where s is the matrix logarithm of the conformation tensor. They are related by the following expression
= R · e s · RT


v1nx −v1ny
=
v1ny v1nx
 λ

e 1
0
=
0 eλ2

A
R
es
A11

2
2
= v1nx
eλ1 + v1ny
eλ2

A12

= v1nx v1ny eλ1 − eλ2

A22

2
2
= v1nx
eλ2 + v1ny
eλ1



where λ1 and λ2 are the eigenvalues of s, while v1nx and v1ny are the normalized eigenvector components
belonging to λ1 . The orthonormal rotation matrix, R holds the components of the eigenvectors, but we
exploit that the second eigenvector is just the transpose of the first. Denoting the components of s by s11 ,
s12 and s22 , the eigenvalues become1


q
1
λ1 = 2 s22 + s11 − (s22 − s11 )2 + 4s212
and


q
1
λ2 = 2 s22 + s11 + (s22 − s11 )2 + 4s212 ,
and the components of the first eigenvector are


λ1 − s22
s12

v1

=



v1x

.q

v1y

.q


.

The eigenvector is normalized
v1n = 

2 + v2
v1x
1y

2 + v2
v1x
1y


.

The special case s12 = 0 is handled by using R = I, when |s12 | ≤  with  = 10−12 .
1 There

also exist explicit expressions for the three-dimensional case.

3

In the following we will express the reaction term Π as R · Ω · RT as described in [5]. The continuous
velocity gradient is rewritten in the principal frame
G̃ =

RT · G · R

G̃11

=

2
2
v1nx
G11 + v1nx v1yx (G12 + G21 ) + v1ny
G22

G̃12

=

2
2
v1nx
G21 + v1nx v1yx (G22 − G11 ) − v1ny
G12

G̃21

=

2
2
v1nx
G12 + v1nx v1yx (G22 − G11 ) − v1ny
G21

G̃22

=

2
2
v1nx
G22 − v1nx v1yx (G12 − G21 ) + v1ny
G11 .

Note that this matrix is not symmetric, and that we have chosen G12 to refer to the derivative of the velocity
vector x-component with respect to the y-coordinate. The diagonal components Ω11 and Ω22 are
Ω11

=

Ω22

=

1 − e−λ1
and
λ
1 − e−λ2
.
2G̃22 −
λ
2G̃11 −

These expressions are specific to the Oldroyd-B model, but it is straight forward to write the expressions for
the Giesekus, Phan-Thien-Tanner, FENE-P or FENE-CR models instead [5].
Giesekus

:

where
PTT(exp.)

:

where

(6)

1 + αGie. (ex − 1) x
(e − 1)
λ
1 ,λ2 )
2G̃11 − kPTT (λ
λ
kPTT (λ1 ,λ2 ) ,
2G̃22 −
λ

Ω11 =
Ω22 =

kPTT (λ1 , λ2 ) = ePTT (λ1 +λ2 −3)
:

where
FENE − CR

2G̃11 − fGie. (λ1 )
2G̃22 − fGie. (λ2 )

fGie. (x) =

where
FENE − P

Ω11 =
Ω22 =

Ω11 = 2G̃11 −
Ω22 = 2G̃22 −
k(λ1 , λ2 ) =

:

(7)

)eλ1 −1

k(λ1 ,λ2
λ
k(λ1 ,λ2 )eλ2 −1
λ

,

1
1 − (λ1 + λ2 )/max2
λ1

2G̃11 − k(λ1 , λ2 ) e λ−1
,
λ2
2G̃22 − k(λ1 , λ2 ) e λ−1
1
k(λ1 , λ2 ) =
1 − (λ1 + λ2 )/a2max

Ω11 =
Ω22 =

In the case of the FENE-P and FENE-CR models, the amax is a maximum extension, i.e. Tr(A) < a2max .
The off-diagonal component, Ω12 , equals


( 
λ1 −λ2
λ1
λ2
e
G̃
+
e
G̃
, |λ1 − λ2 | > 
12
21
λ
λ
e 1 −e 2
Ω12 =
G̃12 + G̃21
, |λ1 − λ2 | ≤ 
where we have reused the numerical parameter for the eigenvectors to handle the special case of identical
eigenvalues. We are now ready to compute the reaction term

Π12

2
2
= v1nx
Ω11 − 2v1nx v1ny Ω12 + v1ny
Ω22

2
2
= v1nx − v1ny Ω12 + v1nx v1ny (Ω11 − Ω22 )

Π22

2
2
= v1nx
Ω22 + 2v1nx v1ny Ω12 + v1ny
Ω11

Π11

To ease modification of the implementation, we have opted to construct functions for each of the above
expressions instead of writing out Π in terms of expressions involving s and G only.
4

With respect to visualization it is worth noting that s and A are positive definite tensors, and one
can thus think of them as ”fields of ellipses”. In other words, the conformation tensor has an extension in
every point given by its eigenvectors and eigenvalues. The trace of the conformation tensor is the sum of
the extensions in the principal directions given by the eigenvectors, while the determinant is the product
of these extensions. The trace and determinant equal the sum and product of the eigenvalues, respectively,
and they are invariant with respect to rotation of the coordinate system, which make them interesting in
the context of post-processing/visualization2 . The last invariant is the angle(s) related to the eigenvectors,
but the eigenvector associated with the largest eigenvalue is normally oriented with the streamlines, so this
angle provides little insight. Furthermore the trace and determinant often look very similar, at least on a
logarithmic scale, and therefore we just plot the trace.

2.1

Non-Dimensional Equations

We define the dimensionless spatial coordinates, velocity, time, pressure and stress (x̃, ṽ, t̃, p̃ and τ̃ e ),
x = Lchar x̃,

v

=

p

=

and τ e

=

Lchar
t̃,
vchar
vchar (ηs + ηp )
p̃
Lchar
vchar
(ηs + ηp )
τ̃ ,
Lchar e
vchar ṽ,

t=

which give rise to the following dimensionless constants
Re

=

β

=

ρvchar Lchar
,
η
ηs
,
ηp + ηs

We = λ

vchar
Lchar

and

The Weissenberg number, We, describes the magnitude of elastic to viscous effects. When either the Weissenberg number goes to zero or the solvent to total viscosity ratio, β goes to unity, the model approaches that
of a Newtonian fluid. Contrarily, the value β = 0 corresponds to the upper convected Maxwell model, which
is characterized by particularly strong viscoelastic effects. Note that we prefer to define a characteristic
pressure in terms of a characteristic velocity rather than the other way around due the fact that we intend
to impose fixed flow rates rather than fixed driving pressures (as in [6]).
We have implemented the above approach to the log-conformation formulation for two dimensions in
COMSOL Multiphysics version 4.3a, and the online version of the manuscript includes a MATLAB script
for computation of the flow of an Oldroyd-B fluid past a confined cylinder using the COMSOL LiveLink
interface. We use default settings for everything, except we force update of the Jacobian every time a time
step is saved.

2.2

The Oldroyd-B model in Weak Form

Converting the strong form of the Oldroyd-B model (2-5) to weak form is trivial for the equations that
relate to the elastic stress (4) and continuity equation (3) as one just multiplies with test functions for the
elastic stress and pressure, respectively. The evolution equation for the conformation tensor (5) requires
SUPG stabilization, which means that the equation is multiplied3 with Atest + (2hmesh /vchar )v · ∇Atest
rather than Atest , where hmesh is a characteristic size for the mesh. In other words a small amount of
stream-wise diffusion is added to allow computation, but since the amount of diffusion is scaled with the
mesh, the implementation still converges towards the solution without any artificial diffusion. Finally there
is the momentum part of the Navier-Stokes equation (2), which is multiplied with a test function for the
P
Q
that det(A) = i eλi = e iλi = eTr(s) , where λi are the eigenvalues of s.
that the test functions for elastic stress and conformation tensor are multiplied element wise (Hadamard product),
denoted with ◦.
2 Note

3 Note

5

velocity and integrated by parts:
Z 
0

=
Ω




Dṽ
˜
· ṽtest dΩ
∇ · τ̃ s + τ̃ e − Re
Dt̃

Z

0

=
+
+

˜ ṽtest dΩ
−
τ̃ s ◦ ∇
Ω

Z 
Dṽ
−1
˜
−
Re
−
Da
α̃ṽ
· ṽtest dΩ
∇τ̃
e
Dt̃
Ω
Z
τ̃ s · n̂ · ṽtest dl

(8)

(9)

∂Ω

The point is that the solvent stress is integrated by parts (Divergence theorem), but the elastic stress is not.
This approach has been followed by other researchers so we have adopted it as well. The advantage is that
the specification of τ̃ s on open boundaries is slightly more simple than τ̃ s + τ̃ e , but in fact we have tested
both without finding any difference. At outlets we thus usually impose normal velocity, zero pressure and
τ̃ s = −Ip̃.
This can also be used for pressure driven inlets, but in that case one should supplement the constraint on
the velocity with one on the conformation tensor. We impose the no-slip(ṽ = 0) at walls, which causes
COMSOL to automatically enforce ṽtest = 0, such that the boundary integral in equation (9) vanishes.

3

Confined Cylinder

The confined cylinder geometry is shown in figure 2, and it is a popular benchmark geometry due to the
absence of geometric singularities. It can be used to verify an implementation by comparing the drag over
the cylinder with a reference. We exploit the symmetry of the problem by modeling only half of the domain.

Figure 2: The confined cylinder geometry with a (dashed) symmetry axis. The fluid comes in from the left,
flows past the obstacle and exits through the boundary on the right (t̂ is the unit tangent vector).
The cylinder radius R, and average inlet velocity are used as characteristic length and velocity, respectively. The blocking ratio is set at 2 (H = 2R), and the outlet length, L1 , is fixed at 10R. In agreement
with references [2] the solvent to total viscosity ratio is taken as 0.59, and inertia is neglected (Re = 0). The
dimensionless drag is computed as
Z
Drag = −2
τ̃ · n̂ · x̂T dl, with
cylinder
h
i
˜ + (∇ṽ)
˜ T + τ̃
τ̃ = −p̃I + ηs ∇ṽ
e
where n̂ is the outward pointing normal vector.
We impose the no-slip boundary condition ṽ = 0 at the walls, and fully developed boundary conditions

6

at the inlet
ṽinlet
A11,inlet
A12,inlet
A22,inlet

3
2


1 − (ỹ/2)2 x̂

2
∂ (ṽ · x̂)
9
= 1+2
We = 1 + 8 ỹ 2 We2
∂ ỹ
∂ (ṽ · x̂)
3
= We
= − 4 ỹWe and
∂ ỹ
= 1
=

(10)

Note that the intermediate expressions can be used in the case of a pressure driven setup. At the outlet we
impose zero pressure together with zero normal viscous stress.
poutlet

=

τ̃ s,outlet · n̂

=

0

and

−Ip · n̂ = 0

(11)

At the symmetry line we impose zero normal velocity.
ṽsym · n̂ = 0

3.1

(12)

Solution Details

In order to find a steady solution one can go through the following steps.
#1 Solve for the case of a Newtonian fluid.
#2 Use #1 as initial condition for a transient simulation running for 10 relaxation times.
#3 Use #2 as initial guess for a non-linear solver to find a steady solution.
We however choose to vary the Weissenberg number slowly using a regularized step function, st, such that
a quasi steady state is achieved.
We
st(t̃)

Westart + (Weend − Westart ) st(t̃), where

,
t̃ < t̃start
 0
3
=
0.5 + 1.5t̄ − 2 (t̄ ) , t̃start ≤ t̃ < t̃end

1
, t̃end ≤ t̃
=

t̄ =

t̃ − (t̃start + t̃end )/2
t̃end − t̃start

(13)

This way we compute a quasi-steady solution to a range of Weissenberg numbers using a single transient
simulation. COMSOL uses a fully implicit transient scheme with adaptive stepsize, so slow variation of the
Weissenberg number does not increase the total computation time [7].

3.2

Results

Since the time over which the Weissenberg number is varied is a numerical parameter, we choose to scale it
with the mesh size, so convergence towards the steady case is achieved, i.e.
Tstep = t̃end − t̃start = 8000(0.07/hmesh )
Figure 3 shows the drag as a function of the Weissenberg number for three different discretizations indicating
convergence up to We = 0.6. We observe good agreement with the reference drag although the simplicity
of the chosen mesh gives rise to a high computational cost. It is well known, that convergence at We = 0.7
is difficult, even with anisotropic mesh adaption [8]. In fact experimental [9] as well as theoretical evidence
[10] exist to suggest, that a transition to a three dimensional flow pattern occurs, in other words the two
7

Figure 3: The drag on a confined cylinder is plotted as a function of the Weissenberg number for three
discretizations together with the data of another study [2].

Figure 4: The trace of the conformation tensor is plotted on a logarithmic scale for the confined cylinder in
the case of We = 0.51 and We = 0.7 (183k DOF). The acceleration of the flow in the wake of the cylinder
generates a strand of highly elongated dumbbells.
dimensional symmetric flow becomes unstable. The three-dimensional solution has the flow alternating
between going above and below the cylinder, and one can think of this solution as having a weaker extensional
character in the wake of the cylinder.
When a viscoelastic fluid enters a region, where the extension rate times the relaxation time reaches
a critical value, the dumbbell extension starts to grow exponentially in time. This kind of non-linearity
explains the significant difference between the wake of extension at We = 0.51 and We = 0.7 illustrated in
figure 4.

4

Bistable Cross-Slot



The FENE-CR model differs from the Oldroyd-B model in the sense, that λ is replaced by λ 1 − Trace(A)/a2max .
In other words the relaxation time decreases as the dumbbell extension increases. This puts an upper bound
a2max , on the trace of the conformation tensor corresponding to a maximum extension. The modification
prevents an unbounded extension in free stagnation points such as the center of the cross-slot geometry
shown in figure 5. This geometry is known to exhibit bistability experimentally [11] as well as numerically
[12].
In this case we keep the average inlet velocity as characteristic velocity and the channel width H as

8

Figure 5: The cross-slot geometry is shown with the fluid coming in from the left and right, diverging at
the center, before exiting through to the upper and lower outlets leaving a free stagnation point (red) in the
center of the geometry.
characteristic length scale. We impose the no-slip boundary condition at walls, but we have been unable
to initialize the transient simulation with the correct expressions for the conformation tensor at the inlets.
Therefore we have used s = I at the inlets rather than values based on the expressions listed in Appendix
I. Furthermore we had to change the damping strategy of the nonlinear method from the default (constant)
to automatic.
At the outlet we impose the same boundary conditions (11) as for the cylinder. We compute the dissipation as an integral over the entire computation domain,
Z
h
i
˜ + (∇ṽ)
˜ T dΩ.
φ̃ =
τ̃ : ∇ṽ
(14)
Ω

where : is the Frobenius product. We choose Re = 0, a2max = 100 and β = 0.2 and the same transient
solution procedure as described in section 3.1. We define the center of the domain,
Ωc

=

x̃ ∈ |x̃|∞ < 0.5,

as we find that the vorticity integrated over this is a good asymmetry parameter. We thus plot the square
of this versus the Weissenberg number in figure 6(a). The system becomes bistable at We = 0.54 , which
gives rise to a kink in the dissipation as shown in figure 6(b). The solution before the kink [figure 7(a)
and (d)] has vertical, horizontal as well as 180-degree rotational symmetry, while the solutions after the
kink [7(b), (c), (e) and (f)] only have 180-degree rotational symmetry, and we will refer to these solutions
with reduced symmetry as ”asymmetric solutions”. The phenomenon can be understood by considering the
velocity magnitude and conformation tensor trace as plotted in figures 7(a-c) and (d-e) respectively. The
dumbbells are extended along the vertical axis, which causes a damping that delays the merging of the
flows, and this gives rise to extra shear and thus also dissipation. When the solution becomes asymmetric,
the damping is moved to the side of the channel, which lowers the dissipation. One can also think of the
asymmetry as a flow pattern with less extension in the stagnation point, similar to the three-dimensional
flow for the confined cylinder.
4 Based on personal correspondence we have concluded that the results of [3] are based on the FENE-MCR model, so no
reference exists for the FENE-CR model.

9

Figure 6: The square of the vorticity integrated over the center of the domain Ωc is plotted as function
of the Weissenberg number in (a). In (b) it is the dissipation that is plotted, which shows that the point
of bistability occurs (We = 0.5) at the point of maximum dissipation and thus also maximum hydraulic
resistance.

5

Three Inlets

In this section we introduce the geometry shown in figure 8. It consists of a hexagon with in- and outlets
attached in an alternating fashion. The geometry always has a vertical symmetry axis, but depending on
whether a certain angle is equal to π/3 or not, it can also have 120-degree rotational symmetry. We define
the distance to the center as in- and outlet lengths, and set equal inlet flow rates as well as outlet pressures.
We set Lin = 6H and Lout = 8L, β = 0.1 and s = I at the inlets, but keep all other parameters identical
to those used for the cross-slot. That is, the channel width and average inlet velocity are still used as
characteristic length and velocity, respectively.
If the angle is equal to π/3 radians, a single stagnation point exists in the center of the geometry regardless
of whether we consider Stokes flow or just the regime of low elasticity. Two stagnation points however appear,
if the angle is different from π/3 radians. Figure 9 shows the solution for angles of π/3, π/3.5 and π/4 radians
at We = 0.37 (low) and We = 0.66 (moderate). For an angle of π/3 radians, the regime of low elasticity (a)
gives a solution with 120-degree rotational- and vertical symmetry, while for moderate elasticity a solution
with only 120-degree rotation symmetry appears (b). This latter solution involves the clockwise flow, but by
symmetry we can argue that there must be an identical solution with counter-clockwise flow. The symmetry
reduction also occurs for an angle of π/3.5 (c-d), in the sense that we go from only vertical symmetry to
no symmetry. Regardless of the Weissenberg number, there are two stagnations points even though the
flow rate between them is small at moderate elasticity (d). In other words turning up the elasticity causes
the flow rate between the stagnation points to decrease significantly. As shown in (e-f) this is different for
an angle of π/4 radians, as the elasticity tends to move the stagnation points apart and increase the flow
10

Figure 7: The magnitude of the velocity vector is plotted for three Weissenberg numbers in the cross geometry
in (a-c), while red stream lines are plotted on top of the trace of the conformation tensor on a logarithmic
scale in (d-f). Note that we plot streamlines as contours of the stream function instead of integrating the
velocity field.

Figure 8: The geometry with three in- and outlets is illustrated with the different length scales and angles
involved. For an angle of π/3 the geometry has a symmetry axis through every channel as well as 120-degree
rotational symmetry, while an angle different from π/3 breaks all symmetries except the one around the
vertical axis.
rate between them in what we call the ”straight flow solution”, that is the solution has identical symmetry
properties at low and moderate elasticity.
Using artificial initial conditions as listed in Appendix II, one can show that the clockwise, counterclockwise and the straight flow all exist as stable solutions at moderate elasticity for angles between π/4 and
π/3.5. This means, that it is the connectivity with the solution at low elasticity that changes. Furthermore
there is an unstable solution with a vertical symmetry axis for the angles π/3 and π/3.5. We calculate the
integral of the vorticity around the center
Ωo

= x ∈ ||x||2 < 0.5,
11

Figure 9: The trace of the conformation tensor is plotted on a logarithmic scale together with streamlines
(red) for two Weissenberg numbers in a geometry with three inlets. The case of an angle equal to π/3 (a-b),
π/3.5 (c-d) and π/4 (d-e) radians is illustrated (see figure 8). Although not apparent from the streamlines,
the upper outlet receives fluid from all three inlets in (d), and the flow between the stagnation points is
reduced with a factor of 4.5 in (d) compared to (c). Any asymmetry around the vertical axis in (a), (c) and
(e) is due to the mesh.
as a measure of asymmetry and plot this as a function of the Weissenberg number for five angles in figure
10 (a), while the dissipation is plotted in figure 10 (b). It is clear that the point of maximum dissipation
coincides with the point of bistability for the three largest angles (π/3, π/3.5, π/3.75), as it was also the
case for the cross-slot in section 4. The dissipation for the two smaller angles also have maxima, but we do
not believe they play the same role as the maxima for the smaller angles. This is due to the fact that the
Weissenberg number at which the maximum occurs seems to decreases for both the rotation and straight
flow solutions. This means that there does not exists an angle for which the two kind of maxima in figure
10 (b) meet, although there must exist an angle where all three solutions are connected to the solution at
low elasticity.

6

Conclusion

We have described the viscoelastic Oldroyd-B model with log-conformation reformulation in detail, and
included a COMSOL implementation in the supplementary material. We have showed, that the code agrees
with a reference in a benchmark geometry. Furthermore we have demonstrated that the code can be modified
to describe the flow of a FENE-CR fluid, and that it can reproduce bistability in the cross-slot geometry.
Finally we computed the flow in a geometry with three in- and outlets, which appears to feature three stable
solutions. The connectivity of the solutions with the solution at low elasticity seems to depend on an angle
in the geometry.

12

Figure 10: The square of the vorticity integrated over the center of the domain Ωo is plotted as functions
of the Weissenberg number for five different angles in (a), while the dissipation is shown (b). The point of
bistability coincides with the point of maximum dissipation for the three larger angles, as it was also the case
for the cross-slot in figure 6. The maximum dissipation of the two smaller angles is indicated with crosses.

References
[1] F. Baaijens, S.H.A. Selen, H.P.W. Baaijens, G.W.M. Peters, and H.E.H. Meijer. Viscoelastic flow past
a confined cylinder of a low density polyethylene melt. Journal of non-newtonian fluid mechanics, 68
(2-3):173–203, 1997.
[2] M.A. Hulsen, R. Fattal, and R. Kupferman. Flow of viscoelastic fluids past a cylinder at high weissenberg
number: stabilized simulations using matrix logarithms. Journal of Non-Newtonian Fluid Mechanics,
127(1):27–39, 2005.
[3] G.N. Rocha, R.J. Poole, M.A. Alves, and P.J. Oliveira. On extensibility effects in the cross-slot flow
bifurcation. Journal of Non-Newtonian Fluid Mechanics, 156(1):58–69, 2009.
[4] A. Fortin, R. Guenétte, and R. Pierre. On the discrete evss method. Computer Methods in Applied
Mechanics and Engineering, 189(1):121–139, 2000.
[5] A. Kane, R. Guénette, and A. Fortin. A comparison of four implementations of the log-conformation
formulation for viscoelastic fluid flows. Journal of Non-Newtonian Fluid Mechanics, 164(1-3):45–50,
2009. ISSN 0377-0257.
[6] K. Ejlebjerg Jensen, P. Szabo, and F. Okkels. Topology optimization of viscoelastic rectifiers. Applied
Physics Letters, 100(23):234102–234102, 2012.

13

[7] COMSOL Multiphysics. 4.3a user’s guide, 2013.
[8] R. Guénette, A. Fortin, A. Kane, and J.F. Hétu. An adaptive remeshing strategy for viscoelastic fluid
flow simulations. Journal of Non-Newtonian Fluid Mechanics, 153(1):34–45, 2008.
[9] G.H. McKinley, R.C. Armstrong, and R.A. Brown. The wake instability in viscoelastic flow past confined
circular cylinders. Philosophical Transactions: Physical Sciences and Engineering, 344(1671):265–304,
1993.
[10] M. Sahin. Parallel large-scale simulation of viscoelastic fluid flow instabilities. In 17th International
Workshop on Numerical Methods for Non-Newtonian Flows, 2012.
[11] P.E. Arratia, CC Thomas, J. Diorio, and JP Gollub. Elastic instabilities of polymer solutions in crosschannel flow. Physical review letters, 96(14):144502, 2006.
[12] RJ Poole, MA Alves, and PJ Oliveira. Purely elastic flow asymmetries. Physical review letters, 99(16):
164503, 2007.

Appendix I: FENE-CR Inlet Conditions
The developed flow conditions for a FENE-CR fluid
3

5

is listed in the following equations

1 − (ỹ/2)2 n̂

ṽinlet

=

−2

A11,inlet

=

a2max − 1


p
a2max a2max − a4max + 8χ2 (a2max − 2)

+

4χ2


A11,inlet + 1
A12,inlet = χ 1 −
a2max
A22,inlet = 1, where
∂ (ṽ · x̂)
= 3We(ỹ/4)(n̂ · x̂)
χ = We
∂ ỹ

Appendix II: Artificial Initial Conditions for the Geometry with
Three In- and Outlets
One can find stable solutions not connected to the solution at low elasticity using artificial initial conditions.
A rotating solution for an angle of π/4 (see figure 8) can in example be found by imposing a rotating force
in the center of the geometry,



F =



(x̂ỹ − ŷx̃)/2
0

, ||r̃||
, ||r̃||

≤ 0.5/ sin(π/6)
,
> 0.5/ sin(π/6)

where r̃ is the position vector and || · · · || is the euclidean norm. This is then multiplied with a step function
as given by equation (13), such that the force vanishes beyond some critical time. Alternatively one might
find a solution going straight through a geometry with an angle of π/3.5 using an upward force,

ŷ(0.25 − x̃2 ) , ||r̃|| ≤ 0.5/ sin(π/6)
↑
F =
.
0
, ||r̃|| > 0.5/ sin(π/6)

5 Note

that one has to set A11,inlet = 1 when χ2 < .

14

