Report

Structure-Preserving B-spline Methods for the Incompressible Navier-Stokes Equations John Andrew Evans Institute for Computational Engineering and Sciences, UT Austin Stabilized and Multiscale Methods in CFD Spring 2013 Motivation So why do we need another flow solver? Incompressibility endows the Navier-Stokes equations with important physical structure: • Mass balance • Momentum balance • Energy balance • Enstrophy balance • Helicity balance However, most methods only satisfy the incompressibility constraint in an approximate sense. Such methods do not preserve structure and lack robustness. Motivation Consider a two-dimensional Taylor-Green vortex. The vortex is a smooth steady state solution of the Euler equations. 1 0.9 0.8 ¶u + Ñ × ( u Ä u ) + Ñp = f ¶t Ñ×u = 0 0.7 y 0.6 0.5 0.4 with Symmetry BCs 0.3 0.2 0.1 0 0 0.2 0.4 x 0.6 0.8 1 Motivation Conservative methods which only weakly satisfy incompressibility do not preserve this steady state and blow up in the absence of artificial numerical dissipation. 2-D Conservative Taylor-Hood Element Q2/Q1 and h = 1/8 Motivation Methods which exactly satisfy incompressibility are stable in the Euler limit and robust with respect to Reynold’s number. Re = ¥ Re = 10 Increasing Reynold’s number 2-D Conservative Structure-preserving B-splines k’ = 1 and h = 1/8 Motivation Due to the preceding discussion, we seek new discretizations that: • Satisfy the divergence-free constraint exactly. • Harbor local stability and approximation properties. • Possess spectral-like stability and approximation properties. • Extend to geometrically complex domains. Structure-Preserving B-splines seem to fit the bill. The Stokes Complex The classical L2 de Rham complex is as follows: ¾¾ ® H ¾¾¾ ® H(curl) ¾¾® H(div) ¾¾® L2 ¾¾ ®0 1 grad curl div From the above complex, we can derive the following smoothed complex with the same cohomology structure: V Q 0 grad : H V : H 1 1 curl { div F := v ÎL2 : curl v ÎH1 Q : L 2 } The Stokes Complex The smoothed complex corresponds to viscous flow, so we henceforth refer to it as the Stokes complex. V Q 0 grad Scalar Potentials curl Vector Potentials grad div Flow Pressures Flow Velocities curl div The Stokes Complex For simply-connected domains with connected boundary, the Stokes complex is exact. • grad operator maps onto space of curl-free functions grad : curl • curl operator maps onto space of div-free functions curl v V : divv • div operator maps onto entire space of flow pressures divV Q The Stokes Complex Gradient Theorem: grad ( ) d s ( (1)) ( (0 )) curl( ) d a d s Curl Theorem: S S div(v ) dV Divergence Theorem: V V grad curl v da div The Discrete Stokes Subcomplex Now, suppose we have a discrete Stokes subcomplex. h h Vh Q h 0 grad Discrete Scalar Potentials Discrete Vector Potentials curl div Discrete Flow Velocities Discrete Flow Pressures Then a Galerkin discretization utilizing the subcomplex: • Does not have spurious pressure modes, and • Returns a divergence-free velocity field. The Discrete Stokes Subcomplex P roposition. If v h V h satisfies div v h , q h L2 0, q h Q h then div v h 0. Proof. Let q h div v h Q h . d iv v h 2 L 2 0 The Discrete Stokes Subcomplex Proposition. For every qh Î Xh3 there exists a v h Î Xh2 such that di v v h = qh and v h 1 £ b h-1 qh 0 where b h only depends on the norm of P2h . The Discrete Stokes Subcomplex Proof. For qh ÎXh3 , let v ÎH1 be such that div v = qh and v 1 £ b -1 qh 0 . Define v h = P2h v. Then div v h = div P2h v = Ph3 div v = qh and v h 1 £ P2h v 1 £ b -1 P2h qh 0 º b h-1 qh 0 . Structure-Preserving B-splines Review of Univariate (1-D) B-splines: Knot vector on (0,1) and k-degree Ξ = {0, 0, 0, 0.2, 0.4, 0.6, 0.8, 0.8, 1, 1, 1}, k = 2 Knots w/ multiplicity B-spline basis on (0,1) by recursion: ì1 if xi £ x < xi +1 , N i,0 (x ) = í î0 otherwise x - xi N i,k (x ) = N (x ) + xi+k - xi i,k-1 x i+k+1 - x N i+1,k-1 (x ) xi+k+1 - xi+1 Start w/ piecewise constants Bootstrap recursively to k Structure-Preserving B-splines Review of Univariate (1-D) B-splines: • Open knot vectors: • Multiplicity of first and last knots is k+1 • Basis is interpolatory at these locations • Non-uniform knot spacing allowed • Continuity at interior knot a function of knot repetition k=2 Structure-Preserving B-splines Review of Univariate (1-D) B-splines: • Derivatives of B-splines are B-splines 1 1 d/dx onto 0 0 1/3 2/3 k=4 1 0 0 1/3 2/3 k=3 1 Structure-Preserving B-splines Review of Univariate (1-D) B-splines: • We form curves in physical space by taking weighted sums. 4 5 3 2 0 “control mesh” 1 - control points Quadratic basis: - knots Structure-Preserving B-splines Review of Multivariate B-splines: • Multivariate B-splines are built through tensor-products • Multivariate B-splines inherit all of the aforementioned properties of univariate B-splines In what follows, we denote the space of n-dimensional tensorproduct B-splines as polynomial degree in direction i k1 ,k2 ,...,ki ,...,kn a1 ,a 2 ,...,a i ,...,a n S i th continuity vector Structure-Preserving B-splines Review of Multivariate B-splines: • We form surfaces and volumes using weighted sums of multivariate B-splines (or rational B-splines) as before. Control mesh Mesh Structure-Preserving B-splines Two-dimensional Structure-Preserving B-splines Define for the unit square: k ,k Sˆ h : S 11 ,22 k ,k 1 k 1,k Rˆ h : S 11 ,22 1 S 11 1,22 k 1,k 1 Wˆ h : S 11 1,22 1 In the context of fluid flow: Sˆ h Sp ace o f Stream Fu n ctio n s Rˆ h Sp ace o f Flo w V elo cities Wˆ h Sp ace o f Pressu res and it is easily shown that: curl ˆ div Sˆ h R Wˆ h 0 h Structure-Preserving B-splines Two-dimensional Structure-Preserving B-splines: Mapped Domains On mapped domains, the Piola transform is utilized to map flow velocities. Pressures are mapped using an integral preserving transform. Structure-Preserving B-splines Two-dimensional Structure-Preserving B-splines We associate the degrees of freedom of structure-preserving B-splines with the control mesh. Notably, we associate: Sˆ h : co n tro l p o in ts Rˆ h : co n tro l faces Wˆ h : co n tro l cells Structure-Preserving B-splines Two-dimensional Structure-Preserving B-splines, k1 = k2 = 2: Structure-Preserving B-splines Two-dimensional Structure-Preserving B-splines, k1 = k2 = 2: Structure-Preserving B-splines Two-dimensional Structure-Preserving B-splines, k1 = k2 = 2: Structure-Preserving B-splines Two-dimensional Structure-Preserving B-splines, k1 = k2 = 2: Structure-Preserving B-splines Two-dimensional Structure-Preserving B-splines, k1 = k2 = 2: Structure-Preserving B-splines Two-dimensional Structure-Preserving B-splines, k1 = k2 = 2: Structure-Preserving B-splines Two-dimensional Structure-Preserving B-splines, k1 = k2 = 2: Structure-Preserving B-splines Two-dimensional Structure-Preserving B-splines, k1 = k2 = 2: Structure-Preserving B-splines Three-dimensional Structure-Preserving B-splines Define for the unit cube: k ,k ,k Sˆ h : S 11 ,22 ,3 3 k 1,k ,k k ,k 1,k k ,k ,k 1 Nˆ h : S 11 1,22 ,3 3 S 11 ,22 1,3 3 S 11 ,22 ,3 3 1 k ,k 1,k 1 k 1,k ,k 1 k 1,k 1,k Rˆ h : S 11 ,22 1,3 3 1 S 11 1,22 ,3 3 1 S 11 1,22 1,3 3 k 1,k 1,k 1 Wˆ h : S 11 1,22 1,3 3 1 It is easily shown that: grad curl div Sˆ h Nˆ h Rˆ h Wˆ h 0 Structure-Preserving B-splines Three-dimensional Structure-Preserving B-splines Flow velocities: map w/ divergence-preserving transformation Flow pressures: map w/ integral-preserving transformation Vector potentials: map w/ curl-conserving transformation h h S h : w : w oF Sˆ h : u N h : v : D F Rh W h : h h T v h oF Nˆ h 1 : det D F D F h u o F Rˆ h h h p : det( D F ) p oF Wˆ h Structure-Preserving B-splines Three-dimensional Structure-Preserving B-splines We associate as before the degrees of freedom with the control mesh. Sˆ h : control points Nˆ h : con trol ed ges Rˆ h : con tro l faces Wˆ h : control cells Structure-Preserving B-splines Three-dimensional Structure-Preserving B-splines Control points: Scalar potential DOF Control edges: Vector potential DOF Structure-Preserving B-splines Weak Enforcement of No-Slip BCs Nitsche’s method is utilized to weakly enforce the no-slip condition in our discretizations. Our motivation is threefold: • Nitsche’s method is consistent and higher-order. • Nitsche’s method preserves symmetry and ellipticity. • Nitsche’s method is a consistent stabilization procedure. Furthermore, with weak no-slip boundary conditions, a conforming discretization of the Euler equations is obtained in the limit of vanishing viscosity. Structure-Preserving B-splines Weak Enforcement of Tangential Continuity Between Patches On multi-patch geometries, tangential continuity is enforced weakly between patches using a combination of the symmetric interior penalty method and upwinding. ! 1 F1 ! F2 "2 ! 3 2 1 x2 F3 0 0 1 "1 x1 Summary of Theoretical Results Steady Navier-Stokes Flow ( ( ) ) Ñ × u Ä u + Ñp = Ñ × 2nÑsu + f in W Ñ ×u = 0 u=0 in W on ¶W • Well-posedness for small data • Optimal velocity error estimates and suboptimal, by one order, pressure error estimates • Conforming discretization of Euler flow obtained in limit of vanishing viscosity (via weak BCs) • Robustness with respect to viscosity for small data Summary of Theoretical Results Unsteady Navier-Stokes Flow ( ) ( ) ¶u/ ¶t + Ñ × u Ä u + Ñp = Ñ × 2nÑsu + f in W ´ (0,T ) Ñ ×u = 0 in W ´ (0,T ) u = 0 on ¶W ´ (0,T ) u = u0 on W • Existence and uniqueness (well-posedness) • Optimal velocity error estimates in terms of the L2 norm for domains satisfying an elliptic regularity condition (local-in-time) • Convergence to suitable weak solutions for periodic domains • Balance laws for momentum, energy, enstrophy, and helicity • Balance law for angular momentum on cylindrical domains Spectrum Analysis Spectrum Analysis Consider the two-dimensional periodic Stokes eigenproblem: 2 u p u in s u 0 in w / P eriodic B C s We compare the discrete spectrum for a specified discretization with the exact spectrum. This analysis sheds light on a given discretization’s resolution properties. Spectrum Analysis Spectrum Analysis: Structure-Preserving B-splines Spectrum Analysis Spectrum Analysis: Taylor-Hood Elements Spectrum Analysis Spectrum Analysis: MAC Scheme Selected Numerical Results Steady Navier-Stokes Flow: Numerical Confirmation of Convergence Rates é ù 2e x (-1+ x)2 x 2 (y 2 - y)(-1+ 2y) u(x,y) = ê x 2 2 ú ë(-e (-1+ x)x(-2 + x(3 + x))(-1+ y) y ) û -424 + 156e + (y 2 - y)(-456 + e x (456 + p(x,y) = x 2 (228 - 5(y 2 - y)) + 2x(-228 + (y 2 - y)) + 2x3 (-36 + (y 2 - y)) + x 4 (12 + (y 2 - y)))) 2-D Manufactured Vortex Solution Selected Numerical Results Steady Navier-Stokes Flow: Numerical Confirmation of Convergence Rates 2-D Manufactured Vortex Solution Selected Numerical Results Steady Navier-Stokes Flow: Numerical Confirmation of Convergence Rates 2-D Manufactured Vortex Solution Selected Numerical Results Steady Navier-Stokes Flow: Numerical Confirmation of Convergence Rates Re 0 1 10 100 1000 10000 Energy 1.40e-2 1.40e-2 1.40e-2 1.40e-2 1.40e-2 1.40e-2 H1 error u 1.40e-2 1.40e-2 1.40e-2 1.40e-2 1.40e-2 1.40e-2 L2 error - u 2.28e-4 2.28e-4 2.28e-4 2.28e-4 2.28e-4 2.28e-4 L2 error - p 3.49e-4 3.49e-4 1.98e-4 1.96e-4 1.96e-4 1.96e-4 Robustness with respect to Reynolds number k’ = 1 and h = 1/16 2-D Manufactured Vortex Solution Selected Numerical Results Steady Navier-Stokes Flow: Numerical Confirmation of Convergence Rates Re 0 1 10 100 1000 10000 H1 error u 6.77e-4 6.77e-4 7.11e-4 2.26e-3 2.16e-2 X L2 error u 6.54e-4 6.54e-6 6.79e-6 1.97e-5 1.86e-4 X L2 error p 1.96e-4 1.96e-4 1.96e-4 1.96e-4 1.96e-4 X Instability of 2-D Taylor-Hood with respect to Reynolds number Q2/Q1 and h = 1/16 2-D Manufactured Vortex Solution Selected Numerical Results Steady Navier-Stokes Flow: Lid-Driven Cavity Flow U u ÎW1,r , p ÎW 0,r , 1< r < 2 H H Selected Numerical Results Steady Navier-Stokes Flow: Lid-Driven Cavity Flow at Re = 1000 Selected Numerical Results Steady Navier-Stokes Flow: Lid-Driven Cavity Flow at Re = 1000 Method umin vmin vmax k’ = 1, h = 1/32 -0.40140 -0.39132 0.54261 k’ = 1, h = 1/64 -0.39399 -0.38229 0.53353 k’ = 1, h = 1/128 -0.39021 -0.37856 0.52884 k’ = 2, h = 1/64 -0.38874 -0.37715 0.52726 k’ = 3, h = 1/64 -0.38857 -0.37698 0.52696 Converged -0.38857 -0.37694 0.52707 Ghia, h = 1/156 -0.38289 -0.37095 0.51550 Selected Numerical Results Unsteady Navier-Stokes Flow: Flow Over a Cylinder at Re = 100 L H D U Lout Selected Numerical Results Unsteady Navier-Stokes Flow: Flow Over a Cylinder at Re = 100 Patch 1 Patch 4 Patch 2 Patch 3 Patch 5 Selected Numerical Results Unsteady Navier-Stokes Flow: Flow Over a Cylinder at Re = 100 Simulation Details: - Full Space-Time Discretization - Method of Subgrid Vortices - Linears in Space and Time - D = 2, H = 32D, L = 64D - Time Step Size: 0.25 - Trilinos Implementation - GMRES w/ ILU Preconditioning - 3000 time-steps (shedding initiated after 1000 time-steps) - Approximately 15,000 DOF/time step Selected Numerical Results Unsteady Navier-Stokes Flow: Flow Over a Cylinder at Re = 100 Selected Numerical Results Unsteady Navier-Stokes Flow: Flow Over a Cylinder at Re = 100 Quantity of Interest: Strouhal Number: St = fD/U Computed Strouhal Number: 0.163 Accepted Strouhal Number: 0.164 Selected Numerical Results Unsteady Navier-Stokes Flow: Three-Dimensional Taylor-Green Vortex Flow 3-D Periodic Flow Simplest Model of Vortex Stretching No External Forcing Selected Numerical Results Unsteady Navier-Stokes Flow: Three-Dimensional Taylor-Green Vortex Flow Time Evolution of Dissipation Rate Reproduced with permission from [Brachet et al. 1983] Selected Numerical Results Unsteady Navier-Stokes Flow: Three-Dimensional Taylor-Green Vortex Flow at Re = 200 Enstrophy isosurface at time corresponding to maximum dissipation Selected Numerical Results Unsteady Navier-Stokes Flow: Three-Dimensional Taylor-Green Vortex Flow at Re = 200 Convergence of dissipation rate time history with mesh refinement (k’ = 1) Selected Numerical Results Unsteady Navier-Stokes Flow: Three-Dimensional Taylor-Green Vortex Flow at Re = 200 Convergence of dissipation rate time history with degree elevation (h = 1/32)