Development of a high fidelity fluid-structure interaction

solver : towards flexible foils simulation

Thomas Fabbri

To cite this version:

Thomas Fabbri. Development of a high fidelity fluid-structure interaction solver : towards flexible
foils simulation. Fluid Dynamics [physics.flu-dyn]. Université Grenoble Alpes [2020-..], 2022. English.
�NNT : 2022GRALI002�. �tel-03623578�

HAL Id: tel-03623578

https://theses.hal.science/tel-03623578
Submitted on 4 Apr 2022

HAL is a multi-disciplinary open access L’archive ouverte pluridisciplinaire HAL, est

archive for the deposit and dissemination of sci- destinée au dépôt et à la diffusion de documents
entific research documents, whether they are pub- scientifiques de niveau recherche, publiés ou non,
lished or not. The documents may come from émanant des établissements d’enseignement et de
teaching and research institutions in France or recherche français ou étrangers, des laboratoires
abroad, or from public or private research centers. publics ou privés.
THÈSE
Pour obtenir le grade de

DOCTEUR DE L'UNIVERSITE GRENOBLE

ALPES
Spécialité : Mécanique des fluides, Procédés, Energétique
Arrêté ministériel : 25 mai 2016

Présentée par

Thomas FABBRI
Thèse dirigée par Guillaume BALARAC, Professeur, Université Grenoble Alpes (UGA), co-
encadrée par Vincent MOUREAU, Chargé de recherche CNRS, CORIA et
co-encadrée par Pierre BENARD, Maître de Conférences, INSA Rouen

préparée au sein du Laboratoire des Ecoulements Géophysiques et Industriels (LEGI)

dans l'Ecole Doctorale I-MEP2 - Ingénierie - Matériaux, Mécanique,

Environnement, Energétique, Procédés, Production

Développement d’un solveur haute

fidélité d’interaction fluide-structure :
vers la simulation de profils flexibles
Development of a high fidelity fluid-
structure interaction solver : towards
flexible foils simulation
Thèse soutenue publiquement le 03/03/2022,
devant le jury composé de :

Franck NICOUD
Professeur des Universités, Université de Montpellier, Rapporteur
Emmanuel LEFRANÇOIS
Professeur des Universités, Université de technologie de Compiègne,
Rapporteur
Rémy DENDIEVEL
Professeur des Universités, Grenoble INP, Président du Jury
Guillaume DE NAYER
Tenured Senior Researcher, Helmut Schmidt Universität, Examinateur
Mathieu OLIVIER
Professeur agrégé, Université Laval, Examinateur
Guillaume BALARAC
Professeur des Universités, Grenoble INP, Directeur de thèse Vincent
MOUREAU
Chargé de Recherche CNRS, Université de Rouen, Membre invité
Pierre BENARD
Maître de Conférences, Université de Rouen, Membre invité
 Remerciements

Je tiens tout d’abord à remercier Emmanuel Lefrançois et Franck Nicoud qui ont accepté
de rapporter cette thèse avec la plus grande attention. Je souhaite également remercier les
autres membres du jury pour le temps et la considération qu’ils ont accordé à mon travail,
ce qui a permit des échanges que j’ai trouvé très enrichissant lors de ma soutenance.
Un immense merci à mes deux co-encadrants, Vincent et Pierre. Vous m’avez donné la
chance de me reconvertir dans un domaine qui était bien éloigné de mes études en m’acceptant
en stage et en m’accordant votre confiance. C’est d’ailleurs grâce au temps et la formation
que vous m’avez fourni lors de mon passage au CORIA que j’ai été en mesure d’effectuer ce
travail.
Bien entendu, un merci tout particulier à mon directeur de thèse, Guillaume. Travailler à
tes côtés durant ces années a été un immense honneur et un grand plaisir. Tu m’as laissé une
grande liberté dans la direction de ce travail, et pour autant tu as toujours été derrière moi
pour me remonter le moral et me remotiver. Malgré ton emploi du temps de ministre, tu as
toujours pris le temps de répondre à mes questions et je m’estime très chanceux d’avoir eu
un directeur de thèse aussi formidable que toi, autant au niveau scientifique que sur le plan
humain.
Cette thèse a été marquée par plusieurs périodes de confinement, ce qui a souligné
l’importance de l’environnement de travail. J’ai pour ma part intégré la team MoST, dans
laquelle règne une atmosphère particulièrement conviviale. Je pense tout d’abord à Patrick
qui, en plus de mettre à notre disposition des moyens de calcul très conséquents et de nous
dépanner régulièrement, illumine les pauses cafés grâce à sa bonne humeur et à ses anecdotes
sur ses (petits) trails de 80km. Je remercie également tous les autres membres de l’équipe
qui participent à créer cette bonne ambiance générale; Giovanni, Manu (même si ses sar-
dines emboucanent la salle de pause), Matthieu (merci encore pour OP), Cyril (l’homme
trompette), Jeremy (qui avait tendance à confondre mon bureau et son salon), Guillaume
(ou Mr “café?” à 19h), François (même s’il s’agit d’un individu détestable), Eliott, Hugo, et
tous les autres y compris les supers stagiaires comme Célia. Enfin un grand merci à mon
irremplaçable co-bureau, j’ai nommé Savinien. Ta rigueur et ton sérieux m’ont énormément
tiré vers la haut durant ces 3 ans (et quelques..), et nos conversations interminables sur des
problèmes physiques et ta Normandie natale vont me manquer. Vous avez tous participé à
rendre ces années bien moins pénibles, et c’était un plaisir de passer tout ce temps avec vous
derrière un café ou une (une?) bière.
Cette thèse s’est déroulée dans le cadre du projet ANR DYNEOL, dont je tiens à re-
mercier tous les acteurs. C’est également l’occasion d’exprimer ma reconnaissance envers le
contribuable français qui finance des beaux projets de recherche comme celui-ci, et je réalise la
chance que j’ai eu de naı̂tre et de grandir dans un système qui permet cela. C’est grâce à

iii
lui et aux bourses mises en place que j’ai pu effectuer mes études post-bac dans des bonnes
conditions, et je pense qu’il est important de se rappeler que cela n’aurait malheureusement
pas été possible partout.
Pour terminer, je tiens à remercier mes proches. Tout d’abord pour le cadre familial qui
m’a fait comprendre que je devais prendre mes études au sérieux. Une pensée pour mon
adorable soeur qui a accepté sans hésiter de relire mon manuscrit pour y corriger les fautes
d’anglais, ce qui n’était pas une mince affaire. Merci aussi à mes amis sur qui je pouvais
toujours compter pour me changer les idées. Enfin, je souhaite exprimer toute ma gratitude
pour celle qui m’accompagne depuis toute ces années et sans qui “rien n’aurait été possible”
(comme ça c’est bon?), Yasmina.

iv
 à vous deux, Gino et Yohann

“La science n’est pas indépendante des femmes et des hommes qui la pratiquent.
La science n’est pas une froide découverte du réel “en lui même“.
Elle est, avant tout, une ”manière de faire un monde“.
Une manière élégante et cohérente mais évidemment non unique et non hégémonique.”
Aurélien Barrau

v
 Contents

Preamble: on the necessity of energy transition 5

1 Introduction 15
1.1 Motivations: experimental research about chordwise flexible blades turbines . . 15
1.2 Numerical methods for FSI simulation . . . . . . . . . . . . . . . . . . . . . . 18
1.2.1 Generalities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
1.2.2 Monolithic and partionned approaches . . . . . . . . . . . . . . . . . . 20
1.2.3 Numerical simulations of chordwise flexible blades turbines . . . . . . . 21
1.3 Objectives and plan . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

2 Arbitrary Lagrangian-Eulerian solver 25

2.1 Presentation of YALES2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
2.2 ALE solver . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
2.2.1 Numerical method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
2.2.1.1 Prediction step in ALE . . . . . . . . . . . . . . . . . . . . . 28
2.2.1.2 Correction step . . . . . . . . . . . . . . . . . . . . . . . . . . 30
2.2.1.3 Large-Eddy Simulation . . . . . . . . . . . . . . . . . . . . . . 30
2.2.1.4 Dynamic Mesh Adaptation . . . . . . . . . . . . . . . . . . . 32
2.2.2 Application case: simulation of a Darrieus turbine . . . . . . . . . . . . 34
2.2.2.1 Case presentation . . . . . . . . . . . . . . . . . . . . . . . . . 34
2.2.2.2 Numerical setup . . . . . . . . . . . . . . . . . . . . . . . . . 36
2.2.2.3 Results and discussion . . . . . . . . . . . . . . . . . . . . . . 39

3 Structural Mechanics Solver 43

3.1 Numerical method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
3.1.1 Generalities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
3.1.1.1 History of FEM . . . . . . . . . . . . . . . . . . . . . . . . . . 44
3.1.1.2 Shape functions . . . . . . . . . . . . . . . . . . . . . . . . . . 44
3.1.2 Stationary problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
3.1.2.1 Stationary linear elasticity problem of one finite element . . . 47
3.1.2.2 Displacement boundary conditions . . . . . . . . . . . . . . . 50
3.1.2.3 Matrices assembly . . . . . . . . . . . . . . . . . . . . . . . . 51
3.1.2.4 Stationary non-linear problem . . . . . . . . . . . . . . . . . . 54
3.1.3 Dynamic problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
3.1.3.1 Generalities . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
3.1.3.2 Newmark scheme . . . . . . . . . . . . . . . . . . . . . . . . . 58

1
Contents

3.1.3.3 Generalised-α scheme . . . . . . . . . . . . . . . . . . . . . . 60

3.1.4 Use of reference space - Quadratic finite elements . . . . . . . . . . . . 62
3.1.4.1 Principle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
3.1.4.2 Implementation in the SMS . . . . . . . . . . . . . . . . . . . 66
3.1.4.3 Forces computation . . . . . . . . . . . . . . . . . . . . . . . . 67
3.1.5 Possible improvements . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
3.2 Verification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
3.2.1 Comparison with Ansys . . . . . . . . . . . . . . . . . . . . . . . . . . 71
3.2.1.1 Stationary problem . . . . . . . . . . . . . . . . . . . . . . . . 71
3.2.1.2 Dynamic problem . . . . . . . . . . . . . . . . . . . . . . . . . 73
3.2.2 Mesh convergence study with the implemented elements . . . . . . . . 81
3.2.2.1 2D finite elements . . . . . . . . . . . . . . . . . . . . . . . . 81
3.2.2.2 3D finite elements . . . . . . . . . . . . . . . . . . . . . . . . 84

4 FSI solver 87
4.1 Mesh Movement Solver . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
4.1.1 Pseudo solid method . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
4.1.2 Implementation in YALES2 . . . . . . . . . . . . . . . . . . . . . . . . 88
4.1.3 Adjustment of pseudo-material properties . . . . . . . . . . . . . . . . 90
4.2 FSI coupling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
4.2.1 Data transfer via CWIPI . . . . . . . . . . . . . . . . . . . . . . . . . . 94
4.2.2 Weak coupling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
4.2.2.1 Case presentation . . . . . . . . . . . . . . . . . . . . . . . . . 97
4.2.2.2 Parallel synchronous scheme . . . . . . . . . . . . . . . . . . . 97
4.2.2.3 Serial staggered scheme . . . . . . . . . . . . . . . . . . . . . 100
4.2.3 Partitioned semi-implicit predictor-corrector coupling scheme . . . . . . 101
4.2.4 FSI coupling with DMA . . . . . . . . . . . . . . . . . . . . . . . . . . 105

5 Validation against 2D numerical laminar cases 107

5.1 Case presentation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107
5.2 CFD tests . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108
5.3 CSM tests . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110
5.4 FSI tests . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112
5.4.1 FSI3 case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112
5.4.2 FSI2 case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116

6 Validation against 3D experimental turbulent cases 119

6.1 Kalmbach experimental case: flexible plate behind a cylinder . . . . . . . . . . 119
6.1.1 Case presentation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119
6.1.2 Numerical setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121
6.1.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123
6.2 Hoerner experimental case: flexible pitching foil . . . . . . . . . . . . . . . . . 128
6.2.1 Case presentation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128

2
 Contents

6.2.2 Rigid foil simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130

6.2.3 Flexible foil simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . 134

7 Conclusion & Perspectives 141

7.1 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141
7.2 Perspectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143

A Shape functions and Gauss points 147

Bibliography 157

3
 Preamble: on the necessity of
energy transition

Although this work focuses primarly on the numerical simulation of fluid-structure inter-
actions, it is interesting to begin with a larger contextualisation. It is thus introduced by a
part aiming to explain the interest of studying and improving wind or marine turbines tech-
nology. It will provide a quick overwiew of causes and consequences of global warming, and
it will show why energy transition, which includes the development of such sources of energy,
has become a necessity. The interested reader is invited to refer to the Intergovernmental
Panel on Climate Change (IPCC) report [1] for further details.
The climate crisis comes from an analysis: global temperature is increasing. This is
illustrated by temperatures curves given in Fig. 0.1. Despites of this undeniable situation,

Fig. 0.1: Left: Change in global surface temperature (decadal average) as reconstructed
(1-2000) and observed (1850-2020). Right: Change in global surface temperature (annual
average) as observed and simulated using human & natural and only natural factors (both
1850-2020). Extracted from [2].

causes may seem unclear. Climate evolution is indeed a complex science and it depends on
a lot of factors. For instance, solar activity variations has an impact on climate but some
studies have shown that in reality this latest is minor [3]. Climate can also change due
to deviation of earth movement and rotation around the sun, called Milankovitch cycles.
Even though this has a strong influence on climate, the time scale (cycles of about 100 000

5
Contents

years) does not fit with the previous observations [4, 5, 6]. The hypothesis supported by a
large majority of climate scientists is that this global warming is due to human activities,
especially greenhouse effect gas emissions. That is why Fig. 0.1 also shows simulations of
climate evolution without human impact and it remains nearly constant.
Greenhouse effect is explained in Fig. 0.2. It involes a difference in transparency of the
atmosphere depending on the radiation that tries to pass through. When incoming solar

Fig. 0.2: Schematic diagram of the global mean energy balance of the Earth. Numbers indi-
cate best estimates for the magnitudes of the globally averaged energy balance components
together with the uncertainty ranges. Extracted from [7].

radiation arrives at the top of the atmosphere, 30% is directly reflected towards space by
clouds (20%), the various layers of the atmosphere (6%), and the surface of the earth (4%),
which includes some places – for example the ice caps – that are particularly highly reflexive
(albedo effect). The remainder of this incoming energy is absorbed by the various components
of our planet (continents, oceans, atmosphere), heating them, and then re-emitted in the form
of infrared energy, as black body radiation. The greenhouse gases here refer to gases that let
pass solar radiation but are opaque to these infrareds. These gases then also heat up because
of infrareds absorption, re-emitting infrareds that goes back to the ground, heating it up a

6
 Contents

second time. The whole process finally reaches an equilibrium, but with higher temperatures
than if these greenhouse gases where not present [7, 8].
Without this greenhouse effect, the average surface temperature would be -18°C; a tem-
perature so low that all water on Earth would freeze, the oceans would turn into ice and life,
as it is known, would not exist [9]. However, gas with this effect are minor ingredients of this
atmosphere. 60% to 70% of the Earth’s greenhouse warming is induced by water vapor, while
carbon dioxide (CO2 ) provides just a few degrees. Moreover, it is important to notice that
Earth atmosphere had known different greenhouse gases concentrations, including CO2 , and
sometimes well beyond today. The natural mechanisms explaining atmosphere CO2 concen-
tration are complex and it exists a lot of positive and negative feedback loops; any interested
reader can refer to [10] for more details. Consequently, as regards the above points, accurate
climate prediction turns out to be really challenging.
Nevertheless, climate scientists keep thinking that the temperature increase is caused by
greenhouse gas emissions due to human activities. These emissions are indeed increasing
since the end of the XIX-th century, as shown by Fig. 0.3. They are in major part caused
by energy production, as it is mainly based on fossil fuels combustion, releasing CO2 . These
energy sources have allowed important economic growth, with counterpart an increase in
global temperatures. That is why temperature increase in 2081-2100 is estimated at around
+3°C in comparison with 1986-2005 according the scenario “business as usual” [11].
Furthermore, for a deeper understanding of the climate crisis, it is essential to explain
why climate change represents a problem. This is detailed in the IPCC report [13] from
where Fig. 0.4 has been extracted. It shows world regions where climate change affects
physical systems, biological systems and human activities. Examples arising include water
stress, disrupting food production, and also sea level rise, flooding densely populated regions.
Projections also show that the number of days per year where temperature and humidity
conditions are lethal for human will increase in several countries [14].
It is then obvious that global warming induces significant mass migration [15], disrupts
human societies, and is causing irreversible damages on biodiversity [16, 17], and that these
effects will increase if greenhouse gas emissions continue to grow [13]. One pathway would be
to take carbon from the atmosphere and to store it apart. This Carbon Capture and Storage
(CCS) strategy could indeed decrease atmosphere CO2 concentration, but this is limited and
demanding to develop at large scale [18, 19, 20, 21]. As a result, a reduction of human
greenhouse gas emissions is mandatory. This reduction is in fact extremely challenging given
that a lot of different sectors are involved, as shown by Fig. 0.5. Therefore, changes have to
be made in most human activities.
As explained above, human greenhouse gas emissions are mainly caused by energy pro-
duced with fossil fuels. Nevertheless, it exists other sources of energy with lower carbon
intensity exist. Figure 0.6 gives greenhouse gas impacts of the different electricity supply
technologies, integrating emissions during their entire lifecycle. It can be seen that nuclear
and renewable energies such as solar, wind or hydropower induce low carbon emissions in
comparison with gas and coal. However, nuclear power requires high initial investment, good
security controls, and is often not accepted by local population. In addition, renewable ener-

7
Contents

Fig. 0.3: Greenhouse gas emissions by sector (top) and primary energy production by source
(bottom). From [12].

8
 Contents

Fig. 0.4: Widespread impacts in a changing world. Impacts are shown at a range of geographic
scales. Symbols indicate categories of attributed impacts, the relative contribution of climate
change (major or minor) to the observed impact, and confidence in attribution. Extracted
from [13].

9
Contents

Fig. 0.5: Global greenhouse gas emissions by sector for the year 2016. Extracted from [22].

10
 Contents

Fig. 0.6: Emissions of selected electricity supply technologies (gCO2 eq.kWh−1 ). Extracted
from [23].

gies as wind or solar are intermittent and so hardly compatible with the necessity to ensure
a stable energy production given that energy storage remains very difficult in large scale. As
for hydropower, it depends on regions topography and most viable sites are already exploited
in Europe. Those reasons explain why fossil fuels have been mainly used during the XX-th
century, and make the energy transition challenging. Research works in this domain are then
encouraged but it may seem unreasonable to wait for new technology given the task scope.
As a matter of fact, it must also be noticed that the problem is not only about electricity
production, as shown by Fig. 0.7. With the example of France, it highlights that half of its
energy consumption comes from fossil energies, despite a low carbon electricity production
system (58.5 gCO2 eq.kWh−1 in 2016, while it was 440.8 gCO2 eq.kWh−1 in Germany [24]).
This is explained with Fig. 0.5; energy is indeed used in different domains under different
forms. For instance, road transports use mainly oil and heating in buildings is often done with
gas combustion. Besides, it has to be reminded that France, just as most occidental countries,
has relocated a part of its industries in other countries (in Asia for instance), so that most
of the consumer goods sold are made elsewhere, and gas emissions due to manufacturing of
these products do not appear in France carbon budget [25]. This pollution offshoring induced
by globalization imposes a common effort of every country, and is meaningless concerning
the global climate given that carbon emissions impact the atmosphere independently of its
location.
Improving the energy efficiency, like developing thermal insulation, and reducing the
energy consumption by changing individual behaviour are necessary actions, but are still
insufficient [26, 27]. A strong global effort is then necessary to produce electricity with
low carbon systems (nuclear and renewable energies), and also to electrify our systems in
agriculture, transports, buildings, etc...

11
Contents

Fig. 0.7: Primary energy consumption in France between 1980 and 2016. From [12].

Even supposing that human activities would not have any consequences on climate, or
climate change would not cause any damage to human society, the situation would still be
worrying because of fossil fuels rarefaction [28], as illustrated by Fig. 0.8. Fossil ressources

Fig. 0.8: Evolution of the world recoverable oil and gas discoveries, except US and Canada.
Extracted from [8].

stocks are indeed in limited quantities on Earth, and are regenerated at timescale incompat-

12
 Contents

ible with the current consumption. Changes in energy production will then be voluntary or
suffered.
In order to limit the climate change and reduce greenhouse gas emissions, several energetic
and social scenarios are studied by IPCC [1], or by national agencies [29]. They predict a
decrease of energy consumption but an increase in electricity production to substitute current
uses of fossil fuels. Consequently, in most of these scenarios, renewable energies have to be
developed to ensure a low carbon electricity production system. It has indeed been shown
that they can reduce significantly carbon emissions, even in France where the nuclear power
is the main source of electricity [30].

It has been explained why developing renewable energies is crucial. The present research
work is part of this energy transition strategy, and is motivated by the need of improving
wind or marine turbines technology. The first Chapter will go over one of these possible
improvements, the use of flexible foils.

13
 Chapter 1

Introduction

This chapter aims at introducing the objectives of this work. It starts by explaining what
motivates the study, providing a literature review of experimental research results involving
chordwise flexible blades. Then, the FSI phenomenon is presented, as the principal numerical
methods used to simulate cases of that kind. State of the art of chordwise flexible foil numerical
simulations is given. As a result, detailed objectives and a work plan can be described.

Contents
1.1 Motivations: experimental research about chordwise flexible blades turbines . . . . . 15
1.2 Numerical methods for FSI simulation . . . . . . . . . . . . . . . . . . . . . . . . . . 18
1.2.1 Generalities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
1.2.2 Monolithic and partionned approaches . . . . . . . . . . . . . . . . . . . . . 20
1.2.3 Numerical simulations of chordwise flexible blades turbines . . . . . . . . . . 21
1.3 Objectives and plan . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

1.1 Motivations: experimental research about

chordwise flexible blades turbines
The necessity to develop and improve wind or marine tubines keeps rising, and today
a majority of these turbines are Horizontal Axis Turbines (HAT). They present indeed a
better efficiency than Vertical Axis Turbines (VAT). However, VAT have several advantages,
as detailed in [31]. At first, they are not sensitive to stream direction changes, which has a
major impact with large scale turbine where flow conditions can be different at the bottom
and at the top of the turbine. Besides, the gravity center is at the base facilitating the
maintenance, and blade shape is the same along the span, implying simpler manufacturing.
Finally, several studies [32, 33] showed that in the case of tubine farm, the power extracted for
a given surface tends to be superior with VAT than HAT, due to turbines wakes for instance.
As illustrated in Fig. 1.1, blades of VAT experience variable flow conditions during a rota-
tion. One challenge is then to find an optimal shape that remains highly efficient regardless
any conditions in order to maximize energy harvesting. Therefore, a static shape can seem
not optimal. A bioinspired approach to control and improve stall characteristics of a foil is
the use of adaptive morphology. The flexible fins or wings of some animals enhance indeed

15
Chapter 1. Introduction

Fig. 1.1: Left and middle: Velocity triangles of a VAT, with ω the rotation speed and θ
the azimuth angle. Only the tangential force generates thrust, causing the turbine rotation,
so that most of the forces results in structural loads. Right: CAD model of a three-bladed
H-Darrieus turbine rotor. Extracted from [34].

their propulsion [35, 36]. The deformation modifies the foil shape, affecting the flow and the
pressure field and then the fluid forces, as shown by Fig. 1.2. The angle of attack α (defined

Fig. 1.2: A rigid and a flexible NACA0018 foil (National Advisory Committee for Aeronau-
tics). The large deformations modify significantly the rotor flow characteristics. Extracted
from [37].

in Fig. 1.2) of the blade is varying along the cycle, reaching angle where the flow has stalled,
which induces important lift reduction. Therefore, even a slight trailing edge deflection to-
wards the flow direction would delay the stall and increase the lift. This may avoid deep and

16
 1.1. Motivations: experimental research about chordwise flexible blades turbines

light dynamic stall at high angles of incidence. Thereby, the dynamic shape modification of
turbine blades adapted to the transient flow conditions have been recently considered as a
promising research subject. It can be done actively or passively.
Active deformations can be performed with a control device. Several studies [38, 39, 40]
showed that turbines with blade pitch control system or imposed blade shape deformations
could see their efficiency improved. However, active blade-shape deformation system imposes
a mechanical system requiring power source and maintenance. This system needs also to be
sufficiently robust to resist to temperature and humidity conditions experienced by a pro-
duction turbine, while being manufactured for a reasonable price for large scale production.
Passive blade shape modification is then preferred, using flexible material relatively light
with respect to the fluid to ensure pressure based deformations. This has been achieved
experimentally by Zeiner-Gundersen [41] with a hydrodynamic VAT, resulting with improve-
ments at low current velocities in comparison with rigid model. Also, MacPhee & Beyene [42]
have designed a three flexible blades HAT. In order to quantify the turbine performance, the
power coefficient, a normalized value of torque, is measured for different values of Tip Speed
Ratio (TSR), which is a ratio between the blade velocity (ωR) and the inflow velocity. Good
improvements were found, as given by Fig. 1.3. More recents studies have been carried out by

Fig. 1.3: Power coefficients comparison for rigid and flexible (morphing) blade VAT, obtained
with the experiment of McPhee & Beyene. Extracted from [42].

Hoerner et al. [34, 37, 43, 44] with regard to an oscillating hydrofoil in a closed water channel
(Fig 1.4). The blade is pitched to replicate the variation of angle of incidence experienced
by a rotating VAT blade. The compound hydrofoil consists of a leading edge milled from
aluminum, a skeleton of carbon and a body of silicone. The rigidity of the foil is adjustable
based on the thickness of the carbon core. Parametric studies of rigidity, oscillation frequency
and tip speed ratio have then been performed. The best results combines a reduction of 25%
in structural loads along with improved blade thrust (see Fig. 1.1) of 20% compared to an
identical rigid turbine [44].

17
Chapter 1. Introduction

Fig. 1.4: Photographies of the experimental set up used by Hoerner et al.. Extracted from [34].

Note that blade deformations can also occur in spanwise direction, even for VAT [45],
but the efficiency seems limited compared to rigid blades. This effect is mainly studied for
HAT [46] or flapping foil [47], so that it will not be discussed in this work.
Giving those promising results, one of the goal of the ANR project DYNEOL is the exper-
imental study of the influence of chordwise flexibility on energy harvesting, with for instance
a four blades VAT. However, the physical phenomena occuring with chordwise flexible blade
remain unclear in some aspects. That is the reason why another objective is to develop a
tool to numerically reproduce those experiments in order to provide a better understanding
of the flow and structure behaviours. The goal of this work is then to develop a solver able
to reproduce cases of that kind, involving Fluid-Structure Interaction (FSI). To explain the
underlaying difficulties, the next sections will present the FSI phenomenon, and then will
give a state of the art of numerical methods for FSI simulations. Those sections will aim to
present the requirements for chordwise flexible blades turbines simulations. On that basis,
detailed objectives of this PhD as well as the plan of this manuscript will be provided.

1.2 Numerical methods for FSI simulation

1.2.1 Generalities
FSI phenomena are very common in nature and human activities. It happens as soon
as a fluid and a solid are in contact and that the solid is deformed significantly under the
effect of the forces applied by the fluid. The structure displacement influences also the flow,
resulting in a complex coupling process [48, 49]. Given that the interaction is instantaneous
and continuous, the study of this phenomenon can be difficult. In nature, FSI occur for flora
with tree deformations caused by wind or submerged aquatic canopies [50]. For fauna, the

18
 1.2. Numerical methods for FSI simulation

fins of fishes [35], wings of birds or insects [36, 51] are flexible so that it has bioinspired several
human designs. Two of these examples are illustrated in Fig. 1.5.

Fig. 1.5: Example of FSI in nature with seagrass meadow as an example of dense submerged
canopy from [50] (top) and bumblebee flight thanks to flexible wings (bottom) from [51].

FSI is also studied for human activites. A lot of applications involve flexible structure in-
teracting with complex flows, such that parachutes [52], aerospatial [53], civil engineering [54]
or even biomedecine [55]. This is the phenomenon causing the foil deformations mentionned
in the last section. More examples are given in [56].
In order to provide a better understanding of those problems, numerical simulations of FSI
cases have been developed substantially in the last decade. This has been made possible by

19
Chapter 1. Introduction

recent increases in the available computational ressources, but it is still very challenging since
it requires advanced numerical methods in different physic fields. FSI simulations consist in
coupling fluid and solid solvers. Among other, problem unknowns are the fluid velocity u,
solid displacement d and fluid forces f . The coupling conditions occur at the interface Γ
between fluid and solid as,
∂d

u= (1.1a)


ˆ∂t 
∂u
f= µ + P n dΓ (1.1b)


Γ ∂n
with µ the dynamic viscosity, P the fluid pressure and n the normal direction to Γ pointing
on the solid to the fluid. In fact, solid imposes its velocity to fluid while fluid applies a force
on solid through viscous shear and pressure.
Besides, FSI computation can be very unstable, especially when solid and fluid densities
are close [57], so that fluid inertial effects are as important as solid ones. The resulting
added-mass effect of the fluid on the structure, i.e. the mass of fluid which is accelerated by
the structure, is known to be source of numerical difficulties [58, 59]. Different approaches to
face this problem are presented in the next sections.

1.2.2 Monolithic and partionned approaches

Monolithic approach consists in writting the entire problem in a single non linear sys-
tem [60, 61]. This system is then solved using a Newton-Raphson type method. Conse-
quently, fluid and solid fields are solved in the same time with a fully implicit formulation,
resulting in an very accurate and unconditionally stable solving [62].
However, this method lacks of adaptability and remains difficult to use on complex FSI
cases (geometry, physics) that require highly advanced solvers for both subtasks. Partionned
schemes, where suited techniques can be used in each domain, are thus preferred as recent
studies can demonstrate.
With this method, two different solvers are used for solid and fluid. They have to exchange
data, fluid forces and structure displacements, at each time step. This technique allows a
large variety of solvers for each physic fields. The fluid is generally solved with an Arbitrary-
Lagrangian-Eulerian (ALE) approach, allowing to compute the flow equations on a moving
grid. In addition, turbulence is often modeled with a Unsteady Reynolds Averaged Navier-
Stokes (URANS) approach, in which none of the scales of the turbulent spectrum are resolved
(see section 2.2.1.3 for more details). For the solid, Finite Element Method (FEM) with
special shell and membrane elements are widespread [63], given that FSI often occurs for
thin solid. In most cases, coupling takes place at the interface Γ between fluid and solid, and
the two meshes complement each other. The fluid mesh then has to move according to the
solid displacement, imposing a mesh movement solving to maintain a good cell quality (see
section 4.1). This step can represent a large part of the computational cost of FSI simulations.
Nonetheless, it must also be noticed that Immersed Boundary Methods (IBM) can also be
used for FSI coupling [55, 64, 65, 66]. With this technique, the Lagrangian solid nodes move

20
 1.2. Numerical methods for FSI simulation

over the fixed fluid grid, resulting in an overlap of the fluid and solid domain. Consequently,
the elastic solid forces are interpolated from the lagrangian solid grid to the eulerian fluid
grid. This results in a body forces source term in the Navier-Stokes equations. Conversely,
the velocity is interpolated from the surrounding fluid to the structural nodes. Therefore,
the fluid mesh remains unchanged with this method, reducing the computational cost, but
the interpolations near the interface result in a loss of accuracy [56]. The interpolation of
the velocities from the incompressible fluid to the structure also implies that the structural
motion is divergence free, although not all solids behave this way.
It exists a lot of different coupling schemes, making the computation more or less accurate
and stable. Two main categories can be distinguished; explicit (or weak) and implicit (or
strong) coupling [56]. With weak coupling partitioned techniques, flow and structural equa-
tions are solved separately and only once per time step. This technique results in a reduced
computational cost, but it has been shown that the added-mass effect causes instability of
explicit coupling with an incompressible fluid [58, 59]. This is highlighted in section 4.2.2.
In order to enforce the equilibrium of the traction and displacement on the fluid-structure
interface for cases with important added-mass effect, a strong coupling scheme is required.
This implies several subiterations at each time step since flow and structural equations have
to be solved in each coupling iteration. Different schemes are possible to minimize the num-
ber of subiterations, with for instance fixed-point method [67]. Examples of scheme will be
given in section 4.2.
Note that another approach is also possible where two systems are solved but the added-
mass effect is computed in matrix form and taken into account for acceleration relaxation [68].
Nevertheless, this solution computational cost is growing very rapidly for solid system with
a large number of degrees of freedom.
A large review of the different possible schemes for partitioned simulation is available
in [56]. Now that the different methods for FSI simulation have been presented, next section
will focus on existing techniques for the simulation of chordwise flexible blades turbines.

1.2.3 Numerical simulations of chordwise flexible blades turbines

Even if the present work is focused on energy harvesting applications, flapping or pitching
foils are also studied to optimize propulsion, as detailed in the review of Wu et al. [69]. The
number of papers in this domain is growing since the last two decades (Fig. 1.6), but just 5% of
them include 3D flexible foil approaches, and it is mainly experimental works. It is indeed still
very difficult to numerically reproduce FSI cases of flexible foil with high fidelity. It requires
a FSI solver able to handle cases with turbulence and complex geometries, for the fluid as
well as for the solid. Consequently, it involves a FSI solver with different characteristics:
• For the flow solving, it is essential to correctly capture the velocity gradients close to
the blade, which can be difficult with techniques based on IBM. As explained earlier,
interpolations near the interface are sources of error. Besides, as a fine description of
boundary layer at high Reynolds number requires small cell sizes, ensure the computa-
tion accuracy in cases with large structure displacements involves a fine grid resolution

21
Chapter 1. Introduction

Fig. 1.6: The number of papers on flapping foils since 1970 (top) and the proportion of paper
involving foil according the nature of structure (bottom). Extracted from [69].

in an extented region. The necessary large amount of cells would induce prohibitively
high computational cost. On the contrary, body fitted techniques do not suffer from
these deficiencies, but require a mesh movement computation.

• In regard to the computation of the structure deformation, several methods exist, gen-
erally based on different types of finite element. However, accurate description of the
foil geometry and internal stress, especially in cases of chordwise flexibility, implies a
solid solver that can use 3D solid elements, where FSI solvers sometimes prefer to use
membrane or shell elements [70, 71, 72]. Those elements are highly adapted to predict
the behaviour of a structure with the corresponding geometry, but do not allow to
reproduce structures with more complex geometry, such as a foil.

• Furthermore, flexible blades are composed of material with density comparable with
water density, implying that the FSI coupling can be very unstable and requires a strong
coupling, as explained in the previous section. A monolithic approach can overcome this
issue [60], but this solution lacks adaptability, so that a partitioned scheme is favored in
most recent studies. This method allows to fully benefit from highly advanced solvers
for both fields of application.

22
 1.2. Numerical methods for FSI simulation

These requirements explain the limited number of numerical studies about flexible foil. How-
ever, some FSI solvers have been developed to face this challenge. For instance, MacPhee
et al. [73] managed to simulate a three flexible blades turbine using the URANS approach
with the k-ω-SST model for turbulence. Results of those simulations are given in Fig. 1.7.
Despite small deformation cases, they found an increase of power coefficients around 9.6%
for the blade with a Young modulus E = 0.5MPa in comparison to the same turbine with
rigid blade. Similar method has been used by Marinić-Kragić et al. [74] to reproduce a

Fig. 1.7: Power coefficient comparisons for each VAWT simulation (left) and displacement
magnitude (in meters) of the blade at E = 0.5MPa with a chord length c = 0.4m (right), by
MacPhee & Beyene. Extracted from [73].

Savonius-type VAWT and they also found an additional 8% increase in the power coefficient.
More recently, other simulations have been performed [75, 76, 77], but always with a RANS
approach.
Nevertheless, RANS-based FSI solvers do not always provide reliable results, because
the pressure distribution is not accurate enough to conveniently determine dynamic stall
separation point position. This is why Hoerner et al. [34] suggest that three-dimensional
Large-Eddy Simulations (LES) approach is required to reproduce with high fidelity a high
Reynolds number case involving a chordwise flexible blade.
Unlike RANS approach, with LES the large scales of the flow are resolved while only the
smaller subgrid scales are modelled. However, as detailed in [78], LES-based FSI solvers face
different issues, because they demand higher grid quality than RANS approach and impose
a fine time and space resolution. The moving grid quality can be difficult to maintain, as
it requires robust and efficient mesh movement method, especially for unstructured grids.
Very few LES-based FSI solvers have been developed [78, 79, 50, 80], but they never use 3D
solid elements, making a chordwise flexible foil simulation impossible. To the best of author’s
knowledge, LES-based FSI solvers meeting all the previous requirements do not exist in the
literature so far. Given all these elements, the next sections will explain the objectives of this
PhD.

23
Chapter 1. Introduction

1.3 Objectives and plan

As explained above, the main motivation of this PhD is the numerical study of VAT with
chordwise flexible blades. Given that the tool required for this kind of simulation does not
seem to exist in the literature so far, the main objective of this work is then to develop a
high fidelity FSI solver, based on LES approach, and able to reproduce a wide variety of FSI
configurations, especially cases presented in section 1.1.
All implementations will be carried out within YALES2, a multi-physics library, initially
designed for fluid mechanics. Thus, a partitioned coupling is chosen to provide genericity,
and benefit from the existing ALE solver of YALES2. This library did not initially possess
a solver for structure dynamics prediction though. In order to ensure compatibility between
solvers and to easily modify the FSI solver for future works, it has been chosen to develop a
Structural Mechanics Solver (SMS) from scratch within YALES2. The resulting FSI solver
will then use LES for the flow solving, and FEM with 3D solid elements for the structure
displacement prediction. Both domain will be discretized with unstructured grids to not be
limited by geometry. Furthermore, for mesh movement, a special method will be developed
to meet all the requirements of these simulations. Once the solver will be functionning and
validated, it will be applied to experimental cases.

This manuscript will then start by presenting the existing library and the fluid ALE solver
in Chapter 2. A numerical study of four rigid blades VAT will also be performed. The FEM,
used for the SMS development, will then be explained in Chapter 3 and tests performed
to ensure good validation of the SMS will be detailed. The FSI solver, which couples the
ALE solver and the SMS, will be discussed in Chapter 4, starting with a presentation of
the original mesh movement method developed on purpose. The coupling scheme used will
then be introduced step by step using a simple test case. The developed FSI solver will be
validated against a 2D numerical benchmark with laminar flow in Chapter 5. It will finally be
used to reproduce 3D experimental turbulent cases in Chapter 6; after a validation against a
first turbulent case, a flexible pitching foil experiment is simulated, confirming the potential
of the developed FSI solver for its intended use.

24
 Chapter 2

Arbitrary Lagrangian-Eulerian
solver

This second chapter presents the fluid solver used in this work. First, it introduces the
library YALES2 used for the entire thesis. An additional focus is then given on the Arbitrary
Lagrangian-Eulerian (ALE) solver, and specific numerical methods used by this lattest are
detailed. The Dynamic Mesh Adaptation (DMA) method description is also provided. The
last part of the chapter presents an application case where an experiment of a four blades
VAT is reproduced numerically thanks to the ALE solver.

Contents
2.1 Presentation of YALES2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
2.2 ALE solver . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
2.2.1 Numerical method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
2.2.1.1 Prediction step in ALE . . . . . . . . . . . . . . . . . . . . . . . . 28
2.2.1.2 Correction step . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
2.2.1.3 Large-Eddy Simulation . . . . . . . . . . . . . . . . . . . . . . . . 30
2.2.1.4 Dynamic Mesh Adaptation . . . . . . . . . . . . . . . . . . . . . . 32
2.2.2 Application case: simulation of a Darrieus turbine . . . . . . . . . . . . . . . 34
2.2.2.1 Case presentation . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
2.2.2.2 Numerical setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
2.2.2.3 Results and discussion . . . . . . . . . . . . . . . . . . . . . . . . . 39

2.1 Presentation of YALES2

The library used in this work is YALES2 (Yet Another LES Solver) [81]. It has been
originally designed to solve low-Mach number Navier-Stokes equations (incompressible and
variable density) using LES and DNS approaches with finite volumes. YALES2 can use
unstructured meshes and adaptive grid refinement with a large number of elements while
maintaining good scalability thanks to a Double Domain Decomposition (DDD) strategy. At
first, the mesh is split into sub-parts that are affected to each computational core. Then, the
sub-parts are divided again into cells group of prescribed size, the element groups (ELGRPs),

25
Chapter 2. Arbitrary Lagrangian-Eulerian solver

as illustrated in Fig. 2.1. This technique allows for easily optimizing the use of processor

Fig. 2.1: Double domain decomposition strategy used in YALES2. Extracted from [81].

memory for cache-aware algorithms and may also be exploited by deflation algorithms [82],
in order to fully benefit from high-performance computing on massively parallel machines.
However, use of DDD requires specific data structure. Each ELGRP represents indeed an
independent mesh block, adding the necessity to dispose of a structure to connect elements at
the interface of several ELGRPs, even on the same processor. This internal communicator has
then a copy of data for nodes, pairs and faces at the interface. In a similar manner, external
communicators handle the communication at the interfaces between processors. Finally, a
last structure is used for element located on boundaries to facilitate their specific treatments.
This architecture is represented in Fig. 2.2. This technique is very effective but it complicates
new implementations within the existing code.
Furthermore, YALES2 uses 4th-order central finite-volume method and 4th-order time
integration called TFV4A [83]. The code, initially developed at CORIA (COmplexe de
Recherche Interprofessionnel en Aerothermochimie), is now used and developed in different
laboratories and also used by industrials to reproduce realistic configurations with high fi-
delity. That is why it has now become a multi-physics library able to reproduce two-phase
flows (Lagrangian particles), spray and atomization (Levelset) and combustion (Tabulated or
finite-rate chemistry) cases. Moreover, it disposes of magneto-hydro-dynamics solver and a
granular flow solver. Examples of applications in various physic fields are shown in Fig. 2.3.
The library can be then used in combustion [84], biomechanics [85], process engineering [86],
geophysical flow [87], hydroelectricity [88, 89] or wind turbine energy [90].
In this work, the Arbitrary Lagrangian-Eulerian (ALE) solver has been used to predict
the flow behaviour in a moving domain. This solver was already existing before the present
work and was used to reproduce rotating turbine for example. Thus, developments performed
in the ALE solver have remained then minor, while on the contrary YALES2 did not have

26
 2.1. Presentation of YALES2

Fig. 2.2: Double domain decomposition strategy used in YALES2. Extracted from [81].

Fig. 2.3: Examples of applications with YALES2 in biomechanic, aerodynamic or hydrody-

namic. More details can be found and animations can be visualized at [91].
.

27
Chapter 2. Arbitrary Lagrangian-Eulerian solver

the Structural Mechanics Solver (SMS) which has been thereby entirely developed during the
PhD (see chapter 3). But first, this chapter will present the ALE solver with the next part
discussing physic theory and numerical methodology, and the last one showing an application
case.

2.2 ALE solver

2.2.1 Numerical method
2.2.1.1 Prediction step in ALE
The equations describing the flow dynamics are the Navier-Stokes equations:

 du = − ∇P + ν∆u + f ,

(2.1a)
v
dt ρf

∇ · u = 0, (2.1b)
where u is the vector field of material velocity, P is pressure, ρf is the fluid density, ν is
the kinematic viscosity, and fv a volumetric force. To take into account moving bodies,
the Arbitrary-Lagrangian-Eulerian (ALE) approach is used. It consists in integrating the
governing equations on deformable control volumes, while nodes of the computational mesh
are moved as Lagrangian points. In order to express the modified Navier-Stokes equation in
the ALE framework, the material time derivative has to be defined. If nodes are moved at a
velocity w, then the material time derivative can be written such as
du ∂u
= + (u − w) · ∇u . (2.2)
dt ∂t
Thus, the governing equations can be rewritten as follows

 ∂u ∇P
+ ((u − w) · ∇)u = − + ν∆u + fv , (2.3a)
∂t ρf

∇ · u = 0. (2.3b)
The time advancement scheme used for explanation is an explicit low-storage four-step
Runge-Kutta (LSRK4) scheme [92] recasted in ALE formalism and coupled with the Chorin’s
projection correction method [93] for the pressure term (more details in next section). Full
methodology is detailed in [94, 95]. This method can be broken down into a prediction step,
where the grid is moved, and a correction step. During the prediction step, Eq.(2.3a) is
then integrated in space on a node-centred control volume ω(t) and in time between tn and
tn+1 = tn + ∆t such as
ˆ tn+1 ˆ ˆ tn+1 ˆ
∂
udωdt + ∇ · ((u − w)u)dωdt = RHS
tn ∂t ω(t) tn ω(t)
ˆ tn+1 ˆ (2.4)
un+1 Vn+1 − un Vn + ∇ · ((u − w)u)dωdt = RHS
tn ω(t)

28
 2.2. ALE solver

where Vn is the control volume at tn , and RHS = ν∆u + fv . For the sake of clarity, the RHS
is omitted in the following. The four sub-steps of the time advancement are computed as:
u0 = un
ˆ
Vn αi ∆t
ui = un − ∇ · ((ui−1 − wn+1 )ui−1 )dω for i = 1, ..., 4 (2.5)
Vi Vi ω(ti )

u∗n+1 = u4
where the index i refers to time of the i-th sub-step ti = tn + αi ∆t, u∗n+1 is the predicted
velocity and αi is a coefficient as α = [1/4, 1/3, 1/2, 1]. It should be noted that for w = 0,
Vn = Vi so that classical LSRK4 scheme is recovered. Also, grid nodal velocity w is considered
constant during the timestep, and consequently, node position x is computed at the beginning
of each sub-step with:
x0 = xn
xi = xi−1 + βi ∆twn+1 for i = 1, ..., 4 (2.6)
xn+1 = x4 + βf ∆twn+1
where βi are coefficients linked to αi equal to βi = [1/8, 1/24, 1/12, 1/4]. They are chosen so
that the grid is at the midpoint configuration, meaning that
xn + αi ∆twn+1
xi = . (2.7)
2
At the end of the fourth sub-step, the grid is at xn+1/2 as confirmed by the βi summation.
It needs then to be replaced to reach its final position with βf = 0.5 to finally get xn+1 =
xn + ∆twn+1 .
The reader’s attention is now drawn to the integration volume ω(ti ) in Eq.(2.5). This
volume has to verify:
ˆ
Vi − Vn = −αi ∆t ∇ · wn+1 dw for i = 1, ..., 4. (2.8)
ω(ti )

This relation, known as the Discrete Geometric Conservation Law (DGCL) [96, 97], states
that for each control volume, the volume change between tn and ti must be equal to the
volume swept by the control volume faces during ti − tn . Considering that w is constant
during the entire time step, the integration volume ω(ti ) cannot be equal to Vi ; in order to
satisfy Eq.(2.8), ω(ti ) must be computed at the midpoint configuration for each sub-step,
which corresponds indeed to the position xi given by Eq.(2.6). The RHS is computed with
values taken at ti−1 but on the midpoint configuration as well. Thus, present scheme allows
to compute velocity prediction u∗n+1 , while respecting the DGCL.
It should also be noticed that, in the original Chorin algorithm, the RHS term is com-
puted without the pressure term. This has been possible since the velocity field u can
be decomposed in a sum of a irrotational field and a solenoidal field, as stipulated by the
Helmholtz-Hodge theorem. The pressure term of Navier-Stokes equation is then taken into
account in the projection step, as detailed in the next section.

29
Chapter 2. Arbitrary Lagrangian-Eulerian solver

2.2.1.2 Correction step

The Chorin algorithm [93] can be explained by writing Eq.(2.3a) with a Euler explicit
scheme such as:
un+1 − un ∇Pn+ 1
= −((un − wn+1 ) · ∇)un + ν∆un + fv − 2
(2.9a)
∆t ρf
un+1 − un ∇Pn+ 1
= −((un − wn+1 ) · ∇)un + RHS − 2
(2.9b)
∆t ρf
∇Pn+ 1
un+1 = un − ∆t((un − wn+1 ) · ∇)un + ∆tRHS − ∆t 2
(2.9c)
ρf
∇Pn+ 1
un+1 = u∗n+1 − ∆t 2
. (2.9d)
ρf
The computation of u∗n+1 , i.e. the prediction step, is explained in the previous section.
Taking the divergence of Eq.(2.9d) and considering that un+1 has to satisfy the divergence
free condition (Eq.(2.3b)), the Poisson equation is obtained:
ρf
∆Pn+ 1 = ∇ · u∗n+1 . (2.10)
2 ∆t
This linear system is solved using a Deflated Preconditionned Conjugate Gradient (DPCG)
solver [82]. This step represents the major part of the computational time. Finally, the
velocity correction step can be performed where the newly obtained pressure field Pn+ 1 is
2
used to rectify u∗n+1 with Eq.(2.9d). It must be precised that in YALES2, the method is
slightly different since the pressure gradient computed from Pn− 1 is added to the RHS
2
before the prediction and removed afterwards, so that it is taken into account for the time
advancement sub-steps [98].

2.2.1.3 Large-Eddy Simulation

The resolution of Navier-Stokes equations for turbulent flows can be achieved thanks
to various modeling approaches leading to different computational cost. First, with Direct
Numerical Simulation (DNS), any scale can be resolved. As a result, very accurate flow
predictions are obtained, however the downside is that cells sufficiently small are required to
capture the smallest eddies. Concretely, computational cost of this method appears to be
prohibitive for highly turbulent flows with current computational capacity.
The opposite approach is the Reynolds Averaged Navier-Stokes (RANS) approach, which
consists in applying a time-average operator to the Navier-Stokes equations. The fluctuat-
ing parts of the flow are then completely modeled, and none of the scales of the turbulent
spectrum are resolved. This approach is very widespread and attractive for industries as it
does not require refined meshes but also results in very short restitution times. Nevertheless,
the accuracy is highly sensitive to the validity of the RANS model which depends of the flow
configuration.

30
 2.2. ALE solver

In the present work, an alternative approach between RANS and DNS is used: the Large-
Eddy Simulation (LES). While the large scales of the flow are resolved, the smaller subgrid
scales (SGS) are modelled using a SGS model. The principles of those three approaches are
illustrated in Fig. 2.4.

Fig. 2.4: Sketch of energy density E vs wavelength k in an homogeneous isotropic turbulence

(log-log scale). Note that the smallest eddies finally vanish into heat. Extracted from [99].
.

With LES formalism, the governing equations are the filtered Navier-Stokes equations,

 ∂ui ∂ui ui 1 ∂P ∂ 2 ui ∂τijSGS
+ =− +ν − + fv i , (2.11a)


∂t ∂xj ρ ∂xi ∂xj ∂xj ∂xj

 ∂ui
= 0, (2.11b)


∂xi


where the · operator is used to denote filtered quantities over scales which are large enough
to be resolved with the computational mesh. Subgrid-scale eddies are then unresolved and
result in the residual-stress tensor τijSGS = ui uj − ui uj . This term has then to be modelled
with a SGS model; in this work, the SGS model are based on eddy viscosity model,
 
1 ∂ui ∂uj

1  S ij = +
τijSGS − δij τkk = −2νt S ij with 2 ∂xj ∂xi (2.12)
3
νt = (Cs ∆)2 D


31
Chapter 2. Arbitrary Lagrangian-Eulerian solver

1
where νt is the turbulent viscosity, ∆ = V 3 the filter length, Cs the model constant and D
the time-scale operator. It exists a lot of different eddy viscosity SGS models [100], but in
this work only the Dynamic Smagorinsky model is considered. It gives
q
D = 2S ij S ij (2.13)

while Cs is computed dynamically [101, 102].

The numerical methodology to solve equations governing the fluid dynamics on a moving
grid with ALE approach has now been presented. The ALE solver has already been validated
in [95]. Next part will then explain how to compute the node velocities w in the entire domain
while maintaining a good grid quality.

2.2.1.4 Dynamic Mesh Adaptation

In general, mesh movement is prescribed on at least one domain boundary, like in the case
of a moving object. A mesh movement solver is then required to compute grid movement
in the rest of the domain while maintaining a good cells quality. The computation of nodes
velocity w can be performed by different methods, as exposed in section 4.1.1. In YALES2, w
can be computed with algorithm based on analogy with spring network [103], or more simply,
it can just be prescribed for basic case. For complex cases, a pseudo-solid based method has
been implemented during the PhD, but as this technique uses concepts presented in the next
chapter, it will be explained in detail in section 4.1.
However, for some cases involving very large deformation, moving or rotating object, mesh
movement method alone cannot conserve a mesh describing the changing geometry. The only
solution is then to produce a completely new mesh, adapted to the new geometry i.e. to use
Dynamic Mesh Adaption (DMA). In YALES2, the re-meshing step is based on MMG, the
sequential anisotropic mesh adaptation library for tetrahedral (3D) and triangle element
(2D) [104, 105]. The parallel mesh adaptation strategy proposed by Benard et al. [106] is
used. The benefit of DMA has been already discussed in various configurations including
simulations of multiphase flows [107, 108] or moving bodies [109].
To prevent from performing the grid adaptation procedure at every step, a criterion on
the maximum allowed rate of element deformation is defined based on the element skewness
S. This latter can be determined quantitatively by computing
Vref − Vc
S= (2.14)
Vref
where Vc is the cell volume and Vref is the volume of the equilateral tetrahedron in 3D (or
triangle in 2D), which fits in the same circumsphere as the one of the element (Fig. 2.5). After
each time step, the mesh quality is assessed by computing the maximum skewness Smax in the
computational domain. If this value exceeds a given threshold Slim , the re-meshing step (grid
adaptation procedure) is triggered, and all data fields are interpolated on the new mesh, as
illustrated in Fig. 2.6. The new mesh is generated according to a given metric field computed
from the old mesh.

32
 2.2. ALE solver

Fig. 2.5: Representation in 2D case of a cell (in purple) and of its reference shape (in green).

Fig. 2.6: Example of case with a rotating plate in a cylindrical domain. The mesh movement is
here prescribed (by the Mesh Movement Solver (MMS)) so that a transition zone is established
between the plate in solid rotation and the fixed walls. When the Smax reaches a threshold
value in the transition zone, DMA is then used to generate a new mesh adapted to the new
geometry.

33
Chapter 2. Arbitrary Lagrangian-Eulerian solver

It should also be noted that, in order to limit the cost of the re-meshing step, the grid
adaptation procedure can be applied only on a part of the flow domain. A mask can in-
deed be applied on the elements which are not deformed during the ALE time step. These
elements are without nodes velocity or located in a region with a uniform movement (solid
rotation or translation). This masking procedure allows also to preserve the mesh quality for
specific regions by ensuring no-mesh deformation, as for the boundary layers, for example.
Furthermore, thanks to this technique, MMG just has to deal with a reduced number of
element, given that high density cells region can usually be masked.
Developments were made during the PhD to establish the mask as a function of the
distance R of cells to previously flagged boundaries. This was usefull for cases with complex
geometries. Also, it had been observed that the successive interpolations of the metric field
could induce diffusion, so that the generated meshes finally did not present the same cell size
distribution than the first mesh. A solution developed to face this issue involves recomputing
the desired metric fied after each re-meshing as function of R, avoiding in this way the
diffusion.
Finally, it should be brought to the attention of the reader a last point about data
interpolation between old and new grids. It has been observed that this step resulted generally
in a pressure spike at the next temporal iteration. This effect rapidly vanishes in the next
iterations, so that it does not affect the flow prediction. Following deeper investigations, it
has been concluded that it might be caused by the velocity interpolation. The hypothesis is
that the interpolation step could not preserve the condition ∇ · u = 0 usually ensured by the
correction step (see section 2.2.1.2). The correction step of the following time step should
then offset and this would affect the pressure computation. This effect can be neglected in
most cases, except FSI cases where the pressure is used to predict the structure displacement.
This challenge and path forwards are discussed further in section 4.2.4.
Numerical methods for flow solving on moving grids has now been explained, as well as
techniques to maintain suitable grid quality. Last part of this chapter presents an experimen-
tal case reproduced with the ALE solver using DMA. It shows that the present methodology
can be applied to case of rigid blades VAT.

2.2.2 Application case: simulation of a Darrieus turbine

2.2.2.1 Case presentation

As explained in section 1.1, the ANR project DYNEOL aims to study VAT performances,
experimentally and numerically, for rigid and flexible blades. Numerical simulations of the
experiments of partner laboratories have then been performed during the PhD. One of them,
run by the PPRIME institute, involved a four blades Darrieus turbine in a water channel.
This experimental setup aims to study the influence of inflow velocity u∞ , water height
and Tip Speed Ratio (T SR = ωR/u∞ with R the radius and ω the rotation speed) on the
turbine performances. Once this first study was performed, turbines blades were replaced
by deformable blades in order to compare power coefficient for the two type of blades as

34
 2.2. ALE solver

function of T SR. To be able to reproduce such cases, a FSI solver was required, hence its
development detailed in the present work. However the first step has been the validation of
the ALE solver by reproducing rigid blades experiment. This study is presented in this part.
The inflow velocity u∞ can here be adjusted until u∞ = 5 m.s−1 , as the water height
Hwater . Pictures of the experimental set-up are given in Fig. 2.7. The water channel di-

Fig. 2.7: Darrieus turbine (left) and water channel (right) used in the experiment.

mensions are given in Fig. 2.8. The turbine has a diameter of D = 0.04 m while blades
profile is a NACA0015 with a chord of c = 0.08 m and a lenght of H = 0.4 m. The
turbine center is at 2.2 m from the inlet, and at h = 0.044 m from the bottom of the chan-
nel. The best efficiency operation point was chosen for the simulation, corresponding to
u∞ = 0.85m/s, Hwater = 0.588 m = 1.47H and T SR = 2.15. The chord based Reynolds
number is Re = 68000.

Fig. 2.8: Water channel dimensions. Note that water height (here Hwater = 0.588m) can be
adjusted.

Next part will present the numerical setup used to reproduce this case.

35
Chapter 2. Arbitrary Lagrangian-Eulerian solver

2.2.2.2 Numerical setup

To reproduce numerically a case with a rotating body, different strategies are possible. A
common approach is to use a rotating coordinate system, where the mesh does not need to
be deformed [110]. Nevertheless, water channel side walls cannot be taken into account with
this method, and considering that the turbine fills half of the channel width, the effect of
blockage cannot be neglected. Another approach involves using two different grids, one for
the turbine, the other one for the channel [111, 112, 89]. In particular, it can then be difficult
to ensure mass conservation in flow solving. Moreover, data are interpolated from one grid
to the other for every temporal iteration, which could be source of error.
The chosen method here relies on the DMA. The node velocity w is prescribed in the
entire domain by defining three zones with two parameters, R1 and R2 :

• The first zone (r < R1 ) is the high density cells zone containing the entire turbine. This
zone moves in solid rotation so that w(r) = ω ∧ r where r is the radial vector between
the turbine rotation axis and the considered point, and r = krk. In this zone, the cells
are not deformed.

• The second region (R1 < r < R2 ) is a transition zone where the radial speed is decreas-
ing linearly from w(R1 ) to zero as function of r. Cells deformations then occurs in this
region, until the maximum skewness Smax reaches a chosen threshold value which will
trigger the DMA.

• In the third zone (r > R2 ), the mesh does not move.

The technique is illustrated in Fig. 2.9. Besides, the masking procedure explained in sec 2.2.1.4
is used here, as illustrated in Fig. 2.10. It can be noticed that only the transition zone could
be not masked, but tests showed that giving to MMG all the remaining domain facilitates
its task, making the adaptation step faster.
It is well known that in numerical simulation, turbine performance prediction is highly
dependent of the close blade mesh resolution [31, 89]. A mesh convergence study has been
then performed to make sure that numerical results were mesh independent. In this study,
the mesh strategy to ensure a fine close blade resolution with a minimum number of elements
involves using prism layers, as illustrated in Fig. 2.11. These prismes need to be then cut into
tetrahedrons to be able to use DMA, even if this region is masked. With this technique, three
meshes are performed with different elements size in the prism layers. ∆s is defined as the
lenght of prism face against the blade. Meshes characteristics are summarized in Tab. 2.1.

18M 50M 166M

−2
∆s/c (×10 ) 1.25 0.625 0.3125
Number of elements (×106 ) 18 50 166

Table 2.1: Characteristics of the three meshes.

36
 2.2. ALE solver

Fig. 2.9: Three different mesh regions for the mesh movement prescription.

Fig. 2.10: Side view (left) and top view (right) of the computational domain where the
masked regions are in red: the DMA does not occur in this region. In the present case, the
masked region corresponds to the solid rotation zone.

37
Chapter 2. Arbitrary Lagrangian-Eulerian solver

Fig. 2.11: Close blades mesh (18M). It is first generated with prisms layers before being cut
into tetrahedrons.

Also, two main simplifications are made here. First, the support attached to the superior
disk shown in Fig. 2.7 is not present in the simulation. This would be numerically possible,
but the hypothesis was made that it would not have a strong influence on the turbine perfor-
mances. Furthermore, slip walls are used on every channel walls because the boundary layers
description would be too costly. It also means that for the top side, the free surface is not
taken into account. Once again, this would be numerically doable but it would require the
development of a two-phase flow ALE solver and would increase the computational cost. The
assumption is then made that the free surface does not affect the flow close to the blades.
However, this has been investigated by testing cases where the channel height above the
turbine was increased to reduce the effect of blockage. The results are discussed in the next
part.

38
 2.2. ALE solver

2.2.2.3 Results and discussion

The first study aims at reaching a mesh convergence to make sure that numerical results
are mesh independent. The case is then reproduce on the three meshes 18M, 50M and 166M.
In order to quantify the power extracted by the turbine from the flow, a power coefficient Cp
is defined such as
ωkFf luid ∧ rk
Cp = (2.15)
ρf rHu3∞
where ρf = 1000kg/m3 and Ff luid is the fluid force applied on the entire turbine surface
Sturbine . This latter is computed as a sum of the viscous shear and pressure forces with
ˆ  
∂u
Ff luid = µ + P n dS (2.16)
Sturbine ∂n

where µ = 10−3 P a.s is the dynamic viscosity, P the pressure and n the unite normal vector.
The results for the three meshes are given in Fig. 2.12, where the power coefficient is plotted
as function of the normalized time. It can be seen very clearly that mean power coefficient

Fig. 2.12: Cp as a function of normalized time. Here τ is characteristic time such as τ =

2π
. As nblades = 4, a complete turbine rotation corresponds to 4τ .
ωnblades

needs to be computed once the flow is established. It is also shown that the power coefficient

39
Chapter 2. Arbitrary Lagrangian-Eulerian solver

presents a τ periodicity. Another periodicity of τ2 ≈ D/u∞ can also appears on this kind of
configurations, corresponding to the travel time of inter blades eddies. In fact, the upstream
blade can generate structures that will collide with downstream blade and disrupt the flow.
This effect is not observed here, possibly because the grid resolution in the turbine center
might not be sufficient to transport these structures. Time averaged power coefficients have
then been computed for each mesh on the three last rotations of the turbine as presented
in Tab. 2.2. These values show that the solution seems converged at the 4th turbine rota-

2nd turn 3rd turn 4th turn

(4τ - 8τ ) (8τ - 12τ ) (12τ - 16τ )
hCp i on 18M 0.068 0.013 0.020
hCp i on 50M 0.234 0.124 0.127
hCp i on 166M 0.331 0.300 0.307
expe 0.356

Table 2.2: Time averaged power coefficient computed on succesive turbine rotations.

tion given that they are closed of those obtained at the 3rd rotation. Despite an error of
13.8% obtained with the finest mesh on the last rotation, it can be deduced that values are
converging towards the experimental values as more refined mesh were used.
However, the first value given by the PPRIME institute was 0.21 instead of 0.356. It
was coming from a technical problem that had been fixed at the end of the PhD. The
overestimation of Cp with the last mesh was then requiring more understanding: that is the
reason why another study was also performed with the 166M mesh to investigate the influence
of the water height Hwater as function of the T SR. For this purpose, two other meshes had
been generated, based on the mesh 166M, but with respectively 1H and 3H added above the
initial channel of 1.47H (H being the blades lenght defined in section 2.2.2.1). As explained
in the previous part, the idea here is to reduce the effect of blockage due to the channel
top wall, considering that the experiment was executed with a free surface. Results are
given in Fig. 2.13. This shows that for T SR < 2.5, considering a slip wall instead of a
free surface seems acceptable. On the contrary, for T SR = 3, the free surface seems to
have a strong influence on the flow. This can be explained by the higher rotation speed ω
used experimentally to increase the T SR. With turbine support also rotating, that resulted
in strong deformation of the free surface so that the upper disk of the turbine was nearly
not immersed anymore. To confirm this assumption, it would be interesting to reproduce the
experiment again with T SR = 3 but with a lower ω and a lower u∞ . Velocity fields computed
with the three different geometries for T SR = 3 are given in Fig. 2.14. It can indeed be seen
that the high velocity region above the turbine computed with Hwater = 1.47H progressively
vanishes when Hwater is increasing. Also, the velocity difference between each side of the

40
 2.2. ALE solver

0.3

0.2
Cp

0.1

166M
0.0 166M+1H
166M+3H
−0.1
expe
1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0
TSR

Fig. 2.13: Power coefficient as function of T SR for different water channel height. The
experiment is performed with a free surface of height correspond to the mesh 166M (Hwater =
1.47H).

Fig. 2.14: Velocity field obtained the three different channel heights Hwater =
[1.47H, 2.47H, 4.47H]), front view just before the turbine (left) and side view (right), for
T SR = 3.

41
Chapter 2. Arbitrary Lagrangian-Eulerian solver

turbine is reduced, explaining the power coefficient drop observed in Fig. 2.13.

The ALE solver used for the FSI coupling has now been presented. As explained at the
beginning of this chapter, this solver was already operational before this work. In addition,
the VAT simulation confirms that this solver is effective and can be used for FSI coupling.
However, implementations in this solver have been successfully performed, especially regard-
ing restarts and lagrangian probes, and few bugs were fixed. Several developments also
concerned the DMA and were mentioned in section 2.2.1.4. Next chapter will now present
the solver used to predict the structure behaviour. Unlike the fluid solver, the SMS has been
entirely developed from scratch during the PhD.

42
 Chapter 3

Structural Mechanics Solver

This chapter presents the solid solver. Given that it has been developed from scratch, the
Structural Mechanics Solver (SMS) constitutes one of the main achievement of this work. The
numerical methods used are introduced step by step, starting by a stationnary linear elasticity
problem of one linear finite element. The problem complexity then progessively increases,
to end with non linear dynamic problem solved with quadratic elements. In the second part
of the chapter, the solver accuracy is verified with several test cases, ensuring the suitable
implementation of the previously described method.

Contents
3.1 Numerical method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
3.1.1 Generalities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
3.1.1.1 History of FEM . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
3.1.1.2 Shape functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
3.1.2 Stationary problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
3.1.2.1 Stationary linear elasticity problem of one finite element . . . . . 47
3.1.2.2 Displacement boundary conditions . . . . . . . . . . . . . . . . . . 50
3.1.2.3 Matrices assembly . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
3.1.2.4 Stationary non-linear problem . . . . . . . . . . . . . . . . . . . . 54
3.1.3 Dynamic problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
3.1.3.1 Generalities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
3.1.3.2 Newmark scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
3.1.3.3 Generalised-α scheme . . . . . . . . . . . . . . . . . . . . . . . . . 60
3.1.4 Use of reference space - Quadratic finite elements . . . . . . . . . . . . . . . 62
3.1.4.1 Principle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
3.1.4.2 Implementation in the SMS . . . . . . . . . . . . . . . . . . . . . . 66
3.1.4.3 Forces computation . . . . . . . . . . . . . . . . . . . . . . . . . . 67
3.1.5 Possible improvements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
3.2 Verification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
3.2.1 Comparison with Ansys . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
3.2.1.1 Stationary problem . . . . . . . . . . . . . . . . . . . . . . . . . . 71

43
Chapter 3. Structural Mechanics Solver

3.2.1.2 Dynamic problem . . . . . . . . . . . . . . . . . . . . . . . . . . . 73

3.2.2 Mesh convergence study with the implemented elements . . . . . . . . . . . 81
3.2.2.1 2D finite elements . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
3.2.2.2 3D finite elements . . . . . . . . . . . . . . . . . . . . . . . . . . . 84

3.1 Numerical method

3.1.1 Generalities
3.1.1.1 History of FEM
Finite Element Method (FEM) is a general numerical method for solving Partial Differ-
ential Equations (PDE) by subdividing a system into smaller parts called finite elements.
The method was originally designed to face complex problems in civil and aeronautical en-
gineering implying structural analysis. The process of evolution leading to the present FEM
is complex, and it involves a lot of mathematicians and physicist such that it is difficult to
quote a precise invention date, but it can be estimated in the 1960s. Actually, it has been
independently developed at different places in the world by pioneers such as Hrennikoff [113]
and Courant [114] in north America, Argyris [115] and Zienkiewicz [116] in Europe and Kang
in China. Mathematicians had back then already developed general techniques to solve PDE
such as finite difference and on the other hand, engineers used to discretize continuum prob-
lems by more simple elastic bars. FEM arised at the frontier of these two fields, as it is shown
in Fig. 3.1 that briefly sums up the historical development of the method. More details can
be found in [117].
The use of FEM has since been generalized in a wide varety of physics fields like heat
transfer, electromagnetism and fluid dynamics. As it is the most widespread method used for
structural mechanics code, it was a fairly obvious choice to use it to develop the Structural
Mechanics Solver (SMS), even though there were no solvers in YALES2 using FEM. The
work presented in this chapter has been mainly based on books from Zienkiewicz [118, 119],
on the documentation of the code ANSYS [120], but also on the work of Debard [121].

3.1.1.2 Shape functions

FSI computations require a solid solver able to solve dynamic problems with finite de-
formations. Nevertheless, this implies several complex numerical methods, and this chapter
aims at providing the entire methodology used. For the sake of clarity, the FEM will then
be explained step by step. This part will start by introducing a crucial concept of FEM,
shape functions. FEM is based on the mesh discretization of a continuous domain into a set
of discrete sub-domains, called elements. At first, only one finite element is considered. The
main idea in order to determine a continuous field u(x, y, z) inside the element consists in

44
 3.1. Numerical method

Fig. 3.1: History of development of FEM. Extracted from [118].

using interpolation functions and values of u at the elements nodes. In a two-dimensional

space, inside an element e of nn nodes, u(x, y) can be approached with:
nn
X
u(x, y) ≈ Nae (x, y)ua , where ua = u(xa , ya ) (3.1)
a=1

and Nae is the interpolation function of node a in element e. These functions are called shape
functions and verify the two following properties:

Nae (xa , ya ) = 1

(3.2a)
Nae (xb , yb ) = 0, a 6= b (3.2b)

with b another node in the element e.

In the case of plane stress, if u(x, y) now represents the displacement of a point at (x, y),
the vector can be written:  
u(x, y)
u(x, y) = , (3.3)
v(x, y)
and shape functions under matrix form become

Nea = Nae I. (3.4)

45
Chapter 3. Structural Mechanics Solver

Eq.(3.1) gives for the displacement u at any point within the element e:
 

 u1  
 v1 
   e  
 
e
u(x, y) N1 (x, y) 0 ... Nnn (x, y) 0

= ... (3.5)
v(x, y) 0 N1e (x, y) ... 0 Nnen (x, y) 
unn 
 

 
 
vnn


and with FEM formulation,

 
 u1 
... u2 = Ne ue .

u = Ne1 Ne2 (3.6)
...
 

The problem consisting in solving the continuous displacement field u(x, y) has then been
discretized since there are now only 2nn unknowns i.e. the components of displacements of
the element nodes.
For the sake of clarity, until section 3.1.4, the method will be explained for linear elements
such as 4-nodes tetrahedron or 3-nodes triangle.

Fig. 3.2: Illustration of the displacement field discretization on a 3-nodes triangle. Extracted
from [121].

For this last element, illustrated in Fig. 3.2, shape functions are linear and can be written:

Na = (aa + ba x + ca y)/(2∆) (3.7)

in which
 a1 = x2 y3 − x3 y2 (3.8a)


b1 = y 2 − y 3 (3.8b)

c1 = x 3 − x 2 (3.8c)

46
 3.1. Numerical method

with other coefficients obtained by cyclic permutation of the subscripts in the order 1, 2, 3,
and where ∆ is the area of triangle obtained with
1 x1 y1
2∆ = det 1 x2 y2 . (3.9)
1 x3 y 3
The space discretization method using finite element has been presented, but the methodology
to compute nodes displacements ue in function of external constraints and material properties
has to be described. The next part will then present the governing equations, which have
to be solved in order to predict the behavior of one element for a stationary linear elasticity
problem.

3.1.2 Stationary problems

3.1.2.1 Stationary linear elasticity problem of one finite element
In this part, only one finite element is considered. Generalization to the whole domain
will be done in section 3.1.2.3. The equilibrium equation governing the solid dynamics can
be written in a Lagrangian frame such as :

σji,j + bi = ρüi (3.10)

where σij are components of Cauchy stress, corresponding to the internal force experienced
by the element, bi are body force components and ρ refers to the solid density. Note that the
partial derivatives are here and in the sequel denoted by
∂f ∂f
f,i = and f˙ = . (3.11)
∂xi ∂t
For stationary problem, Eq.(3.10) gives

σji,j + bi = 0 (3.12)
which corresponds to the balance between external and internal forces. To describe the
deformation state of the solid, a strain-displacement relation needs to be defined; the strain
may be expressed in tensor form as
1
εij = (ui,j + uj,i ) (3.13)
2
for small deformations, i.e.
|εij |  1 and |ωij2 |  kεij k (3.14)
where ωij denotes the small rotations defined by
1
ωij = (ui,j − uj,i ). (3.15)
2

47
Chapter 3. Structural Mechanics Solver

The last equation to describe the behavior of a material is the constituve equation relating
the stress to the strain. The simplest model is the linear elasticity; by neglecting strains
arising from other sources than the displacement, the stress-strain relation gives

σij = Cijkl εkl (3.16)

where Cijkl are elastic moduli which depend on material properties (Eq.(3.26)).
In matrix notation, through the symmetry of the tensors, the three-dimensional forms of
the stress can be written such as
 T
σ = σ11 σ22 σ33 τ12 τ23 τ31 (3.17)

where the shear stress is expressed

τ12 = σ12 = σ21 . (3.18)

In similar manner, strain is given by

 T
ε = ε11 ε22 ε33 γ12 γ23 γ31 (3.19)

where the engineering shear strain is introduced such as

γij = εij + εji = 2εij , i 6= j. (3.20)

Both Eq.(3.13) and Eq.(3.19) then provide for the matrix form of the strain-displacement
relation
ε = Bu (3.21)
where the operator B is defined by
 
∂ 0 0 ∂ 0 ∂
 ∂x1 ∂x2 ∂x3 
BT =  0
 ∂ 0 ∂ ∂ 0 .

(3.22)
∂x2 ∂x1 ∂x3
∂ ∂ ∂
 
0 0 0
∂x3 ∂x2 ∂x1
Using Eq.(3.6), strain-displacement can now be expressed in function of node displacement
ue as
Xnn nn
X
e
ε = Bu = BNa ua = Bea ua = Be ue (3.23)
a=1 a=1

where Bea is the strain matrix of element e at node a given by

 e e e

Na,1 0 0 Na,2 0 Na,3
Bea T =  0 Na,2 e
0 Na,1 e e
Na,3 0 . (3.24)
e e e
0 0 Na,3 0 Na,2 Na,1

48
 3.1. Numerical method

In matrix form, Eq.(3.16) gives

σ = Dε. (3.25)
For an orthotropic elastic material, the elasticty matrix D is given by
Ey
 
Ex (1 − ν 2 Ez ) (ν + ν ν Ez ) Ez (ν + ν ν ) 0 0 0
 h yz E
y h xy xz yz
Ey h xz yz xy

 Ey E E E E E 
2
 (νxy + νxz νyz )
h Ey
z x (1 − ν z
Ex ) z (ν yz + νxz νxy
Ey
z ) 0 0 0 
h xz h 
Ez (ν + ν ν ) Ez (ν + ν ν Ey ) Ez (1 − ν 2 Ey )
 
D=
 h xz yz xy 0 0 0 
 (3.26)
 h yz xz xy
Ex h xy E
x 

 0 0 0 Gxy 0 0 

 0 0 0 0 Gyz 0 
0 0 0 0 0 Gxz

where Exi is the Young modulus material in the xi direction, νxi xj the Poisson’s ratio and Gxi xj
is the shear moduli. The Young modulus corresponds to the stiffness of the material while
the Poisson’s ratio is related to the deformation of the material in directions perpendicular
to the specific direction of loading. In the last equation, h is introduced as

2 Ey 2 Ez 2 Ez Ez
h = 1 − νxy − νyz − νxz − 2νxy νyz νxz . (3.27)
Ex Ey Ex Ex
Besides, Gxi xj is computed by
E xi
Gxi xj = . (3.28)
2(1 + νxi xj )
For plane stress or plane strain formulation of D, please refer to [122].
Now, Eq.(3.12) can be written with FEM formulation:

BT σ + be = 0. (3.29)

By multiplying by Ne T and integrating over the element Ωe , it gives:

ˆ ˆ
eT
B σdΩ + NeT be dΩ = 0. (3.30)
Ωe Ωe

Substituting Eq.(3.25) and Eq.(3.23), it comes

ˆ ˆ
eT
B DεdΩ + NeT be dΩ = 0 (3.31a)
ˆ Ωe
ˆ Ωe
BeT DBe dΩue + NeT be dΩ = 0. (3.31b)
Ωe Ωe

For the linear stationary problem, it is finally obtained:

Ke ue = f e (3.32)

49
Chapter 3. Structural Mechanics Solver

where ˆ
e
K = BeT DBe dΩ (3.33)
Ωe

is the linear stiffness matrix of the element and

ˆ
e
f =− NeT be dΩ (3.34)
Ωe

is the force vector.

At this point, it is interesting to focus on elements near the boundary. Displacement
boundary conditions will be detailed further in section 3.1.2.2, but if an element face Γe is
subjected to a surfacic force t, a loading term on the nodes of the element has to be added
such as ˆ
fe → fe − NeT tdΓ. (3.35)
Γe

Also, it should be noted that for linear elements, derivatives of linear shape functions are
constants, such that the strain matrix Be is constant. It results in a strain and a stress
constant throughout the element. Futhermore, integration required for computation of Ke
(Eq.(3.33)) can be simplified as follows:
ˆ
e
K = BeT DBe dΩ = BeT DBe V e (3.36)
Ωe

with V e the element volume. Finally, the reader can notice that the stress and the strain are
not explicitly present in Eq.(3.32) since they can be considered as post-processing data for
linear elastic problem, where small deformations hypothesis are made.
Most structural problems include displacement boundary conditions, to fix the structure
as well as to impose a deformation. These conditions require a specific method, detailed in
the next section.

3.1.2.2 Displacement boundary conditions

Prescribing a displacement condition to a node can be implemented in several ways.
Taking the example of a one dimensional linear stationary problem, for a 2-nodes segment it
can be written:     
K11 K12 u1 f1
= (3.37)
K21 K22 u2 f2
in which the displacement u1 = u1 is to be imposed. For this purpose, different techniques
exist as listed below.

• The simplest method to implement is a “penalty” approach where the system to solve
is     
α K12 u1 αu1
= (3.38)
K21 K22 u2 f2

50
 3.1. Numerical method

whith α  K11 . This method, sometimes called “Method of large number” can be
easily added in a code, but presents the main drawback to not be perfectly accurate
and lacks of elegance.
• Another method that requires more connectivity, involves modifying the system such
as     
1 0 u1 u1
= . (3.39)
K21 K22 u2 f2
However, if K is symmetric, it has to become
    
1 0 u1 u1
= . (3.40)
0 K22 u2 f2 − K21 u1
This solution, much more rigorous, is more difficult to implement in a code since it
requires enough cells connectivity to be able to retrieve K21 from a boundary node 1,
even if node 2 is not on the boundary.
• A third approach is to modify the system as above and then to remove from the system
all equations where displacement has to be imposed. For a specified displacement of
node i, the i-th column and the i-th line will then be removed. In the example it gives:

K22 u2 = f2 − K21 u1 . (3.41)

This method is optimal since it leads to a system with a minimum number of unknwows,
but requires an important flexibility of the code.
The method is applied to a one-dimensional example, but principle remains the same for
other dimensions; length of vector displacement becomes ndim nn and it is then possible to
impose displacement only for one component. Considering the different methods requirements
and the data structure of YALES2, the second method has been used for the SMS. It has
allowed to implement options such that the user could impose zero displacement, rotation,
vibration or just a given displacement to a chosen domain boundary.
The methodology to formulate a linear stationary structural problem using finite element
has now been presented on one subdivision of the domain, one finite element. The assembly
process has to be explained to be able to write equations for the entire system.

3.1.2.3 Matrices assembly

It is now mandatory to define a system where the unknowns will be the components of
nodal displacements such as
u1 
 


 .. 
  
 .   ui 

 

u= ui where ui = vi (3.42)


 .
.. 

 wi
 

 

unodes
 

51
Chapter 3. Structural Mechanics Solver

is the displacement vector of node i, and ui , vi , wi are its components along each space
dimension. The vector u will then list the displacement of the nnodes nodes of the whole
structure in which all the ne elements participate. To write the corresponding system, both
Equations (3.33) and (3.34) have now to be written for the whole system Ω such as
ˆ
K= BT DBdΩ (3.43)
Ω

and ˆ ˆ
T
f =− N bdΩ − NT tdΓ. (3.44)
Ω Γ

These matrices are in fact composed by submatrices computed by

ne
X ne
X
Kab = Keab and fa = fae . (3.45)
e=1 e=1

This results in the system to solve for a linear stationary problem:

Ku = f . (3.46)

This very simple rule for matrice assembly allows to compute the submatrice on each
element with Eq.(3.33) and Eq.(3.34) while storing results directly in the appropriate global
data. Material properties can then be different for each element, allowing the simulation of
multi-material structure.
Also, it may be noticed that for two nodes a and b that do not belong to the same
element, Kab = 0. Concretely, submatrices Kaa and fa can then be indexed on nodes number,
and submatrices Kab with a 6= b can be indexed on pair numbers. Furthermore, for the
application cases targeted for this PhD, the matrice K will always be symmetric. It has
then been chosen to store only half of submatrices to save computational memory. However,
even if K is symmetric, for two or three dimensionals problem submatrices Kab 6= Kba but
Kab = KTba . This detail complicated a lot matrices computation, linear solving operations
and displacement condition imposition. Nevertheless, data structure of YALES2 have been
proven sufficiently flexible to overcome those difficulties so that the SMS has been written
using those techniques.
The full methodology to write and solve a linear elastic stationary problem using linear
finite elements has been presented. At the start of the PhD, a beam solver has then been
entirely written in Python to confirm the good understanding of the FEM. It used 3-nodes
triangle to solve stationary and dynamic problems, and good agreement were found with
theorical results and Ansys comparison, as shown in Fig. 3.3. Its development constituted a
first step before writting the SMS which had to fit to the specific YALES2 data structure.
Nevertheless, the presented method does not allow to reproduce cases with finite defor-
mation. It is now necessary to explain how a structural problem can be solved whithout
resorting to the small deformation assumption.

52
 3.1. Numerical method

Fig. 3.3: Example of simulation computed with the Python solver, comparison with Ansys
code.

53
Chapter 3. Structural Mechanics Solver

3.1.2.4 Stationary non-linear problem

In section 3.1.2.1, it has been assumed that deformations were remaining small so that
linear relations could be used to represent the body strain. Besides, shape functions were
computed and integrated always on the same geometry. If large deformation are now consid-
ered, it is mandatory to distinguish between initial and deformed configurations, respectively
called reference and current configuration, as illustrated Fig. 3.4.

Fig. 3.4: Reference and current configurations for finite deformation problem. Extracted
from [119].

As much as possible uppercase letters and indices will refer to quantities defined in the
reference configuration and lowercase letters and indices to quantities defined in the current
deformed configuration. The displacement vector can then be used to change between two
coordinate frames such as

xi = δiI (XI + UI ) (3.47)

where Einstein notation and Kronecker symbol δ are used. It should be noted that displace-
ment can be used equally in the two frame since u = U .
It is now possible to define the deformation gradient F by
∂ui
FiI = δiI + = δiI + ui,I . (3.48)
∂XI
Stress and strain need also to be redefined for large deformation. In the reference config-
uration, the Green strain tensor EIJ can be used, expressed in function of the displacement
such as
1h i
EIJ = UI,J + UJ,I + UK,I UK,J (3.49)
2

54
 3.1. Numerical method

where linear part corresponds to the small strain form (Eq.(3.13)).

In a similar manner, the second Piola-Kirchhoff stress SIJ can be related to the previously
used Cauchy stress σij with
det(F )σij = FiI SIJ FjJ . (3.50)
The simplest stress-strain relation is the hyperelastic Saint-Venant-Kirchhoff model given
by:
SIJ = CIJKL EKL , (3.51)
where CIJKL is the same constant elastic moduli already used in Eq.(3.16). This model has
been the only one implemented in YALES2 during the PhD, because it was sufficient to
reproduce all the desired cases. However, data and solver structure had been developed to
facilitate future implementations of other material models. For that matter, given that some
material model results in non-bijective stress-deformation relation (i.e. path dependent), the
possibility to apply progressively forces with successive sub-steps has been implemented, for
stationary or dynamic problem.
Concretely, the computation of the derivatives of displacement is performed by using
Eq.(3.1) such as
nn
X ∂N e a
ui,I = uai ≡ Na,I
e
uai . (3.52)
∂XI
a=1

The strain operator also needs to be changed; using the previous small deformations
b e can be introduced with
problem definition of Bea (Eq.(3.24)), B a

b e = Be + Be N L
B (3.53)
a a a

where the nonlinear part is given by

e e e
 
u1,1 Na,1 u2,1 Na,1 u3,1 Na,1
e e e

 u1,2 Na,2 u2,2 Na,2 u3,2 Na,2 

e e e
u 1,3 N u 2,3 N u 3,3 N
Bea N L
 a,3 a,3 a,3

u1,1 Na,2 + u1,2 Na,1 u2,1 Na,2 + u2,2 Na,1 u3,1 Na,2 + u3,2 Na,1  .
= (3.54)
e e e e e e

 
e
u1,2 Na,3 e e e e e 
+ u1,3 Na,2 u2,2 Na,3 + u2,3 Na,2 u3,2 Na,3 + u3,3 Na,2
e e e e e e
u1,3 Na,1 + u1,1 Na,3 u2,3 Na,1 + u2,1 Na,3 u3,3 Na,1 + u3,1 Na,3

As a result, the strain-displacement relation is now non linear.

Finally, Eq.(3.30) gives for a stationary non linear problem:
ˆ
b T SdΩ = f .
B (3.55)
Ω

To solve this non linear equation, a Newton-Raphson method is used so that the problem
becomes the minimization of the function

55
Chapter 3. Structural Mechanics Solver

ˆ
Ψ(u) = f − b T SdΩ = f − P(u).
B (3.56)
Ω
To the first order, this equation can be approximated as
!i
∂Ψ
Ψi+1 = Ψ(ui+1 ) ≈ Ψ(ui ) + dui = 0 (3.57)
∂u

with i the Newton iteration counter and dui a displacement increment. The tangeant stiffness
matrice is here defined by
∂Ψ
KT = − (3.58)
∂u
for a stationary problem. One Newton iteration consists then in solving

KiT dui = Ψi (3.59)

to find dui knowing that

i
X
i+1 i i 0
u = u + du = u + duk . (3.60)
k=1

The principle of the Newton method is illustrated in Fig. 3.5. The method requires the

Fig. 3.5: Illustration of the Newton-Raphson method.

computation of KT which is given by

ˆ ˆ bT
bTD ∂B ∂f
KT = B b T BdΩ
b + SdΩ − = KM + KG + KL . (3.61)
Ω Ω ∂u ∂u

56
 3.1. Numerical method

The first term, the material tangeant KM , is computed with D b T = D as the Saint-Venant-
Kirchhoff is used. Besides, the tangent term KG is called the geometric stiffness and can be
written ˆ
Kab
G = Gab I where Gab = Na,I SI,J Nb,J dΩ. (3.62)
Ω
The last term KL comes from the loading which depends on deformation state, e.g follow-
ing forces like pressure. In general, its computation results in an asymmetric term; this
complicates greatly the solving of the final linear system, so that this term is neglected for
the computation of KT . However, this does not affect the final result since forces f i are
systematically updated for the computation of Ψi .
The process then repeats until

Ψ(ui+1 ) ∞
≤ εnewton (3.63)

where εnewton is a selected convergence criterion. It should also be noted that displacement
boundary conditions just need to be imposed at the first Newton iteration.
This iterative method is essential when the solution process is path dependent, like for
some non linear constitutive equations of solids. It can be improved with different techniques,
such as quasi-Newton methods or line search procedures [119], to accelerate the convergence
and reduce the number of iterations required. Also, a formulation using only the current
configuration is possible. This formulation was originally used for the SMS, but it leads to
stability issues when it comes to simulations involving a lot of iterations.
The method to solve linear and non linear structural mechanics stationary problem using
FEM has now been detailed. Nonetheless, for FSI applications, the SMS needs to be able
to predict the dynamics of the structure, at any time, under the influence of variable loads.
This last assertion complicates the previous equations, and specific temporal schemes have
to be used for time advancement; this is explained in the following part.

3.1.3 Dynamic problems

3.1.3.1 Generalities
Using the d’Alembert principle, intertial forces through the body can be introduced as:
ˆ
f → f − NT ρNdΩü. (3.64)
Ω

For a dynamic problem, Eq.(3.46) now changes into

Mü + Ku = f (3.65)
where M is the mass matrix defined by
ˆ
Meab = e
Mab I and e
Mab = Nae ρe Nbe dΩe . (3.66)
Ωe

57
Chapter 3. Structural Mechanics Solver

Furthermore, for non linear problem, the function Ψ becomes

ˆ
Ψ(u) = f − B b T SdΩ − Mü. (3.67)
Ω

Concretely, the definition given by Eq.(3.66) corresponds to the so-called “consistent” mass
matrix, but it can be very useful to obtain a diagonal mass matrix for explicit solving for
instance. The mass can then be physically lumped at the element nodes so that
nn ˆ
X
e
Maa = e
Nae ρe Nbe dΩe and Mab = 0 when a 6= b. (3.68)
b Ωe

Consistent or lumped mass matrices both have to respect the mass preservation condition
given by
Xnn X
nn ˆ
e
me = Mab = ρe dΩe (3.69)
a b Ωe

where me is the total mass of the element. It should be noted that as me remains constant
during the deformation, the mass matrix only needs to be computed once at the start of the
simulation whereas KT and f have to be systematically updated.
A damping matrix can also be added to previous equation such as

Mü + Cu̇ + Ku = f . (3.70)

In fact, the determination of the damping matrix C for a given material is difficult, and a
method called Rayleigh damping is often used. This consists in considering that

C = αM + βK (3.71)
where α and β are coefficients determined experimentally [123, 124]. The mass coefficient α
results in decaying damping effects on higher modes while the stiffness coefficient β induces
a damping on lowest modes. Also, the stiffness matrix K used to compute C, becomes KT
for non linear solving and can be updated during the simulation or considered as constant.
This choice has been implemented as an user’s parameter into the SMS.
Explicit solving has not been discussed at this point because timestep constraint results
in smaller timestep than in CFD, making it unusable for FSI applications. Equations cor-
responding to dynamic problems have now been introduced, and different temporal schemes
can be used for solving. The two schemes that have been implemented during the PhD will
now be presented in the two following parts.

3.1.3.2 Newmark scheme

The Newmark method is widely used in dynamic analysis. For a linear dynamic problem
at time tn+1 , system can be written

Mün+1 + Cu̇n+1 + Kun+1 = fn+1 . (3.72)

58
 3.1. Numerical method

This problem with three unknowns ün+1 , u̇n+1 and un+1 needs to be simplified to be solved.
The Newmark method consists in considering that velocity and acceleration can be expressed
such as
u̇n+1 = a1 (un+1 − un ) − a4 u̇n − a5 ün (3.73a)
ün+1 = a0 (un+1 − un ) − a2 u̇n − a3 ün (3.73b)
where: 
1
a0 = (3.74a)


α∆t2





 δ
a1 = (3.74b)


α∆t





 1
 a2 = α∆t (3.74c)



1
 a3 = −1 (3.74d)


 2α
δ


a4 = − 1


 (3.74e)



 α !
δ


 ∆t
 a5 = 2 α − 2 . (3.74f)




α and δ are Newmark’s integration parameters. Eq.(3.72) can then be written

   
a0 M + a1 C + K un+1 = fn+1 + M a0 un + a2 u̇n + a3 ün
  (3.75)
+ C a1 un + a4 u̇n + a5 ün

and once this last is solved to find un+1 , velocity and acceleration can be computed using
Eq.(3.73).
The amount of numerical dissipation of the algorithm can be controlled by δ value. How-
ever, to ensure an unconditionally stability, Newmark’s parameters have to meet the following
requirements: 
1
δ≥


 (3.76a)

 2
!2
1 1
α≥ +δ . (3.76b)




 4 2

In a similar manner than in Ansys, a user parameter, the amplitude decay factor γ, is
introduced such as 
1
δ = +γ (3.77a)


2



1 2
α = 1 + γ (3.77b)
4




γ ≥ 0. (3.77c)


59
Chapter 3. Structural Mechanics Solver

Case with γ = 0 then has no numerical damping.

For non linear dynamic solving, the residual Ψ of Eq.(3.56) is now computed with

Ψin+1 = fn+1
i
− P(uin+1 ) − Müin+1 − Cu̇in+1 (3.78)

and one Newton iteration involves solving

 
a0 M + a1 C + KT in+1 duin+1 = Ψin+1 . (3.79)

Initial displacement u0n+1 , velocity u̇0n+1 and acceleration ü0n+1 are taken as the converged
solution from the last time step.
As explained above, the amout of numerical dissipation can be controlled by γ. However,
in low frequency modes the Newmark method fails to retain the second-order accuracy. To
overcome some drawbacks as identified for the Newmark method, another temporal scheme
has also been implemented and is presented in the next section.

3.1.3.3 Generalised-α scheme

This method was originally developed by Chung & Hulbert [125]. It allows to sufficiently
damps out spurious high-frequency response via introducing controllable numerical dissipa-
tion in higher frequency modes, while maintaining the second-order accuracy. It consists in
solving
Mün+1−αm + Cu̇n+1−αf + Kun+1−αf = fn+1−αf (3.80)

where


 ün+1−αm = (1 − αm )ün+1 + αm ün (3.81a)

 u̇n+1−α = (1 − αf )u̇n+1 + αf u̇n

(3.81b)
f


 un+1−αf = (1 − αf )un+1 + αf un (3.81c)

n+1−αf = (1 − αf )fn+1 + αf fn .
 f

(3.81d)

αm and αf are user’s parameters allowing to control the algorithm properties. Using Eq.(3.73),
the system can be written
 
a0 M + a1 C + (1 − αf )K un+1 = (1 − αf )fn+1 + αf fn − αf Kun
   
+ M a0 un + a2 u̇n + a3 ün + C a1 un + a4 u̇n + a5 ün

(3.82)

60
 3.1. Numerical method

with 
1 − αm
a0 = (3.83a)


α∆t2




(1 − αf )δ


a1 = (3.83b)


α∆t




 a2 = a0 ∆t


 (3.83c)

1 − αm
a3 = −1 (3.83d)


 2α
(1 − αf )δ


a4 = −1 (3.83e)






 α !
∆t δ


 a5 = (1 − αf ) −2 .


 (3.83f)
 2 α

It should be noted that for αm = αf = 0 this scheme is identical to the Newmark method.
To ensure that the scheme is unconditionally stable and second-order accurate, α and δ are
computed with 
1
δ = +γ (3.84a)


2





 1 2
 α= 1+γ (3.84b)
4
γ = αf − αm (3.84c)




1


 αm ≤ αf ≤ .

 (3.84d)
2
This optimal case can be described by defining αm and αf in terms of the spectral radius ρ∞
such as
2ρ∞ − 1 ρ∞
αm = , αf = . (3.85)
ρ∞ + 1 ρ∞ + 1
This ensures an optimal combination of high-frequency and low-frequency dissipation [125].
Thus, by selecting the appropriate values for αm and αf , the user can use different well-
known time integration methods, as illustrated in Fig. 3.6. In a similar manner, non linear
problems give:

Ψin+1 = fn+1−α
i
f
− P(uin+1−αf ) − Müin+1−αm − Cu̇in+1−αf (3.86)

and one Newton iteration consists in solving

 
a0 M + a1 C + (1 − αf )KT in+1 duin+1 = Ψin+1 . (3.87)

Full methodology to reproduce non linear dynamic cases has now been provided. This
work was used to develop a first version of the SMS. However, tests showed that this approach
was not suitable to reproduce 3D large deformation cases. In fact, the code was only able
to use 3-nodes triangles or 4-nodes tetrahedrons, which are linear elements. This results

61
Chapter 3. Structural Mechanics Solver

Fig. 3.6: Different time integration methods according values of αm and αf . Extracted
from [125].

in a prohibitive number of elements necessary to predict accurately complex deformations.

Consequently, second order elements, as well as quadrilateral and hexahedron, had to be
available in the SMS, driving to a lot of mandatory modifications in the data structure.
Their implementation is explained in the next section.

3.1.4 Use of reference space - Quadratic finite elements

3.1.4.1 Principle
As showed in section 3.1.1.2, 3-nodes triangles and 4-nodes tetrahedrons have linear shape
functions. That facilitates considerably any integration and derivation necessary to compute
stiffness matrix by example. Nevertheless, this is not true for other elements, where stress
and strain are not constant in the element. Furthermore, for higher order elements, nodes
are added at the pair centers and/or face centers. For second order elements, faces can then
deform in a quadratic shape, making the metric direct computation of the element difficult,
as illustrated in Fig. 3.7.
The widespread solution used in most FEM solvers to deal with this problem is using

62
 3.1. Numerical method

Fig. 3.7: Geometric deformation of a 6-nodes triangles. Extracted from [121].

a reference space or also called isoparametric map. The idea here is to compute data on
one undeformed element in a suited coordinate system, and then to use an isoparametric
transformation to obtain values of the element of the real space. Concretely, for a real space
(x, y, z), a reference space (ξ, η, ζ) can be defined, as illustrated Fig. 3.8. Shape functions
Niref can be determined in the reference space such as

 nn
 x(ξ, η, ζ) = X N ref (ξ, η, ζ)x

(3.88a)


 i i

 i=1

nn


 X
y(ξ, η, ζ) = Niref (ξ, η, ζ)yi (3.88b)


 i=1



 Xnn
z(ξ, η, ζ) = Niref (ξ, η, ζ)zi



 (3.88c)

i=1

where nn is the number of nodes of the element, x(ξ, η, ζ), y(ξ, η, ζ), z(ξ, η, ζ) are coordinates
of a point inside the element e in the real space and (xi , yi , zi ) are coordinates of the i-th
node of the element e in the real space. These relations allow then to determine the position
of a point in the real space from its position in the reference space. By using these shape
functions, it is now possible to define a Jacobian matrix of the transformation:
 nn nn nn

ref ref ref
  X ∂Ni X ∂Ni X ∂Ni
∂x ∂y ∂z  xi yi zi 
 ∂ξ ∂ξ ∂ξ   i=1 ∂ξ ∂ξ ∂ξ
 
i=1 i=1

   n n n

n n n
ref ref ref
  X 
e
 ∂x ∂y ∂z   ∂Ni X ∂Ni X ∂Ni 
J (ξ, η, ζ) = 
  = xi yi .
zi  (3.89)
 ∂η ∂η ∂η   i=1 ∂η ∂η ∂η

i=1 i=1

   n 
 ∂x ∂y ∂z  X n
ref
nn
ref
nn
ref 

 ∂Ni X ∂Ni X ∂Ni
∂ζ ∂ζ ∂ζ  xi yi zi 
∂ζ ∂ζ ∂ζ
i=1 i=1 i=1

63
Chapter 3. Structural Mechanics Solver

Fig. 3.8: Illustration of elements in real space (x, y, z) and corresponding elements in reference
space (ξ, η, ζ). Extracted from [118].

This transformation has to be bijective: the Jacobian determinant has then to remain positive
within the element. That implies that reference element and real element have to be numbered
in the same way. In addition, usual rules of partial differentiation give

∂Niref ∂Niref ∂x ∂Niref ∂y ∂Niref ∂z

= + +
∂ξ ∂x ∂ξ ∂y ∂ξ ∂z ∂ξ

∂Niref ∂Niref ∂x ∂Niref ∂y ∂Niref ∂z

= + + (3.90)
∂η ∂x ∂η ∂y ∂η ∂z ∂η

∂Niref ∂Niref ∂x ∂Niref ∂y ∂Niref ∂z

= + +
∂ζ ∂x ∂ζ ∂y ∂ζ ∂z ∂ζ

64
 3.1. Numerical method

giving relations between derivatives of the two space:

   
∂Niref ∂Niref ∂Niref ∂Niref
   

 ∂ξ   ∂x   ∂x   ∂ξ 
   
       
ref  ref 
   
ref  ref 
 ∂Ni  ∂Ni  and  ∂Ni  = Je−1  ∂Ni  .
   
  = Je    (3.91)
 ∂η  
 ∂y



 ∂y

  ∂η 
       
ref ref
  ref ref
 
 ∂Ni  
∂Ni
 
∂Ni
  ∂Ni 
∂ζ ∂z ∂z ∂ζ

Jacobian can also be used for integration; the computation of the element volume V e in the
real space can indeed be done from the reference space with
ˆ ˆ
e
V = dxdydz = det(Je (ξ, η, ζ))dξdηdζ. (3.92)
V Vref

This relation allows to compute an integration in the real space of a function f such as
ˆ e ˆ
f (x, y, z)dxdydz = f (x(ξ, η, ζ), y(ξ, η, ζ), z(ξ, η, ζ)) det(Je (ξ, η, ζ))dξdηdζ. (3.93)
V Vref

However, the explicit computation of the right hand side of this last equation can be difficult,
especially when f and J are non linear. That is why it is usually performed numerically
by quadrature; the most accurate methods for polynomial expressions is Gauss-Legendre
quadrature [126]. Gaussian quadrature integrates a function as

ˆ 1 ng
X  d2ng f 
f (ξ)dξ = f (ξj )wj + O (3.94)
−1 j=1
dξ 2ng

where ξj are the points where the function is evaluated, wj is the corresponding weight,
and ng the number of points used for integration. Thus, this formula integrates exactly a
polynomial of order 2ng − 1. Examples of points location and weights are given in Tab. 3.9.
This can be directly used for integration over segment element in the reference space. For
more complex elements, shape functions and gauss points location and weights were found
in Code Aster documentation [127]. Examples are given for quadrilateral and hexahedron in
Appendix A (extracted from [127]).
Now that the principle of the use of a reference space has been explained, its application
in the case of structural mechanics and its implementation in YALES2 will be detailed in the
next part.

65
Chapter 3. Structural Mechanics Solver

Fig. 3.9: Gaussian quadrature points and weights corresponding to Eq.(3.94). Extracted
from [119].

3.1.4.2 Implementation in the SMS

Using a reference space, equations of previous sections have to be adapted. Thus, stiffness
matrix (Eq.(3.33)) computation gives:
ˆ
e
K = BeT DBe dΩ
ˆΩe
= Be (ξ, η, ζ)T DBe (ξ, η, ζ)) det(Je (ξ, η, ζ))dξdηdζ
Ωref (3.95)
ng
X
= Be (ξj , ηj , ζj )T DBe (ξj , ηj , ζj )) det(Je (ξj , ηj , ζj ))wj
j=1

∂Niref ∂Niref ∂Niref

where derivatives of shape functions , and required to determine the
∂x ∂y ∂z
strain operator Be are computed using Eq.(3.91). In a similar manner, mass matrix (Eq.(3.66))

66
 3.1. Numerical method

gives: ˆ
e
M = NeT ρNe dΩ
Ωe
ng (3.96)
X T
ref ref e
= N (ξj , ηj , ζj ) ρN (ξj , ηj , ζj )) det(J (ξj , ηj , ζj ))wj .
j=1

However, in practice, it is difficult to store Je (ξ, η, ζ), Je (ξ, η, ζ)−1 , det(Je (ξ, η, ζ)) or shape
functions derivatives because these matrices are element dependent and include high order
polynomial functions. For Nref , which corresponds to the reference element and is unique for
all the domain, shape function polynomial coefficients were stored but this solution was not
possible for the other matrices. Choice has then been made to store directly det(Je (ξj , ηj , ζj ))
 ∂N ref (ξ , η , ζ ) ∂N ref (ξ , η , ζ ) ∂N ref (ξ , η , ζ ) 
i j j j i j j j i j j j
and , , values on Gauss points locations.
∂x ∂y ∂z
Furthermore, new elements had to be added in YALES2 which is originally a finite-volume
based library. To be able to use the existing data structure of the code, the chosen devel-
opment approach was to add nodes at pairs or face centers once the initial mesh were read.
Consequently, both first order and second order were available for the same mesh file. This
implementation required though a lot of modifications in cells connectivity, particularly to fit
with the massively parallel characterictics of YALES2. Besides, the high number of pairs of
some element (351 for 27-nodes hexahedron) complicated these tasks and required important
memory space for few elements. Consequently, optimal number of elements by CPU and by
element group revealed to be significantly lower for these elements than with linear elements.
All elements finally made available in the SMS are shown in Fig. 3.10 and Fig. 3.11. This cor-
responds to a small case with few elements to highlight the parabolic deformation of second
order elements. Efficiency of these elements is further checked in section 3.2.2.
The last point not mentionned in this section is the computation of volumic and surfacic
forces. Yet, it has to be adapted to the use of reference space. This is the subject of the next
part.

3.1.4.3 Forces computation

Volumic forces require a computation of the element volume, dependent of its deformation
state. Such forces will then impose at each time step a computation of det(Je (ξj , ηj , ζj )). In
fact, the resulting force on an element e subjected to a volumic force fV is given by:
ˆ
e
f = NeT fV dΩ
Ωe
ng (3.97)
X T
= Nref (ξj , ηj , ζj ) fV det(Je (ξj , ηj , ζj ))wj .
j=1

As explained above, Jacobian computation allows integration over the entire element,
but is not adapted to describe the metric of only one face of a deformed element. However,

67
Chapter 3. Structural Mechanics Solver

Fig. 3.10: From top to bottom: 3-nodes triangles, 6-nodes triangles (left) and 4-nodes quadri-
laterals, 8-nodes quadrilaterals, 9-nodes quadrilaterals (right).

Fig. 3.11: From top to bottom: 10-nodes tetrahedrons, 4-nodes tetrahedrons (left) and 27-
nodes hexahedrons, 20-nodes hexahedrons, 8-nodes hexahedrons (right).

68
 3.1. Numerical method

surface or normal vector of such face are mandatory when surfacic forces are applied on one
boundary of the structure, which is precisely what happens in FSI cases. The solution is then
to refer to the reference element of inferior dimension ef ace corresponding to one face of the
entire element e, as illustrated in Fig. 3.12.

Fig. 3.12: Use of reference space to compute faces metric. Here the 3-nodes segment ef ace
(left) corresponds to the reference element of the face of the real 6-nodes triangle e (right).
Extracted from [121].

In this example, shape functions N i of ef ace can be introduced such as:

nn nn
X X
x(ξ) = N i (ξ)xi and y(ξ) = N i (ξ)yi . (3.98)
i=1 i=1

It can then be deduced:

n n
X ∂N i (ξ) 
dx(ξ) = xi dξ = Jx (ξ)dξ (3.99a)
∂ξ
i=1
n n
X ∂N i (ξ) 
dy(ξ) = yi dξ = Jy (ξ)dξ (3.99b)
∂ξ
i=1
p q
ds(ξ) = dx(ξ) + dy(ξ) = Jx (ξ)2 + Jy (ξ)2 dξ = Js (ξ)dξ.
2 2 (3.99c)
The force resulting in the application of a surfacic force fS along ef ace can then be com-
puted with:
ˆ ˆ 1
e eT
f = N fS dS = N(ξ)T fS Js (ξ)dξ
ef ace −1

Xng (3.100)
T
= N(ξj ) fS Js (ξj )wj .
j=1

69
Chapter 3. Structural Mechanics Solver

Normal vector at position ξ and surface S of ef ace can also be computed with:
∂y(ξ) 
 

 
   ˆ ˆ 1 ng
∂ξ Jy (ξ)
 X
n(ξ) = = and S = ds = Js (ξ)dξ = Js (ξj )wj . (3.101)
∂x(ξ)  −Jx (ξ) ef ace −1
−

 
 j=1
∂ξ
In a similar manner, for a 3D real element, normal vector at position (ξ, η) is given by:
∂x(ξ, η)   ∂x(ξ, η) 
   

   



 ∂ξ 




 ∂η 



 ∂y(ξ, η)   ∂y(ξ, η) 
 
n(ξ, η) = ∧ dξdη. (3.102)

 ∂ξ   
 ∂η  
   



 ∂z(ξ, η) 


 


 ∂z(ξ, η) 



∂ξ ∂η
The entire methodology used to develop the SMS during the PhD has now been presented.
Despite the existing data structure of YALES2, many things had to be adapted since the code
has not been originally designed for finite elements. The SMS has been then developed with
the aim of reproducing FSI cases only with YALES2 solvers, but its structure is prepared
to facilitate future implementations, like other material models. Given the limited time of
the PhD, some axes of improvements were not investigated but they might be the subject of
future work. The next section aims to present these axes.

3.1.5 Possible improvements

As shown in next section, the SMS has proven to be able to reproduce accurately structural
mechanic cases. However, one of its main axis of improvements remains the reduction of
computational cost. Moreover, a quick analysis of the computation time repartition revealed
that most of the time was spent in the linear system solving algorithm. In fact, this part
has not been discussed yet, but to solve systems such those of Eq.(3.46) or Eq.(3.87), the
SMS uses a Preconditioned Conjugate Gradient (PCG) algorithm. This method will not be
detailed here, but any interested reader is invited to refer to [128]. For small problem, direct
linear solving is possible but those techniques require large memory storage space and suffer
from poor scalability. Iterative linear solving has then been chosen. Different methods were
already developed in YALES2, however the current routines were not suitable for the specific
form of structural mechanics system so they had to be rewritten to fit with the SMS. In
addition, YALES2 offered only the possibility to use a Jacobi preconditionning (JCG), where
for a system (3.46), the matrix MJacobi used for preconditionning is


MJacobi = diag Kii . (3.103)

Nevertheless, as explained in [129], JCG suits well to CFD problem as they are known to
have a better conditionning that CSM problem. The use of more complex and efficient

70
 3.2. Verification

preconditionner are necessary to be able to simulate efficiently problem with a large number
of elements with the SMS. Issues in linear solving with JCG can also arise from solid composed
by material with very different properties, or from very thin geometries. These limitations
were caught during the SMS development, although implementation of new preconditonners
appeared to be challenging due to parallelism constraints. One relevant choice providing
many opportunities would be to make available the use of the PETSc (Portable, Extensible
Toolkit for Scientific Computation) [130, 131, 132] library by the SMS. This latter is indeed
used by code ASTER [127] for instance. Further work would be needed on this subject, thus
contributing to accelerate significantly the linear system solving.
As mentioned in section 3.1.2.4, another possible improvement would be the implemen-
tation of more advanced techniques than the basic Newton algorithm used presently for non
linear solving [119]. Finally, the SMS does not allow sereral different types of element in
the same mesh. However, this strategy can be useful to reproduce some geometries with a
minimum number of elements. For this purpose, some code modifications are required which
would induce more complex data structure but would allow a greater flexibility.
Now that the numerical methodology used by the SMS has been detailed, next sections
of this chapter will aim to validate this solver, without FSI coupling for the moment. In
this way, comparison with Ansys and mesh convergence studies of the implemented finite
elements are presented in the next parts.

3.2 Verification
3.2.1 Comparison with Ansys
3.2.1.1 Stationary problem
In order to ensure the correct implementation of the methodology presented in the first
part of this chapter, the library Ansys is used for comparison. Given that the structural
solvers of Ansys used FEM, similar results are expected. However, YALES2 data structure
is different and some details of the algorithms used by Ansys remain unclear, e.g for linear
system solving.
One of the main challenge faced for verification is the wide variety of options available
through the SMS. Among them, there are 10 different finite elements, the stationary or
dynamic solving, the linear or non linear solving (i.e. the small deformation assumption, see
section 3.1.2.4), the different boundary conditions, the possibility to add structural damping
(Eq.(3.70)), isotropic or orthotropic material, different temporal schemes, the possibility to
use heterogeneous material properties, etc... Considering the high number of combinations,
this part will give an overview of main tests performed for verification. Cases with 2D
elements will not be treated in this section as a 2D structural test case will be presented in
section 5.3. Furthermore, the choice to use 3D solid elements does not restrain the geometry
selection, which would not be the case with membrane or shell elements. A blade geometry
has then been chosen for tests since the PhD aims at reproducing flexible blade turbine.

71
Chapter 3. Structural Mechanics Solver

A NACA0018 profil with a chord C = 66 mm has been selected, corresponding to the foil
used in the experiments run by Hoerner & al. [34, 44]. For the material, the entire blade is
considered as made of silicon which is a simplification of real case where some parts of the
blade are composed of aluminium and carbon fiber. Silicon properties are E = 0.64 GPa,
ν = 0.48 and ρs = 1 100 kg.m−3 .
The first test introduced in this section is a stationary problem of span-wise flexibility. A
blade with a length L = 4C is meshed with 80 112 tetrahedrons as shown in Fig. 3.13.

Fig. 3.13: Blade mesh used for stationary problem.

The blade is fixed at its extremities along z-axis, and an acceleration field a = 20 m.s−2
is applied in y direction. The problem is solved with linear or non linear method, for 4-
nodes tetrahedrons (4NTET) and 10-nodes tetrahedrons (10NTET). With linear method,
the problem is to solve Eq.(3.46), while with the non linear method, the problem is given
by Eq.(3.55). Results obtained with the SMS are shown in Fig. 3.14. It highlights that
problems where the structure length does not remain constant with large deformation cannot
be solved with a linear approach. On the contrary, if only one blade extremity is fixed,
results between linear and non linear solving are similar. This is explained by the constant
stiffness matrix K used for linear solving, while physically the structure stiffness depends on
its deformation state. This is taken into account with the non linear solving during the step
of the Newton-Raphson method, resulting in a more realistic deformation.
Simulations are also reproduced with Ansys with the exact same mesh, and the normalized
minimum vertical displacements of each case are gathered in Fig. 3.15.
A good agreement is obtained for each cases, especially for linear solving. Differences for
non linear solving are slightly larger but remain neglictible (<2.3 %). This can be explained

72
 3.2. Verification

Fig. 3.14: Stationary problem solved with 10-nodes tetrahedrons and linear (top) or non
linear (bottom) method. Front view (left) and back view (right). Structures are colored by
normalized vertical displacement.

Fig. 3.15: Results obtained (minimum of dy /C) for the four different cases and comparison
with Ansys.

by the complexity of algorithm used for non linear solving, where several linear systems have
to be solved successively during Newton-Raphson method, while only one system is solved
with linear solving. It can also be noticed that higher order elements use results in a larger
bending, indicating that the mesh is not fine enough to use 4NTET.
The SMS has now been validated for stationary problems, and the following section will
detail the validation for dynamic problems.

3.2.1.2 Dynamic problem

The idea is now to consider a case of chordwise flexibility. For that, the blade lenght
has been divided by 4 so that L = C and the number of tetrahedrons is reduced at 20 209.
Besides, the entire leading edge, the first third of the foil, is fixed. A force of Fmax = 5

73
Chapter 3. Structural Mechanics Solver

N is then applied on the top surface of the blade, resulting in a curved deformation state
(Fig. 3.16), similar to those observed during the experiments [44].

Fig. 3.16: Dynamic problem solved with 10-nodes tetrahedrons and non linear method (case
E).

In order to amplify the oscillations amplitude, the density is set to ρs = 55 000 kg.m−3 .
This is not realistic but will induce an oscillating deformation state, more interesting for
comparison. This will also highlight the effect of structural damping on the result. Five
cases are then performed with different combination of options. This is detailed in Fig. 3.17.

Fig. 3.17: Dynamic cases parameters. Structural damping is computed with α = 0.2 and
β = 0.2 (Eq.(3.71)).

Once again, simulations are also carried out using Ansys and the quantity used for compar-
ison is the normalized infinite norm of displacement. The force is applied with an increasing
slope of 2 seconds. It is then maintained during 2 seconds and released within 2 seconds.
After t = 6 s, forces are no longer applied to the blade and its displacement are only caused
by intertial effects. A characterictic time τ = 2 s is defined. The Newmark temporal scheme
is chosen with the amplitude decay factor γ = 0 to avoid numerical damping. The condition
of Eq.(3.76) being respected, the scheme is unconditionally stable so that the time step can
be chosen as ∆t = τ /100 s. Results obtained with YALES2 and Ansys are plotted in function

74
 3.2. Verification

of t/τ and are compared for each case in the following, after which a conclusion is given on
the overall agreement between the two solvers.

75
Chapter 3. Structural Mechanics Solver

• Case A is simulated with 4NTET and linear method; results are given in Fig. 3.18. As no
structural damping is applied, the problem involves solving Eq.(3.65). A good agreement
between simulations is found and the exact same deformation state is predicted. Thanks
to the large value of density ρs , the structure oscillates around its equilibrium position, and
the same results are found between 1t/τ and 2t/τ while minor differences appear for the
oscillating state between 3t/τ and 4t/τ . This can be explained by an accumulation over time
of small discrepancies.

YALES2 LINEAR 4NTET ANSYS LINEAR 4NTET

3.0
0.175
Force
0.150 2.5

0.125
2.0
kdk∞/C

F/Fmax
0.100
1.5
0.075
1.0
0.050

0.025 0.5

0.000 0.0
0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0

t/τ

Fig. 3.18: Comparison of the normalized infinite norm of displacement computed with
YALES2/Ansys for case A.

76
 3.2. Verification

• Case B is simulated with 4NTET, linear method and structural damping; results are given in
Fig. 3.19. The problem involves now solving Eq.(3.70), and the Rayleigh damping parameters
defined in Eq.(3.71) are selected as α = 0.2 and β = 0.2. With Newmark scheme, the system
is given by Eq.(3.75). In this case, no differences on the computed structure displacement
can be observed between the two simulations.

YALES2 LINEAR 4NTET DAMPED ANSYS LINEAR 4NTET DAMPED

0.175 3.0
Force
0.150 2.5
0.125
2.0
kdk∞/C

F/Fmax
0.100
1.5
0.075

0.050 1.0

0.025 0.5

0.000 0.0
0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0

t/τ

Fig. 3.19: Comparison of the normalized infinite norm of displacement computed with
YALES2/Ansys for case B.

77
Chapter 3. Structural Mechanics Solver

• Case C is simulated with 4NTET and non linear method; results are given in Fig. 3.20. With
Newmark scheme, one Newton iteration is to solve Eq.(3.79) with C = 0. The oscillations are
predicted with slight differences on the frequency and the mean value. Also, the amplitude
predicted with Ansys after t/τ = 3 is nearly twice superior to the one computed with the
SMS. As for the stationary case, non linear solving implies a larger number of operations
than linear solving, increasing the accumulation of differences. Nonetheless, this relates to
normalized infinite norm of displacement and the blade global predicted behaviour is the
same with the two solvers.

YALES2 NON LINEAR 4NTET ANSYS NON LINEAR 4NTET

3.0
0.175 Force
0.150 2.5

0.125 2.0
kdk∞/C

F/Fmax
0.100
1.5
0.075
1.0
0.050

0.025 0.5

0.000 0.0
0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0

t/τ

Fig. 3.20: Comparison of the normalized infinite norm of displacement computed with
YALES2/Ansys for case C.

78
 3.2. Verification

• Case D is simulated with 4NTET, non linear method and structural damping; results are
given in Fig. 3.21. The differences observed are of the same order of magnitude than those
observed between equilibrium position of the oscillating state obtained in case C. It remains
then small, so that it can be concluded that structural damping seems to have the same
impact in both simulations.

YALES2 NON LINEAR 4NTET DAMPED ANSYS NON LINEAR 4NTET DAMPED

3.0
0.175
Force
0.150 2.5

0.125
2.0
kdk∞/C

F/Fmax
0.100
1.5
0.075
1.0
0.050

0.025 0.5

0.000 0.0
0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0

t/τ

Fig. 3.21: Comparison of the normalized infinite norm of displacement computed with
YALES2/Ansys for case D.

79
Chapter 3. Structural Mechanics Solver

• Case E is simulated with 10NTET, non linear method and structural damping; results are
given in Fig. 3.21. In comparison with case D, it can be observed that the predicted deforma-
tion is more important, indicating again that the mesh is not fine enought to reproduce with
accuracy this case with 4NTET. Furthermore, the difference is slightly increasing compared
to the previous case, which was expected due to the higher number of nodes and the larger
deflection. Nevertheless, a good agreement is found between the two computations.

YALES2 NON LINEAR 10NTET DAMPED ANSYS NON LINEAR 10NTET DAMPED

0.25 3.0
Force
2.5
0.20

2.0
kdk∞/C

0.15

F/Fmax
1.5
0.10
1.0
0.05
0.5

0.00 0.0
0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0

t/τ

Fig. 3.22: Comparison of the normalized infinite norm of displacement computed with
YALES2/Ansys for case E.

Although larger differences have been obtained for case with non linear solving in com-
parison with linear solving case, discrepancies remain minors so that the overall agreement
can be considered satisfying. This confirms the correct implementation in YALES2 of the
methodology presented in the chapter, while checking on the effective functioning of the SMS.
Another test case in 2D is also proposed in section 5.3.

Next section will present a mesh convergence study performed to confirm the efficiency
of the implemented quadratic finite elements in comparison with the linear elements.

80
 3.2. Verification

3.2.2 Mesh convergence study with the implemented elements

3.2.2.1 2D finite elements

The goal of this part is to verify the correct efficiency of the different finite elements that
can be used with the SMS. With this aim, a basic case of a beam embedded at its extremities
is reproduced. In order to focus on elements precision, the problem is considered stationary
and the linear method is used. Furthermore, the convergence criterion of the PCG is chosen
sufficiently low to not have any influence on the computed result. The beam dimensions are
L × h × w with h = w = 0.5 m and L = 20h such as its cross-section is a square of side h. A
vertical surfacic force of FS = −4.8 N.m−2 is applied on the top face of the beam along y-axis
and the problem consists in computing its stationary deformation state. Material properties
are chosen to ensure a large deformation; E = 7.1 MPa and ν = 0.4. Note that the material
density ρs has no influence here since it is a stationary problem without acceleration field.
This first part presents the test for the 2D case. Five finite elements are studied: 3-
nodes triangle (3NTRI), 6-nodes triangle (6NTRI), 4-nodes quadrilateral (4NQUAD), 8-
nodes quadrilateral (8NQUAD) and 9-nodes quadrilateral (9NQUAD). Six meshes composed
of quadrilaterals and six meshes with triangles are then generated with different cell sizes,
s = [0.500, 0.250, 0.125, 0.100, 0.050, 0.025] m, as shown in Fig. 3.23.

Fig. 3.23: Quadrilateral meshes used for mesh convergence study. The dimension is h = 0.5
m. From top to bottom, cell size s is respectively 0.5 m, 0.25 m, 0.125 m, 0.1 m, 0.05 m and
0.025 m.

81
Chapter 3. Structural Mechanics Solver

Computed solution obtained with 4NQUAD are given in Fig. 3.24. It highlights the
progressive mesh convergence given that the result obtained with the two finest meshes seems
identical and is very different from the solution computed with the coarser mesh.

Fig. 3.24: Displacement norm obtained with different quadrilateral meshes (4NQUAD).

The infinite displacement norms computed with the different meshes are then gathered
in Fig. 3.25. A convergence can also be noticed; every finite element seem to give results
converging on the value obtained with the finest mesh of 9NQUAD. This value is taken as

Fig. 3.25: Results obtained (b = kdk∞ in meter) for meshes composed by different finite
elements of different sizes for the 2D case.

reference bending bref , and an error related to this value can be computed for a mesh with cell

82
 3.2. Verification

of size s as b(s) − bref . The normalized error is finally plotted in function of the normalized
cell size (with sref = 0.025 m) on Fig. 3.26.

3NTRI 8NQUAD
6NTRI 9NQUAD
4NQUAD

10−1
b(s)−bref
bref

10−3

100 101
s/sref

10−1
b(N )−bref
bref

10−3

102 103 104

Nelem

Fig. 3.26: Normalized error obtained with different 2D meshes. Here, bref = b(sref ) is
computed with 9NQUAD, and sref = 0.025 m.

It can clearly be seen that for a same cell size, higher order elements (6NTRI, 8NQUAD
and 9NQUAD) give result closer to the reference than lower order elements (3NTRI and
4NQUAD). Besides, 6NTRI and 8NQUAD seem identically efficient, but it must be noticed
that the same geometry requires more triangles than quadrilaterals to be meshed with the
same cell size. That is why the second graph shows the normalized error in function of
the number of elements Nelem . It highlights that for the same value of Nelem , second order
elements lead to better results than first order ones. Also, 8NQUAD are more precised than
6NTRI, which makes sense giving the larger number of degree of freedom per element.
Therefore, this study confirms the efficiency of the second order 2D finite element imple-

83
Chapter 3. Structural Mechanics Solver

mented in the SMS. These elements allow then to describe correctly the structure behaviour
with a lower Nelem , i.e. for a lower computational cost in comparison with first order ele-
ments. They will be used in Chapter 5. Nevertheless, this study needs also to be carried out
for 3D elements, which is the subject of the next section.

3.2.2.2 3D finite elements

The same case is then extruded with w = h so that the beam cross-section is a square.
The elements to study are: 4-nodes tetrahedron (4NTET), 10-nodes tetrahedron (10NTET),
8-nodes hexahedron (8NHEX), 20-nodes hexahedron (20NHEX) and 27-nodes hexahedron
(27NHEX). Once again, 6 meshes are generated with tetrahedrons and 6 meshes with hexa-
hedrons with the same cell size s than for the 2D case. The results are gathered in Fig. 3.27.

Fig. 3.27: Results obtained (kdk∞ in meter) for meshes composed by different finite elements
of different sizes for the 3D case.

Unlike previous cases, all values seem to converge to the one obtained with the finest
mesh of 10NTET sref . The value computed with this mesh is then taken as bref and the
normalized error is plotted in Fig. 3.28.
It can be noticed that even if for the same cell size s the 10NTET gives the smallest
errors, the number of elements required to describe the beam is way larger than with hexahe-
dron. Therefore, for a same Nelem , accuracy is thereby proven to be better with hexahedron.
Nonetheless, it has to be precised than the computational cost per element is higher with
hexahedron than with tetrahedron due to the higher number of nodes considered. Yet, for
the same s, the cost is still lower with hexahedrons. Besides, tetrahedrons are very useful to
represent complex geometry that cannot be discretized with hexahedrons. Given the error
obtained with 4NTET and 8NQUAD, this study also confirms the efficiency of second order
3D finite element.

Those two studies have demonstrated the benefit to use higher order finite element. That
has completed the verification of the correct implementation of the FEM in YALES2. The
development of the SMS has then been a success, ensuring the availability of a performing

84
 3.2. Verification

4NTET 20NHEX
10NTET 27NHEX
8NHEX

10−1
b(s)−bref

10−2
bref

10−3

10−4
100 101
s/sref

10−1
b(N )−bref

10−2
bref

10−3

10−4
102 103 104 105 106
Nelem

Fig. 3.28: Normalized error obtained with different 3D meshes. Here, bref = b(sref ) is
computed with 10NTET, and sref = 0.025 m.

tool to the YALES2 community to study case involving structural dynamics. In this PhD, it
has been created with the aim of reproducing FSI cases, and consequently, stationnary solving
is not used by the FSI solver. Surprisingly, several parts of the SMS were used to implement
a pseudo-solid method used for mesh movement, allowing considerable time saving benefit.
A FSI solver has then been developed after the SMS, and methods used are presented in the
next chapter.

85
 Chapter 4

FSI solver

This chapter describes the FSI solver, which couples the ALE solver and the SMS pre-
viously presented. It begins with the presentation of the original mesh movement method
developed to deform the fluid mesh according to the solid displacement. The FSI coupling
schemes are then introduced progressively and applied to a simple test case, and their stabil-
ity in function of a given parameter is studied. It finishes with the scheme used in the rest
of this work.

Contents
4.1 Mesh Movement Solver . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
4.1.1 Pseudo solid method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
4.1.2 Implementation in YALES2 . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
4.1.3 Adjustment of pseudo-material properties . . . . . . . . . . . . . . . . . . . 90
4.2 FSI coupling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
4.2.1 Data transfer via CWIPI . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
4.2.2 Weak coupling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
4.2.2.1 Case presentation . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
4.2.2.2 Parallel synchronous scheme . . . . . . . . . . . . . . . . . . . . . 97
4.2.2.3 Serial staggered scheme . . . . . . . . . . . . . . . . . . . . . . . . 100
4.2.3 Partitioned semi-implicit predictor-corrector coupling scheme . . . . . . . . 101
4.2.4 FSI coupling with DMA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105

4.1 Mesh Movement Solver

4.1.1 Pseudo solid method
As discussed in section 1.2.2, immersed boundaries methods do not allow the precise
description of geometries which is critical to predict accurately boundary layers in turbulent
cases. Therefore, body fitted meshes are prefered. Besides, the FSI coupling takes place at
the interface between fluid and structure, so that the fluid domain has to move according
to the solid displacement. This imposes to compute a grid movement for the entire fluid
domain while maintaining suitable cells quality. This is particularly challenging for cases

87
Chapter 4. FSI solver

with large deflection on unstructured grid, because the mesh has no preferential direction
contrary to structured mesh. This facilitates the distribution of the deformation among cells
with algebraic methods such as transfinite interpolation (TFI) [133, 134, 135, 136], tension
spring analogy [103] or torsion spring analogy [137]. These methods are either not appropriate
for unstructured meshes because it implies interpolation along mesh lines, or too expensive
for large scale problems, as explained in [138].
Furthermore, to be able to reproduce numerically fluid structure interaction cases with
turbulent flows, mesh movement method has to be suitable for a grid with a large number
of elements and different mesh size variation, especially for LES which requires fine grid
resolution. However, precision of spatial schemes suffer from mesh distorsion [139], so that
it is preferable to avoid to deform the mesh in regions where strong velocity gradients occur.
Thus, mesh movement requires an algorithm that can be fully parallelized and allowing to
control which regions of the domain will withstand the deformation. These reasons explain
the choice the pseudo-solid method [140, 141] for the present work. Besides, considering that
the SMS had been developed just before, the choice of and algorithm based on the FEM
seemed very natural. Note that alternatives will be discussed as perspective of this work in
section 7.2.
The principle of this method is to consider the fluid domain as a linear elastic solid. The
displacement constrains that the domain undergo then become usual boundaries conditions
of a solid with prescribed displacement, as seen in section 3.1.2.2. Inertial effects are not
advisable here so the problem is considered stationary. The mesh movement consists then in
solving a linear static problem. With a FEM formulation, it gives

Kd = f (4.1)

with K the stiffness matrix, d the nodes displacement and f the forces. With this formulation,
the equilibrium position of the pseudo-solid is considered to be the mesh configuration at
t = 0 = t0 . Here, there is no external forces such as body forces or pressure on boundaries,
however f depends on the displacement conditions. In this work, this method has been
improved to meet FSI simulations requirements. This is detailed in the next two sections.

4.1.2 Implementation in YALES2

As explained above, it is very useful to select regions of the domain that will move without
cell deformation, and the zone that will withstand the deformation. The pseudo-solid method
is very convenient for this purpose because it only requires to change the element pseudo-
Young modulus E, that will directly affect the cell flexibility. In fact, for an heterogeneous
solid, the part that will deform first is the less rigid one. Moreover, in most cases, the smallest
cells are regrouped close to the flexible body because it corresponds to the zone of interest,
and it is essential to maintain the quality of these cells. Therefore, the generic technique
developed here to compute the pseudo-Young modulus field of the pseudo-solid consists in
establishing a E profile as a function of the distance from the object R.

88
 4.1. Mesh Movement Solver

In this work, this distance R is computed with a geometrical method for unstructured
simplicial meshes [108], already used in two phases flow solvers of YALES2 to compute the
distance to the phases interface [107]. The pseudo-Poisson ratio of all elements is 0.2 and
their pseudo-Young modulus E is computed in function of R such as



 E(R < Rmin ) = 100Eint (4.2a)
(R − Rmin )(Eint − Emin )



E(Rmin ≤ R ≤ Rmax ) = Eint − (4.2b)

Rmax − Rmin




 E(R > Rmax ) = Emin (4.2c)
Eint = Emin yr

(4.2d)

where the pseudo-Young modulus ratio yr , Rmin and Rmax are user’s parameters. On the
other hand, the value of Emin has no influence since it is a linear elastic problem. This
method ensures that for R < Rmin , deformation is minimal and cells move mainly according
to the imposed displacement, thanks to the high jump of E at Rmin . On the contrary, the
mesh deformation takes place in the transition zone between Rmin and Rmax . The size of
the different regions can be easily adapted for the different cases with these two parameters.
However the deformation is not uniform in the transition zone, because for yr = 1, only cells
at R = Rmin are deformed, and for high value of yr , most of the deformation will be taken
by cells at R ≈ Rmax . The suitable value for yr then depends on the size of the transition
zone, the displacement amplitude and the cells size distribution. An example of both R and
E fields can be visualized in Fig. 4.1, where the geometry used for the example corresponds
to the case presented in Chapter 5.

Fig. 4.1: Distance field from the object R (left) and pseudo-Young modulus field E (right).

This method has been proven to be efficient but could be improved; the cells quality of the
elements with the lowest E frequently dropped first, but the deformation did not sufficiently
affect cells closer to the object so that the computations stopped and failed even though cells
at Rmin ≤ R ≤ Rmax were still nearly non deformed.

89
Chapter 4. FSI solver

4.1.3 Adjustment of pseudo-material properties

In order to improve the previous method, the idea is to make the cells pseudo-Young
modulus E depend on their quality. The latter can be determined quantitatively by com-
puting the cells skewness S defined by Eq.(2.14). This criterion is then used to adjust the
pseudo-solid elements flexibility during the simulation such as
  
Eint − E(t0 ) S(tn ) − S(t0 )
E(tn+1 , S(tn )) = + E(t0 ). (4.3)
Smax − S(t0 )
where Smax is chosen at 0.90. Note that the pseudo-Young modulus at the beginning of the
computation E(t0 ) is established according to the Eq.(4.2). Equation(4.3) is applied only for
cells where R > Rmin because the properties of cells close to the object has to be homogeneous
as most as possible to guarantee that these elements remain undeformed.
However, this single modification is not enough to considerably improve the method. In
fact, the solving of Eq.(4.1) for a strong displacement will result in a mesh with a deformed
cells zone. That will induce an increase of the pseudo-Young modulus of these cells because of
Eq.(4.3), so that they will become more rigid than their neighbours. The solving at the next
time step will then lead to a mesh where that same zone has not been deformed because of the
previously computed high value of E, but their neighbours on the contrary will withstand
all the deformation. The issue here is that the solved problem is still defined according
to the initial configuration, and therefore the cell skewness S cannot evolve progressively.
To make this technique efficient, the reference configuration used to compute the stiffness
matrix K has to be updated at each time step with the new node position, which means to
update continuously the equilibrium position of the pseudo-solid. The computed displacement
becomes then an increment ddn+1 used to move nodes between their position at tn and the
new one at tn+1 . Besides, to prevent high gradient of E, the pseudo-Young modulus computed
with Eq.(4.3) is filtered once with a scatter-gather algorithm, given that this data is stored
at the element. Tests showed that this strategy improves significantly the efficiency of the
algorithm.
An example of application of the full method is given in Fig. 4.2. It shows skewness
fields of a test case where a domain boundary, the rod behind the cylinder, moves accord-
ing a prescribed sinusoidal displacement at different instants. For this test, a mesh with
inhomogeneous cell size has been used; close to the cylinder and the rod, a zone with small
elements can be identified. The mesh movement strategy has to ensure that these cells will
not be too much deformed as explained in previous section. With the present method, the
good behaviour is obtained by choosing a correct value for Rmin , and the different snapshots
highlight that the skewness in this region remains intact. Moreover, it can also be seen that
the two critical zones in terms of skewness progressively appear; one at R = Rmin and the
other at R = Rmax . That can be explained by the stronger pseudo-Young modulus gradients
of these regions. Nevertheless, the two last snapshots prove that the rigidification of these
cells leads to a smooth cell deformation, resulting in a satisfying skewness distribution.
A comparison of the method of section 4.1.2 and the method of section 4.1.3 is given
in Fig. 4.3. It shows that without the new formulation, only elements with low value of E

90
 4.1. Mesh Movement Solver

Fig. 4.2: Skewness field of a test case with imposed displacement for the moving boundary
at different instants.

tend to deform, while the other cells in the transition zone remain nearly unchanged. With
the adjustment of pseudo-material elements properties, it is clear that the skewness field is
much more homogeneous, thanks to the rigidification of close wall cells. This last technique
appeared consequently more robust to handle strong displacement and maintain a lower
maximum skewness. An example of pseudo-Young modulus field after a large displacement is
given in Fig. 4.4. It emphasizes that regions forced to stretch see their cell skewness decrease,
and on the other hand compressed regions undergo growing skewness. This causes a pseudo-
Young modulus profile quite different from the one at t = 0 given Fig. 4.1. This efficient
mesh deformation method presents then the advantage to be robust, easy to implement and
without parallelization issues. Furthermore, it allows the user to control which mesh region
can be deformed or not.
Nevertheless, for certain cases with very large deformation, moving or rotating object,
DMA is used (see section 2.2.1.4). To combine efficiently this feature with the pseudo-
material method, the stiffness matrix K has to be updated at each re-meshing step. In
fact, a new mesh with a low maximum skewness has to be considered as a new equilibrium
position. Besides, the previous E field resulting from adjusting the pseudo-Young modulus as
a function of skewness becomes irrelevant because the skewness field is completely changed.
The distance from the object R is then recomputed and a new pseudo-Young modulus field
is established with Eq.(4.2). The entire method for mesh movement developed in this work
is summarized in Fig. 4.5.

91
Chapter 4. FSI solver

Fig. 4.3: Comparison of the pseudo-solid method of section 4.1.2 (top) and the method of
section 4.1.3 (bottom).

Fig. 4.4: Pseudo-Young modulus field after large deflection of the rod while adjusting the
pseudo-material properties with Eq.(4.3).

Finally, it must also be noticed that the chosen time step can have a strong influence on the

92
 4.1. Mesh Movement Solver

Fig. 4.5: Scheme of the mesh movement method developed in this work.

reversebility of the algorithm. Figure 4.6 shows the evolution of the maximum skewness of the
domain in function of time for the previously introduced case with a τ periodic prescribed
motion of the rod. It can be seen that for ∆t/τ = 2.5 × 10−2 the maximum skewness is
increasing for each cycle of prescribed displacement, even if the geometry returns to the
initial configuration every τ /2 seconds. This is caused by the incremental formulation used
with the adjustment of properties. On the other hand, with a time step ten times smaller,
this effect vanishes nearly entirely. Note that DMA is not used here to highlight this effect.
In practice, considering the time step necessary for FSI computations, this effect remains
acceptable and is not restrictive because if the maximum skewness finally reach too high
value (kSk∞ < Slim with Slim = 0.99), the DMA can still be used.

93
Chapter 4. FSI solver

0.1
Displacement (m)

0.0 Prescribed tip displacement

−0.1
0 1 2 3 4 5
t/τ
1.0
Skewness

0.8

0 1 2 3 4 5
t/τ
Maximum skewness for ∆t/τ = 2.5 × 10−2 Maximum skewness for ∆t/τ = 2.5 × 10−3

Fig. 4.6: Comparison of the maximum skewness obtained for different timestep ∆t. DMA is
not used here.

4.2 FSI coupling

4.2.1 Data transfer via CWIPI
FSI simulations consist in coupling fluid (ALE) and solid (SMS) solvers. The coupling
conditions occur at the interface Γ between fluid and solid as,

dn+1 − dn

 un+1 = wn+1 = (4.4a)


ˆ  ∆t 
∂un+1
 fn+1 = µ + Pn+1 n dΓ (4.4b)


Γ ∂n

with µ the dynamic viscosity and n the normal direction to Γ pointing on the solid to fluid
direction. In fact, solid imposes its velocity to fluid while fluid applies a force on solid through
viscous shear and pressure. As a partitioned approach is chosen here (see section 1.2), the
interface Γ displacement d has to be sent from the SMS to the ALE solver, and the fluid
forces applied on Γ f has to be sent from the ALE solver to the SMS.
The coupling is then performed by data exchange between solvers at the interface Γ
ensured by the CWIPI (Coupling With Interpolation Parallel Interface) library [142, 143].
This library allows to compute data interpolation based on boundary nodes coordinates,

94
 4.2. FSI coupling

which is mandatory when meshes are not coincident. Nonetheless, the library has not been
designed to support quadratic elements. As a FSI simulation involves surfacic coupling, faces
of the quadratic elements used by the SMS have to be cut into linear elements, as shown in
Fig. 4.7.

Fig. 4.7: Adaptation of quadratic elements faces to make the SMS work with CWIPI.

CWIPI can work with two meshes composed of different element type, with different
cell sizes, as illustrated in Fig 4.8. It also manages the parallel communications allowing to
benefit of the massively parallel performance of each solver. As the processors distribution
on fluid side is completely changed after a re-meshing step, the coupling with CWIPI has to
be entirely reestablished after each DMA.
Now that the fluid solver (Chapter 2), the solid solver (Chapter 3) and the mesh movement
solver have been presented, next part will present different possible coupling scheme, using a
simple test case.

95
Chapter 4. FSI solver

Fig. 4.8: Data transfer via CWIPI from fluid mesh composed with tetrahedrons of size ∆xf
(left) to solid mesh composed with 27-nodes hexahedrons of size ∆xs = 11∆xf (right).

96
 4.2. FSI coupling

4.2.2 Weak coupling

4.2.2.1 Case presentation

This part aims at explaining why a strong coupling scheme is necessary for certain cases.
The coupling scheme used in the rest of this work will then be introduced step by step.
For that, using a simple test case, the density ratio ρs /ρf will be varied to highlight the
instabilities induced by the added mass effect. Note that other parameters, like geometry,
have also an influence on this effect [58, 59], but only ρs /ρf will be discussed here. Thus,
increasingly complex schemes will be tested to study their robustness against more and more
unstable cases. It will finally lead to the actual scheme used to reproduce cases of next
chapters. It should also be precised that more sophisticated schemes exist in the literature,
and that will be discussed as perspective of this work in section 7.2.
A simple 3D case of a flexible beam in a channel has been devised to study the stability
of the different schemes with low cost tests. The case is illustrated in Fig. 4.9. The beam
has a lenght l = 1 m, a square section of side a = 0.1 m and is reproduced with 10 27-
nodes hexahedrons. The material properties are E = 3.6 GPa, νs = 0.4, and ρs will be case
dependent. The fluid domain dimensions are 5m × 4m × 3m and the fluid mesh is composed
by 38 × 103 tetrahedrons. The fluid properties are νf = 10−6 m2 .s−1 and ρf = 103 kg.m−3 .
The inflow velocity is U∞ = 1 m.s−1 . In order to increase the structure deflection, a sinusoidal
displacement of period T = 1 s and amplitude A = 0.1 m is imposed at the basis of the beam.
For the mesh movement, the pseudo-solid method presented in the last section is used and
slip boundary conditions are used for channel walls.
The beam will then bend under the flow effect. Fluid forces applied on the structure and
beam tip displacement will be used for comparison.

4.2.2.2 Parallel synchronous scheme

The first approach proposed is very simple; data for coupling are exchanged at the be-
ginning of the time step and both fluid and solid solving are computed simultaneously. The
scheme is given in Fig. 4.10. The n + 1 time iteration begins with computation of the time
step value ∆tf given the stability condition of the Runge–Kutta method used by fluid solver.
Note that the SMS uses the generalized-α (see section 3.1.3.3) temporal scheme which is un-
conditionally stable. Thus, ∆tf is sent from the fluid to the solid and ∆ts = ∆tf is imposed
for the SMS. In the present case, it results in a time step of about 15 ms. This method
ensures optimal computational cost, but not the consistency between fluid and solid solution
at the end ot the time step. Figure 4.11 compares the results obtained during the first ten
seconds with ρs /ρf = [3.0, 2.0, 1.0]. It can clearly be seen that instabilities on forces (note
that logarithmic scale is used) and dx are growing when ρs /ρf decreases. For ρs /ρf = 1,
the forces diverges so that the computation failed. It must also be noted that important
instabilities on the drag force reappears for ρs /ρf = 2 even after several seconds with small
variations.

97
Chapter 4. FSI solver

Fig. 4.9: Test case configuration (top) and cross-section (bottom).

98
 4.2. FSI coupling

Fig. 4.10: Parallel synchronous scheme. Fluid and solid solving are performed simultaneously.

ρs/ρf = 1.0 ρs/ρf = 3.0

ρs/ρf = 2.0
1063
Lift (N )

10
0
−1036
−10 0 2 4 6 8 10
1063
Drag (N )

10
0
−1036
−10
0 2 4 6 8 10
dx (m)

0.25

0.00
0 2 4 6 8 10
0.2
dy (m)

0.0

0 2 4 6 8 10
time (s)

Fig. 4.11: Result obtained with the parallel synchronous scheme for different density ratios
ρs /ρf .

99
Chapter 4. FSI solver

This low cost scheme can be useful for simple case where a global behaviour is studied,
but lacks of accuracy and does not allow to reproduce cas with important added mass effect.

4.2.2.3 Serial staggered scheme

This second scheme is presented in Fig. 4.12. It is sequential and begins with an estimation
of structure displacement with
d0n+1 = dn + ∆tf vn (4.5)
where v is the solid velocity. This extrapolated displacement is then used by the ALE solver
to compute fluid forces. This latter are sent to the SMS and taken into account to compute
dn+1 . This method ensures a better consistency between solutions than previous scheme, but

Fig. 4.12: Serial staggered scheme.

it is not perfectly accurate as d0n+1 6= dn+1 . Figure 4.13 compares the results obtained during
the first two seconds with ρs /ρf = [1.0, 0.75, 0.5]. With this scheme, case with ρs /ρf = 1
can be reproduced. Nonetheless, for lower density ratios, instabilities appear. They finally
vanish for ρs /ρf = 0.75 but they make the computation fail for ρs /ρf = 0.5. This scheme
is then slower than the previous one, but more precise; however, it does not allow case with
important mass effect. In this test, this problem appears for very low ρs /ρf , but for other
conditions, higher ratios can still result in unstable case, like in section 5.4.2 where ρs /ρf = 10
but the case cannot be reproduced with a weak coupling.
In order to enforce the equilibrium of the traction and displacement on the fluid-structure
interface and ensure consistency between fluid and solid solutions at each time step, a strong
coupling is required. This implies an additional iterative procedure during time advancement.
The chosen strong scheme is explained in the next section.

100
 4.2. FSI coupling

ρs/ρf = 0.5 ρs/ρf = 1

ρs/ρf = 0.75
1063
Lift (N )

10
0
−1036
−100.00 0.25 0.50 0.75 1.00 1.25 1.50 1.75 2.00
1063
Drag (N )

10
0
−1036
−100.00 0.25 0.50 0.75 1.00 1.25 1.50 1.75 2.00
dx (m)

0.25

0.00
0.00 0.25 0.50 0.75 1.00 1.25 1.50 1.75 2.00
0.2
dy (m)

0.0
−0.2
0.00 0.25 0.50 0.75 1.00 1.25 1.50 1.75 2.00
time (s)

Fig. 4.13: Result obtained with the serial staggered scheme for different density ratios ρs /ρf .

4.2.3 Partitioned semi-implicit predictor-corrector coupling

scheme
The FSI coupling scheme applied is given Fig. 4.14 and is similar to the one proposed by
Breuer et al. [78]. This coupling scheme is the one used in the rest of this work, so it will be
explained in details. It can be split into different steps:
1. The n + 1 time iteration begins with computation of the time step value ∆tf given the
stability condition of the Runge–Kutta method used by fluid solver. Note that the SMS
uses the generalized-α (see section 3.1.3.3) temporal scheme which is unconditionally
stable. Thus, ∆tf is sent from the fluid to the solid and ∆ts = ∆tf is imposed for the
SMS.
2. New displacement of the interface Γ is estimated thanks to the following extrapolation

d0n+1 = dn + α0 ∆tf vn + α1 ∆tf (vn−1 − vn ) (4.6)

where v is the solid velocity, and αi are interpolation coefficients. For α0 = 1 and
α1 = 0, Eq.(4.6) corresponds to a linear extrapolation while for α0 = 1 and α1 = 1/2 it
corresponds to quadratic extrapolation. This last option leads to a faster convergence
so it is used in this work.

101
Chapter 4. FSI solver

Fig. 4.14: FSI coupling scheme used for the present work.

3. Mesh movement is solved thanks to the method explained in section 4.1 with the dis-
placement condition d0n+1 . That results in a nodes velocity field wn+1
0
used to perform
the ”velocity prediction” step, as detailed in section 2.2.1.1. The intermediate velocity
u∗n+1 obtained is stored apart.

4. A new pressure field is then determined with Eq.(2.10) and used to perform the ”velocity
correction step” via Eq.(2.9d). With the resulting pressure and velocity fields, the fluid
k
force applied on the solid fn+1 can be computed. Here the upper index k corresponds
to the number of subiterations in the FSI loop. Note that absence of this index refers
to converged values.

k
5. The fn+1 field is sent from the fluid to the solid by using CWIPI library.

6. The displacement dkn+1 in the solid is computed using the methodology detailed in
section 3.1.

7. At this stage, both fluid and solid have been solved at least once. The dynamic equi-
librium is then checked to verify the consistency of the two solutions. This checking is

102
 4.2. FSI coupling

performed with
dkn+1 − d̃k−1
n+1
∞
≤ εF SI (4.7)
dkn+1 − dn ∞

where εF SI is a chosen convergence criterion, and d̃0n+1 = d0n+1 . If the FSI solution is
converged, the solution at the considered time is finally determined and next time step
can be started.

8. Otherwise, the computed displacement has to be underrelaxed on Γ with

d̃kn+1 = ωdkn+1 + (1 − ω)dk−1

n+1 (4.8)

where ω is a constant underrelaxation factor defined by the user. This value can also be
computed at each subiteration by different means [67] but tests for our cases show that
different methods did not allow to guarantee stability. Then, an accurate value of ω has
to be found for the different case; a too large value will make the computation diverge,
or will cause instability on fluid forces. On the other hand, the lower ω is, the greater
will be the total number of subiterations NF SI required before reaching convergence.

9. The underrelaxed displacement d̃kn+1 is then sent to the fluid via CWIPI.
k
10. A new field of wn+1 is computed thanks to the mesh movement solver (MMS). It has to
be precised that before the solving, the domain is put back to its position at tn . A new
increment of displacement has to be computed but the adjustment of pseudo-material
parameters with Eq.(4.3) has to be done only once the convergence is reached.

With the previously stored velocity prediction u∗n+1 and the updated wn+1 k
, new pressure
and velocity fields can be determined so that step 4 can be repeated and updated fluid
k+1
forces fn+1 can be computed, closing so the FSI loop. This method thus ensures stability
with a partitioned coupling and can be easily implemented, but presents the main drawback
that the computational cost will depend on the number of FSI subiterations NF SI necessary
to reach convergence. This number depends on a lot of factors and is initially impossible
to determined. The skipping of the velocity prediction once in the FSI loop was initially
proposed by Fernandez et al. [144]. It allows to reduce computational cost and it has been
shown that it does not affect the final result. Also, the underrelaxation and the convergence
k
checking can be performed on fn+1 instead of dkn+1 .
The case of section 4.2.2.1 is now reproduced with this scheme for ρs /ρf = 0.5. Numerical
parameters are ω = 0.5 and εF SI = 10−3 . Figure 4.15 shows the results. It shows that the
computed solution is very accurate and whithout instabilities. Only the first time iterations
presents a spike of drag, which is due to the non physical initialization of the flow. It can also
be observed on the curve of NF SI that solid state corresponding to the maximums and the
minimums of dy are unstable; at these instants, structure has a small inertia so that it requires
more subiterations to determine the two consistent solutions. This phenomenon will also be

103
Chapter 4. FSI solver

1063
Lift (N )

10
0
−1036
−10 0 2 4 6 8 10
1063
Displacement (m)Drag (N )

10
0
−1036
−10 0 2 4 6 8 10
0.25 dx
0.00 dy

0 2 4 6 8 10
10
NF SI

FSI subiterations
5
0 2 4 6 8 10
time (s)
n+1 ||∞

||dkn+1−dn||∞

100
Solid convergence
||∞ ||dkn+1−dk−1

10−2

0 25 50 75 100 125 150 175

Number of subiteration
100
n||∞
k−1
−fn+1
k −f

10−2
||fn+1
||fn+1

Fluid convergence
k

0 25 50 75 100 125 150 175

Number of subiteration

Number of Subiterations
10
NF SI

0.0 2.5 5.0 7.5 10.0 12.5 15.0 17.5

Number of time iteration

Fig. 4.15: Result obtained with the strong coupling scheme for ρs /ρf = 0.5 (top) and evolu-
tion of solid and force convergence with subiterations for the first 20 time iterations (bottom).

104
 4.2. FSI coupling

observed on cases of the next chapters. Second graph also highlights that displacement and
force have the same convergence behaviour.
This coupling scheme has finally been adopted for the present work. On each case,
tests have to be done to find suitable values of ω and εF SI to guarantee stability but also to
minimize the computational cost. Also, in several cases, it has been noticed that the flow had
to be established before starting the FSI coupling. As a matter of fact, whithout substantial
resulting fluid forces on the structure, it seems that the scheme does not allow to find a
stable solution and the successive computed forces within the FSI loop tend to diverge. This
specific point has been a large source of difficulties in the PhD and took a lot of time to be
identified. An important perspective of this work could be to develop a robust initialization
strategy.

4.2.4 FSI coupling with DMA

It was important to be able to use DMA with the FSI solver to reproduce case with large
deflection or flexible moving object. However, as mentioned in section 2.2.1.4, DMA can
induce pressure spike in some configurations. This strong discontinuity on fluid force can
disturb the FSI coupling and make the computation diverge. Attempts had been made to
correct the pressure spike and reconstruct another pressure field after DMA but it was not
enough to stabilize the coupling.
The solution finally adopted consists in an extrapolation of the fluid forces at the time
step after the re-meshing. It is computed as

fn+2 = 2 fn+1 − fn (4.9)

k
and is sent to the solid instead of the normally computed fn+2 at each subiteration. The
underrelaxation factor is then temporarily set to 1 so that d2n+1 − d̃1n+1 = 0 at the second
subiteration and NF SI is systematically equal to 2 after DMA. This solution was very simple
to implement but revealed quite efficient. Figure 4.16 shows the results obtained for the
previous case where DMA was forced every 100 time iterations. The DMA step can be
identified on the skewness curves by the sudden drops. It can also be noticed that at next
time iteration, NF SI = 2. Other cases with DMA will be presented in next Chapters.

The entire numerical methodology used in the FSI solver has now been presented. The
three schemes studied in this Chapter are still available in the FSI solver of YALES2. These
developments occupied approximatively the first half of the PhD. It was then essential to
correctly validate this new solver. Next Chapter will present the validation against a 2D
laminar reference test case.

105
Chapter 4. FSI solver

Fig. 4.16: Result obtained with the strong coupling scheme for ρs /ρf = 0.5 with DMA.

106
 Chapter 5

Validation against 2D numerical

laminar cases

In this chapter, a complete validation of the previously introduced FSI solver is given,
including mesh convergence studies. For that, a 2D numerical benchmark is reproduced. It
involves the interaction between an elastic rod and a laminar incompressible flow. The chapter
starts with validations of fluid and solid solvers independently, to finally presents two cases
of FSI, with and without Dynamic Mesh Adaptation.

Contents
5.1 Case presentation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107
5.2 CFD tests . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108
5.3 CSM tests . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110
5.4 FSI tests . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112
5.4.1 FSI3 case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112
5.4.2 FSI2 case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116

5.1 Case presentation

In order to assess the present methodology, a 2D laminar FSI test case is first considered:
the benchmark proposed by Turek & Hron [145]. This numerical case consists in the interac-
tion between an elastic rod and a laminar incompressible flow. The benchmark proposes to
validate each solver independently in preliminary tests. Therefore, three different test cases
including mesh convergence studies are presented in this part.
This case is inspired by the older benchmark of an incompressible flow around a cylin-
der [146] except that a flexible rod is attached to the back side of the cylinder. The geometry
and dimensions are given in Fig. 5.1 and Tab. 5.1. By measuring from the left bottom corner
of the channel, the cylinder center position is then C = (0.2, 0.2) while the rod tip is situated
at A = (0.6, 0.2). It should be noticed that the setup is intentionally non-symmetric (with
Hc 6= H/2) to prevent any influence of the computation precision. For the three folowing
tests, reference values have been taken from [145]. They are obtained with a fully implicit
monolithic ALE-FEM method with a fully coupled multigrid solver [60] and are almost grid
independent.

107
Chapter 5. Validation against 2D numerical laminar cases

Fig. 5.1: Sketch of the studied numerical test case.

Geometry parameters Value [m]

Cylinder diameter D = 0.1
Cylinder center x-position Lc = 0.2
Cylinder center y-position Hc = 0.2
Channel lenght L = 2.5
Channel height H = 0.41
Deformable structure length l = 0.35
Deformable structure thickness h = 0.02

Table 5.1: Dimensions of the sketch in Fig. 5.1.

5.2 CFD tests

To first validate the fluid solver alone, the benchmark proposes three cases with different
physical parameters: the case named CFD3 is chosen here. The solid is considered perfectly
rigid so that the test case focuses on the laminar flow description around the cylinder and
the attached rod. Quantities used for comparison will then be the fluid forces applied on the
whole submerged body, computed according Eq. (1.1b), for a fully developed flow and for
one full period of the oscillation. In fact, for this case, a non stationary regime is reached
where pressure distribution fluctuates.
The flow is considered as incompressible, with a density ρf = 1000 kg.m−3 and a kinematic
viscosity of νf = 0.001 m2 .s−1 . A parabolic velocity profile u(x = 0, y) is prescribed at the
inlet as
y(H − y)
u(x = 0, y) = 1.5u∞ . (5.1)
(H/2)2

108
 5.2. CFD tests

In CFD3 case, u∞ = 2 m.s−1 leading to a Reynolds number Re = 200 which leads to vortex
shedding behind the cylinder, as illustrated in Fig. 5.4. No-slip boundary conditions are
applied at channel walls and at the body.
Simulations are performed on three different meshes, M1f , M2f and M3f composed of
triangles. They are characterized by two metric values, ∆x1 and ∆x2 , which correspond to
cell size close to the cylinder and the rod, and to the cell size in the rest of the domain,
respectively. These values are given in Tab. 5.2.

M1f M2f M3f

∆x1 [m] 1 × 10−3 8 × 10−4 6 × 10−4
∆x2 [m] 1.25 × 10−2 1 × 10−2 0.75 × 10−2
Nelem (×103 ) 73 98 173

Table 5.2: Characteristics of the three meshes used for CFD3.

The timestep ∆t used here is computed following the Courant-Friedrichs-Levy (CFL)

convective time step constraint by keeping CFL number CF L = u∆t/∆x smaller than 0.8.
This leads to a time step value around 139 µs, 113 µs and 84 µs for meshes M1f , M2f and
M3f , respectively.
As precised above, fluid forces oscillations are compared here. Amplitudes, mean value
and frequencies have been computed for each mesh and are presented in Fig. 5.2 and in
Tab. 5.3. Despite that the time-steps are about 50 times smaller than in the reference

Mean Drag [N] Amp. Drag [N] Mean Lift [N] Amp. Lift [N] f [Hz]
Ref [145] 439.45 5.6183 -11.893 437.81 4.3956
433.69 6.4297 -7.0625 485.57 4.4326
M1f
(1.31%) (14.44%) (40.62%) (10.91%) (0.84%)
435.35 6.3794 -7.4314 478.94 4.4347
M2f
(0.93%) (13.55%) (37.51%) (9.39%) (0.89%)
435.01 6.0761 -10.505 462.78 4.4373
M3f
(1.00%) (8.15%) (11.67%) (5.70%) (0.95%)

Table 5.3: Results of the CFD3 benchmark and comparison with the reference study [145].

case, amplitudes of drag and lift and mean lift converge towards the reference data when the
grid is refined. Note that the mean drag force and the frequency are less sensitive to mesh
refinement with an error about of 1% even for the coarser mesh. Therefore, the results allow
to consider that the CFD solver is validated.

109
Chapter 5. Validation against 2D numerical laminar cases

M3f Ref
500
Lift (N )

−500
9.0 9.1 9.2 9.3 9.4 9.5 9.6
time (s)
Drag (N )

440

430
9.0 9.1 9.2 9.3 9.4 9.5 9.6
time (s)

Fig. 5.2: Fluid forces integrated on the whole submerged body (cylinder+rod) computed
with M3f and results of the reference study [145]. Time offset is due to different computation
initializations.

5.3 CSM tests

The solid solver needs now to be validated with pure structural test. This test consists in
computing the deformation of the flexible rod in a gravitational field g = (0, 2) m.s−2 without
taking into account the fluid. In this study, the CSM3 test is chosen because it is the only
time dependent case. It starts from the undeformed configuration and as there is no damping,
the structure immediately oscillates periodically. The material is characterized by a Poisson’s
ratio of νs = 0.4, a Young modulus of E = 1.4 MPa and a density of ρs = 1000 kg.m−3 . As
significant deformations are expected, the Saint-Venant-Kirchhoff material model is used.
Three meshes have been made for this test, corresponding only to the deformable part
of the body because the cylinder is always considered perfectly rigid. These meshes are only
composed by 9-nodes quadrilaterals of size ∆x given in Tab. 5.4.
As the structure is clamped to the backside of the cylinder, a non displacement boundary

110
 5.3. CSM tests

M1s M2s M3s

∆x [m] 4 × 10−2 2 × 10−2 1 × 10−2
Nelem 9 18 70

Table 5.4: Characteristics of the three meshes used for CSM3.

condition is applied to the corresponding nodes. The chosen time step is ∆t = 0.005 s (same
than in [145]) and is applied with the generalized-α method with the spectral radius ρ∞ = 0.8
(Eq.(3.85)).
For comparison, displacement of previously defined point A is tracked in time. Once
again, mean values, amplitudes and frequencies are computed for displacements dx and dy
along x-axis and y-axis, respectively. All results are gathered in Tab. 5.5 and in Fig. 5.3.

dx M3s dx Ref
dy M3s dy Ref

0.000

−0.025
Displacement (m)

−0.050

−0.075

−0.100

−0.125
8.00 8.25 8.50 8.75 9.00 9.25 9.50 9.75 10.00
time (s)

Fig. 5.3: Displacement of point A along x-axis (top) and y-axis (bottom) obtained with M3s
and results of the reference study [145].

Differences with reference values are here minor, even though grid used in the reference
study is much finer (5120 elements). Efficiency of 9-nodes quadrilaterals seems here confirmed
since errors are not above 2% for all quantities, except for the amplitude of dy obtained with
M1s , which converges with grid refinement anyway. This close agreement with reference data
allows then to successfully validate the CSM solver.

111
Chapter 5. Validation against 2D numerical laminar cases

Mean dx [mm] Amp. dx [mm] Mean dy [mm] Amp. dy [mm] f [Hz]

Ref [145] -14.305 14.305 -63.607 65.160 1.0995
-14.057 14.057 -63.367 63.367 1.1074
M1s
(1.73%) (1.73%) (0.38%) (2.66%) (0.72%)
-14.398 14.394 -64.113 64.436 1.1053
M2s
(0.65%) (0.63%) (0.80%) (1.10%) (0.53%)
-14.257 14.452 -64.163 64.872 1.0972
M3s
(1.06%) (1.03%) (0.87%) (0.44%) (0.21%)

Table 5.5: Results of the CSM3 benchmark and comparison with the reference study [145].

5.4 FSI tests

5.4.1 FSI3 case
Finally, to validate the coupling between fluid and solid solvers, the FSI3 benchmark has
been reproduced. The fluid forces are now applied to the flexible structure so that it deforms
and starts interacting with the flow. Fluid properties are the same as in CFD3 but for the
solid, the Young modulus is chosen as E = 5.6 MPa. This case is particularly challenging
because ρf = ρs = 1000 k.m−3 , which maximizes the importance of the added-mass effect.
Besides, the channel is narrow compared to the expected structure deflections; that imposes
a very efficient mesh movement algorithm to keep a low maximum skewness. The method
presented in section 4.1 is used here with Rmin = 0.02 m, Rmax = 0.2 m and yr = 100 m.
These values allow the most refined zone close to the body to remain intact while coarser
regions close to the channel walls will withstand the deformation. That can be seen in Fig. 5.4
which shows the resulting grid M1f when the rod is at the maximum deflection. The method
ensures a good grid quality without needs of DMA step. As regards the FSI coupling, the
convergence criterion defined in Eq. (4.7) is set at εF SI = 1 × 10−5 and the underrelaxation
factor is ω = 0.1. Because of the density ratio ρs /ρf = 1, higher values of ω do not allow
to reach convergence within the FSI loop. Nonetheless, it results in a relatively low mean
number of subiterations NF SI of 15.04 while this number averaged 55.21 in [78] for the same
value of ω. However, Breuer et al. [78] managed to reach convergence with ω = 0.5 and
NF SI = 9.38 which was not possible in the present study.
In the original case, the authors propose to progressively establish the flow starting with
a zero velocity field in all the domain. Note that this cannot be done with the present
algorithm because the fluid forces were so low that the importance of the added mass effect
was relatively too important. Such an unstable configuration surely requires a more stable
monolithic approach. Therefore, the coupling has here been started with a fully developed
flow, as computed in CFD3 case.
Figure 5.4 also depicts an instantaneous velocity field. The vortex shedding at the cylinder
is clearly visible and induces oscillating forces applied in the entire body as it has been seen in

112
 5.4. FSI tests

the CFD3 case of section 5.2. These forces cause the flexible structure displacement resulting
in periodic displacement of the rod behind the cylinder.

Fig. 5.4: Velocity field and deformed mesh of FSI3 benchmark reproduced with M 1f & M1s .
Note that the domain has been cropped on the right side to focus on the deformable part.

For comparison, sinusoidal displacement of previously defined point A is measured and

presented in Fig. 5.5 and in Tab. 5.6. Even if numerical methodologies between the present
study and the reference one are different, the results converge toward the reference values.
The only disparity occurs for the mean dy obtained with mesh M2f &M2s . Considering that it
concerns a small value and that error obtained with mesh M3f &M3s drops below the one with
M1f &M1s , it can be considered as a minor deviation. Consequently the overall agreement
with reference data enables to validate the present FSI solver.
More data are given for the computation with M3f &M3s in Fig. 5.7. This has been
performed with 16 CPU for the ALE solver and 1 for the SMS; the mesh movement solving
took about 61% of the computational time, the fluid solving 26% and solid solving 13%.
With the three first graphes, it can be deduced that the computed solution is very stable. As
regards the mesh movement, the nearly constant mean skewness in the domain confirms that
only few cells are problematic. As for the maximum skewness, it slowly increases because of
the non reversibility of the algorithm already highlighted in Fig. 4.6. However, the maximum
skewness finally reaches values where the variation are stable and which allow the computation
to continue despite high skewness. The cycles observed of the number of subiterations NF SI
were already observed in Fig. 4.15. It corresponds to the same oscillations than dy .

113
Chapter 5. Validation against 2D numerical laminar cases

Fig. 5.5: Displacement of point A along x-axis (top) and y-axis (bottom) obtained with
M3f &M3s and results of the reference study [145].

Fig. 5.6: Results of the FSI3 benchmark and comparison with the reference study [145].

114
 5.4. FSI tests

250
Lift (N )

−250
0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5
500
Drag (N )

450

0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5

Displacement (m)

0.025 dx
0.000 dy
−0.025
0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5
1.0
Skewness

Mean skewness
0.5
Maximum skewness

0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5

100
NF SI

0
0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5
time (s)

Fig. 5.7: Results obtained for FSI3 with M 3f & M3s .

115
Chapter 5. Validation against 2D numerical laminar cases

5.4.2 FSI2 case

Another proposed benchmark FSI2 [145] has also been reproduced and statisfying results
have been obtained. In this test case, ρs /ρf = 10 which makes the coupling less unstable,
but it still requires the strong coupling scheme. This allows to apply an underrelaxation of
ω = 0.3 and convergence is reached with a number of FSI subiterations that averages 4.5.
Contrary to the FSI3 case, here the larger deflections impose the use of the DMA, which
took less than 0.1% of the computational time. Animation of this case can be visualized at
https://www.youtube.com/watch?v=qLWOg_EJtTM. With M 2f & M3s , the amplitudes were
predicted with less than 2.3% of differences and frequencies with less than 4.5% in comparison
with [145].

Fig. 5.8: Velocity field and deformed mesh of FSI2 benchmark reproduced with M 2f & M3s .
Note that the domain has been cropped on the right side to focus on the deformable part.

Results are given in Fig. 5.8 and in Fig. 5.9. Here the effect of DMA on pressure men-
tionned in section 4.2.4 can clearly be seen on the two first graphes. At the time iteration
after DMA, the Eq.( 4.9) is then used, allowing the pursuit of the computation. Neverthe-
less, this destabilizes the solution as it can be deduced from the graph of NF SI . The effect

116
 5.4. FSI tests

1000
Lift (N )

3 4 5 6 7 8 9 10 11
1000
Drag (N )

500

3 4 5 6 7 8 9 10 11
Displacement (m)

0.05 dx
0.00 dy
−0.05
3 4 5 6 7 8 9 10 11
Skewness

0.75
Mean skewness
0.50
Maximum skewness
0.25
3 4 5 6 7 8 9 10 11
15
NF SI

10
5

3 4 5 6 7 8 9 10 11
time (s)

Fig. 5.9: Results obtained for FSI2 with M 2f & M3s .

117
Chapter 5. Validation against 2D numerical laminar cases

eventually vanishes but still increases the computational cost.

The entire FSI solver has thus been validated with success on a reference benchmark. It
can be concluded that the strong coupling scheme developed in this work allows to reproduce
accurately cases with important added mass effect. Also, the FSI2 case confirms that the
DMA can be used for FSI coupling, allowing to reproduce case with large deflection on
unstructured grid. Nonetheless, the objective of this work is to reproduce 3D experimental
turbulent cases, as chordwise flexible blade for instance. The FSI solver then needs to be
validated with cases of that kind; this will be presented in the next Chapter.

118
 Chapter 6

Validation against 3D experimental

turbulent cases

In this chapter, the developed FSI solver is used to reproduce turbulent experimental case.
At first, the methodology is validated against a case in which a flexible plate clamped behind
a cylinder immersed in a turbulent flow. Afterwards, a simulation of an experiment with
flexible pitching foil is undertaken in order to verify that the FSI solver is able to reproduce
cases of that kind.

Contents
6.1 Kalmbach experimental case: flexible plate behind a cylinder . . . . . . . . . . . . . 119
6.1.1 Case presentation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119
6.1.2 Numerical setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121
6.1.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123
6.2 Hoerner experimental case: flexible pitching foil . . . . . . . . . . . . . . . . . . . . 128
6.2.1 Case presentation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128
6.2.2 Rigid foil simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130
6.2.3 Flexible foil simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134

6.1 Kalmbach experimental case: flexible plate

behind a cylinder
6.1.1 Case presentation
For completness, the FSI solver is now validated on the experimental case performed by
Kalmbach et al. [147], named FSI-PfS-2a. This case deals with a turbulent flow which leads
to large deformation of the structure.
This test case is derived from the 2D case considered in the previous section; it consists of
a cylinder fixed in a water channel. A flexible rubber structure with an attached steel weight
is clamped behind it. The geometry and dimensions are detailed in Fig. 6.1 and Tab. 6.1.
The inflow velocity uinf low is set at 1.385 m.s−1 resulting in nearly symmetrical, large and
reproducible structural displacement. The measured inflow turbulence level is T uinf low = 0.02

119
Chapter 6. Validation against 3D experimental turbulent cases

Fig. 6.1: Sketch of the studied experimental test case. Extracted from [147].

Geometry parameters Value [m]

Cylinder diameter D = 0.022
Cylinder center x-position Lc = 0.077
Cylinder center y-position Hc = 0.120
Test section lenght L = 0.338
Test section height H = 0.240
Test section width W = 0.180
Deformable structure length l1 = 0.050
Deformable structure height h = 0.002
Deformable structure width w = 0.177
Rear mass length l2 = 0.010
Rear mass height h = 0.002
Rear mass width w = 0.177

Table 6.1: Dimensions of the sketch in Fig. 6.1.

120
 6.1. Kalmbach experimental case: flexible plate behind a cylinder

and is considered sufficiently low to be ignored in the rest of the study, a choice motivated
by previous studies which also did not take it into account [148, 149]. The water density is
ρf = 1000 kg.m−3 and the dynamic viscosity is µf = 0.001 Pa s. The structure is composed
of two different materials; the flexible part is made of para-rubber while the bonded rear
mass is made of steel. All properties can be found in Tab. 6.2. The rear mass is here added
to emphasize the structure motion and trigger the second swiveling mode, unlike the first
case FSI-PfS-1a [148].

Young’s
Density Poisson’s
modulus
[kg/m3] ratio
[MPa]
Flexible structure
1090 4.1 0.48
(para-rubber)
Rear mass
7850 2.1 × 105 0.3
(steel)

Table 6.2: Structure properties for the FSI-PfS-2a case.

The Reynolds number based on cylinder diameter gives Re = 30 470 and it is observed
experimentally that the flow is in a sub-critical regime. The boundary layers around the
cylinder are still laminar but the flow becomes turbulent downstream. At this point, a large
variety of spatial and temporal frequencies appears but only the lowest ones can be found
in the structure displacement. The flexible part thus deforms in the second swiveling mode
with a frequency of fF SI = 11.25 Hz.

6.1.2 Numerical setup

In order to reduce computational cost, a subset case has been executed in [148] where
dimensions along z-axis were reduced. Structure width then became w0 = l1 + l2 , and the gap
between channel walls and the structure was ignored such as the test section width W 0 was
equal to w0 . For the full case, these walls were assumed as slip walls while periodic boundary
conditions were used for the subset case. The comparison of the results in terms of structure
deflection but also flow field between both cases has demonstrated that this simplification
was valid. Only the subset case is then considered in the present study with the difference
that the small gap between the flexible part and the channel walls is taken into account. Test
section height then becomes W 0 = l1 + l2 + (W − w). Therefore, solid nodes on these sides are
free and a boundary condition that imposes a zero z-displacement is not mandatory. For the
fluid, slip walls are assumed for these boundaries to maintain the blocking effect and limit
flow recirculation and thus reproduce the same effect than in the experimental setup.
Considering that the flow relatively away from the body will not affect the structure
motion, fluid mesh has been refined in the zone of influence around the cylinder and the plate.
First cell size close to the body is equal to 9 × 10−3 D and slowly decreases to 2.3 × 10−2 D in

121
Chapter 6. Validation against 3D experimental turbulent cases

the rest of this region, as it can be seen in Fig. 6.2. Consequently, a grid composed of about
40 million tetrahedra is generated. As precised above, all channel walls are considered as

Fig. 6.2: Fluid mesh metric field (log scale)(left) and solid mesh (right) used to reproduce
the FSI-PfS-2a case.

slip walls since the full resolution of their boundary layers would be too costly. Nevertheless,
this simplification is expected to have no influence on solid motion. On the other hand,
the cylinder and the structure are defined as no-slip walls. Once again, the timestep of the
simulation is based on the CF L condition leading to ∆t ≈ 5×10−5 s. For the flow simulation,
LES is performed and the SGS model used is the dynamic Smagorinsky model [101]. The
pseudo-solid mesh movement is used with Rmin = 0.0025 m, Rmax = 0.08 m and yr = 100
and is combined with dynamic mesh adaptation.
Thanks to CWIPI (section 4.2.1), the two meshes do not need to be coincident. Thus, the
solid mesh is composed by 900 27-nodes hexahedrons of size equal to structure height h and
can be visualized in Fig. 6.2. The nodes against the cylinder are fixed but all the other ones
are free to move, building five solid faces coupled with the fluid. For temporal advancement,
the generalized-α method is used with a spectral radius ρ∞ = 1.0 (Eq.( 3.85)). In a similar
manner than in [149], no structural damping are applied.
As explained in the previous section, the flow is first established before starting the FSI
coupling. The underrelaxation factor is ω = 0.2 and the FSI convergence criterion εF SI is
equal to 5 × 104 . It results in a mean number of subiterations NF SI = 22.6. A cluster specifi-
cally designed for memory bound applications (eg: CFD codes such as YALES2) was used. It
is based on AMD Epyc 7302 processors with 16 cores and 128Mo of L3 cache, providing high
memory bandwith for each cores. Nodes with 128GB of RAM and two processors are con-
nected with an InfiniBand HDR100 network (100Gb/s). The computation of the FSI-PfS-2a
case was carried out using 240 cores for the fluid and 16 cores for the solid, so that one phys-
ical second could be predicted in 519 hours wall-clock, i.e. 133kH.CPU. The mesh movement
solving took 28.6% of the computational time, the solid solving 52.2%, the fluid 17.2% and

122
 6.1. Kalmbach experimental case: flexible plate behind a cylinder

dy /D|max dy /D|min Mean dy /D Amp. dy /D f dy /D [Hz]

Exp 0.667 -0.629 0.019 0.648 11.25
0.632 -0.627 0.002 0.629 11.73
This study
(5.27%) (0.37%) (86.37%) (2.89%) (4.23%)

Table 6.3: Numerical results and comparison with the experiment.

the DMA 2.0%. Given that the mesh movement solver and the SMS use the same linear
solver, this computational time distribution highlights the impact of the preconditionner as
mentionned in section 3.1.5.

6.1.3 Results
A visualization of the resulting simulation is given in Fig. 6.3 and 6.4 and an anima-
tion can be watched at https://www.youtube.com/watch?v=X7Zb7y67yAo. It shows that
although the structure displacement is mainly two dimensional, three-dimensional flow struc-
tures appear in the wake of the cylinder. Furthermore, the vorticity field highlights the fact
that the structure motion can deviate the vortex trajectories.
In order to compare with experimental data, a probe is positioned on the solid at 2 mm
from the plate tip. To verify whether the displacement can be considered uniform along the
z-axis, normalized y-displacement dy /D is measured at plate center z = 0.0 and extremities
z ± 0.03. The progression over time for those quantities is plotted in Fig. 6.5. It can be seen
that deviation between the three results is minor, confirming that dy /D can be computed
as a mean value of those three points. This mean value is then averaged in 23 sub-parts of
the phase in order to compare with experimental data. Results are collected over 5 periods
leading to a mean standard deviation of about 0.01. The two phase-averaged dy /D signals
are shown in Fig. 6.6 where standard deviation of each point is also represented. This figure
highlights the good agreement between the numerical results and the experimental data.
Table 6.3 gathers the maximum, minimum, amplitudes, mean and frequency values of the
normalized y-displacement.
Note that the asymmetries of experimental data are unexpected considering the symmetry
of the entire setup, but may be explained in [149]. Following careful analysis, it has been
concluded that errors were related with the assembly of the structure, especially the bond
between the rubber plate and the cylinder and the bond between the cylinder and the test
section. That explains the undesired mean value of dy /D of 0.019 while no mean displacement
is expected in this direction due to symmetry. This experimental issue makes irrelevant any
accurate comparison with the reference value for dy /D|max and dy /D|min . However, the
amplitude of the structure oscillation is still valid and the numerical result deviates only of
about 2.9%. This is also true for the predicted frequency where the error is around of 4.2%.
The overall agreement with experimental values confirms that the present FSI solver is able
to predict turbulent case with large deformation with high fidelity.

123
Chapter 6. Validation against 3D experimental turbulent cases

Fig. 6.3: From top to bottom: velocity field, vorticity field and deformed mesh for the FSI-
PfS-2a case. Note that the fields have been cropped but the entire domain can be seen in
Fig. 6.2.

124
 6.1. Kalmbach experimental case: flexible plate behind a cylinder

Fig. 6.4: Visualization of a Q-criterion volumic rendering for the FSI-PfS-2a case. The
structure is colored by the normalized y-displacement dy /D.

Fig. 6.5: Normalized y-displacement at different z locations.

125
Chapter 6. Validation against 3D experimental turbulent cases

Fig. 6.6: Comparison of the numerical and experimental averaged phase of normalized y-
displacement.

More data are given in Fig. 6.7. It can be seen than despite the use of DMA, the solution
is stable in terms of displacement and forces. Also it should be noticed that the remeshing
step does not always results in low maximum skewness mesh. Nonetheless, the low quality
cells can be within a zone which will not be deformed by the structure motion, as it is shown
by the threshold obtained after 0.6 s. Finally, contrary to the case of section 5.4.2, the DMA
does not induce any spike of the NF SI value.
This simulation confirms the ability of the developed FSI solver to reproduce with high
fidelity 3D turbulent cases. This challenging case had been reproduced only twice in the
literature so far [149, 150], and never with unstructured mesh or 3D solid finite elements.
Now that the FSI solver had been validated on several benchmarks, it can be used to
reproduce a case of a chordwise flexible foil. The objective of the following part is then
to verify is the developed FSI solver is able to reproduce such cases, but also to provide a
deeper understanding of the FSI phenomenon occuring when flexible foil are used. Even if
this key step has been initiated at a late stage of the thesis, some results were obtained and
are presented in the next section.

126
 6.1. Kalmbach experimental case: flexible plate behind a cylinder
Lift (N )

2.5
0.0
−2.5
0.5 0.6 0.7 0.8 0.9
2.0
Drag (N )

1.5

1.0

0.5 0.6 0.7 0.8 0.9

0.5 dy /D z=0.0
dy /D

dy /D z=0.03
0.0
dy /D z=-0.03
−0.5 dy /D mean
0.5 0.6 0.7 0.8 0.9
1.0
Skewness

0.9

0.5 0.6 0.7 0.8 0.9

20
NF SI

10

0.5 0.6 0.7 0.8 0.9

time (s)

Fig. 6.7: Results obtained with present computation. From top to botttom: lift and drag
integrated on the cylinder and the plate, dy /D for different z positions, maximum skewness
and number of subiterations.

127
Chapter 6. Validation against 3D experimental turbulent cases

6.2 Hoerner experimental case: flexible pitching foil

6.2.1 Case presentation
Several experiments involving chordwise flexible blade have been presented in section 1.1.
They have been contemplated to test the present FSI solver, but choice of a case with only one
blade seemed more suitable as it would be part of a good first step before reproducing an entire
VAT. The case of Hoerner et al. [34, 43, 37, 44] has then been chosen. Also, this experiment
was performed in the laboratory of the thesis, so that exchanges with experimentators were
facilitated.
The goal of this experiment is to show how hyperflexible blades effectively adapt dynami-
cally their shape and passively control the flow through deformation [44]. In order to generate
high-resolution flow observation around the hydrofoil, only one blade is studied in a water
channel. With VAT configuration, it can be shown [31] that the angle of attack α (defined
in Fig. 1.2) experienced by one blade is varying according to
 sin θ  ωR
α = arctan where λ = , (6.1)
λ + cos θ u∞
with θ the azimuth angle, R the turbine radius, ω the angular velocity and u∞ the inflow
velocity.
To reproduce this variation, the blade is pitched by a drive system, as illustrated in
Fig. 6.8 (and in Fig 1.4). For that, half of one blade extremity is embedded in a rotating
structure. The blade is disposed at 32 mm from the inlet, and its lenght is L = 173, 5 mm

Fig. 6.8: Scheme of the experimental setup. The water tunnel dimensions are 1m × 0.28m ×
0.175m. Extracted from [43].

so that there is only a gap of 1.5 mm between the free blade extremity and the side wall.

128
 6.2. Hoerner experimental case: flexible pitching foil

Four different hydrofoilds with a symmetric NACA0018 geometry were set up, with a
chord of C = 66 mm. One of them was entirely rigid, while the others were composed of
various materials, as shown in Fig 6.9. The leading edges made of a milled aluminum pieces

Fig. 6.9: Photographies of the bioinspired hydrofoil. Extracted from [43].

reach up to the first quarter of the chord. A carbon fiber plate is embedded in this part and
the rest of the trail is composed by a silicone embodiment. Thus, the profile stiffness could
be adjusted by changing the carbon-fiber blade thickness th between the three hydrofoils,
from 0.3 to 0.7 mm. One part of the blade manufacturing as well as the experimental setup
can be visualized at [151].
Parametric studies of rigidity, oscillation frequency and tip speed ratio have then been
performed [44]. Considering that the goal of this part is to verify if the developed FSI solver
is able to reproduce a flexible pitching foil case, the experiment resulting in the largest foil
deformation seems to be a relevant choice. The corresponding parameters are th = 0.3 mm,
f = 2.63 Hz, u∞ = 3.07 m.s−1 and λ = 2. The chord-based Reynolds number then gives
Re = 202620.
To measure blade deformation, the angle of deformation β is defined in Fig 6.10. There-

Fig. 6.10: Definition of the angle of deformation β.

129
Chapter 6. Validation against 3D experimental turbulent cases

fore, for the rigid blade, it gives α = β. Also, lift and drag forces applied on the blade are
measured, as the velocity field around the foil. Figure 6.11 shows the obtained experimental
results. Before trying to reproduce the flexible blade case, the rigid case is first simulated to

Fig. 6.11: Experimental results. Extracted from [44].

ensure the chosen fluid mesh is able to describe accurately the flow.

6.2.2 Rigid foil simulation

In order to reduce the computational cost, the channel width is reduced from 2.65C to
1C. Preliminary tests have indeed shown that this change has not affected significantly the
fluid forces per unit of length. Also, the gap between the extremity of the blade and the
channel wall is not taken into account, and the channel walls are considered slip. The blade
lenght then becomes L = 1C. The computationnal domain is shown in Fig. 6.12.
For the mesh movement, the grid velocity is prescribed in the whole domain, and DMA
is used. This mesh strategy has already been discussed in section 2.2.2.2, in the case of

130
 6.2. Hoerner experimental case: flexible pitching foil

Fig. 6.12: Fluid geometry. Channel width has been reduced in comparison with experiment
from 2.65C to 1C. Also, the gap between the extremity of the blade and the channel wall
has been ignored.

Fig. 6.13: Prescribed grid velocity field for the rigid foil case. The mesh movement strategy
has already been detailed in section 2.2.2.2.

131
Chapter 6. Validation against 3D experimental turbulent cases

the VAT simulation. An example of grid velocity field is given in Fig. 6.13. The gravity is
taken into account. The water density is ρf = 1000 kg.m−3 and the dynamic viscosity is
µf = 0.001 Pa s.
In order to study mesh influence, two meshes with different cell sizes are generated, one
with 27M tetrahedrons and the other one with 74M tetrahedrons. The metric fields of the
two meshes is given in Fig. 6.14. It can be seen that mesh is refined in a large zone after the

Fig. 6.14: Metric fields of the two meshes, 27M (left) and 74M (right).

foil, in order to accurately predict the dynamic stall phenomenon occuring because of angle
of attack variation. Flow is established before pitching the foil, and the two first oscillations
are computed for each mesh. For accurate comparison, results should be averaged on a large
number of periods, but it cannot be done easily at this point because of the ratio between
the timestep dt and the period T . In fact, for the mesh 74M, dt ≈ 9 × 10−6 s so that one
oscillation requires T /dt ≈ 69500 iterations. However, preliminary tests on a larger amount
of oscillations with the mesh 27M showed that results were similar from the second period,
so that comparison with experimental data for t > 1T seems acceptable.
Results are compared in terms of lift and drag per unit of length. A comparison with
experimental data for the two meshes is provided in Fig. 6.15. Hysteresis effects can be
seen as the second period is slightly different from the first one. It can also be observed
that lift and drag curves present a bump after each stall. This effect was already observed
and analyzed in [31]. Nonetheless, it does not appear in the experimental results; other
tests demonstrated that this phenomenon was caused by the absence of gap at the blade
extremity in the computational domain. Despite important differences in cell sizes, both

132
 6.2. Hoerner experimental case: flexible pitching foil

Exp Rigid 27M

74M
Lift/C (N/m)

50

−50
0.00 0.25 0.50 0.75 1.00 1.25 1.50 1.75 2.00
Drag/C (N/m)

20

0
0.00 0.25 0.50 0.75 1.00 1.25 1.50 1.75 2.00

25
α(◦)

−25
0.00 0.25 0.50 0.75 1.00 1.25 1.50 1.75 2.00
time/T

Fig. 6.15: Results comparison between the two meshes and experimental data for the rigid
case.

133
Chapter 6. Validation against 3D experimental turbulent cases

meshes are giving similar results, in good agreement with experimental results. Those results
confirms that the mesh influence is limited and that mesh 27M allows the accurate fluid forces
computation. Therefore, this mesh will be used onwards for the flexible blade simulation.

6.2.3 Flexible foil simulation

For the flexible blade, reproducing real boundary conditions would imply setting displace-
ment only at one extremity and taking into account the three materials composing the blade.
Nonetheless, tests revealed that this configuration requires prohibitive computational time
due to a lack of efficiency of the linear system solving algorithm of the SMS (see section 3.1.5).
One solution was to consider the first quarter of the blade, entirely made of aluminum, as a
perfectly rigid solid. The computational domain was then reduced to the silicone and carbon
parts, and displacement was set on the solid face which is in contact with the aluminium
part. Figure 6.16 illustrates the remaining solid domain, meshed with 10-nodes tetrahedrons.
A first mesh was generated to take clearly into account the distinction between the silicone

Fig. 6.16: Solid meshes composed by 10-nodes tetrahedrons; about 3000 10-nodes tetrahe-
drons for the mesh with the carbon plate (left) and about 1900 10-nodes tetrahedron without
the plate.

and the 3 mm carbone part. However, the differences in Young modulus values between these
cells induce high condition number of the linear system, so that the solving requires a large
number of iterations (more details in [129]). A similar situation has already been observed
in section 6.1 due the steel weight linked to the rubber. It has finally been decided to not
take into account explicitly the carbon plate in the solid domain, but to use instead a single
material with physical properties between the silicon and the carbon, as given in Tab 6.17.
This simplification made the computation 29 times faster. Several values have been tested
but finally, the Young modulus was chosen as E = 1.28 MPa, the density as ρ = 1250 kg.m−3 ,
and the resulting mesh of 1900 10-nodes tetrahedrons is shown in Fig. 6.16. The poisson ratio
was set to zero to prevent foil distension in spanwise direction.
Furthermore, the linear solving (see section 3.1) of the structure caused instabilities for
the FSI coupling. A possible explanation is that the update of the interface fluid-solid is

134
 6.2. Hoerner experimental case: flexible pitching foil

Fig. 6.17: Physical properties of the different materials.

mandatory for the solid, to ensure consistency with fluid pressure normal to the boundary.
Consequently, non linear solving was necessary.
As regards the FSI coupling, it should be noticed that only the part of the blade corre-
sponding to the structure computational domain was coupled, and as a result, the first quarter
of the blade appeared only in the fluid domain and its movement was prescribed. The un-
derrelaxation factor is ω = 0.5 and the FSI convergence criterion εF SI is equal to 5 × 10−4 .
Moreover, in order to stabilize the FSI coupling and reduce the number of subiterations NF SI
required, structural damping is added with the stiffness proportional damping factor β equal
to 0.05 and the mass proportional damping factor α equal to 0.0 (see Eq.( 3.71)).
It finally results in a mean number of subiterations NF SI = 3.39. The cluster presented in
section 6.1.2 was also used, but this time with 20 cores for the solid and 225 for the fluid, so
that one period of foil oscillation could be predicted in 53.25 hours wall-clock, i.e 13kH.CPU.
It is only 6 times more than the computational cost of the rigid case simulated with the same
fluid mesh. The mesh movement took 68.7% of the computational time, the solid solving 9%,
the fluid solving 19.1% and the DMA 3.2%.
The first thing to consider is to verify that the chosen material properties allow to re-
produce what has been observed experimentally. For that, numerical results are compared
with experimental data in terms of lift, drag and deformation angle β in Fig. 6.18. It can
be seen that the results are in good agreement with the experiment for the forces, and that
the foil deformation is only slightly over estimated. The flow behaviour have then been well
reproduced, allowing physical insights given that the solid problem simplification does not
affect the flow configuration once the blade behaviour is similar. For the analysis of the FSI
phenomenon, different post-processings have been done.
An animation showing the results can be watched at https://www.youtube.com/watch?
v=d1A6ShfmbnM. At first, a visualization of the vorticity volumic rendering for different angle
of attack is given in Fig. 6.19. It compares results computed for the rigid and flexible cases.
For the rigid case, the dynamic stall can clearly be seen. It causes the generation of a first
vortice above the foil. This will create a low pressure region, attracting the surrounding fluid,
especially fluid coming from the other side of the foil. This finally lead to the generation of

135
Chapter 6. Validation against 3D experimental turbulent cases

Exp Flexible Y2

40
Lift/L (N/m)

20

−20

−40
0 1 2 3 4 5 6

20
Drag/L (N/m)

15

10

−50 1 2 3 4 5 6

α
20
Exp β
angle (◦)

Y2 β
0

−20

0 1 2 3 4 5 6
time/T

Fig. 6.18: Comparison of forces and blade deformation with the experimental data for the
flexible case.

136
 6.2. Hoerner experimental case: flexible pitching foil

Fig. 6.19: Vorticity volumic renderings of rigid (left) and flexible (right) foil cases.
137
Chapter 6. Validation against 3D experimental turbulent cases

a second vortice at the trailing edge. That having been said, the analysis of the flexible case
revealed that the generation of this second vortice induces the foil deformation. This second
low pressure region appears indeed where the blade thickness is relatively small, and the
pressure difference is sufficiently strong to deform significantly the foil. It finally results in
completely different flow configuration.
This explanation is confirmed by the analysis of the velocity and pressure fields shown in
Fig. 6.20. The two vortices can be observed as well as the low pressure regions. The pressure

Fig. 6.20: Velocity, skewness and pressure fields computed for the rigid (top) and the flexible
(bottom) foil at α = 29.1◦ .

distribution appplied on the blade is then very different between the rigid and flexible case,
explaining the differences in lift and drag seen in Fig.6.11. Also, the skewness fields highlight
the difference of method used for mesh movement between the two cases, given that it is
more complicated to move the fluid mesh when the blade is deforming.
Foil relative displacement compared to the rigid case for an angle of attack of α = 29.1◦
is given in Fig. 6.21. It can be noticed that at the maximum deflection, deformation reaches
35% of the chord, confirming that the solver can reproduced cases with large foil deformation.
This simulation then confirms that the developed FSI solver is able to reproduce cases
of chordwise flexible blade. These different analysis lead to a better understanding of the
link between the dynamic stall and the flexible foil deformation. This also shows why a
LES approach is crucial to perform accurate analysis of the flow. Further works would be
necessary to complete this study and explain how exactly this phenomenon influences the

138
 6.2. Hoerner experimental case: flexible pitching foil

Fig. 6.21: Relative displacement of the foil compared to the rigid case at α = 29.1◦ .

turbine performances, but this could not be done in this work given the time requirements.

In the computation with the carbon plate, the FSI solver spent more than 90% of the
computational time in the linear solving of the SMS. This highlights the need to enhance
performances of the linear solver used by the SMS, as mentionned in section 3.1.5, to be
able to reproduce cases with multi-materials solid. Nevertheless, the solutions are known
(see section 7.2) and will be implemented in future works. Besides, once the simplification of
the solid problem was made, this test proved that the developed method for the FSI solver
allows indeed to reproduce cases of flexible pitching foil with large deformation, which has
never been done in litterature so far with a 3D LES approach. It will then be used in future
works to perform further physical analysis of this case, as well as the simulation of flexible
blades turbine as part of the ANR project DYNEOL.

139
 Chapter 7

Conclusion & Perspectives

7.1 Conclusion
Considering the necessity of energy transition, the importance to develop and enhance low
carbon power sources as renewable energies keeps rising. In this context, recent experimental
studies tend to show that the use of chordwise flexible blades for vertical axis turbines en-
hances their efficiencies and increases their lifetime. These potential improvements of turbine
efficiency seem to be explained by the impact of the chordwise flexibility on dynamic stall
which implies a drag reduction.
However, the FSI phenomenon occuring in these cases remains unclear in some aspects.
Numerical simulations of these experiments would allow a better understanding of this effect
and would open new perspectives for future improvements in the blade design. As a result,
the main motivation of this thesis is the numerical study of VAT with chordwise flexible
blades.
Although FSI phenomena are very common in nature and human activities, their numer-
ical simulations remain a challenging task. It requires advanced solvers in different physic
fields. Moreover, the resulting coupling can be very unstable, especially when fluid and
solid densities are close. Different approaches are possible to face this problem, but cases
of chordwise flexible blade involve a FSI solver with specific characteristics. For the fluid,
it is necessary to use body fitted meshes to describe accurately the boundary layer at high
Reynolds number. As regards the solid, the foil geometry imposes the use of 3D solid ele-
ment in order to compute structure displacements. Finally, both solvers have to be strongly
coupled to ensure consistency between fluid and solid solutions. These requirements explain
the limited number of numerical studies about flexible foil.
Some FSI solvers have yet been developed to face this challenge, but a RANS approach
is systematically used. This do not always provide reliable results because the pressure
distribution is not accurate enough to conveniently determine dynamic stall separation point
position. This is why Hoerner et al. [34] suggest that three-dimensional LES approach is
required to reproduce with high fidelity a high Reynolds number case involving a chordwise
flexible blade. In fact, very few LES-based FSI solvers have been developed, and they never
used 3D solid elements, making a chordwise flexible foil simulation impossible. To the best
of author’s knowledge, LES-based FSI solvers meeting all the previous requirements do not
exist in the literature so far.
The aim of this work is then to develop a high fidelity FSI solver, able to reproduce a wide
variety of FSI configurations, especially cases involving chordwise flexible blades. This solver

141
Chapter 7. Conclusion & Perspectives

uses LES approach to predict fluid dynamics and 3D solid finite elements for solid dynamics,
all on unstructured meshes (Chapter 1).
All the implementations have been carried out within YALES2, a multi-physics library,
initially designed for fluid mechanics. In order to predict the flow behaviour on moving grids,
specific numerical methods are used by the ALE solver. They have been presented, as well as
a mesh movement strategy using Dynamic Mesh Adaptation (DMA), allowing to reproduce
cases of moving or rotating objects. A simulation of a four rigid blades VAT has then been
performed, confirming that the ALE solver is effective and can be used for FSI coupling
(Chapter 2).
Unlike the fluid solver, the solid solver has been developed from scratch in this work. This
has represented one of the main task of this thesis. The FEM has been explained in details,
and the SMS has been validated carefully with several tests (Chapter 3).
Also, body fitted techniques involve a fluid mesh movement computation. This can be
difficult for large deformations, especially with unstructured grids. An original and efficient
pseudo-solid method has then been proposed, meeting all the underlying requirements. The
FSI coupling scheme has also been introduced step by step using a simple test case. The
benefit of strong coupling has then been highlighted. The potential use of all the previous
elements combination, including DMA, has thus been shown (Chapter 4).
Nonetheless, the developed FSI solver was requiring validation. For this purpose, a refer-
ence 2D numerical benchmark with laminar flow has been reproduced. For each cases, mesh
convergence studies have been performed. After validation of fluid and solid solvers inde-
pendently, the coupling has been validated against two FSI cases, with and without DMA
(Chapter 5).
The methodology has then been applied to 3D experimental turbulent cases. The Kalm-
bach case has been studied. It consists of a flexible rubber structure with an attached steel
weight clamped behind a cylinder. Once the structure is immersed in high Reynolds num-
ber flow, large deformations of the structure are obtained. This case has been succesfully
reproduced even though it had never been simulated with unstructured fluid mesh or 3D
solid finite elements by the past. That confirms the ability of the developed FSI solver to
reproduce with high fidelity 3D turbulent cases.
Simulation of experiment involving chordwise flexible blade has been finally made at the
end of the thesis. Even if simplification has been made for the solid problem, this confirms
that the FSI solver was able to reproduce cases of flexible foil, completing the initial goal
the thesis. The use of the LES approach also allowed physical analysis, providing deeper
understanding of the link between the dynamic stall and the foil deformation. Simulation of
this kind had never been done with a 3D LES approach. This work has demonstrated thus
the potential of the FSI solver for its intended use (Chapter 6).
The development of a high fidelity FSI solver seems a sucess, and despite possible im-
provements detailed in the next sections, it can be used as a tool to reproduce a wide variety
of FSI configurations.
The PhD defense can be watched at https://www.youtube.com/watch?v=Kc7vQDwIQP8&
t=2422s.

142
 7.2. Perspectives

7.2 Perspectives
As mentionned above, the principle axis of improvement is about computational time
reduction. Most possible improvements are actually related to the SMS. In fact, this solver
has been developed from scratch to ensure compatibility between solvers and to easily modify
the FSI solver for future works. However, all features and advanced methods used in the
literature for computational structural dynamics could not be implemented in this work.
That is why the solid solver can be upgraded with well-known methods.
For faster computation, the linear solving algorithm needs to be improved with a better
preconditionner. As detailed in section 3.1.5, One relevant choice providing many opportu-
nities would be to make available the use of the PETSc library by the SMS [130, 131, 132].
Also, it would be useful to implement more advanced techniques than the basic Newton algo-
rithm used presently for non linear solving [119]. Finally, development should be undertaken
to allow the use of different types of finite elements in the computational domain, or also
other types, as shell or membrane elements which are well suited to describe thin structures.
As regards the mesh movement method, the pseudo-solid algorithm could be improved.
Its computational cost can indeed be very important, as in section 6.1. Inspired from existing
techniques [141], implementations to compute mesh deformation on a coarse grid and then
interpolate on the computational grid have been undertaken, but not completed. In addition,
other mesh movement methods for unstructured grid, as techniques based on Radial Basis
Functions (RBF), could be tested and may result in lower computational cost [138].
Another aspect slowing down computations is the large number of FSI subiterations NF SI
required into the FSI loop. Only Aitken method has been tested in this work, but a lot of
other schemes for partitioned simulation exist [56, 152, 153, 154]. A possible solution would
be to combine YALES2 with the preCICE library [155].
Further works could aim at implementing the features to simulate simultaneously several
FSI with different objects in the same domain. For instance, DMA could be used to reproduce
several flexible falling objects, which is very difficult with body fitted techniques. Also, this
would be necessary to simulate an entire flexible blades VAT.

This work will also be used as a base for future studies, considering that the selected
methods can be applied to a wide variety of FSI cases. Further works will be undertaken in
order to complete the flexible foil simulations, and provide new physical insights on the use
of chordwise flexible blade for VAT applications.
On the other hand, at IMAG (Institut Montpelliérain Alexander Grothendieck), a thesis
has been started by B. Thibaud about numerical simulation of hemodynamics in deep vein
valves. In this work, fluid dynamics and structure dynamics are at stake respectively for
blood and valve tissues. The FSI solver is then used to predict the deformations of the
leaflets under real blood pressure condition, as illustrated in Fig. 7.1.

143
Chapter 7. Conclusion & Perspectives

Fig. 7.1: Simulations performed at IMAG with the FSI solver. Top: Leaflets in closed and
open positions, the arrow represents the blood effect. Bottom: Flow velocity field in the
middle place of the vein. Courtesy B. Thibaud.

144
 7.2. Perspectives

The work performed in this thesis has thus provided an advanced tool to the YALES2
community to lead new research studies in the domain of structural mechanics and fluid-
structure interaction.

145
 Appendix A

Shape functions and Gauss points

147
Appendix A. Shape functions and Gauss points
Version
Code_Aster default
Titre : Fonctions de forme et points d'intégration des élé[...] Date : 17/07/2017 Page : 10/28
Responsable : DELMAS Josselin Clé : R3.01.01 Révision :
ad6c457e7f2a

3.2 Quadrangles : ELREFE QU4, QU8, QU9

N4 N7 N3


N8 N6
N9 

N1 N5 N2
Coordonnées des nœuds :
 
N1 -1.0 -1.0
N2 1.0 -1.0
N3 1.0 1.0
N4 -1.0 1.0
N5 0.0 -1.0
N6 1.0 0.0
N7 0.0 1.0
N8 -1.0 0.0
N9 0.0 0.0

Famille Point   Poids

FPG1 1 0 0 4
FPG4 1 −a −a 1.0
2 a −a 1.0
3 a a 1.0
4 −a a 1.0
a=1/  3
FPG9 1 −a −a 25/81
2 a −a 25/81
3 a a 25/81
4 −a a 25/81
5 0.0 −a 40/81
6 a 0.0 40/81
7 0.0 a 40/81
8 −a 0.0 40/81
9 0.0 0.0 64/81
a=0.774596669241483

QU4 : quadrangle à 4 nœuds

nombre de nœuds :4
nombre de nœuds sommets :4

Manuel de référence Fascicule r3.01: Références générales

Document diffusé sous licence GNU FDL (http://www.gnu.org/copyleft/fdl.html)

148
 Version
Code_Aster default
Titre : Fonctions de forme et points d'intégration des élé[...] Date : 17/07/2017 Page : 11/28
Responsable : DELMAS Josselin Clé : R3.01.01 Révision :
ad6c457e7f2a

fonctions de forme, dérivées premières et secondes du quadrangle à 4 nœuds :

{N} { ∂ N /∂  } { ∂ N /∂ }
1−1− / 4 − 1− /4 − 1−/ 4
11− / 4 1− / 4 − 1/ 4
11 / 4 1 / 4 1/4
1−1 / 4 − 1 /4 1−ξ /4

{ ∂2 N /∂ 2 } {∂2 N /∂ ∂ } {∂2 N /∂  2 }

0 1/4 0
0 -1/4 0
0 1/4 0
0 -1/4 0

QU8 : quadrangle à 8 nœuds

nombre de nœuds :8
nombre de nœuds sommets :4

fonctions de forme et dérivées premières du quadrangle à 8 nœuds :

{N} { ∂ N /∂ } { ∂ N /∂ }
1−1− −1−− / 4 1−  2  / 4 1−2  /4
11− −1− / 4 1−  2 − / 4 − 1 −2  / 4
11 −1 / 4 1  2  / 4 12  /4
1−1 −1− / 4 − 1 −2 /4 1−−2  /4
1−2  1− /2 − 1−  − 1−2 / 2
11− 2 /2 1− 2 / 2 − 1
1−2  1 /2 − 1  1− 2 /2
1−1− 2 /2 − 1− 2 /2 − 1−

dérivées secondes du quadrangle à 8 nœuds :

{ ∂2 N /∂ 2 } {∂2 N /∂ ∂ η } {∂2 N /∂  2 }
1−η /2 1−2 ξ −2 η/4 1−/2
1−η /2 −12 ξ −2 η / 4 1/2
1η /2 12 ξ 2 η /4 1/2
1η /2 −1−2 ξ 2 η / 4 1−/2
−1η ξ 0
0 −η −1−
−1−η −ξ 0
0 η −1

QU9 : quadrangle à 9 nœuds

nombre de nœuds :9
nombre de nœuds sommets :4

Manuel de référence Fascicule r3.01: Références générales

Document diffusé sous licence GNU FDL (http://www.gnu.org/copyleft/fdl.html)

149
Appendix A. Shape functions and Gauss points
Version
Code_Aster default
Titre : Fonctions de forme et points d'intégration des élé[...] Date : 17/07/2017 Page : 12/28
Responsable : DELMAS Josselin Clé : R3.01.01 Révision :
ad6c457e7f2a

fonctions de forme et dérivées premières du quadrangle à 9 nœuds :

{N } { ∂ N /∂ ξ } { ∂ N /∂ }
  −1 −1 / 4  2 −1  −1/ 4 −1  2 −1 / 4
  1 −1 / 4  2 1  −1/ 4 1  2 −1 / 4
  1 1 / 4  2 1  1/ 4 1  2 1 / 4
  −1 1 / 4  2 −1  1/ 4 −1  2 1 / 4
2 2
1−  −1/2 −   −1 1−  2 −1 /2
1 1− 2 / 2  2 1 1− 2 /2 −  1 
1−2  1/2 −   1 1− 2  2 1 /2
2 2
−1 1− / 2  2 −1 1− /2 −  −1 
1−  1− 2 
2
−2 1−  2
−2  1−2 

dérivées secondes du quadrangle à 9 nœuds :

{ ∂2 N /∂ 2 } {∂2 N /∂ ∂ } {∂2 N /∂  2 }
  −1 /2 −1 /2  −1 /2  −1 /2
  −1 /2 1 /2  −1 /2  1 /2
  1 /2 1 /2  1 /2  1 /2
  1 /2 −1 /2  1 /2  −1 /2
−  −1 − 2 −1 1− 2
1− 2 −  21 − 1 
2
−  1 − 2η1  1−
2
1− −  2−1 − −1 
−2 1− 2  4  −2 1−2 

Manuel de référence Fascicule r3.01: Références générales

Document diffusé sous licence GNU FDL (http://www.gnu.org/copyleft/fdl.html)

150
 Version
Code_Aster default
Titre : Fonctions de forme et points d'intégration des élé[...] Date : 17/07/2017 Page : 19/28
Responsable : DELMAS Josselin Clé : R3.01.01 Révision :
ad6c457e7f2a

4.3 Hexaèdres : ELREFE HE8 , H20 , H27

N5 N8

N6
N7

N4

N2 N3
x

N5 N20 N8

N17 N26 N19

z = +1

N6 N18 N7

N13 N25 N16

N22 N27 N24

z=0

N14 N23 N15

N1 N12 N4

N9 N21 N11
z = -1

N2 N10 N3

Manuel de référence Fascicule r3.01: Références générales

Document diffusé sous licence GNU FDL (http://www.gnu.org/copyleft/fdl.html)

151
Appendix A. Shape functions and Gauss points
Version
Code_Aster default
Titre : Fonctions de forme et points d'intégration des élé[...] Date : 17/07/2017 Page : 20/28
Responsable : DELMAS Josselin Clé : R3.01.01 Révision :
ad6c457e7f2a

Coordonnées des nœuds :

x y z
N1 -1. -1. -1.
N2 1. -1. -1.
N3 1. 1. -1.
N4 -1. 1. -1.
N5 -1. -1. 1.
N6 1. -1. 1.
N7 1. 1. 1.
N8 -1. 1. 1.
N9 0. -1. -1.
N10 1. 0. -1.
N11 0. 1. -1.
N12 -1. 0. -1.
N13 -1. -1. 0.
N14 1. -1. 0.
N15 1. 1. 0.
N16 -1. 1. 0.
N17 0. -1. 1.
N18 1. 0. 1.
N19 0. 1. 1.
N20 -1. 0. 1.
N21 0. 0. -1.
N22 0. -1. 0.
N23 1. 0. 0.
N24 0. 1. 0.
N25 -1. 0. 0.
N26 0. 0. 1.
N27 0. 0. 0.

Fonctions de forme :

Formule à 8 nœuds
1 1
w 1=  1−x   1− y   1−z  w 5=  1−x   1− y   1z 
8 8
1 1
w 2=  1x   1− y   1−z  w 6=  1x   1− y   1 z 
8 8
1 1
w 3=  1x   1 y   1− z  w 7=  1x   1 y   1 z 
8 8
1 1
w 4=  1−x   1 y   1−z  w 8=  1−x   1 y   1z 
8 8

Manuel de référence Fascicule r3.01: Références générales

Document diffusé sous licence GNU FDL (http://www.gnu.org/copyleft/fdl.html)

152
 Version
Code_Aster default
Titre : Fonctions de forme et points d'intégration des élé[...] Date : 17/07/2017 Page : 21/28
Responsable : DELMAS Josselin Clé : R3.01.01 Révision :
ad6c457e7f2a

Formule à 20 nœuds
1 1
w 1=  1−x   1− y   1−z   −2− x− y−z  w 11=  1− x 2   1 y   1−z 
8 4
1 1
w 2=  1x   1− y   1−z   −2x− y−z  w 12=  1− y 2  1−x   1−z 
8 4
1 1
w 3=  1x   1 y   1− z   −2x y−z  w 13=  1− z 2   1−x   1− y 
8 4
1 1
w 4=  1−x   1 y   1−z   −2−x y− z  w 14=  1−z 2   1x   1− y 
8 4
1 1
w 5=  1−x   1− y   1z   −2− x− yz  w 15=  1− z   1x   1 y 
2
8 4
1 1
w 16=  1− z   1−x   1 y 
2
w 6=  1x   1− y   1 z   −2x− yz 
8 4
1 1
w 7=  1x   1 y   1 z   −2x yz  w 17=  1− x 2   1− y   1z 
8 4
1 1
w 8=  1−x   1 y   1z   −2− x yz  w 18=  1− y 2  1x   1z 
8 4
1 1
w 9=  1−x 2   1− y   1−z  w 19=  1− x 2   1 y   1z 
4 4
1 1
w 10=  1− y   1x   1−z  w 20=  1− y   1−x   1z 
2 2
4 4

Manuel de référence Fascicule r3.01: Références générales

Document diffusé sous licence GNU FDL (http://www.gnu.org/copyleft/fdl.html)

153
Appendix A. Shape functions and Gauss points
Version
Code_Aster default
Titre : Fonctions de forme et points d'intégration des élé[...] Date : 17/07/2017 Page : 22/28
Responsable : DELMAS Josselin Clé : R3.01.01 Révision :
ad6c457e7f2a

Formule à 27 nœuds
1 1
w 1= x  x−1  y  y−1  z  z −1  w 15= x  x1  y  y1   1− z 
2
8 4
1 1
w 2= x  x1  y  y−1  z  z −1  w 16= x  x−1  y  y1   1− z 2 
8 4
1 1
w 3= x  x1  y  y1  z  z−1  w 17=  1−x 2  y  y−1  z  z1 
8 4
1 1
w 4= x  x−1  y  y1  z  z −1  w 18= x  x1   1− y 2  z  z 1 
8 4
1 1
w 5= x  x−1  y  y−1  z  z1  w 19=  1−x  y  y1  z  z1 
2
8 4
1 1
w 20= x  x−1   1− y  z  z1 
2
w 6 = x  x1  y  y−1  z  z 1 
8 4
1 1
w 7 = x  x1  y  y1  z  z 1  w 21=  1−x 2  1− y 2  z  z −1 
8 2
1 1
w 8 = x  x−1  y  y1  z  z 1  w 22=  1−x 2  y  y−1   1−z 2 
8 2
1 1
w 9 =  1− x 2 y  y−1  z  z −1  w 23= x  x1   1− y 2  1−z 2 
4 2
1 1
w 10 = x  x1   1− y  z  z −1  w 24=  1−x  y  y1   1−z 
2 2 2
4 2
1 1
w 11 =  1−x  y  y1  z  z−1  w 25= x  x−1   1− y  1−z 
2 2 2
4 2
1 1
w 12 = x  x−1   1− y 2  z  z−1  w 26=  1−x 2  1− y 2  z  z 1 
4 2
1 w 27= 1−x   1− y  1− z 
2 2 2
w 13 = x  x−1  y  y−1   1−z 
2
4
1
w 14 = x  x1  y  y−1   1− z 
2
4
Formule de quadrature de Gauss à 2 points dans chaque direction (ordre 3) (FPG8)
Point x y z Poids
1 −1 /  3 −1 /  3 −1 /  3 1.
2 −1 /  3 −1 /  3 1/  3 1.
3 −1 /  3 1/  3 −1 /  3 1.
4 −1 /  3 1/  3 1/  3 1.
5 1/  3 −1 /  3 −1 /  3 1.
6 1/  3 −1 /  3 1/  3 1.
7 1/  3 1/  3 −1 /  3 1.
8 1/  3 1/  3 1/  3 1.

Manuel de référence Fascicule r3.01: Références générales

Document diffusé sous licence GNU FDL (http://www.gnu.org/copyleft/fdl.html)

154
 Version
Code_Aster default
Titre : Fonctions de forme et points d'intégration des élé[...] Date : 17/07/2017 Page : 23/28
Responsable : DELMAS Josselin Clé : R3.01.01 Révision :
ad6c457e7f2a

Formule de quadrature de Gauss à 3 points dans chaque direction (ordre 5) : ( FPG27 )

Point x y z Poids
1 − − − c 31
2
2 − − 0. c1 c2
3 − −  c 31
4 − 0. − c 21 c 2
2
5 − 0. 0. c1 c2
2
6 − 0.  c1 c2
7 −  − c 31
8 −  0. c 21 c 2
9 c 31
−  
2
10 0. − − c1 c2
11 0. − 0. c 1 c 22
12 0. −  c 21 c 2
2
13 0. 0. − c1 c2
14 0. 0. 0. c 32
2
15 0. 0.  c1 c2
2
16 0.  − c1 c2
17 0.  0. c 1 c 22
18 c 21 c 2
0.  α
3
19  − − c1
2
20  − 0. c1 c2
21  −  c 31
22  0. − c 21 c 2
2
23  0. 0. c1 c2
2
24  0.  c1 c2
3
25   − c1
2
26 α  0. c1 c2
27    c 31

avec :

 3 5 8
= c 1= c 2=
5 9 9

Manuel de référence Fascicule r3.01: Références générales

Document diffusé sous licence GNU FDL (http://www.gnu.org/copyleft/fdl.html)

155
 Bibliography

[1] “Ipcc, climate change 2021, the physical science basis.” https://www.ipcc.
ch/report/ar6/wg1/downloads/report/IPCC_AR6_WGI_Full_Report.pdf.
2 citations pages 5 and 13

[2] “Ipcc, climate change 2021, the physical science basis, summary for pol-
icymakers.” https://www.ipcc.ch/report/ar6/wg1/downloads/report/IPCC_AR6_
WGI_SPM.pdf. Cited page 5

[3] P. Laut, “Solar activity and terrestrial climate: an analysis of some purported correla-
tions,” Journal of Atmospheric and Solar-Terrestrial Physics, vol. 65, no. 7, pp. 801–
812, 2003. Cited page 5

[4] A. Berger, “Milankovitch theory and climate,” Reviews of geophysics, vol. 26, no. 4,
pp. 624–657, 1988. Cited page 6

[5] P. Huybers and W. Curry, “Links between annual, milankovitch and continuum tem-
perature variability,” Nature, vol. 441, no. 7091, pp. 329–332, 2006. Cited page 6

[6] “Science etonnante website.” https://scienceetonnante.com/2016/07/22/

les-cycles-de-milankovitch-et-les-changements-climatiques. Cited page 6

[7] M. Wild, D. Folini, C. Schär, N. Loeb, E. G. Dutton, and G. König-Langlo, “The

global energy balance from a surface perspective,” Climate dynamics, vol. 40, no. 11-
12, pp. 3107–3134, 2013. 2 citations pages 6 and 7

[8] “Jm jancovici, website.” https://jancovici.com/en. 2 citations pages 7 and 12

[9] K. R. Lang, “Nasa’s cosmos.” https://ase.tufts.edu/cosmos/view_chapter.asp?

id=21&page=1. Cited page 7

[10] “Le reveilleur, website.” https://www.lereveilleur.com/climat-du-passe.

Cited page 7

[11] “Ipcc, climate change 2021, the physical science basis, resume a l intention des
decideurs.” https://archive.ipcc.ch/pdf/assessment-report/ar5/syr/AR5_SYR_
FINAL_SPM_fr.pdf. Cited page 7

[12] “The shift data portal.” https://www.theshiftdataportal.org.

2 citations pages 8 and 12

157
Bibliography

[13] “Ipcc, ar5 climate change 2014: Impacts, adaptation and vulnerability.” https://www.
ipcc.ch/report/ar5/wg2/summary-for-policymakers. 2 citations pages 7 and 9

[14] C. Raymond, T. Matthews, and R. M. Horton, “The emergence of heat and humidity
too severe for human tolerance,” Science Advances, vol. 6, no. 19, p. eaaw1838, 2020.
Cited page 7

[15] J. Stephenson, K. Newman, and S. Mayhew, “Population dynamics and climate change:
what are the links?,” Journal of Public Health, vol. 32, no. 2, pp. 150–156, 2010.
Cited page 7

[16] D. B. Botkin, H. Saxe, M. B. Araujo, R. Betts, R. H. Bradshaw, T. Cedhagen, P. Ches-

son, T. P. Dawson, J. R. Etterson, D. P. Faith, et al., “Forecasting the effects of global
warming on biodiversity,” Bioscience, vol. 57, no. 3, pp. 227–236, 2007. Cited page 7

[17] J. R. Malcolm, C. Liu, R. P. Neilson, L. Hansen, and L. Hannah, “Global warming

and extinctions of endemic species from biodiversity hotspots,” Conservation biology,
vol. 20, no. 2, pp. 538–548, 2006. Cited page 7

[18] M. Bui, C. S. Adjiman, A. Bardow, E. J. Anthony, A. Boston, S. Brown, P. S. Fennell,

S. Fuss, A. Galindo, L. A. Hackett, et al., “Carbon capture and storage (ccs): the
way forward,” Energy & Environmental Science, vol. 11, no. 5, pp. 1062–1176, 2018.
Cited page 7

[19] R. S. Haszeldine, “Carbon capture and storage: how green can black be?,” Science,
vol. 325, no. 5948, pp. 1647–1652, 2009. Cited page 7

[20] J.-F. Bastin, Y. Finegold, C. Garcia, D. Mollicone, M. Rezende, D. Routh, C. M.

Zohner, and T. W. Crowther, “The global tree restoration potential,” Science, vol. 365,
no. 6448, pp. 76–79, 2019. Cited page 7

[21] R. Bettin, “De la compensation du carbone au financement de la neutralité,” in Annales

des Mines-Responsabilite et environnement, no. 2, pp. 52–54, FFE, 2021. Cited page 7

[22] “Our world in data.” https://ourworldindata.org/emissions-by-sector#

energy-electricity-heat-and-transport-73-2. Cited page 10

[23] “Ipcc, technology-specific cost and performance parameters.” https://www.ipcc.ch/

site/assets/uploads/2018/02/ipcc_wg3_ar5_annex-iii.pdf. Cited page 11

[24] “European environment agency, co2 emission intensity.”

3f6dc9e9e92b45b9b829152c4e0e7ade. Cited page 11

[25] D. Dussaux, F. Vona, A. Dechezleprêtre, et al., Carbon Offshoring: Evidence from

French Manufacturing Companies. OFCE, 2020. Cited page 11

158
 Bibliography

[26] J. Webb, “Climate change and society: The chimera of behaviour change technologies,”
Sociology, vol. 46, no. 1, pp. 109–125, 2012. Cited page 11

[27] C. Dugast and A. Soyeux, “Faire sa part,” Pouvoir et responsabilité des individus, des
entreprises et de l’état face Al’urgence climatique. Carbone, vol. 4, 2019. Cited page 11

[28] M. Auzanneau, “L’inexorable déclin du pétrole,” Futuribles, no. 4, pp. 65–74, 2021.
Cited page 12

[29] “Futurs energetiques 2050 - principaux resultats.” https://assets.rte-france.

com/prod/public/2021-10/Futurs-Energetiques-2050-principaux-resultats_
0.pdf. Cited page 13

[30] “Note: precisions sur les bilans co2 etablis dans le bilan previsionnel et les etudes
associees.” https://assets.rte-france.com/prod/public/2020-06. Cited page 13

[31] N. Guillaud, Simulation et optimisation de forme d’hydroliennes à flux

transverse. PhD thesis, Université Grenoble Alpes (ComUE), 2017.
4 citations pages 15, 36, 128, and 132

[32] G. Feng, T. De Troyer, and M. Runacres, “Optimizing land use for wind farms using
vertical axis wind turbines,” in EWEA Annual Event, 2014, p. PO 192, European Wind
Energy Association (EWEA), 2014. Cited page 15

[33] J. O. Dabiri, “Potential order-of-magnitude enhancement of wind farm power density

via counter-rotating vertical-axis wind turbine arrays,” Journal of renewable and sus-
tainable energy, vol. 3, no. 4, p. 043104, 2011. Cited page 15

[34] S. Hoerner, S. Abbaszadeh, T. Maı̂tre, O. Cleynen, and D. Thévenin, “Char-

acteristics of the fluid–structure interaction within darrieus water turbines with
highly flexible blades,” Journal of Fluids and Structures, vol. 88, pp. 13–30, 2019.
7 citations pages 16, 17, 18, 23, 72, 128, and 141

[35] G. V. Lauder and P. G. Madden, “Fish locomotion: kinematics and hydrodynam-

ics of flexible foil-like fins,” Experiments in Fluids, vol. 43, no. 5, pp. 641–653, 2007.
2 citations pages 16 and 19

[36] H. Truong, T. Engels, D. Kolomenskiy, and K. Schneider, “A mass-spring fluid-

structure interaction solver: Application to flexible revolving wings,” Computers &
Fluids, vol. 200, p. 104426, 2020. 2 citations pages 16 and 19

[37] S. Hoerner, C. Bonamy, O. Cleynen, T. Maı̂tre, and D. Thévenin, “Darrieus

vertical-axis water turbines: deformation and force measurements on bioinspired
highly flexible blade profiles,” Experiments in Fluids, vol. 61, pp. 1–17, 2020.
3 citations pages 16, 17, and 128

159
Bibliography

[38] M. T. Bouzaher, B. Guerira, and M. Hadid, “Performance analysis of a vertical axis

tidal turbine with flexible blades,” Journal of Marine Science and Application, vol. 16,
no. 1, pp. 73–80, 2017. Cited page 17

[39] P. Chougule and S. Nielsen, “Overview and design of self-acting pitch control mechanism
for vertical axis wind turbine using multi body simulation approach,” in Journal of
physics: conference series, vol. 524, p. 012055, IOP Publishing, 2014. Cited page 17

[40] Y.-b. Liang, L.-x. Zhang, E.-x. Li, and F.-y. Zhang, “Blade pitch control of straight-
bladed vertical axis wind turbine,” Journal of Central South University, vol. 23, no. 5,
pp. 1106–1114, 2016. Cited page 17

[41] D. H. Zeiner-Gundersen, “A novel flexible foil vertical axis turbine for river, ocean, and
tidal applications,” Applied Energy, vol. 151, pp. 60–66, 2015. Cited page 17

[42] D. W. MacPhee and A. Beyene, “Experimental and fluid structure interaction analysis
of a morphing wind turbine rotor,” Energy, vol. 90, pp. 1055–1065, 2015. Cited page 17

[43] S. J. Hoerner, “Characterization of the fluid-structure interaction on a vertical axis

turbine with deformable blades,” 2020. 3 citations pages 17, 128, and 129

[44] S. Hoerner, S. Abbaszadeh, O. Cleynen, C. Bonamy, T. Maı̂tre, and D. Thévenin,

“Passive flow control mechanisms with bioinspired flexible blades in cross-
flow tidal turbines,” Experiments in Fluids, vol. 62, no. 5, pp. 1–14, 2021.
6 citations pages 17, 72, 74, 128, 129, and 130

[45] W. Liu and Q. Xiao, “Investigation on darrieus type straight blade vertical axis
wind turbine with flexible blade,” Ocean Engineering, vol. 110, pp. 339–356, 2015.
Cited page 18

[46] W. Mo, D. Li, X. Wang, and C. Zhong, “Aeroelastic coupling analysis of the flexible
blade of a wind turbine,” Energy, vol. 89, pp. 1001–1009, 2015. Cited page 18

[47] S. Heathcote, Z. Wang, and I. Gursul, “Effect of spanwise flexibility on flapping

wing propulsion,” Journal of Fluids and Structures, vol. 24, no. 2, pp. 183–199, 2008.
Cited page 18

[48] H.-J. Bungartz and M. Schäfer, Fluid-structure interaction: modelling, simulation, op-
timisation, vol. 53. Springer Science & Business Media, 2006. Cited page 18

[49] T. Belytschko, “Fluid-structure interaction,” Computers & Structures, vol. 12, no. 4,
pp. 459–469, 1980. Cited page 18

[50] S. Tschisgale, B. Löhrer, R. Meller, and J. Fröhlich, “Large eddy simulation of the
fluid–structure interaction in an abstracted aquatic canopy consisting of flexible blades,”
Journal of Fluid Mechanics, vol. 916, 2021. 3 citations pages 18, 19, and 23

160
 Bibliography

[51] T. Engels, D. Kolomenskiy, K. Schneider, M. Farge, F.-O. Lehmann, and J. Sesterhenn,

“Video: Bumblebee flight in turbulence: high resolution numerical simulations,” in 70th
Annual Meeting of the APS Division of Fluid Dynamics, American Physical Society,
2017. Cited page 19

[52] K. Takizawa, S. Wright, C. Moorman, and T. E. Tezduyar, “Fluid–structure interaction

modeling of parachute clusters,” International Journal for Numerical Methods in Fluids,
vol. 65, no. 1-3, pp. 286–307, 2011. Cited page 19

[53] B. Šekutkovski, A. Grbović, I. Todić, and A. Pejčev, “A partitioned solution approach

for the fluid–structure interaction of thin-walled structures and high-reynolds number
flows using rans and hybrid rans–les turbulence models,” Aerospace Science and Tech-
nology, vol. 113, p. 106629, 2021. Cited page 19

[54] A. Seghir, A. Tahakourt, and G. Bonnet, “Coupling fem and symmetric bem for dy-
namic interaction of dam–reservoir systems,” Engineering Analysis with Boundary El-
ements, vol. 33, no. 10, pp. 1201–1210, 2009. Cited page 19

[55] J. Sigüenza, D. Pott, S. Mendez, S. J. Sonntag, T. A. Kaufmann, U. Steinseifer, and

F. Nicoud, “Fluid-structure interaction of a pulsatile flow with an aortic valve model: a
combined experimental and numerical study,” International journal for numerical meth-
ods in biomedical engineering, vol. 34, no. 4, p. e2945, 2018. 2 citations pages 19 and 20

[56] J. Degroote, “Partitioned simulation of fluid-structure interaction,” Archives of

computational methods in engineering, vol. 20, no. 3, pp. 185–238, 2013.
3 citations pages 19, 21, and 143

[57] E. De Langre, Fluides et solides. Editions Ecole Polytechnique, 2001. Cited page 20

[58] P. Causin, J.-F. Gerbeau, and F. Nobile, “Added-mass effect in the design of partitioned
algorithms for fluid–structure problems,” Computer methods in applied mechanics and
engineering, vol. 194, no. 42-44, pp. 4506–4527, 2005. 3 citations pages 20, 21, and 97

[59] C. Förster, W. A. Wall, and E. Ramm, “Artificial added mass instabilities in sequential
staggered coupling of nonlinear structures and incompressible viscous flows,” Computer
methods in applied mechanics and engineering, vol. 196, no. 7, pp. 1278–1293, 2007.
3 citations pages 20, 21, and 97

[60] J. Hron and S. Turek, “A monolithic fem/multigrid solver for an ale formulation of fluid-
structure interaction with applications in biomechanics,” in Fluid-structure interaction,
pp. 146–170, Springer, 2006. 3 citations pages 20, 22, and 107

[61] T. Richter, “A monolithic geometric multigrid solver for fluid-structure interactions in

ale formulation,” International journal for numerical methods in engineering, vol. 104,
no. 5, pp. 372–390, 2015. Cited page 20

161
Bibliography

[62] C. Michler, S. Hulshoff, E. Van Brummelen, and R. De Borst, “A monolithic approach

to fluid–structure interaction,” Computers & fluids, vol. 33, no. 5-6, pp. 839–848, 2004.
Cited page 20

[63] K.-U. Bletzinger, R. Wüchner, and A. Kupzok, “Algorithmic treatment of shells and
free form-membranes in fsi,” in Fluid-structure interaction, pp. 336–355, Springer, 2006.
Cited page 20

[64] F. Sotiropoulos and X. Yang, “Immersed boundary methods for simulating fluid–
structure interaction,” Progress in Aerospace Sciences, vol. 65, pp. 1–21, 2014.
Cited page 20

[65] B. E. Griffith and N. A. Patankar, “Immersed methods for fluid–structure interaction,”

Annual review of fluid mechanics, vol. 52, pp. 421–448, 2020. Cited page 20

[66] W. Kim and H. Choi, “Immersed boundary methods for fluid-structure interaction:
A review,” International Journal of Heat and Fluid Flow, vol. 75, pp. 301–309, 2019.
Cited page 20

[67] U. Küttler and W. A. Wall, “Fixed-point fluid–structure interaction solvers with

dynamic relaxation,” Computational mechanics, vol. 43, no. 1, pp. 61–72, 2008.
2 citations pages 21 and 103

[68] C. Yvin, A. Leroyer, M. Visonneau, and P. Queutey, “Added mass evaluation with
a finite-volume solver for applications in fluid–structure interaction problems solved
with co-simulation,” Journal of Fluids and Structures, vol. 81, pp. 528–546, 2018.
Cited page 21

[69] X. Wu, X. Zhang, X. Tian, X. Li, and W. Lu, “A review on fluid dynamics of flapping
foils,” Ocean Engineering, vol. 195, p. 106712, 2020. 2 citations pages 21 and 22

[70] G. Santo, M. Peeters, W. Van Paepegem, and J. Degroote, “Dynamic load and stress
analysis of a large horizontal axis wind turbine using full scale fluid-structure interaction
simulation,” Renewable energy, vol. 140, pp. 212–226, 2019. Cited page 22

[71] J.-F. Gerbeau, M. Vidrascu, and P. Frey, “Fluid–structure interaction in blood flows
on geometries based on medical imaging,” Computers & Structures, vol. 83, no. 2-3,
pp. 155–165, 2005. Cited page 22

[72] T. Nakata and H. Liu, “A fluid–structure interaction model of insect flight with flex-
ible wings,” Journal of Computational Physics, vol. 231, no. 4, pp. 1822–1847, 2012.
Cited page 22

[73] D. W. MacPhee and A. Beyene, “Fluid–structure interaction analysis of a morphing

vertical axis wind turbine,” Journal of Fluids and Structures, vol. 60, pp. 143–159, 2016.
Cited page 23

162
 Bibliography

[74] I. Marinić-Kragić, D. Vučina, and Z. Milas, “Concept of flexible vertical-axis wind

turbine with numerical simulation and shape optimization,” Energy, vol. 167, pp. 841–
852, 2019. Cited page 23

[75] M. Benaouicha, S. Guillou, A. Santa Cruz, and H. Trigui, “Fluid-structure interaction

approach for numerical investigation of a flexible hydrofoil deformations in turbulent
fluid flow,” in Pressure Vessels and Piping Conference, vol. 51654, p. V004T04A011,
American Society of Mechanical Engineers, 2018. Cited page 23

[76] S. Hosseinzadeh and K. Tabri, “Numerical investigation of hydroelastic response of a

three-dimensional deformable hydrofoil,” in HSMV 2020, pp. 77–86, IOS Press, 2020.
Cited page 23

[77] P.-O. Descoteaux and M. Olivier, “Performances of vertical-axis hydrokinetic turbines

with chordwise-flexible blades,” Journal of Fluids and Structures, vol. 102, p. 103235,
2021. Cited page 23

[78] M. Breuer, G. De Nayer, M. Münsch, T. Gallinger, and R. Wüchner, “Fluid–structure

interaction using a partitioned semi-implicit predictor–corrector coupling scheme for
the application of large-eddy simulation,” Journal of Fluids and Structures, vol. 29,
pp. 107–130, 2012. 3 citations pages 23, 101, and 112

[79] L. Zhang, Y. Guo, and W. Wang, “Large eddy simulation of turbulent flow in a true 3d
francis hydro turbine passage with dynamical fluid–structure interaction,” International
journal for numerical methods in fluids, vol. 54, no. 5, pp. 517–541, 2007. Cited page 23
[80] M. Ilie, “Fluid-structure interaction in turbulent flows; a cfd based aeroelastic algo-
rithm using les,” Applied Mathematics and Computation, vol. 342, pp. 309–321, 2019.
Cited page 23

[81] V. Moureau, P. Domingo, and L. Vervisch, “Design of a massively parallel CFD

code for complex geometries,” C.R. Mecanique, vol. 339(2/3), pp. 141–148, 2011.
3 citations pages 25, 26, and 27

[82] M. Malandain, N. Maheu, and V. Moureau, “Optimization of the deflated con-

jugate gradient algorithm for the solving of elliptic equations on massively par-
allel machines,” Journal of Computational Physics, vol. 238, pp. 32–47, 2013.
2 citations pages 26 and 30

[83] M. Kraushaar, Application of the compressible and low-Mach number approaches

to Large-Eddy Simulation of turbulent flows in aero-engines. PhD thesis, 2011.
Cited page 26

[84] V. Moureau, P. Domingo, and L. Vervisch, “From large-eddy simulation to direct nu-
merical simulation of a lean premixed swirl flame: Filtered laminar flame-pdf modeling,”
Combustion and Flame, vol. 158, no. 7, pp. 1340–1357, 2011. Cited page 26

163
Bibliography

[85] C. Chnafa, S. Mendez, and F. Nicoud, “Image-based large-eddy simulation in a realistic

left heart,” Computers & Fluids, vol. 94, pp. 173–187, 2014. Cited page 26

[86] Y. Dufresne, V. Moureau, G. Lartigue, and O. Simonin, “A massively parallel cfd/dem

approach for reactive gas-solid flows in complex geometries using unstructured meshes,”
Computers & Fluids, vol. 198, p. 104402, 2020. Cited page 26

[87] J. Dagaut, M. Negretti, G. Balarac, and C. Brun, “Linear to turbulent görtler instability
transition,” Physics of Fluids, vol. 33, no. 1, p. 014102, 2021. Cited page 26

[88] F. Doussot, G. Balarac, J. Brammer, Y. Laurant, and O. Métais, “Rans and les sim-
ulations at partial load in francis turbines: Three dimensional topology and dynamic
behaviour of inter blade vortices,” International Journal of Fluid Machinery and Sys-
tems, vol. 13, no. 1, pp. 12–20, 2020. Cited page 26

[89] M. Guilbot, Analyse et optimisation des performances de turbines à axe vertical et flux
transverse par simulations numériques. PhD thesis, Université Grenoble Alpes, 2021.
2 citations pages 26 and 36

[90] P. Bénard, A. Viré, V. Moureau, G. Lartigue, L. Beaudet, P. Deglaire, and L. Bricteux,

“Large-eddy simulation of wind turbines wakes including geometrical effects,” Comput-
ers & Fluids, vol. 173, pp. 133–139, 2018. Cited page 26

[91] “coria cfd youtube channel.” https://www.youtube.com/user/CoriaCFD/videos.

Cited page 27

[92] J. H. Williamson, “Low-storage runge-kutta schemes,” Journal of computational

physics, vol. 35, no. 1, pp. 48–56, 1980. Cited page 28

[93] A. Chorin, “Numerical solution of the navier-stokes equations,” Math. Comp., vol. 22,
pp. 745–762, 1968. 2 citations pages 28 and 30

[94] J. Donea, A. Huerta, J.-P. Ponthot, and A. Rodrı́guez-Ferran, “Arbitrary l agrangian–e

ulerian methods,” Encyclopedia of computational mechanics, 2004. Cited page 28

[95] C. Chnafa, Using image-based large-eddy simulations to investigate the intracardiac flow
and its turbulent nature. PhD thesis, Université Montpellier II-Sciences et Techniques
du Languedoc, 2014. 2 citations pages 28 and 32

[96] C. Farhat, P. Geuzaine, and C. Grandmont, “The discrete geometric conservation

law and the nonlinear stability of ale schemes for the solution of flow problems on
moving grids,” Journal of Computational Physics, vol. 174, no. 2, pp. 669–694, 2001.
Cited page 29

[97] D. Boffi and L. Gastaldi, “Stability and geometric conservation laws for ale formula-
tions,” Computer methods in applied mechanics and engineering, vol. 193, no. 42-44,
pp. 4717–4739, 2004. Cited page 29

164
 Bibliography

[98] J.-L. Guermond, P. Minev, and J. Shen, “An overview of projection methods for in-
compressible flows,” Computer methods in applied mechanics and engineering, vol. 195,
no. 44-47, pp. 6011–6045, 2006. Cited page 30

[99] T. Jaravel, Prediction of pollutants in gas turbines using large eddy simulation. PhD
thesis, 2016. Cited page 31

[100] V. Moureau, Simulation aux grandes échelles de l’aérodynamique interne des moteurs
à piston. PhD thesis, Châtenay-Malabry, Ecole centrale de Paris, 2004. Cited page 32

[101] M. Germano, U. Piomelli, P. Moin, and W. H. Cabot, “A dynamic subgrid-scale eddy

viscosity model,” Physics of Fluids A: Fluid Dynamics, vol. 3, no. 7, pp. 1760–1765,
1991. 2 citations pages 32 and 122

[102] M. Germano, “Turbulence: the filtering approach,” Journal of Fluid Mechanics,

vol. 238, pp. 325–336, 1992. Cited page 32

[103] J. T. Batina, “Unsteady euler airfoil solutions using unstructured dynamic meshes,”
AIAA journal, vol. 28, no. 8, pp. 1381–1388, 1990. 2 citations pages 32 and 88

[104] C. Dobrzynski and P. Frey, Anisotropic Delaunay Mesh Adaptation for Unsteady Sim-
ulations, pp. 177–194. Springer Berlin Heidelberg, 2008. Cited page 32

[105] C. Dapogny, C. Dobrzynski, and P. Frey, “Three-dimensional adaptive domain remesh-

ing, implicit domain meshing, and applications to free and moving boundary problems,”
Journal of computational physics, vol. 262, pp. 358–378, 2014. Cited page 32

[106] P. Benard, G. Balarac, V. Moureau, C. Dobrzynski, G. Lartigue, and Y. D’Angelo,

“Mesh adaptation for large-eddy simulations in complex geometries,” International
Journal for Numerical Methods in Fluids, vol. 81, no. 12, pp. 719–740, 2015. fld.4204.
Cited page 32

[107] S. Pertant, M. Bernard, G. Ghigliotti, and G. Balarac, “A finite-volume method

for simulating contact lines on unstructured meshes in a conservative level-
set framework,” Journal of Computational Physics, vol. 444, p. 110582, 2021.
2 citations pages 32 and 89

[108] R. Janodet, G. Vaudor, G. Lartigue, P. Benard, V. Moureau, and R. Mercier, “An

unstructured conservative level-set algorithm coupled with dynamic mesh adaptation
for the computation of liquid-gas flows,” in 29th European Conference on Liquid Atom-
ization and Spray Systems (ILASS Europe), 2019. 2 citations pages 32 and 89

[109] A. Pushkarev, P. Benard, G. Lartigue, V. Moureau, and G. Balarac, “Numerical ap-

proach for simulation of moving bodies by using the dynamic mesh adaptation method
within ale technique,” in ECCOMAS MSF 2017, 2017. Cited page 32

165
Bibliography

[110] N. Guillaud, G. Balarac, E. Goncalves, and J. Zanette, “Large eddy simulations on

vertical axis hydrokinetic turbines-power coefficient analysis for various solidities,” Re-
newable Energy, vol. 147, pp. 473–486, 2020. Cited page 36

[111] M. Kaiho, M. Ikegawa, and C. Kato, “Parallel overlapping scheme for viscous incom-
pressible flows,” International journal for numerical methods in fluids, vol. 24, no. 12,
pp. 1341–1352, 1997. Cited page 36

[112] C. Kato, M. Kaiho, and A. Manabe, “An overset finite-element large-eddy simula-
tion method with applications to turbomachinery and aeroacoustics,” J. Appl. Mech.,
vol. 70, no. 1, pp. 32–43, 2003. Cited page 36

[113] A. Hrennikoff, “Solution of problems of elasticity by the framework method,” 1941.

Cited page 44

[114] R. Courant, K. Friedrichs, and H. Lewy, “On the partial difference equations of mathe-
matical physics,” IBM journal of Research and Development, vol. 11, no. 2, pp. 215–234,
1967. Cited page 44

[115] J. H. Argyris and S. Kelsey, Energy theorems and structural analysis, vol. 60. Springer,
1960. Cited page 44

[116] O. C. Zienkiewicz, R. L. Taylor, P. Nithiarasu, and J. Zhu, The finite element method,
vol. 3. McGraw-hill London, 1977. Cited page 44

[117] O. Zienkiewicz, “Origins, milestones and directions of the finite element method—a
personal view,” Archives of computational methods in engineering, vol. 2, no. 1, pp. 1–
48, 1995. Cited page 44

[118] O. C. Zienkiewicz, R. L. Taylor, and J. Z. Zhu, The finite element method: its basis
and fundamentals. Elsevier, 2005. 3 citations pages 44, 45, and 64

[119] O. C. Zienkiewicz and R. L. Taylor, The finite element method for solid and structural
mechanics. Elsevier, 2005. 6 citations pages 44, 54, 57, 66, 71, and 143

[120] M. K. Thompson and J. M. Thompson, ANSYS mechanical APDL for finite element
analysis. Butterworth-Heinemann, 2017. Cited page 44

[121] Y. Debard, “Calcul des structures par la methode des elements finis.” http://iut.
univ-lemans.fr/ydlogi/index.html. 4 citations pages 44, 46, 63, and 69

[122] “efunda, engineering fundamentals.” https://www.efunda.com/formulae/solid_

mechanics/mat_mechanics/hooke_plane_stress.cfm. Cited page 49

[123] A. K. Chopra, “Dynamics of structures: Theory and applications to earthquake engi-

neering, prentice hall,” Inc., Upper Saddle River, NJ, 1995. Cited page 58

166
 Bibliography

[124] R. W. Clough and J. Penzien, “Dynamics of structures mcgraw-hill,” Inc Editor, 1993.
Cited page 58

[125] J. Chung and G. M. Hulbert, “A time integration algorithm for structural dy-
namics with improved numerical dissipation: the generalized-α method,” 1993.
3 citations pages 60, 61, and 62

[126] M. Abramowitz and I. A. Stegun, “Handbook of mathematical functions dover publi-

cations,” New York, vol. 361, 1965. Cited page 65

[127] “code-aster documentation.” https://www.code-aster.org/V2/doc/default/fr/

index.php?man=R3. 2 citations pages 65 and 71

[128] J. R. Shewchuk et al., “An introduction to the conjugate gradient method without the
agonizing pain,” 1994. Cited page 70

[129] D. Code Aster, “Généralités sur le gradient conjugué: Gcpc aster et utilisation de
petsc.” 2 citations pages 70 and 134

[130] S. Balay, S. Abhyankar, M. F. Adams, S. Benson, J. Brown, P. Brune, K. Buschelman,

E. M. Constantinescu, L. Dalcin, A. Dener, V. Eijkhout, W. D. Gropp, V. Hapla,
T. Isaac, P. Jolivet, D. Karpeev, D. Kaushik, M. G. Knepley, F. Kong, S. Kruger,
D. A. May, L. C. McInnes, R. T. Mills, L. Mitchell, T. Munson, J. E. Roman, K. Rupp,
P. Sanan, J. Sarich, B. F. Smith, S. Zampini, H. Zhang, H. Zhang, and J. Zhang,
“PETSc Web page,” 2021. https://petsc.org/. 2 citations pages 71 and 143

[131] S. Balay, S. Abhyankar, M. F. Adams, S. Benson, J. Brown, P. Brune, K. Buschelman,

E. Constantinescu, L. Dalcin, A. Dener, V. Eijkhout, W. D. Gropp, V. Hapla, T. Isaac,
P. Jolivet, D. Karpeev, D. Kaushik, M. G. Knepley, F. Kong, S. Kruger, D. A. May,
L. C. McInnes, R. T. Mills, L. Mitchell, T. Munson, J. E. Roman, K. Rupp, P. Sanan,
J. Sarich, B. F. Smith, S. Zampini, H. Zhang, H. Zhang, and J. Zhang, “PETSc/TAO
users manual,” Tech. Rep. ANL-21/39 - Revision 3.16, Argonne National Laboratory,
2021. 2 citations pages 71 and 143

[132] S. Balay, W. D. Gropp, L. C. McInnes, and B. F. Smith, “Efficient management of

parallelism in object oriented numerical software libraries,” in Modern Software Tools
in Scientific Computing (E. Arge, A. M. Bruaset, and H. P. Langtangen, eds.), pp. 163–
202, Birkhauser Press, 1997. 2 citations pages 71 and 143

[133] J. F. Thompson, Z. U. Warsi, and C. W. Mastin, Numerical grid generation: founda-

tions and applications. Elsevier North-Holland, Inc., 1985. Cited page 88

[134] S. P. Spekreijse, “Elliptic grid generation based on laplace equations and algebraic
transformations,” Journal of Computational Physics, vol. 118, no. 1, pp. 38–61, 1995.
Cited page 88

167
Bibliography

[135] S. P. Spekreijse, B. B. Prananta, and J. C. Kok, “A simple, robust and fast algorithm
to compute deformations of multi-block structured grids,” 2002. Cited page 88

[136] M. Potsdam and G. Guruswamy, “A parallel multiblock mesh movement scheme for
complex aeroelastic applications,” in 39th Aerospace Sciences Meeting and Exhibit,
p. 716, 2001. Cited page 88

[137] C. Farhat, C. Degand, B. Koobus, and M. Lesoinne, “Torsional springs for two-
dimensional dynamic unstructured fluid meshes,” Computer methods in applied me-
chanics and engineering, vol. 163, no. 1-4, pp. 231–245, 1998. Cited page 88

[138] C. Sheng and C. B. Allen, “Efficient mesh deformation using radial basis func-
tions on unstructured meshes,” AIAA journal, vol. 51, no. 3, pp. 707–720, 2013.
2 citations pages 88 and 143

[139] M. Bernard, G. Lartigue, G. Balarac, V. Moureau, and G. Puigt, “A framework to

perform high-order deconvolution for finite-volume method on simplicial meshes,” In-
ternational Journal for Numerical Methods in Fluids, vol. 92, no. 11, pp. 1551–1583,
2020. Cited page 88

[140] G. De Nayer, Interaction Fluide-Structure pour les corps élancés. PhD thesis, Ecole
Centrale de Nantes (ECN), 2008. Cited page 88

[141] E. Lefrançois, “A simple mesh deformation technique for fluid–structure interaction

based on a submesh approach,” International Journal for Numerical Methods in Engi-
neering, vol. 75, no. 9, pp. 1085–1101, 2008. 2 citations pages 88 and 143

[142] F. Duchaine, T. Morel, and A. Piacentini, “On a first use of cwipi at cerfacs,” Contract
report TRCMGC-11-3. CERFACS, vol. 213, 2011. Cited page 94

[143] E. Quémerais, “Coupling with interpolation parallel interface,” ONERA web site, 2016.
Cited page 94

[144] M. A. Fernández, J.-F. Gerbeau, and C. Grandmont, “A projection semi-implicit

scheme for the coupling of an elastic structure with an incompressible fluid,” Inter-
national Journal for Numerical Methods in Engineering, vol. 69, no. 4, pp. 794–821,
2007. Cited page 103

[145] S. Turek and J. Hron, “Proposal for numerical benchmarking of fluid-

structure interaction between an elastic object and laminar incompress-
ible flow,” in Fluid-structure interaction, pp. 371–385, Springer, 2006.
7 citations pages 107, 109, 110, 111, 112, 114, and 116

[146] M. Schäfer, S. Turek, F. Durst, E. Krause, and R. Rannacher, “Benchmark compu-

tations of laminar flow around a cylinder,” in Flow simulation with high-performance
computers II, pp. 547–566, Springer, 1996. Cited page 107

168
 Bibliography

[147] A. Kalmbach and M. Breuer, “Experimental piv/v3v measurements of vortex-induced

fluid–structure interaction in turbulent flow—a new benchmark fsi-pfs-2a,” Journal of
Fluids and Structures, vol. 42, pp. 369–387, 2013. 2 citations pages 119 and 120

[148] G. De Nayer, A. Kalmbach, M. Breuer, S. Sicklinger, and R. Wüchner, “Flow past a

cylinder with a flexible splitter plate: A complementary experimental–numerical inves-
tigation and a new fsi test case (fsi-pfs-1a),” Computers & Fluids, vol. 99, pp. 18–43,
2014. Cited page 121

[149] G. De Nayer and M. Breuer, “Numerical fsi investigation based on les: Flow
past a cylinder with a flexible splitter plate involving large deformations (fsi-pfs-
2a),” International Journal of Heat and Fluid Flow, vol. 50, pp. 300–315, 2014.
4 citations pages 121, 122, 123, and 126

[150] A. Kondratyuk, “Investigation of the very large eddy simulation model in the context
of fluid-structure interaction,” 2017. Cited page 126

[151] “Manufacturing and experiment of pitching flexible foil.” https://wasd.urz.

uni-magdeburg.de/hoerner/thesis/videos/fabrication_profil_souple.mp4.
Cited page 129

[152] A. E. Bogaers, S. Kok, B. D. Reddy, and T. Franz, “Quasi-newton methods for implicit
black-box fsi coupling,” Computer Methods in Applied Mechanics and Engineering,
vol. 279, pp. 113–132, 2014. Cited page 143

[153] D. Blom, “Efficient numerical methods for partitioned fluid-structure interaction simu-
lations,” 2017. Cited page 143

[154] A. Naseri, O. Lehmkuhl, I. Gonzalez, E. Bartrons, C. D. Pérez-Segarra, and A. Oliva,

“A semi-implicit coupling technique for fluid–structure interaction problems with strong
added-mass effect,” Journal of Fluids and Structures, vol. 80, pp. 94–112, 2018.
Cited page 143

[155] H.-J. Bungartz, F. Lindner, B. Gatzhammer, M. Mehl, K. Scheufele, A. Shukaev, and

B. Uekermann, “precice–a fully parallel library for multi-physics surface coupling,” Com-
puters & Fluids, vol. 141, pp. 250–258, 2016. Cited page 143

169
Résumé
De récentes études expérimentales tendent à montrer que l’utilisation de profils
flexibles pour des turbines à axe vertical améliore leur rendement et augmente
leur durée de vie. La simulation numérique de telles expériences permettrait une
meilleure compréhension du phénomène d’Interaction Fluide-Structure (IFS) en
jeu. Cependant, cela nécessite un solveur IFS capable de prédire le décrochage
dynamique avec précision, tout en calculant la déformation d’un solide à
géométrie complexe. Le but de la thèse est donc de développer un solveur
capable de reproduire fidèlement des cas de profils flexibles à haut nombre de
Reynolds. Le manuscrit présente ainsi le développement d’un solveur employant
une approche de Simulations aux Grandes Echelles (SGE) pour prédire la
dynamique du fluide et des éléments finis solides 3D pour la dynamique du
solide, le tout sur des maillages non structurés. Les développements ont été
réalisés au sein de la librairie YALES2, qui a initialement été conçue pour la
mécanique des fluides. Un solveur structure a donc été développé lors de ce
travail pour prédire les déplacements du solide. De plus, une méthode originale
du type pseudo-solide a été proposée pour le calcul du mouvement de maillage
fluide. Les solveurs fluide et solide sont fortement couplés avec un schéma
partitionné permettant de reproduire des cas où le solide et le fluide ont des
densités proches. Le solveur IFS résultant de ce travail est donc en mesure de
reproduire une large variété de cas complexes, et peut aussi employer une
méthode d’adaptation dynamique de maillage pour pouvoir prendre en compte de
grands déplacements solides. Une fois les méthodes numériques expliquées en
détail, le solveur IFS est validé. Un cas numérique de référence 2D laminaire est
d'abord reproduit avec succès. Ensuite, une validation 3D est réalisée à partir
d’une expérience d’une plaque flexible accrochée derrière un cylindre disposé
dans un écoulement à haut Reynolds. Enfin, une expérience d’un profil flexible
en oscillation dans un canal est reproduite. Malgré des améliorations possibles
concernant le temps de calcul, cela confirme l’utilisation potentielle du solveur
IFS pour l’usage souhaité.

Mots-clés : Interaction Fluide-Structure, Simulations aux Grandes Echelles,

Méthode des éléments finis, Maillage non structuré, Adaptation Dynamique de
Maillage, Mouvement de maillage, Pseudo-solide, Schéma partitionné

Abstract
Recent experimental studies tend to show that the use of chordwise flexible
blades for vertical axis turbines enhances their efficiency and increases their
lifetime. Numerical simulation of such experiments would provide a better
understanding of the Fluid-Structure Interaction (FSI) phenomenon involved.
However, this requires a FSI solver capable of accurately predicting the stall
dynamics, while computing the deformation of a solid with complex geometry.
Thereby, the goal of the thesis is to develop a solver capable of accurately
reproducing cases of flexible foils at a high Reynolds number. The manuscript
thus presents the development of a solver using the Large-Eddy Simulations (LES)
approach to predict fluid dynamics and 3D solid finite elements for solid
dynamics, all on unstructured meshes. The development work has been carried
out within the YALES2 library, which was initially designed for fluid mechanics.
Consequently, a Structural Mechanics Solver (SMS) has also been developed.
Moreover, an original method based on the pseudo-solid approach has been
proposed for the computation of the fluid mesh movement. Fluid and solid
solvers are strongly coupled with a partitioned scheme, providing the opportunity
to study dynamics between some fluid and solid of close densities. The FSI solver
resulting from this work is therefore able to reproduce a wide variety of complex
cases, and can also use a dynamic mesh adaptation method to take into account
large solid displacements. Following the review of the numerical methods, the FSI
solver can then be validated. To begin, a 2D numerical reference case with
laminar flow is successfully reproduced. After that, a 3D validation is performed
against an experiment where a flexible plate is clamped behind a cylinder
immersed in a high Reynolds number flow. Finally, an experiment of a pitching
flexible foil in a channel is reproduced. Despite possible improvements regarding
computational time reduction, this process confirms the potential of the FSI solver
for its intended use.

Keywords : Fluid-Structure Interaction, Large-Eddy Simulation, Chordwise flexible

blade, Finite Element Method, Unstructured grid, Dynamic Mesh Adaptation,
Mesh movement, Pseudo-solid, Partitioned scheme