**2.1.2 Formulation of problem and boundary conditions**

Inviscid incompressible fluid flow is described by one partial differential equation which satisfy law of fluid continuity (Laplace equation):

$$\frac{\partial^2 \Phi}{\partial \mathbf{x}^2} + \frac{\partial^2 \Phi}{\partial \mathbf{y}^2} = 0$$

Equation is the same for steady and unsteady flow. Hence, methods for steady flow can be applied in a solution procedure for unsteady flow. When is the panel method applied, the problem is reduced on determination of the velocity field which must satisfy boundary condition of zero normal flow at the surface of airfoil.

When the airfoil is moving relative to a steady undisturbed flow, this boundary condition must be expressed in form:

$$
\left[\vec{V}\_{\rm F}\left(t\right) - \vec{V}\_{\rm S}\left(t\right)\right] \cdot \vec{n}\left(t\right) = 0
$$

where *VF* is the fluid velocity, *VS* is the surface velocity and *<sup>n</sup>* is normal at the surface.

When the airfoil is moving relative to a steady incompressible inviscid flow, circulation (t) around airfoil varying with time t. Total circulation around airfoil and vortex wake must be zero (Kelvin theorem), which is possible to achieve only if the circulation with intensity (*/t*)*t* is separated from trailing edge of airfoil in interval dt. This separate vorticity is carried out in downstream direction. There is direct relationship between this and Kutta condition for trailing edge. Kutta condition can be expressed as equality of airfoil upper and lower side pressures on the trailing edge, and on the basis of unsteady Bernoulli equation can be expressed as:

$$\frac{\partial \Gamma}{\partial t} = -\frac{1}{2} \left( V\_{\mathcal{U}}^2 - V\_{\mathcal{L}}^2 \right) = -\frac{\left( V\_{\mathcal{U}} + V\_{\mathcal{L}} \right)}{2} \left( V\_{\mathcal{U}} - V\_{\mathcal{L}} \right)^2$$

where *U* and *L* denote upper and lower side of trailing edge, respectively ( *x l* 1).

Therefore, above mentioned consistence of Kutta condition and vortex shedding model can be clearly seen. Namely, during time above mentioned circulation (t) around airfoil is balanced by vortex shedding with intensity *V V U L* , with mean velocity *V V U L* 2 .

Therefore, unsteadiness of the problem is taken into consideration through unsteady form of kinematics boundary condition, Kelvin theorem, unsteady Bernoulli equation, and unsteady form of Kutta condition for trailing edge.

Numerical method for solution of unsteady flow over airfoil oscillating in incompressible, inviscid flow field is based on well known panel method in form developed for steady flow by Hess and Smith. Due to unsteadiness in addition to boundary conditions, additional conditions are necessary. Kelvin theorem and Kutta condition is to be applied in unsteady

Inviscid incompressible fluid flow is described by one partial differential equation which

2 2 2 2 0 *x y*

 

 

 

Equation is the same for steady and unsteady flow. Hence, methods for steady flow can be applied in a solution procedure for unsteady flow. When is the panel method applied, the problem is reduced on determination of the velocity field which must satisfy boundary

When the airfoil is moving relative to a steady undisturbed flow, this boundary condition

 <sup>0</sup> *V t V t nt F S* 

is the surface velocity and *<sup>n</sup>*

When the airfoil is moving relative to a steady incompressible inviscid flow, circulation (t) around airfoil varying with time t. Total circulation around airfoil and vortex wake must be zero (Kelvin theorem), which is possible to achieve only if the circulation with intensity

> <sup>1</sup> 2 2

Therefore, above mentioned consistence of Kutta condition and vortex shedding model can be clearly seen. Namely, during time above mentioned circulation (t) around airfoil is balanced by vortex shedding with intensity *V V U L* , with mean velocity *V V U L* 2 .

Therefore, unsteadiness of the problem is taken into consideration through unsteady form of kinematics boundary condition, Kelvin theorem, unsteady Bernoulli equation, and unsteady

2 2

where *U* and *L* denote upper and lower side of trailing edge, respectively ( *x l* 1).

*t* is separated from trailing edge of airfoil in interval dt. This separate vorticity is carried out in downstream direction. There is direct relationship between this and Kutta condition for trailing edge. Kutta condition can be expressed as equality of airfoil upper and lower side pressures on the trailing edge, and on the basis of unsteady Bernoulli equation

> *U L U L U L V V V V V V*

is normal at the surface.

**2.1.1 Introduction** 

form. This results with a nonlinear problem.

satisfy law of fluid continuity (Laplace equation):

condition of zero normal flow at the surface of airfoil.

*t* 

form of Kutta condition for trailing edge.

must be expressed in form:

is the fluid velocity, *VS*

where *VF*

can be expressed as:

(*/t*)

**2.1.2 Formulation of problem and boundary conditions** 

#### **2.1.3 Discretization and numerical solution procedure**

Solution of flow over airfoil (moving arbitrary), depending on time and starting from time t=0 is calculated in successive time intervals tk (to= 0, k=1,2,3...). Fig. 1 shows model for the time tk.

Airfoil contour in time tk is replaced by N linear elements. (i)k and k are uniform source and circulation distributions on i-th element (i=1,2,...,N), where (i)k varies from one element to another, k is the same for each element on airfoil and k denotes time tk. Total circulation k is given as <sup>k</sup> (airfoil perimeter).

An elementary vortex wake with length of k and pitch angle of k in regard to the x-axis (to the free stream direction) is attached to the trailing edge. Length k and angle k are arbitrary in first iteration. Their values will be determined as the part of the solution. Circulation in trailing edge vortex wake element is (w)k, where:

$$
\Delta\_k \left( \mathcal{Y}\_w \right)\_k = \Gamma\_k - \Gamma\_{k-1} \tag{1}
$$

Hence, circulation on the element is equal to the difference between circulations around airfoil in times tk-1 and tk, assuming that k-1 has been already determined. Vortex wake consist of concentrated vortices formed by vorticity shed at earlier times, which is assumed to be transformed into discrete vortices. Concentrated vortices is moving with resulting velocity calculated in the center of each vortex at each successive time interval. Therefore, strength and positions of discrete vortices are regarded as known at time tk, and there are N+3 unknowns (i)k (i=1,2,...,N), k, k and k at the time tk,.

Fig. 1. Solution at time tk

They are determined by satisfying following conditions:

 N conditions of zero normal velocity component in external middle point at each segment of airfoil

$$(V\_{nj})\_k = 0\tag{2}$$

Computer program was developed from formulated numerical method. Flow over NACA 0012 airfoil was modeled by using 20 linear elements. Airfoil oscillating at high frequency trough certain angle of attack is considered flow characteristics are computed in time interval from 0 to 30 sec. Fig. 2 shows static and dynamic lift curve resulting from the application of this numerical method. Dynamic lift curve is obtained by airfoil oscillating

Dynamics lift curves for two reduced frequencies k=0.5 and k=0.15 based on obtained results are also shown (Fig. 3). In this case, airfoil oscillating through angle of attach of 2.5, with amplitude of 4. Dynamic lift curve shape depending on reduced frequency can be

25 4 . sin

*t*

Fig. 4. Computed vortex traces as a function of the reduced frequency k,

through angle of attack of 6, with amplitude 5 and reduced frequency k=0.5.

**2.1.4 Discussion of results** 

Fig. 2. Static and dynamic lift curve

Fig. 3. Influence of reduced frequency

seen from the figure.

where ( ) *Vnj <sup>k</sup>* is total normal velocity component at the external middle point of j-th. element in time tk.

 condition of equal pressures in middle points of two elements on an airfoil on both sides of trailing edge (unsteady Kutta condition):

$$(V\_{t1})\_k^2 = (V\_{tN})\_k^2 + 2\left(\Gamma\_k - \Gamma\_{k-1}\right) / \left(t\_k - t\_{k-1}\right) \tag{3}$$

where <sup>1</sup> ( ) *Vt k* is total tangent velocity component at middle point of first element and ( ) *VtN k* is total tangent velocity component at middle point of N-th element in time tk.

 trailing edge vortex wake element length and direction (k and k ) are determined from condition that the element is tangent to the local resultant velocity and that its velocity is proportional to the local resulting velocity. If (Uw)k and (Ww)k are components of total velocity at middle point of trailing edge vortex wake element, excluding influence of the element on itself, then:

$$\tan \theta\_k = \left(\mathcal{W}\_w\right)\_k \Big/ \left(\mathcal{U}\_w\right)\_k \quad , \qquad \Delta\_k = \left[\left(\mathcal{U}\_w\right)\_k^2 + \left(\mathcal{W}\_w\right)\_k^2\right]^{\frac{1}{2}} \left[t\_k - t\_{k-1}\right] \tag{4}$$

By taking into the consideration that the problem considers incompressible flow, formulas for velocities induced from source and vorticity distributions are the same as these for the steady case.

Basic set of equations is nonlinear therefore following iterative procedure for its solution is accepted. Since k and k are assumed it leaves N linear equations from (2) and square equation from (3). N linear equations are solved to give (i)k (i=1,2,...,N) relatively to the k. Then k is determined from square equation (3). Once known (i)k (i=1,2,...,N) and k, (Uw)k and (Ww)k can be determined and replaced into (4) to obtain new values of k and k. Procedure is repeated until k and k are of requested accuracy.

When intensities are determined, source and vortex distribution is known. Pressure coefficient follows from unsteady form of Bernoulli equation:

$$\mathbf{c}\_p = \mathbf{l} - \frac{V^2}{\mathcal{U}\_\alpha^2} - \frac{2}{\mathcal{U}\_\alpha^2} \frac{\partial \mathcal{p}}{\partial t} \tag{5}$$

where V is total velocity on external side of airfoil, and is velocity potential. Forces and moment coefficients are obtained by direct integration of pressure distribution.

Since solution in time tk is determined, model is applied for the time tk+1 in the same way as this for time tk, with vortex vale computed in time tk. Distributed vorticity in trailing edge vortex wake element in time tk is now considered as concentrated in vortex with strength (w)kk at the position in time tk+1 determined by coordinates:

$$\mathbf{x} = \begin{pmatrix} \mathbf{x}\_{\rm TE} \end{pmatrix}\_k + \frac{1}{2} \boldsymbol{\Delta}\_k \cos \theta\_k + \begin{pmatrix} \boldsymbol{\mathcal{U}}\_w \end{pmatrix}\_k \begin{pmatrix} \mathbf{t}\_{k+1} - \mathbf{t}\_k \end{pmatrix} \quad , \quad \mathbf{z} = \begin{pmatrix} \mathbf{z}\_{\rm TE} \end{pmatrix}\_k + \frac{1}{2} \boldsymbol{\Delta}\_k \sin \theta\_k + \begin{pmatrix} \mathbf{V} \mathbf{V}\_w \end{pmatrix}\_k \begin{pmatrix} \mathbf{t}\_{k+1} - \mathbf{t}\_k \end{pmatrix}\_k$$

Resultant velocity in the center of each next concentrated vortex in vortex wake is calculated by solution at time tk, and their position in time tk+1 directly follows.

condition of equal pressures in middle points of two elements on an airfoil on both

2 2

where <sup>1</sup> ( ) *Vt k* is total tangent velocity component at middle point of first element and ( ) *VtN k* is total tangent velocity component at middle point of N-th element in time tk.

condition that the element is tangent to the local resultant velocity and that its velocity is proportional to the local resulting velocity. If (Uw)k and (Ww)k are components of total velocity at middle point of trailing edge vortex wake element, excluding influence of

<sup>1</sup> tan , *k w w k w w kk W U <sup>k</sup> <sup>k</sup> k k*

By taking into the consideration that the problem considers incompressible flow, formulas for velocities induced from source and vorticity distributions are the same as these for the

Basic set of equations is nonlinear therefore following iterative procedure for its solution is

and (Ww)k can be determined and replaced into (4) to obtain new values of k and

When intensities are determined, source and vortex distribution is known. Pressure

2 2 2 <sup>2</sup> <sup>1</sup> *<sup>p</sup> V*

where V is total velocity on external side of airfoil, and is velocity potential. Forces and

Since solution in time tk is determined, model is applied for the time tk+1 in the same way as this for time tk, with vortex vale computed in time tk. Distributed vorticity in trailing edge vortex wake element in time tk is now considered as concentrated in vortex with strength

 1 1 1 1 2 2 *TE k k w k k* cos , *TE k k w k k* sin *k k k k x x*

*U t t zz Wt t*

Resultant velocity in the center of each next concentrated vortex in vortex wake is calculated

k are of requested accuracy.

*U U t*

<sup>1</sup> 1 1 () ( ) *V V t k tN k k k k k* 2 / *t t* (3)

1

i)k (i=1,2,...,N) relatively to the

i)k (i=1,2,...,N) and

k.

k.

k, (Uw)k

2 2 2

k are assumed it leaves N linear equations from (2) and square

(5)

 

*U W tt* (4)

k ) are determined from

element in time tk.

the element on itself, then:

steady case.

Then 

(

accepted. Since k and

Procedure is repeated until k and

equation from (3). N linear equations are solved to give (

coefficient follows from unsteady form of Bernoulli equation:

w)kk at the position in time tk+1 determined by coordinates:

by solution at time tk, and their position in time tk+1 directly follows.

k is determined from square equation (3). Once known (

*c*

moment coefficients are obtained by direct integration of pressure distribution.

sides of trailing edge (unsteady Kutta condition):

trailing edge vortex wake element length and direction (k and

where ( ) *Vnj <sup>k</sup>* is total normal velocity component at the external middle point of j-th.

#### **2.1.4 Discussion of results**

Computer program was developed from formulated numerical method. Flow over NACA 0012 airfoil was modeled by using 20 linear elements. Airfoil oscillating at high frequency trough certain angle of attack is considered flow characteristics are computed in time interval from 0 to 30 sec. Fig. 2 shows static and dynamic lift curve resulting from the application of this numerical method. Dynamic lift curve is obtained by airfoil oscillating through angle of attack of 6, with amplitude 5 and reduced frequency k=0.5.

Fig. 2. Static and dynamic lift curve

Fig. 3. Influence of reduced frequency

Dynamics lift curves for two reduced frequencies k=0.5 and k=0.15 based on obtained results are also shown (Fig. 3). In this case, airfoil oscillating through angle of attach of 2.5, with amplitude of 4. Dynamic lift curve shape depending on reduced frequency can be seen from the figure.

Fig. 4. Computed vortex traces as a function of the reduced frequency k, 25 4 . sin*t*

helicopter rotor planning allows experimental investigations in tunnels and in flight to be final examinations. In that way, a very useful interaction between numerical calculation and

The finale result is obtaining of unsteady lift dependence of angle of attack. Model established in such a way is characteristic for the helicopter rotor blade airflow and it is based on the

The planar potential flow of incompressible fluid can be treated in Cartesian coordinates *x* and *y*. Two dimensional potential incompressible flow is completely defined by the speed

*xy xy*

fulfillment of Cauchy-Reimann conditions enables combining of the velocity potential and

potential flow of incompressible fluid as a function of a complex coordinate. The complex analytical function *w* (*z*), called the complex flow potential, always has a unique value for the first derivative. This derivative of the complex potential is equal to the complex velocity

> *dw u iv V dz*

where *u* is its real part (velocity component in *x*-direction) and *v* is the imaginary part

Circulation and flow are equal to zero for any closed curve in complex plane. Complex

If we assume that the moving vortex of intensity 0 is at the distance *z*0, than the perturbance

 2 2 2

*z z z z a*

<sup>1</sup> <sup>2</sup>

 

2 2

2 2 ( ) ln ln( ) ln *i i a a i i w z V ze V e i a V sin z z z <sup>z</sup>*

*i i a i <sup>i</sup> i az V V e V e a V sin*

 

2 2

*z z*

0 0 2

2

0 0 i

0 0

0 0

0

( )

2

  0

*z z*

( )

(,) (,) *y i x y* . This complex function entirely defines the planar

0 and <sup>2</sup>

0 . The

 

 

influence of the previous lifting surface's wake influence on the next coming blade.

potential and stream function and is presented by Cauchy-Reimann equations:

experimental results is achieved.

stream function *wz x*

(velocity component in *y*-direction).

**2.2.3 Simulation of the moving vortex** 

 

potential of such a flow is:

and it's complex velocity:

at that point, i.e.:

**2.2.2 Foundations of the irrotational 2-D flow** 

 

where *u* and *v* are velocities in *x* and *y* directions respectively.

potential has no singularities except at the stagnation point.

 

> 

The Laplace partial differential equations can also be introduced <sup>2</sup>

 

Vortex wake shape of airfoil oscillating through same angle of attack, with same amplitude, but with different reduced frequency obtained by this method is shown at Fig. 4. It is similar with one obtained by methods of visualization in tunnels.

Numerical method presented in this paper leads to inviscid flow field calculation over airfoil moving arbitrary in time if it is supposed that flow stays attached and that it is separated art trailing edge. Method is completely general although only results of oscillating at high frequencies are presented.
