Structural Analysis Using Openfoam

Fluid-Structure Interaction (FSI) case study of a
cantilever using OpenFOAM and DEAL.II with

application to VIV
Johan Lorentzon
June 17, 2009
ii
.
Thesis for the degree of Master Science in Technical Mathematics.
ISSN 0282-1990
ISRN LUTMDN/TMHP09/5186SE
c Johan Lorentzon, Juni 2009
Division of Fluid Mechanics
Department of Energy Sciencies
Lunds Institute of Technology
Box 118
S-211 00 LUND
Sweden
Abstract
This Master Thesis in Technical Mathematics at LTH, directed towards simulation and com-
putation, has treated the subject of uid-structure interaction (FSI) for incompressible ow
with small vibrations. The open source packages DEAL.II and OpenFOAM have been used
to create a coupling between a nite element formulation for structure and a nite volume
formulation for uid ( gas or liquid ). A staggered solution algorithm for FSI has been im-
plemented in C++ using Aitkens relaxation method together with a Reduced-Order-Model
(ROM). The solution algorithm has been validated by using an application consisting of a
cantilever immersed in a steady ow transversal to its axial direction. Also, vortex-induced
vibrations (VIV) were calculated as a function of the ow velocity and successfully com-
pared to empirical data. The study has demonstrated the usefulness of artical damping to
solve boundary condition problems in incompressible ow. The method presented is general
and has a strong potential in technical applications where the structure is subjected to a
surrounding uid.
iii
iv
Acknowledgement
To my supervisors, Johan Revstedt, Robert Zoltan Szasz, Per Erik Austrell: My thanks for
being supportive in my work. To my friend Mathias Haage, I am truly grateful for your en-
couragement as consultant in the eld of C++. The teachers Per Lidstrom, Matti Ristinmaa,
Mathias Wallin, Claus F uhrer, Niels Saabye Ottosen and Jan Gustavsson my thanks for en-
during my endless questions in the eld of continuum mechanics and mathematics. I would
also like to thank Gunnar Sparr, for the gift of the precious compendium Kontinuerliga
system from year 1983, the very same version i lost on the train last year, a compendium
that truly inspired me in my work. Finally, my thanks to my mentor Andrej Sadlev, for the
years of fruitful discussions and meetings in the topic applied mathematics.
December 2008.
v
vi
Nomenclature
Abbreviations
BC Boundary Condition
CFD Computational Fluid Dynamics
Co Courant number
cv control volume
FD Finite Dierence
FEM Finite Element Method
FIV Flow-Induced Vibrations
FSI Fluid-Structure Interaction
FVM Finite Volume Method
INS Incompressible Navier Stokes Equation
KC Keulegan Carpenter Number
PDE Partial Dierential Equation
Re Reynolds number
ROM Reduced Order Model
St Strouahls number
VIV Vortex-Induced Vibration
Mathematical notation
Tensor notation is adopted to distinguish tensor A
i
from vector notation A, and this implies
the convention [12] of summation index,
A
i
B
i
=
k
A
k
B
k
, (1)
and the dierential convention for tensor,
A
ji,i
=
k
A
jk
. (2)
The dierence between vector and tensor is the applied rule upon coordinate transformation
for a tensor. Although, ambiguity in interpretation of formulas appears, a vector in this thesis
is a representation of a tensor. The origin of these rules comes from the that since physical
laws are invariant to coordinate transformations and when expressed in tensor notation the
physical constants also shall be invariant.
vii
viii
Symbols
A m Amplitude of marker point
A
t
m
2
Area of the tip of cantilever
b Nm
3
Force per volume
C kgs
1
Damping matrix
D m Width of cantilever
E Pa Youngs Modulus
f s
1
VIV frequency
f
S
N Force acting on structure
f
F
N Force acting on uid
f
n
s
1
Natural frequency of structure
I m
4
Moment of inertia
K Nm
1
Stiness matrix
L m Height of cantilever
l m Characteristic length
M kg Mass matrix
n Surface normal
p Pa Pressure
P Pa Pressure dierence
q m Displacement eld of structure
t Pa Traction vector
T
r
s Reference time period
T
c
s Calculated time period
T
v
s Vacuum time period
U ms
1
Flow velocity
U ms
1
Magnitude of velocity at inlet
U
r
ms
1
Relative velocity
v ms
1
Velocity eld of structure
V
r
ms
1
Reduced velocity
Symmetric gradient operator
m
2
s
1
Kinematic viscosity
kgm
1
s
1
Viscosity
f
kgm
3
Density of ow
s
kgm
3
Density of structure
s
1
Forced frequency
Poisson number
Pa Stress
Contents
1 Introduction 3
1.1 The Scenario . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.2 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.3 Outline of the FSI algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.4 The motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.5 Objective of the thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.6 The material used for this thesis . . . . . . . . . . . . . . . . . . . . . . . . . 8
2 Theory 9
2.1 The Continuum hypothesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.2 Numerical approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.3 Governing equations for the incompressible ow . . . . . . . . . . . . . . . . . 10
2.4 The Finite Volume Method applied to the INS . . . . . . . . . . . . . . . . . 11
2.5 Quasi-steady approximation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.6 Weak form for static system of the solid state . . . . . . . . . . . . . . . . . . 14
2.7 The Finite Element Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.8 Governing equation in the state space formalism . . . . . . . . . . . . . . . . 15
2.9 ALE description and imposed BC . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.10 The Fluid-Structure Interaction . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.10.1 Monolithic FSI Problem . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.10.2 Staggered FSI algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.11 Acceleration of the convergence . . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.11.1 Aitkens relaxation method . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.11.2 The ROM method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
3 Implementation 21
3.1 The compilation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
3.2 Main program ICOFSI . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
3.3 OpenFOAM Package . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
3.4 DEAL.II package . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
3.5 InterGridMapping . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
3.5.1 Mapping functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
3.5.2 Transfer step . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
4 The case study 27
4.1 Beam theory of the cantilever . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
4.2 The steady state response for a cantilever. . . . . . . . . . . . . . . . . . . . . 27
4.3 Boundary condition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
ix
CONTENTS 1
4.4 Single cantilever study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
4.5 Multi-block cantilever study . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
4.6 The cantilever mesh . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
4.7 The uid domain mesh . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
5 Validation procedure 33
5.1 Solid state solver . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
5.1.1 Test Case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
5.1.2 Verication of class StaticElastic . . . . . . . . . . . . . . . . . . . . . 34
5.1.3 Verication of class QuasiElasticity . . . . . . . . . . . . . . . . . . . . 34
5.1.4 Verication of the class DynamicElastic . . . . . . . . . . . . . . . . . 35
5.2 icoDyMFoam . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
5.3 setBoundaryIndicator and InterGridMapping . . . . . . . . . . . . . . . . . . 37
5.4 The FSI algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
5.4.1 Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
5.4.2 A note on measurement of the period . . . . . . . . . . . . . . . . . . 38
5.4.3 Result . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
5.5 Error analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
5.5.1 icoDyMFoam . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
5.5.2 DynamicElastic . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
5.5.3 Estimated Error in the period . . . . . . . . . . . . . . . . . . . . . . . 48
5.6 General comments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
5.6.1 The theta method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
5.6.2 Explicit/Implicit FSI . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
5.6.3 Aitkens relaxation and the ROM . . . . . . . . . . . . . . . . . . . . . 53
5.7 Multiple block case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
5.8 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
6 An application to VIV 57
6.1 Problem description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
6.2 Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
6.3 Result . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
7 Discussion 67
8 Future work 71
2 CONTENTS
Chapter 1
Introduction
This chapter describes the scenario and the incentive for this study. The background of uid-
structure interaction and its technical applications is presented with a brief description of
the vortex-induced vibration. This is followed by an outline of the staggered algorithm and
ends with a survey of the research eld on this subject which denes the goals of the thesis
and the material used.
1.1 The Scenario
A cantilever is placed in a domain of a velocity driven uid. The stress acting upon the
structure induces a deformation of the structure to which the uid responds. This mutual
inuence referred to as uid-structure interaction (FSI), is known to cause several interesting
phenomena. Among such is vortex-induced vibration, where the forced movement of a uid
around the structure gives upon point of release from the structure, an angular momentum
manifested as a vortex in the uid. Further, due to no-slip condition between uid and struc-
ture, the structure will experience an extra mass while moving and hence decrease the natural
frequency.
For this simple case there is an analytical formula for the steady FSI but already when in-
cluding the eect of conning walls the task becomes quickly overwhelming while using pencil
and paper. At this point the computer is unhanding us this problem. This requires a theory
to model the FSI and a numerical procedure such that when the problem is implemented,
a computer simulation can gives us an approximate answer to our problem. Although every
step is simple from a mathematical point of view, the number of steps can be an exerting
task for even a skilled scientist. An ecient engineer takes the work of others and reshape
their tools for his own purpose.
The computer can be instructed in several languages, each language fullling their purpose.
For many years Fortran has been considered as the key language for a computational engi-
neer due to its ecient implemented libraries and compilers adapted to high-speed machines
dierent from PC. However, as the architecture of the computer changes the dierence be-
comes more vague and favor Object-Oriented (OO) languages such as C++/Java/Python
which minimizes the development time and spread the usability to far more users due to
its simplicity in implementation [70, 74]. Without further dwelling into details, the actual
question that matters is how to organize the work and how fast it can be implemented, not
how much faster the code will become or easier it is to code on a detailed level.
3
4 CHAPTER 1. INTRODUCTION
1.2 Background
The Fluid-Structure Interaction (FSI) appears as a physical phenomenon in engineering [13]
such as static load, drag, and Flow-Induced-Vibrations (FIV). FIV is further classied into
utter, galloping and vortex-induced vibrations (VIV) [32]. The static load on the structure
originates mainly from the total pressure dierence between front and wake side of the struc-
ture. The dynamic pressure can due to instability in the structure, where the damping is less
than the energy transfered by the vortex-induced wake, create a resonance with the structure
with a negative damping known as utter.
Alternatively it will be manifested as VIV, where the pressure dierence created by the vor-
tices causes an elastic response in the structure with smaller amplitudes inducing fatigue
cycles on the structure. A related issue is bueting, which mainly concern the turbulence
random excitation on the structure. Galloping falls between the regime of utter and VIV.
The classication of FIV is however not precise where overlaps between the areas are a com-
mon thread to misunderstanding. However, the VIV gives the answers to a large number of
questions involved in this phenomenon of FSI,
Estimation of the structural response to a given ow.
The mechanism behind the coupling between uid and structure.
Estimation of the fatigue in structure.
Design structure to minimize the self-induced vibrations.
The questions are related to each other but focused on dierent aspects: prediction, under-
standing and prevention. Typical examples of applications in which utter is considered are
aircraft wings, tall buildings, rotating blades and long-span bridges [51, 47, 48]. The VIV
concerns bridges [59], chimneys [10], heat exchanger tubes [45], power lines and undersea
constructions such as sea cables/risers [18, 9] or biomechanical applications as blood vessels
[53].
A simplied case of self-induced vibration is forced vibrational study where the oscillation
to given frequency of the structure in transversal direction to given ow direction creates
a coupling force with the surrounding uid. Several studies have showed a lock-in frequency
with the uids frequency of vortex shedding, f. In a review by Williamson [58], focused on
cylinder/cantilever, experiments and models explains the phenomenon thoroughly. In this
review, the reduced velocity V
r
(
U
l
) versus
f
show a linear dependency except for a region

around = f. The cylinder/cantilever or n motion behind blu body in enclosed com-
partments are classical applications allowing self-induced vibrations [20, 38, 31]. Further, V
r
versus the reduced amplitude in transversal direction to the ow for self-induced free vibra-
tion
A
l
show a coupling between the form in release of vortex and the magnitude of reduced
amplitude response. The results implicate that a steady state motion turns into a unsteady
motion for a certain frequency irrelevant to the size of the body. At this point it is evi-
dent of the importance in understanding the synchronization of a wake mechanism with the
movement of the structure.
From this knowledge a passive control [6, 32, 10] of FIV can be developed to prevent dam-
aging eects of short term [24] and long term [45] eects on the structure. However, the
characteristic of the vortex shedding depend on the structure [31, 26, 24]. The simplest mode
of vortex shedding is the one studied in the cylinder case, where the alternate release of
vortex is observed as a street of pairs of vortices with opposite spin, denoted as 2S mode.
This is a direct consequence when taking the blu body eect into consideration. This eect
1.3. OUTLINE OF THE FSI ALGORITHM 5
forces the uid to take a curved path which at point of release have a curl with same the spin
as the ow around the body, a Vortex is born, see Figure 1.1. The next type of mode is 2P
appearing when the lock-in frequency cause a phase shift in which the FSI coupling creates
two pair of vortices within each periodic cycle hence four vortices per cycle. The design to
Figure 1.1: A principal sketch of a 2S vortex formation process [42].
prevent these phenomena is however not always the same and to analyze a problem as an
engineer requires skills in empirical estimating the factors of inuence. An infamous example
of such failure in estimating the eect of dynamics of the ow is the Tacoma bridge which
also illustrates the misconception between cause and consequence. The utter is the cause to
the collapse and there was no lock-in between VIV and the natural frequency of the structure
[7, 24]. Still, the pre-loading of the structure could have its source from VIV.
Although the question of dierence seem semantic, since both phenomena are self-induced
vibrations, it is noted that this is an event that has been debated over a half century. In
the process to develop a passive control of VIV, a Direct Numerical Simulation to support
experimental investigation is desirable. This minimize the cost involved with reinforcement
and maintenance.
1.3 Outline of the FSI algorithm
Fluid mechanics is characterized to be a highly chaotic dynamical system and solid mechanics
by its self preserving and elastic behavior. The governing equations expressed in PDE and
their characteristics could therefore be dierent, changing from hyperbolic/parabolic to el-
liptic. By implementing the whole uid-structure interaction in one solver, also known as the
monolith approach [25, 8], requires more computational eort and may lead to less accurate
simulation due to the limitation in the modeling. This loss in accuracy has its origin mainly
in that the numerical method depends on the characteristic of PDE where a system that
changes characteristic during simulation requiring projection techniques solved by introduc-
ing intermediate variables which split the moment equation into sub-equations or by using
semi-discretization techniques and assume ow characteristic such as incompressible ow.
In both cases stability criteria arises and error is introduced, although in dierent senses,
numerical and modeling error respectively. The second aspect in loss of accuracy is the lim-
itation of the computational resources, which requires a cut-o in precision. However, the
traditional way to solve equations of this type is by using the staggered technique [37, 33].
The staggered solver exchanges the interaction term in an iterative scheme by solving each
domain separately. The particular advantage of this method is the possibility to use a black
box solver for each governing equation, and with a small eort merge them into a functional
FSI algorithm. This implicates the possibility to use dierent solution techniques that enables
a more accurate modeling.
The Staggered FSI algorithm
Time Loop
Staggered Loop
Solve Fluid State
Transfer Traction to Solid State Solver
Solve Solid State
Exit Staggered Loop if change of deformation < tolerance
Transfer Deformation to Fluid State Solver
End Staggered Loop
End Time Loop
The accuracy is thus aected by restricting the governing equation in PDE form to a given
characteristic, or rene the constitutive modeling and/or using the staggered techniques which
formulates the linear matrix equations as symmetric, while the monolith approach leads to
asymmetric ones. The problematic part is that the staggered technique for the incompressible
case suer severely by stability issues. This occurs whenever the change in variables under
each iteration induce a change such that when inserted into the equation, the residual is of
comparable size as the coupling term. Improper BC conditions also aect this due to the
domain of dependency for incompressible ow.
1.4 The motivation
ANSYS+CFX [64] have both monolith and staggered techniques available, ABAQUS and
FLUENT [63] have been used in several studies, COMSOLE MULTIPHYSICS [62] has a
module for both techniques. The commercial products give low insight in solution techniques
and their verication is built upon cases which can be downloaded and tested against other
codes using open source, but any manipulation of the algorithm itself is limited although
some has console window option allowing meta-coding. The usage of the staggered technique
with ABAQUS and FLUENT gives elaborate choice of tool box and solution techniques for
each domain separately. However, the COMSOLE MULTIPHYSIC gives the full control of
the solution techniques but compared to ANSYS+CFX 11.0 it is less ecient to a factor of
three using the version 3.4 [62, 66]. Other commercial products of interest are ADINA [4]
and LS-DYNA [67].
For university applications, it is however important apart from the license cost benet, both
from an educational and a research point of view to be able to modify the algorithm. There
are several open source codes for solving FSI. The solvers implemented by individual research
teams often fall within the categories of staggered, monolith or approximative modeling. Fur-
ther, one can distinguish between stochastical approach and deterministic.
One solver of interest is the icofsi solver from OpenFOAM developed by Tukovic and Jasak
[53]. This is a solver based on the Finite Volume Method (FVM) using Green-Cauchy strain
measure that allow arbitrary deformations with application to a cantilever in a rectangular
1.4. THE MOTIVATION 7
domain. This article use the Aitkens under-relaxation technique to accelerate the convergence
of the fully coupled FSI is presented and the dierence is demonstrated between explicit and
implicit FSI. The monolith method suer from the ill-conditioned system matrix and pre-
cision loss due to computational limitation, but there is an approach to solve this [25] by
using a uid Pressure Poisson Equation for INS. The validation of this algorithm is based
upon calculating the added mass and study VIV on a cylinder immersed in a velocity driven
domain.
To improve the convergence of the monolithic approach, an article by M. Razzaq et al pre-
sented a Newton-Raphson scheme using a Finite-Dierence scheme (FD) to evaluate the
Jacobian [43] and a similar approach is performed with a staggered method [37]. A step
further is the reduced-order model (ROM), where the governing equations are replaced by a
local estimate of the PDE by using a FD scheme or a residual method. This is used either as
a preconditioner in the sub-cycle or replaces the solver to save computational eort [54, 16].
In the staggered solution technique, allowing parallel execution is an essential step in order
to maintain a simulation within a reasonable time and to a limited cost solve larger struc-
tures. Several articles have been presented on this subject but [15] highlights the instability
region where one falsely concludes the utter point due to the explicit staggered scheme. This
article also further points out the limitation in fully coupled in precision due to small time
step required and emphasis the need of adaptive time stepping to achieve convergence. In the
eort to stabilize the fully coupled staggered algorithm the Aitkens
2
relaxation method
was implemented with successful result [33, 53].
A topic not included in this thesis is the turbulence, which add the bueting eect in FSI,
this is vital since the random excitations on a structure have a signicant eect on durability
and in aeroelasticity [50]. In all discussed articles above use the Arbitrary Lagrangian-Euler
description essential for FSI and for an excellent review by J.Donea et al on the topic see
[11]. The list of individual contributions in FSI is vast, for further see review by Dowell [13]
and references therein or the Wikipedia link [72].
1.5 Objective of the thesis
This thesis has the goal to implement an implicit FSI algorithm by coupling a CFD and a
FEM solver. Then validate the FSI solver by reproducing the common phenomena observed
such as VIV and frequency/amplitude shifts, then analyze and optimizing its performance.
The code should have generic input allowing dierent types of mesh to be used. Another goal
is to generalize the algorithm to multiple structures that enable dierentiation in modeling.
The code is restricted to elastic modeling with small strain of the structure with limited
adaptive meshing. The programming language is C++ with the goal to create an interface
such that it can easily be developed to include non-linear basis elements, non-linear strain,
under-integrations techniques, accelerated techniques in convergence and adaptive meshing
in each time step. The eort is merely replace/adding a class in the code.
1.6 The material used for this thesis
The open source packages used in this project are OpenFOAM [27, 65] and DEAL.II[61]. This
allows to run the FSI solver on separate machines/threads. The original application is taken
from an article by Sampaoi [50] but using the staggered algorithm and parameter settings
presented by Tukovic and Jasak [53]. An article by M. Dreier et al [55] about the natural
frequency shift of a cantilever immersed in a uid is used in order to verify the solver since
suitable experimental data for the given case and free vibrational study were not be found. In
order to increase the performance, the reduced-order model (ROM) was implemented using
the article by Vierendeels as guidance [54]. Further, to accelerate/stabilize the sub-cycle loop
in the staggered technique, and to locate the quasi-static equilibrium point between the
solvers by the x-point iteration, the Aitkens relaxation method was chosen. This method
is described by K uttler and Wall [33]. An important resource is the forum sites associated
with the homepages for each open source package. The developers and several research teams
report their issue and share knowledge and experience which is very valuable for a beginner
on the subject.
Chapter 2
Theory
The problem can be divided into three domains. First, the physical domain, the container
of our phenomena of interest, formulated in terms of physical laws and variables. Secondly,
the model domain, the set of approximations describing the physical domain to be governed
by equations in PDE. Finally, the numerical domain containing the procedure in how to
implement and solve the PDE in a computer environment. This chapter rst discusses the
limits of the controle volume size where constitutive relation is relevant, then distinguish
the dierence between FVM and FEM. The governing equations describing each domain
is derived and the boundary conditions and explicit formulas used in this thesis will be
discussed.
2.1 The Continuum hypothesis
The concept used in this thesis is the deterministic behavior of physical domain allowing
functional relation between variables. This is when the molecular uctuations becomes in-
nitesimal and statistical mechanics enters the realm of continuum mechanics. The essential
part is the macroscopic model description where constitutive equations interrelate to the
variable representing the physical quantities of interest. The model description of a physical
domain is often expressed in terms of governing equations, a set of state equations whose
state variables are related by PDE in strong form or as state functions. The validity of the
continuum hypothesis can often be estimated by evaluating the size of a given control volume
V when the statistical uctuations of a given parameter becomes negligible. Although this
limit is a gray-scale, less than 1 10
6
m in cell diameter, the statistical description becomes
signicant, and at 1 10
9
m quantum eects cannot be neglected.
2.2 Numerical approach
Given the model expressed in PDE, the problem can be solved in two fundamentally dif-
ferent ways, either one discretize the operators and then solves the parameter relation, or
one discretizes the solution space and uses the exact PDE to dene a residual which is to
be minimized. The dierent approaches can then of course be combined. The mesh is the
topology together with the grid that denes the solution space to a given problem. Discretize
the PDE consistently with ensured stability, then in the limit of innitesimal mesh the Finite
Dierence scheme (FD) will approach the exact solution. Another approach is to integrate
the strong form over the control volume, and even further partial integrate into a weak form
9
10 CHAPTER 2. THEORY
of the system, where a solution to a strong PDE is the solution to a weak form but not
necessary the other way around. It is under this category the Finite Volume Method (FVM)
enters using a conservative form of the PDE. The other approach is to multiply the strong
form with a test function, then from this its weak form is formulated. This leads to a bilinear
form which denes a linear system, whose solution in the limit of innitesimal mesh goes to
the exact solution of the original strong PDE form. A general applied method is the Galerkin
method where the test functions are chosen to be the discretized solution space, the FEM
method.
2.3 Governing equations for the incompressible ow
The continuity equation states that for a given property for a given innitesimal volume
element V with the ux F, the sum of the unsteady change of the and the change in ux
is equal to the source S
t
+ F = S
. (2.1)
By applying the conservation of mass S
= 0, setting =
f
and F = U, where U is the
velocity through given V of the uid, then the continuity equation for mass becomes
f
+ (
f
U) = 0. (2.2)
Dene the material time derivative as
D
Dt
=
t
+U . (2.3)
By identifying the material time derivative in Eqn (2.2) the expression takes the following
form
D
f
Dt
=
f
( U). (2.4)
Incompressible ow is characterized by a zero material derivative which implies the condition
U = 0. (2.5)
Applying the conservation of momentum to the Eqn (2.1), dening =
f
U
i
and the ux
F = U, using the conservation of mass, the Newtons second law can be identied with LHS
and therefore the source term becomes the force f
i
,
f
DU
i
Dt
= S
f
i
. (2.6)
Using continuum mechanics on f
i
, dening b
i
as the body force acting per volume, and stress
tensor
ij
acting upon the boundary, the equation of motion is given by,
f
i
=
ij,j
+ b
i
. (2.7)
At this point, this applies to both uid and solid mechanics and the equation formed by
Eqn (2.6) and Eqn (2.7) is a part of the Incompressible Navier Stokes equations (INS). The
dierence lies upon the constitutive relation dened for the shear stress part of the stress
tensor,
ij
,
ij
= p
ij
+
ij
. (2.8)
For a Newtonian uid using symmetric gradient =
1
2
(+
T
),
ij
= ( U)
ij
+(
j
U
i
+
i
U
j
) = ( U)
ij
+ 2
ij
(U), (2.9)
2.4. THE FINITE VOLUME METHOD APPLIED TO THE INS 11
where
2
3
in order to satisfy entropy condition [52], where the Stokes condition as-
sert equality. Using the incompressibility condition Eqn (2.5) on Eqn (2.9), then Eqn (2.6)
becomes (
2
= )
f
DU
Dt
= p +
2
U+b. (2.10)
The conservation of angular momentum has to be taken into consideration, this gives the
condition of symmetric stress tensor (
ij
=
ji
) as Eqn (2.9) provides. Denote l as the
characteristic length, a the reference speed, dene then U
=
U
a
, t
=
ta
l
, x
i
=
xi
l
, p
=
p
f
a
2
,
f
and
= l,the following dimensionless governing equation is obtained,
= 0, (2.11)
DU
Dt
+
1
Re
2
U
+b
, (2.12)
where Re =
al
is the Reynolds number. Note, OpenFOAM CFD package have not adopted
the convention to scale velocity, only the pressure. The equations originating from the energy
conservation is omitted in this thesis.
2.4 The Finite Volume Method applied to the INS
The center value Finite Volume Method (FVM) uses the conservative integral forms of the
governing equations of the uid characterized by the following [17, 28],
spatial derivatives over volume are converted into integrals over surfaces in terms of
the ux F.
time derivative is semi-discretized.
the grid points dene the faces and the discretization reservoir points are at the center
of the control volume (cv).
the uxes are interpolated at each step.
the integrals are evaluated by the use of the mean value theorem.
In applying the FVM for INS, the strong form of the PDE for the INS is integrated over a
control volume (cv). This control volume is polyhedral, arbitrary as long it is convex and has
planar surface elements. The physical reservoir point P is at the centroid x
P
of the cv, the
interior domain/volume V
P
. Each given face k of the cv has its centroid point x
k
and the
neighbor N.
By Taylor expansion of the tensor to the rst order around the centroid of the cv, the
integration error becomes to the second order accurate with the size of the cv,
=
P
+ (x x
P
) (x
P
) +O(x
2
). (2.13)
The face values
k
of the of the x
P
is estimated by the Central Dierence (CD) approxi-
mation,
()
k

x
k
x
P
x
P
x
N
(x
P
) + (1
x
k
x
N
x
P
x
N
)(x
N
). (2.14)
and dening S(k) to be the face surface normal with the magnitude of the area of the face
k, then S ()
k
is approximated by,
S ()
k
S

N

P
x
N
x
P
. (2.15)
However, CD cause unbounded solutions whenever convective term is the dominating term
[28]. In non-orthogonal meshes, a correction term appears that is omitted here for simplicity
for all ux and face values. In similar fashion for second order temporal discretization the
error is proportional to (t)
2
.
The INS (2.5) and (2.6) become with
i
=
f
U
i
and denitions from previous section,
_
Vp
UdV = 0, (2.16)
_
t+t
t
_

t
_
Vp
i
dV +
_
Vp
(
i
U)dV +
_
Vp
i
pdV
_
Vp
i
dV
_
dt = 0. (2.17)
Dene the F(k) = S ()
k
as the numerical mass ux. Using the divergence theorem on the
convective term, taking
i
to be the value in cv P at x
P
this term then becomes
_
Vp
(U
i
)dV =
_
Vp
U
i
dS
k
F(U
i
)
k
. (2.18)
Taking the transient part of Eqn (2.17) to be semi-discretized by using implicit Euler scheme
where superscript n is an index to time t
n
. Use Eqn (2.18) to discretize the convective term,
and by using divergence theorem on the diusion term together with Eqn (2.15), then the
LHS of Eqn (2.17) becomes
_
t+t
t
_
(
(
n
i
)
P
(
n1
i
)
P
t
)V
P
+
k
S
i
(p)
k
+
k
F(U
n
i
)
k
k
S (
n
i
)
k
_
dt. (2.19)
By assuming no change in the control volume V
P
, the update scheme then becomes,
(
n
i
)
P
= (
n1
i
)
P
+
t
V
P
_
k
S
i
(p)
k
+
k
F(U
n
i
)
k
k
S (
n
i
)
k
_
. (2.20)
By this scheme in Eqn (2.20), a stability condition is introduced, the so called Courant
number Co = |
U
f
(xP xN)
t
| < 1. In the implicit Euler scheme the ux and gradient of the
ux is calculated on the current time scale n in accordance to Eqn (2.14) and Eqn (2.15).
The INS is non-linear in the convective term. For each cv Eqn (2.20) denes a linear relation,
A
cv
cv
= a
P
P
+
N
a
N
N
= R
cv
. (2.21)
The coecients a
P
and a
N
are calculated from values of current time state
n
i
and the R
contain the previous time state. A loop over all cv in a given domain using a topological
index mapping creates a sparse matrix relation,
A[
n
]
n
cv
A
cv
cv
=
cv
R
cv
R[
n1
]. (2.22)
2.5. QUASI-STEADY APPROXIMATION 13
In the case of a structured grid, where the cv become cubes, this scheme coincides with
second order central dierence FD scheme with semi-discretized implicit Euler scheme in
time. Note that A depends on
i
, this non-linearity must be solved, and is termed pressure-
velocity dependency. The procedure to solve Eqn (2.20) implies three such linear equations
to be solved, where the ux interrelates all three. In order to dene the solution procedure
in icoDyMFoam [28, 65], the solver in OpenFOAM used for the uid domain in this thesis,
some denitions are needed. First the H[U] operator,
H[U] =
N
a
N
U
N
+
U
o
P
t
, (2.23)
where U
o
is rest term from Eqn (2.20) minus the pressure term. Then the explicit velocity
correction becomes,
U
P
=
H[U]
a
p
1
a
P
(p)
P

H[U]
a
p
k
1
a
P
S(p)
k
. (2.24)
By using Eqn (2.24) into the incompressible ow condition Eqn (2.16), the pressure equation
is then obtained, a Laplacian solved by using Eqn (2.14) and Eqn (2.15),
k
S
_
(
1
a
P
)
k
(p)
k
_
=
k
S (
H[U]
a
P
)
k
. (2.25)
The BC, physical as well numerical (uxes and interpolated face values) is implemented by
explicit correction in the Eqn (2.22) and projection techniques, the details of these steps are
omitted.
The pressure-velocity dependency is solved by Rhie-Chow interpolation in the following
scheme, known as the PISO algorithm,
Solve the U, by solving the momentum predictor step using ux and pressure from
previous time step in Eqn (2.22).
Calculate H[U] using Eqn (2.23).
Update the ux F using Eqn (2.24).
Solve the pressure equation, Eqn (2.25).
Update U using Eqn (2.24).
Repeat from second step until tolerance is reached in U.
This assumes the coecients in H[U] in terms of previous time steps, hence changing the
dependence of A given in Eqn (2.22) from
n
to
n1
! The correction for this creates the
PIMPLE algorithm which is simply a repetition of the PISO until convergence is achieved.
The issue by this approximation is not further pursued in this thesis but details are given in
the PhD thesis of Jasak [28].
2.5 Quasi-steady approximation
The quasi-steady approximation assumes that the structural response is in equilibrium with
the uid instant, using small strain assumption, with =
1
2
(+
T
) and then the following
constitutive relation is used,

ij
= D
ijkl
kl
( q), (2.26)
where D
ijkl
is the tensor expression for linear elasticity and q the displacement eld. Using
implicit Euler and multiply on both sides the t this becomes at given time step n,
n
ij
=
n1
ij
+D
ijkl
kl
(q
n
), (2.27)
inserting this into equation of motion Eqn (2.7) assuming b
i
= 0,
j
D
ijkl
kl
(q
n
) = f
n
i
+
n1
ij,j
. (2.28)
Rewrite this by using the divergence theorem starting with the LHS, taking the scalar product
with test function over interior = (t
n1
). Dene n
j
the normal vector of the surface
(t
n1
) using that
ij
is symmetric,
(
j
D
ijkl
kl
(q
n
)),
i
)
= (n
j
D
ijkl
kl
(q
n
),
i
)
(D
ijkl
kl
(q
n
),
ij
())
. (2.29)
The RHS,
(f
n
i
,
i
)
+ (
n1
ij,j
,
i
)
= (f
n
i
,
i
)
(
n1
ij
,
ij
())
+ (n
j
n1
ij
,
i
)
. (2.30)
Using the BC, n
j
D
ijkl
kl
(q
n
) = f
i
(t
n
) f
i
(t
n1
) = f
i
f
n1
i
, zero displacement on
D
,
(D
ijkl
kl
(q
n
),
ij
())
= (f
i
(t
n
),
i
)
N
(
n1
ij
,
ij
())
. (2.31)
where
D
is the Dirichlet boundary,
N
is the Neumann boundary. This is the bilinear
form a(, ) for the structural part in Euler description for FEM. The complicated part is the
need for calculate rotation matrix update of the second order stress tensor for each cell. For
innitesmal incremental step q
n
, the method holds for arbitrary deformations q
n
.
2.6 Weak form for static system of the solid state
The equation of motion, Eqn (2.7) applying the Newtons second law on RHS or use Eqn
(2.6) with U=0, in strong form can in tensor notation be formulated as,
ij,j
+b
i
=
s
q
i
. (2.32)
Assume linear elastic constitutive law, using symmetric gradient =
1
2
( +
T
) and dis-
placement eld q,
ij
= D
ijkl
kl
(q). (2.33)
The strong form then becomes
s
q
i
j
D
ijkl
kl
(q) = b
i
. (2.34)
The static system is obtained by putting
s
q
i
to 0. By using isotropic approximation over
the mesh domain with Lame parameter (,), the D
ijkl
can be dened as,
D
ijkl
=
ij
kl
+(
ik
lj
+
il
kj
). (2.35)
The partial dierential operator from Eqn (2.34) of elliptic characteristic, acting on solution
q, can be described as a bilinear operator a(, ) within a nite sub domain V
n
. This is obtained
by multiplying the strong form with a trial function v, integrate over a control volume then
using the divergence theorem,
a(q, v) = (
i
q
i
,
j
v
j
)
+ (
i
q
j
,
i
v
j
)
+ (
i
q
j
,
j
v
i
)
. (2.36)
This relation is used in the FEM approximation as described in the next section.
2.7. THE FINITE ELEMENT METHOD 15
2.7 The Finite Element Method
In the Galerkin approximation, given a nite dimensional subspace V
n
of the solution space
V to the PDE satisfying the BC, then nd the q V
n
such as,
a(q, v) = (f, v), v V
n
. (2.37)
Expanding the displacement eld q by a set of test functions,
q =
k
q
k
q
k
N
q
q. (2.38)
and dene the residual = a(q, v) (f, v), where the v is a trial function in V
n
. The a(, )
depends upon the constitutive modelling. In this thesis Eqn (2.36) is used. The integration
is performed over the cell N, A
N
ij
= a(v
i
, v
j
)
N
, and the summation of all contributing cells
is denoted as the assemble of the stiness matrix,
K
ij
a(v
i
, v
j
) =
N
A
N
ij
. (2.39)
Each cell integration is approximated by Gauss quadrature, and isoparametric mapping.
Given any function g, dene the computational domain V
0
, the current domain by V and
the Gauss quadrature points
i
. This implies a Jacobian for the deformation of the cell from
computational domain to the current domain, and weight factor
i
related to the choice of
quadrature points
i
,
_
V
g(x)dV =
_
V0
g()detJd
i
g(
i
)detJ(
i
)
i
. (2.40)
In adaptive meshing, one has to ensure the condition of conforming elements [60] and in order
to avoid linear dependence in a(, ), a constraint equation is formed to lock the freedom of
the nodes which are on the edges of each element, the sub-cells, referred as to hanging nodes
in the DEAL.II package [61].
The solution procedure is to minimize above residual = 0 which expressed in the basis of
V
n
becomes
Kq = f, (2.41)
where q is the displacement vector in the space spanned by the trial functions v, K is the
stiness matrix and f the force vector. The setting of the trial functions as the displacement
eld functions based upon a given mesh gives the Finite Element Method (FEM). The BC of
Dirichlet type is met by condensing the linear matrix, and the Neumann condition appears
as traction elements in the f vector.
2.8 Governing equation in the state space formalism
The strong form of the PDE for the Eqn (2.7) in state space,
q
i
v
i
= 0, (2.42)
t
v
i
+C
il
v
l
j
D
ijkl
kl
(q) = f
i
. (2.43)
In state space, the structural dynamics for the undamped/[damped] system takes the
following form,
q v = 0, (2.44)
M v + [Cv] +Kq = f. (2.45)
where [60],
K =
_
s
(N
q
)
T
D(N
q
)d
s
, (2.46)
M =
_
s
s
N
T
q
N
q
d
s
. (2.47)
Here K is same as previous section but slightly reformulated due to Eqn (2.36). Assuming
the structural damping matrix C,
C = M+K, (2.48)
known as Rayleigh damping, useful when the natural frequency is known. The damping is
introduced by adding the [Cv] term to the equation in version 0.3 and thereby its distinction.
By semi-discretization of the time variable in Eqn (2.44) and Eqn (2.45), a =
a
n
a
n1
t
, and
introduce blending by using the theta method, a = a
n
+ (1 )a
n1
,
q
n
= q
n1
+ t v, (2.49)
[M+ tC]v
n
= [M(1 )tC]v
n1
tK q + t
f. (2.50)
There are two ways of solving this, updating rst v and then q gives the simplest equation
to solve,
[M+tC+
2
t
2
K]v
n
= [M(1)tC(1)t
2
K]v
n1
tKq
n1
+t
f. (2.51)
However, by multiplying the Eqn (2.44) by [M + tC] and replace the v
n
in Eqn (2.44)
using Eqn (2.45) gives the second variant:
[M+tC+t
2
2
K]q
n
= [M+tCt
2
(1 )K]q
n1
+tMv
n1
+t
2
f, (2.52)
[M+ tC]v
n
= [M(1 )tC]v
n1
tK q + t
f. (2.53)
These are the two forms of Roths equation in matrix form used in this project. In the
stationary case, when velocity goes to zero the two variants of semi-discretized equations
are equivalent to the expected Kq = f. The implemented code using this set of equations
is DynamicElastic for dynamic case, StaticElastic for the static system. Final note, semi-
discretization is used to separate the time domain from the spatial domain in the discretiza-
tion which enables the use of an ecient adaptive meshing, another procedure is to use the
method of lines, which implies a decoupling of the PDE into a set of ODE requiring one mesh
for each equation, hence adaptive meshing becomes more expensive.
2.9. ALE DESCRIPTION AND IMPOSED BC 17
2.9 ALE description and imposed BC
The FEM in this study is expressed in terms of the material description, with displacement
eld as a function of original grid, while the FVM of the uid has the solution with respect
the current grid, the Euler description, with the pressure and velocity expressed in current
grid. The no-slip condition is applied to the moving wall, which gives the additional condi-
tion that the uid cell adjacent at the boundary have the same velocity as the grid points
of solid sharing the same face to given uid cell (U
r
= 0). This implies that the grid points
of the uid move coherently with the points of the solid on the coupling boundary. The
pressure of the adjacent cell has imposed zero gradient condition (p, n) = 0 which implies
pressure static condition, hence the isotropic stress tensor of uid and pressure has the same
direction but with opposite sign. The uid also induce a drag originating from the deviatoric
stress tensor which at given grid point is the applied traction with the included pressure.
The strategy to solve this problem, is by applying the traction vector on the solid state and
to the uid impose the displacement on the boundary to the uid and re-mesh the grid and
re-interpolate U and p. This leads to the Dirichlet-Neumann partition condition in Arbitrary
Lagrange-Euler description for the staggered solution of the FSI problem, for more details see
[33, 11]. The critical aspect, in combining the descriptions, using the benets of both is the
complicating feature of the additional eort to re-mesh the domain described by the Euler
description. The re-mesh procedure in OpenFOAM is based upon the pseudo-solid displace-
ment approximation, which has its direct consequence that it assumes the displacement to
be described by a potential, that is a solution to a Laplace equation.
2.10 The Fluid-Structure Interaction
A monolithic formulation of the FSI in FEM serves the purpose to give a closed expression
of the co-moved mass factor and present the frequency shift for a cantilever with forced
vibration analysis on VIV response. A detailed description of the staggered algorithm using
Dirichlet-Neumann partition technique is presented.
2.10.1 Monolithic FSI Problem
The Laplacian form in strong form of PDE in Eqn (2.17) allows the FEM approximation. For
the purpose of study the added mass eect inducing frequency shift and amplitude response
FSI in this section, assume the uid under FEM approximation.
The mass eect of the uid is to be derived. In the following, for simplicity of the equations,
assume incompressible inviscid ow and ignore the convective term and the gravitational
eect for the uid. By further using incompressible condition to eliminate the velocity then
the Helmholtz governing equation for INS of the pressure is obtained,
2
p = 0. (2.54)
Assert FEM solution to the problem,
p N
p
p, (2.55)
q N
q
q, (2.56)
and dropping the tilde henceforth. For the structure, assume no damping and neglect the
viscous eect. Then the following system equations for the uid-structure coupling is obtained
[60, 46],
Hp +Q
T
q +f
F
= 0, (2.57)
M q +Kq Qp +f
S
= 0, (2.58)
where f
S
and f
F
contains the force boundary terms. The new terms in the dynamical system,
is the coupling term Qp involving traction vector t = pn where n is the uid surface normal
and the pressure term H,
Qp =
_
N
t
q
td, (2.59)
H =
_
f
(N
p
)
T
(N
p
)d
f
. (2.60)
The matrices form the same bilinear functions described in the previous sections. By substi-
tuting the Eqn (2.57) into Eqn (2.58) gives the following structured equation to be solved,
(M+QH
1
Q
T
) q +Kq +QH
1
f
F
+f
S
= 0. (2.61)
The QH
1
Q
T
is known as added mass matrix, or co-moved mass, this is the frequency shift
causing the structure to move slower. For the cantilever, there is an analytical solution [55]
to this added mass matrix, the time period for the natural frequency,
T
r
= (1 +
1
4
L
D
s
)
1
2
T
v
, (2.62)
where T
v
is the natural time period in vacuum. The VIV can be analyzed by forced vibration
analysis. Starting from Eqn (2.58) neglecting the f
S
term, and inserting the forced vibration
assumption q = q
0
e
it
and p = p
0
e
it
,the forced uid structured equations of the VIV
then becomes,
q
0
= (K
2
(M+QH
1
Q
T
))
1
Qp
0
. (2.63)
This is the structural response on forced vibration for the steady incompressible inviscid ow.
As the forced vibration approaches to the natural frequency of the structure, resonance
should be observed. At this point, VIV is obtained as setting = f and it is plausible to
assume Qp
0
p, pressure dierence across the transversal direction sides of the cantilever.
2.10.2 Staggered FSI algorithm
Given the problem dened in ALE description, also known as Dirichlet-Neumann partition
[33]. Dene F the solver for the uid and S the solver for the solid. The sub-cycle in the
staggered algorithm to solve the problem can then be formulated as,
Solve F(q
n+1
, t
n+1
, v
n+1
) = f
F
.
calculate the uids traction vector T on the coupled boundary applied as Neumann
condition on force load f
S
(t
n+1
) . (Neumann transfer step)
Solve S(q
n+1
) = f
S
(t
n+1
).
From the displacement eld q, from solid grid, move the boundary of the uid. (
Dirichlet transfer step)
Move the uid mesh accordingly q = 0. This is the pseudo-solid approximation ap-
plied to the uid mesh.
repeat above procedure until convergence achieved such that predicted displacement
eld is consistent with given pressure dierence, that is

f = F(S
1
(F(f
n+1
))) = f
n+1
or q = S
1
(F(q
n+1
)) = q
n+1
.
The procedure of this kind is also known as a strongly coupled solver [13]. The solver has the
option to use a weakly coupled method, that is only one iteration.
2.11. ACCELERATION OF THE CONVERGENCE 19
2.11 Acceleration of the convergence
This section concern with two basic schemes to accelerate the convergence of the xed-point
iteration in the strongly coupled staggered FSI algorithm.
2.11.1 Aitkens relaxation method
In order to stabilize the convergence in the sub-cycle, a relaxation procedure on the Dirichlet
transfer step was introduced, a generalization of the secant method, where nding the root
of scalar equation f(x
) = 0 is achieved by the algorithm [1],

(x
n+1
x
n
) =
(x
n
x
n1
)
(f(x
n
) f(x
n1
))
f(x
n
). (2.64)
Assume the existence of xed point x
= g(x
), then f(x) = g(x) x becomes the residual

(g(x
n
) = x
n+1
),
x =
_
(x
n
x
n1
)
((x
n+1
x
n
) (x
n
x
n1
))
_
(x
n+1
x
n
). (2.65)
The new solution is then x
n+1
= x
n
+ x, the factor within brackets is the measure of a
scaling of the step. Dene x
n+1
= g(x
n
) and generalize the factor to n-dimensional space,
introducing scalar dots instead of multiplication, and dierence in denominator by norm,
n
=
(x
n
x
n1
)
(x
n+1
x
n
) (x
n
x
n1
)
2
(x
n+1
x
n
). (2.66)
The underlying interpretation is that for an algorithm accurate to the rst order, the dif-
ferential residual gives an estimation of the second order correction. The Aitkens relaxation
then becomes,
x =
n+1
(x
n+1
x
n
) + (1
n+1
)(x
n
x
n1
). (2.67)
The under relaxation factor is then modied for the next step as
n+1
=
n
n
. Then applied
to FSI the
n+1
is limited to an upper bound to achieve an under-relaxed step.
2.11.2 The ROM method
The previous iterative solution can be used both as precondition during the FSI cycle and
to approximate the time steps by interpolated solutions of previous steps. This is known as
reduced order model (ROM) [16]. Using the previous solutions, the local derivatives in the
Taylor expansion are expressed by FD schemes or Gauss-Newton method, giving an estimate
of the next step solution. In this study, since this was not of high priority, the simplest choice
was implemented,
q
n+1
= 3(q
n
q
n1
) +q
n2
. (2.68)
The DynamicElastic have a routine implemented such that it calculates exactly 3 points
then extrapolates a given number of consequent steps, note however within the FSI cycle
this formula is used in all steps. A n-ROM method is referred to as a ROM using n previous
steps.
Chapter 3
Implementation
The functionality of the DEAL.II [61] version 6.1 and OpenFOAM version 1.4.1 [65] packages
will be explained.
3.1 The compilation
Here follows a short instruction how to compile the packages. First download the open source
packages and install them accordingly to the individual help les. Then in constructing the
case, starting with the template discussed in the user guide for the OpenFOAM with wmake
as constructing script. Add under same catalogue the DEAL.II code encapsulated under a
namespace in order to avoid name conicts. Open the option le under the Make catalogue
and add the necessary options/ags to be used in order to compile the DEAL.II. However,
this is done in practice by trial and error using wmake iterative and including all missing
links into the option le from the output of the failed compilation log le.
3.2 Main program ICOFSI
The uid solver used in this thesis and described in the theory section is icoDyMFoam and
the main program icoDyMFoam.C, located under the application section of the OpenFOAM
1.4.1 release, was taken as a template for the staggered algorithm described in section 2.7
and the PISO algorithm was enclosed as header icoDyMFoam.H. The DEAL.II code was
enclosed into a separate class under the namespace DynamicElastic and a namespace share
contained the common data as mapping index between the domains and the traction vector.
The brackets [] in the pseudo code shows to which class in respective step is assigned to in
the FSI algorithm. The evaluation of the traction vector was placed in the main function
in version 0.9. An issue with the gnu g++ compiler prevented to use the template class as
dened by standard C++. The solution to this is to dene the class every time it is used.
Earlier version than 0.7 used only the isotropic part of the tensor, the pressure axial on faces
of cv.
21
22 CHAPTER 3. IMPLEMENTATION
ICOFSI PSEUDO CODE
Load Fluid Mesh
Loop over multiple solid objects
Load Grid [ dealiiGrid class ]
Define boundary condition [dealiiSetBoundaryIndicator]
Load parameters [dealiiParameterHandler]
Mapping between fluid and solid [InterGridMapping]
Loop over Time
Adjustable Timestep
ROM interpolation step [InterGridMapping]
Repeat until convergence FSI cycle
Dirichlet Transfer Step [AutoMesh]
Solve Fluid [icoDyMFoam]
Calculate Traction Vector [Section 3.3]
Refine Mapping functions [InterGridMapping]
Neumann Transfer Step [InterGridMapping]
Solve Solid [DynamicElastic]
Aitkens relaxation [DynamicElastic]
Update
Probes [dealiiProbes]
Fluid mesh [AutoMesh]
A short description of the class and their functionality,
DynamicElastic, FEM solver. The algorithm from section 2.6 with acceleration step
from 2.8.
icoDyMFoam, FVM solver, OpenFOAM INS solver in ALE description, described in
section 2.2.
InterGridMapping, handles the mapping between uid and solid face centers, based
upon proximity distance classier, section 3.5.2.
AutoMesh, assign displacement eld vector from interpolated face center on solid onto
the center of uid.
dealiiParameter, the input handler to DynamicElastic.
dealiiGrid, read the grid from msh format or create hard coded mesh for the Dynam-
icElastic.
dealiiSetBoundaryIndicator, creates the Neumann and Dirichlet BC on the solid mesh.
3.3. OPENFOAM PACKAGE 23
3.3 OpenFOAM Package
The OpenFOAM package developed in object-oriented C++ is a CFD library tool to handle
numerical solution of PDE using FVM. The moment predictor step Eqn (2.22) is evaluated
using class fvV ectorMatrix which discretize and distributes the operators specied within
curled brackets, the code extracted from include le UEq.H in icoDyMFoam,
fvVectorMatrix UEqn
{
fvm::ddt(U)
+fvm::div(phi,U)
-fvm::laplacian(nu,U)
}
solve(UEqn==-fvc::grad(p));
where fvm class refers to the solution method, the Finite Volume Method, fvm :: ddt
the discretization scheme for local derivative in time, fvm :: div as divergence operator
discretized as described by Eqn (2.18), fvm :: laplacian solves the viscosity term using
( ) using Eqn(2.15) and Eqn (2.18). The fvc is referred as the Finite Volume Calculus
class, containing the interpolation schemes, in this case (CD). The discretization scheme is
determined by input le. The following illustrates the transcription between theory and the
OpenFOAM,
A
i
U
i
= R
i
solve(UEqn == -fvc::grad(p))
a
P
AU=UEqn().A()
U =
1
aP
H[U] U=UEqn().H()/AU
= S H[U] phi=vc::interpolate(U)& mesh.Sf()
(
1
aP
p) = () fvm::laplacian(1.0/AU, p)==fvc::div(phi)
U = U
1
aP
p U-=fvc::grad(p)/AU
The package contains data types which enables a compact description of the physical
variables and overloading techniques to allow a direct interpretation between a mathematical
formula in PDE and coding context. This is the benet of Finite Dierence/FVM using collo-
cated grids. The snippet below taken from ICOFSI main code illustrates the implementation
of the second order tensor and interpolation technique, the original source of this code was
extracted from liftDrag.H taken from the OpenFOAM community [65]. This is
ij
part
which acts as the traction to be transfered between the grids, see Eqn (2.8). The viscous
traction follows by the second term in Eqn (2.9).
Foam::tmp<Foam::Field<Foam::Vector<double> > >
ViscousTraction =-nu.value()*
U.boundaryField()[Patch].snGrad();
PressureTraction = p.boundaryField()[OpenFOAM_Patch]*
Density.value*
mesh.Sf().boundaryField()[Patch]/
mesh.magSf().boundaryField()[Patch];
Traction = PressureTraction+ViscousTraction;
The snGrad() function is the member function of U calculating the magnitude of the
viscous part of the traction, the returned value ViscousTraction is the vector eld U with
magnitude of viscous force. Regarding the second term PressureTraction, the Sf() member
function is the normal of the face center, magSf() is the length, the face area. Hence the
normal unit vector of the face cell is computed times the pressure. The overloading tech-
nique of C++ allows the compact formulation of tensor algebra. The original article [27] for
OpenFOAM is recommended for explaining the interface and the structure of the package.
3.4 DEAL.II package
The DEAL.II package, is a C++ implemented FEM library for solving Dierential Algebraic
Equation (DAE) systems, using adaptive mesh and ecient linear algebraic solver libraries as
PetScWrapper and Trilinos both allowing MPI. In a comparative study [36] its strength were
also its weakness. Because of its high modularity, there was no automatic tools, providing
standard solvers, as OpenFOAM provides. Following the notation from the theory section,
the following list creates the connection between the code and theory,
a(, ) system matrix
f system rhs
v
i
(
q
) fe values.shape value(i,q)
q
fe values.quadrature point(q)
|detJ(
i
)|
i
fe values.JxW(q)
The documentation of the code is extensive, with articles on the subject aiding in detail
in the construction of the stiness matrix, the correction for the non-conforming elements
in the adaptive meshing and the use of multi-threading. The tutorial section on this forum
explains the use of the package in executable programs named as step-X, where X = 1, .., 33.
In similar fashion as for OpenFOAM, overloading techniques/member function is used to
give an easy transcription between code and the weak form formulated in the theory section.
The template used for DynamicElastic and StaticElastic is step-18, step-8, step-11 step-10
and step-23.
The DEAL.II functionality is exemplied by the following snippet code of the evaluation of
the surface integral of traction vector dened in section 3.3,
_
s
T v(x)
i
dx
q
1
s
T(
q
)v(
q
)
i
|detJ(
q
)|
i
. (3.1)
The rst loop extends over the faces of a given cell. The boundary indicator() ==1, ensures
the surface face is chosen, the reinit(...) member function initializes the values of all test func-
tions related to given cell. The inner loop is over the quadrature points, which approximates
the integral.
3.4. DEAL.II PACKAGE 25
double invrho=1/rho;
const unsigned int n_face_q = face_quadrature_formula.size();
comp_i = fe.system_to_component_index(i).first;
for (unsigned int face=0;
face<GeometryInfo<dim>::faces_per_cell;++face)
if ( cell->face(face)->boundary_indicator() == 1 ) {
index=cell->face_index(face);
fe_face_values.reinit(cell,face);
for ( unsigned int q=0;q<n_face_q;++q) {
cell_rhs(i)+=
invrho*
shared::db[n_object]->get_traction(comp_i,index)*
fe_face_values.shape_value(i,q)*fe_face_values.JxW(q);
}
}
}
The comp i link the connection to the degree of freedom to which coordinate axis in the
vector index. From this snippet,
system_matrix.copy_from(mass_matrix);
system_matrix.add(theta*time_step,C_matrix);
hanging_node_constraints.condense (system_matrix);
stiffness_matrix.vmult (system_rhs, solution_u_cur);
system_rhs *= -theta * time_step;
stiffness_matrix.vmult (tmp, solution_u_prev);
system_rhs.add (-time_step * (1-theta), tmp);
mass_matrix.vmult (tmp, solution_v_prev);
system_rhs += tmp;
C_matrix.vmult (tmp, solution_v_prev);
system_rhs.add (-(1-theta)*time_step, tmp);
system_rhs += forcing_terms;
the clear dierence between OpenFOAM and DEAL.II is illustrated. The snippet corre-
sponds to the construction of the system matrix and system vector in equation of motion
for the velocity update, Eqn (2.50). The left hand side of the equation is formed by member
function copy from from mass matrix followed by the vector/matrix addition add function.
The vmult is the member function for the matrix-vector multiplication. The arguments within
each (..) of the member function is the input and destination vector/matrix for vmult and
scaling factor and destination vector/matrix for add. The handle of constraints is explicit in
the coding. Careful considerations must be done for all detailed parts of MPI coding, and
is required to be explicit in DEAL.II, while all these are hidden structures in OpenFOAM.
The MPI was due to technical issue removed in version 0.3. The adaptive mesh use the Kelly
Algorithm, described by reference within the tutorials.
3.5 InterGridMapping
3.5.1 Mapping functions
The implementation of the InterGridMapping class contain the key associative mapping
functions
sf
,
fs
between face centers x
f
of uid mesh and the structure x
s
on the cou-
pling surface on respective mapping for the transfer steps. A proximity distance classier
for the denition of these mapping functions is implemented,
(x
s
)
sf
= {x
f
: min
x
f
(x
s
x
f
)}, (3.2)
(x
f
)
fs
= {x
s
: min
x
s
(x
s
x
f
)}. (3.3)
3.5.2 Transfer step
The inverse distance interpolation was applied in Dirichlet transfer step,
(X) =
i
Xx
i
k
Xx
k
i
, (3.4)
which provides an ecient and accurate weight formula. However, due to technical details
omitted for the evaluation of the traction in the Neumann transfer, plain average of the eld
over the uid faces within given solid face was performed,
=
1
N
i
. (3.5)
Chapter 4
The case study
The case is a velocity driven uid domain of rectangular shape with cantilever(s) as a structure
immersed and attached to the wall. All values are in SI units.
4.1 Beam theory of the cantilever
The dimension of the cantilever is (1 1 10) D
3
. Using beam theory, the Bernoullis
kinematic assumption in an one-dimensional dynamical continuous system with cross section
area A
t
and centered path of line is applied [23]. The u
l
states the deection for uniformly
applied traction force over the long-side, and u
t
is the deection by force tangential at the tip
of structure. The formulas for these deections are given as follows, where pressure dierence
between front and wake side of the cantilever is denoted as P ,
u
t
=
FL
3
8EI
= k
t
P
E
, (4.1)
u
l
=
FL
3
3EI
= k
l
P
E
. (4.2)
The corresponding vibrational frequencies,
2
=
(2)
2
T
2
= 12.36
EI
A
t
s
L
4
T = k
s
. (4.3)
The strain energy U
t
associated with deection for u
t
,
U
t
=
1
2
u
t
F =
1
6
F
2
L
3
EI
=
3EI
2L
3
u
2
t
= k
u
(P)
2
E
. (4.4)
4.2 The steady state response for a cantilever.
By rewriting the advection term in Eqn (2.6),
U U = (
U
2
2
) + (U) U. (4.5)
27
28 CHAPTER 4. THE CASE STUDY
Take Eqn (2.6) and drop time derivatives and viscosity term, apply then the Eqn (4.5)
assuming no vorticity (U = 0),
f
(
U U
2
) = p. (4.6)
This is Bernoullis equation for a steady irrotational inviscid ow. By integrating along a
streamline,
p +
f
U U
2
= C, (4.7)
Assume uniform velocity prole over the structure in front of the cantilever, the static pressure
as a function of the ow direction and small deections. The argument is as follows, by
applying the Bernoullis relation on the stagnation point for the cantilever at rest, the uid
element adjacent to the cantilever where all dynamic pressure is transformed into pressure.
Since the cantilever is allowed to move, the no-slip asserts the uid cell shall move with
the speed of the cantilever. The actual energy transfered to the cantilever is reduced by the
amount of energy required to keep the uid cell moving at given velocity, this is given by Eqn
(4.7). This applied to the amplitude response for the cantilever Eqn (4.1) gives the relation
between amplitude A by the ow using the moving uid cell, with average speed of 2Af
n
,
adjacent to the structure,
A = k
f
U
2
1
1 +
_
4U
2
k
2
2
f
f
2
n
+ 1
, (4.8)
where is a function of Re and the structure, a form factor. The dependence in Re arise from
the proportionality between front pressure and the dierential pressure over the structure.
However also depends on the size of the wake, which is a function of the form of the
structure and the Re in the uid. The f
n
can without changing the equation become forced
vibration .
4.3 Boundary condition
In OpenFOAM, BC is denoted as patch. The walls are no-slip, U
r
= 0,
movingWallV elocity is the patch to be used for the coupled boundary and fixed for the other
walls. Static pressure condition is enforced on all patches but the outlet, that is (p, n) = 0,
zeroGradient. The inlet has fixed for U, that is Dirichlet condition. The outlet has fixed p,
that is Dirichlet pressure (p = 0) and U zeroGradient, (U, n) = 0. The setup of BC implies
a velocity driven ow. These were the recommended settings taken from [28, 65]. The BC for
the DEAL.II is automatic, in the sense use the proximity distance classier in calculating the
distance of given solid face center to the coupled boundary uid and denes Neumann if within
given tolerance and Dirichlet with zero displacement otherwise. The reason for not using the
same system as OpenFOAM is that the topology is recalculated at loading/construction of
the grid, an unfortunate feature of DEAL.II. It has been noted that the tolerance setting can
be complicated by the shape of the surface and could create mismatched BC but as far the
uid cells exceed the number of solid this lower the chance this event.
4.4. SINGLE CANTILEVER STUDY 29
4.4 Single cantilever study
The case study was taken from [53], the thickness of the cantilever is D = 0.2 m and the
height is 10D, placed 5D from inlet and 2.5D from the walls while the outlet is placed 20D
from the cantilever. The wire frame of the rectangle domain is thus (26 6 12.5) D
3
.
Figure 4.1: single cantilever case study.
4.5 Multi-block cantilever study
The multi-block case is similar to the single cantilever case, but with four cantilevers where
distance between them is 1D, thus giving a wire framed rectangle of size (28 8 12.5) D
3
.
The case originated from a multi cantilever study [50].
Figure 4.2: The multi cantilever case study.
4.6. THE CANTILEVER MESH 31
4.6 The cantilever mesh
The DEAL.II mesh have some support for generate mesh. A simple structure was desirable
requiring a small code to generate the mesh in DEAL.II and then saved as gmesh format.
The mesh was then to be loaded during the FSI simulation. The simplicity in building this
cantilever mesh was therefor an incentive in chosing this application to the validate the
algorithm. The original mesh is 4x4x16 cells, see left Figure 4.3. The renement procedure is
done by dividing each cell into four cells, where the Kelly algorithm [30] is used in applying
the criteria for this adaptive meshing. In same way coarsening is performed, the rened mesh
contain information on coarser mesh, hence the coarsest mesh cannot be coarser. The loop
over such renement denes the number of steps taken. It has very close connection to the
pseudo solid displacement step in the FSI algorithm. The right Figure 4.3 gives the 8x8x32
cell structure after one renement. A more elaborate adaptive meshing is given in tutorial
step-14. The settings of the Kelly algorithm, is the same as those given in tutorial step-8.
Figure 4.3: left initial mesh, right mesh after one renement. The mesh is taken from simulation at
dierent times
The class which perform this action is coarse and refine, the parameters of tolerance is a
value obtained after a numerous testing within the community and developers of the package
[61]. The Gmesh utility [75] can be used in constructing a mesh of the cantilever, see further
details from the forum in how to dene the grid, note that Gambit also have the option to
construct readable input mesh, however, no eort have been done to investigate this further.
4.7 The uid domain mesh
The mesh is generated by Gambit version 2.4.6 using scaled tetrahedral elements, with a
structured mesh on the coupled boundary of cell size 0.02 m, growth rate 1.1 and 0.1 m as
upper limit on cell diameter. Figure 4.4 shows the mesh around the cantilever. The checkMesh
Figure 4.4: The mesh in XY at Z=1 m and XZ at Y=0.6 m.
post-Processing utility using the quality rules provided by the OpenFOAM community gives
the result that grid is acceptable however there was an indication of stretched elements, within
range of acceptance as long as one chooses a correct discretization and tolerance. There can
be seen a clear change in orientation of cells in an area around the cantilever caused by the
cap and this could implicate an area of conict since ux goes in a dierent direction and
hence introduces a non-orthogonality issue in the correction for the CD formula presented in
the theory section where it is known to cause instability.
Chapter 5
Validation procedure
Each step of the algorithm is validated separately. The solid state is veried against beam
theory. The mapping functions in the transfer steps and their interpolation is veried by direct
an evaluation on test case and a check is performed at each simulation. The FSI algorithm
is veried against frequency shift and amplitude equation.
5.1 Solid state solver
5.1.1 Test Case
The DEAL.II package support quadratic elements, cubes and square cells, although the
number of quadrature points is an optional feature in this code, as well as the order of nite
element basis, allowing under-integration respective nonlinear eects. These are features not
fully exploited in this work, but the thesis [19] implements a technique useful to eliminate
volumetric- and shear-locking as the Poisson approaches its upper bound limit 0.5. The
non-linear trial functions increase the number of elements to such a degree that the current
mapping function from face center of uid to solid is not a viable option, since the size
of number of degrees becomes too large. The test case is the same cantilever used in the
FSI study with prescribed pressure on dedicated areas, long-side and short-side. Only linear
elements are therefore considered and exact integration is used. The adaptive meshing is
performed only once, using the Kelly Algorithm, with the argument that the study concerns
steady PDE. The purpose is to study small deformations and therefore it is plausible to
assume the traction is the same within a small perturbation. The benet from this is that
the stiness matrix and the mass matrix is then calculated only once during a simulation.
The actual update is very ecient, the cost lies in the mapping function between the uid
and solid. The renement is controlled by an iterative update scheme where the cell is rened
if the residual = (Kq f) criteria is not fullled, and even coarsened if the residual is to
small. In this way, it has been noted that the renement allows you to construct a coarse or
ne mesh but both lead to the same mesh after the iterative procedure coarse and renement.
The following code was used to test the class StaticElastic.
33
34 CHAPTER 5. VALIDATION PROCEDURE
#include "InterGridMapping.H"
namespace shared
{
InterGridMapping<3>* icofsi_db_list[10] ;
InterGridMapping<3> icofsi_db;
}
#include "StaticElastic.H"
#include "dealiiGrid.H"
#include "setBoundaryIndicator.H"
#include "dealiiParameterHandler.H"
int main(int argc,char* argv[]) {
shared::icofsi_db_list[0] = new InterGridMapping<3>();
dealiiParameterHandler::Base prm(argv,"object0.prm");
prm.print();
int refinement=0;
dealiiGrid::Grid<3> newGrid;
SmartPointer<Triangulation<3> > tria;
newGrid.create_grid();
tria=newGrid.get_triangulation();
StaticElastic::TopLevel<3> solver(prm,tria,refinement);
solver.set_grid_path(argv,"object");
solver.run();
}
5.1.2 Verication of class StaticElastic
Using the cantilever as case, with coarse grid and rene 0-2 gives a mesh with 16 4 4,
32 8 8 respectively 64 16 16 cells taking linear trial functions and exact integration.
The applied traction is 1Pa, E = 20 10
3
Pa, I =
1
7500
m
4
, and the result u
c
is compared to
reference u
r
from Eqn (4.2) and Eqn (4.1) with relative error |
ucur
ur
|. The reference article
cells 256 2048 16384 Relative Error(%)
u
l
/m 0.1285588 0.1441502 0.1488933 0.73
u
t
/m 0.03430288 0.03844592 0.03968957 0.77
Table 5.1: The deection ut/u
l
for StaticElastic
use 10 10 100 = 10000 elements giving a error of 1.73% [53] for a FVM solver using
Green-Cauchy strain measure.
5.1.3 Verication of class QuasiElasticity
The quasi-steady approximation is implemented as class QuasiElasticity and corresponds
to a rewrite of the step-18 tutorial in DEAL.II with MPI removed and the code with traction
vector taken from step-12. The validation of this routine is excluded from this report due to
its limited use but for certain cases such as those in the VIV section presented, this routine
5.1. SOLID STATE SOLVER 35
is suitable. The reason for its presentation in the thesis is due to the observation of uid
driven motion in VIV where steady state is entered in the rst cycle and that an article on
the subject presented a condition for the applicability of the method in terms of Keulegan
Carpenter Number (KC) [5] which describes the relative importance of the drag force and
the inertia forces. For lower numbers the inertia dominates. In order for this approximation
to be valid, drag force must dominate over inertia, which is from a physical point of view
intuitive. The StaticElastic is the stripped version of QuasiElasticity and the same testing
was performed. The major dierence is that StaticElastic is in Lagrange description, while
QuasiElasticity is in Euler description. This solver was used in version 0.1.
5.1.4 Verication of the class DynamicElastic
The test code used for this class was,
#include "InterGridMapping.H"
namespace shared
{
InterGridMapping<3>* icofsi_db_list[10] ;
InterGridMapping<3> icofsi_db;
}
#include "DynamicElastic.H"
#include "dealiiGrid.H"
#include "setBoundaryIndicator.H"
#include "dealiiParameterHandler.H"
int main(int argc,char* argv[]) {
shared::icofsi_db_list[0] = new InterGridMapping<3>();
dealiiParameterHandler::Base prm(argv,"object0.prm");
prm.print();
SmartPointer<Triangulation<3> > tria;
dealiiGrid::Grid<3> newGrid;
newGrid.set_path(argv,"object");
newGrid.create_grid();
tria=newGrid.get_triangulation();
DynamicElastic::TopLevel<3> solver(prm,tria,i);
solver.set_grid_path(argv,"object");
double time=0;stop=1;
while(time<stop) {
solver.set_current_time(time);
time+=0.001;
solver.run();
solver.next_step();
}
}
The dynamic solution is characterized by the amplitude and frequency. To given material
constants and time steps in this test case, the order of accuracy is 10
8
. After 40-75 cycles
( 1000 iterations in total with xed-point iteration) only the last digit was changed in
energy, this was repeated for dierent densities
s
. The test was performed at version 0.5,
using the more complicated set of update solving two linear equations Eqn (2.52) and Eqn
(2.53). In version 0.8 the rst update scheme was chosen, no dierence is noted. The test
revealed shear-locking when applying the traction on the tip of the cantilever, the longitudinal
waves induced by the impact reects from the base of the cantilever in twisting mode. For
this reason, the traction was applied on the long-side. The same settings as the static system,
the time period, T
, for dierent mesh, can be obtained from the output les directly, as the
turning points, maximum deection or equivalent, when the kinetic energy is at minimum.
The Maximum Displacement VECTOR in the output was introduced as to conrm the
vibrational turning point and at same time conrm the coupling to the kinetic energy. The
result of time periods T
c
is compared to the Eqn (4.3) with relative error |
TcTr
Tr
|, The error
cells 256 2048 16384 Error(%)
T
0.01
/s 0.081 0.086 0.088 <0.5
T
0.1
/s 0.255 0.272 0.277 < 0.5
Table 5.2: Periods from DynamicElastic.
in turning point is 1 ms, due to that the time step is 1 ms. The frequency implicates that
K and M is implemented correctly into Roths equation but it remains to show that the
force contribution to the equation is correct. This is done by introducing damping to the
system. The test is performed on the smallest system 16 4 4, using the same contribution
from mass and kinetic (ie = ) in percent expressed, 5% respective 1% and the time step
changed to 0.01 s. The expected feature is that velocity is damped until it reaches zero then
the solution shall become the same as for static according to the Eqn (2.50) after rearranging
the terms. Figure 5.1 summarize the outcome.
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35
0.02
0.04
0.06
0.08
0.1
0.12
0.14
0.16
0.18
Convergence of DynamicElastic to StaticElastic
Time/sec
A
m
p
l
i
t
u
d
e

/
m

5% =0.1
Static
1% =0.1
1% =0.01
Figure 5.1: A critical damped dynamical system 0-0.30 s
5.2. ICODYMFOAM 37
5.2 icoDyMFoam
The verication of the stand alone solver icoDyMFoam itself is considered as performed within
the community of OpenFOAM, this case is the central case in 2006 release of OpenFOAM
with over two years testing with numerous articles and reports as result. The misconception
however in using a code that is not properly documented have its consequence as [44] set
improper patches, that is BC. This can also been seen in the FSI solver icoStructFoam
distributed within the OpenFOAM forum. The MovingWallV elocity and calculated patch
are vital in order for that the solver to work as intended.
5.3 setBoundaryIndicator and InterGridMapping
The auxiliary class setBoundaryIndicator uses the same proximity distance classier to dene
the BC as used in transfer steps, hence this class is veried in the same step as InterGridMap-
ping is veried. This distance classier is the key function that keeps the grid together. There-
for in every run a sanity check is performed without any extra cost for each run. In the log le
the tags Max distance FOAM2DEAL and Max Distance DEAL2FOAM gives the max-
imum distance in the into maps
fs
and
sf
, that is max
xs
(min
x
f
(x
s
x
f
))
resp max
x
f
(min
x
s
(x
s
x
f
)). For all calculations performed, this was conrmed
valid as giving the distance within half a cell size of respective mesh except for one situation.
This conrms the functionality of the proximity distance classier. The InterGridMapping
was veried further by placing two identical squares adjacent with dierent structured mesh
generated by respective solver and by hand control the results in the printout list. In the
log le there is a printout in the number of faces with Dirichlet and Neumann. The result is
veried as well giving the expected number of faces for each type of condition.
The situation observed was a glitch in the restart of the simulation, discontinuities of to the
order of the cell size due to imperfect match between uid and solid. This is an error that
can be estimated by knowing the value of the gradient of over the cantilever and the cell
size. The problem arise whenever a given point X is on the border line between two cells.
The estimation of the face centers in the Euler description of the solid as the average of
surface vertices gives a discrepancy between uid and solid which can cause a mismatch in
the mapping function and indirectly cause a jump in the displacement from the misplaced
traction vector since the topological mapping is lost and has to be recomputed at the restart,
a unfortunate feature of DEAL.II.
5.4 The FSI algorithm
InterGridMapping is veried as to mapping the traction of the uid to the correct position in
the solid, and the displacement of the solid to correct position of the uid. The same applies for
the topological settings, ie boundary condition, performed by setBoundaryIndicator class.
The icoDyMFoam and DynamicElastic solver are veried. The objective that remains is
to ensure that FSI is implemented correctly. The verication steps are as follows,
The shift in time periods is compared to theoretical values.
THe amplitude is related to the formula postulated in this thesis.
The margin of error is not increased with FSI.
5.4.1 Setup
The single cantilever case, described in section 4.4 was used and using the BC settings from
section 4.3 and mesh generated by Gambit as described in section 4.7. All initial data within
the grid are zero and icoDyMFoam is performed until steady state is reached. This is dened
as ow volume ten times the domain volume according to the User Guide of OpenFOAM
downloaded from [65] with a Co=0.5. However, during the FSI simulation, Co = 0.9. This
is a questionable choice since in the forum of OpenFOAM it has been pointed out that for
stability issue using adjustableTime step, Co< 0.5 is recommended. However, due to the
time limit of this project, the speed of the simulation was increased by using Co=0.9. Three
dierent velocities were chosen at the inlet patch, U= 1, 10, resp 25 ms
1
. The Reynolds
number (
UD
) is set to 400 and hence giving the kinematic viscosity
as 5 10
4
, 5 10
3
,
resp 12.5 10
4
m
2
s
1
. The cantilever use the 32 8 8 mesh with 2048 cells, E = 20 10
6
Pa and Poisson number = 0.3. In this study marker point it is placed at the center of tip of
the free end of the cantilever. The tolerance is set to 1 10
6
for the pressure and velocity in
the icoDyMFoam, three PISO loops, the displacement in AutoMesh to 1 10
6
, the tolerance
for the Aitkens acceleration 1 10
7
and convergence criteria 1 10
12
for the CCG solver of
the linear system in FEM. Finally, the Crank-Nicholson scheme was chosen in the temporal
discretization by setting =
1
2
.
5.4.2 A note on measurement of the period
In the forthcoming tables, Period refers to how many periods a given marker point is evolved
during the simulation, and the reason for this is to present how the time period was measured.
In table 5.3, T
pi
refer to the time of the peak p
i
in total amplitude with respect to original
mesh. T
c
is the estimated time of period. FFT refer to the DFT routine in MATLAB for the
Period T
c
Error(%)
0.5 2 T
p1
<20
0.75 2 (T
p1
T
p2
) <10
1-10 T
pn
T
pn+2
<5
>10 FFT 0
Table 5.3: The measure of Tc.
discrete Fourier analysis. The error is the estimated relative in percent to the FFT obtained
5.4. THE FSI ALGORITHM 39
for >10 period simulations, where it was assumed the DFT gave the exact result and transient
eect was estimated. The in-line driven oscillation have unknown equilibrium point and by
Figure 5.2 there is a transient movement before stationary orbit is entered, showing an initial
movement in the direction of X then due to the VIV the orbit is created as going in an ellipse
path. FFT was necessary in order to study VIV and higher frequencies in the marker points
of the cantilever.
0.5 0 0.5 1 1.5 2 2.5 3
x 10
3
1.5
1
0.5
0
0.5
1
1.5
x 10
4 (X,Y) plot of tje marker point placed at center of topside
X
Y
START
END
Figure 5.2: The trajectory of marker point position in XY space relative its reference point. The
marker rapidly along X with a slow motion in Y direction, implicating the uid lower frequency (Y)
and the in-line motion of the natural frequency (X).
5.4.3 Result
T
v
is the theoretical value of the time period given by Eqn (4.1) in vacuum. T
c
is evaluated
according to table 5.3 in the previous section. The T
r
in the table is evaluated time period
using Eqn (2.62) from [55]. For this reason the frequencies are not given in the tables. The
deviation Dev refers to the relative error of T
c
to this T
r
, that is |
TcTr
Tr
|. The A is the
magnitude of the displacement eld at marker point. The numerical error expressed in Dev
have its origin in the quality of the mesh, the transition, the boundary conditions and the
wall eect [40]. Using table 5.4 making the following optimization using least square,
T
2
c
= (A +B
s
)T
2
v
+C, (5.1)
gives the result A=1.74, B=2.54 and C=-0.004. Higher order in relative mass had no inuence.
The theoretical value is A=1, B=2.5 and C=0 in accordance to Eqn (2.62). Applying this
formula for the tables 5.4-5.6 gave the result in Figure 5.3. The rst eight points from the
left is from table 5.4, the seven to fteen is from table 5.5 and the rest from table 5.6. The
reason to optimize over a limited data is that the study of T
2
require adequate of points
covering the whole range and evenly distributed. This is not possible for U=10 and 25 ms
1
due to limited time of the project to use only one processor per simulation and which in turn
limits the study to short time periods to be used in the optimization, alsowhen taking the
transient eect in account. The grid itself show areas which probably with time will cause
divergence, mostly manifested at higher velocities as Figure 5.14 show. The frequency for the
in-line driven motion should not to given Re depend on the actual ow speed and Figure 5.6
illustrates this with a standard deviation of 2.6 10
3
s, note however the system is damping
dierently between cases and hence some slight dierence is obtained. Furthermore, in the
case of

f
s
1, the proportionality of the shift in frequency between dierent solid mass
should be directly related to the vacuum. Figure 5.5 illustrates the correct scaling of a factor
of
5 which follows by Eqn (4.3) with a relative dierence of 0.8%. The observed damping
at 1 10
2
m region is the area where the current time window of study damp the in-line
and lift up the VIV unless the energetic loss drain the movement, and thus implying the
lower amplitudes to be close early to the steady state and above, more simulation time is
required to notice the eect. One signicant result is the accurate correlation with the in-line
amplitude with the proposed Eqn (4.8) showed in Figure 5.4. Figures 5.7-5.12 is the end time
of the simulation, not the maximum amplitude.
s
/kgm
3
f
/kgm
3
T
v
/s T
c
/s T
r
/s A/m Dev(%) Period
2949 1000 1.0632 1.7442 1.4452 0.20 17 0.5
100 1000 0.1958 1.026 0.9983 0.17 3.7 0.75
50 1000 0.1384 1.016 0.9886 0.17 3 0.75
300 100 0.339 0.5840 0.4591 1.4 10
2
21 >10
300 200 0.339 0.6813 0.5537 2.8 10
2
19 6
10 0.01 0.0619 0.0864 0.0862 2.5 10
6
<1 > 10
50 0.01 0.1384 0.1940 0.1388 2.6 10
6
29 > 10
10 1 0.0619 0.0880 0.0692 1.4 10
4
21 > 10
Table 5.4: U=1 ms
1
Re=400, E=2010
6
Pa,=0.3.
s
/kgm
3
f
/kgm
3
T
v
/s T
c
/s T
r
10 0.1 0.0619 0.0863 0.0627 2.8 10
3
27 >10
10 20 0.0619 0.1509 0.1516 0.33 < 1 0.75
50 0.1 0.1384 0.1902 0.1382 2.8 10
3
27 >10
50 20 0.1384 0.2092 0.1958 0.37 10 1.5
100 20 0.1984 0.2667 0.2398 0.39 19 0.75
10 1 0.0619 0.0914 0.0692 1.8 10
2
24 8
5 20 0.0438 0.1480 0.1452 0.32 2 0.5
Table 5.5: U=10 ms
1
, Re=400, E=2010
6
Pa,=0.3.
s
/kgm
3
f
/kgm
3
T
v
/s T
c
/s T
r
50 0.25 0.0619 0.1927 0.1384 4.3 10
2
28 3
10 0.1 0.0619 0.0861 0.0627 6.4 10
2
23 7
10 3.2 0.0619 0.0908 0.0831 0.38 8 1.25
50 3.2 0.1384 0.1947 0.1491 0.43 23 0.5
100 3.2 0.1984 0.2666 0.2035 0.45 24 0.5
10 1 0.0619 0.0879 0.0692 0.15 23 2.5
10 2 0.0619 0.0923 0.0.0758 0.26 18 0.5
1000 1 0.6191 0.8192 0.6199 0.17 24 0.5
Table 5.6: U=25 ms
1
, Re=400, E=2010
6
Pa,=0.3.
0 5 10 15 20 25
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
Case
P
e
r
i
o
d

T
i
m
e
/
s
e
c
The period time of inline frequency

Eqn
T
c
T
v
T
r
Figure 5.3: Periods for 23 Cases. The upper line is the measured points, the circles is the values
predicted by Eqn (5.1) which is computed by least square to Table 5.4. The middle line () is the
predicted values from Eqn (2.62) with vacuum values using Eqn (4.3) as the bottom line.
14 12 10 8 6 4 2 0
14
12
10
8
6
4
2
0
Eqn
M
e
a
s
u
r
e
d

A
m
p
l
i
t
u
d
e
The loglog plot of correlated observed and theoretical amplitude

VIV + strong damping
inline
Figure 5.4: The correlation between dynamic pressure and deection. The measured total magnitude
amplitude v.s. Eqn (4.8). The slope is the form factor. The deviation is only for case where strong
VIV have large inuence on the data (o).
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35
0
0.5
1
1.5
2
2.5
3
x 10
3
T/sec
A
m
p
l
i
t
u
d
e
/
m
s
=50 T=0.1926 s (r)
s
=10 T=0.0868 (b)

s
=10
s
=50
Figure 5.5: The magnitude amplitude of marker point placed at center of the tip. The scaling
between case with low relative mass for
f
= 0.1 kgm
3
U=10 ms
1
. For

f
s
<< 1, the frequency
should be close to vacuum and then it scales with mass in accordance with Eqn (4.3). Note however
the wall eect scales proportional implicating a force coupling which lowers the frequency.
0.05 0.1 0.15 0.2 0.25 0.3 0.35
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
T/sec
A
/
m
a
x
(
A
)
(g) U=1 (r) U=10 (b) U=25 Re 400 T=0.1914 s T
vac
=0.1384 s

U=25
U=10
U=1
Figure 5.6: The magnitude amplitude of marker point placed at center of the tip. The independence
on U for same relative mass for given frequency. s = 10 kgm
3
,
f
= 1 kgm
3
. This is since for
same Re the ow pattern is the same, and hence with the same force prole over the structure the
same frequency but with a dierent amplitude will be observed.
Figure 5.7: The pressure eld (p) in XY for U=1 ms
1
and Z=1 m, s=300 kgm
3
,
f
=100 kgm
3
and Re=400.
1
and Z=1 m, s=10 kgm
3
,
f
=1 kgm
3
and Re=400.
1
and Z=1 m, s=10 kgm
3
,
f
=0.1 kgm
3
and Re=400.
The pressure in Figure 5.7-5.9 follows closely Bernoullis Eqn (4.7). The size of the wake
is a function of the Re, the form of the blu body as implicated by Eqn (2.12). For the
same Re and shape of blu body one can hence expect the same wake size. However dierent
deformed congurations implicates interaction between vortices which change the size of the
wake as a result. The free end of the cantilever shape the wake and expects to increase the
frequency f
n
with increasing amplitude while the conned wall decreases [40]. Further, the
pressure on respective side is nearly constant over the whole surface except at the corners,
as Figures 5.13-5.15 show. This again support the result of Figure 5.6 that is, the frequency
is independent of U, however according to Figure 5.4 the amplitude play a role, but as long
as the deviation is small, this doesnt change the wake.
Figure 5.10: The magnitude of velocity (U) in XY for inlet U=1 ms
1
and Z=1 m, s=300
kgm
3
,
f
=100 kgm
3
and Re=400.
1
and Z=1 m, s=10
kgm
3
,
f
=1 kgm
3
and Re=400.
1
and Z=1 m, s=10
kgm
3
,
f
=0.1 kgm
3
and Re=400.
Figures 5.10-5.12 are the corresponding velocity plots, the increase in speed follows again
Eqn (4.7) as expected with pressure drop in corridor adjacent to the cantilever and expected
increase in damping due to drag eect. Note the occurrence of areas of possible divergence,
which is directly related to the grid structure, see Figure 5.8 and 5.11 and compare it to
Figure 4.4 on coordinates (0.6,0.2) and (2.3,2.1) in Figures 5.13-5.18. The sinusoidal shape of
the wake gives the frequency of the instability related to the vortex shedding. The size of the
wake can be estimated by Figures 5.13-5.18. The tangential line from the tip as demonstrated
in Figure 5.17 to cutting the wake size, the low velocity inside the wake is however not a
precise measure for the wake. The pressure dierence across the structure has it major source
in change by the change of the dynamical pressure ( U
2
).
Figure 5.13: The pressure eld (p) in XZ for U=1 ms
1
and Y=0.6 m, s=300 kgm
3
,
f
=100
kgm
3
and Re=400.
1
3
,
f
=1 kgm
3
and Re=400.
1
3
,
f
=0.1
kgm
3
and Re=400.
Figure 5.16: The magnitude of velocity (U) in XZ for inlet U=1 ms
1
and Y=0.6 m, s=300
kgm
3
,
f
=100 kgm
3
and Re=400.
1
and Y=0.6 m, s=10
kgm
3
,
f
=1 kgm
3
and Re=400.
1
and Y=0.6 m, s=10
kgm
3
,
f
=0.1 kgm
3
and Re=400.
5.5. ERROR ANALYSIS 47
5.5 Error analysis
Error analysis is performed on each domain separately and then an error estimation of the
time period is presented.
5.5.1 icoDyMFoam
OpenFOAM provides a post-Processing utility to estimate the error, icoEstimateError, a
local residual formula presented in the thesis in [28]. By studying the dierence between Eqn
(2.17) and Eqn (2.19), then discretize to the second order, one obtains a transport equation
for the error ux. This utility implicates no instability is caused by the FSI algorithm, since
the error ux is the same before FSI and after. The static error originates mainly from the
convective term, by studying the case U=1 ms
1
and U=10 ms
1
from Figure 5.19 and
Figure 5.20, is increased by a factor of 10. Note however that this utility doesnt estimate
the error in the transient part which due to the high Co is of signicance. The probes in the
uid revealed high frequency error pulse, with the frequency of the vortex release, which is a
result when the vortices enter the zeroGradient outlet boundary condition, hence the pressure
correction gives a shift in total energy which is reduced when the transient process by initial
movement of the cantilever enters the steady state region. From Figures 5.19-5.21 one obtains
that the major source of error lies in the wake, the vortices and at points of interest due to
possibility of divergence in the grid, as Figure 4.4, 5.8, 6.3 and 5.11 illustrates. There were
three cases which had divergence in the PISO loop, they are connected to the error in the
convective term but there are at least four more cases that showed the same tendency, as
Figure 5.8 implicates. The PISO in icoDyMFoam within the loop doesnt correct for the
orthogonality, since it has been observed to cause instability in convergence. This is done
as a last correction [28]. Figure 5.15 show the connection between the failed convergence
point and the estimated error at the divergence points in pressure eld in the wake, note the
lowest level since that point is the common factor in all diverged cases. The relevant question
mark for Figure 5.16-5.18 is the origin of this regular error at the cantilever at three dierent
heights. An analysis of the axial mode of frequency, indicates a relation to the release of the
axial modes of vorticity.
5.5.2 DynamicElastic
FFT analysis on the marker points for analysis revealed that unknown high frequencies with
an amplitude strength relative to the natural frequency of order of 1:1000, only notable in Z
direction, with no correlation to
f
and varies linearly with the
s
. The frequencies of question
are 582 Hz,732Hz and 1782 Hz. This implicates the possibility that excited axial modes
appear. The Crank-Nicholson error caused by staggered technique in state space creates an
error with high frequency error, and from the verication of DynamicElastic showed no such
feature. This favor the interpretation that the observed frequencies are excitations from the
variation of pressure due to the vorticity acting on the tip of the cantilever and since the large
dierence in frequency between uid and natural frequency for the solid, the mode seen in
the marker point is the natural frequency for longitudinal mode. The error from the probes
with the signal noise caused by the vortices is a constant shift, this doesnt aect the solid
state since its movement is built upon the pressure dierence, hence canceling out the noise.
5.5.3 Estimated Error in the period
By taking the sum of the relative error to the contributing Eqn (4.3), the total estimated
error E for T
c
is the sum of the modeling error of solid (T
n
=
1
fn
), estimated by table 5.1,
the force term in uid U
2
(U) estimated by Figure 5.19-5.21 from the same cases presented
in section 5.4.3 and the measurement error (T
c
) from table 5.3,
E =
U
U
+
T
n
T
n
+
T
c
T
c
. (5.2)
This gives a relative error of the order 5-25%,where on average 3% error arise from the U,
the T
n
have an error around 2% while the measured time period lies between 5 20%. The
degree of inuence of the wall on Dev is a relevant issue but the present error bars are
to large to actually answer this question. The major source of error is that the simulation
is limited in time, only partial periods evolved, and the transient eect as seen by Figure
5.2. The number of period required is not that cannot be estimated, its something that will
be revealed during the simulation. The study indicates although a smaller error due to the
following observations,
All shift by Eqn (2.62) have the correct sign.
The scaling eects in parameters are precise, Figure 5.5.
The amplitude formula Eqn (2.63) show less average error.
The optimization Eqn (5.1) implicates a systematic shift.
It is therefore a system dependent error that is relative and proportional to the quantity
observed, hence implicating either an error in the mesh or its related to the boundary
condition of inlet or the wall eect [40]. The time period is not likely aected by the poor
resolution of the Karman eect, the frequency depends on the dynamics of deection of the
uid around the object, an area that has a good grid resolution, which is the above error
estimate for the uid taken from the Figures 5.19-5.21. The error is mainly generated by
the wake and from plots the axial modes of release on vortices is clearly seen, due to the
interaction with the free end one can see a partial mode A with nodal plane at the midpoint
of the cantilever in Figure 5.22 and 5.24. Figure 5.23 is beginning to diverge and this pattern
with regular error points on the surface is a peculiar feature in divergence in the pressure
correction. The divergence is often a slow process, it gradually increase the pressure until a
point where pressure correction fails. Due to CD there is divergence originating from non-
orthogonal cells note especially the point in XY plane at (0.6,0.2) and compare it to Figure
4.4.
Figure 5.19: The estimated velocity error in XY for U=1 ms
1
3
,
f
=100
kgm
3
and Re=400.
1
and Z=1 m, s=10 kgm
3
,
f
=1
kgm
3
and Re=400.
1
and Z=1 m, s=10 kgm
3
,
f
=0.1
kgm
3
and Re=400.
Figure 5.22: The estimated velocity error in XZ for U=1 ms
1
3
,
f
=100 kgm
3
and Re=400.
1
3
,
f
=1
kgm
3
and Re=400.
1
3
,
f
=0.1 kgm
3
and Re=400.
0 0.5 1 1.5 2
0.25
0.2
0.15
0.1
0.05
0
0.05
0.1
0.15
0.2
p
*
over time at probe 55 U=1
s
=300
f
=100 Re 400
Time/sec
p
*
/
(
p
/
)
Figure 5.25: Probe 55 in the uid in U=1 ms
1
. This plot comes from the VIV study in the
next chapter, the in-line motion is damped which could explain the lowering in size of the white
noise in Figure 5.25. The reason for this behavior is that when the vortex enters outlet, due to the
zeroGradient patch on outlet for velocity, it introduces an error which is showed as noise as long
it passes over the gradient. The whole spectra of error in the probes was not fully studied due to
limited time for the project.
Figure 5.26: The divergence at U=10 ms
1
in XY at Z=1.25 m s=10 kgm
3
,
f
=1 kgm
3
and
Re=400. Several cases provide an error as exemplied in Figure 5.14, which this gure illustrates
close with a XY cut at Z=1.25 m. Exact same type of divergence was observed while using semi-
discretization with =
1
2
for the uid a study omitted in this report. This is well documented within
the OpenFOAM community over the classical oscillation in using CD together with Crank-Nicholson.
5.6 General comments
A brief discussion on the blending in time semi-discretization, the dierence between explicit
and implicit FSI followed by an analysis of the used acceleration techniques.
5.6.1 The theta method
The theta method gives explicit, implicit and Crank-Nicholson semi-discretization of the time
variable. The damping introduced by implicit Euler ( = 1) compared to Crank-Nicholson
( =
1
2
) is notably larger as showed by Figure 5.27. As indicated during the verication
of the DynamicElastic, =
1
2
is energy conservative. The damping observed for =
1
2
in
Figure 5.27 comes from the implicit temporal discretization in the uid and the drag eect.
The explicit Euler method = 0 failed due to the induced counter pressure originating
from the Bernouillis relation (P =
f
1
2
v
2
), that is when velocity v of the adjacent uid
cell to structure changes, a pressure change P will be introduced with opposite sign as a
function of the speed v of the structure. As a rough estimate, when the v approaches the
same velocity as the inlet, the explicit solver will diverge. This implicitly set the upper limit
for displacement taken to be U t which in turn limit the time step. However, the inertia of
the solid structure stabilizes the system, the above relation is for
s

f
, higher stability is
achieved by
s
>
f
.
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0
0.5
1
1.5
2
2.5
x 10
4 U=1 Re 450 , =0.5 (r) =1.0
T/sec
A
m
p
l
i
t
u
d
e
/
m

=0.5
=1
Figure 5.27: The Implicit semi-discretization in time ( = 1) compared to Crank-Nicholson ( =
1
2
)
for U=1 ms
1
, s = 10 kgm
3
and
f
= 1 kgm
3
with Re=400. The signicance of the damping of
the uid is illustrated since with the current mesh for the structure the damping in the structure is
insignicant.
5.6.2 Explicit/Implicit FSI
The FSI implemented oers explicit and implicit coupling solution, also known as weak
and strong coupling. Test calculations have shown that with damping there is no signicant
dierence between the methods using the settings from section 5.4.1, the same convergence
pattern. The explicit case with
s

f
shows same divergence for same reason as the explicit
time step in the semi-discretization in section 5.6.1. The gain in using implicit coupling is
5.6. GENERAL COMMENTS 53
together with ROM step gaining an accuracy of several order and a gain in computational
speed is expected if the tolerance could be chosen adaptively.
5.6.3 Aitkens relaxation and the ROM
For

f
s
> 1 and for exible structures, the Aitkens relaxation method is of signicant impor-
tance, where a gain in speed of convergence by a factor 4 is observed or even more. These
reproduce the values of iterations reported elsewhere [53, 33], although higher order errors
will be introduced for non-linear functions due to the interpolations, indicating the favor
for ROM techniques. As indicated from [33] where a diagram shows that not until higher
tolerances than 1 10
6
a gain was obtained in iterations in the sub-cycle step by using
ROM, but due to the nature of this problem, there was little need for having higher preci-
sion unless studying higher frequencies. The Aitkens relaxation normally converged within
2-7 iterations, exceptions occurred whenever
n
0 resulting in a drift with an increased
error on average 20-30 times the tolerance, with an error of the order 10
5
in present study,
negligible for all practical consideration. In applying the 3-ROM with its FD formula gave no
3.9 4 4.1 4.2 4.3 4.4 4.5
8
6
4
2
0
2
4
6
8
10
x 10
6 U=1 Re=400
s
=10
f
, E=20e6, (b) without damping (r) 1% damping
T/sec
A
m
p
l
i
t
u
d
e
/
m

without damping
1% damping
Figure 5.28: Marker at the cantilever center point on the tip (X component). The regular pattern
appearing as noise, is the high frequency from the longitudinal waves in the cantilever (Z-component)
and the eect of damping at 1%.
gain in computational speed between time steps due to its inherited divergence property with
a scaling factor of 3. By taking the test function e
x
it can be demonstrated that the current
3-ROM worked excellent as long as the interpolation step was exactly known, however, due
to the coupling term, the dynamics change the interpolation and introduces an error and
this made the interpolation collapse since the deviation was increasing creating a counter
pressure eect which made the solution crash within a few cycles. Similar behavior as to
explicit semi-discretization except for three simulation followed by one ROM step. The origin
to the disastrous result with 3-ROM is related to the Courant number and the introduced
signal to noise ratio, see Figure 5.28, a small error in the displacement introduce error in the
uid solver. In other words the FD strategy is not a good choice as interpolation technique.
The reason why the sub-cycle iteration dont diverge is that Aitkens relaxation auto correct
the error and at the convergence point for the given threshold in a sub-cycle, a signicant
improved accuracy is gained during the up to two orders since the formula for FD works
well when the coupling eect can be ignored between the uid and the solid. By introducing
damping the performance can be improved between the time steps making the error decay
during the normal FSI cycles. However, a proper analysis is required to not deteriorate the
physics, as Figure 5.28 illustrates with over-damping. An improved ROM technique involves
either a Proper-Orthogonal-Decomposition technique (POD) [16] or a residual method as
Gauss Newton [54].
5.7 Multiple block case
The multi-cantilever case from section 4.5 was chosen in this study using section 4.7 and
4.3. The main purpose of this study is to reproduce the frequency of the single cantilever
with four cantilever case. However, a complicated issue with the laplacian schemes for the
pseudo-displacement scheme and their pressure correction causes divergences in some of the
cases. But all the converged result shows, that the cantilever is moving towards each other as
to be expected by theory. The lines with larger amplitudes in Figure 5.29 are the two front
0 0.05 0.1 0.15 0.2 0.25
2
0
2
4
x 10
4
T/sec
A
m
p
lit
u
d
e
/
m
Case_IV_A :
s
=10
l
=0.01 E=2010
6
: X direction
0 0.05 0.1 0.15 0.2 0.25
5
0
5
x 10
5
A
m
p
lit
u
d
e
/
m
T/sec
Case_IV_A :
s
=10
l
=0.01 E=2010
6
: Y direction
0 0.05 0.1 0.15 0.2 0.25
1
0
1
2
3
x 10
5
A
m
p
lit
u
d
e
/
m
T/sec
Case_IV_A :
s
=10
l
=0.01 E=2010
6
: Z direction
Figure 5.29: The marker point coordinates for U=1 ms
1
, s = 10 kgm
3
,
f
= 0.01 kgm
3
and
Re=400. Front cantilevers are with larger amplitude and those behind in the wake have smaller.
cantilever and those lines with smaller amplitude are the cantilevers behind, indicates less
deformation for the cantilever in the wake of the front, and all four even phase shifted have
the same time period, comparable to the one cantilever case, T
c
=0.0768 s for U=1 ms
1
and T
c
=0.1053 s at U=10 ms
1
. Eight more cases had slow convergence, such that not even
a half period evolved. The convergence increases linearly with the number of cantilevers,
the Aitkens relaxation converges within 5-20 iterations on average. In total one third of all
calculations diverged.
5.8. CONCLUSION 55
5.8 Conclusion
The FSI reproduces the frequency shift well within the margin of error of the mesh quality
and the measurement of the time period. The amplitude correlation study, the scaling and
the independence of velocity in the frequency increases the condence that the FSI algorithm
works as intended. The error analysis from section 5.5.1 is incomplete since it lack an estimate
of the transient eect and not study the renement eect to conrm the convergence in
parameters is achieved. However, the error lies as expected within the wake, the vortices and
areas of possible divergence in the grid and it increases linearly with the ow speed. Figure
5.25 implicates the presence of transient error which is decaying and hence the dynamical
path in momentum space could not be reliable but the pressure change is more and less
a constant shift and dont aect the pressure dierence. The CD is used in the numerical
boundary conditions and the zeroGradient ux BC is the source of the discretization error
[28]. Although the method is stable an issue appear while using two sequential solvers, see
Figure 5.30. The dilemma arise due to that the solid is more accurate and stable than the
uid. Therefor, the upper bound for the time step is given by the uid solver. However,
due to a combination between temporal discretization and the transfer step where viscous
traction is determined by displacement eld which where the interpolation formula have
higher error than the solid solver, a lower bound is obtained from the solid. This implies that
the total algorithm, seldom achieve the individual solvers optimal performance and there
will be situations that for given tolerance, there is no convergence although each individual
converges.
Figure 5.30: The numerical precision entering the point of oating error. The characteristic upper
and lower limit in grid size for a given tolerance. Taken from the course assignment 1 in FMN110
(2007). [14]
Chapter 6
An application to VIV
The lock-in frequency of the uid driven motion on the structure is studied in this chapter.
By using probes on the structure and in the uid the coherent movement is veried by FFT
and the result is compared to the experimental uid frequency to given Reynolds number.
Eqn (4.4) implicates over 90% of the strain energy in the validation steps originates from
the in-line motion and therefore an articial damping must be used to unmask the VIV
components.
6.1 Problem description
The dimensionless parameter Strouhal number (St) relates to the frequency of the vortex
release caused by the deection of the streamlines around the immersed body. The larger the
body the smaller the frequency, the larger velocity the higher the frequency. The characteristic
St for VIV was observed experimentally to be [68],
St =
fl
U
= 0.198(1
19.7
Re
). (6.1)
The vorticity can be identied by the termU in Eqn (4.5), which describes the rotation of
the uid, giving cause to the Karman vortex street Figure (1.1). In a simplied form, it can be
expressed as the energy that cost to divert the ow is the energetic content in the vortex and
hence the induced strength in the transversal force. The St depends on the velocity (U) and
the characteristic length (l ) of the body. The vortex has a three dimensional structure, and
is released from the cantilever periodically along the centerline in modes A and B discussed
in [21] implying that any point behind the cantilever in the domain of dependence can be
representative to evaluate the St, from the variation of the pressure. The classical description
of the modes 2S,2P and so forth is obtained from experiments using thin layer ow (2D),
however, the mid point of the cantilever due to the nodes of A and B gives the same two
dimensional ow pattern.
6.2 Setup
For the settings of the case, the section 4.3, 4.4 and 4.7 is used. The in-line driven ow
motion of the cantilever is nearly decoupled from the vortices released from the sides of the
cantilever. For this reason the FFT spectra in Y direction gives the VIV frequencies. Table 5.3
requires the correct case to be chosen such that at least 5 vortices released during simulation
57
58 CHAPTER 6. AN APPLICATION TO VIV
to estimate the frequency for the cases which revealed the VIV due to the damping in given
time window, that is an amplitude of order 1 10
2
m. For this reason, the following cases
were chosen, U=1,10 and 25 ms
1
with
s
as 300, 10 and 10 kgm
3
respectively
f
as 100,
1 and 0.1 kgm
3
. One particular feature of simulations from section 5.4.3 is the dominance
in the in-line driven motion, the actual VIV is of an factor 100 or less in displacement length
from starting position. By using the moderate damping such that the VIV becomes more
dominant since energy is continuously introduced by the energy content of the vortex. It is at
this point clear that without damping in uid or structure, the amplitude will have negative
damping. The calculations were therefore performed with 0.1% damping such that the in-
line driven motion is drained and transversal motion remains. The observed frequencies in
the cantilever are those of the marker probe placed at the center of the tip of the free end.
These are compared to the averaged spectra of Fourier analysis for the 120 probes in the
uid aligned from six lines in X direction, from middle height and top cantilevers free end,
3 lines at each level, taking the middle point and corners of the surface, 20 points starting
with 0.2 m distance from the cantilever and 0.2 m between each point. The Fourier analysis
is performed as far as possible the accuracy allows from the transient part of the spectra.
6.3 Result
The study have its origin from the case
s
= 10 kgm
3
and
f
= 0.1 kgm
3
at U=1
ms
1
, Figure 6.1 from which no signicant exchange in energy between uid and solid, show
a particular wave packet in Y direction, which revealed the VIV for the rst time in this
project. There is no observed damping in the system since uid is in steady state ow with no
signicant exchange, its a perturbational response to the forces and the cantilever is loading
up the frequency and damped in time by the inherited energy loss in the discretization. There
will be a damping observed in the marker point data if the coupling term with the uid is
removed. The probes data in the uid data contained a white noise signal in period, see
Figure 5.25. This is with a frequency similar to the release of vortices, which itself is not
a problem since INS is independent of absolute data and the FSI is driven by the pressure
dierence over the cantilever. However a question arise of the signicance of the eect in FSI
as when structure enclose the uid like ow in a pipe, then the response is based on absolute
data. The FFT analysis on marker point in solid for U=1 ms
1
gives a VIV time period of
1.45 s (0.69 Hz) as seen from Figure 6.9 red line (Y) while the blue line (X) is masked by the
in-line movement. As the frequency of VIV can be read either from the white noise signal
mentioned above and in section 5.5.3 it can be compared to the average FFT spectrum of
the uid probes given by Figure 6.11 to be around 1-0.7 s (1-1.3 Hz). The FFT spectra for
the solid Figure 6.13 with U=10ms
1
gives around 0.16 s (8 Hz) time period in the solid
and for the probes in the uid Figure 6.14, which is estimated to be around 0.07-0.14 s (7-14
Hz). From Figure 6.16 the FFT spectra U=25 ms
1
case give VIV at 0.064 s (15 Hz) and
the averaged uid probes FFT analysis gives a value around 0.04-0.06 s (17-25 Hz), Figure
6.17. As clearly seen from the Figures 6.11-6.17, the uncertainty is relatively high for the
uid probes. By the Eqn (6.1) one can conclude the frequency should be 1 s (1Hz), 0.1 s
(10Hz) resp 0.04 s (25Hz). However the validity of the formula is somewhat questionable,
the conning walls and the shape of the blu object have inuence as well as the free end
of the cantilever. Out of the previous 23 simulations it was observed that 12 simulations
were observed with VIV. The analysis of these is complicated by not using damping in these
cases hence the amplitudes have relative large range of error. The time period is supposed
to increase due to the mixing of the wakes [40] but also the free tip wake mechanism will
introduce a shredding of dierent frequency [31]. It was also observed multiples of frequency
6.3. RESULT 59
0 2 4 6 8 10 12 14 16 18
5
0
5
x 10
6
T/sec
A
m
p
l
i
t
u
d
e
/
m
Case_III_A :
s
=10
l
=0.1 E=2010
6
: X direction
0 2 4 6 8 10 12 14 16 18
2
1
0
1
2
x 10
7
A
m
p
l
i
t
u
d
e
/
m
T/sec
Case_III_A :
s
=10
l
=0.1 E=2010
6
: Y direction
0 2 4 6 8 10 12 14 16 18
1
0
1
2
3
x 10
8
A
m
p
l
i
t
u
d
e
/
m
Case_III_A :
s
=10
l
=0.1 E=2010
6
: Z direction
Figure 6.1: A marker point spectrum from validation step, U=1 ms
1
, s=10 kgm
3
,
f
=0.1
kgm
3
and Re=400, with X,Y,Z components over time from top and down. Although the result is
non-physical due to out of range of continuum, it had numerical the purpose to prove the stability
and that the natural frequency was correct. Note however the Y component (middle gure) contains
a wave-packet where the lower frequency is the VIV frequency and higher frequency the natural
frequency.
of VIV, most clearly seen in Figure 6.16 and the origin for this lies partly in the failure of
exactly representing a sinusoidal movement and the instability of the wake, that gives higher
order terms in the spectra as multiples of the lower. Using the 12 cases with VIV with the
diculty in nding a measurable A
y
component show that the study had 7 unique points that
could be used for a plot of reduced velocity versus the reduced amplitude and the quotient
fn
f
[58]. The reduced velocity plots in Figure 6.8 reproduces the linear dependency while
the amplitude response cannot be relied on due to the wall eect, apart from one thing, the
point of reduced velocity where the stationary ow induce an unsteady motion. The plot of
V
r
versus
fn
f
show a linear dependency, but since no true lock-in frequency case been studied,
that is f = f
n
, the attening of the curve is absent in this study [58]. One interesting point
is the dierence between the XY and the XZ plane, the impression is given that the size of
vortices decreases linearly by time in the XZ plane. This is a false picture, the purpose for
showing these gures is to show the importance of the vortices released from the tip.
Figure 6.2: xUy yUx in XY for U=1 ms
1
3
,
f
=100 kgm
3
and
Re=400.
1
and Z=1 m, s=10 kgm
3
,
f
=1 kgm
3
and
Re=400.
1
and Z=1 m, s=10 kgm
3
,
f
=0.1 kgm
3
and
Re=400.
From Figures 6.2-6.4 the Z-component of the vortex show the alternative spin in going
from the wall to the closest point of release on the cantilever and the opposite spin over
the structure and thus shaping the wake. From this it is clear that the wall will give the
impression of stable vortex release. The release is unstable meaning increasing size of the
vortex without the walls in time [68]. The tip have strong vortices release and aects the
wake size as a function of the deection. Figures 6.2-6.4 show same size of wake due to Re,
however Figure 6.3 is smaller as expected due to the free tip for cases with larger amplitudes.
Figures 6.5-6.7 show the inuence of the free end tip of release of vortices, tangential to the
surface, however only an animation gives the dynamics clear picture that the vortices going
from the tip are tangential from the surface of the tip.
6.3. RESULT 61
Figure 6.5: total vorticity in XZ for Y=0.6 m at U=1 ms
1
, s=300 kgm
3
,
f
=100 kgm
3
and
Re=400.
1
, s=10 kgm
3
,
f
=1 kgm
3
and
Re=400.
1
, s=10 kgm
3
,
f
=0.1 kgm
3
and
Re=400.
0 2 4 6 8 10 12
0
0.5
1
1.5
x 10
3
V
r
A
/
D
f
Reduced velocity plots
0 2 4 6 8 10 12
0
0.5
1
1.5
V
r
f
/
f
n
Figure 6.8: Reduced velocity Vr(
U
fnD
) plots, The lower gure is Vr v.s.
f
fn
and upper Vr v.s.
A
. The
gure reproduce the linearity of the lower and the initial branch for the upper, i.e. the region where
the VIV enter resonance with the structure.
The reduced velocity plots cannot due to the wall eect which stabilize the vortices,
reproduce the data for the amplitude response, since it is related to a switch between 2S
and 2P mode which is not plausible due to the conned wall, but it is expected to give
the linearity in the reduced velocity and the reduced frequency since this is independent to
the mode. However, the point at which strong interaction is predicted to occur seem to be
indicated in Figure 6.8. Since the data was insucient for a proper plot, the Eqn (4.8) was
used for a uid with density equal to one, the amplitude was hence divided by uids density.
From Table 6.1 one can see the correct scaling in frequency due to Eqn (6.1), however, in
Inlet velocity (U/ms
1
) Fluid frequency (f/Hz)
1 0.79, 0.74, 1.47, 1.53, 1.53, 1.53
10 6.3, 6.6, 5.2
25 17.07, 17.2
Table 6.1: VIV frequency for 12 case from validation step.
order to prove it is uid driven motion, the uid frequency for the same run must be veried,
since it was unclear at the validation step.
6.3. RESULT 63
0 1 2 3 4 5 6 7 8 9
1
0.8
0.6
0.4
0.2
0
0.2
0.4
0.6
0.8
1
The Marker (X,Y,Z) U=1
s
=300
f
=100 Re=400
Time/sec
A
/
m
a
x
(
a
b
s
(
A
)
)
Figure 6.9: The X,Y,Z components of the marker point at U=1 ms
1
. The lower graph is Y
component.
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5
0.2
0.4
0.6
0.8
1
FFT spectrum U=1
s
=300
f
=100 Re 400
Frequency/Hz
R
e
l
a
t
i
v
e

a
m
p
l
i
t
u
d
e

X
Y
Z
Figure 6.10: The FFT on marker point at U=1 ms
1
. FFT on Y have maximum peak at 0.7 Hz
0 2 4 6 8 10
1
2
3
4
5
6
7
8
FFT average spectrum on Probes for
s
=300
f
=100 Re 400 U=1
frequency/Hx
A
*
Figure 6.11: The averaged FFT spectrum for the uid probes for U=1 ms
1
.
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7
1
0.8
0.6
0.4
0.2
0
0.2
0.4
0.6
0.8
1
s
=10
f
=1 Re=400
Time/sec
A
/
m
a
x
(
a
b
s
(
A
)
)

X
Y
Z
Figure 6.12: The coordinates of the marker points for U=10 ms
1
. The lower graph is Y component
0 5 10 15 20 25 30
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
The FFT spectrum U=10
s
=10
f
=1 Re=400
A
*
/
m
a
x
(
A
*
)
Frequency/Hz

X
Y
Z
Figure 6.13: The FFT on marker point for U=10 ms
1
. FFT on Y have two peaks, the natural
frequency and uid frequency.
0 10 20 30 40 50
100
200
300
400
500
600
700
frequency/Hz
A
*
s
=10
f
=1 Re 400 U=10
Figure 6.14: The averaged FFT spectrum of the uid probes for U=10 ms
1
.
6.3. RESULT 65
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35
1
0.8
0.6
0.4
0.2
0
0.2
0.4
0.6
0.8
1
s
=10
f
=0.1 Re=400
Time/sec
A
/
m
a
x
(
a
b
s
(
A
)
)

X
Y
Z
Figure 6.15: The coordinates of the marker point for U=25 ms
1
. The lower graph is Y component
0 5 10 15 20 25 30 35
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
The FFT spectrum U=25
s
=10
f
=0.1 Re=400
A
*
/
m
a
x
(
A
*
)
Frequency/Hz

X
Y
Z
Figure 6.16: The FFT analysis of the marker point for U=25 ms
1
. FFT on Y gives two peaks, a
failed representation of the sinusoidal shape.
0 50 100 150 200
0
1000
2000
3000
4000
5000
6000
A
*
s
=10
f
=0.1 Re 400 U=25
frequency/Hz
Figure 6.17: The averaged FFT spectrum for the uid probes for U=25 ms
1
.
Chapter 7
Discussion
This report presents a method to resolve the uid-structure interaction (FSI) using a xed-
point iterative scheme with a partitioned Gauss-Seidel technique accelerated with Aitkens
relaxation method. The observed stability of the method is limited by the inherited require-
ment to meet the stability in the OpenFOAM and DEAL.II, a classical feature when using
a coupled solver which has dierent oating point error propagation. The validation of the
solver involves,
reproduced frequency shift in in-line movement to Eqn (2.62).
matched frequency in VIV with probes of uid and solid.
reproduced VIV frequency with regard to Eqn (6.1).
correlation in magnitude of amplitude with regard to Eqn (2.63).
independence in frequency with respect to ow velocity.
scaling in frequency cp to T
v
for

f
s
1.
The open source packages with their libraries and tutorials, DEAL.II and OpenFOAM already
contain the routines to solve the equations in question, and thus the average change in existing
code les necessary to merge the packages is a minor coding eort. One of the goals was to
learn to use the packages as to become tools in future projects. The actual thesis is fullling
that goal. The essential achievement is thus the implementation and the verication of the
FSI algorithm on above points. The Eqn (4.8) is of limited use, restricted to the cantilever
case and the same conditions to fulll the Bernoullis equation. The application to the multi-
cantilever case show promising feature and exibility for a small additional cost of coding.
The application to the VIV is rather a process in verication of the success in merging
the packages than an actual achievement in itself since the mesh is too coarse to have any
scientic value. However, during this verication some insight of the stability of the method
and model eects as boundary conditions and wall, revealed plausible future research elds
which are not covered fully in the literature. In all strong VIV cases, a damping is observed,
more importantly, not observed in other cases as clearly, which is more an coincidence since
the damping of the same order for all cases, reducing it to the steady state level for low
amplitudes (< 1 10
2
m), and in the region of 10
2
m the in-line motion is damped to the
steady state level but the VIV is stronger in energy introduced by the diverted ow around
the body and therefore visible clearly.
67
68 CHAPTER 7. DISCUSSION
The implicit coupling and the explicit scheme in the FSI is relatively similar in convergence
and stability for relative mass

f
s
< 1 and moderate large E, in other words, unless the
pressure change by the xed-point iteration is of comparative order as the original pressure
dierence due to the ow. This instability arise whenever a counter pressure induced by
deformation in one sub-cycle iteration creates a self induced vibration which is diverging for
undamped case. The convergence pattern for the Aitkens relaxation is fully reproduced as
Tukovic [53] and K uttler [33] reported were to given relative mass, same level of iterations
in sub-cycle was observed. The ROM method didnt improve in this region of precision, the
convergence, but a signicant decrease was observed in the nal error in the nal iteration
by one to three order of magnitude than the tolerance. It is expected that the use of a better
interpolation technique would improve the result.
The application itself was pressed to its limits by the Co close to 0.9, the recommended value
is less than 0.5 stated in the manuals. The higher order frequencies in the marker point FFT
spectrum, implicates a Crank-Nicholson error. However from the study of the isolated solid
cases revealed no such bias. Its therefor believed to be a physical interaction described by
the uctuations of the vortices release at the top of the cantilever giving longitudinal waves.
Some of the problematic issues in the simulations encountered can be summarized as,
The drifting in convergence by
i
0 is caused by inadequate setting of tolerance in
AutoMesh and DynamicElastic.
Discontinuity in restart caused by the lost topological mapping between simulations.
An unfortunate feature of DEAL.II.
The coarse uid mesh introduce large error at higher velocities and introduces the
back-ow error.
Vortex induced error in time due to zeroGradient patch.
Mesh deformation of non-convex domain causes divergence in pressure correction.
PISO correction extensive in FSI. One extra loop introduced to gain stability compared
to no FSI.
Error analysis showed that FSI gives no additional instability and error to the ow. The error
is mainly originated by the convective term and it increases linearly with the inlet velocity and
some cells with possible divergence in the grid caused by the scaling from the cantilever gave
points of interest where spontaneous uctuations of velocity creates oscillations of pressure
in similar fashion as the vortex-induced oscillations at the outlet but in this case they are
clear in the spatial discretization.
The reader should be kindly reminded that there is no problem to nd simulations that dont
have the presented divergence properties while discussing numerical error in this thesis, but
the error analysis combined with the limited time in the study justify the presence of these
aws in the result. This as long the error can be estimated and that the FSI algorithm
does not amplify the error. On the other hand, for the purpose of validation, the actual
presence is more a strength than a weakness of the result. The discretization error due to CD
was exemplied showing point of spontaneous uctuations around the corner in the course
material for the project course TME050 at Chalmers [71]. These are similar to those seen in
Figure 6.3 and 6.4 around the cantilever.
The failure in convergence of the Aitkens relaxation occurs whenever the correlation
i
from
Eqn (2.66) goes under a threshold and the iterative scheme is aborted. The eect of this is
69
then a drift of the order of the tolerance times the number of iterations from the previous
converged point. By choosing the tolerance in the uid mesh Automesh compared to the
tolerance of convergence in Aitkens relaxation to one order magnitude less, often resolves
this but not always. This since in the multi-block case it has been observed several cases
of divergence in the pressure correction step. The discontinuity shift observed at the restart
of simulation is due to the auxiliary class setBoundaryIndicator, this is directly related to
the proximity distance classication using a simplied estimation of the displacement, the
average, creating a shift, an o-set error. However the eect is marginal and can in most
cases be corrected by post-processing the data unless the shift itself induces a change in the
oscillation, observed in one case of 23 in the Y direction, one of the strong VIV cases. The
convergence problem in the multi-block case is mostly related to the PISO correction and ALE
deformation using elliptic integral technique, these involve a laplacian to be solved, which
becomes unstable and by increasing the number of PISO loops a signicant improvement is
achieved (>2), or by setting the tolerance for the deformation of mesh point to an adequate
level (trial and error). The extra PISO loop compared to steady ow calculation with no
FSI was also necessary for the single cantilever case, and in one simulation, the pressure
correction failed with divergence as a result. The exact cause to this error is unknown but
as indicated in literature, when the mesh is too ne, and the change is too large, a change
in the nodes which alter the topological mapping is more prominent in non-convex domains.
An observed complication is the back ow at the outlet, this phenomena arise from the zero
gradient boundary condition causing divergence. The largest source of error is however the
limited time in simulation, forcing to high Co and rough estimates of the time periods, before
transient movement enter the steady state movement.
70 CHAPTER 7. DISCUSSION
Chapter 8
Future work
The largest source of error lies in the time required to simulate the periods to a desired
accuracy to estimate the frequency, this implicates the need to use the Message Passing
Interface (MPI). This allows parallell processing, a desirable and feasible feature. The current
major source of model error is the mapping between the mesh and the interpolation of transfer
steps. The accuracy will be improved in the transfer by changing the current transfer step
to a mapping between uid face center and solids face quadrature points. This would resolve
the glitch in the restart of a job and the associated error in the evaluation of the viscous
traction. This due to the articial blown up error by using the displacement eld as source of
boundary eld for velocity for high convective case. An impact of this is that it also allows to
reduce the solid mesh size by allowing higher order elements to be used. As for example, the
cantilever, by including third order elements is sucient, allowing one block is to be exact,
within the Bernoullis kinematic assumption of interaction.
implement MPI for parallell processing
mapping between uid face center to solid quadrature points
change strain measure
implement (POD) ROM
investigate the mechanism in pressure coupling in current FSI
implement BC for the multi-block case
implement the bueting eect on FSI
implement under-integration
adaptive meshing in uid domain using icoEstimateError
This is currently impossible since then only one pressure point will be transfered, the average
pressure on each side. The deviation in frequency in this project was the observed bound-
ary condition eect on the pressure and the wall eect, this is a suitable application of the
algorithm. The use of ROM to reduce the computational cost is a complicated feature since
due to white noise and FD scheme implicate divergence, the work around is to use a modal
analysis and expand the time steps in terms of selected eigenvectors or residual error method
solved by Gauss Newton solution technique. A change in the strain measure, by going from
71
72 CHAPTER 8. FUTURE WORK
the small strain to the Cauchy-Green measure is important in order to raise the code from
master thesis level into research level, especially since small stress tensor requires smaller
deformation than 2% in order to be of scientic value. The issue with the observed twisting
of the cantilever in the verication of the solid solver have its origin from the shear and
volumetric locking phenomena, and hence the solution to this is using under-integration with
a small eort change in creating the stiness matrix [19].
A major improvement in the key associative mapping can be implemented by storing the dis-
tance mapping from previous run and by threshold assign the mapping function and thereby
change the order for
sf
from O(N M) to O(N) where N is the solid center face points and
M the uid and vice versa O(M) for
fs
. The rst version using quasi-static approximation
failed due to improper settings in case, the question is whether this is suitable for cases with
KC 5 as proposed by [5], the code have the benet of allowing arbitrary deformations and
with a proper viscous material modelling it is assumed that modeling exible structures as
rubber materials is a feasible task. The adaptive meshing is an important task for making the
FSI useful for VIV studies. The error lies in the wake and the vortices. By rening the uid
mesh using icoErrorEstimate has been a successful task according to OpenFOAM forum.
Several reports indicates the signicant gain by this procedure in accuracy to given mesh
and thus allowing coarser mesh and in turn faster simulations. A question mark must be set
regarding the oscillation error , whether it is a numerical or a physical consequence of an
inadequate modeling of the solid. The current BC allows only natural Dirichlet condition on
a stationary boundary. There is however a small task to allow the boundary to be attached
to a movable boundary. Finally, allowing turbulence is an important extension in order to
fully estimate bueting, this again a smaller task, a modication in the traction vector.
A complementary literature study revealed the presence of few master thesis on the same
subject, one of particular interest is the one written by Michael St ockli [49] and PhD report
with code on FeniSC, where the simplicity of the code is an appealing feature and the ques-
tion arise whether future development favor this package, however, material modelling often
require signicant change in the inner structure of the assemblage, the question is if this can
be a doable task. Further, if MPI cannot be implemented on this platform the actual bene-
ciary feature of DEAL.II is lessened and if also adaptive meshing was feasible in FeniSC then
a change is an obvious consideration since current version contain only a non-parallelized
adaptive mesh solver[73].
Bibliography
[1] Alexander Aitken, On Bernoullis numerical solution of algebraic equations, Proceed-
ings of the Royal Society of Edinburgh 46 (1926) 289-305.
[2] M.Amabili, Reduced order P.O.D Models for nonlinear vibrations of cylindrical shells
with F.S.I.
http://uid.ippt.gov.pl/ictam04/text/sessions/docs/FSM4/12742/FSM4 12742.pdf
[3] K.J.Bathe , H.Zhang and M.H.Wang, Finite Element Analysis of incompressible and
compressible uid ows with free surfaces and structural interactions. Computers and
Structures 56(2-3) (1995) 193-213 .
[4] Klaus-J urgen Bathe and Hou Zhang, Shanhong Ji, Finite element analysis of ows fully
coupled with structural interactions, Computers and Structures 72 (1999) 1-16.
[5] P.W. Bearman, J.M.R Graham and E.D Obasaju, A model for transverse forces on
cylinders in oscillatory ows. Applied Ocean Research,6(3) (1984) 166-172.
[6] Peter Bearman and M asa Brankovic, Experimental studies of passive control of vortex-
induced vibrations. Europena Journal of Mechanics B/Fluids (2004) 9-15.
[7] K.Yusuf Billah and Robert H.Scanlan , Resonance, Tacoma Narrow bridge failure and
undergraduate physics textbooks. Am. J. Phys. 59(2) (1991) 118-124 .
[8] Frederic J.Blom, A monolithical uid-structure-interaction algorithm applied to the
piston problem. J. Comput. Methods Appl. Mech. Engrg 167 (1998) 369-391.
[9] C.L. Cun ,F. Biolley, E. Fontaine, S.
`
Etienne and M.L. Facchinetti, Vortex-Induced
Vibrations of Risers. Theoretical, Numerical and Experiemental Investigation, Oil, Gas
and Tech, Rev, IFP. 57(1) (2002) 59-69 .
[10] P.Dasdia and S.No, Vortex-Induced Vibration of reinforced concrete chimney: in situ
experimental and numerical prevision. J. Wind Engineering and Industrial Aerodynam-
ics 74-76 (1998) 765-776.
[11] J. Donea Antonio Hueste, J.Ph. Ponthot and A. Rodriguez Ferran, Arbitrary Lagrangian-
Eulerian Methods, Encyclopedia of Computational Mechanics, volume 1 chapter 14.
(2004) John Wiley & Sons.
[12] Einstein, Albert ,The Foundations of the General Theory of Relativity, Annalen der
Physik 6(30) (1916) p.74
[13] Earl H. Dowell and Kenneth C Hall, Modelling of uid-structure interaction. Ann.
Rev Fluid. Mech. 33 (2001) 445-90.
73
74 BIBLIOGRAPHY
[14] Edda Eich-Soellner and Claus F uhrer, Numerical Methods in MultiBody Dynamics.
Teubner-Verlag Stuttgart. Course book for FMN110, my assignement is with own data
by similar context is found in this book.
[15] C.Farhat and M.Lesoinne, Two ecient staggered algoritm for the serial and parallel
solution of three-dimensional nonlinear transient aeroelastic problems. Comput. Meth.
Applied Mech. Engrg. 182 (2000) 499-515.
[16] Zhengkun FENG and Azzeione Soulaimani, Nonlinear aeroelastic modeling using a
reduced order model based on proper orthogonal decomposition, Proceedings of Pro-
ceedings of PVP2007. ASME. Pressure Vessel and Piping Division. July 20-26, 2007.
PVP2007-26006.
[17] J. H. Ferziger, M. Peric, Computational Methods for Fluid Dynamics, Springer, 3rd Ed.,
(2001).
[18] F.Flemming and C.H.K. Williamson, Vortex-induced vibrations of a pivoted cylinder.
J. Fluid. Mech. 522 (2005) 215-252.
[19] Magnus Fredriksson, Doctoral Thesis, Class of accurate low order nite elements. Divi-
sion of Solid Mechanics LTH (2006).
[20] A.L.C. Fujarra, C.Pesce, F.Flemming and C.H.K.Williamson, Vortex-induced vibration
of a exbible cantilever. J. Fluids and Structures 15 (2001) 651-658.
[21] R.D.Gabbai and H. Benaroya, An overview of modeling and experiments of Vortex-
Induced Vibrations of Circular Cylinders. Journal of Sounds and Vibration 282 (2005)
575-616.
[22] Jean-Fr ed eric Gerbeau , Marina Vidrascu , A Quasi-Newton Algorithm Based on a Re-
duced Model for Fluid-Structure Interaction Problems in Blood Flows. ESAIM: Math-
ematical Modeling and Numerical analysis, 137(4) (2003) 631-647.
[23] M.Geradin and D.Rixen, Mechanical vibrations, Theory and Application to Structural
Dynamics (1997) Second Edition WileyBlackwell.
[24] Daniel Green and William G.Unruh, The failure of the Tacoma Bridge: A physical
model. Am. J. Phys. 74 (8) 706-716 (2006).
[25] Daisuke Ishiara and Shinbou Yoshimura, A monolithic approach for the interaction of
incompressible viscous uid and an elastic body based on uid pressure Poisson equa-
tion. Int. J. Numer. Meth. Engng. 64 (2005) 167-203.
[26] Laural Jacquin, Aircraft trailing vortice: An introduction. C.R.Physique 6 (2005) 395-
398.
[27] H.Jasak and H.G. Weller, Application of the nitie volume method and unstructured
meshes to linear elasticity. Int. J. Numer. Meth. Engng. 48 (2000) 267-287.
[28] Hvroje Jasak, Doctoral Thesis, Error Analysis and Estimation for the nite Volume
Method with Application to Fluid Flows, Imperial college, University of London, 1996.
[29] Cristophe Kassiotis, Which strategy to move the mesh in the computational uid dy-
namic code in OpenFOAM
http://powerlab.fsb.hr/ped/kturbo/OpenFOAM/slides/.
BIBLIOGRAPHY 75
[30] D.W. Kelly, J.P. Gago, O.C. Zienkiewicz, I. Babuska, A posteriori error analysis and
adaptive procedure for practical engineering analysis, Int. J. Numer. Methods Eng. 19
(1983) 1593-1619.
[31] T.Kitagawa, Y.Fujino and K.Kimura, Eects of free-end condition on end-cell-induced
vibration. J. Fluid and Structures 13 (1999) 499-518.
[32] Raghavan A. Kumar, Chan-Hyun Sohn and Banglaore HL.Gowda, Passive control of
Vortex-induced vibrations: An overview. Recent Patents on Mechanical Engineering 1
(2008) 1-11.
[33] Ulrich K uttler, Wolfgang A. Wall. Fixed-point uid-structure interaction solvers with
dynamic relaxation. Comput.Mech. 43(1) (2008) 61-72.
[34] Fabian Peng Karrholm, Rhie-Chow interpolation in OpenFOAM, Appendix from Nu-
merical Modelling of Diesel Spray Injection and Turbulence Interaction at Chalmers
University.
http://www.tfd.chalmers.se/ hani/kurser/OS CFD 2007/rhiechow.pdf
[35] Lieve Lanoye, Challenges in realizing vascular uid-structure-interaction using Fluent
and Abaqus software.
http://sympos.elis.ugent.be/archive/symp2006/papers poster/
paper037 Lieve Lanoye.pdf
[36] Anders Logg, Automating Finite element method. Arch. Comput. Methods Eng. 14 (2006)
93-138.
[37] Hermann G.Matthies and J.Sterndorf, Partioned strong uid coupling algorithms for uid-
structure interaction. Computers and Structures 81 (2003) 805-812.
[38] Sanjay Mittal and Sauray Singh, Vortex-induced vibrations at subcritical Re.J. Fluid.Mech.
534 (2005) 185-194.
[39] Prioz Moradnia, Phd course report on IcoDymFoam,
http://www.tfd.chalmers.se/ hani/kurser/
OS CFD 2007/PiroozMoradnia/OpenFOAM-rapport.pdf
[40] Tikeswar Naik , Ellen K. Longmire, S.C. Mantell , Dynamic response of a cantilever in liquid
near a solid wall. Sensor and Actuators A 102 (2003) 240-254.
[41] J.M.T. Penrose,D.R. Hose, C.J. Staples, I.S. Hamill, I.P. Jones and D. Sweeney, Fluid-
structure interactions: coupling of CFD and FE. CAD-FEM User Meeting Internationale
FEM-Technologietage. Sept 2001.
[42] Prandtl and Tietjens, Applied Hydro and Aeromechanics, Figure 60, Dover, ISBN 0-486-60375-
X.
[43] M.Razzaq, Jhron and S.Turek, Numerical simulation of laminar incompressible uid-structure
interaction for elastic material with point constraints. PDF draft document from Institute
of Applied Mathematics, TU Dortmund, Germany, submitted for publication and accepted,
Springer, Berlin (2008).
[44] Vincent Rivola, Comparative Study of the CFD codes Mistral and OpenFOAM - Application
to uid-structure interaction.
http://powerlab.fsb.hr/ped/kturbo/OpenFOAM/slides/
76 BIBLIOGRAPHY
[45] Takaaki Sakai Koji Iwate, Masaki Morishita and Seiji Kitamura, Vortex-Induced Vibrations
of a cantilever circular cylinder in super critical Reynolds Number Flow and its suppression by
Structure Damping. JSME Series B 44(4) (2001) 712-720.
[46] Goran Sandberg, CALFEM: Acoustic and interface elements for structure-acoustic analysis.
(2001) Structural Mechanics LTH, Sweden.
[47] Marianne Shubov, Mathematical modeling and analysis of utter in long-span suspended
bridges and in blood vessel walls, J. Aerospace Engineering 17(2) (2004) p. 70
[48] Marianne Shubov, Mathematical modeling and analysis of utter in bending-torsion, coupling
beam, rotating blades and hard disk drives, J. Aerospace Engineering 17(2) (2004) p. 56
[49] Michael Stockli, A Unied Continuum Fluid-Structure Interaction Solver using an ALE nite
Element Method. An Investigation on how to simulate blood ow. Master of Science Thesis,
KTH, Stockholm, Sweden 2007.
[50] P.A.B. de Sampaio P.H.Hallak, Alvaro. L.G.A Coutinho and Mich elle S. Pfeil , A stabilized
nite element procedure for turbulent uid-structure interaction using adaptive time-space
renement. Int. J. Num. Mech. Fluids, 44 (2004) 673-693.
[51] Robert H.Scanlan,An observation on Low-Speed Aeroelasticity, Amer. Soc. 150th Anniver-
sary Civilengineer paper, Conference paper. http://cedb.asce.org/cgi, by Jerey S. Russell, (ed-
itor) Reston, VA: ASCE, 0-7844-0686-3, (2003), p. 401.
[52] C.Trusedell and K.R. Rajagopal, An Introduction to the Mechanics of Fluids, Birkhauser.
[53] Zeljko Tukovic Hrvoje Jasak, Updated Lagrangian nite volume solver for large deformation
response of elastic body. Transaction of Famena 30(2) (2007) 1-18
[54] Jan Vierendeels, Implicit Coupling of partitioned uid-structure interaction solvers using re-
duced order models. Comp. and Struct. 85, Issue 11-14, (2007) 970-976
[55] Stefan Wagert, Markus Dreier and Martin Hegner, Frequency shifts of cantilevers vibrating in
dierent various media. Appl. Phys. Letter. 69 (19), Nov 4. (1996) 2834-2836.
[56] H.G. Weller, G. Tabor, H.Jasak and C.Fureby, A tensorial approach to computational con-
tinuum mechanics using object-oriented techniques. Computers in Physics 12 No. 6 Nov/dec
(1998) p. 620 .
[57] C.H.K Williamson, Advances in our understand of vortex dynamics in blu body wakes. J.
Wind Engineering and Industral Aerodynamics 69-71 (1997) 3-32.
[58] C.H.K. Williamson and R. Govardhan, Vortex-induced vibrations. Ann. Rev.Fluid.Mech 36
(2004) 443-455.
[59] Li Zhou and Yajun Ge, Wind Tunnel Test for Vortex-Induced Vibration of Vehicle-Bridge
system section model, J. Braz. Soc. Mech. Sci. and Eng. 30(2) (2008) 110-117.
[60] OC Zienkiewicz, RL Taylor, JZ Zhu The Finite Element Method: Its Basis and Fundamen-
tals,(2005). Elsevier Ltd.
[61] The homepage http://www.dealii.org. A collection of papers in how using deal.II at
http://www.dealii.org/developer/reports/.
[62] http://www.comsol.com/products/multiphysics/.
[63] http://www.uent.com/solutions/aerospace/pdfs/nl559.pdf.
[64] http://www.ansys.com/solutions/fsi.asp.
BIBLIOGRAPHY 77
[65] OpenCFD Ltd produce OpenFOAM
TM
, the open source computational uid dynamics (CFD)
toolbox http://www.opencfd.co.uk/openfoam/, Recommended sublink for a collection of papers
and slides on the subject, http://powerlab.fsb.hr/ped/kturbo/OpenFOAM/slides/.
[66] http://femcodes.nscee.edu/src/Benchmark.pdf.
[67] http://www.lstc.com/.
[68] http://en.wikipedia.org/wiki, on search karman vortex street. Note however that is several dif-
ferent formulas for dierent Re they are however closely related and dierence lies in what
source used and experiments performed.
[69] http://www.cs.rpi.edu/szymansk/OOF90/F90 Objects.html.
[70] http://www.amath.washington.edu/lf/software/CompCPP F90SciOOP.html.
[71] http://www.tfd.chalmers.se/hani/kurser/OS CFD 2008/.
[72] http://en.wikipedia.org/wiki/Fluid-structure interaction.
[73] http://fenics.org.
[74] http://www.python.org.
[75] http://caelinux.com.
[ All links are update to the present date of presentation of the thesis. ]

Structural Analysis Using Openfoam

Uploaded by

Copyright:

Available Formats

Structural Analysis Using Openfoam

Uploaded by

Document Information

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Structural Analysis Using Openfoam

Uploaded by

Copyright:

Available Formats

Fluid-Structure Interaction (FSI) case study of a

cantilever using OpenFOAM and DEAL.II with

show a linear dependency except for a region

= l,the following dimensionless governing equation is obtained,

) = 0 is achieved by the algorithm [1],

), then f(x) = g(x) x becomes the residual

) is set to 400 and hence giving the kinematic viscosity

You might also like