An efficient two-dimensional Vortex method with long time accuracy

SIAM J. NUMER. ANAL. Vol. 33, No. 4, pp. 1425-1450, August 1996 () 1996 Society for Industrial and Applied Mathematics 008 AN EFFICIENT TWO-DIMENSIONAL VORTEX METHOD WITH LONG TIME ACCURACY* IBRAHIM BLESS RANERO AND TOM/S CHAC(3N REBOLLO Abstract. This paper deals with efficient techniques for the numerical solution of two-dimensional free-space incompressible Euler equations. We develop an algorithm for fast computation of velocity in a vortex method based upon discretization of vorticity by finite elements. We prove that the method with fast computation of velocity is numerically stable and convergent with second-order accuracy. Some standard numerical tests show that the algorithm with Delaunay regridding bears good stability and accuracy properties for long integration times, with a relatively low computational cost. Moreover, the algorithm is found to be more accurate than high-order vortex-blob methods with regridding for long enough integration times. Key words, vortex methods, finite elements, Delaunay regridding, long time accuracy AMS subject classification. 65N30 1. Introduction. This paper deals with the numerical solution of two-dimensional freespace incompressible Euler equations by means of vortex methods with finite elements. Lagrangian methods, and in particular vortex methods, have many favorable practical properties for the numerical simulation of incompressible flows at a high Reynolds number (cf. Leonard 18], Majda 19] for a detailed bibliography and Anderson 1 for flows in bounded domains). Vortex methods essentially introduce no numerical viscosity and are quite accurate and stable, at least for short times (cf. Beale and Majda [5], Perlman [20]). A complete theory of stability and convergence of vortex methods for free-space twoand three-dimensional Euler equations has been developed since the late seventies. In Anderson and Greengard [2] and Beale and Majda [3], [4] an analysis of convergence in/P-norms, with finite p, may be found. Lately, a theory of convergence in/-norms was developed, mostly by Hou and his collaborators (cf. Hou and Lowengrub [16], Hou [14], Hou [15]). Those works prove that vortex methods are essentially stable in pand/-norms, but stability is conditioned to a relatively high order of consistence. Such "conditioned stability" appears in all convergence proofs of vortex methods, as essentially related to the singularity of the B iot-Savart kernel. In practice, vortex methods compute quite accurate solutions of Euler equations for relatively short integration times. Beale and Majda [5] and Perlman [20] have tested vortex-blob algorithms with high-order kernels. Those experiments confirm the theoretical predictions of order of convergence for moderate times; however, at later times the high-order accuracy progressively deteriorates. This seems to be due to the progressive increase of local grid size that, in general, takes place as the initial grid is deformed by the flow. This leads to an increasing loss of accuracy in the computation of velocities. As stability is conditioned to enough accuracy, a progressive loss of stability also occurs. Beale and Majda conclude that to obtain accurate solutions for long integration times, regridding techniques are usually needed. Introducing regridding techniques decreases the local grid size and allows stability and accuracy for longer times. Unfortunately, this introduces increasing levels of numerical diffusion, counterbalancing the main feature of vortex methods. Our purpose here is essentially to develop a vortex method in which it is possible to reduce the local grid size without introducing numerical diffusion at the grid nodes. Our *Received by the editors May 19, 1993; accepted for publication (in revised form) September 28, 1994. This research was partially supported by Spanish M. E. C. Project DGICYT PB91-0619. Departamento de Matemitica Aplicada, Universidad Complutense de Madrid, Avda. de la Complutense, s/n. 20840 Madrid, Spain ([email protected]). Departamento de Ecuaciones Diferenciales y Amilisis Num6rico, Universidad de Sevilla, C/Tarfia, s/n. 41080 Sevilla, Spain ([email protected]). 1425 1426 IBRAHIM BLESS RANERO AND TOM/S CHACON REBOLLO work is based upon a vortex method introduced in Chacon and Hou 10]. In this method the vorticity is discretized by means of piecewise linear finite elements. Also, the mesh points of the initial triangulation are transported along the streamlines of the discrete flow. Thus, the vorticity is accurately computed at the mesh points, since the vorticity is conserved along streamlines. This method shares many nice properties of vortex methods, being in particular nondissipative. The method is proved to be stable and convergent with second-order accuracy in the uniform norm. However, the stability is conditioned to an order of accuracy bigger than one, much as in classical vortex methods. Numerical experiments show that for short integration times this method is quite accurate, but that as time increases, the triangulations generally become degenerated. This produces a fast loss of accuracy in a short time after degeneration. On the other hand, the method uses a technique of computation of discrete velocities that requires an amount of operations of order O (N2), N being the number of grid nodes in the support of the vorticity. These two drawbacks make the method unfeasible in practical cases. This paper reports some modifications of this method that render it accurate for very long times, with a relatively low computational cost. At first, we develop a fast technique to compute the discrete velocity in the Chacon-Hou algorithm. This is an adaptation of the technique introduced in Greengard and Rokhlin 12] to our context of discretization of vorticity by finite elements. This reduces the amount of computational work required by the method to O(N log 2 N). We prove an estimate of error, in terms of the uniform norm of the vorticity, for the computation of discrete velocities by this technique. Furthermore, we prove that if in the Chacon-Hou algorithm the velocity is computed by this technique with enough accuracy, the numerical solution is still convergent in the uniform norm with second-order accuracy. The proof deals essentially with the fact that our algorithm keeps constant the uniform norm of the discrete vorticity. This allows us to obtain uniformin-time estimates for the error in the computation of the velocities with the fast algorithm and ensure the stability of the algorithm. We also prove that if Delaunay regridding is introduced, the algorithm with fast computation of velocity is still convergent with second-order accuracy. Delaunay regridding constructs the triangulation that maximizes the smallest angle of all possible triangulations supported by a given cloud of points. Introducing this regridding technique in our method produces the only effect of redefining the connections between grid points. The values of the discrete vorticity at the grid points remain unchanged. Thus, the local grid size is reduced, without introducing numerical diffusion. The algorithm is convergent independent of the actual time-stepping strategy used to apply Delaunay regridding. We finally report some numerical experiments dealing with test cases considered by Beale and Majda. We confirm our theoretical expectations on the number of operations in the computation of discrete velocities. We also perform some tests of the modified ChaconHou algorithm, introducing Delaunay regridding. These tests show excellent properties of stability and accuracy, for very long time intervals, in the test cases considered. The errors are conserved close to the initial values, and the theoretical convergence orders are confirmed numerically, even for long integration times. We also observe that this algorithm compares advantageously to a desingularized vortex method of the same order, and also to high-order vortex-blob methods with regridding for long enough times of integration. Our paper is organized as follows. In 2, we describe the algorithm reported in Chacon and Hou 10]. In 3 we develop the technique for fast computation of velocities. Section 4 is devoted to the analysis of the convergence of the modified Chacon-Hou algorithm. Finally, in 5 we report our numerical tests. LONG TIME ACCURATE VORTEX METHOD 1427 2. Description of the base algorithm. In this section we shall describe the vortex method with finite elements introduced in Chacon and Hou 10], as some relevant properties of this method motivate our work. We are interested in the numerical solution of free-space Euler equations in two space dimensions, with a homogeneous condition at infinity. In vorticity (co)-streamfunction (q) formulation, these equations are COt -[- (U V)co 0 in Rex]0, T[, co(x, 0) co0(x) in R e, (1) -Aq co inR 2, lim q(x) 0. Here, u (u l, u2) is the velocity field, defined by the two-dimensional Biot-Savart law, (2) u(x, t) (K co)(x) [ K(x y) co(y, t) dy, where K is the B iot-Savart kernel, (3) g(x) 2zrlxl 2 (-x2, Xl). Equations (1) are equivalent to the usual formulation of Euler equations (cf. Anderson and Greengard [2], Kato [1 ]). The equation for the vorticity in (1) may be integrated exactly on the streamlines X (t; s, associated to the velocity field u. The curve 6 [0, T] --+ X (t; s, or) 6 R 2 describes the trajectory of a fluid particle whose position at time s is the point ot 6 R 2. It satisfies the following ordinary differential equation: -d--(t; s, or) u(X(t; s, ), t) in [0, T], (4) X(s; t,) =. Then, (5) co(X (t; s, or), t) co(or, s) ’ ot R 2 and V s, in [0, T]. Chacon and Hou show that if co is approximated by a piecewise polynomial function on polygons, then the corresponding velocity given by (2) may be computed analytically. This, together with (5), suggests that we approximate the vorticity by finite elements on a moving grid whose nodes describe streamlines of the flow. The algorithm of Chacon and Hou is based upon this idea. Although its description needs rather complex notation, we shall provide it as needed throughout this paper. Consider a triangulation Th of R 2 with triangles {727}i u and nodes {Olj}j6 u. We shall assume that h denotes the length of the longest side of all triangles of Th. Denote by Vh the space of continuous piecewise affine finite elements on triangulation Th, defined by (6) Vh Vh C 0 (RZ) vhl is affine for all triangles r 6 Th }. A function Vh Vh is uniquely determined by the values Vh (otj) ’i iV. 1428 IBRAHIM BLESS RANERO AND TOM/S CHACON REBOLLO Assume that we know an approximation Xj to the point X (tn; 0, otj) for each j 6 N. Define the approximate function X 6 Vh to the flow map X (t; 0, .), respectively, in such a way that it verifies ^n ’j 6N. (c s) x ^n Then, 7 {fi (r/) ’4 r/ Th} defines a triangulation of R e, provided the triangles r do not overlap. If this is the case, we also define the piecewise affine finite element space Q on triangulation 7 in the same way that Vh is defined on Th in (6). We shall denote by/3 (x; {coj }) the canonical interpolation operator on P, defined by We shall also denote by di)j }jEN the canonical base of V h Each (Ioj is uniquely defined by ifj k, j(Otk)-- 0 ifjCk. We are now ready to state the Chacon-Hou algorithm. ALGORITHM A. Suppose that coo is of compact support. Denote wj co0(otj). 1. Initialization (i) Triangulation; (ii) Vorticity; (iii) Velocity; ^0 X =otj and T=Th. 2. Time iteration (i) Update Lagrangian mesh points by the second-order Adams-Bashforth method, ^n+l "n At [ ^n(]) ^n-l(ry-1 ] Xj Xj -- - 3 l,t h --l,t h (ii) Construct a piecewise linear approximation to X (t; 0, .), ^n+l j n+l Z Xj jEN (iii) Update vorticity; ()+1 (X) /3+1 (X, {O)j }). (iv) Update velocity; ^n+l +1 u h =K*d Algorithm A may be viewed as a vortex-blob method with a cutoff function varying in space and time and nonoverlapping smoothing parameter 3 h. Thus, Algorithm A shares with LONG TIME ACCURATE VORTEX METHOD 1429 vortex methods the property of being nondissipative. Furthermore, it is convergent without overlapping the smoothing "blobs." Chacon and Hou prove that Algorithm A is convergent in the uniform norm with secondorder accuracy under some regularity properties of the family of initial triangulations. Accuracy of order higher than one appears to be crucial to ensure the stability of the method. Some numerical tests show that this method is accurate only for relatively short time intervals, much as are vortex-blob methods (see Perlman [20]). However, at later times the triangulations become degenerated and accuracy is progressively lost in a short time. Thus, regridding techniques are needed to increase the interval of accuracy of the method. Chacon and Hou prove that the method is convergent for longer time intervals if some specific regridding techniques are introduced. This occurs in particular if the grid nodes are regrouped to form new triangulations, so that their regularity is conserved. 3. Fast computation of discrete velocities. The method of computation of discrete velocities introduced by Chacon and Hou requires computations of order O(Ne), where N is the number of grid points in the support of w h." o This makes the algorithm slow even for moderately large grid sizes. In this section we develop an adaptation of the algorithm introduced in Greengard and Rokhlin [12] (cf. also Greengard [13]) that computes an approximation to the discrete velocities by means of Taylor expansions. This reduces the computational complexity to O (N log e N). We omit most of the proofs of the results presented in this section, as they are adaptations of the corresponding ones given by Greengard and Rokhlin. Let Th be a triangulation of R e with nodes {otj }jsN. We assume a certain uniformity in the spatial distribution of the grid nodes of Th. This is needed to ensure the convergence of the interpolation by finite elements (cf. Ciarlet 11 ]). We are given a vorticity oh Vh, where Vh denotes the space of piecewise affine elements on Th defined in (6). Moreover, we assume that h has compact support. 3.1. Truncated expansions. Our first goal is to obtain a truncated Laurent expansion that approximates the far-field velocity induced by the vorticity supported by a given subset Q of R 2. This vorticity is given by (7) &Q (i (I)i. oti a Define the set Qh-- {x R2 / d(x’ Q)-- inf lxyl <-h } Then, supp ()Q may include Q, but in any case supp OQ C Qh. Denote by/Q (b/l, b/2) the velocity induced by &O" f K(x x’) OQ(X’) IQ(X) Jsu PPffQ Then, H ul iu2 is a function of complex variables, analytic in C\ supp &a, as the real and imaginary parts of H satisfy the Cauchy-Riemann conditions in C \ supp &a. Moreover, (8) /(Z) ppffQ ZZ Q(xt)dxt’ where Z Xl "3r" ix2, Z X X -qlX2, (Xtl X2). 1430 IBRAHIM BLESS RANERO AND TOMAS CHACON REBOLLO This allows us to expand//in a Laurent series as follows. THEOREM 3.1. Assume supp _Oa ( D(zo, r), where D(zo, r) denotes the complex disc of center zo and radius r > O. Then, (9) 5/(z) (Z Z0) k+l V z C \ D(Zo, r), k=0 where (10) 1 fsu (z’ zo)koQ(X ’) dx’. ak / PPffQ Moreover, given p > 1, let us denote by Ltp (z) the p-term truncated expansion (9). Then, the following error estimate holds: (11) I(z) -p(z)l _< Iz zol-r Iz z01 if Iz z0l > r, where 1 f Ial. (12) A --- PPffa The coefficients a given by (11) may be computed analytically by means of a complex version of Green’s theorem. Indeed, as supp (Q is a reunion of triangles r of Th, for each such triangle there exist br, cr, and dr such that 1_ 1 ()QIr (371’ 372) br xl -b cr x2 qd Pr Z + Pr qdr, where pr = br + cr. Thus, 1 r C supp 5 Q z (z zo) g z dz + - pr (z zo) dz + dr (z zo) k dz The mentioned complex version of Green’s formula is used to compute the integrals in (3.1). It is stated as follows. LEMMA 3.2. Let B be a bounded measurable subset ofR 2 with Lipschitz boundary 0 B. Let f and g be two functions of complex values defined and analytic in some open set containing B U 0 B. Then, f 1 f f (z) g(z) dz. f’z gz dz Let us now consider a triangle r Th included in supp ()Q. By taking f (z) z, e,(z) (z zo) z, we obtain 1 f ) (z zo) z dz z (z zo dz. LONG TIME ACCURATE VORTEX METHOD 1431 Also, by taking Z 2 f(z) g(z) (z Zo) 2’ fo2 f (z zo) tz z (z zo) clz. A similar choice allows us to find the last term in (3.1). Finally, the computation of the coefficients a reduces to that of polynomials on segments of straight lines, which is obtained analytically. Our purpose now is to shift the centers of the truncated Laurent expansion obtained above and to convert these expansions into Taylor expansions. We do this in the next three lemmas. LEMMA 3.3. Assume supp Oa Q D(zo, r). Let zl C be such that Iz0 zi[ > r. Then, bk (13) (z) = (ZZl) k+l k=0 where Y z C \ D(zl, Izl z0l + r), (14) b Z at(z ZO) k-l. /=0 Moreover, if we denote by ld the p-term truncated expansion (13), the following error estimate holds: A [,Zl-ZOl+r] p+I (15) I/d(z) -/d(z)l _< Iz zl- (]Zl zol + r) Iz Zol where The transformation of Laurent expansions into Taylor expansions is made as follows. LEMMA 3.4. Assume supp &a C D(zo, r). Let z C be such that IZl z0[ > (c + 1)r (16) (z) = /k(ZZl) k /Z D(Zl, r), k=0 where (k+l) at (17) (-1) (z zo) + o= (zl zo) l" Moreover, given an integer number q > 1, define the truncated expansion p (18) Upq(Z) /k (Z Zl) k, k=O where the coefficients , are defined by k al (19) = (-1) (zz0) k+l (Zl -zo) l" A 2rr ppQ for some c > 1. Then, 1432 IBRAHIM BLESS RANERO AND TOMAS CHACt3N REBOLLO Then, the following error estimate holds: (20) [(z) LCpq(Z)] r c-1 c-1 +1 where 1 fu IQI, A PPQ Z G D(Zl, r), Lemma 3.4 provides an error estimate for a truncated expansion l/[pq that is defined only with a finite number of data: the p + q coefficients a l, a2 ap+q furnished by either Theorem 3.1 or Lemma 3.3. The translation of the centers of the truncated Taylor expansions is given as follows. LEMMA 3.5. Given the complex numbers Zl, Z2; /91,/91 pp, the following holds: p p (21) & (z z) /Sk (Z Z2) t, k=0 l=0 where (22) p /k /gk (Z2 Zl) k-l. k=l Observe that as formula (21) is exact, the error bound (20) is still true for the truncated Taylor expansion with shifted center if z2 6 D(Zl, r) and z D(z2, r ]Zl z21). Also, assume that the p are the coefficients of the Taylor expansion of/g around z. Then, in general, the/5l are not the coefficients of b/around z2. This would happen only if the sum in (22) were infinite. 3.2. Description of the algorithm. To describe to some extent the algorithm for fast computation of the discrete velocity, we shall need some specific definitions. We assume that the support of the vorticity ffh is included in a square of sides of length H, = 2’2 x 2’2 We refer to this square as the computational domain. We subdivide the computational domain into a family of boxes of decreasing size which will be linked by a hierarchy relation. Let N be the number of nodes Oem in supp (Sh, and define the highest level of refinement: J log 4 N. Given a level of refinement j 0, 1 J, we denote by Q, k 1, 2 4J the squares /-/ We denote by obtained by splitting each side of into 2 j subintervals of length hj 2-7. /Jk the center of box Q. We assume that each box Q at level j is a Cartesian product of intervals of the form Q [a, b[ x [c, d[ H for some a < b and c < d in [-7, [" Then, the boxes at level j do not overlap, but their reunion is the whole computational box LONG TIME ACCURATE VORTEX METHOD 1433 DEFINITION 3.6. Given d > O, we shall say that two boxes Q and Q Jr of the same level j are d-separated if k --i11 " d hj. The fact that two boxes are well separated allows us to approximate uniformly on each one of them the velocity field induced by the vorticity supported by the other one by means of the truncated Taylor development furnished by Lemmas 3.4 and 3.5. To state this result we need to give a formal definition of the aspect ratio of a triangulation, as a measure of its regularity. The aspect ratio of a given triangle is the ratio between the diameter of the smallest circle that can circumscribe the triangle and that of the larger circle that can be inscribed in the triangle. The aspect ratio of a triangulation Th is the largest of all aspect ratios of all triangles of Th. LEMMA 3.7. Given a box Q of level j, denote by co Vh the vorticity supported by QJ (23) Otm E Q Call lt the (p, q)-term Taylor expansion associated to o9 defined by (18). Define Let c > 1 be given. Then, there exists a constant ) > O.depending only on the a.spect ratio of triangulation Th such that if d (c + 1) and box QJ is d-separated from QJ, then C-+-1 (24) xEQ{max lu uJl _< B hj (c 1)2 + c +1 Ilffh where B is a constant depending only on the aspect ratio of triangulation Th. Proof. As h is the longest length of all triangles of Th, then (25) supp coJ C B(flJ, rj) with rj -- hj + h. Also, as h j , there exist two positive constants v and/z depending only on the aspect ratio of T, such that vh<hj<tzh. Thus, 1 1 rj < ----_ hj +-hs <_. )hj, / 2 As the boxes are d-separated, 1 1 where ) + -’v d fl fill> d h j > r. If we take d ) (c + 1), we have the hypotheses of Lemma 3.4 with zo fl J k’ Zl [l r rj. Estimate (24) follows immediately. 1440 IBRAHIM BLESS RANERO AND TOM,S CHACt3N REBOLLO FIG. 5. Numerical convergence order in velocity of Algorithm A’, applied to TC 1, for a short time interval. The lines marked by and + symbols correspond, respectively, to hi 1/12, h: 1/16 and to hi 1/16, h2 1/20 to compute the numerical convergence order. The corresponding theoretical convergence order is p 2 (line marked by 0 symbols). FIG. 6. Numerical convergence order in velocity of Algorithm A’, applied to TC 1, for a short time interval. The lines marked by and + symbols correspond, respectively, to h / 12, h2 / 16 and to h / 16, h2 1/20 to compute the numerical convergence order. The corresponding theoretical convergence order is p 1.5 (line marked by 0 symbols). In Figs. 5 and 6 we represent the behaviour of the numerical convergence order in velocity for TC 1 for short times. Figure 5 corresponds to a theoretical convergence order p 2, obtained by choosing a 2 in Algorithm A’. Figure 6 corresponds to p 1.5, obtained by choosing cr 1.5. The sharp oscillations observed in these curves occur when Deiaunay regridding is performed. In both cases there is a good agreement between the computed convergence order and the theoretical one. Furthermore, the curves corresponding to smaller values of h are globally closer to the theoretical convergence orders. Note that in all cases the computed orders oscillate around the theoretical value. It is interesting to observe that the convergence order of Algorithm A’ may be preset to any value p 6] 1, 2] by simply choosing the parameter cr p at the beginning of the run. However, taking r smaller than 2 would produce a waste of computational work, as the computational complexity of Algorithm A’ is of order N log 2 N for any value of a. In what follows we shall always take tr 2. LONG TIME ACCURATE VORTEX METHOD 1441 rrors FIG. 7. Time evolution of numerical convergence order in vorticity of Algorithm A r, applied to TC2, for the time interval [0, 15]. The lines marked by +, *, and 0 symbols correspond, respectively, to hi 1/12, h2 1/16; hi 1/16, h. 1/20; and to hi 1/16, h2 1/20 to compute the convergence order. A good agreement with the theoretical prediction p 2, which improves for smaller values of h, is observed for the whole time interval. 3"0 1 2.5 Errors FIG. 8. Time evolution of numerical convergence order in velocity of Algorithm A , applied to TC2, for the time interval [0, 15]. The lines marked by +, *, and 0 symbols correspond, respectively, to hi 1/12, h2 1/16; hi 1/16, h2 1/20; and to hi 1/16, h2 1/20 to compute the convergence order. A goodagreementwith the theoretical prediction p 2, which improves for smaller values of h, is observed for the whole time interval. Figures 7 and 8 show the numerical convergence order in vorticity and velocity for TC2 during a relatively long integration time. The sharp oscillations of the former test are still observed. However, again the curves corresponding to smaller values of h are closer to the theoretical order p 2. Note that the numerical orders in velocity are closer to p 2 than those in vorticity. This is probably due to the higher regularity of the velocity field. Note also that the agreement holds even for long times. 5.3. Behaviour for long integration times. Our third set of numerical experiments deals with the analysis of the long time behaviour of errors due to Algorithm A’. At first, we tested the effect of introducing Delaunay regridding in Algorithm A’. In Fig. 9 we represent the time evolution of the percent relative errors in velocity corresponding 1442 IBRAHIM BLESS RANERO AND TOM/S CHACON REBOLLO o Errors FIG. 9. Time evolution of relative errors in vorticity of Algorithm A with and without using Delaunay regridding (marked by and + symbols, respectively). The second one grows exponentially while the first one keeps close to the initial error. to Algorithm A without regridding and to Algorithm A’ using Delaunay regridding, applied to TC with h / 12. We made a run in the time interval [0, 50], which appears as quite a long time for the test case considered. Indeed, in that interval the fastest points in supp 09 completed more than three rotations, while the slowest ones did not complete half a rotation. Note that Algorithm A is still defined when the triangulation TT, becomes degenerated. Indeed, in this case it is still possible to compute directly the exact velocity T. We may observe that in this case the error grows exponentially until attempting relative values of more than 50% at _ 50. When Delaunay regridding is used, there is a fall of error each time it is effectively performed. This is due to the diminishing of the local grid size that renders the linear interpolations on each triangle more accurate. After this, there is a slow exponential increase of errors, which falls again the next time Delaunay regridding is used. Note that regridding is not needed very often. The error at 50 is approximately only two times the initial one. Figures 10 and 11 also show the behaviour of errors in velocity and vorticity corresponding to TC1 and TC2, respectively, with h 1!12, during the time interval [0, 100]. This is a very long time interval for both cases, as at time 100 the unit circle has been dramatically deformed by both flows. The curves present sharp oscillations due not only to the use of Delaunay regridding, but also to the low smoothness of the discrete/-norm used. However, we remark at first that in both cases the errors in velocity and vorticity remain almost constant. Also, in both cases the errors in velocity are substantially smaller than those in vorticity. Again, this is very probably due to the higher smoothness of the velocity fields. Note also that although in TC2 the vorticity is not smooth enough to ensure the convergence of Algorithm A’, in practice second-order convergence is attempted. A possible reason for this fact is that the singularities of the vorticity lie on the curve r 1, while we solve Euler equations only inside the unit circle. Also, the fact that the errors corresponding to TC2 are close to those corresponding to TC 1 is probably due to the radial distribution of the nodes in the triangulation used. Figure 12 represents the triangulation for TC1 at time 99, the last time Delaunay regridding is used. Observe the good quality of the grid, which suggests that our run could continue for longer times with similar error levels. Globally, these tests show that Algorithm A’ with the use of Delaunay regridding is stable and accurate, with second-order accuracy, even for very long integration times. LONG TIME ACCURATE VORTEX METHOD 1443 Errors 2.5 2.0 Z. 5 FIG. 10. Time evolution of errors in velocity (line marked by + symbols) and vorticity (line marked by symbols) for Algorithm A’, applied to TC1 with h 1/12, in the time interval [0, 100]. Both errors remain close to the initial values for the whole time interval. Errors 2.5 2.0 0 9 18 27 36 45 54 63 72 81 90 99 FIG. 11. Time evolution of errors in velocity (line marked by + symbols) and vorticity (line marked by symbols) for Algorithm A r, applied to TC2 with h 1/12, in the time interval [0, 100]. Both errors remain close to the initial values for the whole time interval. 5.4. Comparison to a desingularized vortex method. Our next experiment is to compare the performances of a vortex method on a fixed uniform grid with those of Algorithm A’ using Delaunay regridding. As the vortex method we have used the desingularized point vortex method (DPVM) introduced in Hou [15]. To describe it, let us consider a uniform grid of size h of R 2 with nodes {flj }j EN. Denote coj COO (/j). Then, DPVM computes the discrete velocity at point/3 and time by (43) h(flk, t) g(k l) (COl (-Ok) COk I g( k y) dy, /l Esupp wo Jf2h (t) where ]h (t) is a polygonal approximation of supp co(., t). 1444 IBRAHIM BLESS RANERO AND TOM,/tS CHACON REBOLLO FIG. 12. Triangulation at time 99 corresponding to Algorithm A t, applied to TC2 with h 1/12, in the time interval [0, 100]. This is a stable modification of the point vortex method, uniformly convergent with second-order accuracy. Because of that, it seems to be a good method to compare with ours. In our experiments we have run both algorithms for grids of size h 1/12 and h 1/16. We always take the unit circle to be the set 2h(t). This allows us to compute exactly the integral expression in (43). Also, we solved the equation of characteristics for the DPVM with the Adams-Bashforth second-order scheme, just as in Algorithm A’. In our tests, if no regridding techniques are introduced, the DPVM produces a large increase of errors in a relative time interval. For instance, for TC the relative errors in velocity take values of approximately 80% by time 40. As we pointed out in the Introduction, some convenient regridding technique is needed to obtain accurate solutions for long integration times. Beale and Majda introduced in 1985 a simple, but efficient, regridding technique in the context of a vortex-blob method (VBM). VBMs are based upon the discretization of vorticity as a sum of smooth functions with small supports, called blobs. Given a smooth cutoff function q (i.e., an approximation of the Dirac delta at the origin), the vorticity at a fixed time is approximated by (44) co(x, t) coh(X, t) qa(x Xj(t)) COj h 2, J where 1 %(x) With a discretization of the kind of (44), it is possible to compute the vorticity at any prescribed point. The regridding technique of Beale and Majda consists of reinterpolating the vorticity at the nodes of a uniform grid, whenever the current local grid size is large enough. However, as reported by Beale and Majda, the grid size of the reinterpolating grid should decrease progressively to maintain reasonable error levels. In DPVM the vorticity is discretized as a sum of Dirac masses: (45) co(x, t) -- cob(X, t) Z 6(X Xj(t)) coj h 2. J LONG TIME ACCURATE VORTEX METHOD 1445 4 0 9 18 27 36 45 54 69 72 81 90 99 FIG. 13. Comparison of errors in velocity between the DPVM (line marked by + symbols) and Algorithm A (line marked by symbols), applied to TC1, in the time interval [0, 100]. A progressive increase is observed in the first one, while the second remains almost constant. Consequently, the regridding technique of Beale and Majda cannot be directly applied here. However, it is possible to use this technique after approximating d by a sum of blobs as in (44). In practice, we have used a fourth-order cutoff function : *8(x)=- 2exp -- - This cutoff function was also introduced by Beale and Majda in 1985. For smooth functions f, the error f a, f is of order 64. We have taken 6 of order h, so this accuracy seems to be enough, as the DPVM is of order h 2. The regridding strategy that we have used consists of reinterpolating the vorticity when the smallest angle of the deformed grid is smaller than a preset limit value. Each regridding has been set to produce an increase in the amount of the grid points of approximately 15%. Figures 13 and 14 compare the behaviour of errors due to the DPVM and to the finite element vortex method (FEVM) of Algorithm A’. We represent the errors in velocity and trajectories corresponding to TC1 during the time interval [0, 100] with initial grid size h 1/16. Sharp variations of errors corresponding to the DPVM are observed, probably due to the additional error introduced in the reinterpolation associated to the regridding steps. The errors corresponding to the DPVM increase faster than those corresponding to the FEVM. By time 0, the errors corresponding to both methods take very similar values. By time 100, the later errors are nearly 200 times smaller than the former ones in velocity and nearly 20 times smaller in trajectories. This different growth rate is probably a consequence of the introduction of numerical diffusion in the regridding steps. We may conclude that the FEVM solves more accurately our TC 1 for long times, without introducing numerical diffusion. We should point out that the computational work needed by one time-step with the FEVM is nearly 100 times bigger than the one needed by one timestep with the DPVM. However, this work remains constant in time for the FEVM, while that 1446 IBRAHIM BLESS RANERO AND TOM/S CHAC(3N REBOLLO 2.5 -0. Errors 1(3 27 3(3 45 54 63 72 131 90 99 FIG. 14. Comparison of errors in trajectories of grid points between the DPVM (line marked by + symbols) and Algorithm A (line marked by symbols), applied to TC1, in the time interval [0, 100]. needed by the DPVM increases each time regridding is performed. Thus, for long enough time intervals, both computational efforts will be of the same order. 5.5. Comparison to high-order VBMs. We finally compared the FEVM with the VBM with cutoff functions of orders m 4 and m 6. Specifically, we used those introduced by Beale and Majda, corresponding to p:4, (x)=- 2exp - - exp (- 2@2)] p--6, 1 r 2 1 r 2 J6)(x) - I exp (----)- exp (--2@2)+ - exp (---)]. To ensure the convergence of the method, we took the blob size to be 6 h q, with 0 < q < 1. Thus, for smooth enough initial vorticity, the convergence order of the method is p=mq. In practice, we took q 0.95 in all our experiments, as this value seems to be quasioptimal, as reported by Perlman. The streamline equation has been solved with a fourth-order Runge-Kutta method with very small time-step. We have tested our code for TC1 at 1. Our estimations of computed convergence orders are given in Table 1. They are in very good agreement with those reported by Perlman. We have used the regridding technique described in the preceding subsection, with some minor modifications. Indeed, many possible criteria which can be used to apply regridding are equivalent in practice for the rotating steady solutions we are considering. Either regridding when the smallest angle of the mesh is smaller than a given tolerance, or when the current grid size is long enough, is equivalent to regridding a certain fixed number of time-steps. In any LONG TIME ACCURATE VORTEX METHOD 1447 TABLE Convergence order of the VBM for TC1, estimated at time 1, used to compare to the FEVM. Values of m and h h 0.2 h 0.1 h 0.05 Theoretical orders m 4 2.53 3.32 3.60 3.80 m 6 2.78 4.44 5.16 5.70 O,B 0.7 rIME I FIG. 15. Comparison of errors in velocity between Algorithm A’ (line marked by 0 symbols) and the VBM with m 4 (line marked by + symbols) and m 6 (line marked by symbols)for TC1 and h 1/16 in the time interval [0, 100]. case, the actual tolerance value must be tuned with care to avoid an excessive increase in errors. If regridding is applied too often, we shall progressively introduce high levels of numerical diffusion, but if the grid is excessively distorted when regridding, then the accumulated errors will produce an unrecoverable loss of accuracy. In Figs. 15 and 16 we represent the relative errors in velocity for the FEVM and for the VBM withm 4 andm 6, corresponding to TC1 with h 1/16 and TC2 with h 1/12. We may observe that regridding is applied in all cases an almost constant number of time-steps. In all cases errors are kept almost constant for a short time interval whenever regridding is applied and experience a fast increase when the grid becomes progressively distorted. For m 6, this increase is very fast, and this is probably the reason why regridding produces a decrease in errors. For m 4, the errors do not grow as fast, and regridding produces an increase in them. Also, almost linear growth rates of errors are observed in all cases. These rates are smaller for m 6 than for m 4. For the FEVM, the growth of errors as the grid is distorted is the fastest of all cases considered. Thus, the loss of quality of the grid more dramatically affects the accuracy of the FEVM than that of the VBM. However, applying Delaunay regridding in the FEVM diminishes the errors to values close to the initial ones, counterbalancing almost completely the former increase. In TC 1, which corresponds to a smooth solution, the FEVM presents a better performance at time 100 than the VBM with m 4, while the VBM with m 6 yields a higher accuracy than the FEVM. However, in TC2, which corresponds to a less smooth vorticity, the performance of the FEVM at 100 improves that of the VBM with m 4 and also that of the VBM with m 6. We should also remark that the growth rates of errors for the BVM are 1448 IBRAHIM BLESS RANERO AND TOM/S CHACON REBOLLO (VEL.),_ H=/! 4050. 70, 100. FIG. 16. "Comparison of errors in velocity between Algorithm A’ (line marked by 0 symbols) and the VBM (line marked by + symbols) with m 4 and m 6 (line marked by symbols)for TC2 and h 1/12 in the time interval [0, 100]. in all cases larger than those corresponding to the FEVM. Thus, the comparison is very likely to be even more favourable for the FEVM for later integration times. Finally, we must say that we may not expect to solve two-dimensional Euler equations with any initial condition, simply by using Algorithm A’ combined with Delaunay regridding. It seems clear that, in general, the other regridding rules that we mentioned in 2 are needed to obtain accurate results. The results presented in this paper must be understood in the sense that our FEVM, due to its geometrical adaptability, improves the accuracy of classical vortex methods, without introducing numerical diffusion. Appendix: Proof of Theorem 4.1. Proof. Our proof is an adaptation of the convergence proof for Algorithm A given by Chacon and Hou. The essentials of the proof are as follows. Let us define T* max t. "0 < tn < T, max IlX 311,h < h I+p wherep [min(2 cr)-l] 0<k<n We prove that there exists a separation parameter d, depending only on Coo, T, c, and , such that estimates (33) hold in the time interval [0, T*]. Then, we conclude that there exist two positive numbers h Tand ATsuch that if 0 < h < h Tand 0 < At < AT-, then there must be T*>T. The main innovation in our analysis is that we obtain uniform-in-time estimates for the error in the computation of discrete velocities by Algorithm B. As we shall see, this essentially happens because our algorithm preserves the uniform norm of the discrete solution. Indeed, if 0 _< t, _< T*, following Chacon and Hou we state that if h is small enough, then all triangulations {7}0<_t,<_, are nondegenerated. Moreover, all aspect ratios of these triangulations are uniformly bounded from above by a constant }/ independent of h. Let us define the separation parameter dr )Lr (1 + c), where i z is the parameter associated to ?’ given by Theorem 3.11. Note that d7 depends only on coo, T, c, and the aspect ratio of the initial triangulations. LONG TIME ACCURATE VORTEX METHOD 1449 (46) and We prove now that there exist two positive constants r and C such that supp ^n o) h C B(0, r), 0 < tn < T*, (47) Indeed, assume ^n .n. hr T*. max I[fi--(K.(Oh)]( j)l <Cr 0<tn < O<j<M supp tb-I C B(0, rn-1), supp ^n (-o h C B(O, rn) for some positive numbers rn-1 < rn. Let us denote by the Euclidean norm on R 2 and by II the uniform norm on R 2. Theorem 3.11 yields (48) < C1 rn ^n h f ^n %11 -tIg(] Y)I Coh(Y)dY _< C2 rn Ilco011 JB (O,r.) Note that the last inequality here follows because Algorithm A conserves in time the uniform norm of the discrete vorticity. Consequently, IXj^n+l IXjl -At 3 [u h^n(Jy)[ _+_ [tn-l(j]-l)[ _< IXjl + C3 At, j--1 M. Thus, rn < r ro exp(C3 T), 0 <_ tn < T*. The remainder of the proof is a technical refinement of that of Chacon and Hou. We shall omit it here, as it does not introduce any essential innovation. V] Acknowledgments. The authors wish to thank Macarena G6mez Marmol for her valuable help in obtaining the graphic output. REFERENCES [1] C. ANDERSON, Observations on vorticity creation boundary conditions, in Mathematical Aspects of Vortex Dynamics, R. Caflisch ed., Society for Industrial and Applied Mathematics, Philadelphia, PA, 1988, pp. 144-159. [2] C. ANDERSON AND C. GREENGARD, On vortex methods, SIAM J. Numer. Anal., 22 (1985), pp. 413-439. [3] J. T. BEALE AND A. MAJDA, Vortex methods I: Convergence in three dimensions, Math. Comp., 32 (1982), pp. 1-27. [4] Vortex methods II: High order accuracy in two and three dimensions, Math. Comp., 32 (1982), pp. 29-52. [5] High order accurate vortex methods with explicit vorticity kernels, J. Comp. Phys., 58 (1985), pp. 188-208. [6] M. BERNADOU et al., MODULEF A Modular Library of Finite Elements. INRIA, Rocquencourt, France, 1986. [7] R. BOWYER, Computing Dirichlet tessellations, Comput. J., 24 (1981) pp. 162-166. [8] T. E BUTTI(E, Fast vortex methods in three dimensions, in Vortex Dynamics and Vortex Methods, Lectures in Appl. Math., Vol. 28, K. E. Gustafsson and J. A. Sethian, eds., American Mathematical Society, Providence, RI, 1991, pp. 51-66. [9] A fast adaptive method for patches of constant vorticity in two dimensions, J. Comput. Phys., 89 (1990) p. 161.