2D velocity-vorticity viscous incompressible flows

2D Velocity-vorticity viscous incompressible flows V International Conference on Computational Methods for Coupled Problems in Science and Engineering COUPLED PROBLEMS 2013 S. Idelsohn, M. Papadrakakis and B. Schrefler (Eds) 2D VELOCITY-VORTICITY VISCOUS INCOMPRESSIBLE FLOWS Ra´ul T´ellez, Habersheel Acevedo, and Alfredo Nicol´as∗ Departamento de Matem´aticas, UAM–Iztapalapa,09340 Iztapalapa, M´exico City, M´exico. e-mail: [email protected] Key words: Velocity–vorticity formulation, direct method, moderate Reynolds numbers, steady state flows Abstract. The 2D unsteady Navier-Stokes equations in its velocity-vorticity formulation, after time discretization, leads to a nonlinear elliptic system which may be solved by iterative methods; however, decoupling the nonlinear vorticity equation from the vorticity in the velocity equations by linear interpolation, a direct method is applied, allowing the vorticity equation to be linear. Steady state convergent flows from the un–regularized unit driven cavity problem, which causes recirculation because of the nonzero boundary condition on the top wall, are reported for Reynolds numbers Re, 400 ≤Re ≤4000. 1 1033 Ra´ul T´ellez, Habersheel Acevedo, and Alfredo Nicol´as 1 INTRODUCTION The main goal of this paper is to present numerical results at moderate Reynolds numbers in the range of 400 ≤Re ≤4000. The results are obtained using a simple numerical scheme for the unsteady Navier-Stokes equations in velocity and vorticity variables. The numerical scheme consists in a direct method to solve the nonlinear steady subproblem that results after a convenient time discretization is applied. This is a novelty contribution since usually that kind of problems are solved by iterative techniques: for instance, [1] for isothermal problems and [2] for thermal ones. The flows are obtained from the well known un-regularized driven (or lid-driven) cavity problem which originates recirculation phenomena due to the nonzero velocity boundary condition on the top wall: the recirculation is originated by the fluid flow coming from the upstream top corner, and then hitting the downstream top corner. At moderate Reynolds numbers, say for instance Re ≤7500, the flow approaches to an asymptotic steady state as ttends to ∞. The meshes in this work follow the size dictated by the thickness of the boundary layer (of order of Re−1 2) and no refining on the mesh is required near the boundary. The results clearly show that as the Reynolds number increases the mesh has to be refined and this in turn leads to decrease the time step: numerically, by stability matters and physically, to capture the fast dynamics of the flow. At this stage the results that are shown are obtained with the contour values given by [3] which are easier to obtain than those given by [4]. 2 THE CONTINUOUS PROBLEM AND THE NUMERICAL METHOD Let Ω ⊂RN(N=2,3) be the region of the flow of a viscous incompressible fluid, and Γ its boundary. It is well known that this kind of unsteady flow is governed by the non-dimensional Navier-Stokes equations; in Ω, t>0; given by ut−1 Re∆u+∇p+(u·∇)u=f(1) ∇·u= 0 (2) where u, and pare the velocity and pressure of the flow, respectively. The parameter Re =UL/ν,ν=kinematic viscosity, is the Reynolds number. The momentum equation (1) must be supplemented with appropriate initial and boundary conditions, for instance u(x,0) = u0(x)inΩ(∇·u0= 0) and u=f1on Γ, t≥0( f1·ndΓ = 0) respectively. Taking twice the curl in the momentum equation in the primitive variables formulation of the Navier-Stokes equations given by (1), and restricting to the two-dimensional case, we obtain the velocity-vorticity formulation: ωt−1 Re∇2ω+u·∇ω=f, (3) 2 1034 Ra´ul T´ellez, Habersheel Acevedo, and Alfredo Nicol´as for the vorticity ωand two Poisson equations for the velocity components ∇2u=−∂ω ∂y and ∇2v=∂ω ∂x,(4) which are related, with (1), by ω=∂v ∂x −∂u ∂y; (5) where u=(u, v) is the velocity with components uand v, see [1] for details for the general 3D case. From equation (2) a function, called the streamfunction, is obtained, which using the relations u=∂ψ ∂y ,v=−∂ψ ∂x ,(6) and, from (5), an elliptic equation for ψis obtained: ∇2ψ=−ω. (7) A description of the numerical method follows. The time derivative that appears in the vorticity equation (3) is approximated with the following scheme, which is unconditionally stable when it is implicitly combined with the Laplacian operator, and it is well behaved when t→+∞, [5], ft(x,(n+ 1)∆t)=3fn+1 −4fn+fn−1 2∆t+O(∆t2),n≥1,(8) where fn=f(x,n∆t); then, replacing that derivative in (3) by (8) we have the following totally implicit approximate problem −∇2un+1 =∂ωn+1 ∂y (9) −∇2vn+1 =−∂ωn+1 ∂x ,un+1|Γ=un+1 bc (10) αωn+1 −ν∇2ωn+1 +un+1 ·∇ωn+1 =fω,ω n+1|Γ=ωn+1 bc (11) where α=3 2∆t,fω=4ωn−ωn−1 2∆tand 1 Re has been replaced by the kinematic viscosity coefficient ν, considering U=L= 1. On the other hand, ubc and ωbc denote boundary condition for uand ω. We recall that the novelty here is to apply a direct method decoupling equations (9-10) from (11) with a linear interpolation from the values of ωknown from the two previous time levels; then, in these conditions (9-10) are solved first, then with the values of uand v, (11) is solved as a linear equation. 3 1035 Ra´ul T´ellez, Habersheel Acevedo, and Alfredo Nicol´as Renaming {un+1,vn+1,ωn+1}as {u, v, ω}, and holding the linear interpolation of the right hand sides of uand vmentioned above, one has the following linear, uncoupled elliptic equations −∇2u=2 ∂ω ∂y n −∂ω ∂y n−1 (12) −∇2v=−2∂ω ∂x n +∂ω ∂x n−1 (13) αω −ν∇2ω+u·∇ω=fω,(14) with the boundary condition mentioned in (9-11). To obtain (u1,v1,ω1) in (9-11) or (1214), a first order approximation is applied for the time derivatives in (3), among other options, through a subsequence with smaller time step; a elliptic problem of the form (12-14) is also obtained. For the spatial discretization to solve elliptic problems of the form (12-14) it may be applied either Finite Differences (FD) or Finite Element (FE) if the region Ω is rectangular. For the (FD) case the following classic second order discretizations are used: With h>0, f(x)=f(x−hx)−2f(x)+f(x+hx) h2 x +O(h2 x), f(y)=f(y−hy)−2f(y)+f(y+hy) h2 y +O(h2 y); f(x)=f(x+hx)−f(x−hx) 2hx +O(h2 x),(15) f(y)=f(y+hy)−f(y−hy) 2hy +O(h2 y). Remark 2. We can also approximate f(x) and f(y) with (7) replacing ∆tby hxand hy. Considering the rectangular flow region Ω = (0,a)×(0,b), a>0, b>0, which for the unitary cavity, a=b= 1. With M,N > 0 integers, we denote hx=a M+1 ,hy=b N+1 , xi=i∗hx,i=0,1,··· ,M,M + 1, and xj=j∗hy,j=0,1,··· ,N,N + 1; and by, say ωij =ω(xi,y j). If M=N,hx=hy=h. Then, with M=N, ∂ω ∂y =ωij+1 −ωij−1 2h(16) and 4 1036 Ra´ul T´ellez, Habersheel Acevedo, and Alfredo Nicol´as ∂ω ∂x =ωi+1j−ωi−1j 2h; (17) for αω −ν∇2ωwe have αωn+1 ij −νωn+1 i−1j+ωn+1 i+1j−4ωn+1 ij +ωn+1 ij−1+ωn+1 ij+1 h2.(18) Similarly for −∇2uand −∇2v. Obtaining, from (18), a constant matrix mfor ω. If we keep u·∇ωin the left hand side, discretizating it like u·(∇ω)=un+1 ij ωn+1 i+1j−ωn+1 i−1j 2h−vn+1 ij ωn+1 ij+1−ωn+1 ij−1 2h(19) and adding it to αω −ν∇2ωa matrix depending on time is obtained, i.e., m=m(t). The case M=N= 3 is shown next; it is split in two part for space reasons, the second part shows the last three columns.                 4ν+h2 ∆t−ν+hv1,1 20−ν+hu1,1 200 −ν−hv1,2 24ν+h2 ∆t−ν+hv1,2 20−ν+hu1,2 20 0−ν−hv1,3 24ν+h2 ∆t00−ν+hu1,3 2 −ν−hu2,1 20 04ν+h2 ∆t−ν+hv2,1 20 0−ν−hu2,2 20−ν−hv2,2 24ν+h2 ∆t−ν+hv2,2 2 00−ν−hu2,3 20−ν−hv2,3 24ν+h2 ∆t 000−ν−hu3,1 200 0000−ν−hu3,2 20 00000−ν−hu3,3 2                                000 000 000 −ν+hu2,1 200 0−ν+hu2,2 20 00−ν+hu2,3 2 4ν+h2 ∆t−ν+hv3,1 20 −ν−hv3,2 24ν+h2 ∆t−ν+hv3,2 2 0−ν−hv3,3 24ν+h2 ∆t                (20) 5 1037 Ra´ul T´ellez, Habersheel Acevedo, and Alfredo Nicol´as All calculations are made with m=m(t), and so far it seems better than considering the constant matrix m, mentioned before, as it is usually done. As is already known the matrix is symmetric and pentadiagonal of size N2×N2. For N(or M) large enough it gives rise a sparse matrix. With α= 0 and ν= 1, it is clear that the above matrix is similar to −∇2uand −∇2v, but constants; these latter matrices are diagonally dominant whereas the previous one is strictly diagonally dominant if ∆t, in α, is sufficiently small, [8]. It should be noted that the above constant matrix mis the consequence of taking u·∇ω, in (14), to the right hand side with the linear interpolation of the two previous time levels, i.e., 2(u·∇ω)n−(u·∇ω)n−1, Therefore, in each time step, three algebraic linear equations systems, associated with the semi-discrete system for u,v, and ωin (12 −14), will be solved, that is, algebraic systems of the form Ax=b, with A=m(t)(or m). Up to now these algebraic linear systems have been solved with LinearSolve, storing previously the matrix as a SparseArray, both Mathematica 8 Commands. The results are reported through the contour of ωand ψ; concerning ψ, this is obtained, in the final time in which ωhas been already computed solving the equation for ψgiven by (7). 3 NUMERICAL EXPERIMENTS AND DISCUSSION The spatial discretization must be supplemented with boundary conditions, for each equation, on Γ for all t≥0, in this case for the driven cavity problem, we have (u, v)= (1,0) on the moving top wall y=1and(u, v) = (0,0) elsewhere; this problem, as mentioned above causes recirculation. Then in terms of this boundary condition for u, the one given by u,vand ωreads u=0,v= 0; ω=∂v ∂x on Γx=0 u=0,v= 0; ω=∂v ∂x on Γx=1 u=0,v= 0; ω=−∂u ∂y on Γy=0 u=1,v= 0; ω=−∂u ∂y on Γy=1 (21) Besides, ω(x,0) = ω0(x) denotes the vorticity initial condition which, by (5), has to satisfy ω0=∂v0 ∂x −∂u0 ∂y if u0=(u0,v 0) is the initial velocity. For moderate Reynolds numbers Re, the following values are reported: Re=400, 1000 and 4000; actually, Re = 400 is not shown but the calculation agrees perfectly with the one we are comparing with.. The results, for each case of Re that is considered, correspond to steady state flows, the iso-contours for the vorticity ωare shown first and the streamlines for the streamfunction ψright below. The cases Re=400 and 1000 agree perfectly with the ones reported by [3] who solve the stationary problem using the streamfunction and vorticity formulation; they use the contour values vo ={−5., −3,−1,1,3,5.}for vorticity and 6 1038 Ra´ul T´ellez, Habersheel Acevedo, and Alfredo Nicol´as vz ={−0.11,−0.10,−0.08,−0.06,−0.04,−0.02,−0.01}for the streamfunction, which give rise the streamlines. Other kind of contour values are v1={−3., −2., −1., −.5,0., .5,1. ,2., 3., 4., 5.}for vorticity given by [4] which are harder to obtain than those of [3], as has already already pointed out in earlier works of the third author. The case Re = 4000 does not agree at all with the one [3] but it does with the one in [1], where some explanation is given for this small disagreement. The mesh size is denoted by hand the time step by ∆tand they are specified in each case under study. Figures 1 and 2 pictures de iso-vorticity contours and streamlines for Re = 1000 respectively, obtained with h=1/256 and ∆t=0.005. Figures 3 and 4 those for Re = 4000, obtained with h=1/600 and ∆t=0.0025. 4 CONCLUSIONS We have presented results for moderate Reynolds numbers that are obtained from a direct numerical method applied to the Navier-Stokes equations in its velocity-vorticity formulation, which is not a common formulation to use. Besides the novelty of the direct method the results become from a variable matrix depending on time and just using Mathematica 8. We are working to improve the scheme to be able to handle high Reynolds numbers as well as to extend it to thermal problems like those in [2] and in [7]. REFERENCES [1] Nicol´as, A. and Berm´udez, B., Viscous incompressible flows by the velocity-vorticity Navier-Stokes equations. Computer Modeling in Engineering & Sciences (2007) 20:73–83. [2] Nicol´as, A. and Berm´udez B., 2D mixed convection viscous incomprensible flows with velocity-vorticity variables. Computer Modeling in Engineering & Sciences (2011) 82:163–178. [3] Schreiber R. and Keller H.B., Driven cavity flow by efficient numerical techniques. J. Comp. Phy. (1983) 40:310–310. [4] Ghia U., Ghia K.N. and Shin C.T., Driven cavity flow by efficient numerical techniques. J. Comp. Phy. (1982) 48:387–411. [5] Glowinski, R. Numerical methods for nonLinear variational problems. SpringerVerlag, (1984). [6] Nicol´as, A. and Berm´udez, B., 2D thermal/isothermal incompressible viscous flows. Int. J. for Num. Meth. in Fluids (2005) 48:349–366. [7] Nicol´as, A. and Berm´udez, B., Thermal recirculation of viscous incomprensible flows in enclosures. European Congress on Computacional Methods in Applied Sciences 7 1039 Ra´ul T´ellez, Habersheel Acevedo, and Alfredo Nicol´as and Engineering(ECCOMAS 2012) (2012); J. Eberhardsteiner et al. (eds.), Vienna, Austria. September 10–14. [8] Burden R. L. and Faires J. D. Numerical Analysis. PWS, Boston, USA., (1985) 8 1040 Ra´ul T´ellez, Habersheel Acevedo, and Alfredo Nicol´as Figure 1:ω: Re=1000 (Keller); h=1/256, dt=.005 9 1041