07 Cachan

Fluid-Structure Interaction
Hermann G. Matthies
Institute of Scientic Computing
Technische Universitat Braunschweig
Brunswick, Germany
acknowledgements to: M. Meyer, R. Niekamp T. Srisupattarawanit, J. Steindorf
[email protected]
http://www.wire.tu-bs.de
2
Content
Some of our work on coupled problems,
organisational structures to nance doctoral students.
Partitioned Approach to Coupled Problems
Translation into Software Components
Some Future Research Directions
DFG nanced Graduiertenkolleg
TU Braunschweig Institute of Scientic Computing
3
Our Work on Fluid-Structure Interaction (FSI)
Let us start with an example, multi-physics,
here uid-structure interaction (FSI).
We will look at
Oshore Wind Turbines
Random Surface Waves
Coupling as abstract DAEs
Strong Coupling Algorithms
Software Components and FSI
4
Oshore Windturbine
Multi-Physics / Sub-systems
turbulent wind
blade element aerodynamics
dynamic stall
structure: blades, nacelle, tower
generator, controller
random wavesthree sub-zones
pile, soil
soil far-eld
5
Subsurface
Tower, pile, near-eld soil, and surrounding water
6
Wave Domains
Dierent wave domains to avoid reections
Domain 2
Domain 1
Absorbing zone
Domain 3
7
Subsystems of Oshore-Wind-Turbine I
Wind Stochastic wind model in frequency domain.
Wind/ Blade-Structure Blade element theory for the turbine blades.
Simple force and torque balance in decoupled rotor disk annuli.
Blade/ Nacelle & Tower Manages non-inertial rotating co-ordinate systems,
one for each blade, one each for nacelle, tower. Structure-structure
coupling. (Plus generator, control, etc.).
Random Waves/ Fluid-Fluid Coupling To prevent parasitic reected waves
and get the right random properties three wave domains.
1. Deep water, only linear incident waves2. Shallower water, nonlinear
incident and linearised refracted waves3. Shallow water, fully nonlinear
waves and FSI.
8
Subsystems of Oshore-Wind-Turbine II
Deep Water Linear Airy waves for deep water in spectral description.
Shallower Water Nonlinear incident wave, linear perturbation refracted wave.
Potential ow. FE-model for free surface, fast multipole BEM solver for
water below.
Shallow Water Fully nonlinear wave, both for incident and refracted wave.
Potential ow. FE-model for free surface, fast multipole BEM solver for
water below.
Wave/ Tower Potential FlowFree Surface/ Tower Coupling.
Tower/ Soil Coupling tower/ pile FE-model with a near eld soil FE-model.
Soil/ Soil Near eld FE-model of soil coupled with far-eld SBFEM code.
9
Some Results
Model to be used for Monte Carlo fatigue simulation
10
Some of the OWECS Computational Models
Stochastic / turbulent wind generated by spectral description through (fast)
Fourier transform.
Aeordynamics by blade element theory, simple mass owmomentum and
rotational momentum balance, and mass conservation in each rotor annulus
independently, together with measured / previously calculated prole data
(C
L
, C
D
, C
M
) to determine induced velocities.
Structure description is standard large displacement beam model, in
non-inertial (co-rotational) reference systems. Discretised with corresponding
beam nite element system.
Soil is linear elastic solid, near eld with nite elements, and far eld with
scaled boundary nite element method.
11
Free Surface Conditions I
WavesPotential inviscicd ow with free surface:
Fluid velocity v(x, t) = (x, t) with potential , and for all t:
2
= 0.
Free surface described by level set function F(x, t) = 0.
Eikonal equation for free surface
t
F/|F| movement has to match
normal velocity of uid v n with n = F/|F|:
v n =
F
|F|
!
=

t
F
|F|

dF
dt
:=
t
F +v F = 0.
An evolution equation for F:
t
F + F =
t
F
n
|F| = 0
12
Free Surface Conditions II
Next to this kinematic condition,
a dynamic condition comes from Bernoullis equation,
an evolution equation for on F = 0:
t
= g
1
2
||
2
+
p
a
with free surface elevation , gravitational acceleration g,

uid density , and exterior pressure p
a
.
Dynamics of wave domain purely on boundary,
inside domain Laplaces equation
2
= 0 is a constraint.
13
Fluid-Structure Coupling (FSI)
On uid-structure interface,
normal uid velocity
n
has to match structure velocity
Structure with displacement u and stress tensor
t
(u n) =
n
Next to this kinematic condition,

a dynamic condition comes from Bernoullis equation,
an evolution equation for (with initial depth d)
t
= g(d + +u
z
)
1
2
||
2
+
n + p
a

14
Formulation as PDAE (Water)
On free surface F = 0:
t
F =
n
|F|,
t
= g
1
2
||
2
+
p
a
.
On FSI interface:
n
=
t
(u n)
t
= g(d + +u
z
)
1
2
||
2
+
n + p
a
.
In uid domain:
2
= 0, on other bdrs:
n
= 0.
15
Formulation as PDAE (Structure)
On FSI interface:
t
(u n) =
n
,
n = (
n
p
a
) .
In structure:

2
t
u = div f,
= ConstitMod(u) . . . .
Plus Control, Soil, etc. = more PDAEs.
16
Time-Stepping Algorithms Used
In the structure and for dynamic stall Generalised- method.
On the wave boundary Adams-Bashforth-Moulton
3-step predictor-corrector method.
In the soil Newmarks method.
Global time-stepping and iteration (control)
may each be implemented as separate software components.
17
Coupled Problems
Coupled problems often combine the models of two or more physical systems,
they are multi-physics modells. Dierent approaches:
Monolithical approach, which means one global model of everything
Advantages: All encompassing theoretical and numerical treatment
Disadvantages: Treatment is ever more complex, every time completely
new start, new algorithms, new software, does not scale, not modular
Partitioned approach, which means separate models plus coupling
Advantages: Complexity constrained to one physical domain, theory may
be well worked out, ecient numerical algorithms, existing sophisticated
software for each subsystem, modular, scalable
Disadvantages: Subsystems have to be coupled together, new numerical
and algorithmic problems, coupling software necessary
18
Dierential-Algebraic Coupling
For the sake of simplicity, assume two subsystems:
Assume that subsystems are dierential-algebraic equations (DAEs) with local
dierential variables x
1
X
1
, local algebraic variables y
1
Y
1
:
x
1
= f
1
(x
1
, y
1
, z) ,
0 = g
1
(x
1
, x
2
, y
1
, z) ,
same for the second subsystem with local dierential variables x
2
X
2
, local
algebraic variables y
2
Y
2
:
x
2
= f
2
(x
2
, y
2
, z) ,
0 = g
2
(x
2
, x
1
, y
2
, z) .
Coupling conditions formulated as algebraic constraints with global
algebraic variables z Z:
0 = h(x
1
, x
2
, y
1
, y
2
, z) .
19
Dierential-Algebraic Regularity
Assume that each single subsystem, and also global system is an index-1 DAE.
This means that the operator matrices
D
y
j
g
j
,
_
D
y
j
g
j
D
z
g
j
D
y
j
h D
z
h
_
, (j = 1, 2) ,
_
D
y
g D
z
g
D
y
h D
z
h
_
have to be regular, where D
q
is the partial Frechet derivative w.r.t. q,
and we have set g = (g
1
, g
2
)
T
and y = (y
1
, y
2
)
T
.
After time discretisation of subsystems, dierent algorithms may be used in
partitioned approach for total system,
e.g. simple switching, or basic staggering (weak coupling),
or globally implicit methods with equilibrium iteration (strong coupling),
e.g. by (nonlinear, block) Jacobi, Gauss-Seidel,
or Newton variants for better convergence.
20
General Fluid-Structure-Coupling (FSI)
Fluidincompressible Newtonian uid (i.e. Navier-Stokes eqns.) in ALE
(Arbitrary Lagrangean-Eulerian) formulation: (plus boundary conditions)
in
f
:
f
( v + (v ) v) div +p = r
f
,
2 = (v + (v)
T
) = 2
s
v, div v = 0,
Solidlarge deformation elastic St. Venant material in Lagrangean
formulation: (plus boundary conditions)
in
s
:
s
u DIV(FS) = r
s
, F = I + GRADu
S = (tr E)I + 2E, 2E = (C I), C = F
T
F,
Arbitrary Lagrangean-Eulerian coordinate system:
in
f
: L =
u
21
FSI Interface
Conditions on interface
I
between
f
and
s
:
At spatial location (t) =
0
+ u(
0
, t)
I
continuity of velocities:
v((t), t) = u(
0
, t).
Variational formulation for velocity condition:
_
I
(v((t), t) u(
0
, t)) d
I
= 0
Conservation of momentumbalance of tractions:
( pI) n =
1
J
FSF
T
n , J = det F .
Variational formulation for traction conditiontreat like any other boundary
traction, boundary traction is equal to Lagrange multiplier
I
.
22
Discrete Form of Fluid-Structure-Interaction
The Fluid and the ALE-Domain:
M
f
v +N(v )v +K
f
v +B
f
p = r
f
+T
T
f
I
; ,
B
T
f
v = 0 ,
K
g
= Au .
The Solid:
M
s
u +K
s
(u)u = r
s
T
T
s

I
.
The Interface:
T
f
v = T
s
u .
Equality of interface tractions on coupling interface already included in terms
T
T
f
I
and T
T
s

I
.
23
The DAE Correspondence I
x
1
= f
1
(x
1
, y
1
, z) ,
0 = g
1
(x
1
, x
2
, y
1
, z) ,
The Fluid and the ALE-Domain: (with b := v, the uid accel.) and = :
x
1
:=
_
v
_
, y
1
:=
_
_
b
p
_
, f
1
:=
_
b
_
, z :=
I
,
g
1
:=
_
_
M
f
b +N
f
(v )v +K
f
v +B
f
p r
f
T
T
f
I
B
T
f
M
1
f
(N
f
(v )v B
f
p K
f
v + r
f
+T
T
f
I
)
K
g
Aw
_
_
,
24
The DAE Correspondence II
x
2
= f
2
(x
2
, y
2
, z) ,
0 = g
2
(x
2
, x
1
, y
2
, z) .
The Solid: (with a := w = u, the structural acceleration):
x
2
:=
_
u
w
_
, y
2
:= a , f
2
:=
_
w
a
_
, z :=
I
g
2
:= M
s
a +K
s
(u)u r
s
+T
T
s

I
; ,
0 = h(x
1
, x
2
, y
1
, y
2
, z) .
The Interface: h := T
f
b T
s
a.
25
Use Time-Stepping instead of Iteration
Sought solution x
may be seen as equilibrium point of a dynamical system

with state x
(n)
at time levels n t with the time advance operator
x
(n+1)
= (x
(n)
).
For a convergent iteration the equilibrium x
is an asymptotically stable xed point.

This is actually a method to compute the steady state or equilibrium,
often used in CFD (sometimes called dynamic relaxation).
An iteration method is an abstract dynamical system,
each iteration x
k+1
= (x
k
) corresponds to one time step.
Iteration converges i |
j
| < 1 for all eigenvalues of D(x
).
26
Time-Stepping and Iteration
Time stepping scheme for x
(n+1)
= (x
(n+1)
, x
(n)
) is a
discrete dynamical system x
(n+1)
:= (x
(n+1)
1
).
Convergence equivalent with stability of time-stepping scheme for given t.
Convergence of this iteration and stability of time-stepping may be
investigated with same theory.
Weak coupling/ simple switching corresponds (in simple case is equal to)
strong coupling/ block-Jacobi iteration.
Weak coupling/ basic staggering corresponds (in simple case is equal to)
strong coupling/ block-Gauss-Seidel iteration.
27
Global Equations for Strong Coupling
Coupling condition h = 0 usually adjoined to one subsystem.
Set := (x
1
, y
1
, z)
T
= (v, , b, p, ,
I
)
T
and := (x
2
, y
2
)
T
= (u, w, a)
T
to include interface in rst equation,
otherwise include z =
I
in and not in .
Assume that convergent iterative solvers for subsystems exist:
= F
1
(
1
, ) , and
= F
2
(
1
, ) , = 1, 2, . . . ;
Simplest solution process is nonlinear block-Jacobi,
an additive or parallel Schwarz procedure (corresponds to simple switching):
= F
1
1
(
1
,
1
) ,
= F
2
2
(
1
,
1
) ;
28
Nonlinear block-Gauss-Seidel
Almost as simple is nonlinear block-Gauss-Seidel,
a multiplicative or serial Schwarz procedure (corresponds to basic staggering):
= F
1
1
(
1
,
1
) , and with new
, do
= F
2
2
(
1
,
) .
Theorem:[Arnold, G unther] In block-G-S, let L be Lipschitz-constant of
j
,
and let
= max
t[0,T]
(D
y
2
g
2
)
1
D
z
g
2
_
D
y
1
h(D
y
1
g
1
)
1
D
z
g
1
_
1
D
y
2
h ,
the iteration only converges if < 1, and if at least iterations are performed
so that L
< 1, and the total global time-step error is bounded by

< C(
max{0,2}
(x) +
1
(y)) +
1
(x) +
2
(y),
is extrapolation error,
j
is subsystem integrator error, and
= + O(t) < 1.
29
Dierent Possibilities for block-Gauss-Seidel
For dierent ordering and distribution of constraint, we have
1st uid plus coupling, 2nd solid: = M
1
s
T
T
s
(T
f
M
1
T
T
f
)
1
T
s
,,
where

M
1
= M
1
f
(M
f
B
f
M
p
B
T
f
)M
1
f
is a Schur complement, and
M
p
= (B
T
f
M
1
f
B
f
)
1
. Note
f
/
s
.
1st structure plus coupling, 2nd uid: =
M
1
T
T
f
(T
s
M
1
s
T
T
s
)
1
T
f
.
Note
s
/
f
.
1st uid, 2nd solid plus coupling: = (T
s
M
1
s
T
T
s
)
1
(T
f
M
1
T
T
f
) .
Note
s
/
f
.
1st solid plus coupling, 2nd uid: = (T
f
M
1
T
T
f
)
1
(T
s
M
1
s
T
T
s
) .
Note
f
/
s
.
depends on ratio of
f
and
s
.
30
Block-Newton
Desirable is an iteration scheme which will not depend on
ordering and distribution of constraint: Block-Newton.
In each block-Newton iteration following system has to be solved:
_
I D
F
1
D
F
1
D
F
2
I D
F
2
_ _
_
=
_
F
1
(
F
2
(
)
_
.
Symbolic block-Gauss elimination:
= (I D
F
1
)
1
( F
1
(, )) C ,
with the multiplier matrix C := (I D
F
1
)
1
[D
F
1
].
Further with Schur complement matrix S:
S := (I [D
F
2
] [D
F
2
]C) = r,
with r := ( F
2
(, )) + [D
F
2
]q, q := (I D
F
1
)
1
( F
1
(, )).
31
Solving the Block-Newton System
Solution proceeds by Krylov method (Bi-CGstab):
Solving a system with (I D
F
1
): Apply iterative solver F
1
Same with C, plus nite dierences for [D
F
1
]
Solving the Schur-complement system:
Use Bi-CGstab.
Compute r with iterating subsystem solver F
2
;
compute action of S by nite dierences.
Theorem:[Mackens, Voss] If the single system solvers are quadratically
convergent (or enough iterations are made in the approximative steps), the
global iteration is also quadratically convergent.
32
Quasi-Newton
Quasi-Newton methods are generalisations of secant method:
H
_
=
_
F
1
(
F
2
(
)
_
.
Easy to solve with H
(explicit inverse H
1
).
H
changes by low rank only from step to step.

H
1
= H
1
1
+a
b
T
a rank one updateora rank two update

H
1
= H
1
1
+a
a
T
+b
b
T
, b
are easy to compute from known data

33
A Simple Example
Fester Rand
Fester Rand
Elastische Struktur
Starrer Krper
E
i
n
f
l
u
r
a
n
d
A
u
s
f
l
u
r
a
n
d
1
.
0
1
2
.
0
5.5 14.0
4.0
34
Movement and Pressure Distribution
35
Tip Displacement Response Weak Coupling
1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8
4
3
2
1
0
1
2
3
Simulationszeit t: dt = 0.02
y
V
e
r
s
c
h
ie
b
u
n
g
e
n

a
m

S
t
r
u
k
t
u
r
e
n
d
e
36
Tip Displacement Response Strong Coupling
1 1.5 2 2.5 3 3.5 4
1.5
1
0.5
0
0.5
1
1.5
y
V
e
r
s
c
h
ie
b
u
n
g
e
n

a
m

S
t
r
u
k
t
u
r
e
n
d
e
37
Iteration Count
1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2
0
2
4
6
8
10
12
14
16
18
20
I
t
e
r
a
t
io
n
s
z
a
h
l
Approximatives BlockNewton
BlockGauSeidel
38
Solver Calls
1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2
50
55
60
65
70
75
80
85
90
A
n
z
a
h
l
v
o
n

L
s
e
r
a
u
f
r
u
f
e
n
Approximatives BlockNewton
BlockGauSeidel
39
Another Example
Material St.Venant: = 0.2, E = E
l
Navier-Lame: = 0.2, E = E
r
40
Iteration Count
E
l
= 10
5
, E
r
= 5 10
4
E
l
= 10
5
, E
r
= 10
6
E
l
= 5 10
4
, E
r
= 10
6
iter cpu[t]
Jacobi 18 10.9
Gauss-Seidel 14 7.9
BFGS 12 6.8
Newton 3 13.5
iter cpu[t]
Jacobi -
GS 28 16.0
BFGS 12 4.8
Newton 3 13.7
iter cpu[t]
Jacobi -
GS -
BFGS 25 13.4
Newton 12 52.0
41
Software Components
Back to Oshore Wind Turbine
turbulent winda Matlab routine
aerodynamics, dynamic stall, structure components (blade, nacelle,
tower)separate programs coded in C
stochastic waves as a Matlab routine
each wave domain separate, Matlab for boundary treatment, for Laplaces
equation a fast multipole solver coded in C
soil separate in near and far eld, as C and Fortran codes
42
Some of the Components
Fluid flow
NWT code
Matlab
Soil
Felt code
C
Unbounded Soil
Fortan
Similar code
Laplace Solver
FastLap code
Fortan
Structure
and
WKA code
Matlab
Wind
43
Coupling of Components
Coupling of components is here via Component Template Library (CTL),
software written in C++.
Coupling is via remote method invocation (RMI)like subroutine calls.
No assembler-like message passing book-keeping.
CTL has bindings with C++ (naturally),
and with C, Fortran, Java, Python, Matlab.
Component properties specied in interface denitionlike
an external object in the sense OOP.
CTL supports the following linkage types/ communication protocols:
dynamic linkage, threads, pipes and demons, MPI (message passing interface),
PVM (parallel virtual machine), TCP/IP (via sockets), and some others.
44
Conclusions for FSI
Coupled systems in a partitioned approach
Translates naturally into software components
Coupling may introduce algebraic constraints DAEs
Time-stepping, global iteration, algorithmic control (for strong coupling)
may be software components.
CTL allows naturally looking (and feeling) coupling / binding
Software download and more info at
http://www.wire.tu-bs.de/forschung/projekte/ctl/e ctl.html
45
Some Future Research Directions
Numerical Methods
for breaking surface waves
by coupling BEM-NWT with particle methods
for transition to molecular dynamics (MD)
e.g. ows in capillaries with embedded particles
for uid-particle ows
like debris, mud, sand (coupling owobjects).
46
DFG nanced Graduiertenkolleg
DFG Deutsche Forschungsgemeinschaft
( = German Science Foundation )
nances for up to 9 years
so-called Graduiertenkollege = doctoral research groups
http://www.dfg.de/en/research funding/ ...
coordinated programmes/research training groups/index.html
Intention isbroadly speakingto allow research on a high level
without any teaching or project obligations
At the basis is a common research area
preferably interdisciplinary
carried by a core group of 510
well established researchers.
47
DFG-Graduiertenkolleg Funds
Funds
for stipends for 1015 doctoral students
and 12 postdocs (max. 3 years)
for their travel and conference costs
for short stays abroad
for inviting guests and cooperating researchers
for guest lectureres.
48
DFG-Graduiertenkolleg Study Programme
Funds
for specic study courses
for compact courses
for workshops
for inviting guest-lecturers
for training courses
for co-ordination.
49
DFG-Graduiertenkolleg FSI at TU BS
established October 1998
rst evaluation after 3 years
extended to maximum duration
until September 2007
interdisciplinaryMathematics and Informatics,
Civil Engineering,
and Mechanical Engineering,
as well as DLR (=German AeroSpace Res. Est.) at Braunschweig
http://www.tu-braunschweig.de/gk432/
50
FSI-Subjects at TU BS
Innovative Numerical Methods
e.g. particle methods, model reduction.
Aero-Structure-Interaction
e.g. wind and structures, aero-elasticity of lifting devices,
aero-elastic-optimal design, ight of birds.
Water-Structure-Interaction
e.g. earthquakes and uid tanks, structures in owing water,
waves on structures/ dams, swimming of shes.
Acoustics and FSI
e.g. mainly linearised models of FSI.
51
At Last...
Thank you for your attention !
[email protected]
http://www.wire.tu-bs.de

07 Cachan

Uploaded by

Copyright:

Available Formats

07 Cachan

Uploaded by

Document Information

Original Title

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

07 Cachan

Uploaded by

Copyright:

Available Formats

Fluid-Structure Interaction

with free surface elevation , gravitational acceleration g,

Next to this kinematic condition,

TU Braunschweig Institute of Scientic Computing

may be seen as equilibrium point of a dynamical system

is an asymptotically stable xed point.

< 1, and the total global time-step error is bounded by

changes by low rank only from step to step.

a rank one updateora rank two update

are easy to compute from known data

You might also like