A two-dimensional sail in turbulent flow

in

Fluid Structure Interaction

S.K. Chakrabarti; C.A. Brebbia
in Advances in Fluid Mechanics, Vol.30, pp. 245-254

WIT Press, Southampton, Boston, 2001
ISBN 1-85312-881-3

presented on Fluid Structure Interaction 2001, Chalkidiki, Greece, 2001

besser als postscript-download
dieses HTML-Dokument automatisch erstellt, also evtl. fehlerhaft


U. Bunge, T. Rung, F. Thiele
Hermann-Föttinger-Institute of Fluid Mechanics,
Technical University Berlin, Germany

Abstract

The two-dimensional, incompressible flow at Reynolds number based on chord length of $1.3 \cdot 10^6$ around an inextensible, massless sail (flexible membrane airfoil) with varying excess length is examined solving the Reynolds-Averaged Navier-Stokes (RANS) equations on computational grids deforming according to the sail movement. The grid deformation is computed using the same numerical procedure solving the basic equations of displacements for an isotropic, linearly elastic continuum. Results are presented for fully turbulent conditions employing closure models of different degree of complexity in comparison with experimental and analytical results obtained from literature. The sail is modelled as a chain with frictionless hinges. Good agreement can be found for low angle of attack and small excess length. However, for higher angle of attack, approaching onset of separation and beyond, the predictive accuracy varies significantly with the representation of turbulence in conjunction with strong unsteady phenomena.

1 Introduction

Only recently the turbulent viscous flow equations have been introduced as the fluid dynamic model in sail or membrane airfoil theory. An advanced study of flexible, massless membrane wing aerodynamics incorporating the RANS equations has been presented in Smith, Shyy, [1] Unlike an airfoil which is usually regarded as a rigid 2D structure, the sail shape is a function of both incidence and Reynolds number in which the shape itself interacts with the pressure distribution in a rather complicated manner. A hysteresis in aerodynamic characteristics is observed since there are two equilibria at small angles of attack, i.e., the sail snaps through at a certain angle of attack. However, snapping through could not yet be computed due to the fact that the sail shape becomes highly instable and the grid deformation very large.

There are several references for experimental results, e.g. Greenhalgh, Curtiss, Smith, [2], which will be taken for comparison, and a variety of theoretical investigations concerning membrane airfoil or sail theory exist, cf. [1].

The computation of turbulent flow around a flexible 2D sail poses a challenge for the flow solver as well as the structural-mechanics solver, since the simulation has to be time-accurate on a domain with moving boundaries. The prediction of the aerodynamic forces has to be as exact as possible in order to allow for a realistic representation of the sail shape that deforms according to these loads. As classical inviscid theory only holds for small ranges of angles of incidence with small excess lengths, the application of RANS methods becomes necessary for this task.

A chain model as given in de Matteis, de Socio, [3] is considered for the sail in conjunction with a RANS method instead of the sail equation that is usually used.

2 Description of Method

2.1 Fluid Mechanics

For the computation of the 2D flow around the sail an in-house finite-volume (FV) method is used, which is based on an implicit procedure of second order accuracy in time and space. Due to its fully implicit scheme, true time steps are possible and steady flow simulations are unconditionally stable. Moreover, this permits an unlimited variation of the time-step size. Complex two-dimensional geometries and local mesh refinements can be resolved by semi-structured grids, i.e. using hanging nodes, on the basis of general curvilinear coordinates. Moving meshes are implemented based upon the space-conservation law, such that (forced) unsteadiness caused by moving boundaries can be captured.

The cartesian tensor components as well as all scalar variables are stored in the control-volume centers. Diffusive terms are approximated applying central difference schemes whereas convective terms are approximated by upwind biased, bounded, monotonic functions of higher order, Lien, Leschziner, [4]. The governing equations of the incompressible mean flow field, the continuity equation and the transport equation(s) for turbulent quantities are solved sequentially. The solution is iterated to convergence using a compressible (all-speed) pressure correction scheme of the SIMPLE type, Patankar, [5]. Introducing apparent pressures and viscosities, a generalised Rhie & Chow interpolation avoids an odd-even decoupling of pressure, velocity, and Reynolds stress $\overline{u_i u_j}$, Obi, Peric, Scheurer, [6].

The computational modelling of turbulent effects is based on a variety of models featuring different degrees of complexity to compute the Reynolds stresses arising in the averaging process. Amongst the many implemented modelling practices are standard isotropic Boussinesq-viscosity approaches employing either one or two additional transport equations ( $k-\epsilon$, $k-\omega$). Moreover, their extension to non-linear constitutive stress-strain relations, as well as (non-linear) full Reynolds-stress transport models are available, Rung, Lübcke, Xue, Thiele, Fu, [7]. In the present study, however, attention is restricted to the assessment of one- and two- equation models: the standard (SA) and a modified (SALSA) Spalart-Allmaras one-equation model, and only two different $k-\omega$ two-equation models in comparison, Rung, Thiele, [8], because turbulence models of higher order did not show better results.

Due to the implementation of grid fluxes arising from the formulation of transport equations on moving grids, the step from grids moving in a rigid-body manner to deforming grids does not require any changes in the part of the code solving the equations for the fluid.

2.2 Structure Mechanics

2.2.1 Sail
The model presented in this paper describes the behaviour of an inextensible, massless sail using the structural model of a chain given in de Matteis, de Socio, [3]. This description was preferred in the beginning since an analytical function for the membrane shape as in the sail equation does not easily allow for a computation of the exact movement of the sail's material points. Equilibrium of forces is evaluated at each grid node of the sail, see fig. 1, leading to the following equation:

$$\underline{F}(i-1)=\underline{a} (i) l(i) + \underline{F} (i)\;, \;\;\; \; l(i)=\sqrt{\left(\underline{x}(i)-\underline{x}(i-1)\right)^2}\;,$$ (1)

where $\underline{F}(i)$ stands for the force vector acting on the hinge connecting neighboring sail elements and $l(i)$ for the element length. Assuming that no friction occurs in the hinge and that the aerodynamic load $\underline{a}(i)$ acts in the center between $\underline{x}(i-1)$ and $\underline{x}(i)$, the moment equilibrium is given as:

$$0 = \big[ \frac{1}{2} l(i) \underline{a}(i) + \underline{F} (i) \big] \times \big[\underline{x}(i)-\underline{x}(i-1)\big]\;.$$ (2)

Note that there are no time dependent terms included, thus inertia and damping are completely neglected in the structure mechanics. This simplification has become necessary due to a lack of experimental data concerning mass distribution and structural damping. However, the iterative procedure to converge the solution of these equations has to be started with a relaxation factor of $0.001$ that can slowly be raised each step up to a maximum of 0.5.

Solving the sail equation

$$\frac{y''}{\sqrt{(1+(y')^2)^{3}}}= - \frac {\vert\underline{a}_n\vert}{T} \,\, ,$$ (3)

as also given in [3], on the other hand, a direct (steady) solution for the sail shape as $y=y(x)$ is obtained. Here $T$ stands for the tension in the sail and $\vert\underline{a}_n\vert$ for the amount of aerodynamic force normal to the sail(-element). In order to avoid an unphysical - and thus unfavourable for the grid-deformation algorithm - deformation of the sail elements, the sail points have to be redistributed on this function according to local strain and / or constant length of each element. This and the fact that a direct solution might lead to unstable coupled solutions, especially with bad initial sail shapes far from the final result, motivates the use of the described chain model.

2.2.2 Computational Grid
The starting point for the grid-deformation method is the momentum conservation equation valid for a solid body in steady state without any volume loads:

$$0 = \frac{\partial{T_{ij}}}{\partial{x_{j}}} \, ,$$ (4)

where $\underline{\underline{T}}$ represents the stress tensor. Displacements are imposed on the continuum representing the computational grid via given displacement values on the boundaries.

Assuming ideally linearly elastic material, the stress-strain relation becomes

$$T_{ij}=E_{ijkl} \frac{1}{2} \left( \frac{\partial{v_{k}}}{\partial{x_{l}}} +\frac{\partial{v_{l}}}{\partial{x_{k}}} \right)\;$$ (5)

with only linear terms in the strain-displacement relation where $v_i$ are displacements in the direction of the coordinates and $E_{ijkl}$ is the generalized Hooke-material tensor enabling an arbitrary link of displacement derivatives and stresses (anisotropy). However, the grid is idealized as isotropic material with zero Poisson's ratio reducing the material tensor to an isotropic parameter $E$ aligning stress and strain tensor. This is done in order to avoid transverse contraction of grid cells that could result in a loss of grid quality.

It has to be remarked, that although the material is considered as isotropic, it is not homogeneous, i.e., the scalar $E$, representing Young's modulus or the material stiffness, can vary in space. Thus, the grid stiffness depends on grid-cell volume or deformation undergone (strain energy).

It is desirable to keep the grid orthogonal in near-wall regions, where it is usually very fine, but in order to balance deformation it is necessary to avoid strain weakening, so one possible formulation for $E$ could be:

$$E=\left(\frac {min VOL}{VOL}\right)^\alpha \left(\frac {D_{ii}^2}... ...al{v_{i}}}{\partial{x_{j}}} +\frac{\partial{v_{j}}}{\partial{x_{i}}} \right)\;,$$ (6)

where $min$ stands for the minimum value in the whole computational domain and the exponents $\alpha$ and $\beta$ can be adjusted to requirements: if it is desirable to keep a section of the grid with very fine grid cells nearly undeformed, $\alpha$ has to be raised, whereas a higher $\beta$ tends to distribute the deformation of the grid equally. The best formulation for the cell stiffness will depend on the application, thus there is no general set of parameters $\alpha$ and $\beta$ that constitutes an optimum.

In order to implement this algorithm in the previously described code, the same FV method is applied for the calculation of grid deformation where the grid is treated as a continuum. In the absence of convective and transient terms, diffusive terms are approximated applying central difference schemes and a harmonic interpolation of cell-center values on cell faces is used.

2.3 Coupled Procedure

Loose coupling is achieved by a sequential solution of the governing equations for the fluid domain for every time step based upon the sail shape from the previous step. Applying the loads of the converged flow field, a new sail shape is computed and the iteration is proceeded by one time step.

Though there are more elaborate methods available for coupling, this simple one was chosen since there is no time-dependent term in the sail equation implemented as yet and time-step refinement is more straightforward than any coupling over partitioned time steps. Moreover, the method imposes no constraint on the time-step size consistent with the implicit method for the flow.

The iterative procedure is started with an (arbitrary) initial shape obtained from a polynomial and a grid generated manually around it (or with a grid and sail shape from a previous time step), and the unsteady flow computation is proceeded by one time step. With the resulting force $\underline{a}(i)$ due to the flow solution on each element, eqns (1) and (2) are considered. Eqn (1) gives an explicit expression for the nodal forces that can be calculated depending on the force on the first node which can be arbitrarily initialized. In conjunction with eqn (2), the length of each sail element has to be kept constant. Along with the boundary condition of fixed sail ends, a corrected force can be determined as a constant to start with eqn (1) again. This iterative process is continued until a certain level of convergence is reached.

Afterwards, displacements resulting from the change in sail shape are imposed on the grid and a new grid is computed, and the computation can be proceeded by one time step again.

3 Results

The fully coupled computations are started at $3^o$ angle of attack where the sail shape is awaited to be stable and separation bubbles should be small - if existent at all. After convergence, the angle is raised step by step up to a chosen maximum of $20^o$ far beyond separation, and on another ``computational branch'' decreased step by step to an angle, where the algorithm fails to give a solution, fig. 2. In this case the incidence angle has become negative below a critical value, where tension in the sail is close to zero and luffing occurs. Lowering the angle further, a hysteresis would occur. Thus, the sail shape becomes unstable and the deformations exceed the validity of the method, which cannot yet handle luffing that occurs earlier at higher excess lengths. However, the algorithm gives stable solutions for negative incidence close to $0^o$. The onset of numerical divergence roughly matches the angle where luffing occurs in the experiments.

Figure 2:

Evolution of lift coefficient vs. time for varying angle of attack during computation and schematic drawings of streamlines around sail at different angles.

Computational results are presented as diagrams of lift coefficient vs. angle of attack for three different excess lengths, in fig. 3 for $\varepsilon=0.21\%$ (left) and $\varepsilon=1.14\%$ (right), and in fig. 4 for $\varepsilon=4.22\%$ (left), for a Reynolds number based on chord length of $Re=\frac{U_{\infty}\cdot c }{\nu}=1.3 \cdot 10^6$ obtained on grids with 16.000 control volumes covering a physical domain of $\frac{X}{c}$ from $-10$ to $11$ and $\frac{Y}{c}$ from $-9$ to $9$ at $\frac{X}{c}=-10$ and $-12$ to $12$ at $\frac{X}{c}=11$, fig. 1. This number was found sufficient when evaluating the effects of grid refinement on the prediction of lift for these excess lengths. As low-Reynolds-number boundary conditions were used, the wall-normal distance $y^+$ of the first grid points was kept around $1$. The computations are compared to experimental results and a linear approximation

$$C_l \cong 2 \pi \alpha \pm 0.7 \; (100 \; \varepsilon)^{\frac{1}{2}}\;,$$ (7)

given in [2]. The $\pm$ sign stands for a change in curvature when the sail snaps through, where the negative branch represents the sail bent downwards.

Figure 3:

Lift coefficient versus angle of attack for excess length $\varepsilon=0.21\%$ (left), and $\varepsilon=1.14\%$ (right), different turbulence models.

At low angles of attack good agreement can be found with the linear approximation as well as with the experimental results for small excess lengths. For $\varepsilon=0.21\%$, the shortest sail, no distinctive maximum in computed lift can be seen, whereas for longer sails the loss in lift after separation increases and a maximum lift coefficient originates roughly at $7^o-8^o$, its height dependent on the sail length. For $\varepsilon=4.22\%$ there are only very few experimental data points, so comparability is very poor, moreover, only the one-equation models obtained stable solutions. In this case the loss in lift due to separation is significantly underpredicted and the linear approximation does not hold anymore at all, either.

This reveals that approaching onset of separation, the computational results depend strongly on subtle details pertaining to the representation of turbulence. The one-equation models capture most of the experimental results but fail to render the effects of massive separation, whereas the other models tend to produce more separation resulting in a lower lift coefficient at lower angle of attack. Grid-refinement posed even more problems with stability in these cases, revealing that the effects of inertia and damping on the sail shape computation will have to be studied.

The results for lift with separation are averaged values as the flow becomes unsteady before full separation has developed. Fig. 1 shows the extremely unsteady behaviour of the system at these incidence angles for $\varepsilon=1.14\%$. Unfortunately, this behaviour is not represented in any experimental result, so nothing can be said about the physical reliability of the frequency computed, which is roughly $0.3$Hz, which corresponds to a Strouhal number $St=\frac{f \cdot c }{U_{\infty}}$ of lower than 0.03 in this case. At least, turbulent and mean flow frequencies can thus be estimated to be around three orders of magnitude apart, Rung [9], so that applying RANS methods to these unsteady flows does not contradict the basic assumptions of statistical turbulence modelling.

Figure 4:

Lift coefficient versus angle of attack for excess length $\varepsilon=4.22\%$ (left), different turbulence models. Comparison of sail shapes for different excess lengths and angles of attack (right).

With higher excess length $\varepsilon$ not only the maximum lift but also the sail deformation grows, fig. 4 (right), and the grid-deformation algorithm has to meet these requirements. This can easily be done with the proposed procedure as long as the sail deformation is small. When luffing occurs and the sail snaps through, the loss in grid quality is to high for the procedure to work.

Figure 5:

Streamlines for LLR-$k-\omega$ model for excess length $\varepsilon=1.14\%$, angle of attack: $\alpha=20^o$ (left), $\alpha=9^o$ (right).

In figs. 5 and 6 streamline pictures are shown for angles of attack $\alpha$ ranging from $-0.5^o$ to $20^o$. For high angles of attack beyond maximum lift, fig. 5 (left), the flow is fully separated on the upper side of the sail and steady vortices occur in the wake with a similar streamline pattern and pressure distribution for angles of attack above $14^o$.

Figure 6:

Streamlines for LLR-$k-\omega$ model for excess length $\varepsilon=1.14\%$, angle of attack: $\alpha=2^o$ (left), $\alpha=-0.5^o$ (right, smaller scale).

For lower angles of attack, fig. 5 (right), the streamlines reveal a separation bubble limited in size on the upper side that is also visible in the $C_P$ distribution. This separation bubble grows with increasing angle of attack until it becomes unstable in size. Due to the sharp leading edge - the sail has zero thickness - this separation only vanishes when streamlines and sail are nearly parallel at the front, fig. 6 (left), roughly at $2^o$.

On the other hand, when the angle of attack decreases, a separation bubble occurs on the lower side of the sail, fig. 6 (right), only visible in a magnification or the pressure distribution. In this case the sail is more and more likely to become unstable or snap through.

4 Conclusions and Outlook

The chain model works in conjunction with the RANS-method and all turbulence models investigated are able to capture the experimental results at low angle of attack and the weakly non-linear area is predicted with success for small excess lengths $\varepsilon$. In case of massive separation and for long sails, however, severe deviations and differences between the turbulence models show up and between the experiments and computations.

Hence, more research is considered necessary and additional effects will have to be taken into account, such as transitional effects and three-dimensionality. Moreover, work in progress is directed towards higher-order closure models and obtaining and evaluating more experimental data. Another step will be the introduction of damping and inertia into the sail algorithm to stabilize the procedure which will be necessary for the computation of a full hysteresis including luffing and to be able to handle even higher excess lengths $\varepsilon$. An examination of the effects of grid refinement on the frequency of oscillations and sail shape is also considered necessary.

To further validate the procedure against analytical results for small excess lengths $\varepsilon \ll 1\%$, inviscid computations are aimed at, together with laminar computations as a starting point for introducing fixed and variable transition in a next step. So far, inviscid computations could not be performed due to the loss in stability of the numerical procedure when viscosity is small.

Acknowledgment

This research was partly funded by the Deutsche Forschungsgemeinschaft (DFG-SFB 557 'Beeinflussung komplexer turbulenter Scherströmungen').

References

[1]

Smith, R, Shyy, W. Computation of aerodynamic coefficients for a flexible membrane airfoil in turbulent flow: A comparison with classical theory. The Physics of Fluids, 8(12):3346 - 3353, 1996.

[2]

Greenhalgh, S., Curtiss, H. C., Smith, B. Aerodynamic Properties of a Two-Dimensional Inextensible Flexible Airfoil. AIAA Journal, 22(7):865-870, 1984.

[3]

de Matteis, G., de Socio, L. Nonlinear Aerodynamics of a Two-Dimensional Membrane Airfoil with Seperation. J. Aircraft, 23(11):831-836, 1986.

[4]

Lien, F.S., Leschziner, M.A. Upstream Monotonic Interpolation for Scalar Transport with Application to Complex Turbulent Flows. Int. Journal for Numerical Methods in Fluids, 19:527-548, 1994.

[5]

Patankar, S.V. Numerical Heat Transfer and Fluid Flow. McGraw Hill, New York, 1980.

[6]

Obi, S., Peric, M., Scheurer, G.. Second Moment Calculation Procedure for Turbulent Flows with Collocated Variable Arrangement. AIAA Journal, 29(4):585-590, 1991.

[7]

Rung, T., Lübcke, H., Xue, L., Thiele, F., Fu, S. Assessment of explicit algebraic Reynolds-stress models in transonic flows, volume 4 of Engineering Turbulence Modelling and Experiments. Elsevier, Amsterdam, 659-668, 1999.

[8]

Rung, T., Thiele, F. Computational modelling of complex boundary-layer flows. In Proc. 9th Intern. Symposium on Transport Phenomena in Thermal-Fluid Engineering, pages 321-326, Singapore, 1996.

[9]

Rung, T. Entwicklung anisotroper Wirbelzähigkeitsbeziehungen mit Hilfe von Projektionstechniken. Doctoral Dissertation accepted by: Technical University of Berlin, Berlin, 2000.