SLIDE 1 Cracking the high-Weissenberg number problem
Martien A. Hulsen, Frank P.T. Baaijens Materials Technology, Eindhoven University of Technology, 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
Fourteenth International Congress on Rheology, 22–27 August, 2004, Seoul, Korea
SLIDE 2
Introduction (HWNP) (1)
Observation (last 25 years): All numerical methods break down when the Weissenberg number exceeds a critical value. The precise critical value varies with the viscoelastic model, numerical method and mesh used. ⇒ High Weissenberg number problem (HWNP). ⊲ It is believed to have a numerical origin. ⊲ It is still there, despite: – realistic constitutive models (pom-pom, XPP, FENE, . . . ) – advances in ‘stable’ computational methods (DEVSS, SUPG, DG, . . . )
SLIDE 3
Introduction (HWNP) (2)
For example: L = 30R H = 2R R y x Oldroyd-B: break down at Wi = 0.88 Giesekus: break down at Wi = 1.2
SLIDE 4
Introduction (HWNP) (3)
The situation is quite different than classical CFD at large Reynolds numbers, where: ⊲ stable calculations can be performed, but ⊲ accuracy is lost due to insufficient mesh resolution. Research efforts are directed towards fast and accurate numerical schemes on big computers to solve problems for a better understanding of basic phenomena, such turbulence. ⇒ The situation in Computational Rheology seems hopeless . . .
SLIDE 5
Introduction (HWNP) (4)
. . . until a couple of months ago, that is: Raanan Fattal & Raz Kupferman, JNNFM, in press. where: ⊲ the mechanism for the HWNP has been identified ⊲ a solution has been proposed This talk: a more ‘mechanical’ derivation, application to FEM, convince you that the HWNP has been solved and what’s left to do!
SLIDE 6 Basic equations
Momentum balance and mass balance: ∇p − ∇ · (2ηsD) − ∇ · τ = 0 ∇ · u = 0 with D = 1
2(L + LT) and L = (∇u)T.
Constitutive equation: ˙ c = L · c + c · LT + f(c), τ = η λ(c − I) with f(c) = 1 λ
1 λ
Giesekus
SLIDE 7 Numerical methods
Eulerian setting: ˙ c = ∂c ∂t + u · ∇c FEM: quads with u ∈ (Q2)2, p ∈ P d
1
DEVSS: projected D ∈ (P 1)3 DG: c ∈ (Q1)3 u p c Q2P d
1 Q1
Explicit time integration: Euler forward (1st order) or Adams-Bashforth (2nd order).
SLIDE 8
Analysis of the HWNP in 1D: a numerical instability (1)
Toy problem: find c(x, t) such that ∂c ∂t + a∂c ∂x = bc, x ∈(0, L), t ≥ 0 with ⊲ a > 0, b > 0 ⊲ c(x = 0, t) = c(x, t = 0) = 1 1 exp(bx a ) exp(bt) x = at L c ⇒ exponential growth balanced by convection
SLIDE 9 Analysis of the HWNP in 1D: a numerical instability (2)
Numerical discretization in space (first-order upwind): dci dt = −aci − ci−1 ∆x + bci = (b − a ∆x)ci + a ∆xci−1 ci−2 ci−1 ci ci+1
i − 2 i − 1 i + 1 i
∆x
dc dt = Ac with A = γ β γ ... ... β γ , γ = b − a ∆x, β = a ∆x
SLIDE 10 Analysis of the HWNP in 1D: a numerical instability (3)
To avoid exponential growth in time: γ = b − a ∆x ≤ 0
∆x ≤ a b Define C: C = ∆xb a ,
C = ∆x max(b, 0) |a| then for stability: C ≤ 1 first-order upwind For DG using linear polynomials: C ≤ 2 DG with P1 ⇒ Failure to balance exponential growth with convection when criterion violated.
SLIDE 11
Analysis of the HWNP in 1D: a numerical instability (4)
For Oldroyd-B in extension: a = u, b = 2˙ ǫ − 1 λ, with ˙ ǫ = ∂u ∂x: C = ∆x max(2˙ ǫ − 1 λ, 0) |u| For Giesekus in extension: a = u, b = 2˙ ǫ − 1 λ(1 + 2αcxx), C = ∆x max(2˙ ǫ − 1 λ(1 + 2αcxx), 0) |u|
SLIDE 12
Solution to the problem: log transformation
Solution: solve for s = log c. PDE for s: ˙ s = ˙ c c = bc c = b with ˙ s = ∂s ∂t + a∂s ∂x, and thus: ∂s ∂t + a∂s ∂x = b, x ∈(0, L), t ≥ 0 and c = exp(s). 1 x = at L bt bx a s ⇒ No numerical instability. No restrictions on ∆x
SLIDE 13
Log transformation for tensors
log τij?: No log τ?: No log cij?: No log c: YES! Thus we choose s = log c. BONUS: c = exp(s) always positive definite!.
SLIDE 14 Evolution equation for s = log c (1)
Spectral decomposition: c = c1n1n1 + c2n2n2 + c3n3n3 =
3
cinini, c = c1 c2 c3 s = log c =
3
log(ci) si nini =
3
sinini, n2 n3 n1 s = log c1 log c2 log c3 = s1 s2 s3
SLIDE 15 Evolution equation for s = log c (2)
Substantial derivative of s: ˙ s =
3
˙ sinini +
3
si ˙ nini +
3
sini ˙ ni
si = d(log ci) dt = ˙ ci ci ˙ s =
3
˙ ci ci nini +
3
si ˙ nini +
3
sini ˙ ni ⇒ We need expressions for ˙ ci and ˙ ni
SLIDE 16 Evolution equation for s = log c (3)
We get ˙ ci and ˙ ni from the CE: ˙ c = L · c + c · LT + f(c) From the diagonal components (in the ni-frame): ˙ ci = 2ciLii + fi(c1, c2, c3) From the off-diagonal components (in the ni-frame): ˙ ni = ω · ni =
3
ωjinj ˙ n2 ˙ n3 ˙ n1 n2 n3 n1 with the components ωij of the skewsymmetric tensor ω: ωij = ciLji + cjLij cj − ci , i = j, ci = cj
SLIDE 17 Evolution equation for s = log c (4)
End result: ˙ s = 2
3
Liinini +
3
fi ci nini
+
3
3
i=j
si − sj ci − cj (cjLij + ciLji)ninj
Notes: ⊲ ci = exp(si), c = exp(s) ⊲ lim
ci→cj
si − sj ci − cj (cjLij + ciLji) = Lij + Lji = 2Dij
SLIDE 18
Flow past a cylinder confined between two flat plates
⊲ Oldroyd-B model with ηs/η = 0.59, where η = ηs + ηp ⊲ Giesekus model with α = 0.01. ⊲ Weissenberg number: Wi = λU R , with U average velocity in channel. ⊲ Dimensionless drag force K = Fx ηU L = 30R H = 2R R y x ⊲ Dimensionless scaling: coordinates with R, stresses with ηU R , time with R/U, velocities with U. ⊲ Periodical boundary conditions.
SLIDE 19
Mesh
M0, 120 elements M3 M4 M5 M6 M7 number of elements 1920 4320 7680 17280 30720
SLIDE 20
Value of C on the center line in the wake for Oldroyd-B
0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 1 1.1 1.2 C x Wi=0.7 Wi=0.86 Wi=0.87 Wi=0.88, time=18 Wi=1.0, matrix logarithm
SLIDE 21
Value of C on the center line in the wake for Giesekus
0.5 1 1.5 2 2.5 3 3.5 1 1.1 1.2 C x Wi=1.17 Wi=1.18 Wi=1.19
SLIDE 22
Dimensionless drag coefficient K for Oldroyd-B (1)
Fan Caola Owens Wi M3 M4 M5 M6 M7 1999 2001 2003 0.0 132.358 132.36 132.384 132.357 0.1 130.363 130.36 0.2 126.626 126.62 0.3 123.193 123.19 0.4 120.596 120.59 0.5 118.836 118.83 118.763 118.827 0.6 117.792 117.78 117.775 0.7 117.448 117.340 117.320 117.315 117.315 117.32 117.291 0.8 117.373 117.36 117.237 0.9 117.787 117.80 117.503 1.0 118.675 118.501 118.471 118.49 117.783 118.030 1.1 119.466 118.031 118.786 1.2 120.860 120.650 120.613 119.764 1.4 123.801 123.587 (123.541) 1.6 127.356 127.172 (127.106) 1.8 (131.458) 131.285 2.0 (135.839)
SLIDE 23 Dimensionless drag coefficient K for Oldroyd-B (2)
116 118 120 122 124 126 128 130 132 134 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 dimensionless drag coefficient K Wi DG, matrix logarithm, mesh M4 Fan et al. Caola et al. Owens et al.
SLIDE 24 Convergence of the numerical solution for Oldroyd-B? (1)
20 40 60 80 100 120 1 2 3 4 5 tauxx s Wi=0.6 M3 M4 M5 M6 M7 M4-1D
20 40 60 80 100 120 1 2 3 4 5 s Wi=0.7 Fan et al. M3 M4 M5 M6 M7 M4-1D M5-1D
SLIDE 25 Convergence of the numerical solution for Oldroyd-B? (2)
Wi=1.0
20 40 60 80 100 120 140 160 180 200 1 2 3 4 5 tauxx s M3 M4 M5 2000 4000 6000 8000 10000 12000 3 4 5 6 7 8 9 10 s M3 M4 M5 M4-1D M5-1D
SLIDE 26
Convergence of the numerical solution for Oldroyd-B? (3)
sxx = (log c)xx (= log cxx!) cxx
SLIDE 27 Convergence of the numerical solution for Oldroyd-B? (4)
100 200 300 400 500 600 1 2 3 4 5 tauxx s Wi=1.4 M3 M4 M5, t=32 M5, t=48
100 200 300 400 500 600 700 800 1 2 3 4 5 tauxx s Wi=1.6 M3 M4 M5, t=72
SLIDE 28 Minimum value of det c
1 2 3 4 5 6 7
5 10 15 log(det(c)) x Wi=1.8, M4
SLIDE 29 Giesekus with α = 0.01; cxx
6000 12000 18000
5 10 15 cxx x Wi=100, M3
SLIDE 30
Giesekus with α = 0.01: convergence?
400 800 1200 1600 2000 1 2 3 4 5 cxx s Wi=5.0 M3 M4 M5 M6 M3-1D M4-1D
SLIDE 31
Conclusions
⊲ We believe that the HWNP has been solved. ⊲ The log conformation representation is the essential part here. ⊲ We think Computational Rheology can now become a normal branch of CFD focusing on study of basic phenomena, such as flow instabilities, morphology and mixing of polymer blends, macroscopic stresses of viscoelastic particle suspensions, drag reduction at high Reynolds number flows, . . . ⊲ Remaining issues include: – identifying model problems/artefacts – resolution problems – efficient schemes/solvers for large 3D problems – . . .
