7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Development of Numerical Scheme on FluidStructure Interaction for a Flexible
Object
K. Tsujimoto*, A. Kito*, T. Shakouchi* and T. Ando*
Division of Mechanical Engineering ,Graduate School of Engineering, Mie University
1577 Kurimamachiycho, Tsu, 5148507, Japan
tujimoto@mach.mieu.ac.jp
Keywords: FSI, Flexible object, Immersed boundary method
Abstract
Numerical solving of fluidstructure interaction (FSI) problems are of importance in various scientific fields. In par
ticular, the establishment of stable numerical scheme is needed in order to analyze the details of flow phenomena. In
the present paper, we propose a weakcoupling method for FSI problem: for the flexible object the rigorous equations
of motion are discretized with finite volume method; in the flow computation, the flexible object is reproduced via
immersed boundary method. To demonstrate the performance of proposed scheme, the 3D structure analysis and 3D
FSI problem for a flexible object are solved. The results indicate that the computation is stably conducted using the
proposed method, in spite of the occurrence of fairly large deformation of the object.
Introduction
FluidStructure Interaction (FSI) problem is concerned
with in various research fields such as mechanical,
aerospace, civil and medical engineering. The highly
accurate and stable scheme for the FSI is desired to be
developed. In the view of numerical difficulty, while
the FSI for the structure of rigid body is easiest prob
lem, the large deformation of structure often induces
the divergence of computation due to the distorted grid.
Namely depending on the rate of deformation, the treat
ment of FSI problem is difficult. In general, the major
numerical method for FSI problem is divide into two cat
egory, i.e., strongcoupling method[l] (or monotholic
method[2,3] ) and weakcoupling method[4,5]. As the
merit of strong coupling it has been already mentioned
that the governing equations of both fluid and struc
ture are simultaneously solved, resulting in the stable
and highaccurate computation. On the other hand in
the weakcoupling method, since the both computation
of fluid and structure are performed independently of
each other, the performance for the stability and accu
racy is inferior than that of strongcoupling. However,
the conventional code for the structure and flow compu
tation can be utilized through the slightly modification in
which the information is alternatively exchange between
flow and structure. Recently, since IB method[6,7] in the
flow computation is capable of dealing with any mov
ing structures, a keen issue needing to be solved in the
FSI problem including the large deformation structures,
concentrates to the development of the stable structure
computation. In the present study, we pay attention to
develop the stable scheme for structure computation,
and propose the weakcoupling method in which for
the flexible object the rigorous equations of motion are
discretized with finite volume method (FVM); for the
flow computation, the finite difference method (FDM)is
used and the flexible object is reproduced via immersed
boundary method[7].
Finite volume method for elastic body
Governing equations and discretization
The dynamics of elastic body based on the continuum
model can be written by the conservation law of mass
and momentum:
dp 9ui
P% (1)
dt iaxi
dui Tij (2)
dt Oxj
where p is the density, ui is velocity component.
In the present study, constitutive equations are as
sumed to be a Hookean elastic body:
isj = Sij i,' I (3)
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
di C( (t r()
dt 3
gij = 9 + }ui
K tr(i)
(4) The governing equations are discretized using the hexa
hedron element proposed in a classical FVM[8]. The
velocity and the displacement are defined at the grid
points constructing the hexahedron element, and the
(5) stress is at the center of the element. The temporal dis
cretization is as follows[8];
n+1/2 Ouin n1/2
(6) Ei2 9x
(+ (O )
9x
1/2
At (15)
where sij is the stress deviator tensor, Eij strain tensor
and iij the rate of strain. G is the modulus of elasticity,
K the bulk modulus, / artificial viscosity tensor and
p pressure. tr(i) is the trace of the tensor r:
tr () ( i + + ak) (7)
G and K are written by the Young's modulus E and the
Poisson ratio 7:
E
G 2(1 + K
2(1 + ) '
E
3(1 2y)
In order to consider the large deformation of elastic
dJsij dJsij .
body, the rate of stress d is introduced. d is
called as the Jaumann derivative which means the rate
of stress tensor in a reference frame being corrected ac
cording to the rotational effect:
S= sij + SikQkj + SkjQik
dt
where (2 is an antisymetric rotation tensor:
Qij 09ujdz
(aai
n+l s
+ GE 1/
ot1
1tr 612 At
3 ( )
[Sk2j+ k] At
n+1/2
n+l
x
n+1/2
n 1
n1 2 d+ (Ui
'A t
z dt
xn+ n )
Stabilizing the grid distribution due to hourglass
mode
Among of the numerical difficulties, the hourglass
mode sometimes occurs in the quadrilateral grid as
shown in Fig.1. This mode is called as zeroenergy
mode, the grid continues to deform such that the hour
glass pattern appears. In order to avoid the occurrence of
this mode, we investigate the smoothing filter proposed
by SavitzkyGolay [9](hereinafter SG filter).
O9ui
09x )
If the governing equations are normalized with po, G
and the reference length L, the physical quantity are as
follows:
S Ui P t
ui= P* = L ,'
JG/po GC L ^po/GC
Thus, the nondimensionalized equations are:
dp* p*
dt* xi*
p, du O i8 *
dt* 9x*
Ts = s*i p* i q*
L Figure 1: Hourglass mode
(11)
In order to evaluate the usefulness of SG filter, the de
formation of a cantilever are examined using the above
mentioned FVM. In the case without SG filter, as shown
(12) in fig 2(a), the hourglass mode occures, and then the
computation does not continue due to the numerical di
vergence. On the other hand, in the case with SG filter
(13) the hourglass mode is completely suppressed as shown
in Fig 2(b). Hereinafter, the SG filter is introduced in the
all computaion.
Dynamic of threedimensional cantilever
(14) In order to evaluate the validation of simulation code,
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
K:
AI
(a) without filter
(b) with filter
Figure 2: The effect of SG filter for hourglass mode
t
0(
(a) = 0.2
005
(b) 7 0.3
5
2.0L(n,=22)
.OL(n, =52)
L =
L (n =12)
z Z
Free
Force
005
0 05
(c) 7 0.4
Fixed
Figure 3: Computational domain and coordinate system
vibration analysis of a fixedfree cantilever is examined
as similar to the previous papers[10,11]. In the present
study, the shape of cantilever shown in Fig.3 is dealt
with. The dimension of cantilever is 5L x 2L x L, and
the grid size is n, x x n, 52 x 22 x 12. The
end of cantilever is fixed at x 0 and the another end
is free. The free end of cantilever is displaced 1/20L in
z direction at initial instant, and then the cantilever con
tinues to freely oscillate. The Poisson ratio y varies 0.2,
0,3 and 0.4.
Figures 4 show the time evolution of the tip displace
ment of cantilever in z direction. Since the Hookean
elastic model is assumed, the amplitude of oscillation is
constant as well as the periods of oscillation. Table 1
shows the comparison of the theoretical value and the
Figure 4: Time evolution of tip displacement at different
Poisson ratio
present results as for the natural frequency The theoreti
cal value, f+ is as follows[12]:.
A 1 2 El
27 L2 pA
p; density ,E; Young's modulus
I; moment of inertia in the plane or second moment of
area
L; length of cantilever
A; cross sectional area of cantilever A; Eigen value
In this calculation, A = 3.52 is selected as the first
bending mode[12]. From Table 1, the computational re
sults are good agreement with the theoretical value sug
gesting the computational code is established.
The effect of artificial viscosity
In case of the solid dynamics, finescale numerical os
cillation sometimes occures. In general, to avoid the
r'
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Table 1: Comparison of computational value with theo
retical value
Frequency Error
Theoretical Value (y 0.2) 0.01110
Computational value (7 = 0.2) 0.01070 3.8%
Theoretical Value (y 0.3) 0.01120
Computational value (7y 0.3) 0.01155 3.1%
Theoretical Value (7 = 0.4) 0.01199
Computational value (7 = 0.4) 0.01175 2.0%
01 4
01
(a) CNS = 0
(b) CNS = 0.05
 C 0O
~7)~ jjj( fj
C= 50
Figure 5: Time evolution of tip displacement at different
artificial viscosity
unphysical phenomena, the artificial viscosity is intro
duced. The artificial viscosity tensor is defined as fol
lows[8]:
n 1/2
qij
G, +1/2 1tr (.n 1/2) (21)
G i i; 3 (r 6;y I2)
where Gi is constant.
Gi= CNS(PO) v]
In order to investigate the effect of artificial viscosity,
the coefficient CNS varies 0, 0.05, 0.5 and 5. As the test
case, the previous mentioned cantilever shown in Fig.3
is examined. Figure 5 shows the evolution of the tip dis
placement of cantilever in z direction. At CNS =0 and
0.05, the similar behavior of oscillation is observed, in
dicating that artificial viscosity does not contaminate the
(c) CNS = 0.5
(d) CNS = 5.0
Figure 6: Distribution of normal stress in x direction at
different artificial viscosity
numerical accuracy. At CNS 0.5, the amplitude of
oscillation gradually decreases. At CNS = 5.0, after the
start of simulation, the amplitude rapidly decreases due
to the strong artificial viscosity. Next, the large deforma
tion occures by adding the strong force. Figures 6 show
the distribution of normal stress ,, along the x axis at
z = 0 and y = L. Although at CNS = 0.05, 0.5 and
5.0, the stress smoothly distributes, but the unphysical
oscillation occures at CNS = 0 Thus above mentioned
results suggest that CNS = 0.05 is appropriate to com
pletely suppress the numerical instability.
Analysis of large deformation problem
In order to investigate whether the improved FVM
scheme is able to reproduce the behavior of large defor
mation, the cantilever is deformed by adding the strong
external force. As the test case, previous mentioned
cantilever shown in Fig.3 is deformed by a strong force
adding at the free end of cantilever. Figures 7 show the
surface distribution of normal stress ,xx. From figures,
it is found the fairly large deformation are reproduced
through the stable computation. In the figures, the red
and blue color correspond to the tensile and the compres
sion stress, respectively. In particular, when the large
deformation occures in Fig. 7(c), the strong stress dis
tribute around the root of cantilever, demonstrating the
qualitatively correct. Further when any external force
is imposed, in order to confirm that the improved FVM
scheme is able to reproduce the behavior of large defor
mation, another external force is added at the free end
of cantilever shown in Fig. 8. Figures 9 shows the
surface distribution of normal stress T,,. The behav
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
(a) t
(a) t*=0
(c) t* 2
(e) t* 4
(b)t* 1
(d) t* 3
(f) t* 5
 I
Em
(a) t*= 0
(c) t* 2
(e) t* 4
Figure 7: Contours of normal stress a,, for oneside
forcing case
L (n12) Force
50L(n 52)
Fixed
z L(n. 12)
Figure 8: Computational domain and coordinate system
ior of structure seems to be more complicated compared
to the Fig.7. Also, although the structure does not al
ways moves along the direction of computational grid,
the structure computation is stably conducted.
Numerical scheme for FSI problem
Governing equations and discretization
In the present study, we propose the weakcoupling
method for FSI. So the numerical scheme should be pre
pared for each of structure and flow field. For the flow
computation, the governing equation for unsteady in
compressible viscous flow in Cartesian coordinates sys
Figure 9: Contours of
forcing case
tem are as follows:
+u*
a9t* 0 ,* 3
^^ u=
normal stress ac, for twoside
Dp* 1 C2,,*
0x4 Re x*2
where, Re UoL/v is the Reynolds number. Uo is
reference velocity, v kinematic viscosity. p* is normal
ized with pUs. The momentum equations are discretized
using the CrankNicolson method for the viscous term
and the secondorder AdamsBashforth method for con
vective terms. The discretization is performed with the
fractionalstep method[13]:
u, .'
At
3 Hi
3
2 l
1  V + V2 (i + ,,P24)
2 2Re
where the superscript '1' means the index of time step,
and ui is the intermediate velocity component. Hi repre
sents the convective terms, Introducing Aa, = .a .
the momentum equations can be written as:
(1 At l 2
2Re }
2 1 V2+ 1 (25)
2 Hil R
/
(b) t* 1
I
(f) t* 
(f) t* 5
A{3 1
At Hi
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Structure gnd
SFlud gnd
it*
Qu Ui.+
(a)t*= 0
Figure 10: Approach of velocity near Immersed Bound
ary Method
Sp,
Structure grid
Q
C'
Fluid grid
_
P, I ,
(c)t*= 2
Figure 11: Interpolation of fluid pressure
a1+1 _ = _AtV1l+l
To satisfy the continuity equation, V ..1 = 0 should
be established, thus
V (9 1)
V (v0+ 1
acAtV (V1+1)
1
1lAV a,
alAt
.1+1 = Gi actAtV1 1
p11 11 lAt v2+ 1
2Re
Th spatial discretization is performed with a second
order central difference scheme. The staggered grid sys
tem is employed.
In the flow computaion, the interface between fluid
and structure is tracked using immersed boundary(IB)
method. In the IB method[6], to represent the boundary
on the Eulerian grid, the external force is imposed at the
boundary in the flow field so that the boundary velocity
is a specified value. The original concept was proposed
by Paskin[6]. However, it is wellknown the shortcom
ing, in which the time step should be considerable small
Figure 12: Velocity vectors around the cantilever
compared to the usual one being decided by a CFL con
dition[7].For this problem, Fadulm[7]proposed the im
proved IB method called as 'Direct forcing method'. As
shown Fig. 10, in the direct forcing method, the velocity
at the grid point in the vicinity of the wall, . is approx
imated with both the wall velocity, Uwall and the veloc
ity at the slightly distance from the wall, ui+ 1. The
approximated relation obtained for the velocity is incor
porated into the implicit treatment of viscous diffusion
term[7].
As shown in Fig. 11, the pressure at the structure
grid, pwai is linearly interpolated from the pressure at
the fluid computational grid located around the struc
ture grid. For the structure computation, the external
force for structure is estimated using Pwaii. The above
(b)t* 1
(d)t* 3
I I I
(a)t* 1
p
(b)t* 1
(c)t* 2
(d)t* 3
Figure 13: Velocity vectors around the cantilever
mentioned procedure is summarized as follows:.
1. The flow computaion is conducted using finite dif
ference method to satisfy the interface boundary
condition using IB method (Fig. 10)
2. The pressure at the interface obtained from flow
computaion is interpolated from the flow field to
the structure grid. (Fig. 11)
3. The structure computation is conducted using finite
volume method to satisfy the interface boundary
condition.
Through the step 1 ~ 3 are sequentially repeated, FSI
problem will be solved.
Twodimensional FSI problem
In order to check the IB method introduced in the flow
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Hz2 OL (n,=100)
Hy=0 5L (ny=50
Structure
Hx2 OL (nx=100)
2L/3
Figure 14: Computational domain for fluidstructure in
teraction
computation, the flow around the cantilever which is as
sumed to be a rigid body is simulated. The wall condi
tion is enforced on the both upper and lower boundary.
The cantilever is placed on the lower boundary. The uni
form flow, Uo is imposed on the left side boundary and
a convective outflow boundary is used on the right side
boundary. The grid size is n x = x 1 x 100 in x and
y directions. The dimension of computational domain is
Hx x Hy 2L x 2L. The Reynolds number UoL/v is
3000. The dimension of cantilever is L x 0.1L (the grid
size 30 x 20). Figures 12 show the instantaneous vector
plot and pressure contour. In figl2(b)(c), a large vor
tex is generated around the top of cantilever and moves
downward. The low pressure contour (back color region
) corresponding the vortex core is observed. Not shown
the close up, the velocity vector distribute along the sur
face of cantilever; the isocontour lines of pressure at the
structure surface penetrate perpendicular to the structure
surface. Thus the structure is correctly represented by
IB method. Next, the elastic body of cantilever is simu
lated. The grid size is nx x = 200 x 200 in x and
y directions. The other conditions is similar to that of
rigid case. The Reynolds number Re = UoL/v is 1000.
The Poisson ratio, 0.2, Young's modulus E* 2.4.
Figures 13 show the instantaneous vector plot and pres
sure contour. After the initial time, the cantilever gradu
ally inclines downstream; the vortex is generated around
the top of cantilever and moves away from the top of
cantilever while the structure continues to deform down
stream. The maximum bend occurs in Fig 12(c), then the
cantilever moves slightly back to the upstream. Com
pared to the Fig. 12, because of the large deformation
induces the drag reduction, the size of the issuing vortex
is smaller than that of the rigid case.
FSI for threedimensional problem
The computational volume is shown in Fig.14. The
wall condition is enforced on the both upper and lower
boundary. The cantilever is placed on the lower bound
ary. The uniform flow, Uo is imposed on the left side
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
(a) t*= 0
    :: :
(b)t* 1
(c) t* 2
(d) t* 3
Figure 15: Velocity vectors around the cantilever
boundary and a convective outflow boundary is used on
the right side boundary. Periodicity is enforced in span
wise direction The grid size is n, x, x n, 100 x 50 x
100 in x, y and z directions. The dimension of compu
tational domain is Hx x Hy x Hz 2L x 0.5L x 2L.The
Reynolds number is 3000. The dimension of cantilever
is L x 0.1L x 0.1L, (the grid size 30 x 10 x 10 ), the
Poisson ratio, 0.2, Young's modulus E* 2.4. Fig
ures 15 show the vector plot and 3D view of cantilever.
In the figure, red color denote the centerline through the
cantilever at the initial. It is observed that the cantilever
deforms downstream, at the same time oscillate in the
spanwise direction. Figure 16 show the distribution of
wallnormal vortices, w In the figure red color is posi
tive vorticity, blue color negative one. At the initial vor
ticity distribute around the cantilever, then the sheet like
vortex structure is formed. These structure is obviously
different from twodimensional simulation. Thus three
dimensional simulation is indispensable for FSI prob
(d) t* 3
Figure 16: Vorticity around the cantilever
lem.
Conclusions
We propose the weak coupling method: FVM for the
structure computation, FDM for the flow computation
and both scheme are connected by IB method. The
large deformed elastic body are simulated using the pro
posed scheme.
Conclusions are as follows:
1. The classical FVM is capable of stably simulating
the large deformation of structures by introducing
the appropriate artificial viscosity and spatial filter
which avoid hourglass mode.
2. The large deformation is capable by IB method in
 ~ .
m~~r~~:
: 1
n I 
11 
the flow computation. Thus we confirm that the
proposed weakcoupling method is useful for the
FSI problem.
References
[1] Q. Zhang and T. Hisada, Analysis of fluidstructure
interaction problems with structural buckling and
large domain changes by. ALE finite element
method, Comput. Methods Appl. Mech. Eng.,
Vol.190, pp.63416357, 2001.
[2] B. Hubner, E. Walhor and D. A.Dinkler, Mono
lithic approach to fluidstructure interaction using
spacetime finite elements, Comput. Meth. Appl.
Mech. Eng. Vol.193, pp.20872104, 2004.
[3] F. Blom, A monolithical fluidstructure interac
tion algorithm applied to the piston problem, Com
put. Methods Appl. Mech. Eng. Vol.167, pp.369
391,1998.
[4] S. Piperno, Explicit/implicit fluid/structure stag
gered procedures with a structural predictor and
fluid subcycling for 2D inviscid aeroelastic simu
lations, Int. J. Num. Meth. Fluids Vol.25, pp.1207
1226, 1997.
[5] C. Farhat, M. Lesoinne, Two efficient staggered
algorithms for the serial and parallel solution
of threedimensional nonlinear transient aeroelas
tic problems, Comp. Meth. Appl. Mech. Eng.,
Vol. 182, pp.499515, 2000.
[6] C.S.Peskin, The fluid dynamics of heart valves:
Experi mental, theoretical and computational
methods, Ann. Rev. Fluid Mech. Vol.14, pp.235
259, 1981.
[7] E.A. Fadlun, R. Verzicco, P Orlandi and J. Mohd
Yusof, Combined ImmersedBoundary Finite
Difference Methods for ThreeDimensional Com
plex Flow Simulations, J. Comp. Physics, Vol.161,
pp.3560, 2000.
[8] Mark L. Wilkins, Computer Simulation of Dy
namic Phenomena, Springer Verlag GmbH, 2006.
[9] W. H. Press, B. P. Flannery, S.A. Teukolsky, W.
T. Vetterling, NUMERICAL RECIPES in FOR
TRAN The Art of Scientific Computing, Cam
bridge University Press, pp.644649, 1992.
[10] A.K. Slone K. Pericleous, C. Bailey, M. Cross,
Dynamic fluidstructure interaction using finite
volume unstructured mesh procedures, Computers
& structures, Vol.80, pp.371390, 2002.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
[11] X. Lv, Y. Zhao, X.Y. Huang, G.H. Xia and X.H. Su,
A matrixfree implicit unstructured multigrid finite
volume method for simulating structural dynamics
and fluidstructure interaction, J. Comput. Physics,
Vol.225, pp.120144, 2007.
[12] Cyril M. Harris, Shock and Vibration Handbook,
McgrawHill Handbooks, 1961.
[13] Kim, J. and Moin, P, Application of a fractional
step method to incompressible flows, J. Comp.
Physics, Vol.59, pp. 308323, 1985.
