close

Вход

Забыли?

вход по аккаунту

?

Influence of the Weienberg number on the stability of Oldroyd kind fluids.

код для вставкиСкачать
ASIA-PACIFIC JOURNAL OF CHEMICAL ENGINEERING
Asia-Pac. J. Chem. Eng. 2010; 5: 657–673
Published online 11 September 2009 in Wiley InterScience
(www.interscience.wiley.com) DOI:10.1002/apj.384
Research Article
Influence of the Weißenberg number on the stability
of Oldroyd kind fluids
N. Scurtu1 * and E. Bänsch2
1
Department of Aerodynamics and Fluid Mechanics, Brandenburg University of Technology Cottbus, Siemens-Halske-Ring 14, 03046 Cottbus,
Germany
2
Department of Applied Mathematics III, Friedrich-Alexander University of Erlangen-Nürnberg, Haberstr. 2, 91058 Erlangen, Germany
ABSTRACT: This paper is concerned with nonlinear rheological fluids of Oldroyd type. We present a (formal) stability
analysis of the corresponding system of equations, showing stability limits on the Weißenberg number in certain cases.
To this end, we proceed in several steps, thus separating the possible sources for instabilities. First, a spectral analysis of
the linearized Oldroyd system is presented. Then, the influence of the βa -term on the stability of the constitutive stress
equation and of the full Oldroyd system is examined. Moreover, because this stability analysis is of formal and linear
nature, we complement it by numerical simulations for the system showing that the upper limit of the Weißenberg
number found by the stability analysis is fairly sharp. We thereby try to shed some light on the high Weißenberg
number problem, that is, the problem why in certain cases there seem to exist no solutions to the Oldroyd problem for
large Weißenberg numbers.  2009 Curtin University of Technology and John Wiley & Sons, Ltd.
KEYWORDS: CFD; non-Newtonian fluids; viscoelastic flow simulation; Oldroyd fluid; stability analysis
INTRODUCTION
The equations governing the flow of viscoelastic nonNewtonian fluids have been intensively studied for
instance because of their industrial applications in,
for example, polymer processing. Besides so called
generalized Newtonian fluids, where the viscosity is
a function of the strain, more sophisticated models
introduce a material constitutive equation for taking into
account memory effects of such polymeric materials.
In 1950, Oldroyd proposed a differential equation for
modeling viscoelastic fluids.[1] During the last few
decades, significant progress has been made in the
development of numerical algorithms for stable and
accurate solution of these type of viscoelastic flow
problems, see Refs [2–7].
For incompressible fluids, the constitutive assumptions in combination with the equations of motion taking into account elastic and memory effects lead to
a highly nonlinear system of partial differential equations. For Oldroyd fluids, the flow is described by the
following dimensionless system of equations of hyperbolic–parabolic–elliptic type:
*Correspondence to: N. Scurtu, Department of Aerodynamics and
Fluid Mechanics, Brandenburg University of Technology Cottbus,
Siemens-Halske-Ring 14, 03046 Cottbus, Germany.
E-mail: scurtu@tu-cottbus.de
 2009 Curtin University of Technology and John Wiley & Sons, Ltd.

∂τ + (u · ∇)τ + β (τ, ∇u)

We
a

∂t



+τ
−
2αD(u) = 0,



 ∂u
Re
+ (u · ∇)u



∂t





 −2(1 − α)divD(u) − div τ + ∇p = f ,
div u = 0
(1)
The unknown fields are the velocity u, the pressure
p, and the symmetric stress tensor τ . In order to
fulfill the frame indifference principle, an objective time
derivative is introduced in the constitutive equation
by adding the so called β-term to the material time
derivative:
βa (T , ∇u) := TW − WT − a(DT + TD)
(2)
where D = D(u) = (∇u + ∇ T u)/2 and W =
W (U ) = (∇u − ∇ T u)/2 are the rate of deformation
tensor and the spin tensor, respectively. The Oldroyd
system is to be solved in a domain ⊂ Rd , d ≥ 2,
supplemented by appropriate initial and boundary conditions. For the velocity field, the boundary condition
may be of, for example, Dirichlet type on the whole
boundary = ∂, whereas for the stress field Dirichlet conditions have to impose on the inflow boundary
part − only.
658
N. SCURTU AND E. BÄNSCH
Asia-Pacific Journal of Chemical Engineering
The Oldroyd system is characterized by four dimensionless parameters, the Reynolds number Re > 0, the
Weißenberg number We > 0, defined as the ratio of a
relaxation time of the polymer over a typical time scale
in the flow, the rate of viscoelastic viscosities α ∈ [0, 1]
and a ∈ [−1, 1]. The slip parameter a was interpreted
by Johnson and Segalman[8] as a measure of the nonaffinity of polymer deformation, that is, the fractional
stretch of the polymeric material with respect to the
stretch of the flow field.
In the limit of We = 0, the system reduces to the
Navier–Stokes equations, corresponding to the purely
viscous fluid, whereas the limit case α = 1 corresponds
to the purely elastic Maxwell fluid. The particular cases
of the Oldroyd fluid with a lower and upper convective
derivatives, respectively, are referred to as Oldroyd-A
(a = −1) and Oldroyd-B (a = 1) fluid, see Ref.[1].
The most critical parameter of the system is the
the Weißenberg number. Although small values of the
Weißenberg (describing slightly elastic liquids) can be
handled well, there are unsolved problems for high
Weißenberg numbers corresponding to a large relaxation
and highly elastic fluids.
The problems associated with such cases of highly
elastic viscoelastic fluids, where We is large, are consequently called the high Weißenberg number problem.
An upper limit for Wecr , beyond which no solution
exists, is found in all published works applying numerical techniques to differential and integral models, see
Refs [2,3,7,9,10].
Minor changes in the constitutive equation and/or the
algorithm employed can lead to higher limit values of
We. As We approaches the critical value Wecr , it is often
observed that spurious oscillations appear in the field
variables; the stress components are then more severely
affected than the velocity components, yielding large
and erroneous stress gradients. The spurious oscillations
have no physical background, and their wavelength
depends on the mesh used for the discretization. Mesh
refinement and corner strategies affect the critical conditions for breakdown, but it is difficult to discern an
overall consistent trend in published works. For a number of benchmark problems in case of steady flows, the
reviews[2,10,11] showed that the limits in the maximum
attainable Weißenberg number still exist, independently
of the numerical method used for solving the Oldroyd
problem.
Concerning the existence of solutions for the Oldroyd
system with suitable boundary conditions, Renardy[12]
has summarized the known results that are basically
four types of existence results:
For the time-dependent case, existence of solutions
locally in time and for small data, globally in time
has been proved in Banach spaces by FernándezCara et al .[13] For the Oldroyd-B model without the
convective term in the constitutive equation, existence
of solution in a convex polygonal domain in the
stationary case and for small relaxation times has
been proved by Picasso and Rappaz.[14] Mathematical
existence of a solution in a Banach space was proved
by Bonito et al.[15] for the time-dependent system
but also without the convective term. For the case
of a = 0, Lions and Masmoudi[16] have proved the
global existence of a weak solution that satisfies energy
identities in the two-dimensional case.
The main aim of this paper is to present a (formal)
stability analysis of the Oldroyd system and thereby
shed some light on the question about the mechanism
that leads to an upper limit of We. Upon performing
a linearized stability analysis, we will show that the
system can become unstable at supercritical Weißenberg
numbers.
Because this analysis is of formal nature, we complement it by numerical simulations. We computationally
investigate different situations, where instabilities show
up, and try to classify the origin of the arising instabilities.
Our findings imply that no accurate numerical scheme
should overshoot the limit of Wecr . Simulations that go
much beyond this limit should be estimated with care.
Although we mainly focus on the Oldroyd fluid
model, our results may be extended to other nonNewtonian fluids like the Giesekus model[17] or the PTT
model.[18]
The rest of this paper is organized as follows. Section
Stability Analysis contains a formal spectral analysis of
the linearized Oldroyd system. For the full nonlinear
system, stability limits on the Weißenberg number
will be presented. In Section Numerical Study, we
perform numerical simulations on three test cases for
the validation of the stability limits found in Section
Stability Analysis. We conclude the paper by a summary
in Section Conclusions.
(1) results on existence locally in time for initial value
problems, when the solution is a small perturbation
of the initial data;
(2) results on global time existence and asymptotic
decay if the initial conditions are small perturbations of the rest state;
In this section, we study stability properties of the
Oldroyd system by linear perturbation technique. We
proceed in several steps. First, we examine the case,
when neglecting all nonlinear terms, like the convective
term and the β-term. Then we analyze the stress
transport equation alone, but with the β-term included.
 2009 Curtin University of Technology and John Wiley & Sons, Ltd.
(3) results on existence of steady flows that are small
perturbation of the rest state;
(4) results on existence of steady flows that are small
perturbation of the Newtonian flow.
STABILITY ANALYSIS
Asia-Pac. J. Chem. Eng. 2010; 5: 657–673
DOI: 10.1002/apj
Asia-Pacific Journal of Chemical Engineering
INFLUENCE OF WEIßENBERG NUMBER
Next, we linearize the full system at a given state
(τ 0 , u 0 , p 0 ), but still without convection. Finally, we
study the full system.
It turns out that this procedure is adequate, as we
identify the possible mechanism leading to instability
separately. Our finding is that the possible source for
instabilities is the β-term, as it is a term of 0-th
order with an indefinite sign. Convection may transport
this instabilities. The limitation of our analysis of
course is that we cannot investigate strong nonlinear
interactions. However, numerical examples in the next
section support our theoretical findings.
Let us mention that without loss of generality we
restrict ourselves to the case of two space dimensions.
Stability analysis for the linearized Oldroyd
system
We start by studying stability properties of the Oldroyd
system (1) linearized at u ≡ 0, that is, by neglecting the
convective terms and the β-term:
Problem 1 Find (τ , u, p) fulfilling
summed up. Multiplying the second equation of Oldyfour by ξj and using the third equality of Oldyfour, the
pressure can be written in the form:
p̂ =
τ (x , t) =
R2
exp (−i ξ · x )τ̂ (ξ, t) dξ,
R2

∂ τ̂mj

ˆ

 We ∂t + τ̂mj + i α(ξm ûj + ξj ûm ) = fs ,
∂ û
ξξ ξ
Re ∂tj +i ξm τ̂mj +(1 − α)|ξ |2 ûj −i j m2 s τ̂sm = ĝj , (7)


|ξ |

ξm ûm = 0
where ĝ = fˆ + i ξ · (ξ · fˆ )/|ξ |2 . For simplicity, without loss of generality, in the sequel we will treat the case
fs = 0 and f = 0 only.
Considering now the vector of unknowns U :=
(τ̂xx , τ̂xy , τ̂yy , ûx , ûy ) , system (7) can be put in the
form:
∂U
+ AU = 0, with
∂t
1
1 B
I3
We
We
A =
2
1 C (1 − α)|ξ | I
2
Re
Re
exp (−i ξ · x )û(ξ, t) dξ,
R2
exp (−i ξ · x )p̂(ξ, t) dξ
(8)
where the submatrices B and C are defined as follows
2ξ1 0
B = i α ξ2
ξ1
0 2ξ2
i
ξ1 ξ22 −ξ2 z −ξ1 ξ22
C =
(9)
2
ξ12 ξ2
|ξ |2 −ξ1 ξ2 ξ1 z
p(x , t) =
(6)
Here, |ξ | is the length of the wave vector and z =
ξ12 − ξ22 .
The eigenvalues of the matrix A are given by:
u(x , t) =
|ξ |2
Thus, by eliminating the pressure field in Eqn 5, one
obtains the following system for the unknowns τ̂ and
û

∂τ

 We ∂t + τ − 2αD(u) = fs ,
∂u
(3)

 Re ∂t − 2(1 − α)divD(u) − div τ + ∇p = f ,
div u = 0
with appropriate initial and boundary conditions.
The right-hand side fs in the stress transport equation
is introduced for the matter of generality. We (formally)
perform a Fourier transform of this system and its
unknowns in R2 :
i ξ · fˆ + ξ · (ξ · τ̂ )
λe4,5
This transformation leads to the following system:

∂ τ̂mj


 We ∂t + τ̂mj + i α(ξm ûj + ξj ûm ) = fˆs ,
(5)
∂ û
Re ∂tj + i ξm τ̂mj + (1 − α)|ξ |2 ûj − i ξj p̂ = fˆj ,



ξm ûm = 0
where fˆs and fˆ are the corresponding Fourier transform of the right-hand sides fs and f , m, j ∈ {1, 2},
respectively, and repeated indices are understood to be
 2009 Curtin University of Technology and John Wiley & Sons, Ltd.
1
(1 − α)|ξ |2
, λe3 =
,
We
Re
√
1
1 (1 − α)|ξ |2
+
± e
=
2
Re
We
λe1,2 =
(4)
(10)
with
e =
1
(1 − α)|ξ |2
−
Re
We
2
−
4α|ξ |2
ReWe
(11)
The first three eigenvalues are clearly real and
positive, so for the stability of the system (8) it
remains to show that the other two eigenvalues λ4 , λ5
Asia-Pac. J. Chem. Eng. 2010; 5: 657–673
DOI: 10.1002/apj
659
660
N. SCURTU AND E. BÄNSCH
Asia-Pacific Journal of Chemical Engineering
have positive real parts. To do so, two cases can be
distinguished. First, we observe that for e < 0, the real
part R(λe4,5 ) = (1 − α)|ξ |2 /Re + 1/We is positive and
thus we have stability in this case. Secondly, if e ≥ 0
then λ4,5 ∈ R and one can easily check that these two
eigenvalues are positive. Thus, the linearized Oldroyd
system is unconditionally stable.
Influence of the β-term on the stability of the
constitutive equation
The next step in examining the stability properties of
the Oldroyd system is to study the influence of the βterm. To this end, we first look at the stress problem
alone.
Problem 2 Given a solenoidal velocity field u, solve
in × R+ the following equation for the stress tensor
τ
2α
1
∂τ
+ (u · ∇)τ + βa (τ, ∇u) +
τ =
D
∂t
We
We
(12)
supplemented by the boundary condition τ|− = τ on
the inflow boundary and the initial condition τ|t =0 = τ0 .
Neglecting the convective term (u · ∇)τ , the stress
equation with given divergence-free velocity field
reduces to a simple initial value problem. By representstress tensor in the compact form τ =
T
ing the symmetric
τxx τxy τyy , the stress equation can be written as:
∂τ
+ Aτ = R, A =
∂t
We −1 − 2aD
2(−aDxy + Wyx )
0
xx
−1
−(aDxy + Wyx )
We
−aDxy + Wyx
0
−2(aDxy + Wyx ) We −1 + 2aDxx
(13)
with R some right-hand side depending on u and
D =
Dxx
Dxy
0
Dxy
,W =
−Dxx
Wyx
−Wyx
0
(14)
Because of the incompressibility condition (1)3 , the
component Dyy of the rate of deformation tensor is equal
to −Dxx .
The three eigenvalues of matrix A depend on the
velocity field and the Weißenberg number in the following way:
λs = We −1 ,
2 + D2 ) − W 2
λ± = We −1 ± 2 a 2 (Dxx
xy
yx
In the case of a = 0, this condition leads to an upper
stability limit for the Weißenberg number:
−1
2 + D2 ) − W 2
We cr = R 2 a 2 (Dxx
xy
yx
(16)
which is caused by the β-term. More specifically, in
the case of a = ±1, this stability bound implies the following: if their exists any point x in the computational
domain, where the velocity field satisfies det(∇u(x )) <
0, then the stress equation can be stable only up to
Wecr = (2 | det(∇u)|)−1
(17)
This limit is confirmed by our numerical tests in Section
on Numerical Study, see also Ref.[19]. At this point, we
remark that the restriction imposed in Ref.[14] on the
relaxation time for the test example used there is nothing
else but the stability restriction det(∇u(x )) > 0 found
in this work. In the case of a = 0, the radicand in Eqn
(15) is always negative, so no stability limit for the
Weißenberg number exists. This result is in agreement
with the conclusion of Jourdain et al.[20] , who have
proven formally exponential convergence to equilibrium
for the case a = 0 of Eqn (12).
Constitutive equation with β-term
and convection
Stability limits can also exist for the full stress constitutive equation with convection. However, this situation
is more involved.
In this case, we have to look at the possible development of instabilities along streamlines, on which
the stress equation can be again treated as an initialboundary problem.
It turns out that along streamlines, which leave the
domain in a finite time, instabilities cannot blow up,
but may rather disturb the solution with at most finite
amplitude. This is because these perturbations have only
a finite time to grow during traversing of streamlines
through a region, where the stability condition is
violated.
However, if there exists a stagnation point of the flow,
where det(∇u) < 0 (in the case of a = ±1), or if closed
streamlines intersect such a region, then instabilities
along these streamlines arise for Weißenberg numbers
beyond the critical value.
(15)
Influence of the nonlinearity β on the stability
of the Oldroyd system
The stable regime is again given by the condition that
the real parts of the eigenvalues of matrix A are positive.
In this subsection, the full Oldroyd system but without
the convective terms is considered:
 2009 Curtin University of Technology and John Wiley & Sons, Ltd.
Asia-Pac. J. Chem. Eng. 2010; 5: 657–673
DOI: 10.1002/apj
Asia-Pacific Journal of Chemical Engineering
INFLUENCE OF WEIßENBERG NUMBER
Problem 3 Solve in × R+ the following system of
equations
∂U
+ A5 U = 0, with
∂t
1 B + B̃ A
We
A5 =
1
1 − α |ξ |2 I
2
Re C
Re

1
2α

 ∂τ
∂t + βa (τ, ∇u) + We τ − Re divD = 0
2(1 − α) divD − 1 div τ + 1 ∇p = f , (18)
∂u

 ∂t −
Re
Re
Re
div u = 0
with boundary condition u| = u , and initial conditions
u|t =0 = u0 , τ|t =0 = τ0 .

The matrices B and C are those defined in Eqn (9),
the matrix A is analogous to Eqn (13), but with D and
W replaced by D(u 0 ) and W (u 0 ), and B̃ is defined by:
0
0
ξ1 − (1 + a)τxy
ξ2
−2aτxx
0
(1 − a)τxy
ξ1
λξ =
u(x , t) = u (x ) + u(x , t),
0
(19)
Neglecting higher order terms, this splitting leads
to the following linear system of equations for the
perturbations:
Problem 4 Solve in × R+ for τ , u and p, fulfilling
the following system of equations:
 ∂τ
1 τ
+ βa (τ , ∇u 0 ) + βa (τ 0 , ∇u) + We


∂t


2α D(u) = 0,
− We

 ∂u − 2(1 − α) divD(u) − 1 div τ + 1 ∇p = 0,

 ∂t
Re
Re
Re
div u = 0
(20)
with boundary condition u | = 0 and initial conditions
u |t =0 = u 0 and τ |t =0 = τ 0 .
We would now like to perform a Fourier transform
again like in the previous section. But this is only
feasible, if τ 0 , ∇u 0 , which act like coefficients for the
unknowns in Eqn (20), would be constant in space,
which is an unrealistic assumption in general. However,
in performing a Fourier transform one might think of
not too oscillatory functions τ 0 , ∇u 0 or local properties.
Thus, consider now formally the Fourier transforms
(τ̂ , û, p̂) of the unknown fields (τ , u, p). Removing
again the pressure from the set of equation like in
Eqn (7), Eqn (20) can be written as:
 2009 Curtin University of Technology and John Wiley & Sons, Ltd.
(22)
The case a = 1
In the case of a = 1, for the Oldroyd-B model, one
eigenvalue of matrix A5 is:
τ (x , t) = τ 0 (x ) + τ (x , t),
p(x , t) = p 0 (x ) + p(x , t)

1 [(1 − a)τ 0 − (1 + a)τ 0 ]ξ 
yy
xx 1
2
0
0
−2aτyy
ξ2 − (1 + a)τxy
ξ1
0
0
B̃ = i  1 [(1 − a)τxx
− (1 + a)τyy
]ξ2
2
0
(1 − a)τxy ξ2
Let (τ 0 , u 0 , p 0 ) be a stationary solution of Problem 3, and (τ , u, p) be a small perturbation of the
stationary solution such that
(21)
1−α 2
|ξ |
Re
(23)
and the other four eigenvalues can be obtained from the
solutions λ of the polynomial equation:
1−α 2
|ξ | − λ
Re
2 2
1
1
0
2
× 4 det (∇u ) +
|ξ | +
−λ
−λ
We
We
1
−λ
We
α
|ξ |2
0
0
0
|ξ |2 − ξ1 2 τxx
+ 2ξ1 ξ2 τxy
+ ξ22 τyy
We
Re
1
+
− λ 2Dxx ξ1 ξ2 − Dxy z + Wyx |ξ |2
We
2
0
0
0
× −ξ1 ξ2 τxx
+ z τxy
+ ξ1 ξ2 τyy
Re
4α
−
[Dxx z + 2Dxy ξ1 ξ2 ]2
ReWe
4
0
[Dxx z + 2Dxy ξ1 ξ2 ] [Dxy |ξ |2 − Wyx z ]τxy
+
Re
0
+ [Dxx ξ1 + (Dxy + Wyx )ξ2 ]ξ1 τxx
0
= 0
(24)
+ [−Dxx ξ2 + (Dxy − Wyx )ξ1 ]ξ2 τyy
×
with z = ξ12 − ξ22 as in Eqn (9). It is interesting to
observe that the left-hand side of this equation comprises only terms of order two and four in |ξ |. For
small wave vectors, the fourth-order terms can be
Asia-Pac. J. Chem. Eng. 2010; 5: 657–673
DOI: 10.1002/apj
661
662
N. SCURTU AND E. BÄNSCH
Asia-Pacific Journal of Chemical Engineering
neglected compared with the second-order terms. Thus,
one obtains a reduced equation
1
−λ
λ
We
4 det (∇u ) +
0
1
−λ
We
2 = 0
(25)
which possess the solutions λ0 = 0 and λs , λ± as in
the pure stress constitutive equation (15). This means,
the stability properties of system (20) for case a = ±1
are, as expected, not better than those of the pure stress
constitutive equation (12).
The case a = 0
In the case of a = 0, the matrix A5 has the two
eigenvalues λξ and λs as above. The other three
eigenvalues are solutions of the following equation,
which is of third order with respect to λ:
2 1
1−α 2
−λ
|ξ | − λ
We
Re
1
1−α 2
2
+ 2Wyx
|ξ | − λ +
−λ
Re
We
α |ξ |2
1
z
− τxy ξ1 ξ2 + (τyy − τxx )
We 2
4 Re
1
2
1
+ Wyx τxy z + (τyy − τxx )ξ1 ξ2
=0
Re
(26)
It seems difficult to give explicit solutions for Eqns
(24) and (26) for general wave vectors. However, by
fixing the direction of the wave vector to, for example,
ξ = (ξ1 , 0), ξ = (0, ξ2 ) or ξ = (ξ1 , ±ξ1 ), one can
simplify the analysis of the eigenvalues for a concrete
problem for which the velocity and stress fields are
known. This will be carried out for the examples in
Section on Numerical Study.
Influence of the convective terms on the
stability of the Oldroyd system
Applying the splitting Eqn (19) to the full Oldroyd
system (1), after neglecting higher order terms, the
following linearized problem for the perturbed fields
is obtained:
Problem 5 Solve in × R+ for τ , u and p the
following system of equations

∂τ + (u 0 · ∇)τ + (u · ∇)τ 0 + β (τ , ∇u 0 )

a

∂t



1
2α
0
 + βa (τ , ∇u) +

We τ − We D(u) = 0,
∂u + (u 0 · ∇)u + (u · ∇)u 0 − 2(1 − α) divD(u)


∂t
Re


1 div τ + 1 ∇p = 0,


−

Re
Re
div u = 0
 2009 Curtin University of Technology and John Wiley & Sons, Ltd.
with boundary condition u | = 0 and initial conditions
u |t =0 = u 0 and τ |t =0 = τ 0 .
Introducing again the Fourier transforms of the perturbed fields, Problem 5 can be written as an initial
value problem of type (21) with the matrix A5 now
replaced by
Â5 = A5 − i (ux0 ξ1 + uy0 ξ2 ) I5 + D5
(27)
where A5 is defined by Eqn (21) and

D5 =
0
0
0
∂τxx
 ∂x
 0
B̂ , B̂ =  ∂τxy

D̂
 ∂x
 0
∂τyy
∂x
 0
∂ux
 ∂x
D̂ =  0
∂uy
∂x
0 
∂τxx
∂y 
0 
∂τxy

,
∂y 
0 
∂τyy
∂y

∂ux0
∂y 
 (28)
∂uy0
∂y
The convective terms (u 0 · ∇)τ and (u 0 · ∇)u contribute to the matrix Â5 only through the purely imaginary term i (ux0 ξ1 + uy0 ξ2 ) times identity I5 and thus
does not influence the real parts of the eigenvalues.
The other terms arising from linearizing the convective terms, namely (u · ∇)τ 0 and (u · ∇)u 0 collected in
D5 , however, may influence the real parts of the eigenvalues. Thus, the real parts of the eigenvalues of matrix
Â5 may differ from the real parts of the eigenvalues of
matrix A5 .
The results of the previous stability analysis are
summarized in Table 1. For the case a = ±1, the
critical Weißenberg number
Wecr is (neglecting fourth
order terms) Wecr = (2 | det(∇u 0 )|)−1 , computed from
Eqn (25). By P we denote the set of points x ∈ that
are either stagnation points for the flow or are located
on streamlines that form closed curves in the domain .
For the case a = 0, the critical Weißenberg number
Wecr0 is defined as the minimum value of We for which
Eqn (26) has a solution with negative real part in some
point of the domain for some wave vector.
By We we include the effect of matrix D5 , which
we consider as a perturbation to the eigenvalues without
this term. For a concrete given example one may
compute the exact eigenvalues of matrix Â5 , though.
NUMERICAL STUDY
In this section, we use numerical simulations to validate
our formal stability analysis of the previous section.
To this end we will define three test examples. The
numerical method we are using was developed in
Asia-Pac. J. Chem. Eng. 2010; 5: 657–673
DOI: 10.1002/apj
Asia-Pacific Journal of Chemical Engineering
INFLUENCE OF WEIßENBERG NUMBER
Table 1. Results of the stability analysis for the different cases.P is the set of points that are either stagnation
points for the flow or are located on streamlines that form closed curves in the domain . We0cr is defined as the
minimum value of We for which Eqn (26) has a solution with negative real part. By We we include the effect of
matrix D5 , which we consider as a perturbation to the eigenvalues without this term .
Linear system
Stress eqn./wo conv.
Stress eqn./w conv.
Oldroyd/wo conv.
Full Oldroyd
a=0
Uncond. stable
Uncond. stable
Uncond. stable
Unstable if
We > Wecr0
Unstable if
We > Wecr0 + We
a = ±1
Uncond. stable
Unstable if
det(∇u 0 (x )) < 0,
x ∈ and
We > Wecr
Unstable if
det(∇u 0 (x )) < 0,
x ∈ P and
We > Wecr
Unstable if
det(∇u 0 (x )) < 0,
x ∈ and
We > Wecr
Unstable if
det(∇u 0 (x )) < 0,
x ∈ P and
We > Wecr + We
Ref.[19]. Here, we just sketch the basic ideas. For
further reference, see the citations above.
The method is based on a splitting strategy via the
fractional step θ -scheme for time discretization. The
latter one goes back to Bristeau et al.[21] . By this
algorithm, the three major numerical difficulties of the
Oldroyd system are decoupled: the treatment of the
nonlinearity in the momentum equation (given by the
velocity transport term), the solenoidal condition, and
the stress transport term in the constitutive equation. By
the operator splitting, one reduces the Oldroyd system
to three considerable simpler subproblems: a linear
self-adjoint Stokes like problem, a nonlinear Burgers
problem, and a transport problem for the stress tensor.
The hyperbolic stress constitutive equation is discretized in space by the discontinuous Galerkin method
with piecewise linear functions. For the velocity and
pressure fields, we have chosen the LBB-stable TaylorHood element. This choice of finite elements is justified
by an inf-sup condition, see Ref. [22]. The algorithm
was implemented via the finite element package Alberta,
see Ref. [23].
Numerical tests were performed for different refinement levels of the underlying mesh starting from a
macro triangulation of the domain consisting of two
triangles. For refinement level n, the mesh width h is
then h = 2−n .
All simulations were run as time evolutions. They
were started with an interpolant of an exact, steady
state solution as initial condition. Then we followed the
time evolution of the discrete solution. In case of stable
solutions, they converged to a steady discrete state not
too far off the exact solution. In case of instabilities, the
solutions departed from this exact solution and blew up
for time t → ∞.
To further assess the quality of the numerical solutions, we introduce the experimental order of convergence (EOC), defined by
EOC :=
1 Errh
,
log 2 Errh/2
Errh = lim vex − vh (·, t)
t→∞
(29)
 2009 Curtin University of Technology and John Wiley & Sons, Ltd.
Here, Errh is the error of an exact steady solution vex and
the numerical solution vh corresponding to a mesh with
width h. v may stand for u, p, or τ . The EOC gives
the computational rate of convergence for the steady
discrete solution to the exact one with respect to the
mesh size h.
For the pressure and velocity fields, we chose the
L2 - and H 1 -norm, respectively, and for the stress
components we considered the L2 -norm. Note that for
the choice of elements used in our computations, that
is, globally continuous, piecewise quadratics for the
velocity, globally continuous, piecewise linears for the
pressure, and discontinuous, piecewise linears for the
stresses, the optimal order is 2 (in case of the H 1 -norm)
and 3 (for the L2 -norm) for the velocity, 1 (in case of
the H 1 -norm) and 2 (for the L2 -norm) for the pressure
and 2 for the stresses.
Example 1
We choose the right-hand sides fs and f in such a way
that the exact velocity, stress components, and pressure
are:


uex = sin( π2 x ) sin( π2 y) ,


π x ) cos( π y) ,
cos(

2
2

 xx , xy, yy
= pex = cos( π2 x ) sin( π2 y)
τex
(30)
The reason for this choice is that the eigenvalues
from Section on Stability Analysis, which can be
computed explicitly here, are subcritical in some parts
of the domain, whereas supercritical in others. The
computational domain is the unit square, = [0, 1]2 .
We will now study the behavior of the computational
solutions in a similar order like for the stability analysis.
We start with the first example, where one might
encounter instability.
Asia-Pac. J. Chem. Eng. 2010; 5: 657–673
DOI: 10.1002/apj
663
664
N. SCURTU AND E. BÄNSCH
Stress equation without convective term for
a=1
In the case of a = 1, the eigenvalue λ− from Eqn (15)
for this example is given by
1
We
π
π
π
π
− π cos2 ( x ) sin2 ( y) − sin2 ( x ) cos2 ( y) (31)
2
2
2
2
λ− =
For increasing values of We, this eigenvalue is the
first passing through zero and this takes first place
in the corner point with coordinates (x , y) = (0, 1) at
the critical Weißenberg number Wecr = 1/π ≈ 0.3183.
Beyond this critical state, the part of the domain, where
the eigenvalue λ− is negative, becomes increasingly
larger.
The respective eigenvalues and eigenvectors corresponding to the A-matrix for this example problem in
the corner point (0, 1) are:
1
1
1
, λ+ =
+ π, λ− =
− π,
We
We
We
e0 = [0 1 0], e+ = [0 0 1], e− = [1 0 0] (32)
λ0 =
Consequently, the eigenvector corresponding to the
eigenvalue λ− is the unit vector in the x direction, that
Asia-Pacific Journal of Chemical Engineering
is, λ− affects only the stress component τxx . Thus, when
λ− becomes negative, then the component τxx of the
stress tensor will first be affected by the instability. The
instability in the other stress components τxy and τyy
appear much later.
This behavior is clearly confirmed by the computational results. Figure 1 shows a perturbation emanating
from point (0, 1), and the part of the domain with
instability is broadening with increasing Weißenberg
number. We note that this phenomenon is essentially
independent of the mesh width.
Stress equation with convection for a = 1
As just seen, the stress equation without convection
and a = 1 admits stable solutions only for subcritical
Weißenberg numbers. For the stress transport equation
with convection, the stress is transported along streamlines that are depicted for the present example in Fig. 2.
The inflow boundary for Example 1 is the line y = 0
and x = 1 is the outflow boundary. All interior points
of the domain lie on streamlines that leave the domain,
thus for these points the stress field is unperturbed, see
Fig. 2, which shows a very different behavior than for
the above case without convection, Section on Stress
Equation without Convective Term for a = 1.
Figure 1. Isolines of the solution τxx of the stress equation without convective term for a = 1 and mesh size h = 2−6 .
Stable solution for We = 0.3 and unstable solutions for We ∈ {0.32, 0.4, 0.5} at time t = 10. For t → ∞, the
numerical solutions for the three last cases blew up. This figure is available in colour online at www.apjChemEng.com.
Figure 2. Velocity field (left) and isolines for the steady solution τxx (right) of the stress equation with convection
for a = 1 and mesh size h = 2−6 at We ∈ {0.4, 0.5, 0.6}. The stress field is only slightly perturbed on the top of the
domain. This figure is available in colour online at www.apjChemEng.com.
 2009 Curtin University of Technology and John Wiley & Sons, Ltd.
Asia-Pac. J. Chem. Eng. 2010; 5: 657–673
DOI: 10.1002/apj
Asia-Pacific Journal of Chemical Engineering
INFLUENCE OF WEIßENBERG NUMBER
Oldroyd system without β-term, without
convection
Just like in the stability analysis from Section on
Stability Analysis for the Linearized Oldroyd System,
for the Oldroyd system without the β-term no stability
limits exists. Without showing corresponding results we
note that the computations confirm convergence to the
steady state solution for any value of the Weißenberg
number.
Oldroyd system with β-term, without
convection for a = 1
In the case of the Oldroyd system without convection
for a = 1, the influence of the β-term is similar to the
stress equation without convection for a = 1 with given
velocity field, see Section on Influence of the β-term on
the Stability of the Constitutive Equation. Because for
the pure stress equation the critical state appears at the
corner (0, 1) first, we take a look at the eigenvalues
and eigenvectors of the A5 -matrix from Eqn (21) of the
Oldroyd system at this point. At the corner point (0, 1)
the respective velocity field, the rate of deformation
tensor, the spin tensor, and the stresses calculated from
Eqn (30) are:
ux = uy = 0, Dxx =
π
,
2
π
Dxy = 0, Wyx = 0, det(∇u) = − ,
2
τxx = τxy = τyy = 1
(33)
One eigenvalue of matrix A5 (21) is λξ (defined in
Eqn (23)) and is always positive. Computing the other
four eigenvalues from Eqn (24) explicitly is not a simple
task in general. Therefore, we make a particular, but
hopefully representative choice of the wave vectors, for
example, ξ = (ξ1 , 0) or ξ = (0, ξ2 ), where ξ1 and ξ2
are arbitrary. For this choice, it is easy to compute the
 2009 Curtin University of Technology and John Wiley & Sons, Ltd.
exact eigenvalues of the perturbed solution at the corner
point (0, 1) corresponding to the above wave vectors.
The eigenvalues determined by Eqn (24) are given by:
1
± π,
We
|ξ |2
1
1
(1 − α)
+
=
2
Re
We
λ1,2 =
λ3,4

!
2
"
2
2
"
|ξ |
1
|ξ | 
1
±# (1 − α)
+
−1
−4
.
Re
We
We
Re
The real parts of the eigenvalues λ1 and λ3 are
always positive, the real part of λ4 is positive as long
as We < 1, but the real part of the eigenvalue λ2
becomes negative already at Wecr = 1/π . The eigenvector corresponding to the eigenvalue λ2 is given by
e2 = [1 , 0 , 0 , 0 , 0]. Thus, similar as in the analysis
of the pure stress equation without the convective term,
the stress component τxx is first affected by instabilities,
when λ2 becomes negative.
h
L2
In Fig. 3, the time evolution of the error τxx − τxx
is shown for two subcritical values of the Weißenberg
number We ∈ {0.3, 0.31} and for the supercritical value
We = 0.32. For subcritical values, the time evolution
converged to the stationary (discrete) state, but for
supercritical values no convergence was achieved and
the computations eventually blew up. The stress component τxx is plotted in Fig. 4 for We = 0.32 along
the diagonal y = 1 − x of the square domain. One can
observe the instability appearing in the region near the
corner and its growth in time.
5
x 10−3
We = 0.32
We = 0.31
We = 0.3
4.5
||τxx − τhxx||L2 − error norm
In contrast to the interior points, we see some
perturbations at the line y = 1. The solution of the
stress transport equation with convection is steady
but perturbed. The perturbation exists only along the
streamline y = 1 in the upper element layer and grows
with growing Weißenberg number.
In the following subsections, we describe the results
of simulations for the Oldroyd system both without and
with β-term, as well as without and with convective
terms. The numerical runs were carried out for different
values of the Weißenberg number, while keeping the
other parameters constant at the values Re = 1.0 and
α = 0.89.
4
3.5
3
2.5
2
1.5
1
0.5
0
0
20
40
60
80
100
t − time
h 2
Figure 3. Time evolution of the error τxx − τxx
L for the
Oldroyd-B system without the stress transport term for
h = 2−6 at different We. This figure is available in colour
online at www.apjChemEng.com.
Asia-Pac. J. Chem. Eng. 2010; 5: 657–673
DOI: 10.1002/apj
665
666
N. SCURTU AND E. BÄNSCH
Asia-Pacific Journal of Chemical Engineering
Full Oldroyd-B system, a = 1
In the case of the full Oldroyd system with a = 1, we
analyze the eigenvalues of matrix Â5 (defined in Eqn
(27)) for the particular choice of the wave vectors of
the perturbed solution, ξ = (ξ1 , 0) and ξ = (0, ξ2 ), as
before.
For the wave vector ξ = (ξ1 , 0), the eigenvalues of
matrix Â5 (defined in Eqn (27)) are:
instable regime for each We, whereas for large wave
vectors the stable regime is determined by the condition
$
%
1
π Re
π
2
, 1−
|ξ | > Re and We < min
2
π
2 |ξ |2
(37)
For the wave vector ξ = (0, ξ2 ) of the perturbed
solution, the eigenvalues of matrix Â5 are given by
1−α 2 π
1
|ξ | + , λ1,2 =
± π,
Re
2
We


!
"
2
2
2
2
"
|ξ |
|ξ |
1
π
1
π
|ξ |
π 1 
1
1
+
−
± # (1 − α)
+
−
−1
−
=  (1 − α)
−4
 (34)
2
Re
We
2
Re
We
2
We
Re
2 We
λξ =
λ3,4
1−α 2 π
|ξ | − , λ1,2 =
Re
2

!
"
2
2
"
|ξ |2
|ξ |2
1
1
1
π
1
π
|ξ |
#
(1 − α)
=  (1 − α)
−4
+
+
±
+
+
−1
+
2
Re
We
2
Re
We
2
We
Re
λξ =
λ3,4
1
± π,
We

π 1 
(38)
2 We
The real parts of the eigenvalues λξ , λ1 and λ3 are
always positive, while the real part of λ2 becomes
negative already at We = 1/π , whereas the real part
of the eigenvalue λ4 is positive as long as
2
|ξ |
π 1
1
−1
−
>0
(35)
We
Re
2 We
Again, the eigenvalue λξ is positive only for large
wave vectors satisfying |ξ |2 > π Re/2(1 − α). The real
part of the eigenvalues λ1 and λ3 are always positive,
the real part of λ2 becomes negative already at We =
1/π , but the real part of the eigenvalue λ4 is positive
as long as:
For small wave vectors, if |ξ |2 < Re π/2, the real part
of the eigenvalue λ4 is negative for each value of We,
which interestingly means that we have a long wave
instability. We note that as in our example Re = 1 and
thus Re π/2 ≈ 1.57 and = [0, 1]2 ,
these long wave instabilities do not show up in our
computations.
For large wave vectors, if |ξ |2 > Re π/2, the real part
of the eigenvalue λ4 is positive only for
(39)
We < 1 −
π Re
2 |ξ |2
(36)
Thus, the eigenvalue analysis corresponding to the
wave vector ξ = (ξ1 , 0) shows that for small wave
vectors, |ξ |2 < Re π/2, the Oldroyd-B system is in the
 2009 Curtin University of Technology and John Wiley & Sons, Ltd.
2
|ξ |
π 1
1
−1
+
>0
We
Re
2 We
and accordingly for values of the Weißenberg number
that satisfy
π Re
(40)
We < 1 +
2 |ξ |2
Therefore, the eigenvalue analysis corresponding to
the wave vector ξ = (0, ξ2 ) yields that the stable regime
is given by the condition
$
%
1
π Re
π Re
and We < min
, 1+
|ξ |2 >
2 1−α
π
2 |ξ |2
(41)
From conditions (37) and (41), one can conclude that
the solutions are stable only for large wave vectors and
Asia-Pac. J. Chem. Eng. 2010; 5: 657–673
DOI: 10.1002/apj
Asia-Pacific Journal of Chemical Engineering
INFLUENCE OF WEIßENBERG NUMBER
1.2
4
t = 50
t = 80
exact solution
0.8
0.6
error norm
τxx − stress component
1
0.4
0.2
3.5
||τxx − τhxx||L2
3
||τxy − τhxy||L2
2.5
||τyy − τhyy||L2
||u − uh||L2
2
1.5
1
0
−0.2
x 10-3
0.5
0
0.2
0.4
0.6
0.8
1
0
position along y = 1x
0
5
10
15
20
25
t−time
Numerical values of τxx for the Oldroyd-B
system without the stress transport term along the line
y = 1 − x at We = 0.32. This figure is available in colour
online at www.apjChemEng.com.
Figure 4.
More precisely, for small wave vectors Wecr = 0 and
for large wave vectors Wecr = 1/π if α ≥ 1/π and
Wecr = α if α < 1/π .
Similarly as in the case of the pure stress equation
with given velocity field, the corner point (0, 1) may
cause instability, which is propagated along the streamline traversing this point. The numerical simulations
show stable steady solutions up to We ≈ 0.4, which is
slightly above the critical value Wecr . Beyond this value,
however, for instance for We ∈ {0.5, 0.6}, the perturbation in the τxx -component, emanating from the corner
point, is growing along the line y = 1 and eventually
blowing up.
Because the Oldroyd system is a nonlinear coupled
system for the unknown fields, the instability in τxx
eventually also affects the other components of the
solution, which can be seen in Fig. 5. Here, the growth
of the error norms in time is displayed for We = 0.5
and We = 0.6 for the mesh with size h = 2−6 and time
step size t = 0.01.
The stress component τxx is plotted in Fig. 6 at
different time instants for We = 0.6 along the line
y = 1 − h/2. One can see the deviation of the numerical
values of τxx from the exact solution near the streamline
passing the stagnation point.
Neither further mesh refining nor reduction of the
time step led to better convergence results, confirming
that this phenomenon is inherent in the equation.
To show that the perturbation in the upper element
layer is only caused by the instability arising at the
 2009 Curtin University of Technology and John Wiley & Sons, Ltd.
Oldroyd-B system for h = 2−6 at We = 0.6. This figure is
available in colour online at www.apjChemEng.com.
1
0.5
τxx − stress component
small Weißenberg numbers, which satisfy

 We < 1 ifα ≥ 1
π Re
π
π
and
(42)
|ξ |2 >
 We < 1 − π2 Re2 ifα < π1
2 1−α
|ξ |
Figure 5. Time evolution of the error L2 -norms for the full
0
−0.5
−1
−1.5
−2
t = 15
t = 25
exact solution
−2.5
−3
0
0.2
0.4
0.6
0.8
1
x − position along y = 0.9961
Numerical values of τxx for the OldroydB system along the line y = 1 − h/2 for h = 2−7 at
We = 0.6. This figure is available in colour online at
www.apjChemEng.com.
Figure 6.
stagnation point (x , y) = (0, 1), and is not influenced
by some boundary effect, calculations were performed
also in the domain = [−1, 1] × [0, 2]. In Fig. 7,
the velocity field and the time evolution of the stress
component τxx are presented for We = 0.6 for a mesh
with size h = 2−5 and time step width t = 0.01.
Stress transport problem for a = 0
Now, let us clarify what happens in the case of a = 0.
For the pure stress transport, the eigenvalues (15) are in
the stable regime for each value of the Weißenberg number. Consequently, no stability limit exists in this case.
Asia-Pac. J. Chem. Eng. 2010; 5: 657–673
DOI: 10.1002/apj
667
668
N. SCURTU AND E. BÄNSCH
Asia-Pacific Journal of Chemical Engineering
τxx at t = 30
u at t = 100
τxx at t = 100
Figure 7. Velocity field and propagation of the τxx perturbation at We = 0.6 for the
full Oldroyd-B system; = [−1, 1] × [0, 2]. This figure is available in colour online at
www.apjChemEng.com.
Table 2. Errors and EOCs for Example 1 for a = 0 and We = 1000; ‘Ref.’ is the refinement level of the mesh.
xy
yy
Ref.
τ xx − τhxx L2
EOC
τ xy − τh L2
EOC
τ yy − τh L2
EOC
4
5
6
7
6.272069e-04
1.464079e-04
3.497572e-05
8.504148e-06
2.10
2.07
2.04
–
5.865086e-04
1.410928e-04
3.429843e-05
8.415914e-06
2.06
2.04
2.03
–
5.677649e-04
1.382838e-04
3.388190e-05
8.354178e-06
2.04
2.03
2.02
–
This fact is reflected in Table 2, which contains EOC
tests for We = 1000. The computations converged to a
steady state for any value of the Weißenberg number.
Oldroyd system without convection for a = 0
Although the stress transport equation for a = 0 is
unconditionally stable, in the Oldroyd system without
convection instabilities may arise, as it can be seen
from the eigenvalue analysis of matrix A5 defined in
Eqn (21). Like above, we simplify the explicit analysis
by fixing the wave vector to, for example, ξ = (ξ1 , 0)
and ξ = (0, ξ2 ), where the components ξ1 and ξ2 are
arbitrary.
At the corner (x , y) = (0, 1), where relations (33)
hold, Eqn (26) corresponding to the wave vector ξ =
(ξ1 , 0) or ξ = (0, ξ2 ), admits a solution equal to 1/We
and the other two solutions satisfy
1−α 2
1
1 |ξ |2
2
+
|ξ | +
= 0 (43)
λ −λ
We
Re
We Re
The two solutions of Eqn (43)
&
λ1,2 = b ±
1
b=
2
1 |ξ |2
,
b2 −
We Re
1
1−α 2
+
|ξ |
We
Re
For the wave vectors ξ = (ξ1 , ±ξ1 ), Eqn (26) admits
a solution equal to 1/We and the other two solutions
satisfy the quadratic equation
2
1−α 2
|ξ |
1
1
+
|ξ | +
∓1
= 0
λ −λ
We
Re
We
Re
(45)
with roots
2
&
2
|ξ |
1
∓1
,
=b+
−
We
Re
&
2
|ξ |
1
±
2
λ2 = b − b −
∓1
We
Re
λ±
1
b2
(46)
−
The solutions λ−
1 and λ2 corresponding to the wave
vector ξ = (ξ1 , −ξ1 ) are always positive, but those
corresponding to the wave vector ξ = (ξ1 , ξ1 ), namely,
+
λ+
1 and λ2 , are both positive only as long as We ≤ 1.0.
The eigenvector of matrix A5 corresponding to the
eigenvalue λ+
2 from Eqn (46) for the wave vector
ξ = (ξ1 , ξ1 ), is given by
−1
1
α
+
−1
− λ2
= 2i ξ1
0
We
We
−1
1
α
+
−1
− λ2
−1 1
−2i ξ1
We
We
e2+
(44)
(47)
are always positive.
 2009 Curtin University of Technology and John Wiley & Sons, Ltd.
Asia-Pac. J. Chem. Eng. 2010; 5: 657–673
DOI: 10.1002/apj
Asia-Pacific Journal of Chemical Engineering
Thus, the perturbation will show up first in the stress
componentsτxx and τyy .
By means of the previous analytical considerations
one can explain the instabilities arising in the Oldroyd
system, without and with the convective terms, for
a = 0, when We is increasing. In the case of the
Oldroyd system without convective terms, convergence
to a steady state was obtained computationally up to
Wecr0 = 1.0. At We = 1.1 convergence was achieved
only for large mesh sizes, for example, h ∈ {2−3 , 2−4 },
but not for finer meshes. For supercritical Weißenberg
numbers, the perturbation appeared first in the stress
components τxx and τyy , and slightly later in the
τxy component, which is in agreement with the form
of the eigenvector e2+ in Eqn (47). In Fig. 8, the
growth of the perturbation in the stress component
τxx at We = 1.1 is shown at different times. The perturbation is propagated in all directions of the computational domain, emanating from the corner point
(x , y) = (0, 1).
INFLUENCE OF WEIßENBERG NUMBER
Full Oldroyd system for a = 0
Let us turn to the case of the full Oldroyd system with
a = 0 and Â5 defined by Eqn (27). The eigenvalue
analysis corresponding to the wave vectors ξ = (ξ1 , 0)
and ξ = (0, ξ2 ) shows that only for large wave vectors
satisfying |ξ |2 > π Re/2 the Oldroyd system is in the
stable regime for each We. For the wave vectors ξ =
(ξ1 , ±ξ1 ), two eigenvalues of matrix Â5 are λs (15),
and the other three satisfy the cubic equation
1−α 2 π
1
−λ
|ξ | + − λ
We
Re
2
1−α 2 π
1−α 2
|ξ | − − λ +
|ξ | − λ
Re
2
Re
|ξ |2
α
∓1
=0
(48)
We
Re
The numerical simulations show convergence up to
We ≈ 1.4. Likewise, as for the full Oldroyd-B problem,
Figure 8. Example 1, propagation of the perturbation in τxx for the
Oldroyd system without the convective terms for a = 0 and We = 1.1.
This figure is available in colour online at www.apjChemEng.com.
Figure 9. Example 1, propagation of the perturbation in the pressure for
the full Oldroyd system for a = 0 and We = 1.5, h = 2−6 . This figure is
available in colour online at www.apjChemEng.com.
 2009 Curtin University of Technology and John Wiley & Sons, Ltd.
Asia-Pac. J. Chem. Eng. 2010; 5: 657–673
DOI: 10.1002/apj
669
N. SCURTU AND E. BÄNSCH
Asia-Pacific Journal of Chemical Engineering
200
−5
−6
180
−6 −6
−5
−3
−7
−7
−5
0
0
−3
−5
−6
−6
−5
−7
−7
−3
160
140
−6
−3
−5
−6
0
6
80
10
8
6
4
2
0
−2
We = 0.1
We = 0.109
We = 0.125
We = 0.2
We = 0.5
0
−
− 3
−65
40
As a second example for the full Oldroyd system, we
define = [−1, 1]2 and the following exact solution:

u = (4 (2 − r 2 ) y , −4 (2 − r 2 ) x ),

 ex 2
r = x + y 2,
(49)
xx , xy, yy

= 1 + exp(−2) − exp(−r 2 ),
 τex
pex = exp(−r 2 )
−8
0
7
Example 2
−6
−6
−3
100
0
−5
0
120
60
−4
−7
−6
−7
−5
in the present case, because of the convective terms,
the velocity and stress fields are transported along
the streamlines. Thus again, along the streamlines that
leave the computational domain no instabilities arise
but on the streamline containing the stagnation point,
the negative eigenvalue λ+
2 defined in Eqn (46) and
corresponding to the wave vector ξ = (ξ1 , ξ1 ), gives rise
to instabilities. In Fig. 9, the growth of the perturbation
in the pressure field at We = 1.5 is shown at different
times. It is obvious that the perturbation appears first
in the upper element layer on the boundary streamline
downstream the stagnation point.
real(λ−)
670
0.2
0.4
0.6
0.8
1
1.2
1.4
radial position
Figure 10. Example 2, eigenvalue λ− on the line y = x,
x ≥ 0 at different We. This figure is available in colour online
at www.apjChemEng.com.
20 −
6
−5
−7
20
40
−3
0
−7
60
80
0
−3
−5
6
6
100 120
−5 −6
−7
140 160
−7
−6−5
180 200
Example 2, isolines of Rλ− at We =
0.5. This figure is available in colour online at
www.apjChemEng.com.
Figure 11.
Again the right-hand sides were chosen in such a way
that the above functions are solutions of the Oldroyd–B
system. One key feature of this solution is the fact that
the streamlines consist of circular arcs or closed circles.
The critical Weißenberg number defined by relation
(16) and corresponding to the velocity field (49) is
approximatively Wecr ≈ 0.109. Figure 10 shows the real
part of the eigenvalue λ− on the diagonal line y = x and
x ≥ 0 for different values of We. For We < 0.125 the
points, where Rλ− < 0, are located only on streamlines
that leave the domain. Therefore, only for We > 0.125,
one can expect instabilities to appear. This critical value
of the Weißenberg number is available for the stress
transport equation with a = 1. However, for the full
Oldroyd-B system, like as for Example 1, the stability
limits can be different.
Simulations of the Oldroyd-B system for this example
still converged at We = 0.3, and gave good EOC values,
as it is shown in Table 3.
Table 3. EOC tests for Example 2 applied to the Oldroyd-B system for We = 0.3, t = 0.001.
Ref.
2
3
4
5
Ref.
2
3
4
5
p − ph H 1
2.530e-01
1.149e-01
5.720e-02
2.841e-02
τ xx − τhxx L2
1.355e-02
2.625e-03
6.107e-04
1.531e-04
EOC
1.13
1.00
1.00
–
EOC
2.36
2.10
1.99
–
u − uh L2
9.840e-03
1.211e-03
1.562e-04
1.886e-05
xy
τ xy − τh L2
1.240e-02
2.359e-03
5.562e-04
1.450e-04
 2009 Curtin University of Technology and John Wiley & Sons, Ltd.
EOC
3.02
2.95
3.04
–
EOC
2.39
2.08
1.93
–
u − uh H 1
2.296e-01
5.476e-02
1.402e-02
3.356e-03
yy
τ yy − τh L2
1.427e-02
2.484e-03
5.652e-04
1.439e-04
EOC
2.06
1.96
2.06
–
EOC
2.52
2.13
1.97
–
Asia-Pac. J. Chem. Eng. 2010; 5: 657–673
DOI: 10.1002/apj
Asia-Pacific Journal of Chemical Engineering
INFLUENCE OF WEIßENBERG NUMBER
t = 9.0
u
t = 11.5
Example 2: velocity field and time propagation of the τxx perturbation in
the full Oldroyd-B system for We = 0.5. This figure is available in colour online at
www.apjChemEng.com.
Figure 12.
τxx
u
τxy
Figure 13. Example 3: velocity field and isolines of the stress components τxx and τxy at
We = 100. This figure is available in colour online at www.apjChemEng.com.
Table 4. EOC tests for Example 3 for the full Oldroyd-B system for We = 100 and t = 10−3 .
Ref
3
4
5
6
7
Ref
3
4
5
6
7
p − ph H 1
6.836e+00
2.877e+00
1.020e+00
3.932e-01
1.709e-01
τ xx − τhxx L2
7.378e-01
1.078e-01
2.222e-02
4.158e-03
9.685e-04
EOC
1.25
1.49
1.37
1.20
–
EOC
2.77
2.28
2.41
2.10
–
u − uh L2
6.476e-02
1.009e-02
1.882e-03
2.738e-04
3.444e-05
xy
τ xy − τh L2
4.268e-01
9.977e-02
1.773e-02
4.012e-03
1.009e-03
By increasing the Weißenberg number to 0.5, the
system does not converge any more. In Fig. 12, one
can see the time evolution of the perturbation in the
stress component τxx on mesh with size h = 2−5 .
The perturbation grows and is propagated along the
streamlines, in the region of the negative eigenvalue
(compare with Fig. 11), where the streamlines are closed
curves. On the streamlines that leave the computational
domain, no perturbation appear, even though some of
them traverse the region of unstable regime. But along
the closed streamlines that lie in the unstable region,
the perturbation is propagated and grows in time.
 2009 Curtin University of Technology and John Wiley & Sons, Ltd.
EOC
2.68
2.42
2.78
2.99
–
EOC
2.09
2.49
2.14
2.04
–
u − uh H 1
1.240e+00
4.077e-01
1.071e-01
2.620e-02
6.711e-03
yy
τ yy − τh L2
6.474e-01
1.507e-01
2.148e-02
4.316e-03
9.730e-04
EOC
1.61
1.93
2.03
1.96
–
EOC
2.10
2.81
2.32
2.09
–
Example 3
This third example is taken from Ervin and Ntasin.[24]
The computational domain is an L-shaped domain
given by = [−1, 1]2 \ [0, 1]2 . The velocity, stress and
pressure fields are defined by

y − 0.1 x − 0.1


u=
−
,



r
r
r = (x − 0.1)2 + (y − 0.1)2 ,



 τ = 2 α D(u),

p = (2 − x − y)2
(50)
Asia-Pac. J. Chem. Eng. 2010; 5: 657–673
DOI: 10.1002/apj
671
672
N. SCURTU AND E. BÄNSCH
The velocity and stress fields have singularities just
exterior to the domain near the point (0, 0).
The crucial feature of this example is the fact that
the velocity field satisfies det(∇u) = 0 pointwise everywhere. According to Eqn (15), the real part of the
eigenvalue λ− , corresponding to the stress equation, is
therefore positive in the whole computational domain,
and thus no stability limit for the Weißenberg number exists. Moreover, the streamlines neither comprise
a stagnation point nor are closed curves in the whole
domain . Hence, theoretically, one expects no instabilities for the Oldroyd problem.
Indeed, the simulations converged to the steady state
for every value of the Weißenberg number we tested, see
Fig. 13. In Table 4, corresponding EOCs are presented
for the Oldroyd-B model for We = 100.
Asia-Pacific Journal of Chemical Engineering
for arbitrary values of the Weißenberg number (for
a = ±1).
The eigenvalue analysis showed that even in the case
of a = 0 instabilities can occur. The numerical tests
confirmed this stability limits.
For the full Oldroyd-B system, at least those stability
limits arising in the pure stress constitutive equation,
exist, if the flow domain comprises a stagnation point in
the region where det(∇u) < 0, or if the streamlines are
closed curves that intersect such a region. The numerical
simulations confirmed this stability limit derived by the
eigenvalue analysis in principle, however with a slightly
larger value of the critical Weißenberg number.
The numerical examples given here are only relevant
test cases. Application of the obtained results to real
flows on classical benchmarks, such as the flow in a
4:1 contraction, will be presented in a next paper.
CONCLUSIONS
REFERENCES
In this paper, we studied the stability of the Oldroyd
system by both a formal linear spectral analysis and
numerical simulations.
Our findings are summarized as follows. The pure
stress equation without the β-term as well as the
linearized Oldroyd system are both unconditionally
stable. The possible source for instability is the β-term,
because it is a zero-order term with indefinite sign.
By linear stability analysis, we found conditions
to characterize the region of stability, both for the
pure stress equation and for the full Oldroyd system.
These stability limits were fairly well confirmed by the
simulations.
Considering the pure stress equation without the
stress transport term and for a given stationary velocity
field, in the case of a = 0, a stability limit can exist.
More precisely, if their exists a point or a region in the
computational domain where the velocity field fulfills
det(∇u) < 0, then for√a = ±1 the stress equation is
stable until Wecr = (2 − det(∇u))−1 . In the case of
a = 0, no stability limit was found.
Considering the stress constitutive equation with
given (stationary) velocity field, the stress is transported
along the streamlines. Along streamlines, which leave
the computational domain, no instability arises. However, if there exists a stagnation point of the flow in the
region where det(∇u) < 0 (in the case a = ±1), or if
the streamlines are closed curves that intersect such a
region, then instability of the stress components along
the streamlines arises for Weißenberg numbers above
the critical value. For a = 0 the pure stress equation
is not affected, but in the full Oldroyd system stability
limits exist.
For the Oldroyd system without the stress transport
term, in the case a = 0, one detects the same stability limit Wecr . However, if the velocity field fulfills
det(∇u) ≥ 0 everywhere, the system remains stable
 2009 Curtin University of Technology and John Wiley & Sons, Ltd.
[1] J.G. Oldroyd. Proc. R. Soc. London, 1950; A 200, 523–541.
[2] F.P.T. Baaijens. J. Non-Newtonian Fluid Mech., 1998; 79,
361–385.
[3] C. Bodart, M.J. Crochet. Theoret. Comput. Fluid Dynamics,
1993; 5, 57–75.
[4] R. Guénette, M. Fortin. J. Non-Newtonian Fluid Mech., 1995;
60, 27–52.
[5] D. Sandri. Comput. Meth. Appl. Mech. Eng., 2002; 191,
5045–5065.
[6] P. Saramito. Math. Modell. Num. Anal., 1994; 28, 1–34.
[7] P. Wapperon, M.F. Webster. Comput. Meth. Appl. Mech. Eng.,
1999; 180(3-4), 281–304.
[8] M. Johnson, D. Segalman. J. Non-Newtonian Fluid Mech.,
1977; 2, 255–270.
[9] M.A. Alves, P.J. Oliveira, F.T. Pinho. J. Non-Newtonian Fluid
Mech., 2003; 110, 45–75.
[10] R. Keunings. J. Non-Newtonian Fluid Mech., 1986; 20,
209–226.
[11] B. Caswell. J. Non-Newtonian Fluid Mech., 1996; 62, 99–110.
[12] M. Renardy. Mathematical Analysis of Viscoelastic Flows,
CBMS-NSF Regional Conference Series in Applied Mathematics, SIAM 2000.
[13] E. Fernández-Cara, F. Guillén, R.R. Ortega. Mathematical
modeling and analysis of viscoelastic fluids of the Oldroyd
kind, Handbook of Numerical Analysis, Vol VIII, NorthHolland, Amsterdam 2002; pp.543–661.
[14] M. Picasso, J. Rappaz. ESIAM: Math. Modell. Num. Anal.,
2001; 35(5), 879–897.
[15] A. Bonito, P. Clément, M. Picasso. Numer. Math., 2007; 107,
213–255.
[16] P.L. Lions, N. Masmoudi. China Ann. Math., 2000; 21B(2),
131–146.
[17] H. Giesekus. Phänomenologische Rheologie, Springer, Berlin,
Heidelberg 1994.
[18] N. Phan-Thien, R.I. Tanner. J. Non-Newtonian Fluid Mech.,
1977; 2, 353–365.
[19] N. Scurtu. Stability Analysis and Numerical Simulation of NonNewtonian Fluids of Oldroyd Kind. PhD Thesis. FriedrichAlexander-Universität Erlangen-Nürnberg 2005.
[20] B. Jourdain, C. Le Bris, T. Leliévre, F. Otto. J. Arch. Rational
Mech. Anal., 2006; 181, 97–148.
[21] M.O. Bristeau, R. Glowinski, J. Periaux. Comput. Phys. Rep.,
1987; 6, 73–187.
[22] M. Fortin, A. Fortin. J. Non-Newtonian Fluid Mech., 1989; 32,
295–310.
[23] A. Schmidt, K.G. Siebert.Design of adaptive finite element
software: The finite element toolbox ALBERTA. Springer
LNCSE Series 42 2005.
Asia-Pac. J. Chem. Eng. 2010; 5: 657–673
DOI: 10.1002/apj
Asia-Pacific Journal of Chemical Engineering
[24] V.J. Ervin, L.N. Ntasin. SIAM Numer. Meth. Diff. Equations,
2004; 21, 297–322.
[25] E. Bänsch. Acta Math. Univ. Comenianae, 1998; LXVII,
101–114.
 2009 Curtin University of Technology and John Wiley & Sons, Ltd.
INFLUENCE OF WEIßENBERG NUMBER
[26] R.B. Bird, R.C. Armstrong, O. Hassager. Dynamics of
Polymeric Liquids, John Wiley and Sons, New York 1977.
Asia-Pac. J. Chem. Eng. 2010; 5: 657–673
DOI: 10.1002/apj
673
Документ
Категория
Без категории
Просмотров
1
Размер файла
1 429 Кб
Теги
kinds, number, influence, oldroyd, stability, fluid, weienberg
1/--страниц
Пожаловаться на содержимое документа