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
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, . . . )
Introduction (HWNP) (2) For example: y H = 2 R x R L = 30 R Oldroyd-B: break down at Wi = 0 . 88 Giesekus: break down at Wi = 1 . 2
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 . . .
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!
Basic equations Momentum balance and mass balance: ∇ p − ∇ · (2 η s D ) − ∇ · τ = 0 ∇ · u = 0 2 ( L + L T ) and L = ( ∇ u ) T . with D = 1 Constitutive equation: τ = η c = L · c + c · L T + f ( c ) , ˙ λ ( c − I ) with 1 � � c − I Oldroyd-B λ f ( c ) = 1 � c − I + α ( c − I ) 2 � Giesekus λ
Numerical methods c = ∂ c Eulerian setting: ˙ ∂t + u · ∇ c p u c Q 2 P d 1 Q 1 FEM: quads with u ∈ ( Q 2 ) 2 , p ∈ P d 1 DEVSS: projected D ∈ ( P 1 ) 3 DG: c ∈ ( Q 1 ) 3 Explicit time integration: Euler forward ( 1 st order) or Adams-Bashforth ( 2 nd order).
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 exp( bt ) exp( bx c a ) with ⊲ a > 0 , b > 0 1 ⊲ c ( x = 0 , t ) = c ( x, t = 0) = 1 0 x = at L ⇒ exponential growth balanced by convection
Analysis of the HWNP in 1D: a numerical instability (2) c i +1 Numerical discretization in space (first-order upwind): c i c i − 1 d c i d t = − ac i − c i − 1 c i − 2 + bc i ∆ x = ( b − a ∆ x ) c i + a ∆ xc i − 1 i − 2 i − 1 i i + 1 ∆ x or in matrix form d c d t = Ac with γ 0 γ = b − a β = a β γ A = , ∆ x, ... ... ∆ x 0 β γ
Analysis of the HWNP in 1D: a numerical instability (3) To avoid exponential growth in time: γ = b − a ∆ x ≤ a ∆ x ≤ 0 or b Define C : C = ∆ xb C = ∆ x max( b, 0) a , or for a, b ∈ R : | a | then for stability: C ≤ 1 first-order upwind For DG using linear polynomials: C ≤ 2 DG with P 1 ⇒ Failure to balance exponential growth with convection when criterion violated.
Analysis of the HWNP in 1D: a numerical instability (4) ǫ − 1 ǫ = ∂u For Oldroyd-B in extension: a = u , b = 2˙ λ , with ˙ ∂x : ǫ − 1 ∆ x max(2˙ λ, 0) C = | u | ǫ − 1 For Giesekus in extension: a = u , b = 2˙ λ (1 + 2 αc xx ) , ǫ − 1 ∆ x max(2˙ λ (1 + 2 αc xx ) , 0) C = | u |
Solution to the problem: log transformation Solution: solve for s = log c . PDE for s : s = ˙ c c = bc ˙ c = b s = ∂s ∂t + a∂s with ˙ ∂x , and thus: s bx bt a ∂s ∂t + a∂s ∂x = b, x ∈ (0 , L ) , t ≥ 0 1 and c = exp( s ) . x = at L 0 ⇒ No numerical instability. No restrictions on ∆ x
Log transformation for tensors log τ ij ?: No log τ ?: No log c ij ?: No log c : YES! Thus we choose s = log c . BONUS: c = exp( s ) always positive definite! .
Evolution equation for s = log c (1) Spectral decomposition: 3 � n 3 c = c 1 n 1 n 1 + c 2 n 2 n 2 + c 3 n 3 n 3 = c i n i n i , i =1 c 1 0 0 c = 0 c 2 0 n 2 0 0 c 3 n 1 3 3 � � s = log c = log( c i ) n i n i = s i n i n i , � �� � s i i =1 i =1 log c 1 0 0 s 1 0 0 = s = 0 log c 2 0 0 s 2 0 0 0 log c 3 0 0 s 3
Evolution equation for s = log c (2) Substantial derivative of s : 3 3 3 � � � ˙ s i ˙ s i n i ˙ s = s i n i n i + ˙ n i n i + n i i =1 i =1 i =1 s i = d (log c i ) = ˙ c i or with ˙ dt c i 3 3 3 c i ˙ � � � s = ˙ n i n i + s i ˙ n i n i + s i n i ˙ n i c i i =1 i =1 i =1 ⇒ We need expressions for ˙ c i and ˙ n i
Evolution equation for s = log c (3) We get ˙ c i and ˙ n i from the CE: c = L · c + c · L T + f ( c ) ˙ From the diagonal components (in the n i -frame): ˙ n 3 n 3 c i = 2 c i L ii + f i ( c 1 , c 2 , c 3 ) ˙ ˙ n 2 From the off-diagonal components (in the n i -frame): n 2 ˙ 3 n 1 � n i = ω · n i = ˙ ω ji n j n 1 j =1 with the components ω ij of the skewsymmetric tensor ω : ω ij = c i L ji + c j L ij , i � = j, c i � = c j c j − c i
Evolution equation for s = log c (4) End result: 3 3 3 3 f i s i − s j � � � � ˙ s = 2 L ii n i n i + + ( c j L ij + c i L ji ) n i n j n i n i c i c i − c j i =1 i =1 i =1 j =1 � �� � i � = j diagonal � �� � off-diagonal Notes: ⊲ c i = exp( s i ) , c = exp( s ) s i − s j lim ( c j L ij + c i L ji ) = L ij + L ji = 2 D ij ⊲ c i − c j c i → c j
Flow past a cylinder confined between two flat plates ⊲ Oldroyd-B model with η s /η = 0 . 59 , where η = η s + η p y H = 2 R ⊲ Giesekus model with α = 0 . 01 . x R ⊲ Weissenberg number: Wi = λU R , with U average velocity in channel. L = 30 R ⊲ Dimensionless drag force K = F x ηU ⊲ Dimensionless scaling: coordinates with R , stresses with ηU R , time with R/U , velocities with U . ⊲ Periodical boundary conditions.
Mesh M0, 120 elements M3 M4 M5 M6 M7 number of elements 1920 4320 7680 17280 30720
Value of C on the center line in the wake for Oldroyd-B 5 Wi=0.7 4.5 Wi=0.86 Wi=0.87 4 Wi=0.88, time=18 Wi=1.0, matrix logarithm 3.5 3 2.5 C 2 1.5 1 0.5 0 1 1.1 1.2 x
Value of C on the center line in the wake for Giesekus 3.5 Wi=1.17 Wi=1.18 Wi=1.19 3 2.5 2 C 1.5 1 0.5 0 1 1.1 1.2 x
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)
Dimensionless drag coefficient K for Oldroyd-B (2) 134 DG, matrix logarithm, mesh M4 Fan et al. 132 Caola et al. Owens et al. 130 dimensionless drag coefficient K 128 126 124 122 120 118 116 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 Wi
Convergence of the numerical solution for Oldroyd-B? (1) 120 120 M3 Fan et al. M4 M3 M5 M4 100 100 M6 M5 M7 M6 M4-1D M7 80 80 M4-1D M5-1D 60 60 tauxx Wi=0.7 Wi=0.6 40 40 20 20 0 0 -20 -20 0 1 2 3 4 5 0 1 2 3 4 5 s s
Convergence of the numerical solution for Oldroyd-B? (2) Wi=1.0 200 12000 M3 M3 M4 M4 180 M5 M5 10000 160 M4-1D M5-1D 140 8000 120 100 tauxx 6000 80 60 4000 40 20 2000 0 -20 0 0 1 2 3 4 5 3 4 5 6 7 8 9 10 s s
Convergence of the numerical solution for Oldroyd-B? (3) s xx = (log c ) xx ( � = log c xx !) c xx
Convergence of the numerical solution for Oldroyd-B? (4) 600 800 M3 M3 M4 M4 700 M5, t=32 M5, t=72 500 M5, t=48 600 400 500 300 Wi=1.4 400 tauxx tauxx 300 Wi=1.6 200 200 100 100 0 0 -100 -100 0 1 2 3 4 5 0 1 2 3 4 5 s s
Minimum value of det c 7 6 Wi=1.8, M4 5 log(det(c)) 4 3 2 1 0 -15 -10 -5 0 5 10 15 x
Giesekus with α = 0 . 01 ; c xx 18000 Wi=100, M3 12000 cxx 6000 0 -15 -10 -5 0 5 10 15 x
Giesekus with α = 0 . 01 : convergence? 2000 M3 M4 M5 M6 1600 M3-1D M4-1D 1200 cxx 800 Wi=5.0 400 0 0 1 2 3 4 5 s
Recommend
More recommend
