HraWnrefl 
 
 H 
 
 ■HHHIHI 
 
 ■ »s hi 
 
 
 
 
 
 4H 
 
 ■ 
 
 IjlR; >•■•»«' I 
 
 
 H 
 
 ■I 
 
 feri 
 
 
 ■ 
 
 
 ' v j' 
 
 
 
 IH MOB 
 
 Hh3I8i9w nil 
 HSssiI 
 
 
 HkBhWh 
 
LIBRARY OF THE 
 
 UNIVERSITY OF ILLINOIS 
 
 AT URBANA-CHAMPAICN 
 
 510. 8 + 
 co p. 2 
 
-)/o. re/ 
 
 /ifaJZh 
 
 W" Report No. UIUCDCS-R-?^- 670 
 
 TRANSIENT VISCOUS FLOW AROUND AN ELLIPTIC CYLINDER 
 
 by 
 
 July 197^ 
 
 Gary Arthur ,Roediger 
 
 ^ 
 
 ,\ 
 
 * 
 
 \ in 1 ^ 
 
 
 
 DEPARTMENT OF COMPUTER SCIENCE 
 UNIVERSITY OF ILLINOIS AT URBANA-CHAMPAIGN 
 
 URBANA, ILLINOIS 
 
Digitized by the Internet Archive 
 in 2013 
 
 http://archive.org/details/transientviscous670roed 
 
Report No. UIUCDCS-R-T^- 67O 
 
 TRANSIENT VISCOUS FLOW AROUND AN ELLIPTIC CYLINDER 
 
 by 
 
 Gary Arthur Roediger 
 
 July, 197^ 
 
 DEPARTMENT OF COMPUTER SCIENCE 
 UNIVERSITY OF ILLINOIS AT URBANA -CHAMPAIGN 
 URBANA, ILLINOIS 6l801 
 
Ill 
 
 ACKNOWLEDGMENT 
 
 I wish to thank my thesis advisor, Professor D. Watanabe, for 
 his guidance and valuable discussions in this research and for his many 
 suggestions in preparing the final manuscript. I am indebted to and 
 express special thanks to Professor S. M. Yen for suggesting and 
 supporting this research, for his guidance, encouragement, and many 
 helpful suggestions. 
 
 I would also like to recognize Rich Flood who first wrote the 
 Navier-Stokes equations in conservative form for a fixed elliptical 
 cylindrical system and did initial programming in that respect and thanks 
 to Terry Akai for his many helpful comments. 
 
 I express deep gratitude to the Coordinated Science Laboratory 
 whose facilities made this research and manuscript possible and to 
 Trudy Williams who typed this thesis. 
 
 To my wife, Barbara, I express my thanks for her encouragement 
 and perseverance. 
 
 Most of all I wish to recognize my Lord and Savior, Jesus Christ, 
 as my source of all ability and guide to all opportunities. 
 
IV 
 
 TABLE OF CONTENTS 
 
 Page 
 
 1 . INTRODUCTION 1 
 
 2 . DESCRIPTION OF PROBLEM AND METHOD OF SOLUTION 3 
 
 2.1. Formulation of the Problem 3 
 
 2 .2 . Adaptive Coordinate System 4 
 
 2.3. Coordinate System 7 
 
 2.4. Node Distribution 10 
 
 3 . NUMERICAL METHOD AND IMPLEMENTATION 14 
 
 3.1. Exterior Flow Field 15 
 
 3.2. Flow at the Body 17 
 
 3.3. Computer Simulation 19 
 
 4 . NUMERICAL RESULTS 22 
 
 APPENDIX 40 
 
 A. Navier-Stokes Equations, Accelerated Reference Frame 40 
 
 B. Navier-Stokes Equation in Elliptic Coordinates 44 
 
 LIST OF REFERENCES 47 
 
1. INTRODUCTION 
 
 Computational fluid dynamics is a new branch of physical science 
 which has developed into a powerful tool with the advent of large scale 
 digital computers. The complex gas dynamic equations can now be solved 
 numerically to simulate gas flows under controlled flow conditions. 
 Computers have indeed become "numerical wind-tunnels." Roache [1] presents 
 a survey of previous work in this field. 
 
 This paper presents a study of the transient flow of a compressible 
 viscous gas around an elliptic cylinder which smoothly accelerates to a 
 fixed velocity. The complete time -de pendent Navier-Stokes equations are 
 solved numerically in the body reference frame using a non-orthogonal time- 
 dependent curvilinear coordinate system designed to describe the body 
 exactly and to resolve the important features of the transient flow. The 
 computation is performed in a fixed mesh which maps into a grid in physical 
 space which evolves to accommodate the transient flow. 
 
 Chapter 2 presents the formulation of the problem and a qualitative 
 description of the transient flow, and outlines the development of a flexible 
 quasi-elliptical time -de pendent coordinate system appropriate to the 
 problem. Chapter 3 presents a description of the second order accurate 
 numerical scheme used and a discussion of the computer implementation of 
 the scheme. Finally Chapter 4 presents numerical results for Mach number 
 .169 and local Reynolds number 40 for two cylinders, a nearly circular 
 cylinder of eccentricity .2, and a 2 to 1 elliptic cylinder of eccentricity 
 .866. The results illustrate the development of the boundary layer and 
 
separation and circulation regions and the evolution of the pressure and 
 flow velocity. The flow velocity field is shown in fixed and accelerated 
 
 frames . 
 
2. DESCRIPTION OF PROBLEM AND METHOD OF SOLUTION 
 
 2.1. Formulation of the Problem 
 
 Consider an infinite elliptic cylinder with its major and minor 
 axes aligned with the x and y axes centered at the origin and surrounded 
 by a quiescent compressible viscous gas of infinite volume. The cylinder 
 starts from rest, accelerates smoothly until it reaches a given final 
 velocity, and then continues moving at the final velocity. 
 
 The flow produced by the cylinder should have the following 
 features. As the cylinder starts from rest, a weak wave is produced which 
 propagates at the speed of sound. As the ellipse accelerates, mass 
 accumulates in front of the body and cavitation occurs in the rear. When 
 the final velocity is reached, the mass accumulated in front of the body 
 produces a compression wave. Similarly, the void in the rear produces 
 an expansion wave. This phenomenon is very similar to the one dimensional 
 piston problem[2]. As these waves propagate, they interact in a complex 
 way, and eventually a steady state is reached near the body. Our 
 objective is to simulate this flow by solving the complete time -de pendent 
 Navier-Stokes equations in the body reference frame. 
 
 The smooth acceleration of the cylinder permits the simulation 
 of the transient gas flow as well as the steady flow around the body. 
 In previous studies of exterior flows described by the Navier-Stokes 
 equations, the body was either started impulsively, or the steady inviscid 
 flow was used as the initial flow. These studies were aimed at only the 
 steady flow solution because the transient flow produced by these initial 
 
conditions is physically meaningless. An exception is the study of the 
 transient inviscid flow around a circular cylinder by Moretti[3]. 
 
 2.2. Adaptive Coordinate System 
 
 Cartesian coordinates are best suited for rectangular bodies 
 because the grid can be designed so that the body coincides with grid 
 lines. For curved bodies, however, a rational choice of body nodes leads 
 generally to an irregular Cartesian grid, and the body is only approximated 
 by a piecewise linear curve joining the body nodes. The accuracy of 
 this approximation can be increased by increasing the number of body nodes, 
 but the price paid is a large number of unnecessary nodes and further 
 irregularity in the grid. The situation is further complicated by the 
 existence of a boundary layer which requires additional nodes near the 
 body for resolution. 
 
 These problems can be resolved through the use of a curvilinear 
 coordinate system in which the body coincides with grid lines. The 
 natural system for this study is the orthogonal elliptical cylindrical 
 coordinate system (u,v) in the accelerated body reference frame. This 
 coordinate system has other advantages besides exact body representation. 
 A grid in (u,v) with uniform u and v spacing maps into a grid in 
 physical space whose nodes cluster near the body. The concentration of 
 points is highest where the curvature is greatest. These two characteristics 
 of the coordinate system combine to place more nodes in physical space 
 where the gradients in the flow variables are greatest. A typical grid 
 corresponding to a uniform (u,v) grid for an ellipse of eccentricity .866 
 
is shown in Figure 1. 
 
 This elliptical coordinate system would be ideal for monitoring 
 small perturbations from the steady state solution. However, the use of 
 this coordinate system for the transient flow studied here leads to 
 unresolvable problems. First, as waves move away from the body, they 
 move into the area where the nodes are widely spaced, and hence they 
 cannot be resolved adequately. Second, at the fixed outer boundary, 
 signals arriving from the body cannot be handled properly. The compression 
 and expansion waves exhibit steep gradients which must be resolved by a 
 locally fine grid, and the signals arriving at the outer boundary must be 
 handled carefully to prevent the creation of spurious waves which can 
 destroy the calculation. Adding nodes in the neighborhood of the waves 
 and at the outer boundary would be one solution, but it is impractical 
 because of computer space and time limitations. 
 
 A solution is to introduce a time-dependent non-orthogonal quasi- 
 elliptical coordinate system designed so that the nodes of a fixed 
 rectangular grid in this system map into points in physical space which 
 move to accommodate the transient flow. First u is mapped into a new 
 variable Q by a bilinear transformation so that the interval along any 
 appropriate v hyperbola from the body to the outer boundary is mapped into 
 [0,1]. Then Q is mapped into the final variable u* using a hyperbolic 
 transformation. The outer boundary initially is an ellipse concentric 
 with the body and of the same eccentricity but slightly larger. In the 
 fixed frame, the outer boundary expands uniformly at the speed of sound, 
 
■ • 
 
 « • 
 
 ■ ■ 
 
 ' . * " . « ■' f/f 
 
 • .".»■■ ft i 
 
 ■ . ■ « i a 
 
 . ■ ' ■ - ■ ". ■ -\-AVL II 
 
 * ■ ■ ■ \VsSS 
 ■ • ". ■. ■■..*■■■ 
 
 ■« ■■ 
 
 a * ■ ■ i 
 
 a a 
 
 a 
 a a 
 
 CO 
 
 <u 
 
 CO 
 
 c 
 
 n 
 o 
 o 
 
 CJ 
 
 4J 
 
 13 
 
 x 
 
 u 
 
 o 
 
 a) 
 
 60 
 •H 
 En 
 
asymptotically approaching a circle centered at the origin. Hence the 
 signals created by the moving body remain within the computational space. 
 As the calculation proceeds the transformation parameters adjust the 
 mapping so that the fixed grid nodes map into points in physical space 
 clustered about the steep gradients in the flow. 
 
 2.3. Coordinate System 
 
 The transformation from the fixed frame (x,y) to the accelerating 
 body frame (x*,y*) is given by 
 
 x* = x 
 
 r w dt , 
 
 d x 
 
 t 
 
 y* = y - J* w dt , 
 o 
 
 where w and w are the x and y velocities of the body in the fixed frame, 
 x y 
 
 By varying w and w , the body can move in any direction at any angle 
 attack and at any speed. The conservative form of the Navier-Sokes equations 
 in the accelerating frame with density, momentum, and energy as the 
 primary variables are presented in Appendix A. 
 
 The transformation from the coordinates (x*,y*) to the elliptical 
 coordinates (u,v) is given by 
 
 x* = c cosh(u) cos(v) , 
 y* = c sinh(u) sin(v) 
 
 The u and v level curves are confocal ellipses and hyperbolas with foci 
 
(-c,0) and (c,0). An ellipse of eccentricity e is represented by the u 
 
 level curve for u , where 
 o 
 
 u q = ln[(l+(l-eV)/ e ] . 
 
 For this ellipse the focal length is c, the semi -major axis is c cosh(u ), 
 
 and the semi -minor axis is c sinh(u ). The coordinate u ranges from 
 
 u to °° and v ranges from to 2tt. Hence this transformation maps the 
 o 
 
 exterior of the ellipse into a semi-infinite strip. The region of 
 interest, however, is a rectangle which lengthens in the u direction with 
 time. The Navier-Stokes equations of Appendix A transformed to elliptical 
 coordinates are presented in Appendix B. In the transformed equations, 
 the terms due to the accelerated reference frame are retained in cartesian 
 form for convenience . 
 
 The transformation from u to Q is given by the bilinear mapping 
 
 C = (bu+c)/(u+a) 
 
 where 
 
 a = [(2u o -u 1 )u oo -u 1 u o ]/(2u 1 -u o -uJ , 
 
 b = ( Ul +a)/[2( Ul -u o )] , 
 
 c = - u b 
 o 
 
 This transformation maps the inner boundary u , an intermediate point u, , 
 and the outer boundary u into 0, 1/2, and 1. Hence it maps the region 
 
of interest into a fixed rectangle. 
 
 The final transformation is given by 
 
 * - I + i_ i r i+(2C-l)tanh(g(t)/2)- | 
 u " 2 lot Ll-(2C-l)tanh(a(t)/2)J ' 
 
 v* = v , 
 
 t* = t 
 
 This transformation preserves the fixed rectangular computational space, 
 but is compresses the region in Q near 1 and and stretches the region 
 near 1/2. The computational space is covered by a fixed uniform 
 rectangular grid which maps into a grid in physical space whose nodes are 
 clustered in the boundary layer, particularly in the regions of greatest 
 body curvature, and near the outer boundary. The parameters u. and a 
 are used to control the distribution of the nodes in physical space. 
 
 The Navier-Stokes equations are solved on this uniform grid. 
 All factors in the equations other than the derivatives of the flow 
 variables can be easily evaluated analytically. The partial derivatives 
 of the flow variables, however, must be approximated by finite differences 
 The first partial derivatives of a function f are written in terms of 
 u* and v* by applying the chain rules 
 
 f - f ,. u* r C , 
 
 u u * fa u 
 
 f = f + f , u* £ 
 v v* u* C v 
 
 f t = f t * + f u * (U *C C t + u *a a t' 
 
10 
 
 The second partial derivatives are transformed by applying these equations 
 a second time. Then the derivatives of flow variables with respect to 
 u*, v*, and t* are approximated by finite differences on the uniform grid. 
 
 2.4. Node Distribution 
 
 In the body reference frame, the compression and expansion waves 
 move at different speeds. Hence the outer boundary at u,,, is a function 
 of v* and t*. In the fixed Cartesian frame, the outer boundary is represented 
 
 by 
 
 \x u +c t / + \y,+c t / 1 ' 
 bo b o 
 
 where c is the speed of sound in the gas at rest, and x, and y, are the 
 o b b 
 
 semi -major and semi -minor axes of the boundary. Transforming this representa- 
 tion to the accelerated frame moving only in the x direction and scaling 
 distances with respect to c cosh(u ) yields the following equation for u^, 
 
 u (v*,t*) = ln[A+(A 2 -l)2] , 
 
 where 
 
 CO 
 
 2 i 
 A = (-C+(C -4BD) 2 )/(2B) , 
 
 2 2 2 
 
 B - (cos (v)E+ sin (v)F)/(cosh (u )E) , 
 
 o 
 
 C = (2 x, cos( v))/cosh(u ) , 
 
 G O 
 
 D .= x, 2 - (F sin 2 (v))/(E cosh 2 (u )) - F , 
 a o 
 
 e = (y b +c o t) 2 , 
 
 F = (x, + c t) 2 
 
 D O 
 
 
11 
 
 Here x, is the nondimensiona 1 distance the ellipse has moved. 
 
 The parameters <*(t*) and u (v*,t*) control the node distribution 
 The parameter u.. is chosen so that 
 
 cosh(u ) - cosh(u ) 
 
 .___ 1_ = \ 
 
 cosh(u^) - cosh(u ) 
 
 where X is a constant between and 1. As X increases nodes are pulled 
 away from the outer boundary; as \ decreases, nodes are pulled away 
 from the body. 
 
 Solving this equation for u,(v*,t*) yields 
 
 Ul (v*,t*) = ln[A+(A 2 - 1)2] , 
 
 where 
 
 A = (l-\)cosh(u m ) + A, cosh(u ) 
 
 o 
 
 Since the outer boundary, u , expands while the total number 
 of nodes remains fixed, the nodes tend to move away from the body. However, 
 the grid in physical space must be fine enough to resolve the boundary 
 layer. Hence a maximum length 6 for the distance in physical space 
 between the body and a first u* node from the body at some value of v* 
 is specified for each simulation. This length 6, obviously is inversely 
 proportional to the Reynolds number. 
 
 The parameter a(t*) is used to control the spacing of the first 
 row of u* nodes from the body. The initial value of a is chosen so that 
 the maximum distance, 6, is attained at a given time t . For times 
 
12 
 
 smaller than t or is held constant and the spacing increases until it 
 equals 6 at time t. . Thereafter, a is chosen to keep the spacing equal 
 to 6. 
 
 By carefully choosing X and 6, the node distribution can be 
 tailor to the problem needs. For example, in a low Reynolds number flow 
 containing strong waves, X should be small to concentrate nodes near 
 the waves. In a high Reynolds number flow containing weak waves, X 
 should be large to concentrate nodes near the body and 6 should be small 
 enough to permit adequate resolution of the boundary layer. 
 
 A typical node distribution in the (u*,v*,t*) space for some 
 time, t > 0, is shown in Figure 2. 
 
■ ■ ■ 
 
 13 
 
 — ■ ■ 
 
 ■ « ■ ",1 i 
 
 . "-■« 
 
 > 
 
 3 
 
 6 
 V 
 4J 
 CO 
 
 co 
 
 <D 
 4-1 
 CO 
 
 C 
 
 •r-l 
 TJ 
 M 
 O 
 O 
 
 o 
 
 4-1 
 
 cu 
 
 C 
 
 cu 
 a 
 cu 
 
 I 
 •H 
 4J 
 
 C 
 
 •r-l 
 
 T3 
 •r-t 
 
 r-l 
 
 o 
 
 CM 
 
 cu 
 
 u 
 
 00 
 •H 
 
 fa 
 
14 
 
 3. NUMERICAL METHOD AND IMPLEMENTATION 
 
 The Navier-Stokes equations in the (u*,v*,t*) coordinate system 
 can be written in the form 
 
 f . - G(f .,f „ , ,f ,,f , , ,f , ..) 
 t* u* u*u* v* v*v* u-"v* 
 
 where 
 
 T 
 f = {p,m ,m ,E) , 
 
 p = density , 
 
 m = u component of momentum , 
 u 
 
 m = v component of momentum , 
 v 
 
 E = stagnation internal energy 
 
 These equations were solved numerically using MacCormack's predictor- 
 corrector method [4]. This scheme was chosen because it is second order 
 accurate, is easy to implement, and has worked well on viscous flow 
 problems . 
 
 Because of computer space and time limitations, the angle of 
 attack is set equal to to exploit the symmetry of the flow about the x-axis 
 The body velocity w is described by the cubic spline with knots 0, t f /2, 
 and t satisfying 
 
 w (0) = , 
 x 
 
 w x (t £ /2) = v £ /2 , 
 w x (t) - v £ , t > t £ - . 
 
15 
 
 This spline is continuous, has two continuous derivatives, and is anti- 
 symmetric about (v_/2,t_/2). 
 
 The computational space in (u*,v*) is the rectangle [0,1] x [0,tt] 
 
 The grid used is {u*,v*} where 
 
 i j 
 
 u* - iAu*, 1 ■ 0, 1, ..., N , 
 i 
 
 v* = jAv*, j - 0, 1 N , 
 
 j 
 
 where Au* = l/M and Av* = tt/N. The time is partitioned by the grid {t*} 
 
 n 
 
 where 
 
 e*-0 
 
 t * = t* + At*, n = 0, 1, . . . . 
 n+1 n n 
 
 The value of any function f(u*,v*,t*) at the point (u*,v*) at time t* is 
 
 i j n 
 
 denoted by f . . . 
 
 3.1. Exterior Flow Field 
 
 The exterior flow field is defined by {u*,v*}, i = l, ••• , M, and 
 
 i j 
 
 j=0, ... , N. MacCormack's predictor-corrector method first predicts the 
 
 values of the flow variables at time t , , using 
 
 n+1 
 
 _n+l n n 
 
 f. . = f. , + At G. ., i=l, ... , M-l and j=0, ... , N , 
 i,J i.J i,J 
 
 where 
 
 n n n2"2n n 
 
 G. , s G (A f. ., A v *f. . , 6 ^f. ., 6 f n , A ^ ^f ) , 
 i,j u* i,j' i,j' u* i,j' v* i,j' u*v* i,2 
 
16 
 
 A .f . . = (f, ., , - f. .)/Au* , 
 u* i,j i+l, J t,J 
 
 A * f 4 • = ( f - ,xi " f - -)/ Av * . 
 v* i,j i,j+l 1,J 
 
 6 2 ,f. . * (f 4 . - 2f. . + f. ., .)/Au* 2 , 
 u* i,j l-l, j i,j i+l, j 
 
 6 2 f. . = (f. , , - 2f, . + f. .,.)/A v * 2 , 
 
 A ,. ,f. • = (f- L1 «_, - f-^-. - f- ... + f - .)/Au*A v * . 
 U*V* l,j 1+1, j+1 1+1, j 1,J+1 i,j' 
 
 All the terms in G are evaluated at time t* and the predicted values of the 
 
 n 
 
 flow variables at the outer nodes are set to the freestream conditions at 
 
 time t* . .. 
 n+1 
 
 Then the final values of the flow variables at t ,. are obtained 
 
 n+1 
 
 using the following correcter equation, 
 
 n+1 n _n+l _ n +l , 
 
 f. . = (f. . + f. . + AtG. ,)/2, i=l, ... , M-l and 1=0, ... , N, 
 i.J i-.J i,J i>J 
 
 where 
 
 -n+1 -n+i _n+l 2 -n+1 .2 ~n+l „ ^+l x 
 G. , = G(V f. ., V f. ., 6 .f. ., 6 J.. ., V f. .) , 
 i,j v u* i,j' v* i,j' u* i,j' v« i,j' u*v* i,j 
 
 V f , = (f. - f. . .)/Au* , 
 u* l,j i,j i-l, J 
 
 V f. . = (f. . - f. . 
 
 v* 1,J l,j 1,J 
 
 f. J/Av* > 
 
 6 2 f. . = (f. ., . - 2f. . + f - .)/A u * 2 , 
 u* i,j i+l, j i,j i-l, j 
 
17 
 
 6 2 f . = (f. ... - 2f, . + f. . .)/A v * 2 , 
 v* i,j i,j+l i,j i,J-l 
 
 ? +_*,£. , s (f. ~ f- , i ~ f. n • + f. , J/Au*A v * . 
 u*v* i,j i,j i,j"l i-l, J 1-1, J-l' 
 
 — "k 
 
 All the terms in G are calculated at time t ., and the final values of the 
 
 n+1 
 
 flow variables at the outer nodes are set to freestream conditions at time 
 
 'n+1 
 
 Both the predictor and the corrector require the flow variables 
 
 at v* = - Av* and v . = rr + Av*. These values are obtained by symmetry 
 
 from the values at v* and v„ , as follows 
 
 i N-l 
 
 P i,j = P l,j+2 
 
 (m ). . = (m ), .,„ 
 
 (m ) . . = - (m ) . „ 
 vi,j v i,j+2 
 
 E. . = E. , , 
 i,J i.j+2 
 
 for j = - 1 and N-l. 
 
 3.2. Flow at the Body 
 
 A separate calculation is required to compute the flow variables 
 at the body. The two major conditions imposed on the flow at the body are 
 (1) no-slip at the surface in the accelerated frame, and (2) the wall is 
 adiabatic . 
 
18 
 
 The no-slip or zero particle velocity condition is implemented 
 
 by setting m and m to zero at the wall. For an adiabatic wall, the 
 u v 
 
 internal energy of the gas, e, at the surface has zero gradient normal to 
 the surface. 
 
 u body 
 Now the energy of a perfect gas is 
 
 E = pe + p V-V/2 . 
 Since the velocity vector, V, at the wall is zero, it follows that 
 
 (E ), , = (p e). , = (p E/p). . 
 u body u body u body 
 
 To be consistent with MacCormack's second order scheme, all 
 derivatives must be approximated by a second order accurate difference 
 formula at the wall. Hence the following equation is used 
 
 (f ) . = (l/2Au*)(-3f . +4f. . - f, ,)(u*C ) . • 
 u o,j o,j l,j 2, j q u o,j 
 
 Substituting this expression for E and p in the equation already given 
 and solving for E . yields 
 
 E . = p .(4E. . - E .)/(4p 1 . - p .) . 
 o,J o,j v l,j 2,j l,j r 2,j 
 
 In the predictor p , can be calculated using the continuity 
 
 °,J 
 
 equation and the predicted energy, E . , is obtained from the equation for 
 
 ° > J 
 
 E . above using predicted values at time t* , , on the right hand side. 
 o,j ° n+1 
 
19 
 
 The boundary calculations in the corrector present a problem 
 since nodes (u* ,v*) are needed but do not exist. The required derivatives 
 in the continuity equation must be approximated by a second order accurate 
 difference formula consistent with MacCormack's corrector equation. The 
 following approximation can be used. 
 
 (f ) . = (l/Au*)(-2f . + 3f . . - f , ,)(u*C ) . 
 u o,j o,j l,j 2,j q u o,j 
 
 Using this result in the continuity equation yields the values 
 
 of p .. The energy, E . , is then obtained from the equation for E 
 
 o,J o,j o,j 
 
 above using corrected values at time t* . on the right hand side. 
 
 A second method recommended by Roache[5] for computing the density 
 at the boundary was also used. This method, which employs a staggered 
 mesh, provided results very close to the results of the first method. The 
 method required a large amount of geometric computation because of the 
 time dependent grid and moreover required the solution of a tridiagonal 
 system of linear equations for the density values. Since this method for 
 all its added complexity did not produce significantly different results, 
 it was discarded in favor of the first scheme. 
 
 3.3. Computer Simulation 
 
 The complexity of the coordinate system, Navier-Stokes equations, 
 and numerical scheme, the small time steps required for stability, and 
 the relatively large number of nodes required for adequate resolution, 
 all combine to require a very large amount of computer time for each 
 
20 
 
 simulation. Another important factor is the amount of internal storage 
 available. Since many time dependent variables in the predictor are 
 identical to those of the previous corrector, these variables can be 
 computed at each node once and used both in the corrector at t and the 
 predictor at t if storage is available. Unfortunately, storage was at 
 a premium on the machines used, and many quantities had to be recomputed. 
 This added to the already large computing times for each simulation. 
 
 Initially the outer boundary is a finite distance from the body. 
 Hence calculations are performed only at the nodes within the envelope 
 of the wave created by the initial motion of the cylinder. At the remaining 
 nodes, the flow variables are assigned their freestream values. This 
 procedure saves a large amount of computer time during early times, and 
 works better than the alternative procedure of adding nodes using 
 interpolation. 
 
 The numerical scheme was implemented in a Fortran program. Two 
 versions of the program were written: one for the PDP-10 timesharing 
 system, and one for the IBM 360/75 batch system. The basic programs were 
 identical, but the PDP-10 program contains interactive facilities which 
 allow the user to review results dynamically, set checkpoints, and adjust 
 the time step for stability. Both versions automatically adjust the time 
 step to maintain stability. If instability occurs the time step is 
 decreased and the program restarts using the data at the previous stable 
 checkpoint. Checkpoints for program restarts are created by writing the 
 current values of all time dependent variables into a disk file. This pro- 
 cedure minimizes the effect of a machine failure or numerical instability. 
 
21 
 
 The parameters required for a simulation are the reference 
 Reynolds number, eccentricity, grid parameters M and N, node distribution 
 parameters X, 6, x, , and t , acceleration parameters t. and V f , and the 
 initial time step At*. When a run is restarted using a checkpoint of a 
 previous run, At* may be reset. The results are printed on a line 
 printer and plotted, if desired, on a CALCOMP plotter or a Tektronix 
 4010-1 display. 
 
22 
 
 4. NUMERICAL RESULTS 
 
 Calculations were performed for two cylinders, a nearly circular 
 cylinder of eccentricity .2 and a 2 to 1 cylinder of eccentricity .866. 
 In both calculations Au* is 1/31, Av* is tt/37, \ is .68, 6 at v* = tt is 
 .008, the Prandel number is .75, y is 1.4, the reference Reynolds 
 number is 100, the final Mach number is .169 attained at time .5, and 
 the local Reynolds number, 
 
 Re. = 2c cosh(u )p w (t)/u , 
 
 I O O X 
 
 is 40 at Mach number of .169. 
 
 In both calculations the initial time step is .0001. As the 
 minimum spacing between nodes increases to 6, the time step increases to 
 .001 for the first cylinder and .00095 for the second. Because of the 
 small time steps used, each calculation required a large number of steps 
 and consequently a large amount of time. The first required 15.5 hours 
 on thelBM 360/75 to reach time 11.4 in 13,362 steps, while the second 
 required 18 hours to reach time 13.46 in 15,380 steps'. 
 
 As the body accelerates, mass accumulates in front of the body 
 and cavitation occurs behind the body. This phenomenon is evident in 
 the 3-dimensional pressure plots in (u*,v*) shown in Figures 4 and 5, 
 where the pressure is high in the front and low behind. These two 
 figures also illustrate the effects of viscosity since they are from runs 
 differing only in the reference Reynolds number. The viscosity in Figure 5 
 is 10 times that in Figure 4. Consequently the smearing of the compression 
 
23 
 
 Compression 
 Wave 
 
 Expansion 
 Wave 
 
 Pressure Variation 
 Along the Body 
 
 tt + Av* 
 
 NP-429 
 
 Figure 3. Description of 3-dimensian pressure field plots in 
 (u*,v*) space. 
 
24 
 
 Figure 4. 3 -Dimensional pressure plot, eccentricity is .866, 
 Re is 1000, Mach number is .034, and at time .14. 
 
 Figure 5. 3-Dimensional pressure plot, eccentricity is .866, 
 Re is 100, Mach number is .104, and at time .29. 
 
25 
 
 Figure 6. 3-Dimensiona 1 pressure plot, eccentricity is .866, 
 Re is 100, Mach number is .169, and at time 2.95. 
 
 Figure 7. 3-Dimensiona 1 pressure plot, eccentricity is .866, 
 Re is 100, Mach number is .169, and at time 10.25. 
 
26 
 
 Figure 8. 3 -Dimensional pressure plot, eccentricity is .866, 
 Re is 100, Mach number is .169, and at time 13.46. 
 
 Figure 9. 3-Dimensiona 1 pressure plot, eccentricity is .2, 
 
 Re is 100, Mach number is .169, and at time 10.60, 
 
27 
 
 and expansion waves is greater in Figure 5. 
 
 When the final velocity is attained, compression and expansion 
 waves move out from the front and rear of the body at the speed of sound, 
 Vv. As these waves interact, they create complex waves within the gas. 
 These waves can be seen in Figures 6 and 7. Eventually a steady state 
 is reached near the body. Figures 8 and 9 show the steady state pressure 
 distribution near the body as well as the compression and expansion 
 waves far from the body for both cylinders. 
 
 The pressure coefficient is defined as 
 
 C p = 2(p-l)/„ x 2 , 
 
 where p is the dimensionless pressure and w is the dimensionless velocity 
 
 of the body. The evolution of C at the surface of the bodies is shown 
 
 P 
 
 in Figures 10 through 13. Initially the pressure is high in front and 
 low behind, but as the compression and expansion waves travel around the 
 body and interact, the pressure decreases in front and increases behind 
 and the pressure coefficient reaches the steady state shown in Figure 11 
 at time 10.60 for the near circle and in Figure 13 at time 13.46 for the 
 2 to 1 ellipse . 
 
 The density variation throughout the flow field is 7% for the 
 near circle and 5% for the 2 to 1 ellipse. The temperature variation 
 at the body is 27 for the near circle and 1% for the 2 to 1 ellipse. The 
 variations for the near circle are larger than those for the 2 to 1 
 ellipse because the blunter near circle produces stronger compression 
 
28 
 
 72 108 
 V (degrees) 
 
 180 
 
 NP-428 
 
 Figure 10. Pressure coefficient around body of eccentricity .2. 
 
 2.40 
 
 1.60 
 
 0.80 
 
 C P 
 
 -0.80 
 
 -1.60 
 
 -2.4Q 
 
 1 1 1 1 
 
 1 1 ' 
 Time 
 
 1 I 
 
 = 4.76 
 = 7.25 
 
 "" Time ■ 
 
 ^^v 
 
 Time : 
 
 = 10.60 
 
 ^ 
 
 
 — 
 
 Vv\ 
 
 
 
 Vv> 
 
 
 
 \X\ 
 
 
 — 
 
 \ 
 
 
 X ^-a- — 
 
 \ 
 
 
 
 \\v 
 
 
 
 NN^ 
 
 ^-^* 
 
 
 1- ,1 x -^.. 
 
 — ^ 
 
 u.,- ^ 1 I 
 
 i 
 
 36 
 
 72 108 
 V (degrees) 
 
 144 
 
 180 
 
 NP-426 
 
 Figure 11, Pressure coefficients around body of eccentricity .2 
 for various times. Time 10.60 is steady state curve. 
 
29 
 
 C P 
 
 -1.5 
 
 Time = 0.9 
 
 1 r 
 
 72 108 
 V (degrees) 
 
 144 
 
 180 
 
 NP-432 
 
 Figure 12. Pressure coefficient around body of eccentricity 
 .866. 
 
 2.40 
 
 i 1 r 
 
 T 
 
 Time = 3.04 
 
 Time= 9.12 
 
 Time =13.46 
 
 72 108 144 180 
 
 V (degrees) NP . 427 
 
 Pressure coefficients around body of eccentricity 
 .866 for various times. Given 13.46 is steady 
 state curve. 
 
30 
 
 and expansion waves. This difference in relative strengths of the waves is 
 evident in the 3-dimensiona 1 pressure plots in Figures 8 and 9. 
 
 Figures 14 through 22 show the velocity field development in 
 the body reference frame. They illustrate how the boundary layer 
 evolves and how separation and circulation eventually develop. Figures 18 
 and 22 show that separation occurs further downstream for a slender body 
 than for a blunt body. Let s be the arc length along the body from the 
 
 front stagnation point and let I be half the body perimeter. Figures 18 
 
 s 
 and 22 show separation occurring at y = .63 for the near circle and at 
 
 s 
 
 j = .84 for the 2 to 1 ellipse. 
 
 Figures 23 through 29 show the velocity field development in the 
 fixed frame, They illustrate how the accelerating body affects the gas 
 particles. The particles at the surface move in the same direction and at 
 the same speed as the body. The particles in front of the body are pushed 
 out and around the cylinder, while those behind are pulled in toward the 
 body to fill the void. A circulation region develops near v = tt/2 during 
 acceleration, and then moves slowly downstream. Figures 23 to 29 illustrate 
 the development of these regions. 
 
 These results represent nearly incompressible isothermal flow. 
 In principle this method can also be applied to high subsonic Mach number 
 flows where compressible and thermal effects are more significant. However, 
 finer meshes would be required to resolve the steeper gradients of such 
 flows, and smaller time steps would be required for stability. This might 
 raise the computer storage and processor time to prohibitively high levels 
 on many machines . 
 
31 
 
 Figure 14. Initial boundary layer growth at time .29 
 
 FlHfli VtLOtm -0.3000 nEAOCO 01 TIHt-O. 50000 U-N0OC3 l 
 Hnt-3 04419 *CH.O 169031 L8CR. ntT«L03 MUHBM-40 . 00 
 
 RELRTJVE VELOCITY VECTORS RCCELERRTEO FRPME 
 
 Figure 15. Velocity field development at time 3.04 
 
32 
 
 Figure 16. Velocity field development at time 10.25 
 
 STSTivr^ffiTiff'WTora (catepwrm t 
 
 Figure 17. Circulation and separation at time 13.46 
 
33 
 
 / 
 
 v > 
 
 --'//■ 
 
 £f*>' y 
 
 
 Figure 18 
 
 • f circulation region at time 13.46, 
 Enlarged view of circus 
 
34 
 
 nOX 015 TO 1ST soot' O.Q08 QT l 00 0KACC3 OT. ;.C0130 ml. i|S v-0.60 
 FIfWL VCLOClTr.0.2000 "CflCMCO AT TIMC-O. 50000 U-MJDeS 1 TO 30 
 T1HC-0 38437 MflCM.O. 146091 LOWt *CTNQlD3 Ni_*BM-24.S7 
 
 s£_fl:;'.E velocitt vectors acce.e^'ed fra^e 
 
 Figure 19. Initial boundary layer growth at time .38, 
 
 Figure 20. Velocity field development at time 2.69 
 
35 
 
 ^ 
 
 Figure 21. Circulation and separation at time 10.60, 
 
 s 
 
 Figure 22. Enlarged view of circulation region at time 
 10.60. 
 
36 
 
 MAX 015 TO I3T MODE- 0.008 01 180 OtCflCC? 01" 0.0OOBS —1.20 >-0.68 
 F1MO. VCLKnT-O.aOOO ntOC«0 HI TIW0.SO0OO U-M00f3 I 10 28 
 1IHC0.290SH MflCM-O. I0M89H LOCOL KTKOL05 »»f«CIW>i.B2 
 
 RELATIVE VELOCITY VECTORS FIXED FRAME 
 
 Figure 23. Initial gas particle movement at time .29 
 
 MR* 013 TO 131 MOOt: 0.008 PT 180 OCGflCCS 0T« 00053 .-3.85 X-O.EB 
 FIMPL VCLOCITT-0.2000 «£RCHeO BT TJME-0.50000 U-N00E3 1 TO 17 
 TIK-3.0H413 HPCH-O. 169031 LK«. nTTMOLOS MUHBCT-iJO.OO 
 
 RELATIVE VELOCITY VECTORS FIXED FRAME 
 
 Figure 24. Gas particle movement at time 3.04. 
 
Figure 25. Gas particle movement at time 10.25 
 
 37 
 
 SEJRnSrWEsc'i'l! Wflfco rmmi 
 
 Figure 26. Circulation region of gas particles has moved 
 downstream at time 13.46. 
 
38 
 
 Figure 27. Initial gas particle movement at time .38 
 
 
 \ \ 
 
 
 Figure 28. Gas particle movement at time 2.69. 
 
39 
 
 Figure 29. Circulation region of gas particles has moved down- 
 stream at time 10.60. 
 
40 
 
 APPENDIX A 
 Nayier-Stokes Equati ons, Accelerated Reference Frame 
 
 The two dimensional Navier-Stokes equations for an ideal com- 
 pressible viscous gas in an accelerated Cartesian reference frame (x*,y*) 
 are presented here. 
 
 The transformation from the fixed frame (x,y) to the accelerated 
 frame (x*,y*) is given by 
 
 t 
 
 x* = x - f w dt , 
 •' x 
 o 
 
 y* = y - J w dt , 
 
 o y 
 
 where w and w are the x and y velocities of the accelerated frame relative 
 x y 
 
 to the fixed frame . 
 
 The primary dependent variables chosen are the density p, the x* 
 
 and y* components of the momentum m ^ and m ., and the total energy E. 
 
 The momentum and total energy are defined as follows: 
 
 V = p y y * ' 
 
 E = P(e+(V^+V^)/2) , 
 
 where V . and V . are the x* and y* components of the gas velocity, and e 
 is the internal energy. Other dependent variables are the velocity V, 
 temperature T, stress t, pressure P and heat flux of q. 
 
41 
 
 All dimensional quantities will be nondimensionalized with respect 
 to the following characteristic quantities. 
 
 L = length , 
 
 o ■ velocity , 
 
 where 
 
 Ia/y = 
 
 time , 
 c 
 o 
 
 P = density 
 
 PC 2 
 
 o o , 
 = pressure and energy , 
 
 2 
 c 
 
 -— = temperature , 
 
 c = specific heat at constant volume , 
 v 
 
 c = specific heat at constant pressure , 
 P 
 
 c p 
 Y = ratio of specific heats, — , 
 
 v 
 
 R - gas constant, c - c , 
 
 P v 
 
 c = speed of sound in reference gas 
 o 
 
 The scaling yields the following nondimensional quantities 
 
42 
 
 p c L 
 Re ■ reference Reynolds number 
 
 Pr = Prandtl number = -r* , 
 
 M ■ Mach number = — , 
 
 c 
 o 
 
 where k is the thermal conductivity and \i is the viscosity. 
 
 The Navier-Stokes equations in dimensionless form with the stars 
 dropped for convenience are: 
 Continuity 
 
 - dm oni 
 
 & + — * + — £ = 
 
 ot dx By ' 
 
 x -momentum 
 
 dm -. ^ -, 
 
 TT + =T ( m V ) + f- (m V ) = 5- [t -P ] 
 ot ox x x oy x y ox xx 
 
 St dw 
 xy j 
 
 + o7^" p dT 
 
 y -momentum 
 
 Sm v * a * 
 
 TT- + I" (m V ) + f- (m V ) = ~ [T -P] 
 
 ot ox s y x ay y y oy yy 
 
 oT xy dw y 
 
 + Sx~" ■ P dT 
 
43 
 
 Total energy §£ + ^ (EV,) + ^ (EV y ) - - ^ - Jl 
 
 + §- [V (T " ?) + V T ] 
 
 ox x xx y xy 
 
 + j— [V T + V (T -p)] 
 
 oy x yx y yy 
 
 dw dw 
 
 x _y 
 
 x dt y dt 
 
 where 
 
 T §§ Re Lo? Vd5 + oT7 /3 J » 
 
 T §^] " Re LoTl + d§ J ' 
 
 Y ST 
 
 q ? " (Y-DPrRe &5 ' 
 
 for 5 = x,y and M = y,x. 
 
 Note that w and w are relative to the fixed frame. 
 x y 
 
44 
 
 APPENDIX B 
 Navier-Stokes Equation in Elliptic Coordinates 
 
 This appendix presents the Navier-Stokes equations of Appendix A 
 transformed to elliptic cylindrical coordinates. The transformation from 
 the Cartesian coordinates (x,y) to the elliptic coordinates (u,v) is given by 
 
 x = c cosh(u) cos(v) , 
 y = c sinh(u) sin(v) 
 
 The notation and variable definitions are identical to those of Appendix A. 
 However, the reference length L of Appendix A is now assumed to be the 
 semi -major axis of the cylinder, c cosh(u ), and the momentum is now resolved 
 into u and v components. For convenience, the x and y components of the 
 body velocity w and w are retained in Cartesian form. 
 
 The transformed Navier-Stokes equations [6] are 
 
 Continuity 
 
 ft + ^ [§;*%> + t;< h °v] = • 
 
 u -momentum 
 
 Bm 
 
 ■ST 1 = ^r 1?- [h(T -m V -T + m V)] 
 ot , 2 Lgu uu u u vv v v 
 
 + 2 f- [h(T -m V )] + h 
 ov uv u V 
 
 r- [T - P - m V ] 
 .ou W V V 
 
 a 11 dw x dw y 
 
 §- [T -m V ] \Y + Apl- + BP jr- , 
 ov uv u v JJ dt dt 
 
v -momentum 
 
 45 
 
 Total energy 
 
 3m 
 
 a* 1 - -T if" ft(T -m V -t +m V )] 
 ot ,2 ^-ov w v v uu u u 
 
 5T ~ .2 lav [hl 
 
 + 2 t- [h(T - m 
 ou vu v 
 
 V )] + h||- [T - P - m v 1 
 u Ldv uu u uj 
 
 ou vu V 
 
 m V ]]) - 
 V u JJ 
 
 dw dw 
 
 If ■ h ih < h <V T u U 
 
 - P 
 
 E) + V T 
 
 V vu 
 
 Vi 
 
 + f~ [h(V (t - P - 
 OV V w 
 
 E) + V T 
 
 u uv 
 
 V 1 
 
 where 
 
 ( dw x dw y\ ( dw y dw x\ 
 
 + A V m TT~ + m ~a7~) + B V m 7Z m 77/ ' 
 
 \ u dt v dt / \ u dt v dt ' 
 
 2, . . .2, ,,2 
 
 (sinh (u) + sin (v)) 2 /cosh(u ) , 
 
 h 
 
 A = cos(v) 8inh(u)/(h cosh(u )) , 
 
 B = sin(v) cosh(u)/(h cosh(u )) , 
 
 SV. V-, 
 
 . . 2_ { i Zl + 3 Ml 
 
 §1 Re lh a? u 2 ST| 
 
 77 r lf < h V + Iff < h V]} ■ 
 3n 
 
46 
 
 - 1 \ 
 
 m- fc t& ( Y h) + tn ( V h) ) 
 
 ST 
 
 q § " " (Y-l)PrRe h d§ ' 
 
 for 5 ■ u,v and T] ■ v,u. 
 
47 
 
 LIST OF REFERENCES 
 
 1. Roache, P. J., Computational Fluid Dynamics , pp. 226-228, Hermosa 
 Publishers, Albuquerque, New Mexico, (1972). 
 
 2. Liepmann, H. W. and Roshko, A., Elements of Gas Dynamics , p. 91, 
 Eighth Edition, John Wiley & Sons, Inc., (1967). 
 
 3. Moretti, G. , "Transient and Asymptotically Steady Flow of an Inviscid, 
 Compressible Gas Past a Circular Cylinder", Polytechnic Institute 
 
 of Brooklyn, PIBAL Report Number 70-20, (1970). 
 
 4. MacCormack, R. W. , "The Effect of Viscosity in Hypervelocity Impact 
 Cratering", AIAA 7th Aerospace Sciences Meeting, Paper Number 69-354, 
 (1969). 
 
 5. Roache, op. cit . , pp. 272-275. 
 
 6. Woan, C. J. and Akai, T. J., "Basic Gas Dynamic Equations in a Moving 
 Orthogonal Coordinate System", Report T-8, Coordinated Science Laboratory, 
 Urbana, Illinois, (1973). 
 
 7. Schlichting, H. , Boundary Layer Theory . Sixth Edition, McGraw-Hill 
 Book Co., Inc., New York, (1968). 
 
BIBLIOGRAPHIC DATA 
 SHEET 
 
 1. Report No. 
 
 UIUCDCS-R-7U-670 
 
 "itle and Subtitle 
 
 Transient Viscous Flow Around An Elliptic Cylinder 
 
 3. Recipient's Accession No. 
 
 5. Report Date 
 
 July, 197^ 
 
 7. Author(s) 
 
 Gary Arthur Roediger 
 
 8. Performing Organization Rept. 
 N °- UIUCDCS-R-7^-670 
 
 9. Performing Organization Name and Address 
 
 University of Illinois at Urbana-Champaign 
 Department of Computer Science 
 Urbana, Illinois 6l801 
 
 10. Project/Task/Work Unit No. 
 
 1 1. Contract /Grant No. 
 
 12. Sponsoring Organization Name and Address 
 
 University of Illinois at Urbana-Champaign 
 Department of Computer Science 
 Coordinated Science Laboratory 
 Urbana, Illinois 6l801 
 
 13. Type of Report & Period 
 Covered 
 
 Master of Science Thesib 
 
 14. 
 
 15. Supplementary Notes 
 
 16. Abstracts 
 
 This paper presents a study of the transient flow of a compressible viscous gas 
 around an elliptic cylinder which smoothly accelerates to a fixed velocity. The 
 complete time -dependent Navier- Stokes equations are solved numerically in the body 
 reference frame using a non- orthogonal time -dependent curvilinear coordinate system 
 designed to describe the body exactly and to resolve the important features of the 
 transient flow. The computation is performed in a fixed mesh which maps into a grid 
 in physical space which evolves to accommodate the transient flow. 
 
 17. Key Words and Document Analysis. 17a. Descriptors 
 
 17b. Identifiers/Open-Ended Terms 
 
 17c. COSATI Field/Group 
 
 18. Availability Statement 
 
 Unlimited Release 
 
 PORM NTIS-35 ( 10-70) 
 
 19. Security Class (This 
 Report) 
 
 UNCLASSIFIED 
 
 20. Security Class (This 
 Page 
 
 UNCLASSIFIED 
 
 21- No. of Pages 
 
 50 
 
 22. Price 
 
 USCOMM-DC 40329-P7I 
 

 
 uJ 
 CO 
 
UK 
 
 
 UNIVERSITY OF ILLINOIS-URBANA 
 510MIL6flno C002 no 667 672(1974 
 Riport/ 
 
 3 0112 088401382 
 
 
 
 ■ 
 
 ■ 
 
 ■ 
 
 ' ■ 
 
 ■ 
 
 
 ■-V : '■*• 
 
 
 
 
 ■ ' ■ 
 
 
 
 
 . H 
 
 
 ■ 
 
 
 
 
 
 
 X *i 
 
 *■« I 
 
 
 ! ~*ij 
 
 SK9B 
 
 1 ■ I i H 
 
 HH I H 
 
 mi 
 
 ■ 
 
 v.4*