Jtam.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 47, 4, pp. 779-796, Warsaw 2009 VORTEX IN CELL METHOD FOR EXTERIOR PROBLEMS Henryk Kudela Tomasz Kozłowski Wroclaw University of Technology, Department of Numerical Flow Modelling, Wrocław, Poland e-mail: tomasz.kozlowski@pwr.wroc.pl; henryk.kudela@pwr.wroc.pl The ”vortex in cell” method was used to model the flow past a solid body. Problems connected with flows past a body belong to the cate- gory of exterior problems. The main difficulty here is to establish the boundary condition far from the body, to solve equations of motion on the numerical mesh. The mesh requires a limited calculation area and the boundary condition far from the body. A method of obtaining the accurate boundary condition in such flows was presented. Testing cal- culations were performed for the flow past a cylinder for a wide range of Reynolds numbers. The flow past an ellipse was also performed. A good agreement with experimental results and calculations of other re- searchers was obtained. Key words: vortexmethod, flow past cylinder, external problem 1. Introduction Vortex particle methods have obtained already the state of maturity and pre- cision that one can compare them to the accuracy and efficiency of the classic method like the finite differencemethod or spectralmethod (Cottet etal, 2000, 2002). The importance of the vortex particle method lies in the possibility of analysing more easily and directly the vorticity field due to the fact that in computations particles that carry information about the vorticity field are used. The vorticity is the key to understanding many hydrodynamic pheno- mena (Kudela and Machela, 2007). An attractive feature of the method is also elimination of difficulties connected with modelling of the inertial terms ofNavier-Stokes equations. In the presentwork, theVortex-In-Cell (VIC)me- thod was used to determine the flow past a solid body in an infinite domain. Such a problembelongs to the category of exterior problems. The difficulty in 780 H. Kudela, T. Kozłowski numerical modelling of such a problem is related to the establishment of the boundary condition far from the body. The mesh requires a finite calculation area, and the boundary condition far from the body is needed. In direct vor- texmethods, which use the Biot-Savart law to calculate the velocity field, the condition at infinity is satisfied automatically. However, direct vortexmethods are very time consuming. The VIC method is hundreds times faster than the directmethod (Cottet andPoncet, 2003). TheVICmethod retains the impor- tant features of particlemethods that solve the convective part of fluidmotion. The grid is used to compute the velocity by solution of the Poisson equation for stream function and vorticity. Tomanage the far field boundary condition, we used the method described in Anderson and Reider (1996), Wang (1999). Thatmethod takes advantage of the fact that the domain of non-zero vorticity around the solid body is limited to the domain around the obstacle. In the far distance from the body, where the vorticity is zero, the asymptotic properties of the Laplace equation and its solution represented by Fourier series is used. It enables determination of the proper far field boundary condition. In order to use the fast direct algorithm for the solution of the Poisson equation on a uniform grid we applied the conformal mapping. We tested the VIC method for the flow over a circular cylinder and elongated shape of an ellipse.Wewere especially interested in the ability of theVICmethod for reproduction of very fine vortex structures near the solid body. 2. Vortex particle method and conformal mapping Two-dimensional equations of an incompressible viscous fluid transformed to the Helmhotz equation for stream function and vorticity are ωt+uωx+vωy = ν∆ω ∆Ψ =−ω (2.1) The velocity of the fluid is determined as u= ∂yΨ v=−∂xΨ (2.2) The subscripts denote partial differentiation with respect to the variable that appears in the subscript. Equations (2.1) need to be completed with initial conditions for the vorticity field ω(x,y,0)=ω0 = rot(u0), where u0 denotes the initial velocity field u(x,y,0) = u0 and non-slip boundary condition on the solid wall for the velocity. Vortex in cell method for exterior problems 781 Zeroingof thenormalvelocity component is obtainedby setting Ψ = const on the wall. Zeroing of the tangential velocity is satisfied by introducing a proper portion of vorticity that ensures u ·s=0 were s is a unit tangential vector. To transform the non-rectangular physical region (x,y-variables) to the rectangular one (ξ,η), we used the conformal mapping. In the new variables (ξ,η), equations (2.1) take form (Graham, 1988) ωt+uωξ+vωη = ν J ∆ω ∆Ψ =−Jω (2.3) where J denotes Jacobian of the conformal transformation J =det ∣ ∣ ∣ ∣ ∣ xξ xη yξ yη ∣ ∣ ∣ ∣ ∣ (2.4) and the velocities are expressed by the formulas u= 1 J ∂ηΨ v=− 1 J ∂ξΨ (2.5) In the VIC method, the continuous vorticity field is approximated with a discrete particles distribution. The flow region is covered with the numerical grid h = ∆η = ∆ξ, (ih,jh), (i = 1,2, . . . ,nξ, j = 1,2, . . . ,nη). In the grid nodes, the particles with Γj circulation are placed Γj = ∫ A ω dξdη (2.6) where A=h2, and ω(ξ,η)= ∑ p Γpδ(ξ− ξp)δ(η−ηp) (2.7) The viscous splitting algorithm was used for solution of (2.3) and (2.5). Firstly, the inviscid fluid equation was solved ωt+uωξ +vωη =0 (2.8) From (2.3)1 it stems that the vorticity is constant along the trajectory of fluid particles. According to the Helmholtz theorem (Wu et al., 2005), vortex particles are moving like material fluid particles. Differential equation (2.8) is replaced by the infinite set of ordinary equations dξ dt =u dη dt = v ξ(0,α1)=α1 η(0,α2)=α2 (2.9) 782 H. Kudela, T. Kozłowski where α = (α1,α2) means Lagrangian coordinates of fluid particles. The infinite set ofdifferential equations (2.9) is replacedbyafiniteone.Thenumber of particles is equal to the number of grid nodes. The Lagrangian parameter α takes in each time step the value αi,αj =(ξi,ηj).Thefinite set of equations (2.9) was solved by theRunge-Kuttamethod of the fourth order. The velocity field was obtained by solving Poisson equation (2.3)2 on the numerical grid andmaking use of (2.5). The velocities of particles that are found between the grid nodes were calculated by the interpolation formula u(ξp,ηp)= ∑ j lj(ξp,ηp)uj (2.10) where lj denotes the two-dimensional bilinear interpolation Lagrange base. The viscosity was taken into account by solving the diffusion equation ωt = ν J ∆ω ω(ξ,η,0)=ω0 ω|wall =ωs (2.11) where ωs was calculated on the basis of Poisson equation (2.3)2. Taking in- to account the non-slip condition u = 0, it gives the vorticity on the wall ω=−Ψηη/J. The value of this vorticity was calculated by the Briley formula (Weinan and Liu, 1996) ω(0,j)s =− 1 J 108Ψ1,j −27Ψ2,j +4Ψ3,j 18h2 +O(h4) (2.12) where h denotes the grid step, index 0 refers to the wall and index i (i=1,2,3), to the distance ih from the wall. In order to solve the diffusion equation on the grid, after displacement of particles according to equation (2.9), it was necessary to pass the vorticity from the particles to the grid nodes. In each time step, redistribution of mas- ses of vortex particles onto neighbouring grid nodes was performed (Fig.1), according to the formula ωij = 1 h2 ∑ p Γpϕh(ξ)ϕh(η) (2.13) where ϕh(ξ)=ϕ (ξ− ξi h ) ϕh(η)=ϕ (η−ηj h ) Indexes p, i, j refer to vortex particles andgrid nodes, respectively, and ϕ(·) is the kernel of the interpolation function. The fourth order interpolation kernel, Vortex in cell method for exterior problems 783 known in the literature as the M4, was used (Becerra-Sagredo, 2003; Cottet and Koumoutsakos, 2000) ϕ(x)=            1− 5 2 x2+ 3 2 |x|3 for |x|< 1 1 2 (2−|x|)2(1−|x|) for 1¬ |x| ¬ 2 0 for |x|> 2 (2.14) For particles near the wall, one-sided interpolation functions were used (2.15) derived according to the algorithm presented in Becerra-Sagredo (2003) ϕ(x)=              1+ 1 2 x2− 3 2 |x| for j=0, |x| ¬ 1 −x2+2|x| for j=1, |x| ¬ 1 1 2 x2− 1 2 |x| for j=2, |x| ¬ 1 (2.15) Fig. 1. Redistribution of particle masses onto neighboring grid nodes, (a) for particles inside the computational domain (at least one cell from the wall), (b) for the particles in the vicinity of the wall The interpolation of particle masses onto grid nodes has fundamental me- aning for precision of the VIC method. Both interpolation kernels conserve three first moments (Cottet and Koumoutsakos, 2000) ∑ p xαpϕ (xp−x h ) =xα α=0,1,2 (2.16) Diffusion equation (2.11)1 was solved on the numerical grid with the alterna- ting direction implicit scheme (Thomas, 1995) ωn+ 1 2 =ωn+ ∆t 2J ν(Λξξω n+Ληηω n+1 2) (2.17) ωn+1 =ωn+ 1 2 + ∆t 2J ν(Λξξω n+1+Ληηω n+1 2) 784 H. Kudela, T. Kozłowski where Λ means the three point finite difference quotient that approximates the Laplace operator with respect to the variable that was put in the lower index. The solution of the diffusion equation ends the calculations in the n-th time step of the vortex particles method. The full computational algorithm can be summarised as follows: • The whole computational domain is covered by the numerical grid. In every grid node, the vortex particle is placed with the circulation given by Eq. (2.6). • The Poisson equation for stream function and vorticity is solved, Eq. (2.3)2, the velocities on the grid nodes are computed, Eq. (2.5). • The equation of motion of the inviscid fluid is solved by displacements of the particles, Eq. (2.9). • After displacements of theparticles, the redistributionprocess of particle circulation to the grid nodes takes place, Eqs. (2.13)-(2.15). • The proper amount of vorticity that is needed to fulfil the non-slip con- dition on the wall is calculated, Eq. (2.12). • The diffusion equation is solved, Eq. (2.11)1, with scheme (2.17). This ends one time step. 3. Far field boundary condition for the Poisson equation Here we present amethod how to obtain the accurate boundary condition for the Poisson equation for the exterior problem.We follow the description that was given by Wang (1999). Due to the fact that we used the finite difference method for the solutionof thePoissonequation, our realisation is abitdifferent than inWang (1999). Let us write the boundary value problem for the Poisson equation in an unbounded area ∆Ψ =−ω Ψ(R1,θ)= z(θ) (3.1) where z(θ) is a known function. In the problem of the flow past a solid body, the vorticity ω is different from zero, in some limited area around the body. This area can be enclosed in a finite ring with radius r = a, Fig.2. The condition for Ψ(R2,θ) is unknown. In the area r > a, vorticity equals zero. Since equation (3.1)1 had been solved by the finite difference method, the calculation area was enclosed in the ring R1 < r a. The error Ψ1 has in the ring R1 < r < R2 a Fourier series representation (Farlow, 1993; Wang, 1999) Ψ1(r,θ)= an ln(r)+ ∑ n­1 ( bn+ cn rn ) einθ (3.3) where an, bn, cn are Fourier coefficients which will be determined from boun- dary conditions. Outside the domainwhere ω(r,θ)≡ 0, it is for r >a,Ψ(r,θ) can be presented as (Farlow, 1993; Wang, 1999) Ψ(r,θ)=Γ ln(r)+ ∑ n­1 dn rn einθ (3.4) The Fourier coefficients bn, cn, dn in equations (3.3) and (3.4) can be deter- mined on the basis of the following boundary conditions 786 H. Kudela, T. Kozłowski Ψ(R1,θ)= 0 Ψtrial(R2,θ)=Ψ(R2,θ)+Ψ1(R2,θ) (3.5) ∂Ψtrial ∂r (R2,θ)= ∂Ψ ∂r (R2,θ)+ ∂Ψ1 ∂r (R2,θ) Conditions (3.5)2,3 express the continuity of Ψ and its first derivative with respect to r for the radius r = R2. The derivative of the trial solution ∂Ψr(R2,θ) = g(θ) was calculated by the one-sided finite difference formula of the fourth order g(θ)= 1 12h [3Ψtrial(rn−4,θ)−16Ψtrial(rn−3,θ)+36Ψtrial(rn−2,θ)+ (3.6) −48Ψtrial(rn−1,θ)+25Ψtrial(rn,θ)]+O(h 4) Having the values of g(θ), we expanded this expression into the Fourier series g(θ)= f0(R2)+ ∑ n>1 fn(R2)e inθ (3.7) Series (3.7) approximates the left side of formula (3.5)3. The right side of (3.5)3 was obtained by differentiating series (3.3) and (3.4). Using equations (3.5) and (3.7) the Fourier coefficients an, bn, cn, dn can be expressed as an =0 bn = fn(R2) 2nRn−12 cn =−bnR 2n 1 dn = bn(R 2n 1 −R 2n 2 ) (3.8) One can see from (3.8), that coefficients bn, cn, dn are expresed by the known coefficient fn. Using the trial solution with the fictitious boundary condition Ψtrial(R2,θ), the accurate value of the boundary condition can be determined Ψ(R2,θ)=Ψtrial(R2,θ)− ∑ n­1 fn(R2) 2n R2 ( 1− R2n1 R2n2 ) einθ (3.9) The procedure for getting of the proper boundary value for r = R2 goes as follows: • First, thePoisson equationwas solved using the finite differencemethod with the fictitious boundary condition on the larger radius Ψtrial(R2,θ) • The derivative of the solution obtained on the external radius g(θ) was determined numerically and expanded into the Fourier series Vortex in cell method for exterior problems 787 • The correction of the boundary condition was made according to equ- ation (3.9) • The Poisson equation was solved again with the accurate condition on the radius R2. ThePoisson equationwas solved by a fast direct algorithm of O(h4) order of accuracy. The order of the correction method was the same as the order of the used solver for the Poisson equation. 4. Flow past a circular cylinder At first, we checked the above procedure for the far field boundary condi- tion and the vortex particle ability to reproduce very fine vortex structures. The calculations were carried out for a well documented flow over a circular cylinder. The following conformal transformation was used x+iy=er+iθ (4.1) The Jacobian of the transformation is presented with the following formula J =e2r (4.2) Thephysical flow region inFig.3awas transformedonto the rectangle, Fig.3b. Fig. 3. Transformation of the physical domain by conformal mapping (4.1) for the flow past a cylinder: (a) grid in the physical domain, (b) rectangular grid in the computational domain The calculations were carried out on a regular numerical mesh h=∆θ=∆r, (ih,jh), (i=1,2, . . . ,nr, j=1,2, . . . ,nθ). The Poisson equation for stream function and vorticity was solved in the transformed coordinates (r,θ) which corresponded, according to (4.1), 788 H. Kudela, T. Kozłowski to R1 = 1, R2 = e π in the omputational domain. The Poisson equation was solved with the following boundary conditions ∆Ψ =− ω J Ψ(R1,θ)= 0 Ψ(R2,θ)=Ψc (4.3) where Ψc denotes the corrected stream function value as it was described in the previous section. At the first step it was assumed that Ψ∞ =U∞R2cosθ, U∞ = 1. When treating the problem with a sudden impulsive flow over the cylinder, one encounters the transition problem which has consequences in oscillatory solutions in the initial time steps. To eliminate the transition, a smooth startwas performed. In the first calculation step, the vorticity equaled zero and it a potential solution was assumed. The velocity at infinity of the flow, changed linearly from zero to U∞, in the time interval [0, ts]. Then we restarted the calculation using the obtained solution for t = ts as the new initial condition for the vorticity field. In each time step that followed, on the external radius a correction of the boundary condition Ψc was made as described in Section 3. Calculations were carried out for a wide range of Reynolds numbers Re =U∞2R1/ν. The parameter of the numerical grid, time steps and ts are presented in Table 1. For comparison, the drag force and drag coefficient were determined on the surface of the cylinder (Anderson and Reider, 1996) FD = bµ 2π ∫ 0 ( ω−R ∂ω ∂r ) Rsinθ dθ (4.4) CD = FD ρu2∞Rb where µ is the dynamic coefficient of fluid viscosity, µ = ρν, ρ is density of the fluid, b denotes the unit length of the cylinder and R denotes its radius. We assumed R=1, ρ=1. Table 1.Grid parameters and time steps for the flow past a circular cylinder Re nr nθ ∆t ts 550 256 512 0.01s 0.2s 3000 512 1024 0.005s 0.3s 9500 2048 4096 0.002s 0.4s A characteristic feature of the flow past the cylinder is appearance of re- circulation bubbles behind the cylinder for Re< 50 (Bourd and Coutanceau, Vortex in cell method for exterior problems 789 1980; Van Dake, 1982). For Reynolds numbers Re > 50, the bubbles are no longer stationary andalternately separate fromthe cylinderwall creatingKar- man’s vortex street, Fig.5 (Van Dake, 1982). In Fig.4 it a comparison of the numerical streamline distribution for the Reynolds number Re = 3000 and t = 5 with the experimental one is shown (Van Dake, 1982). The ability of reproduction of a fine vortex structure on the cylinder is clearly visible. The numerical results for Re = 9500 is presented in Fig.6. Qualitatively, a very good agreement was obtained with the visualization data that were published in Van Dake (1982). Fig. 4. Comparison of the streamline distribution for the flow past the cylinder, Re=3000, t=5s, experiment (Van Dake, 1982) (left) and numerical computations (right) Fig. 5. Exemplary picture of the vortex Karman street, Re=550, t=100s In Fig.7 a comparison of the calculated drag coefficient with the results of Koumoutsakos and Leonard (1995) is presented. Some discrepancies are obse- rved for times t < 1. They can be related to the problem how one managed the transition at the beginning of calculations (Anderson and Reider, 1996; Koumoutsakos and Leonard, 1995). In the present work, the initial condition for vorticity was obtained by smooth start. The noted differences at the be- ginning did not influence the results significantly in later time steps. For the Reynolds number Re= 9500 the correct reproduction of vortex structures in the vicinity of the wall requires a very dense numerical mesh (see Table 1). 790 H. Kudela, T. Kozłowski Fig. 6. Comparison of evolution of streamlines with the experimental data. Evolution of vorticity for Re=9500; (a) experimental data (Van Dake, 1982), (b) streamlines from computation, (c) computed isolines of vorticity 5. Flow past an elliptic cylinder The next test was dedicated tomore complicated shape of the solid body.We investigated theflowpasta thin ellipseperpendicular to theflowdirection.The calculations were carried out in elliptic coordinates, which were transformed onto a regular Cartesian mesh using conformal transformations (Fig.8a) x+iy=cosh(ξ+iη) (5.1) Vortex in cell method for exterior problems 791 Fig. 7. Comparison of time evolution of the drag coefficient with the results by Koumoutsakos and Leonard (1995) Fig. 8. (a) Elliptical grid in the physical domain, (b) rectangular grid in the computational domain after transformation The Jacobian of this mapping is J =cosh2ξ− cos2η (5.2) The calculations were carried out for the following parameters: ellipse thickness e = 0.25, ellipse length L = 2, Reynolds number Re = 500 and 10000. On the sharp ends of the ellipse, significant vorticity is generated. A dense grid in this region (Fig.8a), allows one to reproduce flow phenomena generated in this area, however a small time step is required, ∆t = 0.0005. For low Reynolds numbers, behind the ellipse, vortex structures are genera- ted (Fig.9), which after some time become separated from the wall, forming Karman’s vortex street. At high Reynolds numbers, along the line that forms the main vortex structure, smaller secondaryvortex structures appear (Figs.10, 11).Numerical calculations carried out in Koumoutsakos and Shields (1996), Wang et al. (1999) qualitatively reflect the phenomenon discussed. It was found byWang 792 H. Kudela, T. Kozłowski Fig. 9. Vorticity and streamlines for the flow past an ellipse, Re=500 Fig. 10. Visualization of growth of vortices on the accelerated plate VanDake (1982) Fig. 11. Vorticity and streamlines for the flow past the ellipse, Re=10000. Secondary vortices along the primary vortex are visible et al. (1999) that the frequency of the small secondary vortex is controlled by the vortex which appears near the sharp edge of the ellipse (corner vortex), Fig.12b.After sudden start, the corner vortex reaches aquasi-stationary state. The calculations were carried out for the same parameters as in Wang et al. (1999). The vorticity generated at the end of the ellipse very well agrees with the results ofWang et al. (1999). In Fig.10, the experimental visualization of Van Dake (1982) was presented. The similarity of the visualization data and numerical results is visible. Vortex in cell method for exterior problems 793 Fig. 12. Image of the vorticity field and streamlines in the vicinity of the wall (a) Re=500, (b) Re=10000, (c) comparison of the vorticity distribution on the sharp end of the ellipse with the results byWang et al. (1999) 6. Remarks on realisation of boundary conditions The generated portion of the vorticity on the solid boundary to satisfy the no-slip condition on the wall must fulfill some constraints. Let us denote the whole vorticity in the domain as Ω(t)= ∫ D ω(x,y,t) dxdy (6.1) Taking into account that the divergence of ∇· (ωu) = ∇ω ·u, (∇·u = 0, ∇= (∂x,∂y), where n is the unit normal vector, then integrating Helmholtz equation (2.1)1 by sides, and using theGauss divergence theorem, one obtains d dt Ω(t)+ ∫ ∂D (ωu ·n−ν∇ω ·n) dC =0 (6.2) where C = ∂D is the boundary of the domain.On the solid boundarywe have u ·n=0. So it stems from equation (6.2) that d dt Ω(t)= ν ∫ ∂D ∂ω ∂n dC (6.3) It is clear that global vorticity (6.1) in the domain can only change due to the flux of vorticity through the boundary. But the right-hand side of equation (6.3) should be zero to satisfy the single-value property of pressure. Going around the boundary of the profile, we have ν ∫ ∂D ∂ω ∂n dC = ∫ ∂D ∂p ∂s dC (6.4) 794 H. Kudela, T. Kozłowski where smeans the tangent unit vector. So global vorticity (6.1) in the domain should not change, and in our case should be zero. The changes of the global vorticity in the domain during calculations for the flow over the cylinder was presented in Fig.13. One can see that the maximum fluctuations of the total value of vorticity depend on Reynolds number. The maximum fluctuations is on the level ±2 · 10−10, for the Reynolds number Re = 9500. Taking into account the comparison with the experimental and numerical data presented in the previous Section, it seems to be acceptable. For an acceptable outcome of the non-slip boundary conditions at greater Reynolds numbers, the grid near the boundary should be finer. In literature (Koumoutsakos et al., 1994; Cottet andKoumoutsakos, 2000; Kudela andMalecha, 2007), instead of the values of vorticity generated on the wall to satisfy the no-slip condition, one can meet the generation of proper values of the flux of the vorticity. This flux of vorticity is determined by the formula (Koumoutsakos et al., 1994) ∂ω ∂n =− γ ν∆t (6.5) where γ= [[us]]means the intensity of the vortex sheet along the solid bounda- ry that is equal to the jumpof the velocity from zero to the spuriously induced velocity on the solid boundary.We checked bothmethods and obtained nearly the same results. Fig. 13. Changes of global vorticity for the calculated flow over the cylinder 7. Conclusion From the numerical results that are presented in the paper, one can conclude that theVICmethodwith conformalmapping and the algorithm for determi- nation of the far field boundary condition is effective and accurate. The VIC Vortex in cell method for exterior problems 795 method with generation of vorticity on the boundary to satisfy the no-slip condition is also accurate and able to reproduce a very fine vortex structure that arises at the solid boundary. The vortex particle methods are especially well suited for the problem governed by the vorticity and vortex structures. That situation takes place when we study flows over solid bodies. Due to the fundamental meaning of vorticity in fluid mechanics, the importance of the vortex particle method in numerical fluidmechanics arises as well. References 1. Anderson Ch.R., Reider M.B., 1996, AHigh order explicit method for the computation of flow about a circular cylinder, J. Comp. Phys., 125 2. Becerra-Sagredo J.T., 2003, Moment conserving cardinal splines interpo- lation of compact support for arbitrarily spaced data,Research Report No. 10, Zurich, Switzerland 3. Bouard R., Coutanceau M., 1980, The early stage of development of the wakebehind an impulsively started cylinder for 40