A high precision algorithm for axisymmetric flow

We present a new algorithm for highly accurate computation of axisymmetric potential flow. The principal feature of the algorithm is the use of orthogonal curvilinear coordinates. These coordinates are used to write down the equations and to specify quadrilateral elements following the boundary. In particular, boundary conditions for the Stokes’ stream-function are satisfied exactly. The velocity field is determined by differentiating the stream-function. We avoid the use of quadratures in the evaluation of Galerkin integrals, and instead use splining of the boundaries of elements to take the double integrals of the shape functions in closed form. This is very accurate and not time consuming. 1991 Mathematics Subject Classification: 76B99, 76M10, 65N50, 76M20, 65D07


Introduction
Computation of axisymmetric potential flow in water passages of a turbine is an essential step for quasi-three-dimensional solution of problems of analysis and design of the cascade of turbine blade profiles (see e.g.[8,3]).The equation for axisymmetric potential flow must be written in the terms of the Stokes' stream-function, ψ, since the boundary conditions for the velocity potential, ϕ, cannot be specified along the water passages, while the boundary conditions for ψ can be easily specified for the entire boundary (see §2).
This equation may be solved numerically using both finite difference and finite element methods.Finite difference methods cannot yield a sufficiently accurate solution even in the case of a curvilinear orthogonal grid precisely fitting the boundaries.The reason for this is that finite difference solutions satisfy the boundary conditions for ψ only at a finite number of points along the boundary (see §2).
Finite element methods have been applied to the solution of the equation for ψ written in Cartesian coordinates and typically using a rectilinear grid.In this case the outer nodal points of the finite element grid are located on the boundary, however the element sides connecting them are not congruent to the corresponding segments of the boundary.This can lead to large errors in the computation of the Stokes' stream-function ψ, especially where there are sharp changes of the boundary.Curvilinear elements (but not coordinates) have been considered for the above reasons in [6], §2. 4, for the special case of isoparametric triangular or quadrilateral elements with quadratic approximation of serendipity class.
Another difficulty with the finite element method arises from the need to evaluate many double integrals, a costly procedure.Accuracy is often traded off for speed by using a Gaussian quadrature method (see e.g.[6], §2.4.3; [7], §3.8): each one-dimensional integral is approximated by a linear combination of samples at a finite set of points, usually 3, in the domain.For a double integral this usually means 9 samples.Typical relative error for the velocity field found by such a procedure is around 4%-5% and, depending on the circumstances, is sometimes as high as 15%.
To overcome the difficulties with finite element and finite difference methods mentioned above and to achieve high precision, consistency, satisfaction of boundary conditions and high execution speed, we propose a procedure described in detail in subsequent sections.
The main feature of this procedure is the use of an orthogonal curvilinear coordinate system p, q.In our case, the equation for ψ is written in this coordinate system.The usage of the orthogonal curvilinear system enables us to define the element geometry directly in terms of p and q.Furthermore, we define the shape functions in terms of two nondimensional parameters which are linearly related to p and q.
One consequence of this is that the sides of the outer elements follow the boundary.Thus, since the shape functions can specify a quadratic distribution, the solution satisfies the boundary conditions exactly.Specifically, ψ is constant along the solid boundaries of the the water passage.At the remaining boundary ψ is defined as follows: it is a given linear function at the entrance and a given quadratic function at the exit (see §2).
Another consequence of this approach is that we can integrate the shape functions analytically, independently of the geometry of the passage, and express the Galerkin integrals as linear combinations of integrals of the shape functions.This technique is accurate and fast.Our experience shows that it is about 10 3 times faster than Simpson's rule at the same level of precision.In addition, since the integrals of the shape functions do not depend on the geometry of the passage, they do not need to be recalculated when the geometry is changed.
The orthogonal curvilinear coordinate system p, q is constructed by solving a two-dimensional flow problem in the same passage.The p-axes are defined as streamlines and q-axes as equipotentials for the two-dimensional flow.
The stream-function distribution for this flow is obtained using an iterated finite difference scheme.Iteration is used to take care of unpleasant edge effects at the approximate boundary given by the finite difference grid.At each step we find the stream-function and adjust the boundary conditions along the approximate boundary for the next step in order to achieve proper boundary conditions along the genuine smooth boundary.The streamlines are obtained from interpolation of the stream-function distribution.The equipotentials are obtained as a family of curves orthogonal to the streamlines.
The relative error for the stream-function found by this procedure is quite high (around 3%).However, the grid lines, thus obtained, look sufficiently nice and the ratio of Lamé functions, constant along the actual stream lines, had they been found exactly, is close to constant along the grid lines.The latter fact is important for an increase of accuracy of the solution.
Once the grid lines, and thus, curvilinear coordinates, are determined, we apply the Galerkin method to calculate the Stokes' stream-function distri-bution.The grid lines are represented by splines, so that double integrals of the shape functions can be evaluated analytically.The Galerkin integrals are written as linear combinations of the integrals of the shape functions.The velocity field is obtained by differentiation of the Stokes' stream-function and, thus, equals the gradient of the potential.
We have performed a comparison with a closed form solution in a case permitting the latter [5] yielding a relative error in the velocity field of about 0.01%, where the curvilinear coordinates were derived from exact equipotentials and streamlines.In addition, we have performed internal consistency checks in many practical cases, indicating a relative error in velocity ≤ 0.490%.See §5 for details.

Formulation of the problem
Let (r, θ, z) denote the cylindrical coordinates for R 3 .We make the following assumptions about the flow (i) the flow has a potential ϕ, (ii) the flow is axisymmetric with respect to the z-axis, (iii) the flow is meridional, (iv) the flow enters from radial infinity r = +∞, where the passage consists of two planes that are perpendicular to the z-axis, (v) the flow exits at z = −∞, where the passage is either a circular cylinder or a cone.
The problem for a flow with the above assumptions cannot be formulated as a Dirichlet problem for the potential ϕ, because given a potential difference between the entrance and the exit, we do not know a priori the potential distribution along the solid boundaries of the passage.The potential ϕ can be formulated as a solution of a von Neumann or a mixed problem.However, the solution of a von Neumann problem offers less precision unless the method of singularities is used.The method of singularities, while shown to work in certain idealized situations, like flow around a sphere [1], has not been developed for problems that occur in practice.
On the other hand, the problem can be formulated as a Dirichlet problem for the Stokes' stream-function ψ.The boundary conditions are as follows: (i) ψ is constant along the solid walls of the passage, (ii) at the entrance (r = +∞) ψ is a linear function of z, (iii) at the exit (z = −∞) ψ is a quadratic function of r, Boundary condition (i) implies that the normal component of the velocity where is arc length along the boundary.Since a finite difference solution satisfies the boundary condition (i) only at finitely many points, v n may well vanish nowhere at the boundary.In fact, ψ may even be discontinuous.Thus, finite difference techniques cannot provide accurate satisfaction of the boundary condition.
Finite element solutions will satisfy the boundary conditions at the boundary of the grid.If the shape of the elements is arbitrary, then the the boundary of the grid does not coincide with the true boundary of the water passage.Thus, while the boundary conditions are satisfied at the boundary of the grid, they are not necessarily satisfied at the true boundary.In our method we solve this problem by finding an orthogonal curvilinear coordinate system and using it to define elements in such a way that the boundary of the finite element grid coincides with the true boundary of the water passage.
Let (p, q, θ) be such a system, i.e. an orthogonal curvilinear coordinate system, such that θ is the same as in cylindrical coordinates and the solid boundaries of the passage are curves of constant q.Axisymmetry implies that derivatives with respect to θ are 0.For meridional flow v θ = 0. We take this into account and express various differential operators in the (p, q) coordinate system using Lamé metric functions h p , h q , h θ = r (see e.g.10.51(5), 10.611 [4]).Due to the absence of sources we have We can write the components of velocity v as Using the equality of the second mixed partial derivatives of ψ we obtain i.e. the flow is irrotational, so Letting c p = h p /h q , c q = h q /h p we obtain Using the expression for the Laplacian of ψ we transform (1) into 3 Grid generation and splining In this section we describe the construction of an orthogonal curvilinear coordinate system (p, q) for the water passage.We use a finite difference scheme with boundary corrections to solve the two-dimensional Dirichlet problem for the stream-function ψ 2 (x, y).The (x, y) coordinates of the twodimensional problem correspond to the (r, z) coordinates of our axisymmetric setup.The goal here is not to obtain a particularly accurate solution, but to produce sufficiently nice streamlines, which serve as lines of constant q.The orthogonal family of curves is constructed separately.We use interpolation and splining to obtain locally polynomial representations of curves of constant p or q and expressions for various quantities related to the (p, q) coordinates, such as partial derivatives and Lamé functions.In the next section the grid lines become the basis for the partition of the domain into elements.

Domain:
The idealized water passage from r = +∞ to z = −∞ is truncated at the entrance and the exit for practical purposes.At the entrance the passage is truncated by a cylinder and at the exit by a plane, if the passage is cylindrical, or a sphere, if the passage is conical.The position of truncating surfaces was determined empirically by finding where the dependence of the results on the position becomes negligible.
Let Ω be the domain produced by sectioning the truncated passage with the r-z plane.We assume that Ω is a curvilinear solid quadrilateral with smooth edges and use cartesian coordinates (x, y).

Finite differences:
We use a two-dimensional finite difference scheme with a regular cartesian grid.Since we are dealing with the Laplace equation in two dimensions we may solve the finite difference scheme using a widely used iterative procedure rather than solving a linear system.We start by generating a square grid with spacing 2∆h over the domain and setting the values at the boundary nodes.Then we interpolate in the horizontal direction to obtain a first approximation to the values at the nodes.
The iterative procedure for the finite difference scheme is based on the following observations.Consider third order Taylor approximations of ψ 2 at (x, y): (∆y) 3  6 +O (∆y) 4 .
Letting ∆x = ∆y = ∆h, adding the above expressions, and solving for ψ 2 (x, y) gives Since ψ 2 is harmonic, we see that the average of the values at the four corners of a square is a third order approximation for the value at the center of the square.Thus, at each step of the iterative procedure we average the values of four nodes for each square of the grid and obtain approximations for the values at the center of each square.Then we average the results at four neighboring centers, also forming a square indicated by dotted lines in Figure 1, to obtain new approximations for the values at the nodes.

Boundary correction:
If the boundary conditions are applied at the nodes of the finite difference grid near the boundary and not at the boundary itself, then certain unpleasant effects are observed for the streamlines near the boundary.We correct this problem by occasionally adjusting the values at the boundary nodes of the grid using a second order Taylor expansion in the horizontal direction.After the difference scheme converges, we perform the correction and run the scheme again.This is done several times until the overall scheme converges.In the following paragraphs we describe the correction procedure.
Consider a boundary grid node C and nearby nodes A and B with the same y coordinate.Suppose the true boundary intersects the segment CA at the point L (see Figure 1).Then We solve this system for the first and second partial derivatives of ψ 2 with respect to x at L and use these approximate values to obtain The boundary value is substituted for ψ 2 | L in the above formula and the resulting ψ 2 | C is used to set the value at the boundary node C.

Curvilinear coordinate grid:
We construct a curvilinear coordinate grid such that the streamlines, i.e. the lines of constant ψ 2 , are precisely the lines of constant q.Subsequently we construct an orthogonal family of curves with minimal arc length and distributed roughly regularly according to the two-dimensional potential ϕ 2 .We take these curves to be lines of constant p.
At each level y the function ψ 2 (x, y) is represented by cubic splines using Hermite interpolant polynomials.The splining technique is described in detail in [2].The ingredients, besides the values of ψ 2 at the nodes, are the partial derivatives of ψ 2 at the endpoints and the condition that the second derivatives should match at the internal nodes.The partial derivatives of ψ 2 at the boundary are estimated using the values of ψ 2 at a few nodes and Taylor's formula of appropriate order.In this algorithm we used five points.
We use a regular partition of the interval between the boundary values of ψ 2 and at each level of the square grid find the x coordinates of points with a given value of ψ 2 .A linear search finds the neighboring grid nodes such that the given value of ψ 2 is in the interval between the values of ψ 2 at the nodes.Linear interpolation between the two nodes gives a seed to Newton's method.
Once a set of points with the same given value of ψ 2 is found, the corresponding streamline is represented by splines using Hermite interpolant polynomials in terms of the arc length of the curve (see [2]).The partial derivatives at the entrance are At the exit, in the case of a cylinder, we have and, in the case of a cone, the derivative of (x, y) with respect to is normal to the bounding sphere.
Obtaining a set of curves orthogonal to the streamlines is a slightly more intricate procedure.The first goal is to obtain a reasonable partition of a streamline.We start with a regular partition of one of the boundary streamlines.At each point of the partition we draw a perpendicular and calculate the distance along it to the next streamline.The distances and the values of ψ 2 are used to find the distribution of velocity along the streamline and the two-dimensional potential ϕ 2 .The potential interval from the entrance to the exit is partitioned regularly.
The second goal is to construct orthogonal curves from the boundary streamline to the next streamline.This is achieved iteratively to obtain approximately minimal arc length along each orthogonal curve.The procedure is repeated to continue the curves to other streamlines.The derivatives at the boundary used to construct the spline representation are estimated from Taylor's formula and values at several point similarly to what was done previously in this type of situation.
Figure 2 shows a curvilinear coordinate grid generated during a practical application.

Two-dimensional splining:
In the next section we require spline representations for various quantities involving the curvilinear coordinates in all of the domain.This is particularly useful for a fast accurate computation of the Galerkin integrals.Given a coordinate grid (p n , q n ), 1 ≤ n ≤ N , 1 ≤ m ≤ M , a function f is splined from its sample values and those of its partial derivatives in each grid quadrilateral using the following formula: where

Finite element scheme
Finite element methods, also known as Ritz-Galerkin methods, are based on the idea of approximating a functional linear space by a finite dimensional subspace generated by a basis of locally supported functions, called shape functions.The domain is partitioned into subsets, called elements.The shape functions are defined on each element.An approximate solution to a given problem is determined by a linear system for the components with respect to the basis of shape functions for the approximating subspace.
In this section we construct a finite element scheme using curvilinear coordinates (p, q), based on a conventional scheme for cartesian coordinates, and derive the formulae for the coefficients of the linear system for our particular problem.The domain in the r-z plane is partitioned with a grid into curvilinear quadrilaterals, whose sides are lines of constant p and q.We define elements, local interpolating shape functions on them and generate an approximating space G with interpolating global shape functions g n (p, q), 1 ≤ n ≤ K.

Domain:
In the preceding section we truncated the passage at the entrance and the exit and constructed orthogonal curvilinear coordinates (p, q) for the domain Ω produced by sectioning the truncated passage with the r-z plane.
Recall that the edges of Ω are curves of constant p or q with p varying from entrance to exit: Ω = {(p, q): p entrance ≤ p ≤ p exit , q inner ≤ q ≤ q outer } .

Grid:
In the preceding section we partitioned Ω with a grid of curves of constant p or q.For convenience we used the same ∆p and ∆q throughout the grid, although this is by no means mandatory.Given two even integers N and M , let ∆p = p exit − p entrance N , ∆q = q outer − q inner M , p n = p entrance + n∆p, q m = q inner + m∆q.Elements: Each element is defined to be a quadrilateral consisting of four grid quadrilaterals with one common node in the middle of the element (see Figure 3).Thus, an element contains a total of 9 grid nodes.
An element E nm , where n and m are even, is bounded by the curves p = p n , p = p n+2 , q = q m , q = q m+2 .On each such element E nm we define local variables Local shape functions: On each element we define local shape functions where θ i are the Lagrange interpolant polynomials for the nodes t 1 = 0, For each of the 9 nodes in a given element, exactly one of the ϕ ij takes the value 1 at that node and vanishes on all curved grid edges in the element that do not contain the node.

Global shape functions:
For each grid node we construct a globally defined interpolating shape function by patching local shape functions on elements that share the node and extending by zero.The number of elements that can share a node is between one and four.Specifically, a global shape function corresponding to a particular node is defined by its restrictions to the elements.If an element does not contain the node the value is zero.
Otherwise it coincides with the local shape function on that element corresponding to the given node.
It is an easy exercise to see that the global shape functions are continuous.Note that while these shape functions are defined globally, the support of each global shape function is a subset of at most 4 adjacent elements.Moreover each global shape function is 1 at the node it corresponds to and vanishes on all grid edges not containing the node.In particular, it vanishes at all other nodes, thus earning the title interpolating.
Galerkin integrals: A general function in the approximating space G can be written in the form where ψ k are coefficients to be determined, g k are the global shape functions corresponding to each node in the grid after a suitable reindexing, and K = NM.Since g k are interpolating functions, ψ k are the values of ψ at the nodes.
To obtain the desired linear system we require the left hand side of (1) evaluated at ψ to be annihilated by the projection to G: We apply Green's first identity (10.713 [4]) to obtain Since the values ψ k of ψ at the nodes lying on ∂Ω are given by the boundary conditions, we need not project onto the corresponding g k .Let K be the set of all indices k such that g k corresponds to an interior node.The global shape functions g k with k ∈ K vanish on ∂Ω, so the contour integral on the right hand side of ( 4) is zero.Since Since g k are patched from ϕ ij defined on each element, we use additivity of integration and integrate in each element.Suppose that on an element E α = E nαmα , we have g k = θ i αk (t p )θ j αk (t q ).Then
We expand the derivatives of ψ where β is an index enumerating the 9 nodes of each element, to obtain where We use the two-dimensional spline formula (3) to write c p r, c q r, c p ∂r ∂q , c q ∂r ∂p , as linear combinations of products of Hermite interpolant polynomials H ki (τ p ) H j (τ q ); k, = 0, 1; i, j = 1, 2 in each of the four grid quadrilaterals in the element.The Lagrange interpolant polynomials θ ij and their derivatives are polynomials of t p or t q , so J αkβ is easily computed as a sum of four double integrals of polynomials in t p and t q of total degree ≤ 13 .

Verification
The results of the finite element procedure were compared to an analytic solution in a simple situation permitting the latter [5].In this case we used a grid given by exact streamlines and equipotential lines.The highest relative error in velocity was < 0.01%.Here is an excerpt from the results of comparison.In order to obtain accuracy information in practice, we have built internal consistency checks into the program.These include (i) the constancy with respect to q of the potential difference ϕ m from entrance to exit along a p-axis, (ii) the constancy with respect to p of the flow rate Q F along q-axes from the entrance to the exit.
We show an output excerpt for the results of internal consistency checks in a practical case.Internal consistency checks in a practical case.

Conclusions
• The results of the work show that consistent use of curvilinear coordinates in the finite element scheme leads to a significant increase in accuracy over conventional methods.
• Comparison with an exact solution does not reflect well the relative errors obtained in the solution of practical cases.This makes it absolutely essential to develop and perform consistency checks for each practical application.

Figure 1 :
Figure 1: Boundary correction for the finite difference scheme.