NASA Contractor Report 3313 

On One-Dimensional Stretching 

Functions for FiniteDifference Calculations 

Marcel Vinokur 

University of Samtu Clara 

Saizta Clara, Califoriiia 

Prepared for 

Ames Research Center 

under Grant NSG-2086 

National Aeronautics 

and Space Administration 

Scientific and Technical 

Information Branch 

1980 
INDEX 

i 

ii 

I. INTRODUCTION 

I I. TRUNCATION ERROR ANALYSIS FOR ONE-DIMENSIONAL STRETCHING 

FUNCTIONS 

111. A GENERAL TWO-SIDED STRETCHING FUNCTION 

IV. 

V. CONCLUDING REMARKS 

VI. REFERENCES 

A GENERAL INTERIOR STRETCHING FUNCTION 

APPENDIX A: EVALUATION OF ANTI-SYMMETRIC TWO-SIDED STRETCHING 

FUNCTIONS 

INVERSION OF y = sinh x/x 

INVERSION OF y = sin x/x 

FUNCTION DERIVED FROM w = tanz 

APPENDIX B: 

APPENDIX C: 

APPENDIX D: SOLUTION CHARACTERISTICS OF THE TWO-SIDED STRETCHING 

APPENDIX E: EVALUATION OF ANTI-SYMMETRIC INTERIOR STRETCHING 

FUNCTIONS 

1 

4 

9 

18 

22 

27 

28 

34 

38 

40 

48 
Figure 1. Two-sided Stretching Function for So = 100 26 

Figure D. 1. Solution Characteristics of Two-sided Stretching 46 

and Different Values of S1 

Functions Derived from w = tanz 

Figure D.2. Locus of Solution in z-Plane 47 

i 
ABSTRACT 

The class of one-dimens 

difference calculations 

onal stretching functions used in fin 

is studied. For solutions containing 

tea highly 

localized region of rapid variation, simple criteria for a stretching 

function are derived using a truncation error analysis. These criteria 

are used to investigate two types of stretching functions. One is an 

interior stretching function, for which the location and slope of an 

interior clustering region are specified. The simplest such function 

satisfying the criteria is found to be one based on the inverse 

hyperbolic sine. 

type of function is a two-sided stretching function, for which the 

arbitrary slopes at the two ends of the one-dimensional interval are 

specified. The simplestsuch general function is found to be one 

based on the inverse tangent. 

both equal and greater than one was first employed by Roberts. 

general two-sided function has many applications in the construction Of 

finite-difference grids. 

of the references. 

It was first employed by Thomas et al. The other 

The special case where the slopes were 

The 

An example of such an application is found in one 

ii 
I. INTRODUCTION 

Finite-difference calculations of fluid flow problems are best carried 

out using an equispaced grid in a rectangular (or cubic) computational 

domain, with the flow variables and components of the position vector 

as dependent variables, and boundary conditions applied at the edges 

(or faces) of the domain. In order to minimize the number of grid 

points required for a given accuracy, one seeks boundary-fitted coordinate 

transformations that cluster points in regions where the dependent 

variables undergo rapid variation. 

geometry (very large curvatures or corners), compressibility (entropy 

layers, shock waves and contact discontinuities), and viscosity (boundary 

layers and shear layers). 

of such regions of various length scales, and often of unknown location. 

An ideal grid would adjust with each time or iteration step to maintain 

optimum clustering. Such adaptive grid methods, which involve the 

solution of auxiliary equations, have been developed for one-dimensional 

problems (refs. 1-3). Their extension to complex multi-dimensional 

flows is a difficult problem, particularly when the regions requiring 

clustering do not have simple topological properties required by a 

fi ni te-difference grid. 

These regions may be due to body 

A complex flow may thus contain a variety 

There are many practical problems in which the locations and length 

scales of regions of rapid variation can be estimated a priori 

(e.g., known geometry, attached boundary layers, simple shock wave 

configurations). In these cases the clustering can be incorporated 

in automatic grid generators which solve an elliptic boundary-value 

problem (refs. 4-6). The distribution of grid points on the boundaries 

is then normally prescribed algebraically, using one-dimensional 

stretching functions. (Here, stretching function refers to any 

transformation involving stretching or clustering). 

to employ stretching functions to obtain a clustered grid from an 

unclustered grid by applying clustering to one coordinate family only. 

For some simple geometries, one can construct entire clustered grids 

purely algebraically, using only one-dimensional stretching functions. 

It is also possible 

1 
The simplest class of one-dimensional stretching functions is that 

involving two parameters. In interior stretching functions, the 

parameters are the location and slope of a single clustering region. 

In two-sided stretching functions, the slopes at the two ends of the 

one -d imens i onal i n terval a re s peci f i ed . .The anti - symmetri c two - si ded 

stretching function (with the same slope at each end) is of special 

interest, since the portion from the midpoint to either end defines 

a one-sided stretching function (with zero curvature at the other 

end). 

case of the interior stretching function, by locating the clustering 

region at one end. 

case, these two one-sided stretching functions are of a different nature. 

A one-sided stretching function can also be obtained as a special 

Since the clustered end has zero curvature in this 

An interior stretching function based on the inverse hyperbolic sine was 

employed by Thomas, Vinokur, Bastianon, and Conti (ref. 7) in a 

numerical solution of inviscid supersonic flow over a blunt delta wing 

with elliptical cross section. The function was used to cluster points 

on the body at the vertices of very eccentric ellipses. The one-sided 

version clustered points in the flow near the body surface to resolve 

the entropy layer due to the bow shock. No derivation of the stretching 

function was given, and the clustering parameter appearing in it was not 

related to the length scales of regions of rapid variation in physical 

space. 

An anti-symmetric two-sided stretching function of a logarithmic type 

was employed by Roberts (ref. 8) to study boundary layer flows. 

heuristic derivation of the function avoided consideration of the 

truncation errors associated with finite-difference approximations. 

While this function has been used successfully in many flow calculations, 

there is a need for a general two-sided function which allows arbitrary 

stretching or clustering to be specified independently at each end. An 

application would be problems in which the appropriate length scales 

requiring clustering are significantly different at the two ends. 

application is the distribution of grid points on a curve which is 

The 

Another 

2 
defined piecewise, where continuity of grid spacing is desired at the 

ends ofthe piecewise segments. 

a function as a blending or interpolating function to construct two 

and three-dimensional grids using one-dimensional stretching functions 

and shearing transformation. 

ordered grid for wing body flows by Yinokur (ref. 9) is an example of 

these applications. 

A third application is the use of such 

The construction of a single, wellThe present work has two objectives. 

criteria for one-dimensional stretching functions, by considering the 

truncation errors inherent in finite-difference approximations. The 

functions introduced by Thomas et a1 and Roberts will be found to be 

the simplest ones satisfying these criteria. 

derive a simple form of the general two-sided stretching function. 

One is to obtain simple, rational 

The other objective is to 

3 
11. TRUNCATION ERROR ANALYSIS FOR ONE-DIMENSIONAL STRETCHING FUNCTIONS 

An exact analysis of the truncation error inherent in a finite-difference 

calculation would require knowledge of the equation being solved and the 

finite difference approximation that is used. Here we are concerned with 

the special situation where the solution contains a highly localized 

region of rapid variation with respect to some coordinate, and we seek 

approximate criteria for a stretching function that will be independent 

of the equation or difference algorithm. 

approximated are in general non-linear functions of the unknowns and 

their spatial derivatives. The error analysis will be performed in 

terms of the fractional truncation errors for the spatial derivatives. 

The quantities that are 

Let F(t) be the equation describing a 5-coordinate curve, where t is 

any parameter that varies smoothly with arc length. 

the ends of the curve, we introduce the normalized variables 

If A and B denote 

ranging from 0 to 1. For simplicity, all partial derivatives with 

respect to 5 or t will be written as total derivatives. Let @(t) be 

any function of the unknowns. Outside of the region where 

we can define a natural length scale of the variation of $I 

to t as 

d$I/dt = 0, 

with respect 

Since the components of 

calculation of metrics and Jacobians, we similarly define the natural 

length scale of the variation of 

also enter as dependent variables in the 

with respect to t as 

4 
Note that if t is proportional to arc length, then Lrt is precisely 

the radius of curvature, normalized by the length of the curve, 

Assume first that t is used as the normalized computational variable 

(i.e. 5 = t), with At as the uniform grid spacing. 

finite difference approximation to d$/dt. 

write any first order accurate approximation as 

Let 6@/6t denote the 

Using equation (2), one can 

6 %= +$+o 

Similar y, the approximat on to dF/dt can be written as 

- 

-+ sF= x[l dr t 0 At)] . 6t 

If L-’Ot or L- rt became very large in some localized region, then a 

prohibitively small At would be required to obtain a desired fractional 

truncation error. 

of grid points would be wasted. 

and L-’ $E rS computational variable E; for which Lcould be local y very large. though L-l +t or L-lrt 

Outside of the localized region, the excessive number 

The obvious remedy is to seek a new 

remain of 0(1), even 

With the aid of the identities 

and 

,-. 

and definition (2), one can easily show that 
Similarly, using definition (3) , one obtains the inequality 

One criterion for the stretching function C(t) is therefore 

In addition, we require that 

and 

dt 1- L- rt dg = O(1). 

Consider the case where L-l 

follows from equation (11) that 

is very large in some localized region. It 

9t 

in that region. Since L-l 

thickness is of O(L 

interval. Furthermore, since L-l 

region, it follows from equation (11) that dt/dC = 0(1), i.e. that dt/dC 

does not become large anywhere. 

are satisfied if condition (10) is valid. 

remains large over an interval whose v- 

), we require that dt/d< remains small over that 

@t = 0(1) outside of the localized @t 

But these two additional requirements 

Noting that 

it follows upon integration over a finite interval At that 

6 
Applying equation (15) to the localized region, for which At = O(L 

we find,using equation (lo), that 

), $t 

Thus equation (13) is satisfied over the entire localized region. 

Gt = 1 in equation (15), we see that dt/dE = 0(1) is satisfied over the 

complete range of t. In the event tnat L-lrt is larger than L-l +t in the 

localized region, then condition (13) is replaced by 

Letting 

The above analysis is easily extended to higher order finite-difference 

approximations, as well as the treatment of higher derivatives. 

to consider fractional truncation errors due to a second-order accurate 

approximation to d$/dt (and also for regions where d $/dt2 = 0), it is 

appropriate to define a different length scale of the variation of$with 

respect to t as 

In order 

2 

We now require that L-’ remains of 0(1), even though r-’ 

very large. Using the identity 

could be locally $E $t 

we obtain the inequality 

;-2 5 - 1 dt)2 dt -2 

+E - (L- @t dg 

In addition to satisfy 

must sat i sfy 

- = O( tc 

ng equation (lo), the stretching function t;(t) 

1. (21 1 

7 
in the localized region, condition (13) @t Also, if L-' 

must be replaced by 

is larger that L-l 

@t 

A first order f 

i n troduc t i on of 

22 nite difference approximation to d @ / d t r e q u i r e s the 

the length scale 

&$/$I 

dt 

Combining equations (2), (18), and (23), we see that 

Thus conditions (lo), (21), and (13) or (22) are sufficient to guarantee 

that Z-l remains of O(1). The length scale is also the appropriate 

one to use at a point where d@/dt = 0. At such a point, using equations 

46 @t 

(7) and (lsl), we obtain the 

The criteria for <(t) is aga 

by 

'el at i on 

(25) 

t5. + 3 L-l 

n equation (lo), with equation (13) replaced 

if the localized region of rapid variation occurs around the point d@/dt = 0. 

Conditions (22) and (26) are to be replaced by analogous ones containing 

Lrt 

- - 

and Lrt, if these are the more significant length scales. 

In summary, one first defines length scales appropriate to the difference 

approximation and the location of the region of rapid variation. 

criteria for the stretching function c(t) can be stated as: 

The 

8 
1. All the inverse length scales of the variation of t with respect to 5 

must be at most of order one throughout the range of t. 

2, The slope dt/dC must be of the order of the minimum length scale of 

the variation of @or ? with respect to t in the localized region of 

rapid variation. 

These criteria will insure that most of the grid points will be concentrated 

in the localized region of rapid variation, with a sufficient number of 

points left in the remainder of the domain. 

investigate the two-parameter stretching functionsof the next two sections. 

The criteria will be used to 

111. A GENERAL TWO-SIDED STRETCHING FUNCTION 

In this section we derive a general two-sided stretching function 

c(t; so, sl), where€,and t are normalized variables defined by equation (l), 

and the parameters so and s1 are dimensionless slopes defined as 

and 

In order to be useful for contructing finite-difference grids, the 

function must be monotonic, and satisfy conditions (10) and (21) even if 

so or s 

would be desirable for the function to be simple, invertible, and to vary 

continuously over the complete ranges of so and sl. 

becomes very large. For a general range of applications, it 1 

An attractive candidate for such a stretching function is a scaled portion 

of a single, universal function w(z). 

function will be obtained by properly scaling the portion of the universal 

function from corresponding points zo and zl. An additional requirement 

is that the unnormalized function i(i) be independent of the designation 

For a given so and s, , the stretching 

9 
of a particular end as A or B. 

tne universal function to be odd, i.e. 

One can easily show that this restricts 

w(-z) = -w(z). 

The simplest odd, monotonic, invertible functions are sinz and tanz. 

Their hyperbolic relatives are produced by letting z be complex. 

inverse functions are formed by associating z with either t or 5. One 

can determine whether either universal function is suitable as a basis 

for a stretching function by applying conditions (lo) ;.nd (21) for 

very large so or s,. 

simpler anti-symmetric casesO = sl, which corresponded to zo =-zl, with 

z being either real or pure imaginary. 

The 

Actually, the most extreme test occurs for the 

An evaluation of the anti-symmetric two-sided stretching functions 

obtained from sinz and tanz is carried out for the case so = s >1 

in Appendix A. Only tanz produces a stretching function satisfying 

conditions (10) and (21), with the inverse length scales being 

logarithmically of O(1). The stretching function t(t) is a scaled 

portion of the inverse hyperbolic tangent. 

tangent as a logarithm, we obtain exactly the function derived by 

Roberts (ref. 8). 

of t, a property that was used by Roberts to define his function. 

suggests a related stretching function, for which L-l 

linear function of 6. The corresponding universal function is the 

error function erfz. The associated stretching function is also 

analyzed in Appendix A, and found to satisfy conditions (IO) and (21). 

However, it is not invertible, and has larger maximum inverse length 

scales than the former stretching function. 

1 

Expressing the hyperbolic 

It turns out that L-l is a piecewise linear fminction tS 

This 

is a piecewise tS 

On the basis of the above considerations, the universal function w = tanz 

will be used to obtain a stretching function for arbitrary so and sl. 

Introducing the ranges 

10 
and 

Aw = tanzl - tanzO 

we can define the normalized variables 

5 = (z - zo)/Az 

and 

t = (tanz - tanzo)/Aw. 

The slope of the stretching function is then given by 

Using the trigonometric identity 

sin (z, - z,) I v 

0, cos z1 cos z tanz, - tanzO = 

we find for the parameters so and s1 the relations 

sinAz cos zo 

Az cos z1 - 

and 

sin Az cos z, 

= Az cos zo 

This suggests introducing the new parameters 

B =/= 

and 

A = ,/x 

(32) 

(33) 

(34) 

(35) 

11 
The parameters A and B can then be expressed in terms of z 

as 

and z1 0 

sin Az 

Az 

B= ___ 

(38) 

0 

1. 

and cos z 

A= ____ 

cos z 

Using the cosine sum dentity, we can also write equat 

A = cos Az t tanz sin Az 1 

and 

1/A = cos Az - tanzO sin Az. 

on (39) as 

(40a) 

For a given value of B, 

(40) can then be solved to obtain zo and Aw for a given value of A. 

stretching function obtained from equations (30) and (31) can then be 

written as 

Az is obtained by solving equation (38). Equation 

The 

tan ([Az + zo) - tanzO 

t= Aw (41 1 

While equation (41) is a formal expression for the general stretching 

function, it cannot be used for calculations in its present form. 

on the value of B, Az and Aw are either real or pure imaginary. 

ranges of A and B, zo can become complex. 

equation (39), we can eliminate zo from equation (41) and obtain instead 

Depending 

For certain 

Ufing the tangent sum identity and 

tan gAz 

A sin Az + (1 - Acos Az) tan gAz . t= (42) 

This can be further simplified by noting that A = 1 corresponds to the 

anti-symmetric solution which was analyzed in Appendix A. 

this solution as u(<). 

tangent sum identity, we obtain 

Let us denote 

Setting A = 1 in equation (42), and using the 

12 
(43) 

In terms of u, equation (42) then takes the simple form 

U 

it = A t (1 - AJu, 

which can be readily inverted as 

t 

(1/A) + (1 - l/A)t. u= 

(44) 

(45) 

Note that both u(t) and its inverse can be obtained as scaled portions 

of a rectangular hyperbola. 

ends are A and l/A, respectively. 

calculational purposes, equations (44) and (45) are well behaved in the 

neighborhood of A = 1. 

For each function, the slopes at the two 

Finally, we observe that for 

We thus have the remarkable result that one essentially needs to know 

only the anti-symmetric stretching function for the geometric mean of 

the slopes So and S1. The square root of the ratio of those slopes 

determines an additional simple transformation which produces the 

desired stretching function. 

invertible, the resultant stretching function is also invertible. The 

key trigonometric property making this result possible is that the 

tangent of a sum is a rational function of the individual tangents. By 

contrast, the sine of a sum involves the individual sines and cosines, 

and is not expressible as a rational function of sines alone. 

of a stretching function based on w = sinz has been carried out, but is 

not presented here. 

be the arithmetic mean and difference of the two slopes. 

into two functions corresponding to equations (43) and (44) is not 

possible, and one must use the direct form corresponding to equations 

(41) and (42). 

quadratic equation, and the sign of A must be tested in order to choose the 

appropriate root. 

Since both equations (43) and (44) are 

An analysis 

The parameters corresponding to B and A turn out to 

A separation 

For B > 1, the inversion involves the solution of a 

It is indeed serendipity that the tangent function 

13 
dictated by truncation error considerations is also the much simpler one 

for constructing a general two-sided stretching function. 

The calculation of the anti-symmetric function depends on the size of B. 

If B > 1, it follows from equation (38) that z is imaginary and we obtain 

the results (previously obtained in Appendix A) 

and 

sinh Ay 

AY 

B= 

tanh CAY( 5 - 1/2)] 

2 tanh (Ay/2) u = 1/2 -t 

The inversion of equation (47) yields 

tanh-’ [(Zu - 1)tanh (Ay/2)1 

2AY 5 = 1/2 + 

(47) 

Note that the hyperbolic tangent and its inverse can be expressed in 

terms of exponentials and a logarithm, respectively. 

For B < 1, Az is real, and the corresponding results are 

tan-’ [ (2u-1 )tan (Ax/2)1 

2AX and 5 = 1/2 -t 

When B is very near one, both of the above formulations break down, 

since ax and ~y approach zero. 

by expanding equations (49) and (50) in powers of Ax. 

in B-1, one obtains 

The appropriate expressions are obtained 

To first order 

u 6 [l -t 2(B - 1) (c-.5)(1 -<)I 

and 

5 u [1 - 2(B - 1) (U -.5)(1 - u)]. 

(52) 

(53) 

14 
By scaling half of the above functions, one obtains the one-sided 

stretching functions with s given at t = 0 and zero curvature at t = 1. 0 

The results are: 

so > 1 sinh 2Ay 

so = 2Ay 9 

and 

so< 1 

tanh-’ [(t - 1)tanh Ay] 

AY 5=1+ 

= sin 2Ax 

~AX 9 

tan [AX(< - 1.11 

tan-’ [(t - 1 )tan ~xl 

= ’ tanAx ¶ 

AX <=1+ 

(54) 

(55) 

(56) 

(57) 

(58) 

(59) 

The two-sided stretching functions for B>1 and one sided stretching 

function for so>l require the inversion of the function 

y = sinh x/x. (62) 

An approximate analytic representation of the inverse function x = f(y), 

valid over the range of y required by a stretching function, is derived 

in Appendix B. The results are as follows: 

For y < 2.7829681 

(63) 

x = (1 - .15i t .057321429y2 - .024907295i3 

t .0077424461y4 - .0010794123y5), 

15 
where 

y=y-1. 

For y> 2.7829681 

x = v t (1 t l/v) log (2~) - .02041793 + .24902722~ 

+ 1 .9496443w2 - 2.6294547~ 3+ 8.56795911w4, (65) 

where v = log y (66) 

W = (67) l/y - .028527431. and 

The maximum magnitude of the fractional error in y as defined by equation 

(8.2) is .000267732 for 1 < y < 69.64. The magnitude of the error reaches 

.00083 at y = 120.5. These errors are small enough so that the grids 

constructed by the resultant stretching functions wi 11 exhibit slope 

discontinuities that are negligible within the accuracy of any practical 

finite-difference approximation. 

The two-sided stretching function for B < 1 and one-sided stretching function 

for so< 1 require the inversion of the function 

y = sin x/x. (68) 

An approximate analytic representation of the inverse function is derived 

in Appendix C. The results are as follows: 

For y < .26938972 

23 4 

x = ~[l - y t y2 - (1 tn/6)y t 6.794732~ 

- 13.205501y5 + 11 .726095y6]. 

For .26938972 < y < 1 

x = h$ (1 + .15j + .057321429Y2 + .048774238y3 

- . 053337753j4 + .0758451 34y5), 

where 

j=l -y. 

16 
The maximum magnitude of the fractional error in y defined by equation 

(C.2) is .00019717, which is small enough for numerical applications. 

While it is appealing to view the general two-sided stretching function 

as a distortion of an anti-symmetric stretching function via equations 

(43) and (44), there are advantages in looking at the more basic forms 

of the solution given by equations (41) and (42). 

of operations count is important, then the optimum form for the case 

B > 1 is derived from equation (42), as either 

If efficiency in terms 

t= 

or t = 

tanh EA 

sinh A; t (1 - A cosh Ay 

2EAY -1 e 

e' (1 - Ae -Ay) t AeAy -1. 

The most efficient form for the case B< 1 is equation (4'1) with z replaced 

by x throughout. 

It is also instructive to study the general solution (41), and see how it 

reduces to different real representations for various ranges of A and B. 

This is carried out in Appendix i). 

B < 1, the solution is a scaled portion of a tangent, while for B = 1 

it is a scaled portion of a rectangular hyperbola. 

corresponding Ay given by equation (46) , the representation depends on 

the value of A. 

the solution is a scaled portion of a hyperbolic tangent. 

that range the corresponding function is the hyperbolic cotangent, 

while on the boundary it is the exponential function. 

It is easy to demonstrate that for 

For B>1 (and the 

As shown in Appendix il, if exp (-Ay) < A< exp (Ay), 

Outside of 

When A = 1, the stretching function is anti-symmetric, and the solution 

curve contains an inflection point. 

point moves towards one end, and eventually could disappear. 

in Appendix D that an inflection point will be present if l/cosh Ay < A 

< cosh Ay for B>1, and If B<2/n, the 

solution must always contain an inflection point. Figure D.l is a log-log 

plot of the B vs A plane, showing the inflection point regions, as well 

As A departs from one, the inflection 

It is shown 

cos Ax<A<l/cosAx for B<1. 

17 
as the boundaries for the various real representations of the solution. 

Another way to study equation (41) is to follow the path of the solution 

in the complex z-plane. 

results plotted in Figure D.2. The plot shows clearly the transition from 

one representation to another as the parameters A and 6 are varied. 

This is also carried out in Appendix D, with the 

A typical set of curves for the two-sided stretching function t (t) is 

shown in figure 1 for so = 100, and for values of s1 ranging from.1 to 

100. The figure shows the smooth transition from solutions containing 

an inflection point, for large sly to those without an inflection point 

for small sl. 

for s1 = 100 still leaves a sufficient number of points to resolve the 

central region. 

Note that the high concentration of points at the two ends 

IV. A GENERAL INTERIOR STRETCHING FUNCTION 

In this section we derive a general interior stretching function 

<(t; si, ti), where si is the dimensionless slope at the inflection 

point ti, i.e. 

We limit our consideration to si >1, which is the only case of practical 

interest. 

odd universal function w(z). As in Section 111, the simple functions 

sinz and tanz are considered first, and conditions (10) and (21) are 

examined for the anti-symmetric case (ti = 1/2) when si is very large. 

The evaluation is carried out in Appendix E, and sinz is found to produce 

an appropriate function with inverse length scales that are logar thmical 

of O(1). The stretching function t(t) is a scaled portion of the inverse 

hyperbolic sine. 

We again look for a function which is a scaled portion of an 

Y 

18 
1 

The general interior stretching function for arbitrary ti is readily 

obtained from the universal function w = sinz, by letting z = iy. In 

terms of the range Ay = y - yo, and the implicity defined Ei=S(ti), 

the final result can be written as 

I 

1 

I The inverse function is 

1 -1 5 =$ + sinh [(t/ti - 1) sinh CiAy]. 

The relation between$and ti (for a given Ay) is obtained from equation 

(75) by setting t = 5 = 1. 

I 

The result can be written as 

-- - 1 - cosh Ay + sinh Ay coth Sib~. 

ti 

(77) 

Expressing the inverse hyperbolic cotangent in terms of a logarithm, we 

can write the inverse of equation (77) as 

Equations (76) and (78) are precisely the ones given by Thomas et al. 

~ (ref. 7) 

If si and ti are the given parameters, the corresponding value of Ay 

must be calculated. 

and substituting into equation (73). 

This can be done by differentiating equation (75) 

Using equation (77) to eliminate 

Si, we can write the result as 

Ay - 1 + l/ti 

1 2- sinh Ay l2 -1 

(79) 

This is an implicit equation for Ay involving two independent parameters. 

If the interior point is not too close to either end, and the slope si 

is sufficiently large, one can obtain a simplification. 

exp (-2Ay)c 51, we can approximate equation (79) as 

Assuming that 

19 
- 

sinh (Ay/2) 

2Si/ \ ti (1 - ti)' AY/ 2 

Equation (80) i s i n the form of equation (62), and one can use the 

approximate analytic inversion derived in Appendix B. 

The special case of an anti-symmetric solution is obtained by setting 

ti =Si = 1/2. The results (derived in Appendix E) are 

and 5 = 1/2 t 1 sinh-' [(2t-l)sinh(Ay/2)], (82 1 AY 

- sinh (Ay/2) . where si -'1/2 - AY/2 

A one-sided stretching function is obtained by setting ti = ti = 0. 

The results can be written as 

sinh (toy) 

sinh Ay t= 

and <=- 1sinh-'(t sinh Ay) , 

AY 

- sinh Ay 

AY - where Si = So - 

(83) 

It is interesting to compare the one-sided stretching function derived 

from the hyperbolic tangent (equations (54) to (56)) with the above 

function derived from the hyperbolic sine. 

stand for the tangent and sine solutions, we see from equations (54) and 

(86) that for a given so, 

Letting the subscripts T arid S 

2Ay, = Us. 

20 
The maximum inverse length scale L-' 

occurs at t = 0 for the tangent solution, and has the value 

as defined by equation (2) tS 

) ((L- tc max T ) = 2AyT tanh Ay, . 

For the sine solution, the maximum occurs at t = 1, and has the value 

1 ((L- tS)max)S = Ays tanh Ay, . 

Thus, for so large enough so that tanh AyT N tanh Ay, = 1, the two 

solutions have the same maximum inverse length scale. 

(dc/dt) min.occurs at t = 1 for both solutions. 

The minimum slope 

The results are 

tanh Ay, 

( (dS/dt)min)T = 

AYT 

tanh Ay, 

and ( (d</dt)rni n), = 

For large so, we thus obtain 

hyperbolic sine, for identical clustering at t = 0. 

due to the fact that the zero inflection point occurs 

first, and t = 0 for the second. The particular appl 

determine which of these two is preferable. 

The one-sided function derived from the hyperbolic tangent thus has more 

points at the unclustered end ( t = 1) than the one derived from the 

he difference is 

at t = 1 for the 

cation woul d 

21 
V. CONCLUDING REMARKS 

In this work it has been assumed that the metrics and Jacobians that 

arise in the transformed equations are calculated by finite differences. 

If the equation ;(t) of the 5 -coordinate curve is known analytically, 

and the transformation t(€,) is also given analytically, then the nietrics 

and Jacobians can be analytically determined from the derivatives 

d;/dt and dt/dt. The truncation error in the numerical calculation will 

then be due solely to the finite-difference approximations to the 

derivatives of @ with respect to 5 . 

the optimum transformation would be one in which 5 varied linearly with 

@, since this would result in zero truncation errors. 

When @ varies monotonically with t, 

In order to compare transformations for the numerical and analytic 

treatment of Jacobians and metrics, consider a strictly one-dimensional 

case in which the single unknown u varies monotonically with distance x. 

Assume a highly localized interior region of rapid variation whose 

thickness is proportioned to the small parameterv. A simple exatiiple 

of such a solution is 

u - tanh (x/v), (93) 

where x = 0 in the localized region. 

solution of Burgers' equation with fixed end conditions. For 

the analytic treatment of Jacobians and metrics, it follows that c(x) 

should be a scaled protion of the hyperbolic tangent. The analysis of 

Appendix E, which assumes a numerical treatment of Jacobians and metrics, 

shows that this choice is completely unsuitable, and instead favors a 

scaled portion of the inverse hyperbolic sine. 

equation is written so that only derivatives of u appear, then the 

hyperbolic tangent transformation should lead to a numerical solution 

with no truncation errors. But the equation can also be written in a 

form involving derivatives of several functions of u. 

a strong conservation form. 

This is actually the steady-state 

If the differential 

An example is 

In this situation one has several variables 

22 
$(<) to be approximated by finite differences, and no single transformation 

t (5) is optimum for all of them. 

transformation would put all the interior grid points inside the region 

of rapid variation. There would be no points to resolve the boundary of 

this region. By contrast, the inverse hyperbolic sine transformation 

puts a sufficient number of points outside the region of rapid variation 

to resolve the complete one-dimensional region. 

If vis very small, the hyperbolic tangent 

The above discussion indicates that there are special situations and forms 

of the differential equations for which an analytic treatment of Jacobians 

and metrics can provide a desired accuracy with fewer grid points than 

the numerical treatment. These cases appear to be restricted to monotonic 

distributions that can be approximated by simple analytic expressions. 

general applications of one-dimensional stretching functions , these 

special situations will not be met. It is then best to treat the Jacobians 

and metrics numerically, and use the stretching functions derived in 

Sections I11 and IV. 

For 

Another assumption in the derivation of the stretching functions is that 

the dimensionless length scale of the localized region of rapid variation 

could be extremely small. 

transformation dudt to be extremely large. 

encountered, and transformation slopes remain of 0(1), then the form of the 

stretching function is not critical. For example, many authors have used 

a scaled exponential as a one-sided stretching function. This is perfectly 

reasonable as long as the one-sided slope +, is not much larger than one. 

But one can readily show that the maximum inverse length scale is exp(so) 

for large so. 

one-sided stretching function for very large slopes. 

This requires the dimensionless slope of the 

If this condition is not 

Thus a simple exponential does not yield a suitable 

It should be made clear that Roberts was not the first one to use a 

stretching function involving the hyperbolic tangent. 

(ref. lo), in an investigation of flow in a two-dimensional channel with 

a rectangular cavity, used a transformation hased on the hyperbolic 

tangent to transform a semi-infinite region into a finite computational 

region, and to cluster grid points at the corners of the cavity. 

Mehta and Lavan 

The 

23 
maximum dimensionless slope, based on the lengtn of the cavity, had a 

value of .8664 in their calculations. The transformation would have been 

a poor one if they had required a very large slope at the corner, as 

shown in Appendix E. 

based on the inverse hyperbolic tangent in a study of the two-dimensional 

flow around an airfoil. 

region external to the airfoil into a unit circle. The stretching function 

was necessary to cluster points further near the airfoil surface (to 

capture the boundary layer) as well as the free stream (to overcome the 

stretched grid produced by the first transformation). 

the non-dimensional slopes at the two ends were 5.77 and 27.8. In this 

instance, the use of the hyperbolic tangent was both appropriate and 

necessary, as shown by the analysis of Section 111. 

The sanie authors (ref. 11) used a stretching function 

A previous transformation had transformed the 

In their calculation, 

An important criterion in the development of a two-sided stretching 

function is a continuous behavior as So and S, varied from zero to infinity. 

This is necessary to obtain smooth grids constructed algebraically using 

one-dimensional stretching functions. For B > 1 , the required function was 

found to be based on the inverse hyperbolic tangent. 

the same function could be used for B< 1, simply by interchanging cand t 

in the expression. 

of thiswork). 

At first glance, 

(This is what was actually done in the earlier stages 

But this would violate the desired continuous behavior 

in the neighborhood of B = 1. 

hyperbolic tangent is the inverse tangent, which differs from the hyperbolic 

tangent. 

self invertible in this sense,but they do not include the elementary 

functions, and therefore would not be useful as stretching functions. 

Since the inverse tangent is not self invertible, it is necessary to use 

two different representations in calculating the stretching function 

nlrmeri cal ly. 

The analytic continuation of the inverse 

One can actually construct anti-symmetric function which are 

The use of the stretching functions derived in this work requires 

specifying their slopes at one or two points. 

obtained from matching slopes with another function, or estimating the 

length scale of a localized region where an appropriate dependent variable 

The values are either 

24 
undergoes rapid variation. 

derivation of an appropriate length scale is not easy, and the value to 

assign to the slope can be somewhat arbitrary. 

consistent criterion is used in assigning slope values, useful grids for 

numerical calculations can be generated. 

grid which was generated using the general two-sided stretching function 

is found in reference 12. 

Admittedly, in a complex situation, the 

Nevertheless, if a 

A recent example of a complex 

25 
XI 

0 

FIGURE I TWO- SIDED STRETCHING FUNCTION FOR 5i,= 100 

AND DIFFERENT VALUES OF SI. 

26 
VI. REFERENCES 

1. 

2. 

3. 

4. 

5. 

6. 

7. 

8. 

9. 

10. 

11. 

12. 

Gough, D.O., Spiegel, E.A., and Toomre, J., "Highly Stretched Meshes 

as Functionals of Solutions" , Proceedings of the Fourth International 

Conference on Numeri cal Methods i n F1 ui d Dynamics , Spri nger-Verl ag , pp. 191-196, 1975. 

Ablow, C.M. , and Schechter, S., "Campylotropic Coordinates", J. Comp. 

Phys., Vol. 27(3), pp. 351-362, June 1978. 

Pierson, B.L. , and Kutler, P., "Optimal Nodal Point Distribution 

for Improved Accuracy in Computational Fluid Dynamics", AIAA J. , Vol. 18(1), pp. 49-54, January 1980. 

Thompson, J.F. , Thames, F.C. , and Mastin, C.W. , "TOMCAT- A Code for 

Numerical Generation of Boundary-Fi tted Curvi 1 inear Coordinate 

Systems on Fields Containing Any Number of Arbitrary Two-Dimensional 

Bodies", J. Comp. Phys., Vol. 24(3),pp. 274-302, July 1977. 

Middlecoff, J.F., and Thomas, P.D., "Direct Control of the Grid Point 

Distribution in Meshes Generated by Elliptic Equations", Proceedings 

of the AIAA Fourth Computational Fluid Dynamics Conference , Williamsburg, Virginia, pp. 175-179, July 23-25, 1979. 

Steger, J.L. , and Sorenson, R.L. , "Automatic Mesh-Point Clustering 

Near a Boundary in Grid Generation With Elliptic Partial Differential 

Equations", J. Comp. Phys. , Vol. 33(3), pp.405-410, December 1979. 

Thomas, P.D., Vinokur, M. , Bastianon, R.A. , and Conti , R.J. , "Numerical 

Solution for Three-Dimensional Inviscid Supersonic Flow", AIAA J., 

Roberts , G.O. , "Computational Meshes for Boundary Layer Problems" , Proceedings of the Second International Conference on Numerical 

Methods in Fluid Dynamics", Springer-Verlag, pp. 171-177, 1971. 

Vinokur, M., Steger, J.L., and Pulliam, T.H., "On Use of Warped 

Spherical Coordinates to Generate We1 1 Ordered Fini te-Di fference 

Grids for Wing-Body Flows. Open Forum, AIAA 11th Fluid and Plasma 

Dynamics Conference, Seattle, Washington, July 10-12, 1978. 

Mehta, U.B., and Lavan, Z., "Flow in a Two-Dimensional Channel with 

a Rectangular Cavity, J. Appl. Mech., Vol. 36, Series E, No. 4, 

pp. 897-901 , December 1969. 

Mehta, U.B., and Lavan, Z. , "Starting Vortex, Separation Bubbles 

and Stall: A Numerical Study of Laminar Unsteady Flow Around an 

Airfoil". J. Fluid Mech., Vol. 67(2), pp. 227-256, 1975. 

Vol. 10(7), pp. 887-894, July 1972. 

Lombard, C.K. , Davy, W.C. , and Green, M.J. , "Forebody and Base 

Region Real-Gas Flow in Severe Planetary Entry by a Factored 

Implicit Numerical Method- Part I (Computational Fluid Dynamics). 

AIAA Paper No. 80-0065, AIAA 18th Aerospace Sciences Meeting, 

Pasadena, California, January 14-16, 1980. 

27 
APPENDIX A 

EVALUATION OF ANTI-SYMMETRIC TWO-SIDED STRETCHING FUNCTIONS 

In this appendix we evaluate several candidates for an anti-symmetric 

two-sided stretching function <( t) , where 5 and t are normalized variables 

ranging from zero to one. 

designated as 

The common slopes at the ends will be 

The truncation error analysis of Section I1 results in the conditions 

even if B is very large. The corresponding portion of a universal odd 

function w(z) ranges from zo = -A.,-/2 to z1 = Az/2, where Az is the total 

range, and z is either real or pure imaginary. The stretching function 

is therefore defined as either 

or 

L 

2 W (AZ/2) 

Since conditions (A.2) and (A.3) are difficult to satisfy only for very 

large By we restrict the analysis to the case B>l. 

28 
w = sinz 

If we let z be real, so that Az = Ax, the appropriate function B>1 

is obtained from equation (A.5) as 

I t-7 1 - sin [AX (E- T)] 

2 sin AX/^) . 

Using equation (A.l)y we obtain the relation 

tan (A x/2) 

A X/2 B= 

takes the form tS The inverse length scale L-l 

1 L-l t< = AxItan[Ax (5 -91 Iy 

and has a maximum value 

For large B, Ax -f T, and we obtain 

(A.lO) 

Thus equation (A.2) is violated, and the function (A.6) is not suitable. 

If we let z be imaginary, so that Az = iAy, the appropriate function for 

B > 1 is obtained from equation (A.4) as 

- 1 = 2 sinh[Ay(t - i)] 

sinh (Ay/2) 

Applying equation (A.l), we obtain the relation 

(A.ll) 

29 

(A.12) 
The inverse length scale L-l is now 

1 tS = 2 sinh(AY/Z) I sinh CAY (t - 7)II 

L-l tc 1 1 + sinh [ Ay(t - 711 

is readily found to be tS The maximum value of L-l 

= sinh (AY/2). (L- tc max 

For large B, Ay/Z-+B, and we obtain 

(A.13) 

(A.14) 

(A.15) 

Since the maximum inverse length scale becomes exponentially large, 

the function (A.ll) is completely unsuitable. 

w = tanz 

The appropriate function for real z is obtained from equation(A.4)as 

The relation of B to Ax is obtained from equation (A.l) as 

B = Ax/sin Ax. 

The inverse length scale L-’ is found to be tS 

L-~~~ = 2 tan (A x/2) I sin [Ax (2t - 111. 

For B >~/2, the maximum value of L -1 tS is given by 

l) = 2 tan AX/^). (L- tc max 

(A.16) 

(A.17) 

(A.18) 

(A.19) 

30 
For large By Ax + ry and equation (A.17) can be rewritten as 

-+ (r/Z)tan (Ax/.?). Ax 

2 sin (AX/~)COS (a~/2) B= 

We thus find that 

) + 4B/al (L- t< max 

and consequently the funct 

(A.20) 

(A.21) 

on (A.16) is also not suitab e. 

The corresponding function for imaginary z is obtained from equation (A.5) 

as 

1/2 - tanh[Ay[(-l/;)l 

2 tanh Ay/2 . tThe parameter B is now related to Ay by 

sinh Ay B= 

AY * 

The inverse length scale L-l tS becomes 

L-’ = 2Ay I tanh [Ay(<-1/2)] I, t5 

and has the maximum value 

(L-ltS)max = 2Ay tanh (Ay/2). 

For large B, equation (A.23) can be inverted approximately to 

Ay+log (2B log B). 

(A more accurate relation is derived in Appendix B.) For suff 

(A.22) 

(A.23) 

(A.24) 

(A.25) 

yield 

(A.26) 

ciently 

large B, although ~y is logarithmically of 0(1), it is large enough for 

tanh (Ay/2) + 1. Consequently we obtain 

31 
) (L- t< max -+2Ay-2 log (28 log B). (A.27) 

-1 Thus L tS remains logarithmically of 0(1), even when B becomes 

extremely large. 

Condition (A.3) requires the determination of c-1 

(For example, if B = llO1,(L-lt~)max = 20.) 

, which is given tS 

by 

(A.28) 

It follows readily for large B that 

(A.29) (C -l tS)max = (L-l tS 1 max 

and condition (A.3) is also satisfied. 

suitable candidate for a stretching function. 

The function (A.22) is thus a 

w = erfz 

The appropriate function for real z is obtained from equation (A.5) as 

Using equation (A.l), we obtain the relation 

n 

The inverse length sca e L-’’~~ has the simple form 

1/21 Y 

and has the maximum value 

(A.30) 

(A.31) 

(A.32) 

2 ) =Ax. (L- t< max 

(A.33) 

32 
Approximate inversion of equation (A.31) for large B results in 

Ax2+ 4 log [2B(log B/n) 1/21 . (A.34) 

is logarithmically of 0(1), and can also show t5 We again find that L-l 

that(T-’ )max 2 (L-’ )max. The error function is not as simple as the 

hyperbolic tangent, and is not invertible. 

of B, the maximum inverse length scale for the error function is larger 

than that for the hyperbolic tangent. 

candidate for a simple, anti-symmetric, two-sided stretching function. 

t5 t5 

Furthermore, for a given value 

Thus the function (A.22) is the best 

33 
APPENDIX B 

INVERSION OF y = sinh x/x 

Given the function 

y = sinh x/x, 

we seek an approximate inverse function x = f(y) for y>l. 

of the approximation will be measured by calculating the fractional 

error in y, given by 

The accuracy 

error (y) = -1 -1. 

An expansion for f(y) near y = 1 is readily obtained, but it has a 

small radius of convergence. An asymptotic expansion for large y is 

found to converge very slowly. 

is required extends to y -100. A maximum absolute error of .0005 is 

probably sufficient for the numerical applications of this function. 

Hopefully, these criteria can be satisfied by matching appropriate 

expansions for small and large y with the exact solution at some 

intermediate point yl. Since the asymptotic expansion has slow 

convergence, it will also be matched to the exact solution at another 

point y2. 

now be presented. 

The range of y for which a high accuracy 

An outline of the derivation of the two expansions will 

For given values of x1 and x2, the corresponding values of yl, y2, 

(dx/dy)l ,(d x/dy )l ,and (dx/dy)* are obtained from equation (B.l) 

by substitution and implicit differentiation. 

expansion for y<yI. Since y- 1 + x2/6 for small x, f(y)JJ60 for 

y near one. If we introduce the variable 

22 

We first consider the 

y=y-1, 

we can write an expansion for x in the form 

34 
The coefficients A1 and A2 are determined by inverting the expansion of 

equation (B.l) for small x and equating coefficients of 4. The results 

are 

A == -.15 1 (B.5) 

and 

A=- 3x107 

2 32x175. 

22 By differentiating equation (B.4), we can set x, dx/dy, and d x/dy equal 

to their known exact values at yl. 

algebraic equations are readily solved for Aj, A4 and A5. 

The resultant three simultaneous linear 

In order to obtain the expansion for y>yl,we must first determine the 

leading terms in the asymptotic expansion for large x. Equation (B.l) 

can be rewritten as 

x = sinh-' (xy). 

Since sinh-' z - log (2z) asymptotically for large z, we obtain the 

asymp to t i c express i on 

x - v+log (2x) 

where 

v = log y. 

35 
Substituting the zeroth order approximation 

X(0L v 

(B.10) 

into equation (B.8), we obtain the next approximation 

x(? v + log (2v). 

The further substitution of x(l) into equation (B.:?,) results t~ 

v + log [2v(l + log(2v)/v)] 

-v + log(2v) + log (2v)/v 

or x(2)-v+ (1 + l/v) log (2v). 

(B.11) 

6.12) 

The asymptotic approximation (B. 12) agrees very we1 1 with the exact 

solution for yr 1, giving a maximum absolute error of ,fj2 occurring at 

y- 200. The next higher order approximation gives poorer results in 

our range of interest. In order to obtain better accuracy, and join 

the solution with the expansion for small y, we add a polynominal in 

inverse powers of y. 

y2, it is more reasonable to define instead the variable 

Since we will also match exact conditions at 

w = l/y - 1/Y2. (B.13) 

The expansion for x is thus taken to be 

x = v + (l+l/v) log (2v) + Bo + Blw + B2w 2 + B3w3 + B4w 4 - (B. 14) 

The coefficients Bo and B1 are determined by equating x and dxldy to 

their known exact values at y2. 

B4 are obtained by equating x, dxldy, and d x/dy to their known values at 

Similarly, the coefficients B2, B3 and 

22 

Yl. 

36 
For fixed y1 and y2, as y increase above one, the error determined by 

equations (B.2) and(B.4) becomes positive, reaches a maximum, and 

decreases to zero at yl. 

becomes positive again beyond yl, reaches a maximum, decreases through 

zero, reaches a minimum (which is negative), and increases to zero 

at y2. Beyond y2 the error grows monotonically negative. 

absolute error reaches a local niaximum three times in the range 

1 <y<y2. 

three maxima equal. By a trial and error procedure, this condition was 

found for the values x1 = 2.722567 and x2 = 6.05012, corresponding to 

y, = 2.7829681 and y2 = 35.053980. 

error is .000267732, which is more than sufficient. 

reached at y = 69.64. 

reaching .0006 at y- 100 and .00083 at y = 120.5. Thus, the desired 

criteria are essentially satisfied by the two approximations. The 

numerical values for the final coefficients are given by equations 

(65) and (67). 

Using equation (B.14), we find that the error 

Thus the 

The optimum choice of y1 and y2 is that which makes these 

The corresponding maximum absolute 

This value is again 

The magnitude of the error increases beyond this, 

37 
APPENDIX C 

INVERSION OF y = sin x/x 

Given the function 

y = sin x/x, (C.1) 

we seek an approximate inverse function x = g(y) in the range O<y< 1. 

The accuracy of the approximation will be measured by calculating the 

fractional error in y, given by 

The approximation is derived by matching appropriate expansions for y 

near zero and y near one with the exact solution at some intermediate 

point yl. The derivation is briefly outlined below. 

22 For a given xl, the corresponding values of yl, (dx/dy)l ,and (d x/dy )1, 

are obtained from equations (C.l) by substitution and implicit 

differentiation. 

equation (C.1) near y = 0, as a polynomical in powers of T-x. 

the inversion of this expansion, we take for an approximation the 

expression 

The approximation for y<yl utilizes the expansion of 

Based on 

(C.3) 2 2 3 4 5 6 x = ~[l - y+y -(l+ ~/6) y-t B4y + B5y -t B6y 3. 

22 

The resultant three simultaneous 

By differentiating equation (C.3), we can set x, dx/dy, and d x/dy 

to their known exact values at y. 

linear algebraic equations are readily solved for B4, B5 and B6. 

equal 

The approximation for y>yl similarly utilizes the expansion of equation 

(C.1) near y = 1. Introducing the variable 

y=1-y, 

3c 
we obtain an approximation in the form 

X =,m (1 + AIY + A2y2 + A3y3 + A4j4 + A5j5). 

The coefficients A1 and A2 obtained from the expansion near y = 1 are 

A, = .15 

and 

3x107 

2 32x175 . A =- 

The coefficients A3, A4 and A5 are obtained by equating x, dxldy, and 

d xldy to their known values at yl. 22 

The error defined by equation (C.2) reaches a local maximum between y = 0 

and y = y 

y = 1. 

equal. 

(corresponding to y1 = .26938972). 

which is more than sufficient for numerical calculations. The numerical 

values for the final coefficients are given by equations (71) and (72). 

and a local minimum (which is negative) between y = y1 and 

The optimum choice of y1 will make the magnitudes of these extrema 

This was found by trial and error to be produced by x1 = 2.428464 

1’ 

The maximum absol Ute error is .000197170, 

39 
APPENDIX D 

SOLUTION CHARACTERISTICS OF THE TWO-SIDED STRETCHING FUNCTION 

DERIVED FROM w = tanz 

In this appendix we study the two-sided stretching function derived 

from 

w = tanz, 

where z is the complex variable 

z = x + iy. 

If zo and zl are the two end points of the solution in the z-plane 

which define the range 

0, AZ = Z1 - Z 

it is shown in Section I11 that the appropriate variables for the 

stretching function are defined as 

6 = (Z - z0)/Az 

(D.3) 

(D.4) 

and 

t = (tanz - tanzo)/(tanzl - tanzo) 9 0.5) 

while the governing parameters A and B are related to zo and zl, through 

B = sinAz/Az (D.6) 

and 

coszo 

COSZl 

A=- . 

40 

(D.7) 
~ 

Using equation (D.3), A can also be expressed in the form 

1/A = cosAz - tanzO sinAz . (D.8) 

It is also useful to express tanzO as a complex number in the form ~ 

sin 2x0 + isinh 2yo 

tanzO = cos 2x0 + cosh 2yo 

We will first study the real representations of the solution for 

various ranges of A and B. 

Az is real. 

that tanzo, z 

scaled portion of 

If B< 1, it follows from equation (D.6) that 

Using equations (D.8), (D.9), and (D.4), one can prove 

and z are all real. 0’ The stretching function is thus a 

w = tanx. 

(0.10) 

This representation breaks down as B-tl, since Az-tO. 

from equation (D.8) and (D.7) that z++IT/~. 

representation in the limit by introducing 

For AZ1, it follows 

We can obtain the correct 

I 

- 

x = x TIT/2, (D.ll) 

where X30. Since tan (X+IT/~)= -cot X=-l/x as X-tO, it follows that 

for B = 1 the stretching function is a scaled portion of one branch of 

the rectangular hyperbol a 

I 

w = -l/L (D.12) 

If B>1, it follows from equation (D.6) that Az is pure imaginary, 

i.e. Ax = 0. 

the real part of z is constant. 

to be pure imaginary, one can deduce from equation (D.9) that x I IT/^, 0, 

or +n/2. 

is x = 0, and equation (D.l) becomes 

w = itanhy. 

From equation (D.4) it then follows that x = xo, so that 

Since equation (D.8) required tanzO 

If A is sufficiently close to one, the appropriate solution 

i 

i 

(D.13) 

41 
Thus the stretching function is a scaled portion of the hyperbolic 

tangent. As A departs from one, this representation must break down. 

One can easily show from equation (D.8) that as A-+exp(-Ay), yo++ 03 

and as A+exp(Ay), yo + - 

equation (D.13) takes the limiting form 

. Since lyl+ain these two cases, 

The stretching function then becomes a scaled portion of an exponential. 

This limit is reached when B attains the critical value 

B* = sinh llog AL 

log A 

(D.15) 

For A near one, equation (D.15) has the approximate form 

log B* = 1/6 (log A) 2 , (D.16) 

while for A >> 1 and l/A>> 1 we obtain the asymtotic form 

For a fixed A 1, as B decreases below B? the appropriate solution for 

z is 

(D.18) 

The sign of the correct branch follows from equation (D.ll) by taking 

the limit B+l. Equation (D.l) then becomes 

w = icothy. (D.19) 

Thus, for 1 <B<B*, the stretching function is a scaled portion of one 

branch of the hyperbolic cotangent. The boundaries for the various 

real representations of the solution art. shown in figure D.l, which is 

a log-log plot of the B vs A plane. 

42 
Another interesting property is to determine the conditions under 

which the solution curve contains an inflection point. This is certainly 

true for the anti-symmetric curve (A = 1). As A departs from one, the 

inflection point moves towards one end, and eventually could disappear. 

The critical point is reached when either zl= 0 or zo = 0. From equations 

(D.3) and (D.7), it follows that the former condition occurs when 

A = cosAz. (D.20) 

The inflection point thus occurs at the right end of the curve when B 

attains the critical value 

and 

Bt =/= / cosh-’ A (A> 1) 

B+ =,/z / cos-’ A (A< 1). 

(D. 21 a) 

(D.21b) 

For A near one, equation (D.21) has the approximate form 

+ 

~OCJ B = 1/3 log A, 

while for A>>1 we obtain from equation (D.21a) the asymtotic form 

+ log B +log 2A - log (2 log 2A). 

+ 

It also follows from equation (D.21b) that A+O, B 

minimum value 

attains the 

Bmin = 2/~ . 

The condition zo = 0 is found to occur when 

(D.22) 

(D.23) 

(D.24) 

A = l/cosAZ. (D.25) 

43 
Thus the inflection point occurs at the left end of the curve when B 

attains the critical value 

2 6- =,/l - 1/(A ) / cos-’ (1/A) (A>1) 

and 

6- =/sl / cosh-I (1/A) (A<l). 

For A near one, equation (D.26) has the approximate form 

while for A << 1, we obtain from equation (D.26b) the asymtotic form 

(D. 26a) 

(D. 26b) 

(D.27) 

Condition (D.24) is again reached when A+ m. 

From the above relations, one can conclude that the solution curve 

contains an inflection point unless B lies between B and B-. 

Alternatively, an inflection point will be present if l/coshAy < A< coshAy 

for B > 1, and cosAx < A< l/cosAx for B< 1. 

always occur if B < 2/~ . 

has inflection points are also shown in figure D.l. 

regions do not cover completely the regions where solutions are based on 

the tangent and hyperbolic tangent. 

which a solution curve based on the tangent and hyperbolic tangent does 

not contain an inflection point. 

t 

An inflection point must 

The regions of the B-A Plane where the solution curve 

Note that these 

Thus there are narrow ranges for 

Still another way to study the solution is to follow its path in the 

complex z-plane. 

for constant B, as A increases from zero to infinity. It is easy to 

show that as A+O, XO+-r/2 for all solutions. 

z+~/2 for all solutions. The solution paths are shown in figure D.2, 

This is best done by finding the locus of the solution 

Similarly, as A+m, 

44 
with a dashed l i n e f o r B < 1 and solid l i n e f o r B > 1. The details of 

the solution histories are presented below. 

Since z is real, the solution remains on the real axis. When A = 0, 

xu = -.rr/2, and x1 = -.rr/2 t Ax. (If B > &'T, x1 < 0, and the solution 

has no inflection point). As A increases, both xo and x1 increase. 

(If B>~/TT, x1 = 0 and an inflection point appears when A = cosAx). 

When A = 1, xo = -Ax/2, x1 = Ax/2, and we have the anti-symmetric 

soluiion. (If B>Z/.rr, xo = 0 and the inflection point disappears 

when A = l/cosAx). When A reaches +a, x1 = v/Z, and xo = v/2 - Ax. 

With Az = iAy, the solution at A = 0 has zo = -.rr/2, and z1 = -~/2 + iAy. 

As A increases, the solution stays on the line x = -v/2, and both yo 

and y1 increase. When A = exp(-Ay), yo = y1 = +. 

beyond exp(-Ay), the solution returns along the imaginary axis, and 

both yo and y1 decrease from +a. When A = l/coshAy, yo = 0, and an 

inflection point first appears. 

we have the anti-symmetric solution. When A = coshAy, y1 = 0, and the 

inflection point disappears. 

increases beyond exp(Ay), the solution returns along the line x = n/2, 

and both yo and y1 increase from -m. 

zo = .rr/2 - iAy. 

As A increases 

When A = 1, yo = -Ay/2, y1 = Ay/Z, and 

When - A = exp(Ay), yo = y1 - -00. As A 

When A reaches +a, z1 = .rr/2, and 

B=l 

When A = 0, zo = z1 = -.rr/Z. 

point. 

jump to +.rr/2, where they remain as A increases to +a. 

the boundary between B > 1 and B < 1, the solution can occur only at these 

points where the B> 1 and B< 1 solutions intersect. 

D.2, the solid line (B>l) and dashed line (B(1) intersect precisely 

at the three points defined above. 

As A increases, the solution remains at that 

When A = 1 , zo and z1 jump to the origin. When A > 1 , zo and z1 

Since B = 1 is 

As shown in figure 

45 
N 

z 

Ia 

3 

z 

0 

LL 

ll 

a 

46 
X 

> 

a 

(0 

II 

a 

> I 

a, I 

Q 

-I 

II V 

1 )8 a a,o > 

>" 8< 

a I 

Q, I 

I 

II 

a 

.. 

cu 

> 

a 

I 

Q) 0 

V II 

d a 

8< 

I 

47 
APPENDIX E 

EVALUATION OF ANTI-SYMMETRIC INTERIOR STRETCHING FUNCTIONS 

In this appendix we evaluate several candidates for an anti-symmetric 

interior stretching function <(t), where 5 and t are normalized 

variables ranging from zero to one. The slope at the midpoint will 

be designated as 

and only the case B>1 will be considered. The criteria to be 

satisfied are that the inverse length scales L-l and L (defined 

in Appendix A, equations (A.2) and (A.3)) be of order one, even if B 

is very large. We seek a solution that is a scaled portion of a 

universal odd function w(z). In terms of the range Az, expressions 

foreor t are given by equations (A.4) or (A.5). 

--1 

tS tS 

w = tanz 

If we let z be real, so that AZ = AX, the appropriate function for B>1 

is obtained from equation (A.5) as 

Using equation (B.l) we obtain the relation 

takes the form tE The inverse length scale L-' 

L-l t( = 2Ax Itan [Ax((-l/Z)]l, 

and has a maximum value 

(L-' tS ) max = ZAX tan (A~/z). 

(E.4) 

(E.5) 

43 
For large B, Ax-, TI and we obtain 

Since the inverse length scale is proportioned to B, the function (E.2) 

is not suitable. 

If we let z be imaginary, so that Az = iAy, the appropriate function 

for B>1 is obtained from equation (A.4) as 

5- 1/2 = y$p$-#J . 

Applying equation (E.l), we obtain the relation 

The inverse length scale L-l is now tS 

L- ' tS= 2 tanh (Ay/2) lsinh [Ay(Zt-l)]I, 

and has a maximum value 

(L-lt.)max = 2 tanh (Ay/2) sinh Ay. 

For large B, Ay/2+B, tanh(A.y/2) +l, and we obtain 

26 ) + 2 sinh 2B-e . (L- tS max 

(E.7) 

(E.9) 

(E.lO) 

(E.ll) 

Since the inverse length scale grows exponentially with B, the function 

(E.7) is completely unsuitable. 

w = sinz 

The appropriate function for real z is obtained from equation (A.4) 

as 

49 
The relation of B to Ax is obtained from equation (E.l) as 

B== Ax/2 

sin Ax/2) . 

(E.12) 

(E.13) 

Since Ax 5 if function (E.12) is to remain monotonic, it follows from 

equation (E.13) that B 5 ~/2. Therefore a monotonic stretching function 

for large B cannot be obtained using function (E.12). 

The corresponding function for imaginary z is obtained from equation (A.4) 

as 

The parameter B is now related to Ay by 

s i n h ( Ay / 2 ) 

AY/2 - 

B= 

The inverse length scale L-l tS becomes 

L-ltS = Ay tanh [Ay(S-1/2)], 

and has the maximum value 

) = Ay tanh (Ay/2). (L- t< max 

For large 9, equation (E.15) can be inverted approximately to yield 

Ay+2 log (29 log 9). 

For sufficiently large B, ay is large enough for tanh (ay/2) +l. 

Consequently we obtain 

1 +Ay-2 log (29 log B). ( L-t<)max 

50 

(E.14) 

(E.15) 

(E.16) 

(E.17) 

(E.18) 

(E.19) 
Thus L-ltS remains logarithmically of 0(1), even for very large B. The 

inverse length scale I-' tS becomes 

It follows that for large B 

(E.20) 

(E.21) 

Since both inverse length scales remain of 0(1), the function (E.14) is 

a suitable candidate for a stretching function. 

51 
1. Report No. 

NASA CR-3313 

7. Authoris) 

2. Government Accession No. 3. Recipient's Catalog No. 

~ ~~ ~ ~~ ~ I 8. Performing Organization Report No 

4. Title and Subtitle 

On One-Dimensional Stretching Functions for FiniteDifference Calculations 

5. Report Date 

October 1980 

6. Performing Organization Code 

16. Abstract 

Marcel Vinokur, Principal Investigator 

School of Engineering 

Santa Clara, California 

12. Sponsoring Agency Name and Address 

National Aeronautics and Space Administration 

Washington, D. C. 20546 

, 

9. Performing Organization Name and Address 

Engineering and Applied Science Research 

University of Santa Clara 

The class of one-dimensional stretching functions used in finite-difference calculations is 

studied. 

criteria for a stretching function are derived using a truncation error analysis. 

criteria are used to investigate two types of stretching functions. 

stretching function, for which the location and slope of an interior clustering region are 

specified. 

on the inverse hyperbolic sine. The other type of 

function is a two-sided stretching function, for which the arbitrary slopes at the two ends 

of the one-dimensional interval are specified. The simplest such general function is found 

to be one based on the inverse tangent. 

and greater than one was first employed by Roberts. The general two-sided function has 

many applications in the construction of finite-difference grids. 

application is found in one of the references. 

For solutions containing a highly localized region of rapid variation, simple 

One is an interior 

The simplest such function satisfying the criteria is found to be one based 

These 

It was first employed by Thomas et al. 

The special case where the slopes were both equal 

An example of such an 

10. Work Unit No. 

505-31-11-06 

11. Contract or Grant No. 

NSG- 2086 

13. Type of Report and Period Covered 

Contractor Report 

July 1, 1978 to June 30, 1979 

14. Sponsoring Agency Code 

7. Key Words (Suggested by Author(s) J 18. Distribution Statement 

9. Security Classif. (of this report) 

Unclassified 

Stretching Functions, 

Grid Clustering, and 

Finite-Difference Grids 

20. Security Classif. (of this page) 21. NO. of Pages 22. Price" 

Unclassified 55 A04 

Unclassified - Unlimited 

Star Category - 64 

'For sale by the National Technical Information Service, Springfield, Virginia 22161 

NASA-Langley, 1980