Chapter 1

2 downloads 0 Views 12MB Size Report
Rearranging the above equation as,. (10c). The solution of is then expressed by,. (10d) ...... Appl. Math. ..... tering of a laser light sheet projected through the flow.
he rs

,I

nc .

PHYSICS RESEARCH AND TECHNOLOGY

TURBULENT FLOWS

N

ov

a

Sc i

en

ce

Pu bl

is

PREDICTION, MODELING AND ANALYSIS

PHYSICS RESEARCH AND TECHNOLOGY

,I

he rs

Additional e-books in this series can be found on Nova’s website under the e-books tab.

nc .

Additional books in this series can be found on Nova’s website under the Series tab.

ENGINEERING TOOLS, TECHNIQUES AND TABLES

is

Additional books in this series can be found on Nova’s website under the Series tab.

N

ov

a

Sc i

en

ce

Pu bl

Additional e-books in this series can be found on Nova’s website under the e-books tab.

he rs

,I

nc .

PHYSICS RESEARCH AND TECHNOLOGY

TURBULENT FLOWS

Pu bl

is

PREDICTION, MODELING AND ANALYSIS

ce

ZIED DRISS

N

ov

a

Sc i

en

EDITOR

New York

Copyright © 2013 by Nova Science Publishers, Inc.

nc .

All rights reserved. No part of this book may be reproduced, stored in a retrieval system or transmitted in any form or by any means: electronic, electrostatic, magnetic, tape, mechanical photocopying, recording or otherwise without the written permission of the Publisher.

,I

For permission to use material from this book please contact us: Telephone 631-231-7269; Fax 631-231-8175 Web Site: http://www.novapublishers.com

is

he rs

NOTICE TO THE READER The Publisher has taken reasonable care in the preparation of this book, but makes no expressed or implied warranty of any kind and assumes no responsibility for any errors or omissions. No liability is assumed for incidental or consequential damages in connection with or arising out of information contained in this book. The Publisher shall not be liable for any special, consequential, or exemplary damages resulting, in whole or in part, from the readers’ use of, or reliance upon, this material. Any parts of this book based on government reports are so indicated and copyright is claimed for those parts to the extent applicable to compilations of such works.

Pu bl

Independent verification should be sought for any data, advice or recommendations contained in this book. In addition, no responsibility is assumed by the publisher for any injury and/or damage to persons or property arising from any methods, products, instructions, ideas or otherwise contained in this publication.

en

ce

This publication is designed to provide accurate and authoritative information with regard to the subject matter covered herein. It is sold with the clear understanding that the Publisher is not engaged in rendering legal or any other professional services. If legal or any other expert assistance is required, the services of a competent person should be sought. FROM A DECLARATION OF PARTICIPANTS JOINTLY ADOPTED BY A COMMITTEE OF THE AMERICAN BAR ASSOCIATION AND A COMMITTEE OF PUBLISHERS. Additional color graphics may be available in the e-book version of this book.

Sc i

Library of Congress Cataloging-in-Publication Data

N

ov

a

ISBN: 978-1-62417-742-2

Published by Nova Science Publishers, Inc. † New York

,I

nc . Preface

he rs

CONTENTS

vii

Chapter 3

Chapter 4

is

27

Development of a New Curvature Law of the Wall for Internal Swirling Axial Flows Xiuhua A. Si and Jinxiang Xi

71

On Noise Prediction from a Compact Region of Turbulence in an Infinite Circular Hard-Walled Duct A. O. Borisyuk

93

Computer Simulation of Turbulent Flow Generated by a Deformed Anchor Impeller S. Karray, Z. Driss, H. Kchaou and M. S. Abid

en

Chapter 5

Turbulent Flow Structures for Different Roughness Conditions of Channel Walls: Results of Experimental Investigation in Laboratory Flumes D. Termini

Sc i

Chapter 6

a

Chapter 7

N

ov

Chapter 8

1

Large-Eddy Simulation of Turbulence-induced Aero-Optical Effects in Free Shear and Wall Bounded Flows K. Volkov

Pu bl

Chapter 2

On Analytical Solutions to the Three-Dimensional Incompressible Navier-Stokes Equations with General Forcing Functions and Their Relation to Turbulence Gunawan Nugroho

ce

Chapter 1

Study of Turbulent Flow on an Open Circuit Wind Tunnel Z. Driss and M. S. Abid Color Doppler Ultrasound (C. D. U. S.) Analysis of Turbulent Flow in Three Different Animal Models Jorge Elias Jr., Karina M. Mata, Cleverson R. Fernandes and Simone G. Ramos

119

147 161

185

vi

Contents

Chapter 10

Entrainment of Coarse Solid Particles by Energetic Turbulent Flow Events M. Valyrakis Elastic Effects on the Inviscid Instability of Shear Flows Ahmed Kaffel

219 251

N

ov

a

Sc i

en

ce

Pu bl

is

he rs

,I

Index

201

nc .

Chapter 9

nc . ,I he rs

PREFACE

"TURBULENT FLOWS: PREDICTION, MODELING AND ANALYSIS"

N

ov

a

Sc i

en

ce

Pu bl

is

Turbulent flow means fluid flow in which the fluid undergoes irregular fluctuations. Understanding the turbulent behavior in flowing fluids is one of the most intriguing and important problems that have been the focus of research for decades due to its great importance in a variety of engineering applications. Common examples of turbulent flow include atmospheric and ocean currents, flow through turbines and pumps, blood flow in arteries, oil transport in pipelines, lava flow, and the flow in boat wakes and around aircraft wing tips, etc. Moreover, studying and understanding of turbulence is necessary to comprehend the flow of blood in the heart, especially in the left ventricle, where the movement is particularly swift. A solid grasp of turbulence, can also reduce the aerodynamic drag on an automobile or a commercial airliner, increase the maneuverability of a jet fighter or improve the fuel efficiency of an engine. In the present book, we focus on areas of current turbulence research. Recent progress on modeling and analysis of turbulent flows is reviewed, and likely directions for future research on these topics are indicated. This text is unusual in as much as it provides both general commentaries as well as recent specialized developments in the field of turbulence modeling. As such it provides access to both historically validated and accepted results and newer ideas and approaches to the problem of modeling turbulence. Also of interest, we chose to devote a fraction of this book to recent developments in stability of parallel shear flows which will be useful to study eventual transition to turbulence. This book is unique in that it provides a balanced perspective with an emphasis on both rigorous mathematical developments and a focus on engineering problems in a field traditionally dominated by empiricism. We sincerely hope that the reader will find our choice of material interesting and stimulating. The book is foremost intended for researchers and graduate students with a basic knowledge of fundamental fluid dynamics. We particularly hope it will help young researchers at the beginning of their scientific careers to gain an overview as well as detailed knowledge of the recent developments in the field. Many colleagues have generously provided comments and material from their past and current research. We especially wish to thank all of the authors who submitted chapters at our requests. Editors

a

ov

N ce

en

Sc i he rs

is

Pu bl

,I

nc .

nc .

In: Turbulent Flows: Prediction, Modeling and Analysis ISBN: 978-1-62417-742-2 Editor: Zied Driss © 2013 Nova Science Publishers, Inc.

,I

Chapter 1

Pu bl

is

he rs

ON ANALYTICAL SOLUTIONS TO THE THREE-DIMENSIONAL INCOMPRESSIBLE NAVIER-STOKES EQUATIONS WITH GENERAL FORCING FUNCTIONS AND THEIR RELATION TO TURBULENCE Gunawan Nugroho

Department of Engineering Physics, Institut Teknologi Sepuluh Nopember Jl. Arief Rahman Hakim, Sukolilo, Surabaya, Indonesia

ce

ABSTRACT

en

The three-dimensional incompressible viscous flows are solved in this work with the application of the transformed coordinate which defines as a set of functionals, hi    ki x  li y  mi z  ci t . The solution is proposed from the base of general Riccati

a

Sc i

equation, which is firstly substituted into the Navier-Stokes equations to produce the polynomial equation with variable coefficients. The resultant solutions from the system of Riccati and polynomial are then evaluated by the proposed method of integral evaluation. The simulation shows the fluctuation of the decaying velocity through time due to energy dissipation. The rapid velocity accumulations are also detected providing that the solutions may produce singularity in the small scale of turbulent flows. The bifurcation is then detected which revealed the strong nonlinearity in the small scale of turbulent flows.

N

ov

Keywords: continuity equation, the navier-stokes equations, analytical solutions, partial differential equations, integral evaluation

MSC Number (2010): 35A09, 35C05, 76D05



Email address: [email protected], [email protected]

2

Gunawan Nugroho

1. INTRODUCTION

N

ov

a

Sc i

en

ce

Pu bl

is

he rs

,I

nc .

The problem of searching for the classes of analytical solutions of the full Navier-Stokes equations is highly demanding from both theoretical and practical viewpoints, as has been described in the literature [1]. Analytical solutions often play a special role in the theory of nonlinear equations and they are found able to describe the detailed behaviour of the concerning systems [2]. Analytical solutions also facilitate a theoretical understanding, paving the way to global solutions. They may help explain the issue of global smoothness in time [3]. Moreover, the solutions may be examined as models for turbulence [4]. Also, some particular solutions such as vortex solutions play a significant role in the development of turbulence theories [5]. The main difficulty of analytical solution of the Navier-Stokes equations is the contribution of the nonlinear terms representing fluid inertia which then troubled the conventional analysis in general cases. However, it is not surprising that they attract many efforts until recently, since the full solution of the three-dimensional Navier-Stokes equations still remains as one of the open problems in mathematical physics, especially in the sense of functional analysis [6]. The most important developments can be found in the review papers or in the monograph as in [7,8]. As one of the most interesting is a numerically based existence theorem which then can be established the properties of the solution [9]. By taking that numerical solutions converge to analytical solutions, the similar procedure can then be applied. Concerning to the analytical solutions, there are some works have already been conducted in the literatures [10,11,12]. As in the most cases, analytical solutions are examined only in special conditions in which the nonlinearity are weakened or even removed from the analysis. The type of the simplified analysis is applicable for steady and unsteady Couette and Poiseuille flows in which the nonlinear terms are removed permanently [13]. The other less known example is applicable to Beltrami flows in which the nonlinear terms are nonzero in the Navier-Stokes equations but they fade in the vorticity equations [14]. However, more sophisticated analysis of the Navier-Stokes equations is also conducted and gives more insight to the problems. One of them is the transformation of the NavierStokes equations to the Schrodinger equation, performed by application of the Riccati equation [15]. It has good prospects since the Schrodinger equation is linear and has well defined solutions. The method of Lie group theory is also applied in order to transform the original partial differential equations into ordinary differential systems [16]. It is concluded that an approximate series solution is obtained. The same route is taken by Meleshko [17] and by Thailert [18], in transforming the Navier-Stokes equations to solvable linear systems. Furthermore, less popular methods, such as the Hodograph-Legendre transformation, have also been applied to reduce the original problem to one more tractable, and thus closer to the goal of obtaining analytical solutions [19]. The reduction of the full set of Navier-Stokes equations to be a class of nonlinear ordinary differential equation is also performed [20]. The solution applied to both zero and constant pressure gradient cases. The method of introducing special solutions for velocity has also been investigated in [21,22]. The connection of the solutions with turbulence is discussed intensively by investigating that global in time continuation is still unsolvable for three-dimensional flows [25,26].

On Analytical Solutions to the Three-Dimensional Incompressible Navier-Stokes ...

3

ce

Pu bl

is

he rs

,I

nc .

The question arises as if such singularities exist, they might be related to turbulence by invoking that we have global smooth solution for two-dimensional flows, and turbulence is three-dimensional phenomena. This hypothesis, however, has serious difficulties as the observed phenomena are so far bounded in nature. In that case, different viewpoint is proposed. The argument that turbulent solutions should have no singularities is supported based on the triviality of the solution for simple energy equation. Analysis of global trivial solutions is important from mathematical and physical aspects, it has wide application and may be correlated to many areas where some hierarchical solutions are arranged. Classical procedure of vector identities is implemented for producing trivial solutions. Violation from trivial solutions is also investigated. Investigation of nontrivial solutions is related to the rate of energy generated or destroyed. The nontrivial solutions from the Navier-Stokes equations are then established as turbulent solutions due to energy accumulations. Apart from the simplified analysis, the current research method is also inspired by the importance of the Burger equation which is able to describe one-dimensional turbulence. The Burger equation can be viewed as a Riccati equation in one-dimensional case which the solution is mathematically well-posed and physically relevant for showing the collection of shocks, decaying velocity and its profile in nozzle flows even though it is not very accurate for three-dimensional cases. However, the equation provided us with the mathematical insight and physical understanding on the nature of intermittency of turbulence. In this work, the mathematical formulations by implementing the Riccati equation are developed more. The procedure for obtaining the solutions are derived in a quite simple way which is based on the implementation of coordinate transformation of a set of traveling wave ansatzs. The resulting expression is then produced polynomial equation which then combined with the solution of the Riccati to form final solutions.

en

2. PROBLEM FORMULATION

Sc i

The problem of searching the analytical solutions to the Navier-Stokes equations is a challenging task. The equations are usually simplified and solved concerning to the specific problems. However, our contribution is to evaluate the scheme for general incompressible flow cases in cartesian coordinate. Consider the continuity and the three-dimensional incompressible Navier-Stokes equations with external forces as in the following form, (1a)

u u u u 1 p  2u  2u  2u u v w    2   2   2  F1 t x y z  x x y z

(1b)

v v v v 1 p  2v  2v  2v u v  w     2   2   2  F2 t x y z  y x y z

(1c)

N

ov

a

u v w   0 x y z

Gunawan Nugroho w w w w 1 p 2 w 2 w 2w u v w    2   2   2  F3 t x y z  z x y z

where



is static pressure,

p

is fluid density, 

(1d)

nc .

4

is kinematic viscosity and

    x ,  y ,  z  . All solutions describe the three velocity components in the three

spatial directions, i.e., u  u  x, y, z, t  , v  v  x, y, z, t  and w  w  x, y, z, t  .

u

 v



 wz dx  K1  y, z , t 

y

(2)

is

x

he rs

,I

The presence of forcing functions Fi is physically relevant and can be widely observed in atmospheric flows, which the flow under the presence of the coriolis and electromagnetic forces. Beginning the process, the following relation is produced from the continuity equation (1a),

Pu bl

2.1. The Formulation for v Velocity and the Pressure Relation Equation (2) can be substituted into (1d) to give, wt    

 v x

y







1  wz dx  K1  wx  vwy  wwz   pz   wxx  wyy  wzz  F3  

(3a)

p

 F    w 3

xx



 wyy  wzz  wx  

 v x

y

en

z

ce

Then the pressure relation can be determined as,





 wz dx  K1   wt  vwy  wwz dz  K 2  x, y, t  

(3b)

The next step is implementing equation (2) and (3b) into (1c),

 v x

  y 

y





 z





 wz dx  K1vx  vt  vv y  wvz   vxx  v yy  vzz  K 2 y

Sc i

vx





F3   wxx  wyy  wzz  wx  

 x



 v y  wz dx  K1   wt  vwy  wwz dz   



(4a)

Suppose that the coordinate also satisfy the following set of traveling wave ansatz,

N

ov

a

hi    ki x  li y  mi z  ci t

(4b)

where ki , li , mi and ci are constants. The step is now transforming the cartesian coordinate into  -coordinate by setting the

subscript i  1, 2,3, 4 such that we have four equations for the coordinate transformation. The results of x  1   , y   2   , z   3   and t   4   are then determined. Applying one

On Analytical Solutions to the Three-Dimensional Incompressible Navier-Stokes ...

5

of the traveling wave ansatzs, i.e. h1    k1 x  l1 y  m1 z  c1t into (4a), the considered equation then transformed into,

,I

nc .

 k 2 l 2 m2  k1 k c l m l l v  klv  kmw   1 K1v   1 v  1 vv  1 wv    12  12  21  v  1 F3  F2  1 K 2  h1 h1 h1  h1 h1 h1 h1 h1 m1 h1   2 3 2 2  k l  l  cl  l lm  k l2 l   21 1  21  12 1  w  1 w  k1 1 v  k1l1w   K1   1 1 w  1 vw  1 ww    h1 m1 h1 m1 h1 m1  h1 m h m h m h 1 1 1 1    1 1  

Rearranging (4c) to give, a1v  a2 vv   a3  a4 w  v  a5 w v  a6 w  a7 w  a8 ww 

he rs

(4c)

l1 l F3  F2  1 K 2 m1 h1

(4d)

Pu bl

is

where ai are variables that depend on  . Lemma 1: Equation (4d) is transformable into the system of Riccati and third order polynomial equations. Proof: Consider the Riccati equation, v  b1v 2  b2 v  b3

(5a)

such that the following relation is satisfied, v  b1 v 2  2b1vv  b2 v  b2 v  b3



ce







v  2b12v3  b1  3b1b2 v2  b2  2b1b3  b22 v  b3  b2b3

(5b)

en

Substituting into (4d) to get the following polynomial equation,

 2a b

2 1 1







 a2b1 v3  a1b2  2a1b1b3  a1b22  a2b3  a3b2  a5 w  a4b2 w v 

Sc i

 a1b1  3a1b1b2  a2b2  a3b1  a4b1w v 2  a6 w  a7 w  a8 ww  a4b3 w  a1b3  a1b2b3  a3b3 

(5c) l1 l F3  F2  1 K 2 m1 h1

N

ov

a

Therefore, equation (4d) is transformed to the system of (5a) and (5c). This proves lemma 1.





Set 2a1b12  a2b1  0 and  a1b1  3a1b1b2  a2b2  a3b1  a4b1w  0 , thus the expression

for b1 and b2 are defined by,

6

Gunawan Nugroho

b2  

a2 2a1

(5d)

 a1b1  a3b1  a4b1w

(5e)

3a1b1  a2

a6 w  a7 w  a8 ww  a4b3 w  a1b3  a1b2b3  a3b3 

l1 l F3  F2  1 K 2 m1 h1

(5f)

he rs

v

,I

Therefore, equation (5c) is reduced into,

nc .

b1  

a1b2  2a1b1b3  a1b22  a2b3  a3b2  a5 w  a4b2 w

is

The expression for b3 will be determined later as a requirement of unique solutions under general initial-boundary values. Also note that the solution for the system (5a) and (5f) will be similar to that of the velocity in z direction as derived in the subsequent paragraph.

Pu bl

2.2. The Formulation for w Velocity

The step now is performing equations (2) and (3b) into (1b),   t 

 v x

y



 wz dx  K1      

 v x

y

  wz dx  K1  v y  wz  v   y 





 



 v x

y



 wz dx  K1  



 F   w  w  w    xx yy zz    3   (6a)  v y  wz dx  K1     dz     x  x  z  w v y  wz dx  K1  wt  vwy  wwz   x       x 2 2      K 2 x    vxy  wxz  2  v y  wz dx  K1   2  v y  wz dx  K1    F1       z  x   y  x 

















en





ce

 w  z 

Performing the coordinate transformation (4b),  l  c1 m l k1l1v  k1m1w  K1     k1l1v  k1m1w   K1   1 v  1 w   v 1 k1l1v  k1m1w  K1    h1 h1 h1  h1  kl  m k km l2 w 1 k1l1v  k1m1w  K1  1 K 2    121 v  1 2 1 w    12 k1l1v  k1m1w  K1   h1  h1 h1 h1 h1     k3  k 2 l  kc  k1  k l2 km  k kl k 1  21 1  1 2 1  w  1 w  1 1 v  k12 w   K1   1 1 w  1 1 vw  1 ww    F3    2    h1 m1 h1 m1 h1  m1  h1 h1 m1 h1    m1  h1 m1   









N

ov

a

Sc i





m12 h12









 k1l1v  k1m1w  K1   F1 (6b)

Thus by grouping the above equation, the momentum in z direction is given by,

On Analytical Solutions to the Three-Dimensional Incompressible Navier-Stokes ...

 c k1 k l2 m2  F3  F1  1 K 2   1   12   21  K1   h1 m1 h1 h1 h1  (6c)   a13vv   a14  a15 w  v  a16 w v

7

a12 v

nc .

a9 w  a10 ww  a11w 

Since equation (6c) is similar to that of (4d), it is reasonable to state that the solution for

,I

v is also similar to (5f).

he rs

Repeating the procedure described by (5a – f) with different variables, b4 , b5 and b6 yielding,  c k1 k l2 m2  F3  F1  1 K 2   1   12   21  K1  h1 m1 h1 h1 h1   2 a12b5  2a12b4b6  a12b5  a13b6  a14 b5  a16 w  a15b5 w

a9 w  a10 ww  a11w  a15b6 w  a12b6  a12b5b6  a14b6  v

(6d)

is

Equating (5f) and (6d) to get,

l1 l F3  F2  1 K 2 m1 h1

Pu bl

a6 w  a7 w  a8 ww  a4b3 w  a1b3  a1b2b3  a3b3  v

a1b2  2a1b1b3  a1b22  a2b3  a3b2  a5 w  a4b2 w  c k1 k l2 m2 F3  F1  1 K 2   1   12   21  h1 m1 h1 h1 h1  a12b5  2a12b4b6  a12b52  a13b6  a14b5  a16 w  a15b5 w

a9 w  a10 ww  a11w  a15b6 w  a12b6  a12b5b6  a14b6 

(6e)

ce



  K1  

en

which then can be performed algebraically to form a single expression.

3. THE ANALYTICAL SOLUTIONS

Sc i

In this research, we are interested to the class of solutions that is physically important [23,24], ih  w  r   e 1  

(7a)

N

ov

a

Substituting (7a) into (6e) to get,



 

r a17 r  a18 rr  a19 r  a20 r  a21  r a22 r  a23rr  a24 r  a25 r  a26  a27 r  a28 rr  a29 r  a30 r  a31  0 By setting,



(7b)

8

Gunawan Nugroho r  b7 r 2  b8 r  b9

(7c)









r  2b72 r 3  b7  3b7b8 r 2  b8  2b7b9  b82 r  b9  b8b9

(7d)

nc .

such that the following relation is fulfilled,

,I

Thus, applying (7c) and (7d) into (7b), the resulting polynomial equation is defined by, a32 r 5  a33r 4  a34 r 3  a35 r 2  a36 r  a37  0

where a37 is given by,









he rs

(7e)

a37  a17 b9 b9  b8b92  a19b92  a21b9  a27 b9  b8b9  a29b9  a31

(7f)

Pu bl

is

Set a37  0 then b8 can be obtain as a function of b9 . This procedure will reduced equation (7e) to the fourth order polynomial which the root is solvable by radicals, say be written as r     . Meanwhile, b7 will be defined later by the similar condition applied to b3 .

3.1. The Solution for the Riccati Equation

In order to solve the system of (7c) and (7e), the following step is necessary [25] f 2 f2

to generate,

b7 Z 2  f1Z  C1 f 2b9 C1 f 2

en

Z 

ce

Lemma 2: Consider equation (7c) and set b8  f1 

where Z  C1 f2 r . There exists a function  such that f2b9   Z , which generate the

Sc i

Bernoulli equation. The Riccati equation then has a closed-form exact solution when  is solvable.

N

ov

a

Proof: Set b8  f1 

r 

f 2 f2

to rearrange equation (7c) as,

f2 x 1 r  C1 f 2 r   b7 r 2  f1r  b9 f2 C1 f 2

(8a)

where C1 is a constant. Suppose that Z  C1 f2 r , then the following equation is produced,

On Analytical Solutions to the Three-Dimensional Incompressible Navier-Stokes ... Z 

b7 Z 2  f1Z  C1 f 2b9 C1 f 2

9

(8b)

b7 Z 2   f1  C1  Z C1 f 2

(8c)

,I

Z 

The solution of (8c) is expressed as,



  f1 C1 d

(8d)

b7   f1 C1 d e d  C1 f 2



b7   f1 C1  d   f1 C1 d  C1 f 2 D . e d   D , then e    C1 f 2 b7



Pu bl

Solving for D ,

 b7   f1 C1 d e d  e   C1 f 2



Set  

b7 b9



d

(8e)

A and differentiate (8e) once to give, B 

A

ce

D

he rs

Let



e

is

Z

C1 f 2b9

nc .

Set f2b9   Z , the original problem is transformed into,

b7 b9 B d A

g

en

 A   f1  C1 B  d e  b9 Be  C1 f 2

(8f)

Thus, A and B are given by, 

 f1d

Sc i

A

f 2 ge

C1



1  f 2 ge  B 

f1d 

and B  

d

b9



g 1 gb7 d A

(8g)

N

ov

a

The function  can be determined as, A    B



f 2b9 e C1



 f1d

 A gb d 1

7

(8h)

  f1d  1 f 2 ge  d B



Without loss of generality suppose that, f 2 ge

 f1d  A to produce,

10

Gunawan Nugroho

A  B

 f1d

b7  f1d d  e   f2 A d  C1  B

 

(8i)

nc .



f 2b9 e

  f1d    f 2 b9 e   





1 2  f1d   b7  e  d  d    f2    

(8j)

he rs

 A C2      B    

,I

where  will be defined later. The solution for  is then,

is

The above equation is combined with (8d) to form the solution of the general Riccati equation. This proves lemma 2.

3.2. The Complete Solution

Pu bl

The system which is represented by (7c) and (7e) will have a simultaneous solution as it is driven by equation (8) and equation (7e). The claim is concluded in the following theorem, Theorem 1: Consider the solution of the general Riccati equation as described by (8d) and (8j). By combining with the root of (7e), r     , then the expressions of f1 and f 2 can be determined. The resulting expressions thus complete the solution of the system defined by (7c) and (7e).

ce

Proof of theorem 1: Equating the results from (7e) and lemma 2 as follows, b9

     

  f1d    f 2 b9 e   



en

  or b     C2 9  



1 f1d   2  b7   e d  d     f2    

(9a)

Sc i

Integrate the above equation and perform the algebraic calculations to give,

a



with  

ov

N



f1d

b7  f1d d  C e  e  3 f2 f 2b9  C3 b9

  



b9

 d   

  



b9

 d   

2

  e f2  



f1d 

2

 .  

Differentiate (9b) once to get the relations of f1 and f 2 as in the following,

(9b)

On Analytical Solutions to the Three-Dimensional Incompressible Navier-Stokes ...   1  1  b7    f1       f2    f 2  

d

(9d)

Equating the above equation with b8  f1 

f 2 f2

,I



to get the following expression,

he rs

b7

 b   b8  f1   7   f1     

(9e)

  



b9

 d   

2

 b8 d    e   

 b8d d

Pu bl

C3 b9

is

which has also solved  as in the following,



nc .

(9c)

The solution for f 2 is then,

 f1d e  f2  e 

11

 b e 7

(9f)

It is important to mention that the solution for f1 does not exist. This condition is



b9

 d    

 b e 7

 b8 d d 

(9g)

en

2C3  b8 d  e   

ce

consistent with equation (9d), since f1 will vanish when (9d) is substituted into (8j). The above equation can be rearranged as,

Equation (9g) forms the linear first order equation in  , the solution is then expressed

Sc i

as,

N

ov

a

  2C3  2b8C3           b9      b7  2C3      b9

    exp  

  2C3  2b8C3           b    7  2C3      b9



    d  

Therefore, the explicit solution of the system (7c) and (7e) is given by,

(9h)

Gunawan Nugroho

r

b9



 b7

  f1 C1 d e

 f2



b7 e f2

  f1 C1 d

e

 d



 



  C1  d  

(9i)

 b7   C1  d      d

b7 

  e

nc .

12

C2 

   

     

  f1d    f 2 b9 e   



1





 1  f1d   2  b7 C2     d  d    e   f2        



2 b9 d   

  

b7  d   b9 e    



b7

  e



b7



d

1  2   d  d       .   

he rs

C  2 

,I

where the function  is determined by,

(9j)

w  w1  w2  ........ v  v1  v2  ..........

Pu bl

is

This completes the proof of theorem 1. Note that the procedure that explained by equation (8) and equation (9) can also be applied to system (5a) and (5f). Therefore, by substituting (9i) into (7a) the solution of mass and momentum equations, w, v, p and u are defined by (7a), (5g) or (6d), (3b) and (2). Therefore, based on the obtained solutions, the result then can be generalised as in the following,

(9k)

ce

u  u1  u2  ..........

p  p1  p2  .........

en

which can also be solved by the same procedure.

Sc i

4. INTEGRAL EVALUATION OF ANALYTICAL SOLUTIONS

a

It is important to note that the integrals which appear in the analytical solutions are usually approximated in series forms [26], by which the solution is then no longer exact. In order to resolve the problem, now the following integral is considered [27],

N

ov

A

 e

 fd d 

(10a)

By setting, A





e 

fd 

d   C3  R  Q  e

 gd

(10b)

On Analytical Solutions to the Three-Dimensional Incompressible Navier-Stokes ...

13

  fd e

C3





 R  Q  e

 gd  R  Q  e  gd  R  Q  ge  gd     

Rearranging the above equation as,          f  gd  R    g  R  e  Q    g  Q  C3      

1





e

 gd  e gd   e f  gd   C3



,I

(10c)

he rs

The solution of R is then expressed by, R

nc .

where C3 is a constant, equation (10b) can be differentiated once to give,

        Q    g  Q  d       

is

Let,

Pu bl

      f  gd  e  Q    g  Q   f3 C3    

(10d)

(10e)

Then, R is evaluated in the following, 1





e

 gd   



  gd f3 d  e   

  

 

  gd  f3 d  ge  d    

(10f)

ce

R

Suppose that from equation (10e),

en

  f  gd e  C4 C3

Sc i

where C4 is also a constant. The expression for e



 fd  e  gd

(11a)

a

C3C4

e

 gd  is written as,

N

ov

Thus, equation (10f) will become,

R

  gd  1 e  C3C4

   



   fd f3 d   e   

  

 

    fd f3 d    e    

    d    

(11b)

14

Gunawan Nugroho

Without loss of generality, set

   fd   , and the expression of f3 is 

 f  d  ln   e 3

1     fd ln  e    

     

(11c)

,I

f3 

nc .

obtained as,

Q

1





e

 gd

 C

4

 f3  e

he rs

The solution for Q is consequently obtained from (10e) as in the following relation,

 gd d

Substituting (11a) to get,

 C

4

 f3   e

 fd d 

is

  gd  1 e  C3C4

Pu bl

Q

(11d)

Equations (10b), (11b) and (11d) will give the evaluation as,





e 

fd 

d   C3  R  Q  e

 gd  1  e  fd  1 C4 

C4

  C4  f3  e



fd 

d

(11e)

 fd d 

1    fd  e 3 

 f

en

 e

ce

where f3 is determined by (11c). Equation (11e) can be differentiated once and rearranged to be,   d  

(12a)

1     fd     fd n e  L and f3  ln  e     L , with n is an arbitrary       

Sc i Now suppose that

constant. The relation of  is then given by,

N

ov

a

 ne

n

 fd1n        f    

 

(12b)



1

Let       1 n , equation (12b) will then produce,    n  1  n e

n

 fd  2  1  n    f   2 n en  fd     1  n    f    n  1  n en  fd  2                



 

 



(12c)

On Analytical Solutions to the Three-Dimensional Incompressible Navier-Stokes ...

15

The solution for  is similar to that of formulation described by equation (8),  

n  fd    f     n  1  n e   2      

1 2      n fd      1  n     f     n  1  n e   2                   d    n n  fd n n  fd  d fd  2   e  d   n 1  n 1  fd 2   e     e  e  e e d          



,I

     C5    



he rs



nc .

1  n  

(12d)

2   n fd       n  1 3    f  x   1    f    n e    4    2  

  

1

2 n  fd     ne    d     

Pu bl

  C6    

is

 n  fd   1   Without loss of rigor, set      f    n e  to produce a more simplified 2   solution,



(12e)

 fd  e   n

n

fd 

d 

1 1 n 1    fd L   e 1 n 1  n  

en





e 

ce

The step is now performing the integration of (12a) to give, 1 n

  

(13a)

1

Substitute the relation       1 n into (13a) and rearranging the result to produce,

Sc i

1   2          n fd  n fd     C7     1    f    n e    n e     d           2              2   3   1       n  1 4    f  x  2     f        

N

ov

a



Rearrange and differentiate the above equation once to get,





e 

fd 

d 

 1 n



e 

fd 

16

Gunawan Nugroho   fd        e    fd  2 fd     2   1           e 1  e  D  C8     f       f   2   fd   2   fd   2        e  fd d          e d     e d              





nc .



,I

(13b) where D is expressed by 2

 n n  fd    e  

he rs

 n  fd   1     3   D   n  1   f      f   n e  2  4     



   fd   e fd     e  d            fd      e  d      2   2 fd       fd   1    f    2      e  d          e        2   



Pu bl

  

3

   fd   e  d  D  C8  E  F    





is

It is not hard to see that equation (13b) will produce a polynomial equation as in the following,



 1   and F     f  . Thus, the integral of 2    defined by the root of equation (13c). By taking one root of (13c) as,

en 

e 

fd 

d 

  fd   e  

C8 2 3       n  1 3    f     1    2  4        1     1        f      f     2     2     

a ov

N

  G  H 3  G 2 



1



1 2

3    G  H 3  G 2  

where G and H are defined by,



  e

 fd d  can be

  1      1     fd  f       f     e            2     2  

Sc i



fd 

ce



With E   e 

(13c)

 n  fd   f   n e  

  

 n n  fd    e           

1



1 2

3  

(13d)

On Analytical Solutions to the Three-Dimensional Incompressible Navier-Stokes ...

17

 fd    e    

  1      1     fd   f       f     e    2     2         

and

 fd    e   

  fd  e    

2  fd    1     f     2 e         2    

ce

C82 2  2  n fd        n  1 3    f     1    f    n e    4    2  

a ov

N

 f  d d 

    

2

 n n  fd 1    1            f       f      e 2       2       

Theorem 2: Consider the following integral equation,

    e

3

 n n  fd  1     1      f      f     e           2       2  

which then solves the integral in (10a). Therefore, the following theorem is just proved,

A

3

(13e)

  1      1     fd     f       f     e      2       2           

en

 fd    e   

     

 n n  fd 1    1            f       f      e 2       2       

C8 2 3   n fd        n  1 3    f     1    f    n e  2   4      

Sc i



 n n  fd 1    1            f       f      e 2       2       

C83 27  2  n fd        n  1 3    f     1    f    n e    4    2  

H

he rs

C8 2 2   n fd        n  1 3    f     1    f    n e  2   4      

  

is



2  fd    2 e    

Pu bl



  fd   e  

,I

nc .

    fd   1    f         e          2         fd   fd   2  fd              1      e    e           f     2 e       2         fd          1                 f    e          2      C82  G  2 6  2 n fd   n fd                  n  1 3    f     1    f    n e    n e     1     f   1     f            4 2  2  2                         

2

(13f)

18

Gunawan Nugroho There exists a functional      , with  and  are given by,  n  fd   1    f   n e  , 2  



 ne

n

 fd 

Such that the integral A can be evaluated as,

 fd d  C5

  1      1     fd  f       f     e             2     2  

  

 n n  fd    e           

is

2 3   n  fd   1       3  n  1   f      f   n e   4     2    1     1        f      f      2     2     

  G  H 3  G 2 



3    G  H 3  G 2  



1



1 2

3  

ce



1

1 2

Pu bl

 e

  fd   e  

2  n n  fd     d    e     

,I

1

2   n fd       n  1 3    f  x   1    f    n e    4    2  

he rs

  C4    

nc .

  

en

where C4 , C5 and n are arbitrary constants and G and H are defined by (13e) and (13f).

5. THE ONSET OF TURBULENCE

Sc i

Let  be the region of interest. It is supposed that the associated problem is a connected, bounded region in three-dimensional domain. Moreover, let i , i  1,..., m , be sub regions characterised by simple boundaries and ai , i  1,..., n be the sub regions where boundaries

a

are not simple. This means that the considered boundary i and ai are defined as regular and irregular surfaces respectively. Thus, the boundary-value problem is determined as follows; given density   0 such

N

ov

that velocity U is real vector field consists of  u, v, w components and p is real scalar field defined in every region i , i  1,..., m , and ai , i  1,..., n which fulfill the continuity and incompressible Navier-Stokes equations (1) with boundary and initial conditions U  U  x, y, z,0 on 1 or a1 . The problem also denotes n is the unit vector normal to the surfaces S parallel to the velocity. It will be demonstrated that if the nonzero divergence

On Analytical Solutions to the Three-Dimensional Incompressible Navier-Stokes ...

19

he rs

,I

nc .

condition in (1a) due to the production term in the control volume is implemented, the boundary value problem does not have any trivial solutions. It is very important consequence since the nonzero divergence in continuity equation has wide applications, for example, combustion problems. It is noted that the previous boundaries are characterized by the Cartesian coordinates and only applied to few simple applications. For this reason it is important to prove that the associated problem can be extended to general coordinates for non-simple boundaries. The general coordinates are also defined in i , i  1,..., m and ai , i  1,..., n such that general solution is mapped to the general coordinate in the same regions. The mathematical formulations of the associated general coordinates are,     x, y , z , t      x, y , z , t 

(14)

    x, y , z , t 

is

   t 

Pu bl

such that  ,  ,  and  are continuous in i , i  1,..., m and ai , i  1,..., n . It is also important to note that the boundary conditions considered in (1) are also applied to the general coordinates.

5.1. The Non Existence of Trivial Solutions

ce

In this section, some important properties of the solution of (1) are investigated in order to prove the triviality of the solutions with respect to the boundary values. The energy rate is defined as a product of static pressure p and flow rate U across the control surface. The

en

direction of velocity is parallel to the unit normal control surface n . It is supposed that the region  be the associated problem and consists of i parts of region with simple and non simple boundaries. Lemma 3: There exist p and U as solutions to the continuity and incompressible Navier-

Sc i

Stokes equations such that the rate of energy in the whole domain of i , i  1,..., m and ai , i  1,..., n with respect to the initial and boundary conditions U  U  x, y, z,0 on 1 or

a1 , is equal to zero. The solutions are, p = constant and U  constant

N

ov

a

in  i and ai . Proof: The assumed zero rate energy can be written as, E  pAU  12 U 2 AU  0 t

(15a)

20

Gunawan Nugroho Consequently, the divergence theorem can be applied to the whole region of interest as,

S

1 2



U 2U  n dS 

    pU  

1 2



U 3 d   const

(15b)

nc .

  pU  n 

  pU  n  S

1 2



U 2U  n dS 

  p U  U p  U  U  U U dV 2

2

V

,I

By using vector identity,   fF   f  F  F f , hence, it is identified,

(15c)

he rs

The above equation is then applied into all regions. Summing up the simple and non simple regions to give,

   p.U  U .p  U .U  U .U d      p.U  U .p  U .U  U .U d   const m

n

2

i 1

2

i

2

i 1

ai

2

is

(15d)

1  is set to unity. Suppose that the rate of energy is zero everywhere, then, equation 2 (15d) is automatically satisfied. The first and third terms of energy rate are always zero according to continuity. Therefore, since p and U are nonzero, then p and U must be

Pu bl

with

zero. Thus, the following trivial solution is defined,

(15e)

ce

p = constant and U  constant

Sc i

en

and we are done. It is obvious that the above trivial solution will only valid in special cases. The violation of the condition described in lemma 3 is strongly related to the Navier-Stokes equations which can be correlated to the onset of turbulence as explained in next section. In this case, the analyticity is implemented in order to generalize solutions to the limit of or near to triviality in some sharing regions. Hence, another interesting result is produced, Proposition: Any nontrivial solutions are analytic in regions  i and ai .

Furthermore, the function f  x  is considered in the intersection region i  ai .

Suppose that f  x  is at least twice differentiable and converge to a certain constant value A

N

ov

a

in a certain location which is related to lemma 1. According to the property of analytic function, it is reasonable to assume that if f  x  is convergent at a location y which is close enough to x will generate a certain value A ' such that A '  A . Therefore, for p and U with associated boundary conditions, it is eligible to consider a

sub region  in i  ai such that x    y where   x i and   y ai .

On Analytical Solutions to the Three-Dimensional Incompressible Navier-Stokes ...

21

Then, according to the analytic property of p and U , it is reasonable to conclude that p and U in  are equal to B ' and C ' where B '  B and C '  C in i , i  1,..., m and

nc .

ai , i  1,..., n . Therewith, the following statement is produced,

Lemma 4: Given any nontrivial solutions such that p  0 and U  0, can be applied to  , where  is a boundary of  i and also become at least one sub region or part of ai ,

,I

where p and U are trivial. The corresponding solution can be interchanged at the boundary such that there exist

he rs

p  constant and U  constant

in  .

is

It is noted that equation (15e) satisfies continuity and not the incompressible NavierStokes equations because of the presence of the forcing functions. Therefore, the following non existence theorem of trivial solutions is just already proved. Theorem 3: Any solution for initial-boundary value problems of the continuity and incompressible Navier-Stokes equations that satisfies equation (1) is not trivial, i.e. p ≠

en

ce

Pu bl

constant and U ≠ constant. Note that the investigated problem is strictly obeys Navier-Stokes equations which the condition of p  constant and U  constant is a very special case, even without forcing. Suppose that energy rate is produced if the condition violated, it is plausible that the excess energy will be distributed to the whole domain, then the observed variables will also deviate from trivial conditions. Furthermore, the case considered here strictly admits continuity in the form of velocity divergence, meaning that if in some cases velocity divergence is not zero, triviality in p and U will be more difficult to obtain, therefore theorem 3 imply, Corollary: Any non trivial solutions for boundary value problems of the continuity and incompressible Navier-Stokes equations as a possible onset of turbulence.

6. SOLUTION EXAMPLES AND DISCUSSIONS

N

ov

a

Sc i

In this section, the coefficients of Riccati equation will be determined to produce the known solutions. The collections of analytical solutions to the simplified Navier-Stokes equations can be investigated in the literatures [28]. Many solutions of the Navier-Stokes equations are unstable and can hardly have physical meanings. However, the method which is proposed in this work may handle physically reasonable solutions. Figure 1 shows the fluctuation velocity through time with different time samplings in the small scale. It is also shown that the velocity is decaying in the small scale, this is due to the dissipation mechanism of the turbulent flows which convert energy from large scale into heat at the molecular level. The other aspect of the solution is depicted in figure 2, the velocity accumulates rapidly with the extended time. The it is clear from the integral evaluation that the solution become bounded if the coefficient appeared in the Riccati equation is also bounded. This striking features will be the subject of further investigation of the mathematical properties of the solutions.

22

Gunawan Nugroho 5

nc .

3

,I

2

1

0 0.01

0.015

0.02

0.025

0.03

0.035

Time (s)

5

0.045

0.05

Pu bl

4 Velocity (m/s)

0.04

is

6

he rs

Velocity (m/s)

4

3

2

ce

1

0 0.01

0.015

0.02

0.025

0.03

0.035

0.04

0.045

0.05

Time (s)

en

Figure 1. The fluctuations of velocity through time. 4e+023

Sc i

3.5e+023

3e+023

2.5e+023 2e+023

N

ov

a

1.5e+023 1e+023 5e+022 0 -5e+022 -1e+023

0

0.1

0.2

0.3 Time(s)

Figure 2. (continued).

0.4

0.5

On Analytical Solutions to the Three-Dimensional Incompressible Navier-Stokes ...

23

4e+023 3.5e+023

nc .

3e+023 2.5e+023 2e+023 1.5e+023 1e+023

,I

5e+022 0

-1e+023

0

0.1

0.2

he rs

-5e+022

0.3

0.4

Time(s)

Figure 2. The rapidly increasing velocity through time. 1.5

is

0.5

0

-0.5

-1

0

Pu bl

imaginary

1

0.5

1

2

3

4

5

ce

real 2

en

1.5

a

Sc i

imaginary

1

0.5

0

-0.5

-1

0

1

2

3

4

5

6

real

N

ov

Figure 3. The bifurcation of solutions.

Figure 3 finally shows the bifurcation phenomena which appeared in the solutions. The phenomena are detected by plotting the real and imaginary parts of the solutions, which explicitly shows the strong nonlinearity of the governing equations. This special feature is known for the incompressible Navier-Stokes equations by the qualitative analysis of dynamical systems which deserves further investigation in explicit solutions.

24

Gunawan Nugroho

he rs

,I

nc .

The above example shows that the physically simple and useful solutions can be produced within a specific initial-value if the coefficients are carefully chosen. However, the construction of analytical solution in this work may lead to the sensitivity with respect to initial condition because the method is usually performed forward. Moreover, it is found that u is forced in order to be unique and p is not unique which can yield the unstable and fluctuating solutions and leads to turbulent solutions. Therefore, it is noted that the method may also valid for turbulent flow problems. According to the proposition which was stated in the previous work that any non trivial solutions can be regarded as turbulent solution [29,30], one can investigate the nature of the solution to describe the structure of turbulent flows. It is well accepted that turbulence structures are represented by the condition of non zero

mean fluctuation [31]. Based on the decomposition scheme, w  w  w ' and (16a) it is not

is

hard to verify that the relation w  w is hold such that w  w  w ' . Therefore, w '  0 , which is also applied to v, u and p as well.

Pu bl

CONCLUSION

a

[2]

N

ov

[3]

[4] [5]

REFERENCES

Kao S.K., An Analytical Solution for Three-Dimensional Stationary Flows in the Atmospheric Boundary Layer over Terrain, Journal of Applied Meteorology, Vol. 20, 1980. Galaktionov V.A. and Svirshchevskii A.R., Exact Solutions and Invariant Subspaces of Nonlinear Partial Differential Equations in Mechanics and Physics, Taylor and Francis Group, Boca Raton, 2007. Zhou Y., On a Regularity Criterion in Terms of the Gradient Pressure for the NavierStokes Equations in R N , Z. angew. Math. Phys 57, 2006, pp. 384 – 392. Ohkitani K., A Survey on a Class of Exact Solutions of the Navier-Stokes Equations and a Model for Turbulence, Publ. RIMS, Kyoto Univ., 40, 2004, pp. 1267 – 1290. Lyberg M.D. and Tryggeson H., An Analytical Solution of the Navier-Stokes Equations for Internal Flows, J.Phys.A: Math. Theor. 40, No. 24, 2007, pp. 465 – 471.

Sc i

[1]

en

ce

The method for producing the analytical solutions of the three-dimensional incompressible Navier-Stokes equations is introduced in this chapter. The solution is constructed under the implementation of the Riccati equation which then produces polynomial equation. The solution that is generated from the system of Riccati and polynomial equations is then expanded by the utilization of the proposed integral evaluation. The results show the fluctuation of the decaying velocity through time due to energy dissipation. It also reveals the rapid velocity accumulations providing that the solutions may produce singularity in the small scale of turbulent flows. The bifurcation is then detected which shows the strong nonlinearity in the small scale of turbulent flows.

On Analytical Solutions to the Three-Dimensional Incompressible Navier-Stokes ...

[12] [13] [14]

[15] [16]

[17] [18] [19]

[20]

nc .

Sc i

[21]

,I

[11]

he rs

[10]

is

[9]

Pu bl

[8]

ce

[7]

Ohkitani K., A Miscellany of Basic Issues on Incompressible Fluid Equations, Nonlinearity 21, 2008, pp. 255 – 271. Galdi G.P., Heywood J.G. and Rannacher R., Contribution to Current Challenges in Mathematical Fluid Mechanics, Birkhauser Verlag, Basel, Switzerland, 2004. Foias C., Manley O., Rosa R. And Temam R., Navier-Stokes Equations and Turbulence, Cambridge University Press, UK, 2004. Heywood J.G., Nagata W. and Xie W., A Numerically Based Existence Theorem for the Navier-Stokes Equations, J. Math. Fluid Mech. 1, 1999, pp. 5 – 23. Wang C. Y., Exact Solutions of the Unsteady Navier-Stokes Equations, Appl.Mech. Rev .42, 1989, pp. 269 – 282. Wang C. Y., Exact Solutions of the Steady Navier-Stokes Equations, Annu. Rev.Fluid. Mech.23, 1991, pp. 159 – 177. Okamoto H., Exact Solutions of the Navier-Stokes Equations via Leray’s Scheme, Japan J. Indust. Appl. Math. 14, 1997, pp. 169 – 197. White F.M., Fluid Mechanics 5th edition, McGraw-Hill, 2003. Shapiro A., The Use of an Exact Solution of the Navier-Stokes Equations in a Validation Test of a Three-Dimensional Nonhydrostatic Numerical Model, Monthly Weather Review, Vol. 121, 1993. Christianto V. and Smarandache F., An Exact Mapping from Navier-Stokes Equation to Schrodinger Equation via Riccati Equation, Progress in Physics, Vol. 1, 2008. Kamran F., Zu-Chi C., Xiaoda j., Cheng Y., Similarity Reduction of a (3+1) NavierStokes System, Engineering Computations: International Journal for Computer-Aided Engineering and Software, Vol. 23, No. 6, 2006, pp. 632 – 643. Meleshko S.V., A Particular Class of Partially Invariant Solutions of the Navier-Stokes Equations, Nonlinear Dynamics 36, 2004, pp. 47 – 68. Thailert K., One Class of Regular Partially Invariant Solutions of the Navier-Stokes Equations, Nonlinear Dynamics, 2005. Mohyuddin M.R., Siddiqui A.M., Hayat T., Siddiqui J. and Asghar S., Exact Solutions of Time-Dependent Navier-Stokes Equations by Hodograph-Legendre Transformation Method, Tamsui Oxford Journal of Mathematical Sciences 24 (3), 2008, pp. 257 – 268. Nugroho G., Ali A.M.S., and Abdul Karim Z.A. On a Special Class of Analytical Solutions to the Three-Dimensional Incompressible Navier-Stokes Equations, Applied Mathematics Letters 22, 2009, pp. 1639 – 1644. Sidorov, A. F., On Two Classes of Solutions of Fluid and Gas Dynamics Equations and their Connection with Theory of Travelling Waves, Journal of Applied Mechanics and Technical Physics 30 (2), 1989, pp. 34 – 40. Rajagopal, K. R., A class of exact solutions to the Navier–Stokes equations, International Journal of Engineering and Science 22 (4), 1984, pp. 451 – 458. Moore R.L., Exact Non-Linear Forced Periodic Solutions of the Navier-Stokes Equation, Physica D. 52, 1991, pp. 179 – 190. Green J., Atmospheric Dynamics, Cambridge University Press, UK, 1999. Nugroho G., Reduction of Higher Order Linear Ordinary Differential Equations into the Second Order and Integral Evaluation of Exact Solutions, Tamsui Oxford Journal of Informatics and Mathematical Sciences, 2011 (Submitted). Wong R., Asymptotic Approximations of Integrals, Society for Industrial and Applied Mathematics, 2001.

en

[6]

[22]

a

[23]

N

ov

[24] [25]

[26]

25

26

Gunawan Nugroho

N

ov

a

Sc i

en

ce

Pu bl

is

he rs

,I

nc .

[27] Nugroho G. Ali A.M.S. and Abdul Karim Z.A., The Navier-Stokes Equations: On Special Classes and Properties of Exact Solutions of 3D Incompressible Flows, Lambert Academic Publisher GmbH and Co.KG, Dudweiler Landstr, Germany, 2011. [28] Drazin P. and Riley N., The Navier-Stokes Equations:A Classification of Flows and Exact Solutions, Cambridge University Press (Reprinted), UK, 2007. [29] Nugroho G., Ali A.M.S. and Abdul Karim Z.A., A Class of Exact Solutions to the Three-Dimensional Incompressible Navier-Stokes Equations, Applied Mathematics Letters 23, 2010, pp. 1388 – 1396. [30] Nugroho G., Ali A.M.S., and Abdul Karim Z.A., Properties of V   and V     Classes of Solution to the Three-Dimensional Incompressible NavierStokes Equations, International Journal of Applied Mathematics and Mechanics 6 (11), 2010, pp. 98 – 117. [31] Holmes P., Lumley J.L. and Berkooz G., Turbulence, Coherent Structures, Dynamical Systems and Symmetry, Cambridge University Press, UK, 1996.

nc .

In: Turbulent Flows: Prediction, Modeling and Analysis ISBN: 978-1-62417-742-2 c 2013 Nova Science Publishers, Inc. Editor: Zied Driss

,I

Chapter 2

is

he rs

L ARGE -E DDY S IMULATION OF T URBULENCE -I NDUCED A ERO - OPTICAL E FFECTS IN F REE S HEAR AND WALL B OUNDED F LOWS

Pu bl

K. Volkov∗ School of Mechanical and Automotive Engineering, Faculty of Science, Engineering and Computing, Kingston University, London, UK

Abstract

N

ov

a

Sc i

en

ce

Optical aberrations induced by turbulent flows are a serious concern in airborne communication, imaging and optical systems because the quality of the beam degrades due to variations of the index of refraction along its path. The aero-optical aberrations produced by compressible flows are a direct consequence of the detailed physics in the flow. Time-accurate computational study of aero-optical distortions is performed based on large-eddy simulation (LES) technique because of its ability to account for the optical phase errors due to anisotropy of the flowfield and to resolve large-scale coherent structures at a reasonable computational cost. The re-normalization group (RNG) sub-grid scale model is used to account for the effect of the small turbulent eddies on the flowfield. The results of numerical simulation of turbulent flows in flat plate boundary layer, free mixing layer and free turbulent jet provide general statistics of optical distortions by examining the time-averaged intensity pattern and instantaneous flowfield. The flow-induced optical aberrations are studied for different inflow conditions. The computational tools and the results are beneficial in design and optimization of optical systems, where the control of laser beam propagation is an important problem, and coherent optical adaptive equipment and techniques based on the optical detection of light-scattering properties of the flow. Mitigation strategy is to use some sort of shear-layer control to regularize the vortex roll-up, thereby making it amenable to some form of adaptive-optic correction.

Keywords: turbulence, aero-optic, flow, compressibility, large-eddy simulation, boundary layer, mixing layer, jet, finite volume method ∗

E-mail address: [email protected]

28

K. Volkov

en

ce

,I he rs

Pu bl

is

Latin Symbols a Speed of sound A Aperture c Constant cp Specific heat capacity at constant pressure e Specific total energy f Damping function I Intensity k Turbulent kinetic energy l Length scale l Length M Mach number n Index of refraction p Pressure Pr Prandtl number q Heat flux R Radius Re Reynolds number Ri Richardson number t Time T Temperature S Strain rate tensor u, v, w Velocity components uτ Dynamic velocity vx , vr , vθ Cylindrical velocity components vx , vy , vz Cartesian velocity components x, r, θ Cylindrical coordinates x, y, z Cartesian coordinates

nc .

Nomenclature

N

ov

a

Sc i

Greek Symbols γ Specific heat capacity ratio Γ Circulation δ Thickness δ Filter width ε Dissipation rate of turbulent kinetic energy λ Thermal conductivity µ Dynamic viscosity ν Kinematic viscosity ρ Density τ Shear stress tensor ϕ Phase function ω Frequency Ω Rotation tensor

Large-Eddy Simulation of Turbulence-induced Aero-optical Effects ...

Superscripts n Time layer + Non-dimensional near-wall value

Introduction

en

1.

ce

Pu bl

Abbreviations CFD Computational Fluid Dynamics DNS Direct Numerical Simulation LES Large Eddy Simulation OPD Optical Path Difference OPL Optical Path Length RANS Reynolds-Averaged Navier–Stokes RNG Re-Normalization Group SGS Sub-Grid Scale SR Strehl Ratio TVD Total Variation Diminishing

he rs

,I

nc .

Effective Control volume Sub-grid scale Turbulent Wall Total Reference Free stream

is

Subscripts eff i, j, k sgs t w 0 ∗ ∞

29

N

ov

a

Sc i

Coherent vortical structures in free shear flows and wall bounded flows, separation of flow and unsteady wake impose a number of optically delirious effects on the laser beam projected from the turbulent flow. These optical effects are referred to as aero-optical effects, but in fact, this classification is itself broken into three distinctive effects: line-ofsight jitter, wavefront aberrations (usually reported with tip/tilt removed) and shock effects. Aero-optical distortions induced by turbulent flow are a serious concern of the airborne optical system, because the quality of the beam degrades due to variations of the index of refraction along its path. The refractive index of the optical window has the gradient distribution in non-uniform compressible flow. These randomly distributed variations of refractive index induce phase perturbations in a transmitted optical field that may result in bore-sight and centroid errors for tracking systems, blur and identification problems for imaging systems, and defocus and jitter for directed energy systems, any of which can substantially impact mission effectiveness. Amplitude and phase of the beam are subjected to fluctuations corresponding to the structural changes (broadening, deflection, and splitting) of the beam [30]. When an ini-

30

K. Volkov

z

Flow

x

ρ(x,y,z,t)

Coherent beam

he rs

Eddies

,I

θ

nc .

Disturbed wave-front

Undisturbed wave-front

is

Figure 1. Propagation of coherent beam through a turbulent flow.

N

ov

a

Sc i

en

ce

Pu bl

tially planar optical wavefront passes a compressible flow or thermal environment, different parts of the wavefront experience different density in the medium and hence have different propagation speeds, and the wavefront becomes deformed (Figure 1). The consequences of such deformations include optical beam deflection (bore-sight error) and jitter, beam spread, and loss of intensity. Wavefront distortions also cause reductions of resolution, contrast, effective range and sensitivity for airborne electro-optical sensors and imaging systems. A small initial deformation of the wavefront can lead to large errors on a distant target. An airborne optical beam generally encounters two distinct turbulent flow regimes [49]: the turbulence in the vicinity of the aperture produced by the presence of solid boundaries, and atmospheric turbulence. The effects induced by turbulent, variant-index flowfields is quantified by analysis of the far-field irradiance pattern. The short-path (near-field, y  A) aberration is usually termed aero-optics, and the long-path (far-field, y  A) aberrations are usually termed atmospheric propagation. In atmospheric applications an optical beam is required to be transmitted through a relatively long distance, over which the quality of the beam degrades due to variations of the index of refraction along its path. For air and many fluids, the refractive index is linearly related to the density of the fluid through the Gladstone–Dale equation, and therefore density fluctuations due to flow turbulence are the root cause of optical aberrations. The required spatial and temporal frequencies associated with aero-optical problems are at least an order of magnitude greater than those presently correctable by adaptive-optic systems for the atmospheric propagation case [28]. The solution of the problem is divided into three levels of understanding, arranged in increasing order of complexity [30]. The first level of understanding allows one to estimate the statistical level of optical distortion that is likely to be encountered for a particular optical propagation scenario through a flowfield. At the second level, one can predict not only the statistically averaged distortion but also the spatial and temporal frequencies associated with the time-varying optical distortion. The third level is reached when actual time

Large-Eddy Simulation of Turbulence-induced Aero-optical Effects ...

31

N

ov

a

Sc i

en

ce

Pu bl

is

he rs

,I

nc .

histories associated with the optical distortions can be calculated and linked to the specific structure of the flowfield. The recent development of aero-optical research are summarized in [30,47,48]. A comparison study of atmospheric propagation and aero-optics, which distinguished the aerooptical transmission from atmospheric propagation in terms of the optically distorting characteristic parameters of flowfields, is performed in [42, 43]. The aero-optical interactions along the propagation path in turbulent compressible separated shear layers through experiments on direct imaging of the refractive index field, which investigate the quantities of the interaction optical path difference and its root-mean-square value as functions of propagation distance along the laser beam direction as well as functions of the laser aperture size, is examined in [74]. Approaches at investigating aero-optical effects employ both computational fluid dynamics (CFD) and experiments. As turbulence remains an unresolved problem, CFD solutions and the consequent aero-optical methods generally employ empirical turbulence models and ideal assumptions, for instance, the statistical Sutton model [47, 49] and the turbulent model consisting of spherical vortices with different sizes and distributions [51]. CFD provides three options, including direct numerical simulation (DNS), large-eddy simulation (LES) and Reynolds-averaged Navier–Stokes (RANS) equations. DNS is a simulation in which the instantaneous Navier–Stokes equations are numerically solved. It resolves the whole range of spatial and temporal scales of the turbulence. As a result, the computational cost of DNS is very high and it is seldom used in practical calculations. LES calculates the large scale flow motions while models the smaller ones with the subgrid scale model (SGS), but it is also a simulation that employs the instantaneous (filtered) Navier–Stokes equations to obtain temporal flow characteristics. Most of the early computational studies were based on statistical analysis with simplifying assumptions such as homogeneous and isotropic turbulence, and therefore are not directly applicable in realistic aero-optical flowfields. The mesh resolutions were poor in the cases involving complex geometries, and the numerical schemes employed were either total variation diminishing (TVD) or upwinding techniques that are highly dissipative. The commonly used models of fluctuations of index of refraction are reviewed, and the aero-optical distortions produced by the largest-scale flow structures in the boundary layer on a flat plate, free mixing layer and free turbulent jet are simulated. The semi-empirical correlations, analytical model and numerical model are developed to predict the dispersion of phase function. Time-accurate computational study of aero-optical distortions is performed based on large-eddy simulation because of its ability to account for the optical phase errors due to anisotropy of the flow field and to resolve large-scale coherent structures at a reasonable computational costs. The results obtained are compared with available computational and experimental data.

2.

Modelling and Simulation

Historical development of aero-optical computational and experimental studies is reviewed in [30]. Research in the area of turbulent distortions of optical waves is traced back to the 1950s and 1960s [8,50] when the scattering of acoustic and electromagnetic waves due to random

32

K. Volkov

N

ov

a

Sc i

en

ce

Pu bl

is

he rs

,I

nc .

fluctuations of refractive index were studied, mostly in the context of atmosphere propagation. The first users of aero-optics information were fluid mechanics researchers who took advantage of the variable index of refraction that accompanies variable density flows (shadowgraph systems, interferometric techniques, airborne telescopes, ground-based astronomical systems). In this case, shear and boundary layers over the viewing aperture degrade the system performance. Most of the early studies were based on statistical analysis with simplifying assumptions such as homogeneous and isotropic turbulence, and therefore were not directly applicable in realistic aero-optical flowfields. Optical and flow parameters for the case of homogeneous and isotropic turbulence were discussed in [49], where statistical model of prediction of far-field optical aberrations was developed. It was in the late 1980s when aero-optics in the modern sense, i.e. the study of optical distortions due to near-aperture turbulence, came into consideration. Many experimental studies have been performed to develop high-speed wavefront measurement tools [9, 30], study the refractive index structures [5, 14, 18], apply distortion scaling laws [23, 26], and devise control techniques to suppress or modify optically important turbulence structures [22, 44]. Despite advances in wavefront sensor technology, significant limitations still exist in terms of spatial and temporal resolutions. Computational studies, when performed adequately, allow to probe the flow and optical fields in greater detail and hence complement experiments to further understanding in this area. A time series of realizations of the shear layer density field was performed in [39, 65] for temporally developing shear layers and in [33] for spatially developing shear layers. These approaches yield a solution for the tested conditions only, providing little insight into how the density field varies with changing initial and boundary conditions or the actual mechanism causing the density changes. Different interferometric methods of high-speed flow analysis (direct interferometry, pulsed interferometry, shearing interferometry) provide a time-averaged assessment of the optical phase variance over the aperture by quantifying the phase variance of individual interferograms and averaging the variance over the number of interferograms used to quantify the aperture aberrations. However, these methods provide no information concerning temporal frequencies [4, 6, 14, 70]. Indirect methods of assessing the near-field optical phase variance depend on a quantification of the density fluctuations which are generally experimentally measured using hot-wire methods. Significant progress has been made in both measuring the wavefront dynamics for the atmospheric propagation problem and using this information for designing and applying adaptive optic equipment and techniques. Aero-optical distortions created by shear layers with convective Mach numbers of 0.15, 0.54, and 0.96 were studied experimentally in [13]. The shear layers were produced with helium or nitrogen as the higher speed flow mixing with a lower speed ethylene flow. While the ethylene and nitrogen were density matched, the ethylene produced good Rayleigh scattering of a laser light sheet projected through the flow. Photographs of the resulting flow visualization effectively gave the instantaneous number density of ethylene in the flowfield. From the number density field, the index field was deduced, and the resulting wavefront distortion was computed. Unfortunately, the use of two fluids of such disparate index, while

Large-Eddy Simulation of Turbulence-induced Aero-optical Effects ...

33

N

ov

a

Sc i

en

ce

Pu bl

is

he rs

,I

nc .

clearly incorporating Mach effects on mixing, effectively masked any effects of compressibility alone on wavefront distortion. Low-order descriptions of turbulence include proper orthogonal decomposition (POD) analysis, linear stochastic estimation (LSE), time series of multiple-view wavefronts. The key issue in modeling aero-optical effects is the choice of the method used to calculate the field of the fluctuating parameters of the turbulent flow, in particular, the field of density fluctuations [14, 30] (the density is related to the refractive index by the Gladstone– Dale law). A number of semi-empirical models of different levels of complexity were developed for numerical simulations. The most popular models are the model of homogeneous turbulence [47] and the model of inhomogeneous turbulence [48]. In models [47, 48, 52], the dispersion of density fluctuations is taken as the difference between the medium densities near the wall and at the boundary layer edge. In the approach used in [55], the field of density fluctuations is related to temperature fluctuations by the equation of state (pressure fluctuations are neglected). The temperature fluctuations are found with the use of the Reynolds analogy and Prandtl mixing length model. Semi-empirical models based on the assumption of flow quasi-steadiness were developed in [23, 25, 42]. For the quasi-steadiness condition to be satisfied, the response time of the system and the signal propagation time should not exceed the frozen turbulence time, which is of the order of 10−3 ÷ 10−2 s. The analytical and numerical models developed in [16] had two distinct parts. First, the time-varying velocity flowfield was modeled using a discrete vortex method (DVM). The DVM provided a common velocity-field time series from which different compressibility mechanisms were evaluated. Using each velocity flowfield time series as an input, a second numerical code implemented the density or index of refraction model. In order to account for the effect of pressure fluctuations on the instantaneous density field, a biquadraticsurface-fitting method was used to iteratively solve the unsteady two-dimensional Euler equations. A mathematical model of the optical window under the thermal environment, and analysis of the optical path difference (OPD) caused by the aero-optical effects were developed in [73]. An analytical method to calculate the optical characteristics of the window was performed. A numerical simulation of the optical window was presented, which evaluated temperature field, deformation field, stress field and OPD distributing. The OPD is influenced mostly by the temperature (thermal optical effect). The influence of the deformation and the stress (elastic optical effect) is much smaller. The study of [69] employs CFD results as the source data to discuss the problem of sparse measurement in aero-optics. For this purpose, proper orthogonal decomposition (POD) approach has been employed in aero-optics to reduce the amount of information required to compute or reconstruct the aberrated wavefront, so that the system temporal frequency could be increased. However, POD approach also requires empirical information of the reconstructed mode number to define its represented energy. A neural network (NN) model for high speed turbulent boundary layer inducing aero-optical distortions was proposed to solve the spatial frequency problem in sparse measurements. The temperature fluctuation, a measurable flow state, is set as the NN model input, and the OPD is outputted. Decomposition of the instantaneous velocity and pressure fields into a time-averaged quantity, a random (incoherent) component that describes small-scale motion, and a peri-

34

K. Volkov

N

ov

a

Sc i

en

ce

Pu bl

is

he rs

,I

nc .

odic (coherent) component corresponding to large vortices, was used in [32]. The vortex viscosity hypothesis was applied to close the system of the governing equations. The majority of researches on aero optics have concentrated on the computational approximations using geometrical optics, mainly about the computation of optical path differences or optical path length. Earlier computational approaches for aero-optics typically involved RANS calculations using a turbulence model such as the k–ε model, from which crude algebraic relationships were used to obtain the root-mean-square of intensity and length scales of the index of refraction field. A more sophisticated approach was used in [46], in which a model transport equation for the root-mean-square refractive index fluctuations was solved. The turbulence information was then fed into an optics model based on geometric optics, which predicts properties of the beam such as amplitude loss and spreading [30]. Time-accurate computational studies of aero-optical distortions were performed in [55, 56], where DNS of a homogeneous shear flow and turbulent channel flow was performed to study the induced optical wavefront errors. The simulations were based on incompressible flow equations at relatively low Reynolds numbers, and the fluctuating refractive index was modeled by a passive scalar. The anisotropy of the optical phase errors due to anisotropy of the flow was noted. Because of its ability to resolve large-scale motions at a reasonable computational cost, LES has recently been used for aero-optics. LES of a compressible turbulent mixing layer and ray tracing through it was performed in [11]. LES of aero-optical distortions in a fuselage/turret configuration was reported in [29], where ray tracing was used to obtain the wavefront error and represented it in terms of Zernike polynomials, and in [24]. Also, instantaneous far-field intensity patterns using the Fraunhofer approximation was studied. LES of aero-optics of a supersonic boundary layer was used in [52, 54]. Experiments and LES to investigate control of flowfields to mitigate the distortion of a laser beam passing through a cavity shear layer were used in [44]. In these works, the major contribution to optical distortions was assumed to be from the large resolved scales of the flow. The grid resolutions were poor in the cases involving complex geometries, and the numerical schemes employed were either TVD or upwinding techniques that are highly dissipative. This can severely impact the effectiveness of SGS model and artificially damp the small resolved scales of the flow, which can be important [34]. LES of turbulence-induced aero-optical effects in mixing layers and free jets was performed in [58–60, 64] based on unstructured finite volume CFD code. The paper [24] summarizes recent efforts to investigate the fluid dynamics and aerooptical environments around wall-mounted hemisphere on cylinder turrets with either conformal or flat windows at moderate Mach numbers between 0.3 and 0.5. In addition, results for a hemi sphere turret were given. A short discussion of beam jitter effects and flow topologies around turrets at transonic and supersonic flows was also provided. The density fluctuations are important to the aero-optical calculations, however, the conventionally employed Faver filtering equation does not contain the related term of density fluctuations and the physical meanings of filtered parameters are not clear. LES results provide a set of three-dimensional flowfields in grids at discrete time steps [2, 53]. These fields are utilized to calculate the relevant OPD [29]. Diffraction is always present and is a dominant effect especially near sharp edges where an edge can be either an amplitude edge

Large-Eddy Simulation of Turbulence-induced Aero-optical Effects ...

35

3.

ce

Pu bl

is

he rs

,I

nc .

(i.e. an aperture and the edge of a vehicle body in the flow) or a phase edge (high refractive index gradients) [38], however, this study focused on the near-field aero-optical effects of the boundary layer over a side-mounted window. The non-uniform flowfield around high speed flying vehicle was studied in [71]. The computations were performed on a typical conical-head flying vehicle. The mean flow density distribution, the imaging deviations, the propagation path distances in the non-uniform flowfield, and the density distributions along propagation paths for two different flying cases were computed. The RANS approach was adopted for the density computation based on realizable k–ε model. Three flowfield related factors were considered in order to reduce the deviation: the propagation path distance in the non-uniform flowfield, the angle of incidence at the freestream boundary, the density distribution along propagation path in the non-uniform flowfield. High-order compact differences were used in [68] to solve the parabolic beam equations for the Laplacians and Runge–Kutta integration along the beam path in the realm of aerooptics. The grid-based model was used in [66, 73] to study the aero-thermal optical effect and the aerodynamic optical effect near side-mounted optical window. A computational model for studying aero-optical imaging through the supersonic flowfields was presented in [66]. This model integrates the CFD grids with the model of angular spectrum propagation to construct serially connected optical sub-systems for representing the supersonic flowfields. The CFD grid model developed in [66] was thought as the foundation of the numerical simulation. The aero-optical propagation was viewed as the optic angular spectrum of plane wave transmitting grid by grid, and the total optical transfer function of such flowfields was derived and further digital image processing method was utilized to simulate the aerooptical imaging through supersonic flowfields.

Optical Flow Visualization

N

ov

a

Sc i

en

Interferometry, schlieren imaging and shadowgraph are techniques which are widely used in optical flow visualization. These optical techniques use the principle that as light passes through a flowfield its phase and direction are changed due to variations of the refractive index induced by non-uniform density distribution in the flowfield. Calculating light refracting through a flow presents a number of challenges. Light paths are recomputed whenever the viewpoint changes thus an interactive method for determining them at each frame is presented. This makes it possible to analyze physical phenomena that manifest themselves as density changes in the flowfields. Shadowgraphs provide an intuitive way of looking at small changes in flow dynamics through caustic effects while schlieren cutoffs introduce an intensity gradation for observing large scale directional changes in the flow. The idea of shadowgraph technique is that small changes do not scatter light to a large degree but it was noticed that shining a bright light through them produces a clear image of the flow by looking at the shadows formed from light refraction. In a shadowgraph system refracted light is imaged on a film plane. Shadowgraphs look at changes in the second derivative and are a poor indicator of the amount or direction of refraction. If all rays were

36

K. Volkov

I∼∇

ZL0

nc .

refracted the same amount in the same direction then the resulting image would be identical to a translated image of no refraction at all. For shadowgrams, the image intensity is proportional to the gradient of the integration of ∇ρ(x, y, z) in the direction perpendicular to the light ray, which is calculated by

∇ρ(x, y, z)dl.

,I

0

en

ce

Pu bl

is

he rs

Schlieren photographic techniques provide additional information by introducing a onedimensional cutoff that shifts intensity values based on the amount and direction of displacement at the focused cutoff region. In the schlieren system the light source is then refocused in a small area and a cutoff is inserted to reduce light from the light source. A vertical knife-edge is then inserted at the center of the re-focused light source. If no light is refracted then the knife-edge reduces the light source by half resulting in a gray image. Refracted light causes shifts in the focused image of the original light source resulting in more or less of the focused light being blocked by the cutoff. If the focused image is shifted down the resulting region is darker and if shifted up then more of the original light gets through to the film plane. A knife-edge cutoff provides information about the amount of light shifted along a single axis. Another common type of cutoff is a circular cutoff that shades the image based on the amount of displacement without the directional information of the knife-edge. Color filters are also used as a cutoff to produce colors based on the direction of displacement. The intensity of each point in schlieren photographs is proportional to the density gradient perpendicular to the knife edge because the ray deflection is proportional to the density gradient. When the knife is set to be perpendicular to the x-axis in the physical space the intensity, I, is proportional to the integration of the density gradient along the ray ∂ρ(x, y, z) dl. ∂x

0

a

Sc i

Since both the schlieren and shadowgraph methods provide only qualitative information, and these methods are not used as frequently as interferometry in applications to the CFD validation. Interferometry tracks changes in phase-shift resulting in bands appearing. Interferometry allows a different view of data by tracking phase-shift through the flow producing visual bands. Custom color filters allow for exploring specific regions in the flow. Multi-field data present a difficult problem and an interesting exploration of schlieren visualization. Interferometry differs from schlieren and shadowgraph images by looking at phase shift instead of refraction. When light travels through a disturbance and encounters a change in refractive index the speed of that light changes resulting in a phase shift. The idea behind interferometry is to directly measure this phase shift and display it providing a picture of changes in refractive index. In holographic interferometry, double exposure interferograms are generated by exposing the film to the object and the reference beam. For infinite-fringe interferometry, the object beam passes through the flowfield, and its phase changes due

ov

N

I∼

ZL0

Large-Eddy Simulation of Turbulence-induced Aero-optical Effects ...

37

ZL0h i ρ(x, y, z) − ρ0 dl, 0

,I

2π G ∆φ(xi , yi ) = λ

nc .

to variations of the refractive index caused by density changes between exposures but the reference beam does not change. The phase shift of the object beam relative to the reference beam between exposures is calculated by integrating density distribution

he rs

where λ is the wavelength of the light, ρ0 the density of the undisturbed flow, G is the Gladstone–Dale constant that changes depending upon gas species and varies slightly with the light wavelength, and L0 the length of the light path through the test section. The image intensity of the infinite-fringe interferograms I is calculated by h i I = 1 + cos ∆φ(xi, yi ) + φ0 ,

Pu bl

is

where (xi , yi) denotes the image plane, and φ0 is an initial phase shift to compensate for any phase shift between two exposures and is taken as zero in most of the cases. In the case of finite-fringe interferograms achieved by tilting the reference beam between exposures, the image intensity is given by h i I = 1 + cos ∆φ(xi , yi ) + 2πvxxi + 2πvy yi , where vx and vy are the special frequency components of the unperturbed fringes.

Performance of Aero-optical System

ce

4.

Sc i

en

A measure of the effect of optical aberrations caused by the flow is to examine its effect on the focusality of the laser in the far-field. This measure is usually expressed in terms of a Strehl ratio, SR, which is the tilt-removed, on-axis, far-field intensity of the aberrated beam divided by that of the unaberrated beam (diffraction-limited beam). The Strehl ratio is defined as SR(t) =

I(t) , I0

N

ov

a

where I(t) is the maximum light irradiance in the far-field pattern for an aberrated system, and I0 is the maximum light irradiance for the optical system without aberrations (diffraction-limited on-axis irradiance). According to the large-aperture approximation, an estimation of the time-averaged Strehl ratio for a given time-averaged optical-phase variance, σϕ2 , is described by  SR = exp −σϕ2 . The phase variance is the normalized OPD variance   2πOPDrms 2 2 σϕ = . λ

38

K. Volkov

,I

nc .

where λ is the wavelength of the optical beam, and OPDrms is the spatial root-mean-square of OPD(x, y, z, t). Aberration-free optical systems are still limited in the amount they focus by the wave nature of light. The smallest realizable focal spot is called the diffraction-limited spot. The irradiance profile for this spot is shown using physical optics techniques, and it is given by the analytical relation for a circular aperture   2J1 (r) 2 I(r) = , r

I0 =

P P d2 ∼ , πY 2 λ2 /d2 Y 2 λ2

he rs

where J1 (r) is the first-order Bessel function, and r is the relative distance from the center of the beam. The maximum irradiance of this central lobe is given by

Pu bl

is

where P is the laser output power, Y is the propagation distance (the focal length of the optical system), and d is the aperture diameter. The maximum irradiance increases with decreasing wavelength as shown in the Figure 2, where the peak irradiance is nondimensionalized by the peak irradiance for a wavelength of 1 mm. This increase in I0 with decreasing λ encourages the use of lasers operating at visible wavelengths.

I0 (λ)/I0 (1 µm)

1

N

ov

a

0

0

CO 2 (λ = 10.6 µm)

0.2

DH (λ = 3.8 µm)

Sc i

0.4

HF (λ = 2.8 µm)

en

HeNe (λ = 0.63 µm)

0.6

COIL (λ = 1.315 µm)

ce

0.8

2

4

6

8

12

10

λ, µm Visible/Near visible

Infrared

Figure 2. Peak irradiance of diffraction-limited spot.

The variation in SR with wavelength is illustrated in the Figure 3, which shows that whether a given flowfield (constant OPDrms) presents an aero-optical problem is a strong function of the laser’s wavelength. An aberration that is insignificant (SR = 0.95 for the

Large-Eddy Simulation of Turbulence-induced Aero-optical Effects ...

39

CO2 laser) is suddenly become a large problem for a shorter wavelength laser (SR = 0.04 for the chemical oxygen-iodine laser, COIL).

nc .

I/I0

1

,I

0.8

he rs CO 2 (λ = 10.6 µm)

is 2

0

Pu bl

0

DH (λ = 3.8 µm)

0.2

HF (λ = 2.8 µm)

0.4

COIL (λ = 1.315 µm)

HeNe (λ = 0.63 µm)

0.6

4

6

8

10

12

λ, µm

Visible/Near visible

Infrared

5.

ce

Figure 3. Effect of wavelength on aberrated-system performance.

Phase Function of Wave Front

Sc i

en

The propagation of electromagnetic wave is described by the Maxwell equations. It is assumed that the time scale related to the wave propagation is much smaller than the turbulence time scale, the refractive index is time-independent, and the medium is nonconducting and has a constant magnetic susceptibility. Propagation of the wave front is described by the equation ∇2 E −

n2 ∂ 2 E = 0, c20 ∂t2

(1)

N

ov

a

where E is the electric field strength. The refractive index is given by n = c0 /c, where c and c0 are the speeds of light in the medium and vacuum, respectively. For a monochromatic sinusoidal wave with frequency ω, the equation (1) has an solution h i E(r, t) = E 0 exp ϕ(k, r) − ωt , (2) where E is the field strength at point r at time moment t, and E 0 is the field strength at the point r = 0. The wave number, k, frequency, ω, and wavelength, λ, are coupled by the relationship k = ωc = 2π/λ.

40

K. Volkov

α(x, t) =

dϕ(x, t) , dt

ϕ(x) =

Zt

,I

where ϕ(x, t) is the phase function of the wave front. Integration over time gives

nc .

The angle describing the direction of the beam that passed through a test section of the variable-density medium is determined as the time derivative of the optical path of the beam

α(x, t) dt.

he rs

0

The main parameter of interest for practical applications is the difference between the instantaneous and mean values σϕ (x) = ϕ(x) − hϕ(x)i .

Pu bl

is

The transition from the time formulation of the problem to the space formulation is performed with the Taylor hypothesis of frozen turbulence. The optical path length is determined by integrating the distribution of the refractive index along the beam propagation direction [47] Z h i ϕ(x) = 1 + n(x) dl. (3)

ce

The density field computed by CFD is converted into the refractive index field. The dependence of the refractive index on fluid density is described by the Gladstone–Dale law [30] h i n(x) = 1 + G ρ(x) − ρ∞ , (4)

Sc i

en

where ρ∞ is the density of the unperturbed medium, and G is the Gladstone–Dale constant. Gladstone–Dale coefficient is weakly dependent on the wavelength and is approximately expressed by   7.52 · 10−3 −4 G(λ) = 2.23 × 10 1+ (m3 /kg), λ2

N

ov

a

where λ is wavelength of the infrared light, and its units is µm (G = 2.23 × 10−4 m3 /kg for air, and the weak dependence on the wavelength is neglected). Substitution of equation (4) into equation (3), gives Z n o ϕ(x) = 1 + G [ρ(x) − ρ∞ ] dl. The wave front is defined as the surface normal to the direction of wave propagation ∇ϕ =

∂ϕ ∂ϕ i+ j. ∂x ∂y

Large-Eddy Simulation of Turbulence-induced Aero-optical Effects ...

nc .

The phase function of the wave front is given by Z ϕ = k(x, t) · dx,

41

∆n(x, y, z, t)dz,

is

∆ϕ(x, y, t) = k

Zz2

he rs

,I

where k = 2πl/λ is the wave vector, and l is a unit vector specifying the wave propagation direction. When a plane monochromatic wave passes through a test section with a varying refractive index, its amplitude remains almost constant and the phase is changed, so that ϕ = ϕ0 + ∆ϕ, where ∆ϕ is the phase shift due to the inhomogeneity of the medium. The field strength of a perturbed wave is obtained by multiplying equation (2) by factor exp (i∆ϕ), where ∆ϕ(x, y) denotes the phase shift along the optical path length. The phase distribution in plane (x, y) normal to the direction of wave propagation at time t is found by integrating the distribution of the refractive index over the test section thickness (the beam deflections are neglected)

z1

Pu bl

where ∆n(x, y, z, t) is the variation of the refractive index along the beam propagation direction at time t, and L = z2 − z1 is the thickness of test section. In practice, rather than measuring an absolute optical path length (OPL), only the relative difference in OPL across the aperture or optical path difference (OPD) is important, and its is defined as OPD(t, x) = OPL(t, x) − OPD(t),

Dispersion of Phase Function

en

6.

ce

where the overbar denotes the spatial averaging over the aperture.

Sc i

Dispersion of small-scale fluctuations of density, σρ2 , and the corresponding correlation length scale, lρ, are related to the dispersion of wave phase, σϕ2 , by the expression [47] σϕ2

= αβ

2

ZL

σρ2 lρ dy,

β=

2π dn = kG(λ), λ dρ

0

a

where L is the optical path length, and the integral is taken across the boundary layer. The constant, α, depends on the form of the correlation function assumed for the density fluctuations (α = 2 for the exponential correlation function and α = π for the Gaussian correlation function). According to the equation (5), this leads to an 11% difference in the dispersion of phase fluctuations. Length scale, lρ , is found by integrating the correlation function

ov

N

(5)

+∞ Z lρ = Rρρ (y)dy. −∞

42

K. Volkov

nc .

In the conditions of local equilibrium, the correlation length scale coincides with the correlation scale of velocity fluctuations, lρ ∼ lu ∼ k3/2 /ε. This assumption is not applicable near the wall, where velocity fluctuations turn to zero, therefore k = 0 and lu = 0, but lρ 6= 0. To estimate the phase fluctuations in the near-wall region, the equation (5) is usually replaced with a semi-empirical equation [47] σϕ2 = β 2 ly δσρ2 ,

,I

(6)

σρ2 = A2 (ρw − ρ∞ )2 ,

he rs

where ly is the integral turbulence length scale in the direction normal to the wall. It is assumed that ly ∼ 0.1δ, where δ is the boundary layer thickness. The dispersion of the density fluctuations is estimated as (7)

σρ2 (y)ly (y)dy.

Pu bl

σϕ2 = β 2

ZL

is

where ρw and ρ∞ are the density at the wall and the density in the free flow, and A =0.1– 0.2. At ly  δ, formula (6) is refined by integrating across the boundary layer [47, 52]

(8)

0

Sc i

en

ce

In models of [47, 48, 52] based on equations (6) and (8), the dispersion of density fluctuations is determined as the difference between the density near the wall (inner region) and the density in the outer region of the boundary layer. The applications of the models of [47,48] are discussed in [52]. In the models of [55,67], the density fluctuations are found by solving the transport equation of a passive scalar. If pressure fluctuations are disregarded, density fluctuations are related to temperature fluctuations based on the equation of state [55, 67]. The temperature fluctuations are determined using the Reynolds analogy between the velocity and temperature transport and the Prandtl mixing length model [27], whilst the accounting for the fact that the corresponding linear scales satisfy the relationship lu /lT = Prt , where the turbulent Prandtl number is taken to be constant.

7.

Vortex Model

N

ov

a

To study the nature of optical aberrations, a vortex model of the flow in the boundary layer and in the mixing layer is considered, which is based on the exact solution of the plane Lamb vortex problem. The pressure inside the vortex is not constant, which is responsible for changes in density and refractive index. The pressure and velocity are related by the equation (the viscous effects are neglected) u2θ 1 ∂p =− . r ρ ∂r

Large-Eddy Simulation of Turbulence-induced Aero-optical Effects ...

43

nc .

In the limiting cases corresponding to forced and free vortices, the distribution of tangential velocity is described by the relation ( U r/R if r  R, uθ (r) = (9) U R/r if r  R,

he rs

,I

where R is the vortex radius and U is the maximum velocity expressed in terms of the circulation and vortex radius, U = Γ/(4πR). As the vorticity takes a constant value in the forced vortex and is equal to zero in the free vortex ( Ω if r  R, ω= 0 if r  R.

r/R , 1 + (r/R)2

(10)

Pu bl

uθ (r) = 2U

is

A combined vortex consists from the inner core rotating as a solid with a constant vorticity, Ω, and the outer area where the vorticity tends to zero. The velocity distribution described by the equation (9) is not differentiable at the point r = R. To construct a smooth velocity distribution, the interpolation formula is used

which satisfies the limiting relations (9). The difference between the velocity profiles described by equations (9) and (10) decreases with decreasing vortex radius. Under the assumption of validity of the isentropic conditions   p ρ γ = , p∞ ρ∞

ce

the density distribution is found from the relation  1/(γ−1) (U/a∞ )2 ρ(r) = ρ∞ 1 − 2(γ − 1) , 1 + (r/R)2

en

satisfying the limiting cases

ρ(r) =

(

Sc i

where the density at the vortex center is "

ρmin

if r  R,

ρ∞

if r  R,

ρmin = ρ∞ 1 − 2(γ − 1)



U a∞

(11)

2 #1/(γ−1)

N

ov

a

At ∆ρ/ρ∞  1, the density distribution is described by the formula ρ − ρ∞ 2 (U/a∞ )2 . =− ρ∞ γ 1 + (r/R)2

(12)

where a∞ is the speed of sound at free stream. As the density at the vortex center is smaller than that at the vortex periphery, optical aberrations appear. The flow is replaced by an infinite chain of vortices located on one line at an identical distance from each other and having a circulation Γ. The resultant velocity distribution is described by the relations derived in [40].

44

8.

K. Volkov

Governing Equations

nc .

In Cartesian coordinates (x, y, z), an unsteady three-dimensional flow of a viscous compressible fluid is described by the following equation written relative to filtered quantities ∂Q ∂Fx ∂Fy ∂Fz + + + = 0. ∂t ∂x ∂y ∂z

(13)

he rs

,I

Equation (13) is complemented with the equation of state of ideal gas    1 2 2 2 p = (γ − 1)ρ e − v + vy + vz . 2 x



ρvx ρvx vx + p − τxx ρvxvy − τxy ρvx vz − τxz (ρe + p)vx − vx τxx − vy τxy − vz τxz + qx

ce

  Fx =   

Pu bl

is

The vector of conservative variables, Q, and the flux vectors, Fx , Fy and Fz , have the following form   ρ  ρvx     Q=  ρvy  ,  ρvz  ρe



Sc i

en

  Fy =   

ρvy ρvy vx − τyx ρvy vy + p − τyy ρvy vz − τyz (ρe + p)vy − vx τyx − vy τyy − vz τyz + qy

N

ov

a



  Fz =   

ρvz ρvz vx − τzx ρvz vy − τzy ρvz vz + p − τzz (ρe + p)vz − vx τzx − vy τzy − vz τzz + qz



  ,   

  ,   

  .  

The components of viscous stress tensor and components of heat flux vector are found from the relations   ∂vj ∂vi ∂T 2 ∂vk τij = µeff + − δij , qi = −λeff . ∂xj ∂xi 3 ∂xk ∂xi

Large-Eddy Simulation of Turbulence-induced Aero-optical Effects ...

45

he rs

,I

nc .

Here, t is the time, ρ is the density, vx , vy , and vz are the velocity components in the coordinate directions x, y, and z, p is the pressure, e is the total energy per unit mass, T is the temperature, and γ is the specific heat ratio. The equation (13) is suitable for both laminar and turbulent flows, and it formally coincides with unsteady RANS equations. The effective viscosity, µeff , is calculated as a sum of molecular viscosity, µ, and eddy viscosity, µsgs , and the effective thermal conductivity, λeff , is expressed in terms of viscosity and Prandtl number   µsgs µ µeff = µ + µsgs , λeff = cp , + Pr Prsgs where cp is the specific heat capacity at constant pressure. Molecular and sub-grid Prandtl numbers are Pr = 0.72 and Prsgs = 0.9 for air. The Sutherland’s law is used to obtain molecular viscosity as a function of temperature 

T T∗

3/2

T ∗ + S0 , T + S0

is

µ = µ∗

ce

Pu bl

where µ∗ = 1.68 × 10−5 kg/(m s), T∗ = 273 K and S0 = 110.5 K for air. Modelling of unsteady subsonic flows faces the problem of specification of the boundary conditions for the fluid outflow boundary containing intense vortex structures. Possible non-physical effects of generation and reflection of acoustic waves on the outlet boundary distort the real pattern of the flow. The non-reflecting boundary conditions are applied to all variables on the outlet boundary. Cylindrical coordinates (x, r, θ) are used to specify the boundary conditions instead of Cartesian coordinates (x, y, z). The radial velocity, vr , and the tangential velocity, vθ , are related to Cartesian velocities by the relations vy y + vz z

(y 2

+ z 2 )1/2

,

vθ =

vz y − vy z (y 2 + z 2 )1/2

.

en

vr =

Sc i

In view of the symmetry of computational domain and boundary conditions, the timeaveraged tangential velocity is zero.

9.

Sub-grid Scale Model

N

ov

a

To close the filtered Navier–Stokes equations, written in the form (13), the eddy viscosity model is used. The inclusion of compressibility in the SGS model slightly affects the results of numerical simulation [35]. In the Smagorinsky model [45], the eddy viscosity is calculated as µsgs = ρ (CS ∆)2 |S| ,

(14)

where |S| = (2Sij Sij )

1/2

,

1 Sij = 2



∂vi ∂vj + ∂xj ∂xi



.

46

K. Volkov

nc .

It is usually assumed that CS ∼ 0.1. To take into account the influence of the wall on the mixing length, CS ∆, the equation (14) is supplemented by the Van Driest damping function "   # 3 y+ fµ = 1 − exp − , 25

X=

µ2sgs µeff , µ3

(15)

is

h i1/3 µeff = µ 1 + H(X − C) ,

he rs

,I

where y + = ρuτ y/µ, uτ = (τw /ρ)1/2, and τw is the wall shear stress. The RNG model is simple and efficient from the computational point of view and storage requirements, and in contrast to the Smagorinsky model correctly predicts the behavior of eddy viscosity in the laminar flow region without any additional modifications [72] (µsgs → 0 if Re  1). In the RNG model, the calculation of effective viscosity reduces to the solution of nonlinear equation [72]

Pu bl

where H(X) is the Heaviside function, and C = 100. The eddy viscosity is found from the relation µsgs = ρ (CR ∆) |S|2 ,

(16)

where CR = 0.157. At X  C, the equation (16) yields the Smagorinsky model (14), and the model constant is (2CR )1/4 = 0.119. 2π

ce

CS =

en

To incorporate the effects of curvature of streamlines, the eddy viscosity is multiplied by the damping function depending on the Richardson number [41] f (Ri) = (1 − αRi)1/2

(α ∼ 0.1).

Sc i

The Richardson number is found by the formula   |Ω| |Ω| Ri = 1+ , |S| |S|

N

ov

a

where

|Ω| = (2Ωij Ωij )1/2 .

The filter width, ∆, is related to the control volume ∆ = V 1/3 = (∆x∆y∆z)1/3,

where V is the volume of cell, ∆x, ∆y, and ∆z are mesh steps in the coordinate directions x, y, and z, respectively.

Large-Eddy Simulation of Turbulence-induced Aero-optical Effects ...

47

nc .

In the boundary layer, the mesh step in the direction normal to the wall, ∆x, is replaced b and the filter width is found as by the ∆x,  1/3 b ∆ = ∆x∆y∆z .

he rs

,I

b = ∆x near the wall, and ∆x b = ∆x The modified representation of the filter width gives ∆x far from the wall. The mesh step size, ∆x, smoothly changes between the thresholds. The modified mesh step size, ∆x, is the averaged of ∆x in the near-wall region, and the modified b is found from the formula mesh step size, ∆x,     −1/α 1 α 1 α b ∆x = (α = 3). + ∆x ∆x

is

A smaller width of the filter allows reproduction of a wider frequency range of fluctuations of flow parameters, whereas an increase in filter width facilitates smoothing of the solution (LES transforms to DNS as ∆ → 0).

Pu bl

10. Numerical Method

Sc i

en

ce

Numerical calculations are performed with unstructured in-house finite volume CFD code [61–63]. The code uses an edge-based data structure to give the flexibility to run on meshes composed of a variety of cell types. The fluxes are calculated on the basis of flow variables at nodes at either end of an edge, and an area associated with that edge (edge weight). The edge weights are pre-computed and take account of the geometry of the cell. The non-linear CFD solver works in an explicit time marching fashion, based on a three-step Runge–Kutta stepping procedure and piecewise parabolic method [63]. The governing equations are solved with Chakravarthy–Osher scheme for inviscid fluxes [7], and the central difference scheme of the 2nd order for viscous fluxes. Convergence to a steady state is accelerated by the use of multigrid techniques [61], and by the application of block-Jacobi preconditioning for high-speed flows, with a separate low-Mach number preconditioning method for use with low-speed flows [62]. The sequence of meshes is created using an edge-collapsing algorithm. The re-normalization group (RNG) sub-grid scale model is used to account for the effect of the small turbulent eddies on the flowfield. Equation (13) is written in the form dQni,j,k dt

+ L(Qni,j,k ) = 0,

N

ov

a

where L(Qni,j,k )

= +

(Fx )ni+1/2,j,k − (Fx )ni−1/2,j,k ∆xi,j,k n (Fz )i,j,k+1/2 − (Fz )ni,j,k−1/2 ∆zi,j,k

+

(Fy )ni,j+1/2,k − (Fy )ni,j−1/2,k ∆yi,j,k

+

.

The subscripts i, j and k refer to the control volume, and the superscript n refers to the time layer.

48

K. Volkov The three-step Runge–Kutta method is used for discretization in time [36] (1)

(n)

(n)

i 3 (n) 1 h (1) (1) Qi,j,k + Qi,j,k + ∆tL(Qi,j,k ) , 4 4 i 1 (n) 2 h (2) (2) Qi,j,k + ∆tL(Qi,j,k ) . = Qi,j,k + 3 3

(n+1)

Qi,j,k

,I

(2)

Qi,j,k =

nc .

Qi,j,k = Qi,j,k + ∆tL(Qi,j,k ),

x + iy = r exp (i θ) ,

rc = min r(θ), θ

π 3π 6θ6 . 2 2

is

whose radius is found from the relation

he rs

The domain of stability has the form of a circle of radius rc (the inequality ∆t < rc should be satisfied for the difference scheme to be stable). On a complex plane, the domain of stability is presented as a circle

en

ce

Pu bl

In particular, rc = 1.25 for the three-step Runge–Kutta method. An advantage of the Runge–Kutta method is that it ensures positiveness of the difference scheme. If the solution and the operator L(Q) are positive at the time tn , they also remain positive at the time tn+1 . The flux vector is split into the inviscid and viscous components. The use of centered difference schemes for discretization of convective terms at high Reynolds numbers leads to computational instability and non-physical oscillations of the solution. The method of piecewise-parabolic reconstruction and the Chakravarthy–Osher scheme are used for discretization of inviscid fluxes, and the 2nd order centered finite-difference schemes are used for discretization of viscous fluxes. The inviscid flux is found from the relation i 1h F (QL , QR ) = F (QL ) + F (QR ) − |A| (QR − QL ) . 2

N

ov

a

Sc i

where the subscripts L and R refer to cells on the left and on the right edges of the control volume. The matrix A is presented in the form A = R|Λ|L, where Λ is the diagonal matrix composed from the Jacobian eigenvalues, and R and L are the matrices composed from its right and left eigenvectors, respectively. The Chakravarthy–Osher scheme implies a piecewise-parabolic distribution of the flow variables within the control volume. At an intermediate step, the auxiliary variables are determined as i αi1,m+1/2 = lm+1/2 (Qm − Qm−1 ) , i αi2,m+1/2 = lm+1/2 (Qm+1 − Qm ) , i αi3,m+1/2 = lm+1/2 (Qm+2 − Qm+1 ) .

Here, l = {l1 , l2 , . . .} is the set of left eigenvectors. The superscript i corresponds to the ith eigenvalue and to the ith eigenvector. The subscripts 1, 2 and 3 are used to numerate the auxiliary quantities.

Large-Eddy Simulation of Turbulence-induced Aero-optical Effects ...

49

X 1 + κ

 1 − κ ei ri − α e2,m+1/2 λi+ m+1/2 m+1/2 4 4 i  X 1 + κ i 1−κ i e α e2,m+1/2 + α e3,m+1/2 λi+ − ri , m+1/2 m+1/2 4 4 i

Fm+1/2 = Hm+1/2 +

α ei2,m+1/2 +

nc .

The family of schemes satisfying the TVD condition has the form

he rs

,I

where λ = {λ1 , λ2 , . . .} and r = {r1 , r2 , . . .} are the set of eigenvalues and the set of right eigenvectors, respectively. The first term in the right hand side determines the flux by the first-order scheme i 1h Hm+1/2 = F (Qm+1/2 ) + F (Qm ) − 2  1 X  i+ i − αi2,m+1/2 rm+1/2 . λm+1/2 − λi− m+1/2 2 i

Pu bl

is

The coefficients of the finite difference scheme are found from the relations n o i e α e1,m+1/2 = minmod αi1,m+1/2 , bαi2,m+1/2 , n o α ei2,m+1/2 = minmod αi2,m+1/2 , bαi1,m+1/2 , n o i e α e2,m+1/2 = minmod αi2,m+1/2 , bαi3,m+1/2 , n o α ei3,m+1/2 = minmod αi3,m+1/2 , bαi2,m+1/2 .

ce

Here, b = (3 − κ)/(1 − κ). The parameter κ determines schemes with different levels of accuracy. The flux limiter, minmod, has the form

en

minmod(x, y) = sign(x) max {0, min [ |x|, y sign(x)]} .

Sc i

The time step is found from the estimation of the inviscid and viscous fluxes   1 1 1 β = max , , ∆ti CFL ∆t1i ∆t2i

N

ov

a

where CFL is the Courant–Friedrichs–Lewy number, and β ∼ 0.5. The time step ∆t1i is calculated on the basis of the spectral radius of the Jacobian of the discrete inviscid operator, and the time step ∆t2i is found on the basis of the quasi-linear form of viscous fluxes written in primitive variables and the spectral radius of the Jacobian of the discrete viscous operator. The computational procedure is implemented as a sequence of the following steps. 1. Reconstruction of the solution in each control volume and extrapolation of the unknowns to find the state of the flow on the control-volume faces from values prescribed in the center. 2. Solving the Riemann problem for each face of the control volume with allowance for the local flow direction (normal to the control-volume face). 3. Evolution of the time step.

50

K. Volkov

σρ2 =

N i2 1 Xh ρ(x, y, z, ti) − hρ(x, y, z)i , N i=1

nc .

The root-mean square value of density fluctuations is found by processing the calculated results by the formula

11. Results and Discussion

he rs

,I

where N is the number of time steps, and the angle brackets mean time-averaging flow variables. To obtain a statistically reliable average pattern, 5×104 time steps are performed. The computational procedure is implemented as a computer code in C/C++ programming language. Parallelization of the computational procedure is performed by a message passing interface (MPI).

Pu bl

is

The results of numerical simulation of turbulent flows in flat plate boundary layer, free mixing layer and free round turbulent jet provide general statistics of optical distortions by examining the time-averaged intensity pattern and instantaneous flowfield. The flowinduced optical aberrations are studied for different inflow conditions.

111. Boundary Layer on a Flat Plate

N

ov

a

Sc i

en

ce

A coherent beam propagating through a randomly inhomogeneous medium (e.g. turbulent boundary layer or shear flow) feels optical aberrations. While the outer region of the boundary layer produces relatively small aberrations of a regular character, processes in the inner region of the boundary layer determine the electromagnetic wave propagation. Optical aberrations in the outer region of the boundary layer have a regular character and are relatively small. The x coordinate is directed along the plate, and the y coordinate is normal to it. The origin is located on the leading edge of the flat plate. The velocity profile specified on the inlet boundary of the computational domain is obtained by solving the Blasius problem with small random perturbations imposed onto this profile (white noise) [40]. No-slip and no-penetration boundary conditions are specified on the wall. The plate has a constant temperature Tw = 300 K. Non-reflecting boundary conditions are specified on the boundaries through which the gas leaves the computational domain. Periodic boundary conditions are specified in z-direction. The computations are performed at Reδ = 2.6 × 104 (Reynolds number based on the boundary-layer thickness, δ = 18 mm) and Reδm = 2.8 × 103 (Reynolds number based on the momentum thickness, δm = 2.2 mm). + The mesh contains 180 × 90 × 90 cells. For this mesh, x+ = 45, ymin = 1, z + = 14, where      1/2 uτ µw ∂ hui + y = y, uτ = . νw ρw ∂y w The time step is ∆t = 2.1 × 10−5 s, and 40 000 time steps are performed to obtain time-averaged pattern of the flow.

Large-Eddy Simulation of Turbulence-induced Aero-optical Effects ...

u+

nc .

20

51

16

,I

12

8

he rs

LES Reihardt law Linear profile Log profile

4

0 0 10

2

1

10

10

is

y+

Pu bl

Figure 4. Velocity profile in the boundary layer compared to Reihardt law.

Sc i

en

ce

The velocity profile in the near-wall region of the boundary layer based on LES calculations is in good agreement with the Reihardt law [40] based on the experimental data and covering the viscous sub-layer, the buffer region and logarithmic region of the boundary layer (Figure 4). The resulting density fluctuations are by 10–20% higher than the experimental values, presented in [52] at M = 0.88 (Figure 5). The maximum of the y axis distribution of the density fluctuations is located well close to the wall. At 0.2 < y/δ < 0.6, the root-meansquare density remains nearly constant, and the dispersion of phase fluctuations is virtually uniform in this region. Comparison of the length scales lu and lρ shows that they are equal only in the interval 0.1 < y/δ < 0.22. In the remaining part of the boundary layer, the correlation scale of the velocity fluctuations is considerably smaller than that of the density fluctuations. Using the vortex model, assuming that the beam passes through the vortex center, and integrating with respect to r from −∞ to +∞, the following relation is obtained G ∆ϕ = −2π γ



U a∞

2

ρ0 R.

(17)

N

ov

a

The subscript 0 refers to parameters at the stagnation point. An uniform chain of vortices located on one line and moving with a constant velocity is considered. Each vortex inserts optical perturbations whose maximum intensity is determined in accordance with the equation (17). Assuming these perturbations to have an approximately sinusoidal character, the dependence of the intensity of optical aberrations on the vortex density, velocity, intensity, and vortex radius is obtained σϕ ≈ 4.48Gθ2 M2∞ ρ0 R,

(18)

52

K. Volkov

0.20

ρ', kg/m 3

0.15

nc .

2

,I

1

1 LES 2 Calculations, Tromeur et al (2006) Experiment, Tromeur et al (2006)

0 0

0.3

0.6

is

0.005

he rs

0.01

0.9

1.2

1.5

Pu bl

y/δ

Figure 5. Density fluctuations in the boundary layer.

ce

where θ = U/U∞ is the turbulence intensity, M∞ is the free-stream Mach number. Under the assumption that R ∼ δ ∗ /2 (where δ ∗ is the momentum thickness) and that the turbulence intensity is 10%, the equation (18) yields σϕ ≈ 1.48 × 10−5 Φ,

a

Sc i

en

where Φ = δ ∗ ρ0 M2∞ /ρ∞ . According to the estimation of [23] performed on the basis of these measurements, the constant factor in the equation (19) is 2.4 × 10−5 . The equation (19) yields a linear dependence of the optical path length of the beam on dynamic pressure, which agrees with the data of [17], but contradicts the data of [20], which 1/2 state that σϕ ∼ p0 . The results of calculations concerning the boundary layer on a flat plate are shown in the Figure 6 (it is assumed that α = 2). The plate has a constant temperature of Tw = 300 K. The dispersions predicted with the equation (8), developed in [52] (line 4), are overestimated as compared to those calculated with the equation (6) at A = 0.1 and A = 0.2, developed in [47] (line 2 and line 3). Unlike the predictions based on correlation, developed in [52] (line 4), the dispersion of the wave phase provided by the LES calculations (line 1) has an inflection point. Following the Reynolds analogy, it is assumed that the correlation factor between the fluctuations of velocity and temperature is Ru0 T 0 = −1. The results computed for large eddies indicate that in the interval 0.12 < y/δ < 0.86 the correlation factor is nearly constant, Ru0 T 0 = −0.5, which is in line with the data of [52]. Figure 7 shows the results calculated for the boundary layer on a flat plate by the equa-

ov

N

(19)

Large-Eddy Simulation of Turbulence-induced Aero-optical Effects ...

0.024

1 2 3 4

LES Sutton (1985), A = 0.1 Sutton (1985), A = 0.2 Tromeur et al (2006)

nc .

σϕ2 , µm 2

4

,I

0.032

53

he rs

0.016 3

1

0.008

00

0.6

y/δ

0.9

1.2

1.5

Pu bl

0.3

is

2

Figure 6. Dispersion of phase fluctuations in the boundary layer.

ce

tion (19) (solid line), the data of [23] (symbols ◦ and ), and the LES results for x = 0.11 (symbols •) and 0.16 (symbols ). It is seen that numerical calculations yield lower values of the dispersion of the wave-front phase function than the results obtained in [23] for both values of the x coordinate.

en

112. Free Mixing Layer

Sc i

Two semi-infinite gas flows move in one direction with velocities u1 and u2 along the plane x < 0, y = 0 (Figure 8). The flows contact each other at the point x = 0 and further on (with x > 0) the boundary between them is turbulised. The flow in the free mixing layer is characterized by the effective Mach number Mc =

u1 − u2 . a1 + a2

N

ov

a

Momentum thickness is

1 δm (x) = (∆u)2

+∞ Z ρ (u − u1 ) (u2 − u) dy.

−∞

Reynolds number based on velocity difference and Reynolds number based on effective velocity are as follows Rem =

∆uδm0 , ν

Rec =

uc δm0 . ν

54

K. Volkov

nc .

σϕ , µm

0.04 0.03

0.01

0.5

1

1.5

2

2.5

Pu bl

Φ, µm

is

0 0

he rs

Correlation Gordeyev et al (2003), x = 0.11 cm Gordeyev et al (2003), x = 0.16 cm LES, x = 0.11 cm LES, x = 0.16 cm

0.02

,I

0.05

Figure 7. Root-mean-square values of optical aberrations in the boundary layer at x = 0.11 cm and x = 0.16 cm. Vorticity thickness is

umax − umin . max(∂u/∂y)

ce

δω (x) =

y

Reynolds number based on vorticity thickness is

en

∆uδω0 Reω = , ν

∆u ∂u/∂y



,

x=0

a

Sc i

where δω0 is the vorticity thickness at x = 0. The velocities of mixing flows are u1 = 78 m/s and u2 = 200 m/s. The working fluid is air. For these parameters, δω0 = 0.021 m and Reω = 21600, which corresponds to fully developed turbulent flow. The computations are performed on a 250 × 80 × 80 mesh at Rem = 2 × 105 (the Reynolds number is based on the momentum thickness) within the range Mc =0.15–0.8. At Mc = 0.15, compressibility has no effect on the properties of the flow (ρ2 /ρ1 = 1), while at Mc = 0.9 the density ratio is ρ2 /ρ1 = 4. The computational domain has a length 25δω0 and width 6δω0 (the y coordinate varies in the interval −3δω0 < y < +3δω0 ). The thickness of the computational domain is 5δω0 (periodic boundary conditions are used in z direction). The mesh is refined to the edge of the plate, near the flow separation plane and by the remote boundaries of the computational domain. The coordinates of nodes in y direction

ov

N

δω0 =



Large-Eddy Simulation of Turbulence-induced Aero-optical Effects ...

55

y u1

,I

nc .

u1

x

he rs

O

u2

is

u2

Pu bl

Figure 8. Free mixing layer.

are obtained from the correlation [31] (  ∆y0 aj−1 − 1 /(a − 1) for j ≤ jp , y=  yp + ∆yp bj−jp /(b − 1) for j > jp ,

ce

where ∆yp = ∆y0 ajp −1 (j = 1, 2, . . ., N ). The node with yp coordinate is located on the boundary between the internal and external sub-regions of the mixing layer, and its index is jp = 1 + log [1 + yp (a − 1)/∆y0 ] / log a.

en

The constants, a and b, determine the closeness of mesh nodes with y < yp and y > yp respectively, and they are ∆yj+1 /∆yj . The maximum aspect ratio of cell is

Sc i

Am =

min

j=0,...,N

∆yj

=

∆yN . ∆y0

a

In computations, ∆x = 0.25, yp /Ly = 0.125, ∆y0 = 0.12, a = 1.0, b = 1.02 and Am = 4.8. The stretch function is applied to the nodes located above the flow separation plane, up to the upper boundary of the computational domain. The node distribution obtained is then mirror-reflected from the flow separation plane. The time step is ∆t = 1.5 × 10−3 s, and 50 000 time steps are performed to obtain time-averaged pattern of flowfield. Free shear flows are unstable, and oscillation in these flows arises in the absence of external sources of disturbances. The profile of streamwise velocity is defined as [31]   1 1 y ∆u u(y) = (u1 + u2 ) + (u1 − u2 ) tanh ζ, δm (0) = , 2 2 2δm (0) (du/dy)y=0

ov

N

max ∆yj

j=0,...,N

56

K. Volkov

v = va sin ωt,

he rs

,I

nc .

where δm (0) is momentum thickness at the point x = 0, and ζ is a random number chosen from the normal distribution of probability with a zero mathematical expectation and a unit variance. In order to specify boundary conditions on the inlet boundary, the profiles of velocity and other flow variables are be defined, which were obtained from the calculation of turbulent boundary layer on a flat plate (with the same mesh resolution). This way seems quite expensive from the computational point of view, and other methods of formulation of unsteady boundary conditions on the inlet boundary are applied. Boundary conditions at the inlet are specified using results for turbulent boundary layer on a flat plate. The flow interface is disturbed at random with given amplitude. The transverse velocity is changed according to the law

N

ov

a

Sc i

en

ce

Pu bl

is

where va is the amplitude, and ω is the frequency (va ∼ 2 m/s and ω ∼ 400 rad/s in calculations). Non-reflecting boundary conditions are applied to the outlet boundary. Slip boundary conditions are specified on the remote boundaries of fluid domain in y direction. Periodic boundary conditions are specified in the z direction. Momentum thickness decreases with increasing effective Mach number as shown in the Figure 9. Symbols • correspond to the computational data, symbols  correspond to data of [37], symbols M, O and ♦ correspond to data of [15], and symbols ◦ correspond to data of [21]. The profiles of the density fluctuations in the mixing layer are shown in the Figure 10 at Mc = 0.8. It is seen that the distributions of the density and pressure fluctuations in the cross sections of the mixing layer are similar with a maximum at the line separating the mixing flows. The profiles of pressure fluctuations are filled to a greater extent, and the maximal amplitude of the pressure fluctuations is close to a linear dependence on the streamwise coordinate. Small deviations from the linear dependence occur only at x/L > 0.68. The dependence of the maximal density fluctuations on the x coordinate is nonmonotonic (Figure 11). At x/L < 0.6, it is close to a linear one, and shows a peak and a smooth decay. The calculated results show that lρ ∼ 4lu , where lu = 0.2k3/2/ε. Because of the Kelvin–Helmholtz instability of the free mixing layer, large-scale coherent vortex structures are formed. The downstream increase in the size of these structures is close to a linear dependence. Formation of such structures leads to the emergence of optical aberrations, which follow a sinusoidal dependence over the space and time coordinates [42]. For isentropic flows, the change in density is directly proportional to the change in pressure and inversely proportional to the change in temperature. Assuming that the change in pressure in the mixing layer is approximately proportional to the squared characteristic velocity Uc = (u2 + u1 )/2, the density ratio is ∆ρ 1 ∆p 1 ρ0 2 = ≈ U = M2c , ρ0 γ p0 γ p0 c

where Mc = (u1 −u2 )/(a1 +a2 ) is the effective Mach number, and u1 and u2 are the velocities of the mixing flows. The root-mean-square value of the optical aberrations is related

Large-Eddy Simulation of Turbulence-induced Aero-optical Effects ...

δm/δm0

nc .

1

57

,I

0.8

he rs

0.6

LES Goebel et al (1991) Elliott et al (1990) Papamoschou et al (1988)

0.2

0 0

0.4

0.6

0.8

Pu bl

0.2

is

0.4

1

Mc

Figure 9. Dependence of momentum thickness on effective Mach number.

ce

to the integral of the change in density. Choosing the vorticity thickness as a characteristic scale, we obtain σϕ ≈ G∆ρδω ≈ Gρ0 M2c δω .

(20)

In equation (20), the vorticity thickness is found as [37]

en

dδω (1 − λ)(1 + s1/2 ) , = 0.17 dx 1 + λs1/2

z1

N

ov

a

Sc i

where λ = u2 /u1 , s = ρ2 /ρ1 . Taking into account that δω ∼ x, a linear dependence of the optical aberrations on the coordinate, σϕ ∼ ρ0 M2c x, is obtained. Phase perturbations in the case of the wave passing through a layer of thickness L = |z2 − z1 | are found by integrating the distribution of the refractive index over the layer thickness Zz2 ∆ϕ(x, y, t) = k ∆n(x, y, z, t)dz, where ∆n(x, y, z, t) is the change in the refractive index along the beam path. At small perturbations of the wave vector, the integral takes the form Zz2 ϕ(x, y, t) = k0 n(x, y, z, t)dz. (21) z1

58

K. Volkov

y/δ x/L 1 0.25 2 0.5 3 0.75

,I

0.12

he rs

0 1

3

0

0.006

is

2

-0.12

-0.24

nc .

0.24

0.012

0.018

0.024

0.03

Pu bl

ρ' max, kg/m3

Figure 10. Density fluctuations in the cross section of the mixing layer. Using the equation (21), the calculated results are presented in the form [14]

ce

ϕ(x, y, t) 1 ϕ(x, e y, t) = ' k0 L∆n L∆n

Zz2h

z1

i n(x, y, z, t) − n∞ dz.

(22)

Sc i

en

It follows from the equation (22) that the following relation is valid in the direction of beam propagation ∂ϕ e 1 ψ(x) = L = ∂x ∆n

Zz2

z1

∂n dz, ∂x

where L is the size of the turbulent mixing zone. The spectra follow the dependence

a

For the mixing layer, it is assumed that ∆n = |n1 − n2 | and L = 2δ, where δ is the local width of the mixing zone, the subscripts 1 and 2 refer to the mixing flows. The results presented in the Figure 12 are calculated for the mixing layer by the semiempirical dependence (20) at Mc = 0.24 (u1 = 260 m/s and M1 = 0.77, u2 = 0.06 m/s and M2 = 0.06, in this case dδω /dx ≈ 0.25) are in reasonable agreement with the data of the physical experiment [42] (symbols ◦ correspond to A = 20 cm, and symbols  correspond to A = 30 cm) in the case of rather large aperture sizes (A > 20 cm) and with the LES results (symbols •).

ov

N

Sψ (kx L) ∼ (kxL)2 Sϕ (kx L).

Large-Eddy Simulation of Turbulence-induced Aero-optical Effects ...

ρ'max , kg/m3

nc .

0.032

59

,I

0.024

he rs

0.016

0.008

0

0.2

0.4

x/L

0.6

0.8

1.0

Pu bl

0

is

LES

Figure 11. Dependence of the maximal density fluctuations in the mixing layer on streamwise coordinate.

Sc i

en

ce

Dashed line in the Figure 12 shows the results obtained in [23] with the use of a sinusoidal law for the deflection angle α = sin(2πf t), where f = Uc /Λ is the frequency of vortex formation, Λ is the vortex size, and Uc is the characteristic velocity. Such an approach allows obtaining the final dependence of the magnitude of optical aberrations on the streamwise coordinate. The agreement between the predicted and measured data, however, becomes worse in the downstream direction. In the intermediate range of wave numbers, the spectrum in the mixing layer shows the power dependence Sϕ (kxL) ∼ (kx L)−q ,

a

where q ∼ 2 (Figure 13). The compressibility only has a weak influence on the spectrum behavior (this influence is largely displayed at large wave numbers). The density fluctuations grow with an increase in the effective Mach number. Calculations give lρ ∼ 4lu , where lu = 0.2k3/2/ε.

N

ov

113. Free Turbulent Jet

A free turbulent jet is considered, which flows out of a round nozzle into submerged space. The origin of point is located at the nozzle outlet. The positive reading of x coordinate is made in the direction of propagation of the jet. The radius of nozzle, ra, is the characteristic scale for variables with the dimensionality of length. The velocity ua and

60

K. Volkov

nc .

σ , µm 0.3 ϕ

2

0 0

10

20

30

40

Pu bl

x, cm

he rs

1

is

0.1

Correlation Gordeyev et al (2003) LES Siegenthaler et al (2005), A = 20 cm Siegenthaler et al (2005), A = 30 cm

,I

0.2

Figure 12. Root-mean-square values of optical aberrations in the mixing layer.

10 -2



Mc 1 0.2 2 0.8

en

10 -4

ce

2

10 -6

Sc i

1

N

ov

a

10 -8

10 -1 0 10

101

kx L

10 2

10 3

Figure 13. Spectrum of phase fluctuations in the free mixing layer.

Large-Eddy Simulation of Turbulence-induced Aero-optical Effects ...

61

ce

Pu bl

is

he rs

,I

nc .

temperature Ta at the nozzle outlet are the characteristic scales for variables with the dimensionality of velocity and temperature. The temperature of surrounding fluid is T∞ . The flow in the jet is characterized by the preheating parameter, ϑ = Ta/T∞ (for isothermal jet, ϑ = 1) and by the intensity of turbulence, θ = u0 /ua , where u0 is the fluctuation velocity. The calculations are performed for ra = 5 mm, ua = 80 m/s, ρa = 0.58–1.26 kg/m3 , Ta = 280–600 K, ρ∞ = 0.58–1.26 kg/m3 , T∞ = 280–600 K. The parameters at the nozzle outlet correspond to Reynolds number Re = 1.2 × 105 , which is maintained constant with relevant variation of dynamic viscosity, and to the range of variation of preheating parameter ϑ = 0.48–2.15. The calculations are performed in the domain [0, Lx] × [−Ly , Ly ] × [−Lz , Lz ], where the length of domain is Lx = 100ra, and its width and height in the inlet and outlet sections are Ly = Lz = 10ra and Ly = Lz = 40ra. The mesh contains 350 × 150 × 150 nodes. In the initial region of the jet, the mesh is uniform up to x ∼ 10ra, then mesh is refined, so ∆xmin = 0.08ra and ∆xmax = 0.15ra. In the cross section, the mesh nodes are concentrated towards the nozzle outlet, so ∆ymin = ∆zmin = 0.03ra and ∆ymax = ∆zmax = 0.09ra. The time step is ∆t = 0.08ra/ua = 5.8 × 10−5 s. A statistically reliable averaged pattern of flow is obtained by making 105 time steps. Correct specification of boundary conditions requires the calculation of flow in the pipe and boundary layer on the external surface of the nozzle [57] (a part of the nozzle is included in the computational domain). No calculation of flow in pipe is performed, and a velocity profile is preassigned at the nozzle outlet boundary at |r| ≤ ra, on which random sinusoidal oscillations are imposed [3]    ua 0.5 − |r| vx (r, t) = [1 + α sin(Sh t)] , 1 + tanh 3 2δ

en

where δ is the momentum thickness, α is the amplitude of disturbances, and Sh is the Strouhal number. It is assumed that δ/ra ∼ 0.1, α = 2.1 × 10−3 , and Sh = 0.42. Small random oscillations are further imposed on the radial distribution of circumferential velocity   vθ (r, t) = 0.025 exp −3(1 − |r|)2 φ,

N

ov

a

Sc i

where φ is a random number from uniform distribution in [−0.5, 0.5]. The radial velocity at the nozzle exit is taken to be zero vr (r, t) = 0. Boundary conditions at a large distance from the jet flowing into the submerged space (at the lower and upper boundaries) depend on the ejection properties of the jet. The boundary conditions away from a jet flowing out into submerged space are defined by the ejection properties of the jet, and an induced potential flow directed towards the jet exists away from the jet exists. The properties and parameters of this flow are not known in advance and are defined by the jet proper. The calculations of steady-state flows reveal that the best results are produced by the boundary conditions based on the exact solution which describes a potential flow outside of a round turbulent jet. Non-reflecting boundary conditions are used on the outflow boundary. Large-scale eddy structures present in the shear layer and shown in the Figure 14 have the form of toroidal axially symmetric eddies merging at a certain distance apart from the nozzle outlet (of the order of one to two section diameters).

a)

K. Volkov Midsection

b)

x/ra = 10

c)

x/ra = 80

,I

Figure 14. Contours of vorticity magnitude at time moment t = 2.32 s.

nc .

62

N

ov

a

Sc i

en

ce

Pu bl

is

he rs

Eddy structures existing in the initial part of the jet are fairly small. The typical eddy size gradually increases from the start point downstream, intensifying the exchange of momentum between the layer and the surrounding fluid. The ellipsoidal shape of the coherent structure points to the anisotropy of turbulent pulsations in the region where large-scale eddies exist. The generation of eddies is due to the Kelvin–Helmholtz instability of the shear layer. The maximum and minimum of the vorticity approximately correspond to the eddy centers. At low Reynolds numbers, Re ∼ 103 , the jet at the nozzle outlet is nearly axially symmetrical. As the Reynolds number increases with distance from the nozzle up to Re ∼ 104 , a weak sinusoidal mode appears. The results computed demonstrate that the velocity and temperature profiles in the flow cross sections exhibit a typical jet pattern, as well as a kink in the middle part. The profiles become wider with increasing axial coordinate, which points to the increase in the thickness of the zone of jet mixing with the surrounding fluid. Near the jet boundary, the profiles become sloping. As the degree of preheating of the jet increases, the maximum of fluctuations of axial velocity shifts towards the nozzle outlet. A further increase in the degree of preheating causes a qualitative variation of the pattern of this impact, and the maximum of fluctuations of axial velocity decreases. A similar impact is made by the non-isothermality of flow on the distribution of intensity of temperature fluctuations, the maximum of which is of the order of 13% of the temperature difference at the nozzle outlet. The maximum of fluctuations of axial velocity occurs at x/ra ∼ 16. In the region of the core of constant velocity, the intensity of turbulence, while remaining almost constant in the transverse direction, increases from 1.5–2% at the nozzle outlet to a value of the order of 6% at the end of the initial region with increasing distance from the nozzle. Because no turbulent mixing is present within the jet core, the effect of the increase in the degree of turbulence is associated with the penetration of pressure fluctuations from the mixing zone into the core. The fluctuations of velocity on the jet axis reach a maximum at a distance equal to two lengths of initial region. This value remains almost constant for both isothermal and non-isothermal jets and depends relatively weakly on the conditions of outflow. The radial coordinate, at which the correlation moments of velocity and temperature reach a maximum, almost coincide in sections x/ra = 10 and x/ra = 30, and correspond to the position of maximum of turbulent kinetic energy (r/ru = 0.8 according to the data of [10] and r/ru = 0.7 according to the data of [12]). In the near field of jet at x/ra ∼ 10, the heat transfer occurs along axial coordinate, because hvx0 T 0 i / hvr0 T 0 i ∼ 2 at 0.2 < r/ru < 1.4. Downstream at x/ra ∼ 30, the

Large-Eddy Simulation of Turbulence-induced Aero-optical Effects ...

63

Sϕ (kxL) ∼ (2rakx )−q ,

he rs

,I

nc .

contributions made by the axial and radial heat fluxes become approximately equal, so that hvx0 T 0 i / hvr0 T 0 i ∼ 1.1 at 0.5 < r/ru < 1.8. In a plane jet, hvx0 T 0 i / hvr0 T 0 i ∼ 2 at x/ra = 80 and 0.4 < r/ru < 1.3 [1] (Re = 7.9 × 103 and ∆T = 25 K). In an axisymmetric jet, the measurements of [1] give hvx0 T 0 i / hvr0 T 0 i ∼ 1 downstream of section x/ra ∼ 100 and at 0.3 < r/ru < 1.3. The profiles of correlation moments of density and velocity along the jet axis, which characterize the anisotropy of mass flow, are in qualitative agreement with the data of [19] based on k–ε turbulence model which does not include the anisotropy of turbulent fluctuations of velocity but enables one to estimate the anisotropy of diffusion flow of mass. The anisotropy is observed at 4 < x/ra < 40, and hρ0 vr0 i / hρ0 vx0 i ∼ 1 at x/ra > 40 downstream. In the free jet at Re = (2–8) × 104 and high wave numbers, 2rakx > 1, the wave spectrum is described by a

Re 1 5.0×103 2 1.5×104

1

en

10-2



ce

100

Pu bl

is

where q ∼ 2.5 (Figure 15). The spectrum for the jet is steeper than that for the mixing layer because of a stronger turbulent mixing, and is almost independent on the Reynolds number at high wavenumbers (2rakx > 1, where ra is the nozzle-exit radius). The values calculated for the jet are normalized to ∆n = |na − n∞ | and L = 2ra, where the subscripts a and ∞ refer to the parameters at the nozzle outlet and in submerged space.

2

Sc i

10-4

N

ov

a

10-6

10-8 -1 10

100

kx ra

101

102

Figure 15. Spectrum of phase fluctuations in the round turbulent jet.

64

K. Volkov

12. Conclusion

Sc i

en

ce

Pu bl

is

he rs

,I

nc .

Laser beams are widely used in communication and optical measurement systems. Optical effects are produced by density and temperature fluctuations which result in fluctuations of the refractive index of medium. Laser beams experience fluctuations of the beam size, beam position and intensity distribution within the beam due to refractive turbulence. Computational modeling of flows is vital for understanding and controlling processes relevant to aeronautical, environmental and numerous engineering problems. Most flows are turbulent, this dynamic complexity generally being treated using inexpensive deterministic mathematical models (e.g. RANS) requiring calibration and validation against measurements. The increasing power of computers is allowing much less model dependent approaches, in particular LES. An analysis of results calculated by solving RANS equations and a comparison of these results with LES data show that semi-empirical models of turbulence underpredict the dispersion of fluctuations of the wave phase, which can differ from LES data severalfold (because of uncertainty of constant factors involved into the corresponding relations). Semiempirical relations derived from the vortex model of the flow ensure reasonably accurate calculations of the root-mean-square values of optical aberrations in the boundary layer on a flat plate and in the free mixing layer. The dependence of the root-mean-square fluctuations of the wave phase on the dynamic pressure is linear for the boundary layer on a flat plate. For the mixing layer, this dependence is linear beginning from the coordinate counted along the flow direction. The results obtained provide information for designing and applying adaptive-optic equipment and techniques, and they are beneficial in design and optimization of optical systems, where the control of laser beam propagation is an important problem, e.g. non-intrusive measurement methods (laser two focus, laser Doppler anemometry, Doppler global velocimetry, particle image velocimetry, interferometry, laser-induced fluorescence), and coherent optical adaptive equipment and techniques (interferometric techniques, shadowgraph systems, optical windows, airborne telescopes, ground-based astronomical systems) based on the optical detection of light-scattering properties of the flow. Mitigation strategy is to use some sort of shear-layer control to regularize the vortex roll-up, thereby making it amenable to some form of adaptive-optic correction.

References

N

ov

a

[1] Antonia R.A., Prabhu A., Stephenson S.E. Conditionally sampled measurements in a heated turbulent jet. Journal of Fluid Mechanics, 1975, 72, pp. 455–480. [2] Ayyalasomayajula H., Arunajatesan S., Kannepalli C., Sinha N. Large eddy simulation of a supersonic flow over a backward-facing step for aero-optical analysis. AIAA Paper 2006-1416. [3] Boersma B.J., Lele S.K. Large eddy simulation of compressible turbulent jets. Annual Research Brief, Center for Turbulence Research, 1999, pp. 365–377.

Large-Eddy Simulation of Turbulence-induced Aero-optical Effects ...

65

[4] Carroll B.F., Boulos E., Sytsma M., Cattafesta L.N., Hubner J.P., Sheplak M. Aerooptic measurement facility characterization. AIAA Paper 2004-0936.

nc .

[5] Catrakis H.J., Aguirre R.C. New interfacial fluid thickness approach in aero-optics with applications to compressible turbulence. AIAA Journal, 2004, 42(10), pp. 1973– 1981.

,I

[6] Catrakis H.J., Aguirre R.C., Nathman J.C., Garcia P.J. Large-scale refractive turbulent interfaces and aero-optical interactions in high Reynolds number compressible separated shear layers. Journal of Turbulence, 2006, 7(54), pp. 1–21.

he rs

[7] Chakravarthy S.R., Osher S. A new class of high-accuracy TVD schemes for hyperbolic conservation laws. AIAA Paper 85-0363. [8] Chernov L.A. Wave propagation in a random medium. New York, McGraw-Hill, 1960.

Pu bl

is

[9] Cheung K., Jumper E.J. High temporal bandwidth optical wavefront sensor technologies. Proceedings of 2004 AFOSR Contractors Meeting in Turbulence and Rotating Flows, 4–5 August 2004, Denver, Colorado, USA. [10] Chevray R., Tutu N.K. Intermittency and preferential transport of heat in a round jet. Journal of Fluid Mechanics, 1978, 88, pp. 133–160. [11] Childs R.E. Prediction and control of turbulent aero-optical distortion using large eddy simulation. AIAA Paper 93-2670.

ce

[12] Chua L.P., Antonia R.A. Turbulent Prandtl number in a circular jet. International Journal of Heat and Mass Transfer, 1990, 33(2), pp. 331–339.

en

[13] Dimotakis P.E., Catrakis H.J., Fourguette D.C. Beam propagation and wavefrontphase integrals in high Reynolds number shear layers and jets. AIAA Paper 98-2833.

Sc i

[14] Dimotakis P.E., Catrakis H.J., Fourguette D.C. Flow structure and optical beam propagation in high-Reynolds number gas-phase shear layer and jets. Journal of Fluid Mechanics, 2001, 433, pp. 105–134. [15] Elliott G.S., Samimy M. Compressibility effects in free shear layers. Physics of Fluids, 1990, 2(7), pp. 1231–1240.

a

[16] Fitzgerald E.J., Jumper E.J. Further consideration of compressibility effects on shearlayer optical distortion. AIAA Paper 99-3617.

N

ov

[17] Fitzgerald E.J., Jumper E.J. Scaling aero-optic aberrations produced by high-subsonicMach shear layer. AIAA Journal, 2002, 40(7), pp. 1373–1381.

[18] Fitzgerald E.J., Jumper E.J. The optical distortion mechanism in a nearly incompressible free shear layer. Journal of Fluid Mechanics, 2004, 512, pp. 153–189.

66

K. Volkov

nc .

[19] Gazzah M.H., Sassi M., Sarh B., Gokalp I. Numerical simulation of variable density subsonic turbulent jets by using the k–ε model. International Journal of Thermal Sciences, 2002, 41(1), pp. 51–62. [20] Gilbert K.G. KC-135 aero-optical boundary-layer/shear-layer experiments. In: AeroOptical Phenomena, Eds K.G. Gilbert and L.J. Otten. New York, AIAA, 1982, 80, pp. 306–324.

,I

[21] Goebel S.G., Dutton J.C. Experimental study of compressible turbulent mixing layers. AIAA Journal, 1991, 29(4), pp. 538–546.

he rs

[22] Gordeyev S., Jumper E.J. Aero-optical and hot-wire measurements of the flow around the hemispherical turret with a flat window. AIAA Paper 2004-2450. [23] Gordeyev S., Jumper E.J. Aero-optical characteristics of compressible, subsonic turbulent boundary layers. AIAA Paper 2003-3606.

is

[24] Gordeyev S., Jumper E. Fluid dynamics and aero-optics of turrets. Progress in Aerospace Sciences, 2010, 46(8), pp. 388–400.

Pu bl

[25] Gordeyev S., Jumper E.J. The optical environment of a cylindrical turret with a flat window and the impact of passive control devices. AIAA Paper 2005-4657. [26] Gordeyev S., Jumper E.J., Ng T.T., Cain A.B. Aero-optical characteristics of compressible, subsonic turbulent boundary layers. AIAA Paper 2003-3606.

ce

[27] Huang P.G., Coleman G.N., Bradshaw P. Compressible turbulent channel flows: DNS results and modeling. Journal of Fluid Mechanics, 1995, 305, pp. 185–218.

en

[28] Hugo R.J., Jumper E.J. Experimental measurement of a time-varying optical path difference by the small-aperture beam technique. Applied Optics, 1996, 35(22), pp. 4436–4447. [29] Jones M.I., Bender E.E. CFD-based computer simulation of optical turbulence through aircraft flowfields and wakes. AIAA Paper 2001-2798.

Sc i

[30] Jumper E.J., Fitzgerald E.J. Recent advances in aero-optics. Progress in Aerospace Sciences, 2001, 37(3), pp. 299–339.

[31] Li Q., Fu S. Numerical simulation of high-speed planar mixing layer. Computers and Fluids, 2003, 32(10), pp. 1357–1377.

N

ov

a

[32] Lifshitz Y., Degani D., Tumin A. On the interaction of turbulent shear layers with harmonic perturbations. Flow, Turbulence and Combustion, 2008, 80(1), pp. 61–80. [33] Liou T., Lien W., Hwang P. Compressibility effects and mixing enhancement in turbulent free shear flows. AIAA Journal, 1995, 33(12), pp. 2332–2338. [34] Mani A., Wangy M., Moin P. Computational study of aero-optical distortions by a turbulent wake. Annual Research Brief, Center for Turbulence Research, 2006, pp. 79–91.

Large-Eddy Simulation of Turbulence-induced Aero-optical Effects ...

67

nc .

[35] Martin M.P., Piomelli U., Candler G.V. Subgrid-scale models for compressible largeeddy simulations. Theoretical and Computational Fluid Dynamics, 2000, 13(5), pp. 361–376. [36] Osher S. Riemann solvers, the entropy condition, and difference approximation. SIAM Journal of Numerical Analysis, 1984, 21(2), pp. 217–235.

,I

[37] Papamoschou D., Roshko A. The compressible turbulent shear layer: an experimental study. Journal of Fluid Mechanics, 1988, 197, pp. 453–477.

he rs

[38] Peters B., Brown D., Cole T., Haight J., Havener G. CFI-aero-optic images: issues for wave optics modeling. AIAA Paper 94-2622.

[39] Sandham N.D., Reynolds W.C. Three-dimensional simulations of large eddies in the compressible mixing layer. Journal of Fluid Mechanics, 1991, 224, pp. 133–158. [40] Schlichting H., Gersten K. Boundary layer theory. Springer Verlag, Berlin, 2000.

Pu bl

is

[41] Shen S., Ding F., Han J., Lin Y.-L., Arya S.P., Proctor F.H. Numerical modeling studies of wake vortices: real case simulation. AIAA Paper 99-0755. [42] Siegenthaler J.P., Gordeyev S., Jumper E. Shear layers and aperture effects for aerooptics. AIAA Paper 2005-4772.

[43] Siegenthaler J.P., Jumper E.J., Gordeyev S. Atmospheric propagation vs. aero-optics. AIAA Paper 2008-1076.

ce

[44] Sinha N., Arunajatesan S., Seiner J.M., Ukeiley L.S. Large-eddy simulations of aerooptic flow fields and control application. IAA Paper 2004-2448.

en

[45] Smagorinsky J. General circulation experiments with the primitive equations. Monthly Weather Review, 1963, 91(3), pp. 99–164.

Sc i

[46] Smith R., Truman C. Prediction of optical phase degradation using a turbulent transport equation for the variance of index-of-refraction fluctuations. AIAA Paper 900250. [47] Sutton G.W. Aero-optical foundations and applications. AIAA Journal, 1985, 23(10), pp. 1525–1537.

a

[48] Sutton G.W. Effect of inhomogeneous turbulence on imaging through turbulent layers. Applied Optics, 1994, 33(18), pp. 3972–3976.

N

ov

[49] Sutton G.W. Effect of turbulence fluctuations in an optically active fluid medium. AIAA Journal, 1969, 7(9), pp. 1737–1743. [50] Tatarski V.I. Wave propagation in a turbulent medium. New York, McGraw-Hill, 1961. [51] Trolinger J.D., Rose W.C. Technique for simulating and evaluating aero-optical effects in optical systems. AIAA Paper 2004-471.

68

K. Volkov

nc .

[52] Tromeur E., Garnier E., Sagaut P. Analysis of the Sutton model for aero-optical properties of compressible boundary layers. Journal of Fluids Engineering, 2006, 128(2), pp. 239–246. [53] Tromeur E., Garnier E., Sagaut P. Large-eddy simulation of aero-optical effects in a spatially developing turbulent boundary layer. Journal of Turbulence, 2006, 7(paper 1), pp. 1–28.

he rs

,I

[54] Tromeur E., Garnier E., Sagaut P., Basdevant C. Large eddy simulations of aerooptical effects in a turbulent boundary layer. Journal of Turbulence, 2003, 4(paper 5), pp. 1–16.

[55] Truman C.R. The influence of turbulent structure on optical phase distortion through turbulent shear flows, AIAA Paper 92-2817.

is

[56] Truman C.R., Lee M.J. Effects of organized turbulence structures on the phase distortion in a coherent optical beam propagating through a turbulent shear flow. Physics of Fluids, 1990, 2(5), pp. 851–857.

Pu bl

[57] Tucker P.G. Novel MILES computations for jet flows and noise. International Journal of Heat and Fluid Flow, 2004, 25(4), pp. 625–635. [58] Volkov K. Aero-optic effects in turbulent wall-bounded and free shear flows. Engineering Letters, 2010, 18(3), pp. 285–290.

ce

[59] Volkov K.N. Effect of turbulence on propagation of a coherent beam in the boundary layer and mixing layer. Journal of Applied Mechanics and Technical Physics, 2010, 51(6), pp. 827–838.

en

[60] Volkov K.N. Large-eddy simulation of free shear and wall-bounded turbulent flows. In: Atmospheric Turbulence, Meteorological Modelling and Aerodynamics, Eds. P.R. Lang and F.S. Lombargo. USA, Nova Science, 2010, pp. 505–574.

Sc i

[61] Volkov K.N. Multigrid techniques as applied to gas dynamic simulation on unstructured meshes. Computational Mathematics and Mathematical Physics, 2010, 50(11), pp. 1837–1850. [62] Volkov K.N. Preconditioning of the Euler and Navier–Stokes equations in lowvelocity flow simulation on unstructured grids. Computational Mathematics and Mathematical Physics, 2009, 49(10), pp. 1789–1803.

N

ov

a

[63] Volkov K.N. Unstructured-grid finite-volume discretization of the Navier–Stokes equations based on high-resolution difference schemes. Computational Mathematics and Mathematical Physics, 2008, 48(7), pp. 1181–1202.

[64] Volkov K.N., Emelyanov V.N. Aero-optic effects in turbulent flow and their simulation. Technical Physics, 2008, 53(2), pp. 217–223. [65] Vreman B., Geurts B., Kuerten H. Large-eddy simulation of the turbulent mixing layer. Journal of Fluid Mechanics, 1997, 339, pp. 357–390.

Large-Eddy Simulation of Turbulence-induced Aero-optical Effects ...

69

[66] Wang T., Zhao Y., Xu D. Numerical study of evaluation optical quality of supersonic flow fields. Applied Optics, 2007, 46(23), pp. 5545–5551.

nc .

[67] Wei H., Chen C.P. A nonequilibrium algebraic model for turbulent density fluctuations. International Journal of Heat and Mass Transfer, 1996, 39(18), pp. 3989–3991.

,I

[68] White M.D. High-order parabolic beam approximation for aero-optics. Journal of Computational Physics, 2010, 229(15), pp. 5465–5485.

he rs

[69] Wu L., Fang J., Yang Z., Wu S. Study on a neural network model for high speed turbulent boundary layer inducing optical distortions. International Journal for Light and Electron Optics, 2011, 122(17), pp. 1572–1575.

[70] Wyckham C.M., Zaidi S.H., Miles R.B., Smits A.J. Characterization of optical wavefront distortions due to a boundary layer at hypersonic speeds. AIAA Paper 2003-4308.

is

[71] Xu L., Cai Y. Imaging deviation through non-uniform flow fields around high-speed flying vehicles. International Journal for Light and Electron Optics, 2012, 123(13), pp. 1177–1182.

Pu bl

[72] Yakhot A., Orszag S.A., Yakhot V., Israeli M. Renormalization group formulation of large-eddy simulation. Journal of Scientific Computing, 1986, 1, pp. 1–51. [73] Zhang Y.-P., Fan Z.-G. Study on the optical path difference of aero-optical window. International Journal for Light and Electron Optics, 2007, 118(12), pp. 557–560.

ce

[74] Zubair F.R., Catrakis H.J. Aero-optical interactions along laser beam propagation paths in compressible turbulence. AIAA Journal, 2007, 45(7), pp. 1663–1674.

N

ov

a

Sc i

en

Reviewed by Professor Vladislav N. Emelyanov, Head of Department of Gas and Plasma Dynamics, Faculty of Power Engineering, Baltic State Technical University, ul. 1-ya Krasnoarmeyskaya, 1, 190005 Saint Petersburg, Russia, e-mail: [email protected].

a

ov

N ce

en

Sc i he rs

is

Pu bl

,I

nc .

nc .

In: Turbulent Flows: Prediction, Modeling and Analysis ISBN: 978-1-62417-742-2 Editor: Zied Driss © 2013 Nova Science Publishers, Inc.

,I

Chapter 3

he rs

DEVELOPMENT OF A NEW CURVATURE LAW OF THE WALL FOR INTERNAL SWIRLING AXIAL FLOWS Xiuhua A. Si1, and Jinxiang Xi2, 1

is

Calvin College, Grand Rapids, MI, US Central Michigan University, Mount Pleasant, MI, US

Pu bl

2

ABSTRACT

N

ov

a

Sc i

en

ce

There has been considerable controversy during the past decade concerning the validity of the classical log law that describes the mean-velocity profile in the overlap region of turbulent wall-bounded flows. In this chapter, a new law of the wall accounting for curvature effects for internal swirling axial flows was derived and validated again benchmark experiments. The influence of the curvature on the turbulence mixing lengths in both axial and tangential directions was analyzed theoretically based on the Reynolds stress transport equations. For equilibrium flows with weak curvature, identical mixing lengths were derived for the axial and tangential directions. Additionally, the effect of finite local curvature and shear stress ratio on the near-wall velocities was systematically examined. It was observed that the curvature effect in swirling axial flows is suppressed by a factor of 1/(1+σw2)compared to that in curved channel flows, where σw is the ratio of the axial to swirl shear stress. For a given curvature radius, the maximum velocity deviation occurs when the axial-to-swirl shear stress ratio is zero. Finally, the performance of the new curvature law was evaluated by implementing it as a wall function in a well-established CFD code. The new wall function provides improved agreement for swirl velocity distributions inside labyrinth cavities in comparison with existing experimental laser Doppler anemometry measurements.

Keywords: turbulent wall-bounded flows, streamline curvature, law of the wall, overlap region, swirling axial flow

 

Department of Engineering, Calvin College, 3201 Burton SE, Grand Rapids, MI 49546, US. Department of Mechanical Engineering, Central Michigan Univresity, 2425 S. Mission Ave, Mount Pleasant, MI 48859, US. E-mail: [email protected].

72

Xiuhua A. Si and Jinxiang Xi

is

Greek symbols. α δ δij ε ε κ σw ν τ φ

en

ce

Pu bl

empirical parameter in mixing length correction. turbulent boundary layer thickness. Kronecker Delta. kinetic energy dissipation rate. perturbation parameter. von Kármán constant. axial-to-swirl wall shear stress ratio. kinematic viscosity. shear stress. energy redistribution.

Superscripts

fluctuating quantity. dimensionless value in wall coordinate.

Sc i

' +

1. INTRODUCTION

a

The universal log law of wall has been derived independently by Prandtl and Taylor using mixing length arguments, by von K´arm´an using dimensional reasoning, and by Millikan using asymptotic analysis [1]. Even though this law of wall has been playing a vital role in turbulent modeling, its validity has been frequently questioned when the scenario is skewed away from the low-Reynolds-plate-flows such as a flow with high-Reynolds number, rotation, high pressure gradient, or dramatic heat and mass transfer. In turbomachinery industries, swirling axial flow along curved solid walls is a fundamental phenomenon; the key question hereof is what the curvature effects are on the mean and statistical turbulence quantities. In the past decades, this question has been experimentally studied by Bradshaw [1, 2], Patel et al. [3], Meroney et al. [4], So and Mellor [5,6], Hunt et al. [7], Gillis and

ov

N

,I

turbulence kinetic energy. turbulence mixing length in axial direction. turbulence mixing length in circumferential direction. kinetic energy production. total wall shear velocity. radius of the wall curvature. Reynolds number (=2 uin Cr /ν ). Richardson number. velocity in axial direction [m/s]. wall shear velocity in axial direction [m/s]. velocity in circumferential direction [m/s]. wall shear velocity in circumferential direction [m/s]. coordinate normal to the curved wall.

he rs

k lrx lrθ P qτ R Re Ri u uτ w uτ y

nc .

NOMENCLATURE

Development of a New Curvature Law of the Wall for Internal Swirling Axial Flows 73 Johnston [8], Muck et al. [9], Hoffmann et al. [10], Gibson, et al. [11], and Barlow et al. [12]. As a result of their efforts, it is now generally acknowledged that:

,I

nc .

(i) Curvature influences both mean flow and turbulence structures; (ii) Even small curvature could induce large changes in Reynolds stresses; and (iii) Turbulent shear stresses are dampened on convex surfaces and amplified on concave surfaces.

N

ov

a

Sc i

en

ce

Pu bl

is

he rs

Theoretical efforts to understand the effects of curvature on the law of the wall have been undertaken by Bradshaw [1,2,13], So [14–16], Moser and Moin [17], Patel and Sotiropoulos [18], among others. The analogy between the curvature effect and buoyancy was first recognized by Prandtl in 1929 and was later employed in a model by Bradshaw [1]. In this model, the curvature effect was described for two-dimensional flows by defining a dimensionless parameter known as the gradient Richardson number Ri. Patel [3] studied secondary flows in a curved channel and showed a similar effect between curvature and pressure gradients, with the influence of convex curvature being similar to that of an adverse pressure gradient. In contrast, the effect of concave curvature resembled that of a favorable pressure gradient. In a similar manner, permeable walls [19], heat flux [20] and fluid compressibility [21] can all similarly alter the near-wall velocity gradients, shifting the velocity profiles above the classical log-law on a wall with blowing or cooling (in liquids), while shifting the profile beneath the classical log-law on a wall with suction, heating (in liquids), or compressibility. This is schematically illustrated in Figure 1. Despite its inadequacies, the classical logarithmic law was still widely used for curved walls by both experimentalists and modelers with the qualification that its validity region is closer to the wall than that on a flat plate. Attempts to account for the curvature effects in turbulent modeling were performed by Wilcox and Chambers [22], who adjusted for curvature by adding extra source terms in the k–ω model equations. Rodi [23] investigated the curved shear layers using two-equation turbulence models and showed a considerable influence on the solution accuracy by employing a curvature wall function. Galperin and Kantha [24] used an algebraic Reynolds stress model to obtain a curvature law by assuming equilibrium and constant flux inside the boundary layer.

Figure 1. Deviation of the universal law of wall.

74

Xiuhua A. Si and Jinxiang Xi

N

ov

a

Sc i

en

ce

Pu bl

is

he rs

,I

nc .

Assuming a non-constant-flux shear layer, Kind et al. [25] proposed a correction to the axial and tangential mixing length components based on the shear stress variations in the direction normal to the wall. However, their mixing length correction neglected the direct effect of curvature and assumed it had only an indirect effect on near-wall velocity distributions. This was improved later by Kim and Rhode [26] who proposed a curvature correction to the circumferential component of the mixing length using the local gradient Richardson number [1]. In Kim and Rhode’s work, the streamline curvature effect was limited only to the circumferential direction. A curvature effect was not assumed in the axial direction, a fact which is not supported by experimental observations. The goal of this paper is to assess curvature effects in both the axial and circumferential directions by means of a theoretical analysis using the Reynolds stress model and to propose a generalized curvature law for swirling axial flows based on this analysis. In this chapter, a new law of the wall accounting for curvature effects in swirling axial flows is derived. The curvature parameter that enters this problem is y/(y+R), where y is the normal coordinate, and R is the radius of curvature. The difficulty lies in the subtle way in which curvature influences the mean flow and the turbulence structures. Kind et al. [25] suggested that the near-wall flow was determined by wall shear stress and that curvature only indirectly affected the flow structure through its influence on the wall shear stress. Schwarz et al. [32] attempted to explain the curvature effects in terms of the centrifugal forces and the counteracting cross-stream pressure gradients. Experiments have shown large changes in Reynolds stresses even on small curvatures. However, these changes cannot be inferred from the relative magnitude of the curvature related terms y/(y+R) in the mean-flow equations, which is often an order of magnitude less than the resultant shear stress variation. Bradshaw [2] attributed this apparent incapacity of the mean-flow equations to high-order fluctuating interaction terms such as pressure redistribution and turbulence diffusion. Because of this, Patel et al. [18] suggested that even the Reynolds stress model requires further development of the high-order interaction terms in order to capture the curvature effects.

Figure 2. Typical curvature flows.

Development of a New Curvature Law of the Wall for Internal Swirling Axial Flows 75

nc .

The present study aims to better describe the near-wall behaviors of the turbulent swirling axial flows due to wall curvature in rotating scenarios. The goal is to improve the prediction accuracy of the swirl velocities in labyrinth seals by formulating a new curvature wall function. This is significant because the swirl variation along the rotor surface has a major influence on rotordynamic instabilities [32–36].

,I

2. METHODS

he rs

2. 1. Curvature Law Derivation



u 'i u 'k u j , k  u ' j u 'k ui ,k  -

q'  1 q '2 2 2 u ' u '  q '  ij ij  i j  3lt  3   3

ce



Pu bl

is

It is generally accepted that streamline curvature has a significant effect on flow properties in the near-wall region for 2-D curved channel flows as illustrated in Figure 2a. These results, however, cannot be readily extended to a three-dimensional swirling axial flow model because of the additional degree of freedom, i.e., the axial direction. Specifically, the influence of the circumferential curvature on tangential shear stress in swirling axial flows was found to be different from that in curved channel flows [16, 26]. Furthermore, behavior of the axial shear stresses under the influence of circumferential curvature has rarely received serious attention and remains ambiguous. This section theoretically examines whether circumferential curvature changes the shear stress in the axial direction, and how this change further alters the flow conditions in the circumferential direction. We begin with the Reynolds stress transport equations. Based on the assumption of local equilibrium and small advection and diffusion, the Reynolds stress transport equation inside a turbulent boundary layer is formulated as follows [16], (1)

where q '  u 'i u 'i  2 , lt is a turbulent length scale, and Λ is the Kolmogorov scale. The terms

en

1

Sc i

in Eq. (1) account for, respectively, production, energy redistribution, and dissipation. Only the turbulence return-to-isotropy term was retained in the redistribution term. The mean-strain rapid-distortion term was neglected because of its small effect in shear flows under small external body forces [33]. Expanding Eq. (1) in cylindrical coordinates, the turbulent shear stresses in both axial and circumferential directions were derived:

 u   w w   u     uv   lrx2         rx  r   r r   r 

(2)

 u   w w   w w     vw   lr2          r  r   r r   r r 

(3)

2

2

N

ov

a



and 2



2

76

Xiuhua A. Si and Jinxiang Xi

3 2  lr   y 1   Ri    Ri   2  

1

2

(4)

1

2

(5)

he rs

lrx  w/ r   1   Ri  lr  w / r 

,I

and

nc .

where lrx and lrθ are the turbulent mixing lengths in the axial and tangential directions, respectively.

u / y w / y  w / r

; Ri 

w / r  w / r  w / r 

 u / r 

2

  w / r - w / r 

2

(6)

ce



Pu bl

is

In the above derivation, the only approximation introduced is the assumption that w / r  w / r , which is generally true for both weak curvatures (e.g., δ/R3), the anchor impeller becomes deformed, and the velocity field is preserved.

N

ov

a

Sc i

en

ce

Pu bl

4. 1. 2. Flow Patterns in the R-Z Plane Figures 7, 8 and 9 show the secondary flow (U,W) in the different r-z planes defined by the angular coordinates equal to θ=35°, θ=45° and θ=55° respectively. In the upstream of the anchor impeller (Figure 7), the hydrodynamic flow was characterized by the presence of a circulation zone near the bottom of the tank. Near the walls, it created an ascending vertical jet that moves towards the second circulation zone. This jet was also transformed into a descending jet near to the axis of the anchor impeller. If the anchor geometry undergoes a deformation due to the flexion of the arm, the circulation zones would change their position from one iteration to an other.

Figure 6. Flow patterns induced in the r-θ plane defined by z=1.

Pu bl

is

he rs

,I

nc .

Computer Simulation of Turbulent Flow Generated by a Deformed Anchor Impeller 129

N

ov

a

Sc i

en

ce

Figure 7. Flow patterns induced in the r-z plane defined by θ=35°.

Figure 8. Flow patterns induced in the r-z plane defined by θ=45°.

S. Karray, Z. Driss, H. Kchaou et al.

Pu bl

is

he rs

,I

nc .

130

Figure 9. Flow patterns induced in the r-z plane defined by θ=55°. 2nd iteration

4th iteration

N

ov

a

Sc i

en

ce

Initial

(a) Radial position r=0.5

(b) Radial position r=0.8

(c) Radial position r=0.9

Figure 10. Axial profiles of the radial velocity.

Moreover, the centres of the lower circulation zones move towards the walls of the tank due to the flexion effect of the arm. At the same time, the centres of the higher circulation

Computer Simulation of Turbulent Flow Generated by a Deformed Anchor Impeller 131

ce

Pu bl

is

he rs

,I

nc .

zones come near to the anchor impeller axis. The positions of these centres were stabilized in the second iteration (Figure 7.c). The plane defined by the angular position θ=45° corresponds to the plane which coincides with the anchor impeller plane in its initial state (Figure 8.a). This figure shows correctly the no sliding condition which assumes that the velocity field is null. In the following iterations, this zone disappears due to the deformation of the anchor impeller arm. Also, we observe a strong penetration of the fluid towards the tank interior (Figure 8). Moreover, the phenomenon of brewing is clearly observed with the flexion of the arm. This phenomenon is very clear during the second iteration. This result is already expected, since the anchor impeller arms undergo a large deformation. In the downstream of the anchor impeller (Figure 9), the velocity field converges towards the axis of the anchor impeller because of the depression which is located behind the anchor impeller arms.

N

ov

a

Sc i

en

Figure 11. Turbulent kinetic energy k in the r-θ plane defined by z=1.

Figure 12. Turbulent kinetic energy k in the r-θ plane defined by z=1.5.

132

S. Karray, Z. Driss, H. Kchaou et al.

nc .

The brewing phenomenon was detected along the arm. This phenomenon was developed into two opposite directions. The first moves towards the walls of the tank while the second moves towards the anchor axis. In the third iteration, the flow has only a lower circulation zone (Figures 9.d, 9.e and 9.f).

he rs

,I

4. 1. 3. Axial Profiles of the Radial Velocity Component Figure 10 illustrates the predicted axial profiles of the dimensionless radial velocity component U(z). The dimensionless radial coordinates used are equal to r=0.5, r=0.8 and r=0.9. The superimposed profiles correspond to the second and the fourth iteration of the coupling algorithm. These profiles show a parabolic shape function which develops with the level of the horizontal arm of the anchor impeller. The extrema of these functions is respectively defined by the axial positions equal to z=0.33, z=0.37 and z=0.35. During these various iterations, the intensity of the radial jet decreases with the increase of the horizontal arm retreated. This result is already confirmed within a retreated blade paddle [28].

N

ov

a

Sc i

en

ce

Pu bl

is

4. 1. 4. Distribution of the Turbulent Kinetic Energy in the R-Θ Plane Figures 11, 12 and 13 presented respectively the evolution of the distribution of the turbulent kinetic energy k in the horizontal planes located in the middle height of the tank, in the top and the bottom of the anchor impeller. These planes are respectively defined by the axial positions equal to z=1, z=1.5 and z=0.16. These results show that the area defined by the maximum values is located in the wake which develops at the arm end of the anchor. These maximum values are reached during the first iteration and are about k=0.35, k=0.25 and k=0.1 respectively. In addition, the wake zone is extended at the initial state in the mid-height of the tank (Figure 11.a). In these conditions, the angular sector is equal to αin=90°. This sector decreases from one iteration to another. Indeed, for the first iteration, the value of the angular sector is equal to α1it=45° (Figure 11.b).

Figure 13. Turbulent kinetic energy k in the r-θ plane defined by z=0.16.

is

he rs

,I

nc .

Computer Simulation of Turbulent Flow Generated by a Deformed Anchor Impeller 133

a

Sc i

en

ce

Pu bl

Figure 14. Turbulent kinetic energy k in the r-z plane defined by θ=30°.

N

ov

Figure 15. Turbulent kinetic energy k in the r-z plane defined by θ=45°.

This value decreases further during the various iterations, and stabilizes at the third iteration. In these conditions, the angular sector is equal to α3it=23° (Figures 11.d, 11.e, 11.f). In the top of the anchor impeller (Figure 12), it’s noted a progress in the training of the area defined by the maximum values. Indeed, in the initial state (Figure 12.a), we observe a

134

S. Karray, Z. Driss, H. Kchaou et al.

Pu bl

is

he rs

,I

nc .

complete disappearance of this area. However, we find in the first iteration a more reduced shape which widened more and more during the other iterations (Figure 12.b). In the entire domain swept by the vertical anchor arms, the turbulent kinetic energy remains high. Outside this domain, the turbulent kinetic energy quickly becomes very weak. In the anchor impeller bottom, it is shown that through the field swept by the vertical arms of the anchor, the turbulent kinetic energy remains rather high. Outside this area, the turbulent kinetic energy becomes gradually very weak (Figure 13).

N

ov

a

Sc i

en

ce

Figure 16. Turbulent kinetic energy k in the r-z plane defined by θ=45°.

Figure 17. Dissipation rate of the turbulent kinetic energy ε in the r-θ plane defined by z=1.

Computer Simulation of Turbulent Flow Generated by a Deformed Anchor Impeller 135

Pu bl

is

he rs

,I

nc .

4. 1. 5. Distribution of the Turbulent Kinetic Energy in the r-Z Plane Figure 14, 15 and 16 present the distribution of the turbulent kinetic energy k in the r-z plane located in the upstream, at the level and in the downstream of the anchor impeller. These planes are defined by the angular positions equal to θ=30°, θ=45° and θ=60°. These results show that the turbulent kinetic energy is high between the walls and the anchor impeller. Outside this area, k decreases. Therefore, we can easily note that these results roughly reproduce the geometrical shape of the impeller. Indeed, in the upstream of the anchor impeller, the maximal value of the turbulent kinetic energy is obtained during the first iteration (Figure 14.b). Also, it’s noted that the end of the vertical arm deformed reaches the representation plane. This is displayed by the presence of a small zone, where the turbulent kinetic energy of the fluid is null. This zone changes shape during the iterations. At the level of the anchor impeller, where the arms of the anchor have undergone a deformation due to the flexion, the maximum value changes during these various iterations (Figure 15). Indeed, it is noted that the zone, where the turbulent kinetic energy of the fluid is null, change after the first iteration. These results are obviously due to the flexion of the arm. In the downstream of the anchor impeller, we registered an important turbulent kinetic energy production because of the inhalation of the fluid behind arms. This production decreases during the iterations (Figure 16).

Sc i

en

ce

4. 1. 6. Distribution of Dissipation Rate of the Turbulent Kinetic Energy in the r-Θ Plane Figures 17, 18 and 19 show the distribution of the dissipation rate of the turbulent kinetic energy ε in the horizontal planes located in the mid-height of the tank, in the top and in the bottom of the anchor impeller during the iterations of the coupling algorithm. These planes correspond to various axial positions defined respectively by z=1, z=1.5 and z=0.16. According to these results, we observe a distribution similar to the turbulent kinetic energy. Indeed, it is noted that the area of maximum values is located in the wake which develops near the end of the arm. However, the dissipation rate of the turbulent kinetic energy becomes very weak outside of the field swept by the arm. The greatest rate is reached in the fourth and the fifth iteration and it is equal to =0.6 (Figure 17). In the both planes located just in the bottom (z=0.16) and in the top (z=1.5) of the anchor impeller, a weakening of the dissipation rate of the turbulent kinetic energy takes place. In the same way, the maximum value of ε decreases from iteration to another (Figure 18 and 19).

N

ov

a

4. 1. 7. Distribution of Dissipation Rate of the Turbulent Kinetic Energy in the r-z Plane Figures 20, 21 and 22 show the distribution of the dissipation rate of the turbulent kinetic energy in the r-z planes for an angular positions equal to θ=30°, θ=45° and θ=60°. These planes are located in the upstream, at the level and in the downstream of the anchor impeller. These results show that the dissipation rate of the turbulent kinetic energy is high between the walls and the anchor impeller. It undergoes a very fast drop outside this field. During these iterations, the maximum value changes due to the flexion of the anchor impeller arm. In upstream of the anchor impeller, the maximal value of the dissipation rate of the turbulent kinetic energy is obtained during the first iteration and is equal to 1.1 (Figure 20.b). Also, the deformed arms reach the representation plane.

136

S. Karray, Z. Driss, H. Kchaou et al.

ce

Pu bl

is

he rs

,I

nc .

This is visualized by the presence of a small zone where the dissipation rate of the turbulent kinetic energy of the fluid decreases. At the level of the anchor impeller, the maximum value of the dissipation rate of the turbulent kinetic energy is reached before the third iteration and is equal to 0.3 (Figure 21.d). Indeed, it is noted that the end of the deformed arms reaches the representation plane. This result is visualized by the presence of a small zone where the dissipation rate of the turbulent kinetic energy of the fluid is null. This zone changes in the first iteration, due obviously to the flexion of the arms (Figure 21.b). In downstream of the anchor impeller, the dissipation rate of the turbulent kinetic energy is important because of the inhalation of the fluid behind blades. This production decreases during the iterations. The maximum value of the dissipation rate of the turbulent kinetic energy is obtained during the initial state and is equal 0.14 (Figure 22.a).

N

ov

a

Sc i

en

Figure 18. Dissipation rate of the turbulent kinetic energy ε in the r-θ plane defined by z=1.5.

Figure 19. Dissipation rate of the turbulent kinetic energy ε in the r-θ plane defined by z=0.16.

Pu bl

is

he rs

,I

nc .

Computer Simulation of Turbulent Flow Generated by a Deformed Anchor Impeller 137

N

ov

a

Sc i

en

ce

Figure 20. Dissipation rate of the turbulent kinetic energy ε in the r-z plane defined by θ=30°.

Figure 21. Dissipation rate of the turbulent kinetic energy ε in the r-z plane defined by θ=45°.

S. Karray, Z. Driss, H. Kchaou et al.

Pu bl

is

he rs

,I

nc .

138

N

ov

a

Sc i

en

ce

Figure 22. Dissipation rate of the turbulent kinetic energy ε in the r-z plane defined by θ=60°.

Figure 23. Turbulent viscosity in the r-θ plane defined by z=0.16.

is

Figure 24. Turbulent viscosity in the r-θ plane defined by z=1.

he rs

,I

nc .

Computer Simulation of Turbulent Flow Generated by a Deformed Anchor Impeller 139

N

ov

a

Sc i

en

ce

Pu bl

4. 1. 8. Distribution of the Turbulent Viscosity in the r-Θ Plane Figures 23 and 24 show the distribution of the turbulent viscosity in the horizontal planes located in the bottom of the anchor impeller and in the mid-height of the tank during the iterations of the coupling algorithm. These presentations planes correspond to the axial positions equal respectively to z=0.16 and z=1. In the field swept by the arms of the anchor, the turbulent viscosity remains rather high. Outside this area, it is noted that the turbulent viscosity undergoes a very fast drop (Figure 24), due to the deceleration of the flow.

Figure 25. Turbulent viscosity in the r-z plane defined by θ=30°.

S. Karray, Z. Driss, H. Kchaou et al.

is

Figure 26. Turbulent viscosity in the r-z plane defined by θ=45°.

he rs

,I

nc .

140

Pu bl

In the bottom of the anchor, the turbulent viscosity undergoes a progressive reduction outside this area swept by the arms of the anchor impeller (Figure 23). By comparing these results from iteration to another, it’s noted that the maximum value is reached during the second iteration and it is equal to υt = 1100(Figure 24.b).

N

ov

a

Sc i

en

ce

4. 1. 9. Distribution of the Turbulent Viscosity in the r-z Plane Figures 25, 26 and 27 present the distribution of the turbulent viscosity for an angular positions equal to θ=30°, θ=45° and θ=60°.

Figure 27. Turbulent viscosity in the r-z plane defined by θ=60°.

Pu bl

is

Figure 28. Radial profiles of the angular velocity.

he rs

,I

nc .

Computer Simulation of Turbulent Flow Generated by a Deformed Anchor Impeller 141

Figure 29. Meshing.

N

ov

a

Sc i

en

ce

These positions correspond to planes located in the upstream, at the level and in the downstream of the anchor impeller. In upstream of the anchor impeller, the turbulent viscosity is high between the walls and the anchor impeller (Figure 25). Outside this domain, these values decrease gradually until reach very weak values at the level of the anchor axis. Also, it’s noted that the deformed vertical arm reaches the representation plane. This is confirmed by the presence of a zone where the turbulent viscosity of the fluid is null. This zone changes during the iterations (Figure 25.b, 25.c). At the level of the anchor impeller (Figure 26), the vertical arm deformation reaches the representation plane. This is visualized by the presence of a zone where the turbulent viscosity of the fluid is null. This zone corresponds to the arm of the anchor in the initial state (Figure 26.a). Then, it decreases from iteration to another. The maximum values remain high inside the tank. They are reached in the second iteration and is equal to υt=1100 (Figure 26.c). In downstream of the anchor impeller (Figure 27), the turbulent viscosity reaches very big values from the second iteration. These values are equal to υt =1100 (Figure 27.c).

4. 1. 10. Comparison with Anterior Results The power number calculated within an anchor impeller is equal to Np=0.52. This value, obtained by a volume integration of the field of the dissipation energy, is compared to the experimental results found in the literature. In fact, the values found by Pollard and Kantyka

142

S. Karray, Z. Driss, H. Kchaou et al.

,I

nc .

[29], and Nagata [19] range between 0.5 and 0.6. The comparison of the power number with the experimental results shows a good agreement. Although in figure 28, the radial profiles of the dimensionless angular velocity component V(r) are presented in the laminar flow. In these conditions, the angular coordinate are equal to θ= 90°. The results founded by Murthy and Jayanti [30] are superposed over an average error of 8 %. The good agreement confirms the validity of the numerical method.

4. 2. Static Studies

en

ce

Pu bl

is

he rs

In the static studies, we are interested in the structural part of the mixing system. Particularly, the static behaviour of the anchor impeller and the deformation of the anchor impeller arms of under the flow effect have been studied.

N

ov

a

Sc i

Figure 30. Displacement field of the anchor impeller.

Figure 31. Points positions.

ce

Pu bl

is

Figure 32. Variation of the radial displacement Δr.

he rs

,I

nc .

Computer Simulation of Turbulent Flow Generated by a Deformed Anchor Impeller 143

en

Figure 33. Variation of the angular displacement Δθ.

N

ov

a

Sc i

The discretisation is ensured by finite elements method. In these conditions, 5163 nodes and 6580 elements are obtained. 2784 elements are quadratic type on the level of the anchor arms and 3796 elements are tetrahedral type on the anchor axis (Figure 29). Figure 30 presents the results from the CSD code. Particularly, we are interested in the study of the static behaviour of the anchor impeller and the evolution of the displacement field of the arms during various iterations of our coupling algorithm. Particularly, we are interested in the nodal displacement of the anchor arm to allow us to mesh the fluid domain. This operation is repeated several times until equilibrium is obtained. Under these conditions, it is supposed that the deformation of the anchor axis due to the torsion is negligible. To study the displacements variations, six points placed exactly in the anchor impeller are especially monitored. Their positions can be seen in figure 31. We can follow the radial and angular displacement of these points. The time history of these points is shown in figures 32 and 33. In the beginning of the simulation, all displacements reach their absolute maxima. Through the different iterations, they start oscillating around some values. Obviously, the biggest radial variation displacement Δr is at the point 1. It is fluctuating around Δr=0.03 m. On the

144

S. Karray, Z. Driss, H. Kchaou et al.

nc .

other hand, the smallest displacements are at the point 5 and is equal to Δr=0.0005 m. The angular variation displacement is biggest at the point 2 and fluctuates around Δθ=14°. Also, the value of points 5 and 6 are rather small and is equal to Δθ=3.5°.

CONCLUSION

ce

Pu bl

is

he rs

,I

In this chapter, coupling codes are used for solving fluid-structure interaction (FSI) problems in the stirred tank equipped by an anchor impeller. A computational fluid dynamics (CFD) code and a computational structural dynamics (CSD) code have been coupled by using an efficient coupling interface. Particularly, the fluid field in stirred tank was solved by the CFD code in turbulent regime. The deformation of the anchor impeller was solved by the CSD code. Specific techniques of rearrangement of grid and treatment of the boundary conditions are used to follow the behaviour of the system. This method takes advantage of the parallel process involved within each analysis code. This allows both parts of the fluidstructure interaction problem to be solved in the best possible way: a finite volume method for the fluid dynamics and a finite element method for the structure. The CFD results obtained allow a visualizing of the velocity field, the turbulent kinetic energy, the dissipation rate of the turbulent kinetic energy, the turbulent viscosity and the mechanical deformation. These results proof that the fluid flow can have a significant effect on the deformation of the anchor impeller arms. So, we can conclude that it’s fundamental to take consideration of this phenomenon in the design of the industrial process. The comparison of the power number and the velocity profile with the results founded in the literature shows a good comparison. This proves the validity of the numerical method.

REFERENCES

Sc i

[2]

F. J. Blom. A monolithical fluid–structure interaction algorithm applied to the piston problem. Computer Methods in Applied Mechanics and Engineering, 167, 369-391, 1998. E. H. Van Brummelen, S. J. Hulshoff and R. De Borst. Energy conservation under incompatibility for fluid–structure interaction problems, Computer Methods in Applied Mechanics and Engineering, 192, 2727-2748, 2003. C. Michler, S. J. Hulshoff, E. H. Van Brummelen, and R. De Borst. A monolithic approach to fluid–structure interaction, Computers and Fluids, 33, 839-848, 2004. S. Piperno and C. Farhat. Partitioned procedures for the transient solution of coupled aeroelastic problems—Part II: energy transfer analysis and three-dimensional applications, Computer Methods in Applied Mechanics and Engineering, 190, 31473170, 2001. D. C. Sternel, M. Schäfer, M. Heck and S. Yigit. Efficiency and accuracy of fluidstructure interaction simulations using an implicit partitioned approach, Comput. Mech. 43, 103-113, 2008.

en

[1]

[3]

N

ov

a

[4]

[5]

Computer Simulation of Turbulent Flow Generated by a Deformed Anchor Impeller 145

[11] [12]

[13]

[14]

[15]

[16]

nc .

,I

Sc i

[17]

he rs

[10]

is

[9]

Pu bl

[8]

ce

[7]

S. Karray, Z. Driss, H. Kchaou, M. S. Abid, Numerical simulation of fluid-structure interaction in a stirred vessel equipped with an anchor impeller, Journal of Mechanical Science and Technology, 25 (7), 1-12, 2011. C. Cebral, S. Piperno and B. Larrouturou. Partitioned procedures for the transient solution of coupled aeroelastic problems. Part 1: Model problem, theory and twodimensional application, Journal of Computer Methods in Applied Mechanics and Engineering, 124, 79-112, 1995. K. C. Park and C. Felippa. A Partitioned analysis of coupled systems. Computational Methods for Transient Analysis in Computational Methods in Mechanics, 1, 157-218, 1983. W. M. Wood. Practical Time-stepping Schemes, Oxford University Press, New York, 1990. T. J. Pedley and K. D. Stephanoff. Flow along a channel with a time-dependent indentation in one wall: the generation of vorticity waves, Journal of Fluid Mechanics, 160, 337-367, 1985. M. E. Ralph and T. J. Pedley. Viscous and inviscid flows in a channel with a moving indentation, Journal of Fluid Mechanics, 209, 543-566, 1989. S. Natarajan and M. R. Mokhtarzadeh-Dehghan. Numerical prediction of a (potential) soft acting peristaltic blood pump. International Journal for Numerical Methods in Fluids, 32, 711-724, 2000. C. Cossu and L. Morino. On the instability of a spring-mounted circular cylinder in a viscous flow at low reynolds numbers, Journal of Fluids and Structures, 14, 183-196, 2000. A. Beckert and H. Wendland. Multivariate interpolation for fluid-structure interaction problems using radial basis functions, Aerospace, Science and Technology, 5, 125-134, 2001. G. Sieber. Numerical Simulation of Fluid-Structure Interaction Using Loose Coupling Methods, PhD thesis, at the Department of Numerical Methods in Mechanical Engineering, Darmstadt University of Technology, 2002. M. Glück, M. Breuer, F. Durst, A. Hlfmann and E. Rank. Computation of wind-induced vibrations of flexible shells and membranous structures. Journal of Fluids and Structures, 17, 739-765, 2003. E. Bucchignani, F. Stella and F. Paglia. A partition method for the solution of a coupled liquid-structure interaction problem, Applied Numerical Mathematics, 51, 463-475, 2004. Y. Wang. Combination of C. F. D. and C. S. D. packages for fluid-structure interaction, Journal of Hydrodynamics, 20, 756-761, 2008. S. Nagata. Mixing principles and applications. Halstead press, Japan, 1975. M. Ammar, Z. Driss, W. Chtourou, and M. S. Abid, Effects of baffle length on turbulent flows generated in stirred vessels, Cent. Eur. J. Eng., 1, 401-412, 2011. M. Ammar, Z. Driss, W. Chtourou, M. S. Abid. Effect of the Tank Design on the Flow Pattern Generated with a Pitched Blade Turbine, International Journal of Mechanics and Applications, Vol. 2, N. 1, pp. 12-19, 2012. Z. Driss, S. Karray, H. Kchaou and M. S. Abid. C. F. D. simulation of the laminar flow in stirred tanks generated by double helical ribbons and double helical screw ribbons impellers, Cent. Eur. J. Eng., 1, 413-422, 2011.

en

[6]

[18]

a

[19] [20]

N

ov

[21]

[22]

146

S. Karray, Z. Driss, H. Kchaou et al.

N

ov

a

Sc i

en

ce

Pu bl

is

he rs

,I

nc .

[23] Z. Driss, S. Karray, H. Kchaou, M. S. Abid. Computer Simulations of Fluid-Structure Interaction Generated by a Flat-Blade Paddle in a Vessel Tank, International Review of Mechanical Engineering (I. R. E. M. E.), 1, 608-617, 2007. [24] S. Karray, Z. Driss, H. Kchaou, M. S. Abid. Hydromechanics characterization of the turbulent flow generated by anchors impellers. Engineering Applications of Computational Fluid Mechanics, 5 (3), 315-328, 2011. [25] S. V. Patankar. Numerical heat transfer and fluid flow, Series in Computational Methods in Mechanics and Thermal Sciences, McGraw Hill, New York, 1980. [26] Z. Driss, G. Bouzgarrou, W. Chtourou, H. Kchaou, M. S. Abid, Computational studies of the pitched blade turbines design effect on the stirred tank flow characteristics, European Journal of Mechanics B/Fluids, 29, 236-245, 2010. [27] W. Chtourou, M. Ammar, Z. Driss and M. S. Abid, Effect of the turbulence models on Rushton turbine generated flow in a stirred vessel, Cent. Eur. J. Eng., 2011, 1, 380389. [28] Z. Driss, H. Kchaou, M. Baccar and M. S. Abid. Numerical investigation of internal laminar flow generated by a retreated-blade paddle and a flat-blade paddle in a vessel tank, International Journal of Engineering Simulation, 6, 10-16, 2005. [29] J. Pollard and T. A. Kantyka. Heat transfer to agitated non-Newtonian Fluids, T. Instn. Chem. Eng., 47, 21-27, 1969. [30] S. S. Murthy, S. Jayanti. Mixing of power-law fluids using anchors: Metzner-Otto concept revisited, AIChE Journal, 1, 30-40, 2003.

nc .

In: Turbulent Flows: Prediction, Modeling and Analysis ISBN: 978-1-62417-742-2 Editor: Zied Driss © 2013 Nova Science Publishers, Inc.

,I

Chapter 6

is

he rs

TURBULENT FLOW STRUCTURES FOR DIFFERENT ROUGHNESS CONDITIONS OF CHANNEL WALLS: RESULTS OF EXPERIMENTAL INVESTIGATION IN LABORATORY FLUMES

Pu bl

D. Termini

University of Palermo, Palermo, Italy

ABSTRACT

a

Sc i

en

ce

The present chapter synthesizes the main findings of experimental studies conducted at DICAM in order to analyze the effect of roughness conditions of the channel walls on the velocity field and turbulent structure in a straight open channel flow. In particular, the research concerns on the dynamics of the horizontal turbulence and the coherent motion. The experiments were carried out both over flat bed and for two different roughness conditions of the side-walls and over fixed bed-forms. The measured velocity data have been used to identify the formation of the “horizontal” burst sequence along the examined channel and to compare the dynamics of turbulent events between the different considered roughness conditions of the channel’s walls. Results show that the bed forms lead to the regularization of the horizontal burst sequence; as effect of the side-walls roughness, the spatial lag of events of the same type increases moving towards the channel axis.

N

ov

Keywords: open channel flow, flow turbulence structure, burst cycle, walls roughness, laboratory experiments

Dipartimento di Ingegneria Civile, Ambientale e Aerospaziale, dei Materiali (D. I. C. A. M.), University of Palermo (Italy), Viale delle Scienze, 90128 Palermo, Italy, E-mail: [email protected].

D. Termini

148

NOMENCLATURE

– – – – – – –

u* – z – z/h –

nc .

ks ks + q RIV S S0 t

,I



he rs

k

median sediment diameter. sediment diameter – 16% finer by weight. sediment diameter – 84% finer by weight. occurrence frequency of k-th turbulent event. water depth. discriminating function. index (k = I, II, III, IV) indicating the quadrant where the generic turbulent event falls. equivalent sand roughness. no-dimensional equivalent sand roughness. threshold level of the excluding hole. ratio fIV / fII. initial section. reference section. time. measured velocity component (i = x stream–wise direction, i = y transverse direction). instantaneous turbulent fluctuation component in the i–th direction. instantaneous filtered turbulent fluctuation component in the i–th direction. shear stress velocity. distance from the bed. relative water depth.

is

– – – – –

Pu bl

D50 D16 D84 fk h



x _ –

spatial lag. time lag. water kinematic viscosity. correlation coefficient.

Sc i

en

 – xy –

ce

Greek symbols

a

As it is well-known, the generation, development and break up process of coherent structures in open-channel flows is an important mechanism which determine the main energy exchange in turbulence [1-5]. However, the investigations on the dynamics of coherent structures are still not enough. The turbulent coherent structures evolve periodically as part of the so-called burst sequence that includes four events (inward interaction, ejection, sweep, burst ) in which a small-scale eddy develops into a large-scale one [6]. The potential energy occurring during the growing of the eddy near the wall (inward interaction) is converted during the ejection event into kinetic turbulent energy which is further produced during the subsequent sweep event; the dissipation of turbulent kinetic energy occurs during the burst event so that the energy is transferred to smaller eddies in the well-known cascade of energy process [7].

ov

N

1. INTRODUCTION

Turbulent Flow Structures for Different Roughness Conditions of Channel Walls

149

Sc i

en

ce

Pu bl

is

he rs

,I

nc .

The coherent structures can be formed either by vertical eddies (which rotate in the vertical plane x, z) or by horizontal eddies (which rotate in the horizontal plane x, y) [8]. In the past, most studies on bursting phenomena were conducted in laboratory channels with smooth side-walls and bed [9-12]. But the point is that the occurrence frequency of each turbulent event might be affected by the roughness condition of the channel walls. Investigations of bursting phenomena over rough beds have been limited in comparison with those over smooth beds. Most of these studies have principally concerned well sorted bed, i.e. composed of homogenous particles. It has been verified that (see as an example in [13-15]) the effect of the bed roughness (the bed roughness is defined by the equivalent to sand roughness, ks, and it is analyzed at scale ks + = ks u*/with u*=shear velocity, Kinematic water viscosity) on turbulence intensity is evident for the relative water depth z/h < 0.30 (h=water depth). As ks+ increases the viscous sub-layer disappears and the separation of lowspeed fluid occurs in the cavity between roughness elements (see as an example [14]). Furthermore, experiments have shown that sweep motions become stronger than ejection as ks+ increases [16]. However, the mechanism of turbulence production over a rough bed has not yet been sufficiently clarified [17]. There have been few studies of turbulence in flows over rough natural gravel surfaces, which are characterized by no-regular dimension of particles on the bed (see as an example [18-19]). These studies have highlighted that vortices scaling with the grain size form in the downstream wake of protruding clusters; larger structures may be created by amalgamation of numerous smaller scale structures. Very few studies have been conducted in order to analyze the effect of the side-walls roughness. Some researchers [20-21] have highlighted the importance of the effect of the side-walls roughness in turbulence anisotropy and in determining secondary currents, especially in narrow channels [22]. Furthermore, recent experimental literature findings [8, 23, 24] suggest that the side-walls roughness influences especially the evolution of the horizontal turbulent structures. The point is that still today the existing information on the dynamics of horizontal bursts remains rather vague. In such a context, experimental program has been recently carried out in a straight laboratory flume constructed at the laboratory of Dipartimento di Ingegneria Civile, Ambientale, Aerospaziale, dei Materiali (DICAM) – University of Palermo (Italy) to give a contribution on understanding on the evolution of horizontal turbulence and the influence of the roughness of the channel’s walls. This chapter synthesizes the main results obtained until today.

2. EXPERIMENTAL APPARATUS

N

ov

a

The laboratory flume is 0.4 m wide and 7.0 m long. The bed of the channel (longitudinal slope of 0.45%), consists of quartz sand (D50 = 0.65 mm, D16= 0.55 mm, D84 = 0.90 mm). The bed is fixed in a first reach 1.0 m long and it is mobile in the remaining 6.0 m reach. The channel’s side-walls are rigid and made of Plexiglas. Three experimental runs were conducted in a turbulent subcritical flow with flow discharge of 0.013 m3/s and initial longitudinal bed slope of 0.45%: for the first run (SB run) the side-walls of the channel were smooth (Plexiglas strips) and the bed surface was flat and rigid; for the second run (RB run) rough side-walls

D. Termini

150

( zi   )( zi  k   ) ( zi   )2( zi  k   )2

i = 1, …, Nb

Pu bl

r( k ) 

is

he rs

,I

nc .

were used (the Plexiglas strips were covered by the same sand as the bed); the third run (DB run) was conducted with smooth side-walls (Plexiglas strips) as the SB run, but over fixed deformed bed. The deformed bed configuration was obtained as result of a previous mobilebed run. During the mobile-bed experiment downstream migrating bed-forms occurred. The position of these bed-forms through time was measured both by using a graduated rule fixed to the Plexiglas walls of the channel and by a PV09 profile indicator (by Delft Hydraulics - precision of 0.1 mm) mounted on a carriage running along slides parallel to the channel side-walls. At the end of the experiment, 33 longitudinal bed profiles were measured in five transverse abscissas (see Figure 1). Then, the measured data were interpolated by applying the linear interpolation method based on optimal Delaunay triangulation [25]. The instantaneous 3D surfaces were used to reproduce the bed-form evolution through time. Both by the aforementioned analysis and by direct observations during the mobile bed experiment, the bed-form wavelength was estimated of about 80 cm. Furthermore, it was verified that the positions of pools and crests of the bed-forms were alternated passing from the side-walls to the channel axis. The spatial description of bed forms evolution was also evaluated by the autocorrelation of the instantaneous bed level data in the spatial domain. The autocorrelation function was estimated as [26]:

(1)

en

ce

where zi is the bed level measured at point i, k is the autocorrelation lag, Nb is the number of data considered and  is the arithmetic mean of the data series considered. Equation (1) was applied for different values of the lag (k=0 to N-1) and assuming Nb=N (being N the number of measured data). Figure 2 shows the autocorrelation function estimated along axis 1 (i.e. near the side-wall) and axis 3 (i.e. at channel axis). It can be observed that the relative width of peaks corresponds to sections distant each other of 0.8-0.9 m. Furthermore, positive peaks at axis 3 correspond to negative peaks near the side-wall. This behavior confirms that the position of the bed-forms was alternated passing from the side-walls to the channel axis. During each experiment (SB, RB and DB - duration of about 7 hours) the instantaneous longitudinal ( u x (z,t ) - with z=distance from the bed, t=time) and transversal ( u y (z,t) )

N

ov

a

Sc i

velocity components were measured in points, spaced of 3 mm apart, along five verticals of 17 cross-sections (cross-sections A-S of Figure 1). As discussed in previous works [27, 28], data analysis has been restricted to the channel reach between sections D and H, where the bed forms wavelength was clearly distinguishable.

Figure 1. Longitudinal axes and measurement channel reach.

151

he rs

,I

nc .

Turbulent Flow Structures for Different Roughness Conditions of Channel Walls

Figure 2. Autocorrelation function.

is

For each measurement point, the instantaneous velocity components were measured by using the 2D (side-looking) Acoustic Doppler Velocimeter (ADV) (sampling frequency of 25 Hz). More details on experiments and measurement conditions can be found in [26-28].

Pu bl

3. HORIZONTAL TURBULENT STRUCTURES 3. 1. Occurrence Frequency of Events

Sc i

en

ce

In order to identify the evolution of coherent horizontal vortices, the quadrant analysis [29] has been applied. Quadrant analysis, as well as the POD (proper Orthogonal decomposition) method, is a technique frequently used in detecting bursting phenomena. The quadrant analysis is based on conditional single-point measurements [17]; the POD is based on the two-point velocity correlation [30] and, thus, requires the knowledge of a two-point correlation tensor over a great number of points. According with quadrant analysis [29], the events can be classified on the basis of the sign of the components of the turbulent fluctuations. Thus, the instantaneous turbulent fluctuation components have been estimated as deviations of the instantaneous flow velocity components ( ui (z,t) - with i = x, y) from the corresponding time-averaged velocity values. The instantaneous plane [ ux ( z,t) ; uy ( z,t) ] has been divided into four quadrants so that points falling in quadrant I correspond to inward interactions events ( ux ( z,t ) >0, uy ( z,t ) >0); points falling in quadrant II correspond to sweep events ( ux ( z,t ) >0, uy ( z,t )